Spectral characteristic evolution: a new algorithm for gravitational wave propagation

As an international network of gravitational wave observatories come online, the race to the first direct detection of gravitational waves is expected to herald the beginning of gravitational wave astronomy. Detectors such as Advanced LIGO, VIRGO, GEO, and KAGRA aim for strain sensitivities approaching 10⁻²⁴ [1–4]. Nevertheless, signal candidates from compact binaries or supernovae will be on the cusp of detectability, with very poor signal to noise ratios requiring matched filtering for detection [5]. Matched filtering requires a comprehensive template bank, the generation of which has been a primary goal of the field of numerical relativity. These templates cover a range of expected astronomical phenomena, and are generated by a variety of numerical codes, including the Spectral Einstein Code (SpEC) [6]. Filling out the template bank requires a balance of numerical relativity and analytic waveforms (post-Newtonian [7] and effective-one-body [8]), with waveform models often calibrated using numerical results [9–11].

One technical challenge facing the construction of a large template bank is extraction of gauge invariant waveforms from simulations. In general, large computationally intensive simulations are needed to describe the physics of events such as supernovae and compact object binary inspirals. While a waveform-like signal can readily be extracted from anywhere in the computational domain, waveforms are only rigorously defined at future null infinity. Finite radii waveform approximations are universally contaminated by coordinate system dynamics, or gauge effects, which are poorly understood and nearly impossible to remove or even quantify. Comparison with Cauchy characteristic extraction, or CCE, an alternative which applies characteristic evolution to enable waveform computation at future null infinity, suggests that extrapolation gauge errors could dominate the global error [12]. Characteristic evolution has been previously implemented at up to 4th order radial accuracy [13], while complete extraction has been achieved with finite difference/volume methods up to 2nd order [14, 15].

Here we implement inner boundary extraction and evolution. This must be combined with an appropriate outer boundary algorithm to generate the gauge-invariant news at future null infinity. In a hypothetical complete system illustrated in figure 1(d), Cauchy characteristic matching allows inflowing energy to cross the inner boundary (injection) and enter the Cauchy evolution, ameliorating issues with boundary reflections. Matching, however, faces many technical challenges beyond the scope of this paper, and has only been implemented in the linearized limit [16, 17].

Zoom In Zoom Out Reset image size

Another method is extrapolation [12], which takes metric quantities at a series of increasing radii, then fits them to a function in $l=1/r$. In this coordinate, l = 0 corresponds to infinitely far from the origin, where the extrapolated metric quantities are converted to waveforms. Because each sampling point is subject to an unknown degree of coordinate effect contamination, the extrapolated waveform is itself not gauge invariant.

While characteristic evolution is computationally and conceptually much more involved than extrapolation, it is able to provide gauge invariant waveforms unaffected by coordinate effects. These waveforms are unique modulo the supertranslations, a four parameter subgroup of the BMS group, corresponding to arbitrary inertial observer initial conditions at future null infinity [18, 19]. Removing gauge effects through characteristic evolution is essential for obtaining accurate and useful waveforms. Essentially, characteristic evolution takes metric data at a topologically spherical worldtube Γ enclosing the relativistic Cauchy evolution and evolves it, as well as initial data, outward to future null infinity, or ${{\mathcal{I}}^{+}}$. At the outer boundary, the metric quantities can be read off and, in combination with an inertial conformal coordinate system, used to calculate the true waveform.

The calculation is performed in the Bondi system, in which radial coordinates are outgoing null rays: normal to the worldtube and to each time slice. The spherical coordinates and time-like foliation is adapted from the Cauchy evolution via the worldtube boundary data, illustrated in figure 1(c).

While this simplifies aspects of the evolution, the domain is still infinite, and must be compactified. We use the compactification $r=R\rho /(1-\rho )$, where r is the Bondi radius, ρ is the compactified null radial coordinate and R the compactification parameter. Setting R to the Bondi radius of the worldtube $R={{r}_{|\Gamma }}$ confines $\rho \in [1/2,1]$, as shown in figure 2(c). While the worldtube is defined by a constant coordinate radius sphere, the Bondi radius at this surface is a variable function of time and angles, giving rise to a wobbly, non-spherical shape in the Bondi coordinate system. This variable compactification parameter requires additional terms in the equations.

Zoom In Zoom Out Reset image size

Fixing $\rho \in [1/2,1]$ and having a variable compactification parameter is different to the approach used in earlier incarnations of CCE, which either interpolated the boundaryʼs variable position or used a different compactification scheme [12, 20]. While the Pitt null code [14] assumed that a coordinate sphere of constant Cartesian radius formed the worldtube, we lift that constraint here. A fixed computational domain for characteristic evolution enables conceptually simple radial integration and dynamically variable extraction radii.

Figure 3 shows a series of diagrams illustrating the derivation of the computational domain.

Zoom In Zoom Out Reset image size

An implementation of characteristic evolution was developed by the group of Winicour during the mid 1990s [14, 16, 20], and is now part of the publicly available Einstein tool kit (PITTNullCode) [21]. In its original form it uses finite differences, achieving 2nd order accuracy. The code has been updated and adapted many times since, and in its current form takes $\mathcal{O}$(days) to produce a waveform at an accuracy that matches that of the Cauchy evolution for $\mathcal{O}(1000\;M)$ SpEC runs. Five years ago an $\mathcal{O}(1000\;M)$ run was sufficiently challenging that there was no CCE bottleneck. Today, with hundreds of runs exceeding dozens of orbits and a few exceeding ${{10}^{5}}\;M$, a faster algorithm is needed. The goal is to have the SpEC characteristic code run alongside the Cauchy evolution at an insigificant additional cost.

This algorithm would build on the formalism of the Pitt null code, but exploit SpEC capabilities to use spectral methods and a much coarser grid to achieve equivalent accuracy. The finite difference algorithm marches outwards one null parallelogram at a time, shown in figure 1(b). A spectral algorithm would rapidly calculate high-accuracy radial data in a single step of numerical integration, as shown in figure 1(c).

The difficulty of a spectral approach lies in consistently treating divergences. As written, the hypersurface evolution equations' source terms are linear or greater order in r and diverge at infinity. In a finite domain, large but finite terms on either side would numerically cancel, leaving a valid result with no further complications. The Pitt null code uses an asymptotic form to solve the final step to ${{\mathcal{I}}^{+}}$. To include a point at infinity within a single domain spectral scheme, however, this divergence has to be understood and pre-emptively cancelled. Here, we present a novel approach to regularizing the full nonlinear system, enabling a fully spectral treatment.

3.1. General information

The Bondi metric can be expressed as

$$\begin{eqnarray}\begin{array}{rcl} {\rm d}{{s}^{2}} & = & -\left( {{{\rm e}}^{2\beta }}(rW+1)-{{r}^{2}}{{h}_{AB}}{{U}^{A}}{{U}^{B}} \right){\rm d}{{u}^{2}} \\ {} & {} & -2{{{\rm e}}^{2\beta }}{\rm d}u{\rm d}r-2{{r}^{2}}{{h}_{AB}}{{U}^{B}}{\rm d}u{\rm d}{{y}^{A}}+{{r}^{2}}{{h}_{AB}}{\rm d}{{y}^{A}}{\rm d}{{y}^{B}}. \\ \end{array}\end{eqnarray} \tag{ 3.1 }$$

y^A are angular coordinates, where the uppercase Latin indices ^A, ^B, and so on range from 2 to 3. The quantities $W,{{h}_{AB}},{{U}^{A}},{{Q}_{A}},\beta$ parametrize the metric, representing respectively the mass aspect, the conformal two-metric, the shift and its radial derivative, and the lapse. Where non-scalar, they are related to complex spin-weighted quantities by contracting them with the appropriate dyad. The dyad q^A is a complex field on the unit sphere satisfying ${{q}^{A}}{{q}_{A}}=0$, ${{q}^{A}}{{\bar{q}}_{A}}=2$, ${{q}^{A}}={{q}^{AB}}{{q}_{B}}$, with ${{q}^{AB}}{{q}_{BC}}=\delta _{C}^{A}$ and ${{q}_{AB}}=({{q}_{A}}{{\bar{q}}_{B}}+{{\bar{q}}_{A}}{{q}_{B}})/2$, the unit sphere metric. Under this convention, the spin-weighted functions $U={{U}^{A}}{{q}_{A}}$ and $Q={{Q}_{A}}{{q}^{A}},$ while $J={{h}_{AB}}{{q}^{A}}{{q}^{B}}/2$ uniquely determines the spherical conformal two-metric component of the general four-metric [20]. We chose a dyad consistent with our formulation of the eth operator ð [22], given by ${{q}^{A}}=(-1,-{\rm i}/{\rm sin} \;\theta )$. This is regular everywhere except the poles, which we can avoid through careful choice of grid points. It is worth noting that any choice of angular coordinates is possible. Other conventions use multiple patches to avoid singularities at the poles, in which the phase dislocation due to spin-weight when moving from patch to patch is explicit.

The key to the Bondi formulation is that all the spin-weighted metric quantities have a heirarchical structure, enabling their natural ordering as a nested series of self-referential equations on the outgoing null hypersurface. Their derivation is discussed in Bishop et al [16, 20]. In a relatively simple form they are

$$\begin{eqnarray}&&\beta {{,}_{r}}={{N}_{\beta }},\end{eqnarray} \tag{ 3.2 }$$

$$\begin{eqnarray}&&({{r}^{2}}Q){{,}_{r}}=-{{r}^{2}}(\bar{\eth }J+\eth K){{,}_{r}}+2{{r}^{4}}\eth ({{r}^{-2}}\beta ){{,}_{r}}+{{N}_{Q}},\end{eqnarray} \tag{ 3.3 }$$

