## Extracting Gravitational Waves and Other Quantities from Numerical Spacetimes

June 14, 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

and

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:

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:

• Project the required metric components, and radial derivatives of metric components, onto spheres of constant coordinate radius (these spheres are chosen via parameters).

• Transform the metric components and there derivatives on the 2-spheres from Cartesian coordinates into a spherical coordinate system.

• Calculate the physical metric on these spheres if a conformal factor is being used.

• Calculate the transformation from the coordinate radius to an areal radius for each sphere.

• Calculate the $$S$$ factor on each sphere. Combined with the areal radius This also produces an estimate of the mass.

• Calculate the six Regge-Wheeler variables, and required radial derivatives, on these spheres by integration of combinations of the metric components over each sphere.

• Contruct the gauge invariant quantities from these Regge-Wheeler variables.

#### 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.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 (l) is used for the first column.

• rsch_R?.[tu]l

The extracted areal radius on each 2-sphere.

• mass_R?.[tu]l

Mass estimate calculated from $$g_{rr}$$ on each 2-sphere.

• Qeven_R?_??.[tu]l

The even parity gauge invariate variable (waveform) on each 2-sphere. This is a complex quantity, the 2nd column is the real part, and the third column the imaginary part.

• Qodd_R?_??.[tu]l

The odd parity gauge invariate variable (waveform) on each 2-sphere. This is a complex quantity, the 2nd column is the real part, and the third column the imaginary part.

• ADMmass_R?.[tu]l

Estimate of ADM mass enclosed within each 2-sphere. (To produce this set doADMmass = ‘‘yes’’).

• momentum_[xyz]_R?.[tu]l

Estimate of momentum at each 2-sphere. (To produce this set do_momentum = ‘‘yes’’).

• spin_[xyz]_R?.[tu]l

Estimate of momentum at each 2-sphere. (To produce this set do_spin = ‘‘yes’’).

### 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(th eta)))” dr_gtt ”drsch_dri=int(dr_gt t)” dr_gpp ”drsch_dri=int(dr_gp p/sin(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) see [1] below weight is one for all points, error is O(1/N) open extended ”weight is 23/12,13/12,4/3,2/4, ...0,27/12, error O(1/N)”

[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) Range Default: aerial radius aerial radius ”calculate invariant aerial radius: r=int(sqrt(gtt*gpp -gtp))” see [1] below ”assume Schwarzschild coordinates: r=int(1/2 (gtt+gpp/sin(theta ))” Schwarzschild gtt assume Schwarzschild coordinates: r=int(gtt) Schwarzschild gpp ”assume Schwarzschild coordinates: r=int(gpp/sin(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

Implements:

waveextractl

Inherits:

grid

admbase

staticconformal

io

sphericalsurface

#### Grid Variables

##### 10.0.1 PRIVATE GROUPS
 Group Names Variable Names Details gridsizes_group compact 0 int_nphi description gridsizes stored as variables int_ntheta dimensions 1 distribution CONSTANT group type ARRAY size 100 tags checkpoint=”no” timelevels 1 variable type INT handles_group compact 0 sum_handle description handles for reduction operators dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT do_nothing_group compact 0 do_nothing description if equal to 1 description then WaveExtract won’t do anything (example: all detectors are out of range) dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT sym_factor_group compact 0 sym_factor description symmmetry factor for integrals (depends on domain description sym_factor=2 for bitant for example) dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type REAL interp_metric_arrays compact 0 gxxi description 2D (theta description phi) Arrays for holding the metric and conformal factor and their first derivatives interpolated onto the extraction coordinate sphere gxyi dimensions 2 gxzi distribution DEFAULT gyyi group type ARRAY gyzi size MAXNTHETA size MAXNPHI gzzi tags checkpoint=”no” psii timelevels 1 dx_gxxi variable type REAL surface_arrays compact 0 interp_x description 2D (theta description phi) grid arrays for points on the sphere interp_y dimensions 2 interp_z distribution DEFAULT psi_ext_deriv group type ARRAY ctheta size MAXNTHETA size MAXNPHI cphi tags checkpoint=”no” sintheta timelevels 1 costheta variable type REAL
 Group Names Variable Names Details surface_integrands compact 0 weights description weights and temporary integrands thetaweights dimensions 2 phiweights distribution DEFAULT int_tmp1 group type ARRAY int_tmp2 size MAXNTHETA size MAXNPHI int_tmp3 tags checkpoint=”no” int_tmp4 timelevels 1 int_tmp5 variable type REAL schwarzschild_mass_radius_group compact 0 dtau_dt description Schwarzschild radius description mass and assorted spherical background pieces sph_grr dimensions 0 sph_gtt distribution CONSTANT sph_dr_gtt group type SCALAR sph_gpp tags checkpoint=”no” rsch2 timelevels 1 rsch variable type REAL schwarzschild_mass_radius_results_group compact 0 Schw_Masses description contains Schwarzschild mass/radius from all detectors Schw_Radii dimensions 1 distribution CONSTANT group type ARRAY size MAXIMUM_DETECTOR_NUMBER tags checkpoint=”no” timelevels 1 variable type REAL moncriefq_results_group compact 0 Qodd_Re_Array description contains Moncrief Qeven description Qodd wave indicators from all detectors Qeven_Re_Array dimensions 3 Qodd_Im_Array distribution CONSTANT Qeven_Im_Array group type ARRAY size MAXIMUM_DETECTOR_NUMBER size L_MODE size 2*M_MODE+1 tags checkpoint=”no” timelevels 1 variable type REAL l_m_modes_info_group compact 0 l_min description Information about the modes used for extraction l_max dimensions 0 l_step distribution CONSTANT m_min group type SCALAR m_max tags checkpoint=”no” m_step timelevels 1 max_det_no_param variable type INT moncriefq compact 0 Qodd_Re description Moncrief Qeven description Qodd wave indicators Qodd_Re description real & imaginary part Qodd_Im dimensions 2 Qeven_Re distribution CONSTANT Qeven_Im group type ARRAY size L_MODE size M_MODE+1 tags checkpoint=”no” timelevels 1 variable type REAL
 Group Names Variable Names Details metric_tmp compact 0 gxx_tmp description temp metric for 3d rotation gxy_tmp dimensions 3 gxz_tmp distribution DEFAULT gyy_tmp group type GF gyz_tmp tags tensortypealias=”DD_sym” checkpoint=”no” prolongation=”none” gzz_tmp timelevels 1 variable type REAL current_detector_group compact 0 current_detector description the index number of the current detector dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT current_detector_radius_group compact 0 current_detector_radius description coordinate radius of the current detector dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type REAL my_out_every_det my_out_every_det compact 0 description output frequency dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 vararray_size 100 variable type INT

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: MoncriefQ metric_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