Extracting Gravitational Waves and Other Quantities from Numerical Spacetimes

Gabrielle Allen

March 27, 2024

Abstract

NB: This documentation is taken from the Extract thorn on which WaveExtractL is based. There may be some differences between WaveExtractL and Extract, which are not documented here.

1 Introduction

Thorn Extract calculates first order gauge invariant waveforms from a numerical spacetime, under the basic assumption that, at the spheres of extract the spacetime is approximately Schwarzschild. In addition, other quantities such as mass, angular momentum and spin can be determined.

This thorn should not be used blindly, it will always return some waveform, however it is up to the user to determine whether this is the appropriate expected first order gauge invariant waveform.

2 Physical System

2.1 Wave Forms

Assume a spacetime \(g_{\alpha \beta }\) which can be written as a Schwarzschild background \(g_{\alpha \beta }^{Schwarz}\) with perturbations \(h_{\alpha \beta }\): \begin {equation} g_{\alpha \beta } = g^{Schwarz}_{\alpha \beta } + h_{\alpha \beta } \end {equation} with \begin {equation} \{g^{Schwarz}_{\alpha \beta }\}(t,r,\theta ,\phi ) = \left ( \begin {array}{cccc} -S & 0 & 0 & 0 \\ 0 & S^{-1} & 0 & 0 \\ 0 & 0 & r^2 & 0 \\ 0 & 0 & 0 & r^2 \sin ^2\theta \end {array}\right ) \qquad S(r)=1-\frac {2M}{r} \end {equation} The 3-metric perturbations \(\gamma _{ij}\) can be decomposed using tensor harmonics into \(\gamma _{ij}^{lm}(t,r)\) where

            ∑∞  ∑l
γij(t,r,𝜃,ϕ) =        γlmij (t,r)
            l=0m= −l

and

            ∑6
γij(t,r,𝜃,ϕ) =   pk(t,r)Vk (𝜃,ϕ )
            k=0

where \(\{{\bf V}_k\}\) is some basis for tensors on a 2-sphere in 3-D Euclidean space. Working with the Regge-Wheeler basis (see Section 6) the 3-metric is then expanded in terms of the (six) standard Regge-Wheeler functions \(\{c_1^{\times lm}, c_2^{\times lm}, h_1^{+lm}, H_2^{+lm}, K^{+lm}, G^{+lm}\}\) [19][16]. Where each of the functions is either odd (\(\times \)) or even (\(+\)) parity. The decomposition is then written

\begin {eqnarray} \gamma _{ij}^{lm} & = & c_1^{\times lm}(\hat {e}_1)_{ij}^{lm} + c_2^{\times lm}(\hat {e}_2)_{ij}^{lm} \nonumber \\ & + & h_1^{+lm}(\hat {f}_1)_{ij}^{lm} + A^2 H_2^{+lm}(\hat {f}_2)_{ij}^{lm} + R^2 K^{+lm}(\hat {f}_3)_{ij}^{lm} + R^2 G^{+lm}(\hat {f}_4)_{ij}^{lm} \end {eqnarray}

which we can write in an expanded form as

\begin {eqnarray} \gamma _{rr}^{lm} & = & A^2 H_2^{+lm} \Y \\ \gamma _{r\t }^{lm} & = & - c_1^{\times lm} \frac {1}{\s } \Yp +h_1^{+lm}\Yt \\ \gamma _{r\p }^{lm} & = & c_1^{\times lm} \s \Yt + h_1^{+lm}\Yp \\ \gamma _{\t \t }^{lm} & = & c_2^{\times lm}\frac {1}{\s }(\Ytp -\cot \t \Yp ) + R^2 K^{+lm}\Y + R^2 G^{+lm} \Ytt \\ \gamma _{\t \p }^{lm} & = & -c_2^{\times lm}\s \frac {1}{2} \left ( \Ytt -\cot \t \Yt -\frac {1}{\sin ^2\theta }\Y \right ) + R^2 G^{+lm}(\Ytp -\cot \t \Yp ) \\ \gamma _{\p \p }^{lm} & = & -\s c_2^{\times lm} (\Ytp - \cot \t \Yp ) +R^2 K^{+lm}\sin ^2\t \Y +R^2 G^{+lm} (\Ypp +\s \c \Yt ) \end {eqnarray}

A similar decomposition allows the four gauge components of the 4-metric to be written in terms of three even-parity variables \(\{H_0,H_1,h_0\}\) and the one odd-parity variable \(\{c_0\}\)

\begin {eqnarray} g_{tt}^{lm} & = & N^2 H_0^{+lm} \Y \\ g_{tr}^{lm} & = & H_1^{+lm} \Y \\ g_{t\t }^{lm} & = & h_0^{+lm} \Yt - c_0^{\times lm}\frac {1}{\s }\Yp \\ g_{t\p }^{lm} & = & h_0^{+lm} \Yp + c_0^{\times lm} \s \Yt \end {eqnarray}

Also from \(g_{tt}=-\alpha ^2+\beta _i\beta ^i\) we have \begin {equation} \alpha ^{lm} = -\frac {1}{2}NH_0^{+lm}Y_{lm} \end {equation} It is useful to also write this with the perturbation split into even and odd parity parts:

                 ∑          ∑
gαβ = gbaαcβkground+    glαmβ,odd+    glαmβ,even
                 l,m          l,m

where (dropping some superscripts)

\begin {eqnarray*} \{g_{\alpha \beta }^{odd}\} &=& \left ( \begin {array}{cccc} 0 & 0 & - c_0\frac {1}{\s }\Yp & c_0 \s \Yt \\ . & 0 & - c_1\frac {1}{\s } \Yp & c_1 \s \Yt \\ . & . & c_2\frac {1}{\s }(\Ytp -\cot \t \Yp ) & c_2\frac {1}{2} \left (\frac {1}{\s } \Ypp +\c \Yt -\s \Ytt \right ) \\ .&.&.&c_2 (-\s \Ytp +\c \Yp ) \end {array} \right ) \\ \{g_{\alpha \beta }^{even}\} &=& \left ( \begin {array}{cccc} N^2 H_0\Y & H_1\Y & h_0\Yt & h_0 \Yp \\ . & A^2H_2\Y & h_1\Yt & h_1 \Yp \\ . & . & R^2K\Y +r^2G\Ytt & R^2(\Ytp -\cot \t \Yp ) \\ . & . & . & R^2 K\sin ^2\t \Y +R^2G(\Ypp +\s \c \Yt ) \end {array} \right ) \end {eqnarray*}