$$\begin{eqnarray}&&U{{,}_{r}}={{r}^{-2}}{{{\rm e}}^{2\beta }}Q+{{N}_{U}},\end{eqnarray} \tag{ 3.4 }$$

$$\begin{eqnarray}&&({{r}^{2}}W){{,}_{r}}=\frac{1}{2}{{{\rm e}}^{2\beta }}\mathcal{R}-1-{{{\rm e}}^{\beta }}\eth \bar{\eth }{{{\rm e}}^{\beta }}+\frac{1}{4}{{r}^{-2}}({{r}^{4}}(\eth \bar{U}+\bar{\eth }U)){{,}_{r}}+{{N}_{W}},\end{eqnarray} \tag{ 3.5 }$$

$$\begin{eqnarray}&&{\rm where}\;{\rm the}\;{\rm scalar}\;{\rm curvature}\end{eqnarray} \tag{ 3.6 }$$

$$\begin{eqnarray}&&\mathcal{R}=2\;K-\eth \bar{\eth }K+\frac{1}{2}({{\bar{\eth }}^{2}}\;J+{{\eth }^{2}}\bar{J})+\frac{1}{4\;K}(\bar{\eth }\bar{J}\eth J-\bar{\eth }J\eth \bar{J}),\end{eqnarray} \tag{ 3.7 }$$

$$\begin{eqnarray}&&{\rm and}\;{\rm the}\;{\rm time}\;{\rm derivative}\;{\rm term}\end{eqnarray} \tag{ 3.8 }$$

$$\begin{eqnarray}&&2(rJ){{,}_{ur}}=\left( (1+rW)(rJ){{,}_{r}} \right){{,}_{r}}-{{r}^{-1}}({{r}^{2}}\eth U){{,}_{r}}+2{{r}^{-1}}{{{\rm e}}^{\beta }}{{\eth }^{2}}{{{\rm e}}^{\beta }}-(rW){{,}_{r}}J+{{N}_{J}},\end{eqnarray} \tag{ 3.9 }$$

$$\begin{eqnarray}&&{\rm where}\end{eqnarray} \tag{ 3.10 }$$

$$\begin{eqnarray}&&1={{K}^{2}}-J\bar{J}.\end{eqnarray} \tag{ 3.11 }$$

The nonlinear terms ${{N}_{\beta }},{{N}_{Q}},{{N}_{U}},{{N}_{W}},{{N}_{J}}$ are given in appendix D ¹ . The radial compactification is given by $r=R\rho /(1-\rho )$, where the compactification parameter $R(u,\;\theta ,\phi )$ is the (not necessarily constant) Bondi areal radius of the worldtube.

On each hypersurface in the spacetime foliation, each equation is solved in turn. Given J, β is solved, then U, Q, and W in turn, enabling the computation of $J{{,}_{u}}$. $J{{,}_{u}}$ permits a step forward in time and J is thus defined on the next hypersurface.

Spherical derivatives are implemented using the eth operator ð on spin-weighted spherical harmonics. ð is given by the contraction of the dyad with the derivative operator [22].

In spherical coordinates, this can be written

$$\begin{eqnarray}&&\eth \eta =-({{{\rm sin} }^{s}}\;\theta )\left( \frac{\partial }{\partial \;\theta }+\frac{{\rm i}}{{\rm sin} \;\theta }\frac{\partial }{\partial \phi } \right)({{{\rm sin} }^{-s}}\;\theta \;\eta ),\end{eqnarray} \tag{ 3.12 }$$

$$\begin{eqnarray}&&\bar{\eth }\eta =-({{{\rm sin} }^{-s}}\;\theta )\left( \frac{\partial }{\partial \;\theta }-\frac{{\rm i}}{{\rm sin} \;\theta }\frac{\partial }{\partial \phi } \right)({{{\rm sin} }^{s}}\;\theta \;\eta ).\end{eqnarray} \tag{ 3.13 }$$

$\eth$ and $\bar{\eth }$ increment and decrement respectively the spin-weight s of the quantity they act upon. Details are given in appendix C.

In the Pitt null code approach to CCE, a finite difference-based algorithm is used to solve the hypersurface equations, radially marching from the inner boundary outward, one point at a time. In our algorithm, we use spectral methods to calculate the radial indefinite integral on all collocation points along a spherical radial spoke in a single calculation. A spectral approach is faster and more accurate, but requires the finessing of a few technical difficulties.

The most obvious of these appears in the hypersurface equations for W and Q. To calculate a numerical integral, it is necessary to express the integrand as a bounded function on all radial collocation points, including those at ${{\mathcal{I}}^{+}}$. Given that Q is regular and $\mathcal{O}(1)$ at ${{\mathcal{I}}^{+}},{{r}^{2}}Q$ is clearly divergent. The integral of the right-hand side is similarly divergent for the outermost collocation point, where $r\to \infty$ or $\rho \to 1$. Here, we solve this problem by expressing the right-hand side as a Laurent series around the pole at ρ = 1 and then repeatedly integrating by parts. Mathematically, if the equation is written

$$\begin{eqnarray}&&{{({{r}^{2}}Q)}_{,\rho }}=A+\frac{B}{1-\rho }+\frac{C}{{{(1-\rho )}^{2}}}+\frac{D}{{{(1-\rho )}^{3}}},\end{eqnarray} \tag{ 3.14 }$$

then, as shown in appendix A

$$\begin{eqnarray}\begin{array}{rcl} Q & = & {{Q}_{0}}+\frac{D}{2{{\rho }^{2}}{{R}^{2}}}-(-C+D^{\prime} /2)\frac{1-\rho }{{{\rho }^{2}}{{R}^{2}}}-(B-C^{\prime} +D^{\prime\prime} /2)\frac{{{(1-\rho )}^{2}}{\rm log} (1-\rho )}{{{\rho }^{2}}{{R}^{2}}} \\ {} & {} & -(B^{\prime} -C^{\prime\prime} +D^{\prime\prime\prime} /2)\frac{(\rho +(1-\rho ){\rm log} (1-\rho )){{(1-\rho )}^{2}}}{{{\rho }^{2}}{{R}^{2}}} \\ {} & {} & +\frac{{{(1-\rho )}^{2}}}{{{\rho }^{2}}{{R}^{2}}}\int A+(B^{\prime\prime} -C^{\prime\prime\prime} +D/2)(\rho +(1-\rho ){\rm log} (1-\rho )), \\ \end{array}\end{eqnarray} \tag{ 3.15 }$$

where $D^{\prime} =\frac{\partial D}{\partial \rho }$ and so on. In this context, the $\int$ operator refers to a radial numerical integral where the integration constant is fixed such that the inner boundary value (at the worldtube) vanishes. Note that all terms within the integral are bounded within the domain $\rho \in [1/2,1]$, including at ${{\mathcal{I}}^{+}}$. On the outer boundary, terms in ${\rm log} (1-\rho )$ must cancel out to preserve ${{C}^{\infty }}$ differentiability.

A second, less obvious, issue is that the right-hand side of the $J{{,}_{u\rho }}$ hypersurface equation has nonlinear terms with the desired quantity $J{{,}_{u}}$ in them (as seen in the equation for P₁, equation (A.6), Bishop et al [16]). In order to perform the radial integral and solve for $J{{,}_{u}}$ in a single step, these nonlinear terms have to be somehow removed. One approach is to factorize using an integrating factor and, conceptually, that is what is done. We shall illustrate this first with a simple example.

Given

$$\begin{eqnarray}&&J{{,}_{u\rho }}+AJ{{,}_{u}}=B,\;{\rm then}\end{eqnarray} \tag{ 3.16 }$$

$$\begin{eqnarray}&&\left( {{{\rm e}}^{\int A}}J{{,}_{u}} \right){{,}_{\rho }}={{{\rm e}}^{\int A}}B,\;{\rm and}\end{eqnarray} \tag{ 3.17 }$$

$$\begin{eqnarray}&&{{J}_{,u}}={{{\rm e}}^{-\int A}}\int {{{\rm e}}^{\int A}}B.\end{eqnarray} \tag{ 3.18 }$$

The actual equations are, however, more difficult as the variable is complex, and the nonlinear term also includes the complex conjugate $\bar{J}{{,}_{u}}$ term, which makes a simple factorization impossible. Writing $J{{,}_{u}}=\Phi$, the actual equation can be written

$$\begin{eqnarray}&&{{(r\Phi )}_{,\rho }}-(rJ)(\Phi \bar{\Gamma }+\bar{\Phi }\Gamma )=A+\frac{B}{1-\rho }+\frac{C}{{{(1-\rho )}^{2}}},\;{\rm where}\end{eqnarray} \tag{ 3.19 }$$

$$\begin{eqnarray}&&\Gamma =\left( {{J}_{,\rho }}-J\frac{{{K}_{,\rho }}}{K} \right).\end{eqnarray} \tag{ 3.20 }$$

The right-hand side is a Laurent expansion analogous to the equations for Q and W. The part $(\Phi \bar{\Gamma }+\bar{\Phi }\Gamma )$ is a quantity plus its conjugate and is thus wholly real. This leads to the insight that a real–imaginary split of the equation is productive. Writing $J={{J}_{\mathcal{R}}}+{\rm i}{{J}_{\mathcal{I}}}$, we can exploit the isomorphism between complex numbers and non-singular 2 × 2 matrices by writing

