Extracting Gravitational Waves and Other Quantities from Numerical Spacetimes

Gabrielle Allen

June 21, 2025

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 \[ \gamma _{ij}(t,r,\theta ,\phi )=\sum _{l=0}^\infty \sum _{m=-l}^l \gamma _{ij}^{lm}(t,r) \] and \[ \gamma _{ij}(t,r,\t ,\p ) = \sum _{k=0}^6 p_k(t,r) {\bf V}_k(\t ,\p ) \] 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_{\alpha \beta } = {g}^{background}_{\alpha \beta } + \sum _{l,m} g^{lm,odd}_{\alpha \beta } +\sum _{l,m} g^{lm,even}_{\alpha \beta } \] 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/eNˆ )
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 SPHERICALSURFACEINT

nghoststheta
 Scope: shared from SPHERICALSURFACEINT

nsurfaces
 Scope: shared from SPHERICALSURFACEINT

sfpar_maxnphi
 Scope: shared from SPHERICALSURFACEINT

sfpar_maxntheta
 Scope: shared from SPHERICALSURFACEINT

symmetric_x
 Scope: shared from SPHERICALSURFACEBOOLEAN

symmetric_y
 Scope: shared from SPHERICALSURFACEBOOLEAN

symmetric_z
 Scope: shared from SPHERICALSURFACEBOOLEAN

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 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