Now, for such a Schwarzschild background we can define two (and only two) unconstrained gauge invariant quantities \(Q^{\times }_{lm}=Q^{\times }_{lm}(c_1^{\times lm},c_2^{\times lm})\) and \(Q^{+}_{lm}=Q^{+}_{lm}(K^{+ lm},G^{+ lm},H_2^{+lm},h_1^{+lm})\), which from [3] are

\begin {eqnarray} Q^{\times }_{lm} & = & \sqrt {\frac {2(l+2)!}{(l-2)!}}\left [c_1^{\times lm} + \frac {1}{2}\left (\partial _r c_2^{\times lm} - \frac {2}{r} c_2^{\times lm}\right )\right ] \frac {S}{r} \\ Q^{+}_{lm} & = & \frac {1}{\Lambda }\sqrt {\frac {2(l-1)(l+2)}{l(l+1)}} (4rS^2 k_2+l(l+1)r k_1) \\ & \equiv & \frac {1}{\Lambda }\sqrt {\frac {2(l-1)(l+2)}{l(l+1)}} \left (l(l+1)S(r^2\partial _r G^{+lm}-2h_1^{+lm})+ 2rS(H_2^{+lm}-r\partial _r K^{+lm})+\Lambda r K^{+lm}\right ) \end {eqnarray}

where

\begin {eqnarray} k_1 & = & K^{+lm} + \frac {S}{r}(r^2\partial _r G^{+lm} - 2h^{+lm}_1) \\ k_2 & = & \frac {1}{2S} \left [H^{+lm}_2-r\partial _r k_1-\left (1-\frac {M}{rS}\right ) k_1 + S^{1/2}\partial _r (r^2 S^{1/2} \partial _r G^{+lm}-2S^{1/2}h_1^{+lm})\right ] \\ &\equiv & \frac {1}{2S}\left [H_2-rK_{,r}-\frac {r-3M}{r-2M}K\right ] \end {eqnarray}

NOTE: These quantities compare with those in Moncrief [16] by

\begin {eqnarray*} \mbox {Moncriefs odd parity Q: }\qquad Q^\times _{lm} &=& \sqrt {\frac {2(l+2)!}{(l-2)!}}Q \\ \mbox {Moncriefs even parity Q: } \qquad Q^+_{lm} &=& \sqrt {\frac {2(l-1)(l+2)}{l(l+1)}}Q \end {eqnarray*}

Note that these quantities only depend on the purely spatial Regge-Wheeler functions, and not the gauge parts. (In the Regge-Wheeler and Zerilli gauges, these are just respectively (up to a rescaling) the Regge-Wheeler and Zerilli functions). These quantities satisfy the wave equations

\begin {eqnarray*} &&(\partial ^2_t-\partial ^2_{r^*})Q^\times _{lm}+S\left [\frac {l(l+1)}{r^2}-\frac {6M}{r^3} \right ]Q^{\times }_{lm} = 0 \\ &&(\partial ^2_t-\partial ^2_{r^*})Q^+_{lm}+S\left [ \frac {1}{\Lambda ^2}\left (\frac {72M^3}{r^5}-\frac {12M}{r^3}(l-1)(l+2)\left (1-\frac {3M}{r}\right ) \right )+\frac {l(l-1)(l+1)(l+2)}{r^2\Lambda }\right ]Q^+_{lm}=0 \end {eqnarray*}

where

\begin {eqnarray*} \Lambda &=& (l-1)(l+2)+6M/r \\ r^* &=& r+2M\ln (r/2M-1) \end {eqnarray*}

3 Numerical Implementation

The implementation assumes that the numerical solution, on a Cartesian grid, is approximately Schwarzshild on the spheres of constant \(r=\sqrt (x^2+y^2+z^2)\) where the waveforms are extracted. The general procedure is then:

3.1 Project onto Spheres of Constant Radius

This is performed by interpolating the metric components, and if needed the conformal factor, onto the spheres. Although 2-spheres are hardcoded, the source code could easily be changed here to project onto e.g. 2-ellipsoids.

3.2 Calculate Radial Transformation

The areal coordinate \(\hat {r}\) of each sphere is calculated by \begin {equation} \hat {r} = \hat {r}(r) = \left [ \frac {1}{4\pi } \int \sqrt {\gamma _{\t \t } \gamma _{\p \p }}d\t d\p \right ]^{1/2} \end {equation} from which \begin {equation} \frac {d\hat {r}}{d\eta } = \frac {1}{16\pi \hat {r}} \int \frac {\gamma _{\t \t ,\eta }\gamma _{\p \p }+\gamma _{\t \t }\gamma _{\p \p ,\eta }} {\sqrt {\gamma _{\t \t }\gamma _{\p \p }}} \ d\t d\p \end {equation} Note that this is not the only way to combine metric components to get the areal radius, but this one was used because it gave better values for extracting close to the event horizon for perturbations of black holes.

3.3 Calculate \(S\) factor and Mass Estimate

\begin {equation} S(\hat {r}) = \left (\frac {\partial \hat {r}}{\partial r}\right )^2 \int \gamma _{rr} \ d\t d\p \end {equation}

\begin {equation} M(\hat {r}) = \hat {r}\frac {1-S}{2} \end {equation}

3.4 Calculate Regge-Wheeler Variables