$$\begin{eqnarray}&&\left( \begin{array}{ccccccccccccccc} r{{\Phi }_{\mathcal{R}}} \\ r{{\Phi }_{\mathcal{I}}} \\ \end{array} \right),\rho -\left( \begin{array}{ccccccccccccccc} {{J}_{\mathcal{R}}}{{\Gamma }_{\mathcal{R}}} & {{J}_{\mathcal{R}}}{{\Gamma }_{\mathcal{I}}} \\ {{J}_{\mathcal{I}}}{{\Gamma }_{\mathcal{R}}} & {{J}_{\mathcal{I}}}{{\Gamma }_{\mathcal{I}}} \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} r{{\Phi }_{\mathcal{R}}} \\ r{{\Phi }_{\mathcal{I}}} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} {\rm rh}{{{\rm s}}_{\mathcal{R}}} \\ {\rm rh}{{{\rm s}}_{\mathcal{I}}} \\ \end{array} \right).\end{eqnarray} \tag{ 3.21 }$$

We have restored the integrating factor form of the equation, only now in matrix form. The use of non-constant matrices for an integrating factor requires the calculation of commutators, which in all but the trivial case are non-zero. The required formalism is the Magnus expansion [23, 24], in which the usual integrating factor is supplemented by integrals of progressively higher order commutators.

Let

$$\begin{eqnarray}F=\left( \begin{array}{ccccccccccccccc} {{J}_{\mathcal{R}}}{{\Gamma }_{\mathcal{R}}} & {{J}_{\mathcal{R}}}{{\Gamma }_{\mathcal{I}}} \\ {{J}_{\mathcal{I}}}{{\Gamma }_{\mathcal{R}}} & {{J}_{\mathcal{I}}}{{\Gamma }_{\mathcal{I}}} \\ \end{array} \right),\end{eqnarray} \tag{ 3.22 }$$

and

$$\begin{eqnarray}&&{\boldsymbol{ \Phi }} =\left( \begin{array}{ccccccccccccccc} {{\Phi }_{\mathcal{R}}} \\ {{\Phi }_{\mathcal{I}}} \\ \end{array} \right).\end{eqnarray} \tag{ 3.23 }$$

Then

$$\begin{eqnarray}&&{{(r{\boldsymbol{ \Phi }} )}_{,\rho }}-F.(r{\boldsymbol{ \Phi }} )={\rm exp} \left( \Omega (\rho ) \right).{{\left( {\rm exp} \left( -\Omega (\rho ) \right).(r\Phi ) \right)}_{,\rho }},\end{eqnarray} \tag{ 3.24 }$$

where

$$\begin{eqnarray}&&\Omega (\rho )=\mathop{\sum }\limits_{k=1}^{\infty }{{\Omega }_{k}}(\rho ).\end{eqnarray} \tag{ 3.25 }$$

Write $F({{\rho }_{1}})={{F}_{1}}$. Then ${{\Omega }_{k}}$ forms a series called the Magnus expansion

$$\begin{eqnarray}&&{{\Omega }_{1}}(\rho )=\int _{{{\rho }_{0}}}^{\rho }{\rm d}{{\rho }_{1}}{{F}_{1}},\end{eqnarray} \tag{ 3.26 }$$

$$\begin{eqnarray}&&{{\Omega }_{2}}(\rho )=\frac{1}{2}\int _{{{\rho }_{0}}}^{\rho }{\rm d}{{\rho }_{1}}\int _{{{\rho }_{0}}}^{\rho }{\rm d}{{\rho }_{2}}[{{F}_{1}},{{F}_{2}}],\end{eqnarray} \tag{ 3.27 }$$

$$\begin{eqnarray}&&{{\Omega }_{3}}(\rho )=\frac{1}{6}\int _{{{\rho }_{0}}}^{\rho }{\rm d}{{\rho }_{1}}\int _{{{\rho }_{0}}}^{\rho }{\rm d}{{\rho }_{2}}\int _{{{\rho }_{0}}}^{\rho }{\rm d}{{\rho }_{3}}[{{F}_{1}},[{{F}_{2}},{{F}_{3}}]]+[{{F}_{3}}[{{F}_{2}},{{F}_{1}}]],\end{eqnarray} \tag{ 3.28 }$$

and so on.

In practice, this series must be truncated. Fortunately, $|F|\approx {{10}^{-3}}\lt \lt 1$ in practical cases, so the series converges rapidly, necessitating calculation of only the first three terms. Although a formalism exists [19] to deal with the nonlinear terms without resorting to the Magnus expansion, it requires transformation in terms of dyads, whose issues around the poles we have been careful to avoid.

With the Magnus expansion, the troublesome nonlinear terms are readily dealt with and the equation can be expressed in the form

$$\begin{eqnarray}&&(r\Phi ){{,}_{\rho }}=\tilde{A}+\frac{\tilde{B}}{1-\rho }+\frac{\tilde{C}}{{{(1-\rho )}^{2}}}.\end{eqnarray} \tag{ 3.29 }$$

This is solved analogously to the radial hypersurface equations for W and Q.

3.2. Details specific to our implementation

The spectral characteristic evolution algorithm is analytically derived, but how does it actually perform? We tested the stability and convergence of the code under a variety of circumstances designed to far exceed the demands of any actual CCE run [12].

For a stability test, we ran the algorithm with a variety of settings for a million steps with white noise initial and boundary conditions of magnitude 10⁻⁶. The linear setting truncates all nonlinear terms, and represents a baseline condition. The nonlinear setting restores all nonlinear terms in the equations. The most general setting includes a variable inner boundary position, encoded in the magnitude of the compactification parameter R. Each of these three conditions was run, and in all three cases, the norm of J was stable. In particular, all three runs do not grow exponentially, as seen in figure 4.

Zoom In Zoom Out Reset image size

The purpose of a spectral algorithm is to obtain faster convergence, particularly with respect to radial integration. We ran a series of tests while varying ${{N}_{\rho }}$ between 6 and 46. The test was run on the generic run discussed in section 6 between $t=1000\;M$ and $t=1100\;M$, and the results averaged between $t=1050\;M$ and $t=1100\;M$ to remove any transient contamination. While our algorithm can use any one-dimensional spectral representation for the radial direction, these tests were conducted using Gauss–Chebyshev–Lobatto polynomial basis functions, as discussed in the previous section. Figure 5 shows that local relative error in J converges exponentially as the number of radial points increases up to around 24, at which point convergence becomes sub-exponential due to roundoff error introduced in the integration algorithm.

Zoom In Zoom Out Reset image size

Convergence with angular resolution was calculated identically to that for figure 5, and shows rapid exponential convergence, as seen in figure 6. Our implementation uses the Spinsfast package [27] for spin-weighted spherical harmonic computation in a manner analogous to Spherepack.

Zoom In Zoom Out Reset image size

Finally, timestep convergence was analyzed. For the purposes of this test, the minimum grid spacing to timestep was held constant, while the timestep was varied over more than an order of magnitude. The code displays 4th order time convergence, consistent with the chosen integrator RK4, as shown in figure 7.

Zoom In Zoom Out Reset image size

For longer run comparisons we used the generic precessing binary black hole simulation detailed in table 1 (case 4) of Taylor et al [12]. Its parameters are mass ratio q = 3, black hole spin ${{\chi }_{1}}=(0.7,0,0.7)/\sqrt{2}$ and ${{\chi }_{2}}=(-0.3,0,0.3)/\sqrt{2}$, number of orbits 26, total time $T=7509\;M$, initial eccentricity 10⁻³, initial frequency ${{\omega }_{{\rm ini}}}=0.032/M$, and extraction (coordinate) radius $R=100\;M$. We performed three runs using the Pitt null code for baseline comparisons with three runs using the new spectral characteristic evolution code, using parameters shown in table 1. The spectral code converges rapidly to within the error implied by the Pitt null code, as shown in figure 8. Parameters for the SpEC runs were specifically chosen for comparible levels of error with the Pitt null code. Note that despite SpEC code spatial resolution parameters being chosen for consistent minimum grid spacing to time step ratio, rather than equal numbers of points, the global error comparison is not adversely affected. Both codes were run on the same cluster with dense output.

Zoom In Zoom Out Reset image size

Table 1.  Parameters used for code comparisons. The Pitt null code uses two stereographic patches of size N²_stereo. For both codes, the total number of angular points is given by $2N_{{\rm stereo}}^{2}$ and $2{{L}^{2}}$ respectively.

Run	Pitt 1	Pitt 2	Pitt 3	SpEC 1	SpEC 2	SpEC 3
Nr	100	150	200	10	12	14
Nstereo or L	40	60	80	12	14	17
$\Delta t/M$	0.1	$0.0666...$	0.05	1.0	$0.666...$	0.5
T (CPU hours)	2688	5760	6912	12	31	52

Note also that the resolution of the SpEC runs is an order of magnitude lower, the time steps an order of magnitude longer, and the run time more than two orders of magnitude faster for equivalent or superior accuracy. In this case, accuracy is ultimately limited by the order of the time stepper (RK4) and the order of the time interpolation on the inner boundary (also 4th order).

Global accuracy and convergence was assessed by graphing the relative complex difference between runs of adjacent accuracy. The errors are computed according to

$$\begin{eqnarray*}&&{{E}_{{\rm Pitt}\;{\rm low}\;{\rm res}}}=\frac{|{{j}_{{{22}_{{\rm Pitt}\ 1}}}}-{{j}_{{{22}_{{\rm Pitt}\ 2}}}}|}{|{{j}_{{{22}_{{\rm Pitt}\ 2}}}}|},{{E}_{{\rm Pitt}\;{\rm high}\;{\rm res}}}=\frac{|{{j}_{{{22}_{{\rm Pitt}\ 2}}}}-{{j}_{{{22}_{{\rm Pitt}\ 3}}}}|}{|{{j}_{{{22}_{{\rm Pitt}\ 3}}}}|},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{{E}_{{\rm SpEC}\;{\rm low}\;{\rm res}}}=\frac{|{{j}_{{{22}_{{\rm SpEC}\ 1}}}}-{{j}_{{{22}_{{\rm SpEC}\ 2}}}}|}{|{{j}_{{{22}_{{\rm SpEC}\ 2}}}}|},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{{E}_{{\rm SpEC}\;{\rm high}\;{\rm res}}}=\frac{|{{j}_{{{22}_{{\rm SpEC}\ 2}}}}-{{j}_{{{22}_{{\rm SpEC}\ 3}}}}|}{|{{j}_{{{22}_{{\rm SpEC}\ 3}}}}|},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{{E}_{{\rm Pitt}\;{\rm versus}\;{\rm SpEC}}}=\frac{|{{j}_{{{22}_{{\rm SpEC}\ 3}}}}-{{j}_{{{22}_{{\rm Pitt}\ 3}}}}|}{|{{j}_{{{22}_{{\rm SpEC}\ 3}}}}|},\end{eqnarray*}$$

where j₂₂ is the ₂₂ spherical harmonic coefficient of J.

Figure 8 shows strong and consistent convergence of the SpEC runs over the full $7509\;M,$ as well as consistent agreement between the SpEC and Pitt runs. Figure 9 shows convergence during the junk radiation part of the run, where partial loss of agreement between the Pitt and SpEC runs is caused by their respective differences in unphysical junk radiation resolution. The error nevertheless remains well bounded for the entire run.

Zoom In Zoom Out Reset image size

A new algorithm for spectral characteristic evolution has been developed, implemented, and demonstrated within the SpEC framework. It exploits analytic regularization of the vacuum hypersurface equations and the accuracy and speed of spectral methods. Stability, self-convergence, and convergent agreement with the existing finite difference characteristic evolution code are demonstrated. This algorithm will form the basis for a complete extraction and matching methodology that will enable gauge invariant waveforms and junk radiation-free initial conditions to be computed on the fly.

We thank Jeffrey Winicour for his invaluable insight and encyclopædic knowledge of all things extraction. We thank Nicholas Taylor for his generic spin BBH run we used to test and baseline code performance. We thank Mark Scheel, Yanbei Chen, and Christian Reisswig for their advice, support, and technical expertise. This research used the Spectral Einstein Code (SpEC) [6]. The Caltech cluster zwicky.cacr.caltech.edu is an essential resource for SpEC related research, supported by the Sherman Fairchild Foundation and by NSF award PHY-0960291. This research also used the Extreme Science and Engineering Discovery Environment (XSEDE) under grant TG-PHY990002. The UCSD cluster ccom-boom.ucsd.edu was used during code development. This project was supported by the Sherman Fairchild Foundation, and by NSF Grants PHY-1068881, AST-1333520, and CAREER Grant PHY-0956189.

This appendix describes the process of future null infinity regularization using integration by parts. A Laurent series is an expansion about a pole of (in this case) finite order. It is the logical extension of a Taylor series to functions that diverge in well defined ways.

We wish to radially integrate the following equation

$$\begin{eqnarray}&&{{({{r}^{2}}Q)}_{,\rho }}=A+\frac{B}{1-\rho }+\frac{C}{{{(1-\rho )}^{2}}}+\frac{D}{{{(1-\rho )}^{3}}}.\end{eqnarray} \tag{ A.1 }$$

Note that $r=R\rho /(1-\rho )$, so both sides are infinite at ρ = 1. The most divergent term is the D term—we integrate by parts

$$\begin{eqnarray}&&\int \frac{D}{{{(1-\rho )}^{3}}}{\rm d}\rho =\frac{D/2}{{{(1-\rho )}^{2}}}-\int \frac{D^{\prime} /2}{{{(1-\rho )}^{2}}}{\rm d}\rho .\end{eqnarray} \tag{ A.2 }$$

In practice, limits of integration are $\rho \in [1/2,1]$ where $\rho =1/2$ is the worldtube inner boundary of the domain.

This term is now an integral of the same order as the term in C, as well as a term external to the integral that is of the same order. This is cancelled out through division by ${{r}^{2}}={{R}^{2}}{{\rho }^{2}}/{{(1-\rho )}^{2}}$. In this way, all the terms are regularized, i.e. finite throughout the domain. We include the $D^{\prime} /2$ term in the C integral, and so on

$$\begin{eqnarray}\begin{array}{rcl} \int \frac{-C+D^{\prime} /2}{{{(1-\rho )}^{2}}}{\rm d}\rho & = & \int \frac{\partial }{\partial \rho }\left( \frac{1}{1-\rho } \right)(-C+D^{\prime} /2){\rm d}\rho \\ {} & = & \frac{-C+D^{\prime} /2}{1-\rho }-\int \frac{-C^{\prime} +D^{\prime\prime} /2}{1-\rho }{\rm d}\rho , \\ \end{array}\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&\begin{array}{l} \int \frac{B-C^{\prime} +D^{\prime\prime} /2}{1-\rho }{\rm d}\rho \\ \quad =-(B-C^{\prime} +D^{\prime\prime} /2){\rm log} (1-\rho )+\int (B^{\prime} -C^{\prime\prime} +D^{\prime\prime\prime} /2){\rm log} (1-\rho ){\rm d}\rho , \\ \end{array}\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&\begin{array}{l} \int (B^{\prime} -C^{\prime\prime} +D^{\prime\prime\prime} /2){\rm log} (1-\rho ){\rm d}\rho = \\ \quad -(B^{\prime} -C^{\prime\prime} +D^{\prime\prime\prime} /2)(\rho +(1-\rho ){\rm log} (1-\rho )) \\ \quad +\int (B^{\prime\prime} -C^{\prime\prime\prime} +D/2)(\rho +(1-\rho ){\rm log} (1-\rho )){\rm d}\rho . \\ \end{array}\end{eqnarray} \tag{ A.5 }$$

Crucially, the integral term is now bounded in the domain. This means it can be computed numerically. Computationally, the limits of integration are enforced by subtracting the value of the function on the inner boundary. Q₀ is the boundary value of Q. Combining terms, equation (3.11) can be expressed:

$$\begin{eqnarray}\begin{array}{rcl} Q & = & {{Q}_{0}}+\frac{D}{2{{\rho }^{2}}{{R}^{2}}}-(-C+D^{\prime} /2)\frac{1-\rho }{{{\rho }^{2}}{{R}^{2}}}-(B-C^{\prime} +D^{\prime\prime} /2)\frac{{{(1-\rho )}^{2}}{\rm log} (1-\rho )}{{{\rho }^{2}}{{R}^{2}}} \\ {} & {} & -(B^{\prime} -C^{\prime\prime} +D^{\prime\prime\prime} /2)\frac{(\rho +(1-\rho ){\rm log} (1-\rho )){{(1-\rho )}^{2}}}{{{\rho }^{2}}{{R}^{2}}} \\ {} & {} & +\frac{{{(1-\rho )}^{2}}}{{{\rho }^{2}}{{R}^{2}}}\int A+(B^{\prime\prime} -C^{\prime\prime\prime} +D/2)(\rho +(1-\rho ){\rm log} (1-\rho )). \\ \end{array}\end{eqnarray} \tag{ A.6 }$$

The boundary algorithm formalism is drawn from Bishop et al [16]. It begins with metric quantities on the worldtube forming the boundary between the space-like foliated numerical GR simulation and the null-foliated CCE domain. It calculates several intermediate helper quantities to simplify the computational complexity. Finally, it produces boundary values for each of the hypersurface or Bondi metric quantities. Note that all quantities in this section refer to inner boundary values only.

B.1. Initial metric quantities

The metric quantities are extracted directly from the simulation. They are the covariant three-metric g_ij, the contravariant three-metric ${{g}^{ij}}={{({{g}_{ij}})}^{-1}}$, the co- and contravariant three-metric derivatives ${{g}_{ij}}{{,}_{\alpha }}$ and ${{g}^{ij}}{{,}_{\alpha }}$ (calculated using ${{g}^{ij}}{{,}_{\gamma }}=-{{g}^{ik}}{{g}^{jl}}{{g}_{kl}}{{,}_{\gamma }}$), and their four-metric counterparts. The time components are calculated from the lapse $\alpha =\sqrt{{{g}_{i0}}{{\beta }^{i}}-{{g}_{00}}}$ and shift ${{\beta }^{i}}={{g}^{ij}}{{g}_{j0}}$ as follows

$$\begin{eqnarray}&&{{g}_{k0}}={{g}_{0k}}={{g}_{ki}}{{\beta }^{i}},\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray}&&{{g}_{00}}={{g}_{i0}}{{\beta }^{i}}-{{\alpha }^{2}}={{g}_{ik}}{{\beta }^{i}}{{\beta }^{k}}-{{\alpha }^{2}},\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&{{g}_{k0}}{{,}_{\gamma }}={{g}_{0k}}{{,}_{\gamma }}=({{g}_{ki}}{{\beta }^{i}}){{,}_{\gamma }}={{g}_{ki}}{{,}_{\gamma }}{{\beta }^{i}}+{{g}_{ki}}{{\beta }^{i}}{{,}_{\gamma }},\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&{{g}_{00}}{{,}_{\gamma }}={{g}_{i0}}{{,}_{\gamma }}{{\beta }^{i}}+{{g}_{i0}}{{\beta }^{i}}{{,}_{\gamma }}-2\alpha \alpha {{,}_{\gamma }},\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&{{g}^{\alpha \beta }}{{,}_{\gamma }}=-{{g}^{\alpha \delta }}{{g}^{\beta \epsilon }}{{g}_{\delta \epsilon }}{{,}_{\gamma }}.\end{eqnarray} \tag{ B.5 }$$