\begin {eqnarray*} c_1^{\times lm} &=& \frac {1}{l(l+1)} \int \frac {\gamma _{\hat {r}\p }Y^*_{lm,\t } -\gamma _{\hat {r}\t } Y^*_{lm,\p } } {\s }d\Omega \\ c_2^{\times lm} & = & -\frac {2}{l(l+1)(l-1)(l+2)} \int \left \{ \left (-\frac {1}{\sin ^2\t }\gamma _{\t \t }+\frac {1} {\sin ^4\t }\gamma _{\p \p }\right ) (\s Y^*_{lm,\t \p }-\c Y^*_{lm,\p }) \right . \\ &&\left . + \frac {1}{\s } \gamma _{\t \p } (Y^*_{lm,\t \t }-\cot \t Y^*_{lm,\t } -\frac {1}{\sin ^2\t }Y^*_{lm,\p \p }) \right \}d\Omega \\ h_1^{+lm} &=& \frac {1}{l(l+1)} \int \left \{ \gamma _{\hat {r}\t } Y^*_{lm,\t } + \frac {1}{\sin ^2\t } \gamma _{\hat {r}\p }Y^*_{lm,\p }\right \} d\Omega \\ H_2^{+lm} &=& S \int \gamma _{\hat {r}\hat {r}} \Ys d\Omega \\ K^{+lm} &=& \frac {1}{2\hat {r}^2} \int \left (\gamma _{\t \t }+ \frac {1}{\sin ^2\t }\gamma _{\p \p }\right )\Ys d\Omega \\ &&+\frac {1}{2\hat {r}^2(l-1)(l+2)} \int \left \{ \left (\gamma _{\t \t }-\frac {\gamma _{\p \p }}{\sin ^2\t }\right ) \left (Y^*_{lm,\t \t }-\cot \t Y^*_{lm,\t }-\frac {1}{\sin ^2\t } Y^*_{lm,\p \p }\right ) \right . \\ &&\left . + \frac {4}{\sin ^2\t }\gamma _{\t \p }(Y^*_{lm,\t \p }-\cot \t Y^*_{lm,\p }) \right \} d\Omega \\ G^{+lm} &=& \frac {1}{\hat {r}^2 l(l+1)(l-1)(l+2)} \int \left \{ \left (\gamma _{\t \t }-\frac {\gamma _{\p \p }}{\sin ^2\t }\right ) \left (Y^*_{lm,\t \t }-\cot \t Y^*_{lm,\t }-\frac {1}{\sin ^2\t } Y^*_{lm,\p \p }\right ) \right . \\ &&\left . +\frac {4}{\sin ^2\t }\gamma _{\t \p }(Y^*_{lm,\t \p }-\cot \t Y^*_{lm,\p }) \right \}d\Omega \end {eqnarray*}

where

\begin {eqnarray} \gamma _{\hat {r}\hat {r}} & = & \frac {\partial r}{\partial \hat {r}} \frac {\partial r}{\partial \hat {r}} \gamma _{rr} \\ \gamma _{\hat {r}\t } & = & \frac {\partial r}{\partial \hat {r}} \gamma _{r\t } \\ \gamma _{\hat {r}\p } & = & \frac {\partial r}{\partial \hat {r}} \gamma _{r\p } \end {eqnarray}

3.5 Calculate Gauge Invariant Quantities

\begin {eqnarray} Q^{\times }_{lm} & = & \sqrt {\frac {2(l+2)!}{(l-2)!}}\left [c_1^{\times lm} + \frac {1}{2}\left (\partial _{\hat {r}} c_2^{\times lm} - \frac {2}{\hat {r}} c_2^{\times lm}\right )\right ] \frac {S}{\hat {r}} \\ Q^{+}_{lm} & = & \frac {1}{(l-1)(l+2)+6M/\hat {r}}\sqrt {\frac {2(l-1)(l+2)}{l(l+1)}} (4\hat {r}S^2 k_2+l(l+1)\hat {r} k_1) \end {eqnarray}

where