B.2. Intermediate quantities

Additional spherical derivates are calculated by evaluating the angular derivatives of the spherical harmonic expansions of the quantities on the worldtube.

B.2.1. Jacobians

To perform coordinate transformations, Jacobians and their derivatives were derived. They included the spherical-to-Cartesian $\Lambda _{i}^{(r,A)}$ and its derivative $\Lambda _{i}^{(r,A)}{{,}_{\alpha }}$, the Cartesian-to-spherical $\Lambda _{(r,A)}^{i}$ and its derivative $\Lambda _{(r,A)}^{i}{{,}_{\alpha }}$.

B.2.2. Null generator

We calculate the unit normal vectors to the time slice $({{n}^{\gamma }})$ and the sphere $({{s}^{\gamma }})$ to calculate the null generator $({{l}^{\gamma }})$

$$\begin{eqnarray}&&{{n}^{\gamma }}=\left( 1/\alpha ,-{{\beta }^{i}}/\alpha \right),\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&{{s}^{\gamma }}=\left( 0,\frac{{{g}^{ij}}{{x}_{j}}}{\sqrt{{{g}^{ij}}{{x}_{i}}{{x}_{j}}}} \right),\end{eqnarray} \tag{ B.7 }$$

$$\begin{eqnarray}&&{{l}^{\gamma }}=\frac{{{n}^{\gamma }}+{{s}^{\gamma }}}{\alpha -{{\beta }^{i}}{{s}^{j}}{{g}_{ij}}}.\end{eqnarray} \tag{ B.8 }$$

Their time derivatives are

$$\begin{eqnarray}&&{{n}^{\alpha }}{{,}_{t}}=\frac{1}{{{\alpha }^{2}}}(-\alpha {{,}_{t}},\alpha {{,}_{t}}{{\beta }^{i}}-\alpha {{\beta }^{i}}{{,}_{t}}),\end{eqnarray} \tag{ B.9 }$$

$$\begin{eqnarray}&&{{s}^{i}}{{,}_{t}}=(-{{g}^{im}}+{{s}^{i}}{{s}^{m}}/2){{g}_{mn}}{{,}_{t}}{{s}^{n}},\end{eqnarray} \tag{ B.10 }$$

$$\begin{eqnarray}&&{{s}^{\alpha }}{{,}_{t}}=\left( 0,{{s}^{i}}{{,}_{t}} \right),\end{eqnarray} \tag{ B.11 }$$

$$\begin{eqnarray}&&{{l}^{\alpha }}{{,}_{t}}=\frac{{{n}^{\alpha }}{{,}_{t}}+{{s}^{\alpha }}{{,}_{t}}-{{l}^{\alpha }}\alpha {{,}_{t}}+{{l}^{\alpha }}({{g}_{ij}}{{,}_{t}}{{\beta }^{i}}{{s}^{j}}+{{g}_{ij}}{{\beta }^{i}}{{,}_{t}}{{s}^{j}}+{{g}_{ij}}{{\beta }^{i}}{{s}^{j}}{{,}_{t}})}{\alpha -{{g}_{ij}}{{\beta }^{i}}{{s}^{j}}}.\end{eqnarray} \tag{ B.12 }$$

B.2.3. (Affine) spherical metric quantities

Dramatic simplification is possible by calculating a number of auxiliary metric terms in a spherical coordinate system. Here, $\lambda$ denotes a derivative in the null direction, whereas A, B, C etc sub- and superscripts denote indices across the spherical components of the coordinate system, i.e. θ and ϕ

$$\begin{eqnarray}&&{{g}_{\alpha \beta }}{{,}_{\lambda }}={{l}^{\gamma }}{{g}_{\alpha \beta }}{{,}_{\gamma }},\end{eqnarray} \tag{ B.13 }$$

$$\begin{eqnarray}&&{{\eta }_{AB}}=\Lambda _{A}^{i}\Lambda _{B}^{j}{{g}_{ij}},\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&{{\eta }_{AB}}{{,}_{\lambda }}=\Lambda _{A}^{k}\Lambda _{B}^{l}{{g}_{kl}}{{,}_{\lambda }}+({{l}^{\mu }}{{,}_{A}}\Lambda _{B}^{l}+{{l}^{\mu }}{{,}_{B}}\Lambda _{A}^{l}){{g}_{\mu l}},\end{eqnarray} \tag{ B.15 }$$

$$\begin{eqnarray}&&{{\eta }_{AB}}{{,}_{t}}=\Lambda _{A}^{i}\Lambda _{B}^{j}{{g}_{ij}}{{,}_{t}},\end{eqnarray} \tag{ B.16 }$$

$$\begin{eqnarray}&&{{\eta }_{uA}}{{,}_{\lambda }}={{l}^{t}}{{,}_{A}}{{g}_{tt}}+\Lambda _{A}^{k}{{g}_{tk}}{{,}_{\lambda }}+{{l}^{k}}{{,}_{A}}{{g}_{tk}}+\Lambda _{A}^{k}{{l}^{\mu }}{{,}_{t}}{{g}_{\mu k}}.\end{eqnarray} \tag{ B.17 }$$

The contravariant quantities are similarly defined

$$\begin{eqnarray}&&{{\eta }^{AB}}{{\eta }_{BC}}=\delta _{C}^{A},\end{eqnarray} \tag{ B.18 }$$

$$\begin{eqnarray}&&{{\eta }^{A\lambda }}={{\eta }^{AB}}\Lambda _{B}^{k}{{g}_{tk}},\end{eqnarray} \tag{ B.19 }$$

$$\begin{eqnarray}&&{{\eta }^{\lambda \lambda }}=-{{g}_{tt}}+{{\eta }^{A\lambda }}\Lambda _{A}^{k}{{g}_{tk}},\end{eqnarray} \tag{ B.20 }$$

$$\begin{eqnarray}&&{{\eta }^{AB}}{{,}_{\lambda }}=-{{\eta }^{AC}}{{\eta }^{BD}}{{\eta }_{CD}}{{,}_{\lambda }},\end{eqnarray} \tag{ B.21 }$$

$$\begin{eqnarray}&&{{\eta }^{A\lambda }}{{,}_{\lambda }}={{\eta }^{AB}}{{\eta }_{uA}}{{,}_{\lambda }}-{{\eta }^{AB}}{{\eta }^{C\lambda }}{{\eta }_{CB}}{{,}_{\lambda }}.\end{eqnarray} \tag{ B.22 }$$

B.2.4. Dyad quantities

The dyad allows construction of spin-weighted scalars on the sphere. Given $r=\sqrt{{{x}^{2}}+{{y}^{2}}+{{z}^{2}}},$

$$\begin{eqnarray}&&{{q}^{i}}=\frac{1}{r\sqrt{{{x}^{2}}+{{y}^{2}}}}\left( -xz+{\rm i}yr,-yz-{\rm i}xr,{{x}^{2}}+{{y}^{2}} \right),\end{eqnarray} \tag{ B.23 }$$

$$\begin{eqnarray}&&{{q}_{A}}=\Lambda _{A}^{i}{{q}^{i}}/r\ ,\end{eqnarray} \tag{ B.24 }$$

since ${{q}_{i}}={{q}^{i}}$. Explicitly

$$\begin{eqnarray}&&{{q}_{A}}=\left( -1,-{\rm i}\frac{{{r}_{c}}}{r} \right),\end{eqnarray} \tag{ B.25 }$$

$$\begin{eqnarray}&&{{q}^{A}}=\left( -1,-{\rm i}\frac{r}{{{r}_{c}}} \right),\end{eqnarray} \tag{ B.26 }$$

where ${{r}_{c}}=\sqrt{{{x}^{2}}+{{y}^{2}}},$ and its derivatives are defined as follows

$$\begin{eqnarray}\begin{array}{rcl} q_{,t}^{i} & = & \frac{1}{{{r}^{3}}r_{c}^{3}}\left( z(({{x}^{4}}-{{y}^{2}}({{y}^{2}}+{{z}^{2}})){{x}_{,t}}+xy(2r_{c}^{2}+{{z}^{2}}){{y}_{,t}}) \right. \\ {} & {} & -xr_{c}^{4}{{z}_{,t}}+{\rm i}{{r}^{3}}(x(-y{{x}_{,t}}+x{{y}_{,t}})),z(xy(2r_{c}^{2}+{{z}^{2}}){{x}_{,t}}-({{x}^{4}}-{{y}^{4}}+{{x}^{2}}{{z}^{2}}){{y}_{,t}}) \\ {} & {} & -yr_{c}^{4}{{z}_{,t}}+{\rm i}{{r}^{3}}(y(-y{{x}_{,t}}+x{{y}_{,t}})),\left. r_{c}^{2}(z(z(x{{x}_{,t}}+y{{y}_{,t}})-r_{c}^{2}{{z}_{,t}})) \right), \\ \end{array}\end{eqnarray} \tag{ B.27 }$$

$$\begin{eqnarray}&&q_{,x}^{i}=\frac{1}{{{r}^{3}}r_{c}^{3}}\left( z({{x}^{4}}-{{y}^{2}}({{y}^{2}}+{{z}^{2}}))-{\rm i}{{r}^{3}}xy,xyz(2r_{c}^{2}+{{z}^{2}})-{\rm i}{{r}^{3}}{{y}^{2}},x{{z}^{2}}r_{c}^{2} \right),\end{eqnarray} \tag{ B.28 }$$

$$\begin{eqnarray}&&q_{,y}^{i}=\frac{1}{{{r}^{3}}r_{c}^{3}}\left( xyz(2r_{c}^{2}+{{z}^{2}})+{\rm i}{{x}^{2}}{{r}^{3}},-z({{x}^{4}}-{{y}^{4}}+{{x}^{2}}{{z}^{2}})+{\rm i}xy{{r}^{3}},y{{z}^{2}}r_{c}^{2} \right),\end{eqnarray} \tag{ B.29 }$$

$$\begin{eqnarray}&&q_{,z}^{i}=\frac{{{r}_{c}}}{{{R}^{3}}}\left( -x,-y,-z \right).\end{eqnarray} \tag{ B.30 }$$

Then ${{q}^{i}}{{,}_{\lambda }}={{l}^{\alpha }}{{q}^{i}}{{,}_{\alpha }}$. Note that the qⁱ are not necessarily constant $(q_{,t}^{i}\ne 0)$ as the properties of the worldtube are time-dependent. This approach differs from the Pitt null code.

B.2.5. The Bondi r

The Bondi radius r_b is an areal radius, not a coordinate radius. The value of the Bondi radius r_b at every point on the worldtube provides the compactification parameter. Given $|{{q}_{AB}}|={{{\rm sin} }^{2}}\;\theta =\frac{{{x}^{2}}+{{y}^{2}}}{{{x}^{2}}+{{y}^{2}}+{{z}^{2}}},$

$$\begin{eqnarray}&&{{r}_{b}}={{\left( \frac{|{{\eta }_{AB}}|}{|{{q}_{AB}}|} \right)}^{\frac{1}{4}}},\end{eqnarray} \tag{ B.31 }$$

$$\begin{eqnarray}&&{{r}_{b}}{{,}_{\lambda }}=\frac{{{r}_{b}}}{4}{{\eta }^{AB}}{{\eta }_{AB}}{{,}_{\lambda }},\end{eqnarray} \tag{ B.32 }$$

$$\begin{eqnarray}&&{{r}_{b}}{{,}_{t}}=\frac{{{r}_{b}}}{4}{{\eta }^{AB}}{{\eta }_{AB}}{{,}_{t}}.\end{eqnarray} \tag{ B.33 }$$

In contrast, derivatives of the coordinate radius r are

$$\begin{eqnarray}&&r{{,}_{\lambda }}={{l}^{\alpha }}r{{,}_{\alpha }},\end{eqnarray} \tag{ B.34 }$$

$$\begin{eqnarray}&&r{{,}_{t}}=\frac{xx{{,}_{t}}+yy{{,}_{t}}+zz{{,}_{t}}}{r},\end{eqnarray} \tag{ B.35 }$$

$$\begin{eqnarray}&&r{{,}_{\alpha }}=\frac{1}{r}\left( {{x}^{i}}{{x}_{i,t}},{{x}^{i}} \right).\end{eqnarray} \tag{ B.36 }$$

In the present algorithm the Cartesian extraction radius r is held constant, but derivatives are included for the general case of a variable extraction coordinate radius.

B.2.6. Derivatives of hypersurface quantities

Calculation of Q and $\Phi =J{{,}_{u}}$ requires the null-derivatives of several other hypersurface quantities

$$\begin{eqnarray}&&\beta {{,}_{\lambda }}=\frac{{{r}_{b}}}{8{{r}_{b}}{{,}_{\lambda }}}\left( J{{,}_{\lambda }}\bar{J}{{,}_{\lambda }}-\frac{1}{1+J\bar{J}}\mathbb{R}{{(\bar{J}J{{,}_{\lambda }})}^{2}} \right),\end{eqnarray} \tag{ B.37 }$$

$$\begin{eqnarray}&&J{{,}_{\lambda }}=2\;J\left( \frac{r{{,}_{\lambda }}}{r}-\frac{{{r}_{b}}{{,}_{\lambda }}}{{{r}_{b}}} \right)+\frac{{{r}^{2}}}{r_{b}^{2}}({{g}_{ij}}{{q}^{i}}{{,}_{\lambda }}{{q}^{j}}+\frac{1}{2}{{g}_{ij}}{{,}_{\lambda }}{{q}^{i}}{{q}^{j}}),\end{eqnarray} \tag{ B.38 }$$

$$\begin{eqnarray}&&U{{,}_{\lambda }}=2\beta {{,}_{\lambda }}U+2\beta {{,}_{\lambda }}{{\eta }^{A\lambda }}{{q}_{A}}-{{\eta }^{A\lambda }}{{,}_{\lambda }}{{q}_{A}}-{{r}_{b}}{{,}_{\lambda B}}{{\eta }^{AB}}{{q}_{A}}+{{r}_{b}}{{,}_{B}}{{\eta }^{AB}}{{,}_{\lambda }}{{q}_{A}}.\end{eqnarray} \tag{ B.39 }$$

B.3. Hypersurface quantities

The hypersurface quanties are J, K, β, U, Q and W, as well as the time derivative $J{{,}_{u}}$

$$\begin{eqnarray}&&J=\frac{{{r}^{2}}}{2r_{b}^{2}}{{g}_{ij}}{{q}^{i}}{{q}^{j}},\end{eqnarray} \tag{ B.40 }$$

$$\begin{eqnarray}&&\Phi =J{{,}_{u}}=2\;J\left( \frac{r{{,}_{t}}}{r}-\frac{{{r}_{b}}{{,}_{t}}}{{{r}_{b}}} \right)+\frac{{{r}^{2}}}{r_{b}^{2}}({{g}_{ij}}{{q}^{i}}{{,}_{t}}{{q}^{j}}+\frac{1}{2}{{g}_{ij}}{{,}_{t}}{{q}^{i}}{{q}^{j}}),\end{eqnarray} \tag{ B.41 }$$

$$\begin{eqnarray}&&K=\sqrt{1+J\bar{J}},\end{eqnarray} \tag{ B.42 }$$

$$\begin{eqnarray}&&\beta =-\frac{1}{2}{\rm log} ({{r}_{b}}{{,}_{\lambda }}),\end{eqnarray} \tag{ B.43 }$$

$$\begin{eqnarray}&&U=-\left( {{\eta }^{\lambda A}}+\frac{{{r}_{b}}{{,}_{B}}}{{{r}_{b}}{{,}_{\lambda }}}{{\eta }^{AB}} \right){{q}_{A}},\end{eqnarray} \tag{ B.44 }$$

$$\begin{eqnarray}&&Q=r_{b}^{2}(J\bar{U}{{,}_{\lambda }}+\sqrt{1+J\bar{J}}U{{,}_{\lambda }}),\end{eqnarray} \tag{ B.45 }$$

$$\begin{eqnarray}&&W=\frac{{{r}_{b}}{{,}_{\lambda }}{{\eta }^{\lambda \lambda }}}{{{r}_{b}}}-\frac{2{{r}_{b}}{{,}_{u}}}{{{r}_{b}}}-\frac{1}{{{r}_{b}}}+\frac{2{{r}_{b}}{{,}_{A}}{{\eta }^{A\lambda }}}{{{r}_{b}}}+\frac{{{r}_{b}}{{,}_{A}}{{r}_{b}}{{,}_{B}}{{\eta }^{AB}}}{{{r}_{b}}{{,}_{\lambda }}{{r}_{b}}}.\end{eqnarray} \tag{ B.46 }$$

B.4. Bondi metric reconstruction

J, K, β, U, Q and W can be combined to reconstruct the Bondi metric, with standard coordinate 4-vector ordering $(u,r,\;\theta ,\phi )$ [20]

$$\begin{eqnarray}{{\eta }^{\alpha \beta }}=\left( \begin{array}{ccccccccccccccc} 0 & -{{{\rm e}}^{-2\beta }} & 0 & 0 \\ -{{{\rm e}}^{-2\beta }} & (rW+1){{{\rm e}}^{-2\beta }} & \frac{1}{2}(U+\bar{U}){{{\rm e}}^{-2\beta }} & -\frac{{\rm i}r}{2{{r}_{c}}}(U-\bar{U}){{{\rm e}}^{-2\beta }} \\ 0 & \frac{1}{2}(U+\bar{U}){{{\rm e}}^{-2\beta }} & \frac{1}{2{{r}^{2}}}(2\;K-J-\bar{J}) & \frac{{\rm i}}{2r{{r}_{c}}}(J-\bar{J}) \\ 0 & -\frac{{\rm i}r}{2{{r}_{c}}}(U-\bar{U}){{{\rm e}}^{-2\beta }} & \frac{{\rm i}}{2r{{r}_{c}}}(J-\bar{J}) & \frac{1}{2r_{c}^{2}}(2\;K+J+\bar{J}) \\ \end{array} \right),\end{eqnarray} \tag{ B.47 }$$

$$\begin{eqnarray}&&\begin{array}{l} {{\eta }_{\alpha \beta }}\;= \\ \left( \begin{array}{ccccccccccccccc} -(rW+1){{{\rm e}}^{2\beta }} & -{{{\rm e}}^{2\beta }} & \frac{{{r}^{2}}}{2}\left( (J+K)\bar{U} \right. & \frac{{\rm i}r{{r}_{c}}}{2}\left( \bar{J}U-J\bar{U} \right. \\ +\frac{{{r}^{2}}}{2}\left( 2\;KU\bar{U}+J{{{\bar{U}}}^{2}}+\bar{J}{{U}^{2}} \right) & {} & +\left. (\bar{J}+K)U \right) & +\left. K(\bar{U}-U) \right) \\ {} & {} & {} & {} \\ -{{{\rm e}}^{2\beta }} & 0 & 0 & 0 \\ \frac{{{r}^{2}}}{2}\left( (J+K)\bar{U}+(\bar{J}+K)U \right) & 0 & \frac{{{r}^{2}}}{2}(J+\bar{J}+2\;K) & \frac{{\rm i}r{{r}_{c}}}{2}(\bar{J}-J) \\ \frac{{\rm i}r{{r}_{c}}}{2}\left( \bar{J}U-J\bar{U}+K(\bar{U}-U) \right) & 0 & \frac{{\rm i}r{{r}_{c}}}{2}(\bar{J}-J) & -\frac{r_{c}^{2}}{2}(J+\bar{J}-2\;K) \\ \end{array} \right). \\ \end{array}\end{eqnarray} \tag{ B.48 }$$