\begin {eqnarray} k_1 & = & K^{+lm} + \frac {S}{\hat {r}}(\hat {r}^2\partial _{\hat {r}} G^{+lm} - 2h^{+lm}_1) \\ k_2 & = & \frac {1}{2S} [H^{+lm}_2-\hat {r}\partial _{\hat {r}} k_1-(1-\frac {M}{\hat {r}S}) k_1 + S^{1/2}\partial _{\hat {r}} (\hat {r}^2 S^{1/2} \partial _{\hat {r}} G^{+lm}-2S^{1/2}h_1^{+lm} \end {eqnarray}

4 Using This Thorn

Use this thorn very carefully. Check the validity of the waveforms by running tests with different resolutions, different outer boundary conditions, etc to check that the waveforms are consistent.

4.1 Basic Usage

4.2 Output Files

Although Extract is really an ANALYSIS thorn, at the moment it is scheduled at POSTSTEP, with the iterations at which output is performed determined by the parameter itout. Output files from Extract are always placed in the main output directory defined by CactusBase/IOUtil.

Output files are generated for each detector (2-sphere) used, and these detectors are identified in the name of each output file by R1, R2, ….

The extension denotes whether coordinate time (ṫl) or proper time (˙u  l) is used for the first column.

5 History

Much of the source code for Extract comes from a code written outside of Cactus for extracting waveforms from data generated by the NCSA G-Code for compare with linear evolutions of waveforms extracted from the Cauchy initial data. This work was carried out in collaboration with Karen Camarda and Ed Seidel.

6 Appendix: Regge-Wheeler Harmonics

\begin {eqnarray*} (\hat {e}_1)^{lm} &=& \left ( \begin {array}{ccc} 0 & -\frac {1}{\s }\Yp & \s \Yt \\ . & 0 & 0 \\ . & 0 & 0 \end {array}\right ) \\ (\hat {e}_2)^{lm} &=& \left ( \begin {array}{ccc} 0 & 0 & 0 \\ 0 & \frac {1}{\s }(\Ytp -\cot \t \Yp ) & . \\ 0 & -\frac {\s }{2}[\Ytt -\cot \t \Yt -\frac {1}{\sin ^2\t }\Ypp ] & -\s [\Ytp -\cot \t \Yp ] \end {array}\right ) \\ (\hat {f}_1)^{lm} &=& \left ( \begin {array}{ccc} 0 & \Yt & \Yp \\ . & 0 & 0 \\ . & 0 & 0 \end {array}\right ) \\ (\hat {f}_2)^{lm} &=& \left ( \begin {array}{ccc} \Y & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end {array}\right ) \\ (\hat {f}_3)^{lm} &=& \left ( \begin {array}{ccc} 0 & 0 & 0 \\ 0 & \Y & 0 \\ 0 & 0 & \sin ^2\t \Y \end {array}\right ) \\ (\hat {f}_4)^{lm} &=& \left ( \begin {array}{ccc} 0 & 0 & 0 \\ 0 & \Ytt & . \\ 0 & \Ytp -\cot \t \Yp & \Ypp + \s \c \Yt \end {array}\right ) \end {eqnarray*}

7 Appendix: Transformation Between Cartesian and Spherical Coordinates

First, the transformations between metric components in \((x,y,z)\) and \((r,\t ,\p )\) coordinates. Here, \(\rho =\sqrt {x^2+y^2}=r\s \),

\begin {eqnarray*} \frac {\partial x}{\partial r} &=& \sin \t \cos \p = \frac {x}{r} \\ \frac {\partial y}{\partial r} &=& \sin \t \sin \p = \frac {y}{r} \\ \frac {\partial z}{\partial r} &=& \cos \t = \frac {z}{r} \\ \frac {\partial x}{\partial \t } &=& r\cos \t \cos \p = \frac {xz}{\rho } \\ \frac {\partial y}{\partial \t } &=& r\cos \t \sin \p = \frac {yz}{\rho } \\ \frac {\partial z}{\partial \t } &=& -r\sin \t = -\rho \\ \frac {\partial x}{\partial \p } &=& -r\sin \t \sin \p = -y \\ \frac {\partial y}{\partial \p } &=& r\sin \t \cos \p = x \\ \frac {\partial z}{\partial \p } &=& 0 \end {eqnarray*}
\begin {eqnarray*} \gamma _{rr} &=& \frac {1}{r^2} (x^2\gamma _{xx}+ y^2\gamma _{yy}+ z^2\gamma _{zz}+ 2xy\gamma _{xy}+ 2xz\gamma _{xz}+ 2yz\gamma _{yz}) \\ \gamma _{r\t } &=& \frac {1}{r\rho } (x^2 z \gamma _{xx} +y^2 z \gamma _{yy} -z \rho ^2 \gamma _{zz} +2xyz \gamma _{xy} +x(z^2-\rho ^2)\gamma _{xz} +y(z^2-\rho ^2)\gamma _{yz}) \\ \gamma _{r\p } &=& \frac {1}{r} (-xy\gamma _{xx} +xy\gamma _{yy} +(x^2-y^2)\gamma _{xy} -yz \gamma _{xz} +xz\gamma _{yz}) \\ \gamma _{\t \t } &=& \frac {1}{\rho ^2} (x^2z^2\gamma _{xx} +2xyz^2\gamma _{xy} -2xz\rho ^2\gamma _{xz} +y^2z^2\gamma _{yy} -2yz\rho ^2\gamma _{yz} +\rho ^4\gamma _{zz}) \\ \gamma _{\t \p } &=& \frac {1}{\rho } (-xyz\gamma _{xx} +(x^2-y^2)z\gamma _{xy} +\rho ^2 y \gamma _{xz} +xyz\gamma _{yy} -\rho ^2 x \gamma _{yz}) \\ \gamma _{\p \p } &=& y^2\gamma _{xx} -2xy\gamma _{xy} +x^2\gamma _{yy} \end {eqnarray*}

or,

\begin {eqnarray*} \gamma _{rr}&=& \sin ^2\t \cos ^2\p \gamma _{xx} +\sin ^2\t \sin ^2\p \gamma _{yy} +\cos ^2\t \gamma _{zz} +2\sin ^2\theta \cos \p \sin \p \gamma _{xy} +2\sin \t \cos \t \cos \p \gamma _{xz} \\ && +2\s \c \sin \p \gamma _{yz} \\ \gamma _{r\t }&=& r(\s \c \cos ^2\phi \gamma _{xx} +2*\s \c \sin \p \cos \p \gamma _{xy} +(\cos ^2\t -\sin ^2\t )\cos \p \gamma _{xz} +\s \c \sin ^2\p \gamma _{yy} \\ && +(\cos ^2\t -\sin ^2\t )\sin \p \gamma _{yz} -\s \c \gamma _{zz}) \\ \gamma _{r\p }&=& r\s (-\s \sin \p \cos \p \gamma _{xx} -\s (\sin ^2\p -\cos ^2\p )\gamma _{xy} -\c \sin \p \gamma _{xz} +\s \sin \p \cos \p \gamma _{yy} \\ && +\c \cos \p \gamma _{yz}) \\ \gamma _{\t \t }&=& r^2(\cos ^2\t \cos ^2\p \gamma _{xx} +2\cos ^2\t \sin \p \cos \p \gamma _{xy} -2\s \c \cos \p \gamma _{xz} +\cos ^2\t \sin ^2\p \gamma _{yy} \\ && -2\s \c \sin \p \gamma _{yz} +\sin ^2\t \gamma _{zz}) \\ \gamma _{\t \p }&=& r^2\s (-\c \sin \p \cos \p \gamma _{xx} -\c (\sin ^2\p -\cos ^2\p )\gamma _{xy} +\s \sin \p \gamma _{xz} +\c \sin \p \cos \p \gamma _{yy} \\ && -\s \cos \p \gamma _{yz}) \\ \gamma _{\p \p }&=& r^2\sin ^2\t (\sin ^2\p \gamma _{xx} -2\sin \p \cos \p \gamma _{xy} +\cos ^2\phi \gamma _{yy}) \end {eqnarray*}

We also need the transformation for the radial derivative of the metric components:

\begin {eqnarray*} \gamma _{rr,\eta }&=& \sin ^2\t \cos ^2\p \gamma _{xx,\eta } +\sin ^2\t \sin ^2\p \gamma _{yy,\eta } +\cos ^2\t \gamma _{zz,\eta } +2\sin ^2\theta \cos \p \sin \p \gamma _{xy,\eta } \\ && +2\sin \t \cos \t \cos \p \gamma _{xz,\eta } +2\s \c \sin \p \gamma _{yz,\eta } \\ \gamma _{r\t ,\eta }&=& \frac {1}{r}\gamma _{r\t }+ r(\s \c \cos ^2\phi \gamma _{xx,\eta } +\s \c \sin \p \cos \p \gamma _{xy,\eta } +(\cos ^2\t -\sin ^2\t )\cos \p \gamma _{xz,\eta } \\ && +\s \c \sin ^2\p \gamma _{yy,\eta } +(\cos ^2\t -\sin ^2\t )\sin \p \gamma _{yz,\eta } -\s \c \gamma _{zz,\eta }) \\ \gamma _{r\p ,\eta }&=& \frac {1}{r}\gamma _{r\p }+ r\s (-\s \sin \p \cos \p \gamma _{xx,\eta } -\s (\sin ^2\p -\cos ^2\p )\gamma _{xy,\eta } -\c \sin \p \gamma _{xz,\eta } \\ && +\s \sin \p \cos \p \gamma _{yy,\eta } +\c \cos \p \gamma _{yz,\eta }) \\ \gamma _{\t \t ,\eta }&=& \frac {2}{r}\gamma _{\t \t }+ r^2(\cos ^2\t \cos ^2\p \gamma _{xx,\eta } +2\cos ^2\t \sin \p \cos \p \gamma _{xy,\eta } -2\s \c \cos \p \gamma _{xz,\eta } \\ && +\cos ^2\t \sin ^2\p \gamma _{yy,\eta } -2\s \c \sin \p \gamma _{yz,\eta } +\sin ^2\t \gamma _{zz,\eta }) \\ \gamma _{\t \p ,\eta }&=& \frac {2}{r}\gamma _{\t \p }+ r^2\s (-\c \sin \p \cos \p \gamma _{xx,\eta } -\c (\sin ^2\p -\cos ^2\p )\gamma _{xy,\eta } +\s \sin \p \gamma _{xz,\eta } \\ && +\c \sin \p \cos \p \gamma _{yy,\eta } -\s \cos \p \gamma _{yz,\eta }) \\ \gamma _{\p \p ,\eta }&=& \frac {2}{r}\gamma _{\p \p }+ r^2\sin ^2\t (\sin ^2\p \gamma _{xx,\eta } -2\sin \p \cos \p \gamma _{xy,\eta } +\cos ^2\phi \gamma _{yy,\eta }) \end {eqnarray*}

8 Appendix: Integrations Over the 2-Spheres

This is done by using Simpson’s rule twice. Once in each coordinate direction. Simpson’s rule is \begin {equation} \int ^{x_2}_{x_1} f(x) dx = \frac {h}{3} [f_1+4f_2+2f_3+4f_4+\ldots +2f_{N-2}+4 f_{N-1}+f_N] +O(1/N^4) \end {equation} \(N\) must be an odd number.

References

[1]   Abrahams A.M. & Cook G.B. “Collisions of boosted black holes: Perturbation theory predictions of gravitational radiation” Phys. Rev. D 50 R2364-R2367 (1994).

[2]   Abrahams A.M., Shapiro S.L. & Teukolsky S.A. “Calculation of gravitational wave forms from black hole collisions and disk collapse: Applying perturbation theory to numerical spacetimes” Phys. Rev. D. 51 4295 (1995).

[3]   Abrahams A.M. & Price R.H. “Applying black hole perturbation theory to numerically generated spacetimes” Phys. Rev. D. 53 1963 (1996).

[4]   Abrahams A.M. & Price R.H. “Black-hole collisions from Brill-Lindquist initial data: Predictions of perturbation theory” Phys. Rev. D. 53 1972 (1996).

[5]   Abramowitz, M. & Stegun A. “Pocket Book of Mathematical Functions (Abridged Handbook of Mathematical Functions”, Verlag Harri Deutsch (1984).

[6]   Andrade Z., & Price R.H. “Head-on collisions of unequal mass black holes: Close-limit predictions”, preprint (1996).

[7]   Anninos P., Price R.H., Pullin J., Seidel E., and Suen W-M. “Head-on collision of two black holes: Comparison of different approaches” Phys. Rev. D. 52 4462 (1995).

[8]   Arfken, G. “Mathematical Methods for Physicists”, Academic Press (1985).

[9]   Baker J., Abrahams A., Anninos P., Brant S., Price R., Pullin J. & Seidel E. “The collision of boosted black holes” (preprint) (1996).

[10]   Baker J. & Li C.B. “The two-phase approximation for black hole collisions: Is it robust” preprint (gr-qc/9701035), (1997).

[11]   Brandt S.R. & Seidel E. “The evolution of distorted rotating black holes III: Initial data” (preprint) (1996).

[12]   Cunningham C.T., Price R.H., Moncrief V., “Radiation from collapsing relativistic stars. I. Linearized Odd-Parity Radiation” Ap. J. 224 543-667 (1978).

[13]   Cunningham C.T., Price R.H., Moncrief V., “Radiation from collapsing relativistic stars. I. Linearized Even-Parity Radiation” Ap. J. 230 870-892 (1979).

[14]   Landau L.D. & Lifschitz E.M., “The Classical Theory of Fields” (4th Edition), Pergamon Press (1980).

[15]   Mathews J. “”, J. Soc. Ind. Appl. Math. 10 768 (1962).

[16]   Moncrief V. “Gravitational perturbations of spherically symmetric systems. I. The exterior problem” Annals of Physics 88 323-342 (1974).

[17]   Press W.H., Flannery B.P., Teukolsky S.A., & Vetterling W.T., “Numerical Recipes, The Art of Scientific Computing” Cambridge University Press (1989).

[18]   Price R.H. & Pullin J. “Colliding black holes: The close limit”, Phys. Rev. Lett. 72 3297-3300 (1994).

[19]   Regge T., & Wheeler J.A. “Stability of a Schwarzschild Singularity”, Phys. Rev. D 108 1063 (1957).

[20]   Seidel E. Phys Rev D. 42 1884 (1990).

[21]   Thorne K.S., “Multipole expansions of gravitational radiation”, Rev. Mod. Phys. 52 299 (1980).

[22]   Vishveshwara C.V., “Stability of the Schwarzschild metric”, Phys. Rev. D. 1 2870, (1970).

[23]   Zerilli F.J., “Tensor harmonics in canonical form for gravitational radiation and other applications”, J. Math. Phys. 11 2203, (1970).

[24]   Zerilli F.J., “Gravitational field of a particle falling in a Schwarzschild geometry analysed in tensor harmonics”, Phys. Rev. D. 2 2141, (1970).

9 Parameters




active
Scope: private  BOOLEAN



Description: Should we run at all?



  Default: yes






calc_when_necessary
Scope: private  BOOLEAN



Description: Start calculation at T = R_detector - 50 and stop after T_merger+Ringdown_margin+R_detector



  Default: no






cartoon
Scope: private  BOOLEAN



Description: use cartoon



  Default: no






cartoon_grid_acc
Scope: private  REAL



Description: accuracy to use for dx=dy=dz check



Range   Default: 1e-2
0:*
any non-negative value






cauchy_radius_dr
Scope: private  REAL



Description: seperation of 2 detectors for the extraction



Range   Default: -1
-1:*
-1 means: let the code figure it out






cauchy_radius_end_coord
Scope: private  REAL



Description: where to end Cauchy line in coordinates



Range   Default: -2
-2:*
-2 means deactive (ie use percentage), -1 means go up to maximum






cauchy_radius_end_factor
Scope: private  REAL



Description: where to end Cauchy line as factor of grid size



Range   Default: 1
-1:1
-1 means deactive (ie use coordinate), 1 means go up to maximum






cauchy_radius_extreme
Scope: private  BOOLEAN



Description: use min. radius along axis, dr=dx



  Default: (none)






cauchy_radius_start_coord
Scope: private  REAL



Description: where to start Cauchy line in coordinates



Range   Default: -2
-2:*
-2 means deactive (ie use percentage), -1 means start from first grid point






cauchy_radius_start_factor
Scope: private  REAL



Description: where to start Cauchy line as factor of grid size



Range   Default: (none)
-1:1
-1 means deactive (ie use coordinate), 0 means first grid point






cauchy_time_id
Scope: private  BOOLEAN



Description: compute initial data (ID) for dt_Zerilli?



  Default: (none)






check_rmax
Scope: private  BOOLEAN



Description: Whether to check for rmax or not (must be set to no for Llama



  Default: no






corotate
Scope: private  BOOLEAN



Description: do we corotate? give omega if so



  Default: no






corotate3d
Scope: private  BOOLEAN



Description: do we corotate in 3D (easy slow way)? give omega if so



  Default: no






detection_mode
Scope: private  KEYWORD



Description: Give Specific locations or let the code place the detectors



Range   Default: specific detectors
Cauchy extraction
a lot of detectors in some specified range
specific detectors
give some values






detector_metric
Scope: private  STRING



Description: Metric w.r.t which we extract



Range   Default: admbase::metric
.*
Any group with storage on. Must be in order xx,xy,xz,yy,yz,zz






detector_radius
Scope: private  REAL



Description: Coordinate radius of detectors



Range   Default: (none)
0:*
maximum radius is checked at runtime, detector should be in the grid






drsch_dri_computation
Scope: private  KEYWORD



Description: How to calculate dr_Schwarzschild/dr_isotropic



Range   Default: average dr_gtt dr_gpp
see [1] below
”average using drsch_dri=int(1/2(dr _gtt+dr_gpp/sinˆ2  (th eta)))”
dr_gtt
”drsch_dri=int(dr_gt t)”
dr_gpp
”drsch_dri=int(dr_gp p/sinˆ
2  (theta))”



[1]

average dr\_gtt dr\_gpp




integration_method
Scope: private  KEYWORD



Description: which method to use for the integration



Range   Default: Gauss
Gauss
error is O(1/eˆN  )
see [1] below
weight is one for all points, error is O(1/Nˆ2  )
open extended
”weight is 23/12,13/12,4/3,2/4, ...0,27/12, error O(1/Nˆ4  )”



[1]

extended midpoint rule




interpolation_operator
Scope: private  STRING



Description: What interpolator should be used to interpolate onto the surface



Range   Default: Hermite polynomial interpolation
.+
A valid interpolator name, see Thorn AEILocalInterp






interpolation_order
Scope: private  INT



Description: interpolation order for projection onto the sphere



Range   Default: 2
1:*
Positive Please






interpolation_stencil
Scope: private  INT



Description: interpolation stencil, this is needed to work out how far out you can place the detectors. It depends on the interpolator and the interpolation order used.



Range   Default: 4
1:*
Positive Please






l_mode
Scope: private  INT



Description: all modes: Up to which l mode to calculate/ spefic mode: which l mode to extract



Range   Default: 2
2:*
Positive Please, note that l=0,1 are gauge dependent and not implemented






m_mode
Scope: private  INT



Description: all modes: Up to which m mode to calculate/ specific mode: which m mode to extract



Range   Default: (none)
0:*
Positive Please






make_interpolator_warnings_fatal
Scope: private  BOOLEAN



Description: Report interpolator warnings as level-0 error messages (and abort the run) ?



  Default: no






maximum_detector_number
Scope: private  INT



Description: How many detectors do you want. NOTE: This also fixes the number of detectors in Cauchy mode!



Range   Default: 9
0:100
Positive Please, less than hard limit FIXME






maxnphi
Scope: private  INT



Description: how many points in phi direction at most



Range   Default: 100
0:*
Positive Please - even or odd depends on integration scheme used






maxntheta
Scope: private  INT



Description: how many points in theta direction at most



Range   Default: 100
0:*
Positive Please - even or odd depends on integration scheme used






mode_type
Scope: private  KEYWORD



Description: Which type of mode extraction do we have



Range   Default: all modes
all modes
Extract all modes up to (l_mode, m_mode).
specific mode
Select one specific (l_mode, m_mode) mode






nphi
Scope: private  INT



Description: how many points in phi direction



Range   Default: (none)
0:*
Positive Please - even or odd depends on integration scheme used






ntheta
Scope: private  INT



Description: how many points in theta direction



Range   Default: (none)
0:*
Positive Please - even or odd depends on integration scheme used






observers
Scope: private  BOOLEAN



Description: do we use the observer thorn



  Default: no






origin_x
Scope: private  REAL



Description: origin in x direction



Range   Default: (none)
*:*






origin_y
Scope: private  REAL



Description: origin in y direction



Range   Default: (none)
*:*






origin_z
Scope: private  REAL



Description: origin in z direction



Range   Default: (none)
*:*






out_dir
Scope: private  STRING



Description: Output directory for Extract’s files, overrides IO::out_dir



Range   Default: (none)
.+
A valid directory name
ˆ
$
An empty string to choose the default from IO::out_dir






out_every
Scope: private  INT



Description: At which iterations should we be called



Range   Default: 1
*:*
negative means set via out_every_det






out_every_det
Scope: private  INT



Description: At which iterations should this detector be evaluated



Range   Default: 1
*:*






out_format
Scope: private  STRING



Description: Which format for Scalar floating-point number output



Range   Default: .13f
see [1] below
output with given precision in exponential / floating point notation



[1]

\^({\textbackslash}.[1]?[0-9])?[EGefg]\$




out_style
Scope: private  KEYWORD



Description: Which style for Scalar output



Range   Default: gnuplot
gnuplot
1D output readable by gnuplot
xgraph
1D output readable by xgraph






phicorotate
Scope: private  BOOLEAN



Description: undo corotation on phi itself



  Default: no






ringdown_margin
Scope: private  REAL



Description: The expected length of ringdown. If calc_when_necessary=yes, then T_merger+ringdown_margin+R_detector is the time when calculations are switched off



Range   Default: 150.0
0:*
Any positive number






rotation_omega
Scope: private  REAL



Description: omega from driftcorrect



Range   Default: 0.0
:






rsch2_computation
Scope: private  KEYWORD



Description: How to calculate (Schwarzschild_Radius) ˆ
2



Range   Default: aerial radius
aerial radius
”calculate invariant aerial radius: rˆ2  =int(sqrt(gtt*gpp -gtpˆ2  ))”
see [1] below
”assume Schwarzschild coordinates: rˆ2   =int(1/2 (gtt+gpp/sinˆ2
(theta ))”
Schwarzschild gtt
assume Schwarzschild coordinates: rˆ2  =int(gtt)
Schwarzschild gpp
”assume Schwarzschild coordinates: rˆ2  =int(gpp/sinˆ2  (th eta))”



[1]

average Schwarzschild metric




start_iteration
Scope: private  INT



Description: First iteration when we should be called



Range   Default: (none)
*:*






start_time
Scope: private  INT



Description: First time when we should be called



Range   Default: -1
*:*






subtract_spherical_background
Scope: private  BOOLEAN



Description: Subtract spherical background before calculation of Extraction quantities



  Default: yes






surface_index
Scope: private  INT



Description: surface that contains extraction sphere



Range   Default: -1
-1:






switch_output_format
Scope: private  INT



Description: How many output files do you suffer in your directory, when you have more detectors than this number, the file format is switched



Range   Default: 10
1:*
positive






use_conformal_factor
Scope: private  BOOLEAN



Description: Do we use the conformal factor (if possible) with this metric?



  Default: yes






use_spherical_surface
Scope: private  BOOLEAN



Description: use spherical surface thorn provided spheres



  Default: no






verbose
Scope: private  INT



Description: how verbose should the output be?



Range   Default: 1
0:10
higher number means more output






write_timer_info
Scope: private  BOOLEAN



Description: write timing information to stdout



  Default: no






nghostsphi
Scope: shared from SPHERICALSURFACE INT






nghoststheta
Scope: shared from SPHERICALSURFACE INT






nsurfaces
Scope: shared from SPHERICALSURFACE INT






sfpar_maxnphi
Scope: shared from SPHERICALSURFACE INT






sfpar_maxntheta
Scope: shared from SPHERICALSURFACE INT






symmetric_x
Scope: shared from SPHERICALSURFACE BOOLEAN






symmetric_y
Scope: shared from SPHERICALSURFACE BOOLEAN






symmetric_z
Scope: shared from SPHERICALSURFACE BOOLEAN



10 Interfaces

General

Implements:

waveextractl

Inherits:

grid

admbase

staticconformal

io

sphericalsurface

Grid Variables

10.0.1 PRIVATE GROUPS




  Group Names    Variable Names    Details   




gridsizes_group   compact0
int_nphi   descriptiongridsizes stored as variables
int_ntheta   dimensions1
  distributionCONSTANT
  group typeARRAY
  size100
  tagscheckpoint=”no”
  timelevels1
 variable typeINT




handles_group   compact0
sum_handle   descriptionhandles for reduction operators
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  tagscheckpoint=”no”
  timelevels1
 variable typeINT




do_nothing_group   compact0
do_nothing   descriptionif equal to 1
    descriptionthen WaveExtract won’t do anything (example: all detectors are out of range)
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  tagscheckpoint=”no”
  timelevels1
 variable typeINT




sym_factor_group   compact0
sym_factor   descriptionsymmmetry factor for integrals (depends on domain
    descriptionsym_factor=2 for bitant for example)
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  tagscheckpoint=”no”
  timelevels1
 variable typeREAL




interp_metric_arrays   compact0
gxxi   description2D (theta
    descriptionphi) Arrays for holding the metric and conformal factor and their first derivatives interpolated onto the extraction coordinate sphere
gxyi   dimensions2
gxzi   distributionDEFAULT
gyyi   group typeARRAY
gyzi   sizeMAXNTHETA
    sizeMAXNPHI
gzzi   tagscheckpoint=”no”
psii   timelevels1
dx_gxxi  variable typeREAL




surface_arrays   compact0
interp_x   description2D (theta
    descriptionphi) grid arrays for points on the sphere
interp_y   dimensions2
interp_z   distributionDEFAULT
psi_ext_deriv   group typeARRAY
ctheta   sizeMAXNTHETA
    sizeMAXNPHI
cphi   tagscheckpoint=”no”
sintheta   timelevels1
costheta  variable typeREAL








  Group Names     Variable Names    Details   




surface_integrands   compact0
weights   descriptionweights and temporary integrands
thetaweights   dimensions2
phiweights   distributionDEFAULT
int_tmp1   group typeARRAY
int_tmp2   sizeMAXNTHETA
    sizeMAXNPHI
int_tmp3   tagscheckpoint=”no”
int_tmp4   timelevels1
int_tmp5  variable typeREAL




schwarzschild_mass_radius_group   compact0
dtau_dt   descriptionSchwarzschild radius
    descriptionmass and assorted spherical background pieces
sph_grr   dimensions0
sph_gtt   distributionCONSTANT
sph_dr_gtt   group typeSCALAR
sph_gpp   tagscheckpoint=”no”
rsch2   timelevels1
rsch  variable typeREAL




schwarzschild_mass_radius_results_group   compact0
Schw_Masses   descriptioncontains Schwarzschild mass/radius from all detectors
Schw_Radii   dimensions1
  distributionCONSTANT
  group typeARRAY
  sizeMAXIMUM_DETECTOR_NUMBER
  tagscheckpoint=”no”
  timelevels1
 variable typeREAL




moncriefq_results_group   compact0
Qodd_Re_Array   descriptioncontains Moncrief Qeven
    descriptionQodd wave indicators from all detectors
Qeven_Re_Array   dimensions3
Qodd_Im_Array   distributionCONSTANT
Qeven_Im_Array   group typeARRAY
  sizeMAXIMUM_DETECTOR_NUMBER
    sizeL_MODE
  size2*M_MODE+1
  tagscheckpoint=”no”
  timelevels1
 variable typeREAL




l_m_modes_info_group   compact0
l_min   descriptionInformation about the modes used for extraction
l_max   dimensions0
l_step   distributionCONSTANT
m_min   group typeSCALAR
m_max   tagscheckpoint=”no”
m_step   timelevels1
max_det_no_param  variable typeINT




moncriefq   compact0
Qodd_Re   descriptionMoncrief Qeven
    descriptionQodd wave indicators
Qodd_Re   descriptionreal & imaginary part
Qodd_Im   dimensions2
Qeven_Re   distributionCONSTANT
Qeven_Im   group typeARRAY
  sizeL_MODE
    sizeM_MODE+1
  tagscheckpoint=”no”
  timelevels1
 variable typeREAL








  Group Names     Variable Names     Details   




metric_tmp   compact0
gxx_tmp   descriptiontemp metric for 3d rotation
gxy_tmp   dimensions3
gxz_tmp   distributionDEFAULT
gyy_tmp   group typeGF
gyz_tmp   tagstensortypealias=”DD_sym” checkpoint=”no” prolongation=”none”
gzz_tmp   timelevels1
 variable typeREAL




current_detector_group   compact0
current_detector   descriptionthe index number of the current detector
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  tagscheckpoint=”no”
  timelevels1
 variable typeINT




current_detector_radius_group   compact0
current_detector_radius  descriptioncoordinate radius of the current detector
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  tagscheckpoint=”no”
  timelevels1
 variable typeREAL




my_out_every_det my_out_every_det   compact0
  descriptionoutput frequency
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  tagscheckpoint=”no”
  timelevels1
 vararray_size100
 variable typeINT




Uses header:

carpetinterp2.hh

11 Schedule

This section lists all the variables which are assigned storage by thorn Llama/WaveExtractL. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage

 

Always: Conditional:
MoncriefQmetric_tmp
  my_out_every_det
  l_m_modes_info_group
  current_detector_group
  current_detector_radius_group
  sym_factor_group handles_group do_nothing_group
  gridsizes_group
  surface_arrays
  surface_integrands
  interp_metric_arrays
  Schwarzschild_Mass_Radius_group
  MoncriefQ_Results_group Schwarzschild_Mass_Radius_Results_group
  current_detector_group current_detector_radius_group
   

Scheduled Functions

CCTK_PARAMCHECK (conditional)

  wavextrl_paramcheck

  check parameters for waveextract

 

 Language:c
 Type: function

CCTK_STARTUP (conditional)

  wavextrl_startup

  register waveextract as an io method

 

 After: ioutil_startup
 Language:c
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_subtrsphermetric

  substract spherical background part of metric from metric

 

 After: wavextrl_schwarzmassrad
  Language:fortran
 Options: global
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_moncriefq

  compute moncrief qeven, qodd from regge wheeler quantities

 

 After: wavextrl_subtrsphermetric
 Language:fortran
 Options: global
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_writedata

  write out results to disk and stdout - one file for each detector and (l,m) mode

 

 After: wavextrl_moncriefq
 Language:c
 Options: global
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_adjustdetector

  decrease current_detector, go the the next detector

 

 After: wavextrl_writedata
  Language:c
 Options: global
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_adjustdetector

  decrease current_detector, go the the next detector

 

 After: wavextrl_moncriefq
 Language:c
 Options: global
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_timerinfo

  output timer info if requested

 

 After: wavextrl_adjustdetector
 Language:c
 Options: global
 Type: function

CCTK_ANALYSIS (conditional)

  wavextrl_writelotsofdata

  output one file per (l,m) mode, all detectors in one file

 

 After: wavextrl_calcsatdetector
  Language:c
 Options: global
 Type: function

CCTK_BASEGRID (conditional)

  wavextrl_init

  setup weights for integration

 

 After: spatialcoordinates
 Language:fortran
 Type: function

CCTK_BASEGRID (conditional)

  wavextrl_setup_detectors

  initial setup of all detectors

 

 After: spatialcoordinates
 Language:fortran
 Options: level
 Type: function

CCTK_BASEGRID (conditional)

  wavextrl_setup_sphericalsurface

  setup detectors from spherical surface

 

 After: sphericalsurface_setup
   wavextrl_resetcurrdet
  Language:fortran
 Options: local
 Type: function

CCTK_ANALYSIS (conditional)

  wavextrl_resetcurrdet

  reset the value of the current_detector, needed for the while loop next

 

 Before: wavextrl_calcsatdetector
  Language:c
 Options: global
 Type: function

CCTK_ANALYSIS (conditional)

  wavextrl_calcsatdetector

  calculations done for each detector, we loop over the detectors

 

 After: wavextrl_resetcurrdet
  Type: group
 While:waveextractl::current_detector

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_setupsphere

  setup sintheta, sinphi arrays

 

 Language:fortran
 Options: global
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_projectsphere

  interpolate 3d quantities into 2d grid arrays (on the sphere), project onto sphere

 

 After: wavextrl_setupsphere
 Language:c
 Options: global
 Type: function

WavExtrL_CalcsAtDetector (conditional)

  wavextrl_schwarzmassrad

  calculate schwarzschild radius and mass and spherical background

 

 After: wavextrl_projectsphere
 Language:fortran
 Options: global
 Type: function