The hypersurface equations are written making use of the ð formalism [22], which simplifies equations with spherical symmetry. In our implementation, we use a basis capable of multiple values on the poles, so only a single patch is necessary to cover the entire sphere.

Rank 2 tensors on the sphere can be broken down and expressed as a sum of spin-weighted spherical harmonics when contracted with an appropriate dyad q^A. Following [22](equation (8)), a generic tensor can be written as the sum of respectively a symmetric trace free part, a trace part, and an antisymmetric part:

$$\begin{eqnarray}&&{{w}_{AB}}={{t}_{AB}}+\frac{p}{2}{{q}_{AB}}+{\rm i}\frac{u}{4}\left( {{q}_{A}}\overline{{{q}_{B}}}-\overline{{{q}_{A}}}{{q}_{B}} \right).\end{eqnarray} \tag{ C.1 }$$

$p={{q}^{CD}}{{w}_{CD}}$ and $u={\rm i}({{q}^{C}}{{\bar{q}}^{D}}-{{\bar{q}}^{C}}{{q}^{D}}){{w}_{CD}}/2$ are both real scalar fields with spin-weight 0. Writing $t={{t}_{AB}}{{q}^{A}}{{q}^{B}}$ yields a complex scalar field with spin-weight 2. Together, these three scalars p, u and t completely specify the tensor field independent of choice of basis.

ð is a spherical derivative operator on spin-weighted spherical harmonics. In spherical coordinates

$$\begin{eqnarray}&&\eth \eta =-({{{\rm sin} }^{s}}\;\theta )\left( \frac{\partial }{\partial \;\theta }+\frac{{\rm i}}{{\rm sin} \;\theta }\frac{\partial }{\partial \phi } \right)({{{\rm sin} }^{-s}}\;\theta \;\eta ),\end{eqnarray} \tag{ C.2 }$$

$$\begin{eqnarray}&&\bar{\eth }\eta =-({{{\rm sin} }^{-s}}\;\theta )\left( \frac{\partial }{\partial \;\theta }-\frac{{\rm i}}{{\rm sin} \;\theta }\frac{\partial }{\partial \phi } \right)({{{\rm sin} }^{s}}\;\theta \;\eta ).\end{eqnarray} \tag{ C.3 }$$

This can be thought of as a contraction of the dyad with the spherical derivatives operator.

In terms of spin-weighted spherical harmonics

$$\begin{eqnarray}&&\eth \left( _{s}{{Y}_{\operatorname{lm}}} \right)=+{{\sqrt{(l-s)(l+s+1)}}_{s+1}}{{Y}_{\operatorname{lm}}},\end{eqnarray} \tag{ C.4 }$$

$$\begin{eqnarray}&&\bar{\eth }\left( _{s}{{Y}_{\operatorname{lm}}} \right)=-{{\sqrt{(l+s)(l-s+1)}}_{s+1}}{{Y}_{\operatorname{lm}}}.\end{eqnarray} \tag{ C.5 }$$

Used in combination, these definitions allow the ð formalism to be used in the spectral domain. We used the package Spinsfast [27] to implement spin-weighted spherical harmonics.

ð assumes a spherical domain. However, in general the compactification parameter $R={{r}_{b|\Gamma }}$ is not constant on the sphere. The corrected operation is

$$\begin{eqnarray}&&\eth \eta =\tilde{\eth }\eta -\eta {{,}_{\rho }}\frac{\rho (1-\rho )}{R}\tilde{\eth }R,\end{eqnarray} \tag{ C.6 }$$

where $\tilde{\eth }$ denotes the naive operator in the affine coordinate system.

This correction is the same that operates at the end of the calculation of $J{{,}_{u}}=\Phi$ in appendix D, only in reverse [28].

These are largely drawn from Bishop et al [20], with some re-arrangement of terms to ensure internal consistency with the Magnus expansion formalism, and the compactification transformation $\rho =r/(R+r)$. Additionally, terms are grouped by their order in the Laurent expansion, where relevant.

D.1. β terms

$$\begin{eqnarray}&&{{\beta }_{,r}}={{N}_{\beta }},\end{eqnarray} \tag{ D.1 }$$

$$\begin{eqnarray}&&{{N}_{\beta }}=\frac{\rho {{(1-\rho )}^{3}}}{8R}(J{{,}_{\rho }}\bar{J}{{,}_{\rho }}-{{K}^{2}}{{,}_{\rho }}).\end{eqnarray} \tag{ D.2 }$$

D.2. Q terms

$$\begin{eqnarray}&&{{({{r}^{2}}Q)}_{,\rho }}=\frac{{{Q}_{C}}+{{Q}_{{\rm CNL}}}}{{{(1-\rho )}^{2}}}+\frac{{{Q}_{D}}}{{{(1-\rho )}^{3}}}.\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&{{Q}_{C}}=-{{R}^{2}}{{\rho }^{2}}\bar{\eth }J{{,}_{\rho }}-{{R}^{2}}{{\rho }^{2}}\eth K{{,}_{\rho }}+2{{R}^{2}}{{\rho }^{2}}\eth \beta {{,}_{\rho }},\end{eqnarray} \tag{ D.4 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{Q}_{{\rm CNL}}} & = & {{R}^{2}}{{\rho }^{2}}(1-K)(\eth K{{,}_{\rho }}+\bar{\eth }J{{,}_{\rho }})+\bar{J}\eth J{{,}_{\rho }}+J{{,}_{\rho }}\bar{J}{{,}_{\rho }}+J\bar{\eth }K{{,}_{\rho }} \\ {} & {} & +K{{,}_{\rho }}\bar{\eth }J-J{{,}_{\rho }}\bar{\eth }K+\frac{\eth \bar{J}(J{{,}_{\rho }}-{{J}^{2}}\bar{J}{{,}_{\rho }})+\eth J(\bar{J}{{,}_{\rho }}-{{{\bar{J}}}^{2}}J{{,}_{\rho }})}{2\;{{K}^{2}}}, \\ \end{array}\end{eqnarray} \tag{ D.5 }$$

$$\begin{eqnarray}&&{{Q}_{D}}=-4{{R}^{2}}\rho \eth \beta .\end{eqnarray} \tag{ D.6 }$$

D.3. U terms

$$\begin{eqnarray}&&U{{,}_{\rho }}={{U}_{A}}+{{U}_{{\rm ANL}}}.\end{eqnarray} \tag{ D.7 }$$

$$\begin{eqnarray}&&{{U}_{A}}=\frac{{{{\rm e}}^{2\beta }}}{R{{\rho }^{2}}}Q,\end{eqnarray} \tag{ D.8 }$$

$$\begin{eqnarray}&&{{U}_{{\rm ANL}}}=\frac{{{{\rm e}}^{2\beta }}}{R{{\rho }^{2}}}(KQ-Q-J\bar{Q}).\end{eqnarray} \tag{ D.9 }$$

D.4. W terms

$$\begin{eqnarray}&&({{r}^{2}}W){{,}_{\rho }}=\frac{{{W}_{C}}+{{W}_{{\rm CNL}}}}{{{(1-\rho )}^{2}}}+\frac{{{W}_{D}}}{{{(1-\rho )}^{3}}}.\end{eqnarray} \tag{ D.10 }$$

$$\begin{eqnarray}&&{{W}_{C}}=\mathcal{R}-{{R}^{2}}{{{\rm e}}^{2\beta }}(\eth \beta \bar{\eth }\beta +\eth \bar{\eth }\beta )+\frac{{{R}^{2}}{{\rho }^{2}}}{4}(\eth \bar{U}{{,}_{\rho }}+\bar{\eth }U{{,}_{\rho }}),\;{\rm where}\end{eqnarray} \tag{ D.11 }$$

$$\begin{eqnarray}&&\mathcal{R}=\frac{R{{{\rm e}}^{2\beta }}}{2}\left( 2\;K-\eth \bar{\eth }K+\frac{1}{2}({{{\bar{\eth }}}^{2}}\;J+{{\eth }^{2}}\bar{J})+\frac{1}{4\;K}(\bar{\eth }\bar{J}\eth J-\bar{\eth }J\eth \bar{J}) \right),\end{eqnarray} \tag{ D.12 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{W}_{{\rm CNL}}} & = & R{{{\rm e}}^{2\beta }}\left( (1-K)(\eth \bar{\eth }\beta +\eth \beta \bar{\eth }\beta )+\frac{1}{2}(J\bar{\eth }{{\beta }^{2}}+\bar{J}\eth {{\beta }^{2}}) \right. \\ {} & {} & -\left. \frac{1}{2}(\eth \beta (\bar{\eth }K-\eth \bar{J})+\bar{\eth }\beta (\eth K-\bar{\eth }J))+\frac{1}{2}(J{{{\bar{\eth }}}^{2}}\beta +\bar{J}{{\eth }^{2}}\beta ) \right) \\ {} & {} & -{{{\rm e}}^{2\beta }}\frac{{{R}^{3}}{{\rho }^{4}}}{8}\left( 2\;KU{{,}_{\rho }}\bar{U}{{,}_{\rho }}+J{{{\bar{U}}}^{2}}{{,}_{\rho }}+\bar{J}{{U}^{2}}{{,}_{\rho }} \right), \\ \end{array}\end{eqnarray} \tag{ D.13 }$$

$$\begin{eqnarray}&&{{W}_{D}}={{R}^{2}}\rho (\eth \bar{U}+\bar{\eth }U).\end{eqnarray} \tag{ D.14 }$$

D.5.  $J{{,}_{u}}=\Phi$ terms

$$\begin{eqnarray}&&(r\Phi ){{,}_{\rho }}-(rJ)(\Phi \bar{\Gamma }+\bar{\Phi }\Gamma )={{\Phi }_{A}}+{{\Phi }_{{\rm ANL}}}+\frac{{{\Phi }_{B}}+{{\Phi }_{{\rm BNL}}}}{1-\rho }+\frac{{{\Phi }_{C}}+{{\Phi }_{{\rm CNL}}}}{{{(1-\rho )}^{2}}}.\end{eqnarray} \tag{ D.15 }$$

$$\begin{eqnarray}&&\Gamma =\left( J{{,}_{\rho }}-J\frac{K{{,}_{\rho }}}{K} \right)\ .\end{eqnarray} \tag{ D.16 }$$

$$\begin{eqnarray}&&{{\Phi }_{A}}=(1-\rho )J{{,}_{\rho }}+\frac{1}{2}R{{\rho }^{2}}W{{,}_{\rho }}J{{,}_{\rho }}+\frac{\rho }{2}(1-\rho +R\rho W)J{{,}_{\rho \rho }},\end{eqnarray} \tag{ D.17 }$$

$$\begin{eqnarray}&&{{\Phi }_{B}}=\left( \frac{3}{2}R\rho -R{{\rho }^{2}} \right)WJ{{,}_{\rho }}-\frac{1}{2}R\rho \eth U{{,}_{\rho }}+\frac{{{{\rm e}}^{2\beta }}}{\rho }({{\eth }^{2}}\beta +\eth {{\beta }^{2}}),\end{eqnarray} \tag{ D.18 }$$

$$\begin{eqnarray}&&{{\Phi }_{C}}=-R\eth U,\end{eqnarray} \tag{ D.19 }$$

$$\begin{eqnarray}&&{{\Phi }_{{\rm ANL}}}=-4J\beta {{,}_{\rho }},\end{eqnarray} \tag{ D.20 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\Phi }_{{\rm BNL}}} & = & {{N}_{1B}}+{{N}_{2B}}+{{N}_{3B}}+{{N}_{4B}}+{{N}_{5B}}+{{N}_{6B}}+{{N}_{7B}} \\ {} & {} & +{{P}_{1B}}+{{P}_{2B}}+{{P}_{3B}}+{{P}_{4B}}, \\ \end{array}\end{eqnarray} \tag{ D.21 }$$

$$\begin{eqnarray}&&{{\Phi }_{{\rm CNL}}}={{N}_{2C}}+{{N}_{3C}}+{{P}_{3C}},\end{eqnarray} \tag{ D.22 }$$

N terms are as in Bishop et al [20], with a prefactor $\frac{R}{2{{(1-\rho )}^{2}}}$ and the usual compactification transformation

$$\begin{eqnarray}\begin{array}{rcl} {{N}_{1B}} & = & -\frac{{{{\rm e}}^{2\beta }}}{2\rho }\left( K(\eth J\bar{\eth }\beta +2\eth K\eth \beta -\bar{\eth }J\eth \beta ) \right. \\ {} & {} & +\left. J(\bar{\eth }J\bar{\eth }\beta -2\eth K\bar{\eth }\beta )-\bar{J}\eth J\eth \beta \right), \\ \end{array}\end{eqnarray} \tag{ D.23 }$$

$$\begin{eqnarray}&&{{N}_{2B}}=-\frac{R\rho }{4}\left( \eth J\bar{U}{{,}_{\rho }}+\bar{\eth }JU{{,}_{\rho }} \right),\end{eqnarray} \tag{ D.24 }$$

$$\begin{eqnarray}&&{{N}_{2C}}=-\frac{R}{2}\left( \eth J\bar{U}+\bar{\eth }JU \right),\end{eqnarray} \tag{ D.25 }$$

$$\begin{eqnarray}&&{{N}_{3B}}=\frac{R\rho }{2}\left( (1-K)\eth U{{,}_{\rho }}-J\eth \bar{U}{{,}_{\rho }} \right),\end{eqnarray} \tag{ D.26 }$$

$$\begin{eqnarray}&&{{N}_{3C}}=R(1-K)\eth U-RJ\eth \bar{U},\end{eqnarray} \tag{ D.27 }$$

$$\begin{eqnarray}&&{{N}_{4B}}=\frac{{{{\rm e}}^{-2\beta }}}{4}{{R}^{2}}{{\rho }^{3}}\left( {{K}^{2}}{{U}^{2}}{{,}_{\rho }}+2JKU{{,}_{\rho }}\bar{U}{{,}_{\rho }}+{{J}^{2}}{{{\bar{U}}}^{2}}{{,}_{\rho }} \right),\end{eqnarray} \tag{ D.28 }$$

$$\begin{eqnarray}&&{{N}_{5B}}=-\frac{R\rho }{4}J{{,}_{\rho }}(\eth \bar{U}+\bar{\eth }U),\end{eqnarray} \tag{ D.29 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{N}_{6B}} & = & \frac{R\rho }{2}\left( \frac{1}{2}(\bar{U}\eth J+U\bar{\eth }J)(J\bar{J}{{,}_{\rho }}-\bar{J}J{{,}_{\rho }})+\bar{U}\bar{\eth }J(JK{{,}_{\rho }}-KJ{{,}_{\rho }}) \right. \\ {} & {} & -\left. \bar{U}(\eth J{{,}_{\rho }}-2\;K\eth KJ{{,}_{\rho }}+2\;J\eth KK{{,}_{\rho }})-U(\bar{\eth }J{{,}_{\rho }}-K\eth \bar{J}J{{,}_{\rho }}+J\eth \bar{J}K{{,}_{\rho }}) \right), \\ \end{array}\end{eqnarray} \tag{ D.30 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{N}_{7B}} & = & \frac{R\rho }{2}(KJ{{,}_{\rho }}-JK{{,}_{\rho }})\left( \bar{U}(\bar{\eth }J-\eth K)+U(\bar{\eth }K-\eth \bar{J}) \right. \\ {} & {} & +\left. K(\bar{\eth }U-\eth \bar{U})+(J\bar{\eth }\bar{U}-\bar{J}\eth U) \right). \\ \end{array}\end{eqnarray} \tag{ D.31 }$$

P terms are as in Bishop et al [20] except for the terms in $J{{,}_{u}}$, which have been moved to the left-hand side of the equation. They have a prefactor $\frac{J}{2\rho (1-\rho )}$ and the usual compactification transformation

$$\begin{eqnarray}&&{{P}_{1A}}=-4J\beta {{,}_{\rho }},\end{eqnarray} \tag{ D.32 }$$

$$\begin{eqnarray}&&{{P}_{1B}}=-4R\rho J\beta {{,}_{\rho }}W,\end{eqnarray} \tag{ D.33 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{P}_{2B}} & = & \frac{{{{\rm e}}^{2\beta }}}{2\rho }J\left( -2K(\eth \bar{\eth }\beta +\bar{\eth }\beta \eth \beta )-(\bar{\eth }\beta \eth K+\eth \beta \bar{\eth }K)+J({{{\bar{\eth }}}^{2}}\beta +\bar{\eth }{{\beta }^{2}}) \right. \\ {} & {} & +\left. \bar{J}({{\eth }^{2}}\beta +\eth {{\beta }^{2}})+(\bar{\eth }J\bar{\eth }\beta +\eth \bar{J}\eth \beta ) \right), \\ \end{array}\end{eqnarray} \tag{ D.34 }$$

$$\begin{eqnarray}&&{{P}_{3B}}=\frac{R\rho }{4}J(\bar{\eth }U{{,}_{\rho }}+\eth \bar{U}{{,}_{\rho }}),\end{eqnarray} \tag{ D.35 }$$

$$\begin{eqnarray}&&{{P}_{3C}}=\frac{R}{2}J(\bar{\eth }U+\eth \bar{U}),\end{eqnarray} \tag{ D.36 }$$

$$\begin{eqnarray}&&{{P}_{4B}}=-\frac{{{{\rm e}}^{-2\beta }}}{2}{{R}^{2}}{{\rho }^{3}}J\left( 2KU{{,}_{\rho }}\bar{U}{{,}_{\rho }}+J{{{\bar{U}}}^{2}}{{,}_{\rho }}+\bar{J}{{U}^{2}}{{,}_{\rho }} \right).\end{eqnarray} \tag{ D.37 }$$

The non-spherical and non-constant inner boundary creates a discrepancy between the Bondi and affine coordinates. The corrective factor is given by

$$\begin{eqnarray}&&f{{,}_{\tilde{u}}}=f{{,}_{u}}+f{{,}_{r}}r{{,}_{u}}=f{{,}_{u}}+\rho (1-\rho )\frac{R{{,}_{u}}}{R}f{{,}_{\rho }},\end{eqnarray} \tag{ D.38 }$$

where once again ${{,}_{\tilde{u}}}$ denotes a derivative performed in the affine coordinate system, whereas ${{,}_{u}}$ is the regular derivative in the Bondi coordinate system. Time steps are performed in the affine coordinate system, whereas the hypersurface equations are calculated and solved in the Bondi coordinate system.

This term must be added at the end of the calculation, in any instance where $R{{,}_{u}}$ is non-zero.