General relativistic null-cone evolutions with a high-order scheme

# General relativistic null-cone evolutions with a high-order scheme

Christian Reisswig111Einstein Fellow Theoretical Astrophysics Including Relativity, California Institute of Technology, Pasadena, CA 91125, USA    Nigel Bishop Department of Mathematics, Rhodes University, Grahamstown, 6139 South Africa    Denis Pollney Department of Mathematics, Rhodes University, Grahamstown, 6139 South Africa
July 6, 2019
###### Abstract

We present a high-order scheme for solving the full non-linear Einstein equations on characteristic null hypersurfaces using the framework established by Bondi and Sachs. This formalism allows asymptotically flat spaces to be represented on a finite, compactified grid, and is thus ideal for far-field studies of gravitational radiation. We have designed an algorithm based on 4th-order radial integration and finite differencing, and a spectral representation of angular components. The scheme can offer significantly more accuracy with relatively low computational cost compared to previous methods as a result of the higher-order discretization. Based on a newly implemented code, we show that the new numerical scheme remains stable and is convergent at the expected order of accuracy.

###### pacs:
04.25.dg, 04.30.Db, 04.30.Tv, 04.30.Nk

## 1 Introduction

Characteristic formulations of the Einstein equations have proven to be an important tool for numerical relativity. Most recently, they have been employed for the practical problem of measuring gravitational waves in a gauge-invariant and unambiguous way from numerically evolved spacetimes of binary black hole mergers [1, 2, 3, 4], rotating stellar core collapse [5], and collapsar formation [6]. The technique, called Cauchy-characteristic extraction (CCE) (see [7] for a review) takes metric boundary data produced by a 3+1 evolution on a worldtube of finite radius, and uses null-cone evolutions of the Einstein equations to transport that metric data to future null infinity, , the conformal outer boundary of spacetime where gravitational radiation is invariantly defined and interpreted in the Bondi gauge [8, 9, 10]. Thus, this technique removes the influence of near-zone and coordinate effects [1, 2, 3, 4, 5].

The null formulation is extremely efficient for evolving fields in the wave zone, where the null coordinates are well-behaved and caustics along geodesics are unlikely to be an issue. On relatively small computational grids, by comparison with standard 3+1 methods, it is possible to achieve an accuracy which has proven to be sufficient for practical applications (gravitational wave measurement from compact bodies). In the case of characteristic gravitational-wave extraction, the inner boundary data for characteristic evolution is constructed on a worldtube at some distance form the source where the curvature gradients are already rather small compared to those close to the source. Thus, characteristic extraction requires comparatively little numerical resolution and is therefore less computationally demanding than 3+1 evolutions in the near-zone of a dynamical source.

Although the computational effort is smaller, it is still non-trivial. To yield sufficient accuracy, for example to extract the gravitational radiation from a binary black hole evolution, the characteristic computation by current methods requires several days on up to a dozen processors on a workstation or small cluster to complete a time-series encompassing a dozen orbits of the binary (where is the dimensionless mass of the spacetime). The application of higher-order discretization schemes, on the other hand, may deliver sufficient accuracy at much lower computational cost so that the additional effort of characteristic extraction would become negligible. If the characteristic code could be run concurrently with the underlying 3+1 simulation, this would allow for on-the-fly extraction of waveforms, as well as the interesting potential of using characteristic methods to provide exact non-linear boundary conditions for 3+1 codes.

Characteristic codes have a long history in numerical relativity. A prominent result was the first stable dynamical evolution of a black hole spacetime in three spatial dimensions achieved by the Pittsburgh group [11, 12]. Since then, the Pittsburgh null code (or PITTNullCode) has become the main building block for current implementations of characteristic extraction used in numerical relativity simulations [1, 2, 3, 4, 5], and is now part of the publicly available Einstein Toolkit [13]. The code employs a single-null coordinate system, and is formulated in terms of spin-weighted variables that are related to the original variables defined by Bondi and collaborators [8, 9, 10]. It is built around a numerical scheme based on points located at the corners of null parallelograms, which was originally shown to be stable [14] in the context of the scalar wave equations. The PITTNullCode in its original form is 2nd-order accurate in space and time discretization, meaning that the error decreases as as the resolution is increased.

Over the years, there have been several improvements to the original algorithm as implemented in [11, 12]. Stereographic coordinates were replaced by more uniform angular grids [15, 16]. Fourth-order accurate angular derivatives have also been introduced, though retaining the 2nd-order parallelogram-based radial integration scheme [15, 3].

An alternative to the 2nd-order accurate null parallelogram was attempted by Bartnik and Norton [17, 18]. They developed an algorithm based on a null quasi-spherical gauge, using method-of-lines integrators and 4th-order in time accuracy. They also used a spherical harmonic decomposition of variables on angular shells, and thus their code was pseudo-spectral in the angular directions. Their formulation of the null evolution equations (and coordinates) lead to certain numerical complications, including the need for high-accuracy interpolation operators to compute radial derivatives, and an elliptic gauge equation. Ultimately, the code did not demonstrate long-term stability.

In this paper, we present a new high-order integration scheme for characteristic evolutions in general relativity. The evolution equations are written in the Bondi form, following the prescription of [12]. Time and outward radial integrations are performed by method-of-lines schemes. In particular, for the time integration, we use a classical 4th-order Runge-Kutta integrator. For the radial direction, we use a modified Adams-Moulton multi-step method. A multistep method is required by the lack of information at points between grid spacings, which is needed by the intermediate steps of the Runge-Kutta schemes. A similar method, for the case of 2nd-order accuracy and axisymmetry, was previously used in [19]. In the present context of full three-dimensional characteristic evolutions, we additionally need to discretize the angular direction. For this, we use spectral expansions in terms of spin-weighted real-valued spherical harmonics. The equations are solved on an angular collocation grid using the pseudo-spectral method. Radial and time integrations are performed on the angular spectral coefficients of the evolution variables.

By referring to a simplified linear model problem in Section 2, we argue analytically that the proposed method is stable. We review the Bondi evolution equations in Section 3. In Section 4, we describe the numerical methods, including time and radial integrators, and pseudo-spectral derivatives. Finally, we test a newly implemented three-dimensional high-order code, demonstrating the expected order of convergence and accuracy of the method. We find that the new scheme is significantly more efficient than the old scheme, reaching the same level of accuracy using only very few radial and angular points.

## 2 Integration schemes and stability

Characteristic evolution schemes based on null or double-null coordinates have a significantly different character than conventional 3+1 (Cauchy) evolutions. In 3+1 evolutions, spacetime is foliated along a timelike vector field by spacelike hypersurfaces . In so-called 2+1+1 characteristic evolutions (the characteristic initial boundary value problem [20]), spacetime is foliated along a timelike vector field by null hypersurfaces, which are characteristic surfaces of the Einstein field equations. In this section, we describe a simplified model characteristic problem which exhibits the important features of the Einstein system, which we outline explicitly in Section 3.

The characteristic method solves for the values of a field, , which obeys a hyperbolic equation of the form

 J,ur=F(J,J,r,J,rr,J,A). (1)

Here is a retarded time coordinate labeling individual null slices, is a radial coordinate along each null slice, and the index labels angular coordinates. We allow the function to depend on , and its partial derivatives (to 2nd-order in ), which we label by subscripts

 ∂2J∂u∂r⇒J,ur. (2)

We consider the problem on a domain bounded on the interior by a timelike worldtube, , at some finite areal radius from the centre of our coordinate system, and on the exterior by . Boundary data consists of the variables required to evaluate on , as well as initial data for along a single null cone at .

We introduce an intermediate variable , which allows us to recast Eq. (1) as a pair of 1st-order equations

 Φ,r = F(J,J,r,J,rr,J,A), (3a) J,u = Φ. (3b)

Equation (3a) does not involve time derivatives and can thus be solved on a null slice by radial integration. Given data for on a slice either from the last time step or by appropriate initial data, we propagate the system forward in time by first solving Eq. (3a) along the radial direction, and then integrating Eq. (3b) forward in time to determine on the next slice.

The schemes which we use for both radial and null integrations fall within the broad class of method-of-lines integrators for partial differential equations. These schemes assume that we have been able to evaluate the RHS at a point, so that standard ordinary differential equation (ODE) solvers can be applied to evolve the functions either forward in time, , or (in the case of ) outward in .

For the time direction, we can use an explicit integrator, such as a standard 4th-order Runge-Kutta scheme. Recall, however, that in taking the solution from to , the Runge-Kutta method involves calculating the right-hand side (RHS) at a number of intermediate steps. In the case of the integration the RHS is , and thus we need to compute via radial integration of Eq. (3a) at intermediate substeps between our timesteps of size .

Applying such a scheme in the radial direction for is more problematic. In that case, the RHS is given by , which is a function of data which is only known on discrete spheres (separated by a fixed distance, ), and so we cannot evaluate the intermediate substeps required by a Runge-Kutta-type integrator. Alternatively, we could make use of a multi-step algorithm, such as Adams-Bashforth or Adams-Moulton. These methods evaluate the RHS over some number of previous points, which in the case of the radial integration correspond to a set of the radial spheres on which and its derivatives can been evaluated.

We have examined a number of numerical schemes for carrying out the radial integration required by Eq. (3a) for the Einstein system, but almost universally found them to be unstable in empirical tests. To investigate the stability of different numerical methods, we turn to a simpler model which embodies the main features of the Einstein system, on which we can carry out a von Neumann analysis.

By setting individual terms in the Einstein equations to zero and examining the subsequent numerical evolution, we came to the conclusion that the key terms determining the stability are those which involve the variables and . Thus we constructed two simplified systems which consist exclusively of these terms in the radial integration of a variable . That is, we considered the systems

 Φ,r=J,r,J,u=Φ, (4)

and

 Φ,r=J,rr,J,u=Φ, (5)

individually. In practice, we replace the areal radius by a compactified coordinate , defined by

 x=1rΓ+r (6)

where is a non-zero parameter corresponding to the radius of the inner boundary (see Eq. (25), below), and have found this transformation to be important in the stability analysis. In terms of , Eqs. (4) and (5) become

 Φ,x+(r+rΓ)2rrΓΦ=rΓr(r+rΓ)J,x,J,u=Φ, (7)

and

 Φ,x+(r+rΓ)2rrΓΦ=rΓ2(r+rΓ)2J,xx,J,u=Φ, (8)

respectively.

The von Neumann analysis corresponds to assuming the following form for the model variables:

 J=ewueikx,Φ=EJ, (9a) Δx=h,Δu=μh,ν=kh2. (9b)

The equations are evolved using a 4th-order Runge-Kutta integration, for which the general stability analysis is quite involved. However, to leading order in , it is the same as the Euler method (and this was also found for other evolution algorithms such as the Adams-Bashforth methods). Thus,

 ewμh=1+μhE, (10)

and the stability is determined by the sign of :

 R(E)⎧⎨⎩>0,\text{Unstable},=0:\text{Stability unknown},<0:Stable. (11)

The quantity is determined by the particular finite difference algorithm used to evaluate Eq. (7) or Eq. (8), followed by substitution of Eqs. (9a, 9b). We investigated both second and 4th-order explicit (Adams-Bashforth) and implicit (Adams-Moulton) multi-step methods, with finite differences evaluated using centred, forward and backward methods of the appropriate accuracy. The calculations are somewhat lengthy and were done using a computer algebra script (see the file vN_comp3.map in the online supplement).

We found that the case Eq. (7) involving always leads to values of for which so that the system is always unstable, independent of the integration scheme for . However, there are cases for which Eq. (8), involving , is stable. Since a second derivative term in the stability analysis is divided by compared to division by for the first derivative term, to leading order in the stability of the second derivative term dominates that of the first derivative term.

The stability analysis of Eq. (8) indicates that among the tested methods, only the following cases are stable:

where by “forward” difference operator we indicate that radial derivatives of are calculated using a stencil which involves points in the positive radial direction. We were able to confirm these results by empirical tests with the simplified system, Eq. (8). Thus, for the full Einstein equations, described in the next section, we implemented a scheme in which radial integrations are carried out using a 4th-order Adams-Moulton method with upwinded radial derivatives.

## 3 The Einstein equations in the Bondi-Sachs framework

### 3.1 Coordinates

The Bondi formulation writes the Einstein equations in terms of a null foliation of an asymptotically flat spacetime. We introduce coordinates . The coordinate is a radial surface area coordinate, and is a retarded time coordinate which replaces the time, , of 3+1 formulations. The are angular coordinates, labeled by uppercase indices that take the values and . The coordinates label individual null geodesics extending from a worldtube . The worldtube is chosen so that on . In these coordinates, the general spacetime line element is

where satisfies

 hABhBC=δAC,det(hAB)=det(qAB), (13)

and is the unit sphere metric.

The spacetime is described by , , , and , which are functions of the coordinates. It is convenient to write quantities in terms of spin-weighted scalars in order to remove explicit angular tensor components. This simplifies the expression of the field equations in a way that is independent of the choice of angular coordinates. To this end, we introduce a complex dyad satisfying , , . 222An explicit form for will not be needed here, since we will represent angular dependence in terms of spin-weighted spherical harmonic basis functions constructed using a particular dyad representation adapted to the choice of angular coordinates.. By projecting the angular variables onto this dyad, we define the complex valued scalars

 J=12hABqAqB,K=12hABqA¯qB,U=UAqA, (14)

of spin weights 2, 0, and 1, respectively. The components of are uniquely determined by due to the determinant condition, Eq. (13), thus fixing as a function of via

 K=√1+J¯J. (15)

Corresponding to the complex dyad, we introduce complex angular covariant differential operators and which maintain the property of spin-weight when acting on a scalar of spin-weight  [21]. The action of the and operators is restricted to transformations of our spin-weighted spherical harmonic basis functions (see Section 4.3).

Before writing out the Einstein equations, we note that it is convenient to introduce an additional intermediate variable defined by

 Q:=r2e−2βhABUB,rqA. (16)

This spin-weight 1 variable, which is the first radial derivatives of , will allow us to write the equations in 1st-order form. Also, we re-express in terms of a new variable

 ^W:=V−rr2, (17)

which has a regular limit as .

### 3.2 Einstein equations

In Bondi coordinates, the vacuum Einstein equations

 Rab=0 (18)

give rise to a hierarchy of equations which we can characterize as (i) hypersurface equations, (ii) evolution equations, and (iii) constraints.

Hypersurface equations do not depend on -derivatives, and thus can be evaluated within a slice. They are determined by the components , , and and lead to the following hierarchy of equations:

 β,r = Nβ, (19a) (r2Q),r = (19b) U,r = r−2e2βQ+NU, (19c) (r2^W),r = 12e2βR−1−eβð¯ðeβ (19d) +14r−2(r4(ð¯U+¯ðU)),r+NW,

where the Ricci scalar is given explicitly by

 R=2K−ð¯ðK+12(¯ð2J+ð2¯J)+14K(¯ð¯JðJ−¯ðJð¯J), (20)

and , , , and are non-linear aspherical terms given explicitly in A. The equations are solved in succession, assuming available data on a slice and constraint satisfying inner boundary data at the worldtube for each of the hypersurface variables. This allows us to solve for , which in turn provides data for the equation for . Given and , we can then solve for , and finally for .

The component of the Einstein equations determines the evolution equation for :

 2(rJ),ur−((1+r^W)(rJ),r),r= −r−1(r2ðU),r+2r−1eβð2eβ−(r^W),rJ+NJ, (21)

where the non-linear aspherical terms have been gathered in the quantity (specified in A). We introduce an intermediate variable

 Φ:=J,u. (22)

In terms of ,  Eq. (21) becomes a new hypersurface equation

 2(rΦ),r−((1+r^W)(rJ),r),r= −r−1(r2ðU),r+2r−1eβð2eβ−(r^W),rJ+NJ. (23)

which is integrated radially from using known values of the hypersurface variables determined in Eq. (19a)-(19d). Then, is determined by a timelike integration of

 J,u=Φ. (24)

We have expressed the Einstein system in a form analogous to the simplified model described in Eq. (3). The source for is complicated, but determined entirely by radial integration. Note the presence of the in the second term of Eq. (23), which we have highlighted in Section 2 as key to determining the stability of numerical evolution schemes.

Finally, we take advantage of the nature of null geodesics in asymptotically flat spacetimes to compactify the radial direction so that is a boundary point of a closed domain. We replace the areal radius by a new coordinate via the invertible coordinate transformation

 x(r)=rrΓ+r,r(x)=rΓx1−x, (25)

where is a constant, which we choose to be the radius of the worldtube, . In this coordinate, the equations have a regular limit as , and furthermore, we are able to set terms of order for , to zero at (see Section 3.3). Throughout the domain, derivatives are evaluated numerically in terms of the new coordinate, , and then transformed into -derivatives using the standard Jacobian transformations. For Eq. (25), these are

 dxdr=rΓ(r+rΓ)2,d2xdr2=−2rΓ(r+rΓ)3. (26)

### 3.3 Form of the equations at J+

In the limit of , corresponding to , the numerical treatment of the equations requires special care. The problematic terms are those involving the Jacobian and the coordinate function , which are not regular as . If the coordinate function and the Jacobian are explicitly inserted into the equations, using their form given by Eq. (25) and Eq. (26), respectively, divergent terms are seen to cancel. However, since we do not explicitly impose a specific compactification — we have formulated the problem in terms of generic Jacobians rather than the specific formulas of Eqs. (25) and (26) — we need to be careful to avoid irregular terms at .

It is sufficient to require that in the limit , the compactified coordinate transformation and its Jacobian approach the explicit forms given in Eq. (26) and Eq. (25), respectively. For instance, consider the equation

 U,x=(r2(x)dxdr)−1\rme2β(Q+r2e−2βNU). (27)

According to Eq. (26) and Eq. (25), we have

 (r2(x)dxdr)−1→r−1Γ (28)

as . Hence, at ,

 U,x=r−1Γ\rme2β(Q+r2e−2βNU)=r−1Γ\rme2β(KQ−J¯Q). (29)

We proceed in a similar manner for the other hypersurface equations. The specific form of the equations at is given in B. Note additionally that for , and it is possible to directly evaluate the respective quantity without radial integration at . For instance, as ,

 Q=−2ðβ. (30)

## 4 Numerical methods

### 4.1 Discrete representation of the evolution variables

The evolution algorithm is a hybrid of finite-difference (for radial and time integration) and pseudo-spectral (for angular directions) methods. In the compactified radial direction , fields are evaluated on a uniform grid of points, , at points . The inner coordinate radius, , is that of the world-tube , where we need to specify appropriate boundary data at any given time to carry out a radially outward hypersurface integration. As we will see in Section 4.4, our radial integration scheme actually requires radial points to start the algorithm. We therefore need to provide boundary data on the first radial points so that our worldtube spans three radial points. Boundary data is required for

 {β,Q,U,^W,Φ}for allyAi=0,1,2∈Γ|u. (31)

The outer boundary of the compactified radial grid is placed at the outermost gridpoint corresponding to future null infinity .

At each radial point , we represent angular dependence as a spectral expansion in terms of real-valued spin-weighted spherical harmonics, according to

 Φi(yA)=∞∑ℓ=sm=+ℓ∑m=−ℓΦℓm(xi)sZℓm(yA), (32)

where the are spin real-valued spherical harmonics, defined in terms of the standard basis [22] by

 sZℓm=⎧⎪ ⎪⎨⎪ ⎪⎩1√2(sYℓm+(−1)msYℓ−m),m>0\,,sYℓm,m=0\,,i√2((−1)msYℓm−sYℓ−m),m<0\,. (33)

We use the to accommodate the real-valued spin-0 quantities, which naturally yield real valued coefficients.

We store a finite number of harmonic coefficients for each variable, terminating the sum according to the maximum number of measurable gravitational wave modes contained in the solution. For our test case with linearized solutions as discussed in Section 5, this is . For the realistic case of current binary black hole merger simulations the number of resolved gravitational-wave modes in the 3+1 evolution is typically , beyond which their amplitude is below the level of discretization error. Although this case is not considered here and is left for future work, we do present results of a stability test in which . We store the spectral coefficients of the expansion of each evolution variable at each point of the radial grid.

Radial and time integration is performed entirely on the spherical harmonic coefficients of the evolution variables, with the one exception being Eq. (56), to be introduced below. This equation contains non-linear terms of the form where and are both functions of angular coordinates . To compute non-linear terms occurring either in Eq. (56) or in the RHS for a given hypersurface equation, we first need to recompose the involved variables on a collocation grid. After the terms have been evaluated on the collocation grid, we decompose them back into real-valued spin-weighted spherical harmonics.

We construct a pseudo-spectral collocation grid for spherical harmonics by defining a set of grid points on a spherical shell using a spherical-polar coordinate system with constant grid spacing in and direction

 ⎧⎨⎩(θj,ϕk)=⎛⎝πj+12Nθ,2πkNϕ⎞⎠:j,k∈N;0≤j

Recomposing quantities on the collocation grid is easily done by evaluating the sum of the spherical harmonic expansion via

 f(θj,ϕk)=∑ℓmfℓmsZℓm(θj,ϕk). (35)

Decomposing a quantity into spherical harmonics requires surface integration over . The expansion coefficients are computed according to a discrete version of

 fℓm=∫ΩdΩf(θ,ϕ)s¯Zℓm(θ,ϕ), (36)

where is the surface element on the collocation grid. A numerical integration algorithm which is exact for spherical harmonics up to order is given by Gauss-Chebyshev quadratures using points on ( e.g. [23]). This algorithm makes use of coordinate dependent weights in direction. In , the weights are simply since is a periodic coordinate direction. Since we use equally spaced points in (equivalent to Chebyshev nodes in ), the weights in direction are [24]

 wj=4NθNθ/2−1∑ℓ=012ℓ+1sin((2ℓ+1)θj). (37)

The integral Eq. (36) reduces to

 fℓm=πNθ2πNϕNθ∑j=0Nϕ∑k=0f(θj,ϕk)s¯Zℓm(θj,ϕk)sinθjwj, (38)

provided we have

 (Nθ,Nϕ)≤(2(ℓ+1),2(ℓ+1)) (39)

points on .

We have thus established an exact mapping between the representation in terms of spherical harmonic coefficients and the representation on the collocation grid. To speed up the computation, we precompute the , and the product .

As an example for our algorithm, consider the linearized version of the hypersurface equation for ,

 U,x=(r2(x)dxdr)−1\rme2βQ. (40)

We recompose and on the sphere using their spectral expansion coefficients in order to carry out the required multiplications. The complete procedure for integrating Eq. (40) can be summarized as follows.

1. Define grid points on a spherical shell according to Eq. (34).

2. Recompose

 β(xi,θj,ϕk)=∑ℓmβℓm(xi)0Zℓm(θj,ϕk), (41)

and

 Q(xi,θj,ϕk)=∑ℓmQℓm(xi)1Zℓm(θj,ϕk). (42)
3. Loop over all angular points and . For each angular point, compute using  Eq. (40).

4. Decompose via

 (Uℓm),x(xi)=πNθ2πNϕNθ∑j=0Nϕ∑k=0U,x(θj,ϕk)s¯Zℓm(θj,ϕk)sinθjwj. (43)
2. Radially integrate to obtain .

The last step, radial integration, is described in more detail in Section 4.4.

### 4.2 Radial derivatives and dissipation

Radial derivatives of all hypersurface quantities are generally obtained from the RHS of their corresponding radial ODE integrations and hence do not need to be recomputed by means of finite difference operators. However, the metric variable (and also ) itself is not directly obtained via radial integration and hence must be computed everywhere. We approximate and by means of finite difference operators of 4th-order. The radial derivative of can be obtained by using Eq. (15).

According to stability analysis and empirical findings (see Section 2), we apply fully side-winded derivatives with the stencil points in the direction of . We use 4th-order first and second derivatives

 ∂fi = 1Δx(−2512fi+4fi+1−3fi+2+43fi+3−14fi+4), (44) ∂2fi = 1Δx2(154fi−776fi+1+1076fi+2−13fi+3+6112fi+4−56fi+5), (45)

where is the grid spacing in the compactified radial coordinate direction.

Close to when , we switch to 4th-order centred stencils

 ∂fi = 1Δx(+112fi−2−23fi−2+23fi+1−112fi+2), (46) ∂2fi = 1Δx2(−112fi−2+43fi−1−52fi+43fi+1−112fi+2), (47)

and when , we switch to side-winded stencils pointing towards the inner boundary

 ∂fi = 1Δx(2512fi−4fi−1+3fi−2−43fi−3+14fi−4), (48) ∂2fi = 1Δx2(154fi−776fi−1+1076fi−2−13fi−3+6112fi−4−56fi−5). (49)

In addition, we apply a numerical dissipation operator to . We use a 5th-order Kreiss-Oliger dissipation operator of the form

 Dfi=ϵdiss64Δx(fi−3−6fi−2+15fi−1−20fi+15fi+1−6fi+2+fi+3), (50)

where controls the strength of the applied dissipation operator . At the outer boundary (at ), where we do not have enough points to compute the dissipation operator, we use one-sided stencil derived for an overall 4th-order accurate summation-by-parts (SBP) operator (though we do not make explicit use of the SBP property). The particular stencil coefficients are derived in [25]. We explicitly state the stencil coefficients in C.

Empirical tests have shown that radial dissipation applied to is crucial to improve the stability properties of our scheme.

### 4.3 Angular derivatives

Numerical derivatives in the angular direction are obtained via analytic angular derivatives of the spin-weighted real-valued spherical harmonic spectral basis functions. The action of the derivative on the real-valued spherical harmonics is given by [22]

 ðsZℓm=√(ℓ+s+1)(ℓ−s)s+1Zℓm. (51)

The action of on a spectrally expanded function is

 ðf=∑ℓmfℓm√(ℓ+s+1)(ℓ−s)s+1Zℓm. (52)

Similarly, the action of is given by

 ¯ðsZℓm=−√(ℓ−s+1)(ℓ+s)s−1Zℓm. (53)

The hypersurface equations are integrated in the radial direction using a multistep method. The classes of methods that we have studied for this problem are either the explicit Adams-Bashforth methods, the implicit Adams-Moulton methods, and a combination of the two in the form of a predictor-corrector scheme. However, as discussed in Section 2, at 4th-order, explicit methods are unstable for our particular set of equations.

Thus, the radial integration uses a fully implicit method. Fortunately, the Einstein equations in Bondi-Sachs form are particularly convenient for this purpose, as they form a hierarchy (as outlined in Section 3) with only one unknown function in each equation. Furthermore, the equations are linear in this unknown. Schematically, we write each equation in the form

 dydx+yg(x)=f(x), (54)

and the 4th-order fully implicit Adams-Moulton scheme can be written in explicit form as

 yi+1 (1+38hgi+1)=yi+h(38fi+1+1924(fi−yigi) (55) −524(fi−1−yi−1gi−1)+124(fi−2−yi−2gi−2)).

We again note that the quantities we work with are the spherical harmonic coefficients of each variable, which are which are functions of the compactified radius , according to the procedure outlined at the end of Section 4.1.

The radial integration requires special treatment due to the nonlinear term, in Eq. (21) (given explicitly in A). It is a function of both and , and thus, according to Eq. (22), both and . We write Eq. (21) in the form

 Φ,x+Φa+¯Φb=RΦ, (56)

where

 a = (r(x)dxdr)−1+J2K(¯JK,x−¯J,xK) b = J2K(JK,x−J,xK). (57)

and contains the remaining terms but does not involve or or their derivatives. This is integrated using the scheme of Eq. (55) to obtain an equation that involves both and . Taking the complex conjugate leads to a second equation in the two unknowns, and solving the system gives

 Φi+1 ([1+3h8¯ai+1][1+3h8ai+1]−(3h8)2¯bi+1bi+1) (58) =Ti(1+3h8¯ai+1)−¯Ti3h8bi+1,

where

 Ti=Φi + 3h8RΦ,i+1 (59) + 19h24(RΦ,i−Φiai−¯Φibi) − 5h24(RΦ,i−1−Φi−1ai−1−¯Φi−1bi−1) + h24(RΦ,i−2−Φi−2ai−2−¯Φi−2bi−2).

The scheme above does not allow us to work directly with the spherical harmonic coefficients of , , , , , and , due to the non-linear terms and which must be evaluated on the collocation grid. To compute Eq. (59), we therefore recompose , , , , and on the collocation grid defined by Eq. (34) to perform the required multiplications. Having evaluated according to Eq. (59) on the collocation grid, we decompose to obtain its spectral coefficients in terms of spin real-valued spherical harmonics for each radial point .

Note that the radial integration schemes Eq. (55) and Eq. (59) both require data on radial points to start the algorithm. These must be provided as boundary data on the worldtube .

### 4.5 Time integration

The evolution equation for has the form

 J,u=Φ. (60)

This equation can be straightforwardly integrated via a 4th-order Runge-Kutta scheme using the spectral coefficients of . In addition, we add numerical dissipation to Eq. (60). To be explicit, we solve

 (Jℓmi),u=Φℓmi+DJℓmi,∀ℓ,m,i, (61)

where is a dissipation operator defined in Section 4.2. Since it is necessary to solve the hypersurface equations to obtain the , the hypersurface equations must be solved for each intermediate Runge-Kutta step.

### 4.6 Summary of algorithm

1. Assume data for in the form of spectral coefficients at (intermediate) timestep for each radial shell . If is the first timestep, the are given by initial data.

2. Compute , as well as radial derivatives and from by means of Eq. (44) and Eq. (45), respectively. Since is related non-linearly to , we need to recompose from to evaluate on the collocation grid. Afterwards, we decompose to obtain .

3. Provide inner worldtube boundary data for , , and at (intermediate) timestep in terms of spectral coefficients on the first radial points.

4. Integrate hypersurface equations in the order (i) , (ii) , (iii) , (iv) , and (v) by using the steps described in Section 4 and Section 4.4.

5. Evaluate next Runge-Kutta step for to obtain at next (intermediate) step as described in Section 4.5.

### 4.7 Remarks on the computational implementation

We have implemented a new code within the Cactus computational toolkit [26, 27]. The underlying grid array structures are provided by Carpet [28, 29]. Memory handling of the collocation grid, and recomposition/decomposition in terms of spin-weighted spherical harmonics is provided by SphericalSlice [30].

Since the memory consumption of the implemented code is rather low, and since the computational efficiency is high, we do not currently decompose the domain to distribute the work load across multiple processing units. We do, however, make use of multi-threading via OpenMP to enable faster processing on shared memory units. In particular, we use multi-threading in the following two kinds of loops: (i) when looping over spectral coefficients to perform, for instance, radial integration for each separate mode, and (ii) when looping over points on the collocation grid to evaluate non-linear terms in the equations (such as RHS evaluation). Depending on the number of spectral coefficients and the number of available cores within one shared memory unit, the observed scaling can be close to the optimum (though we remark that we have done a rather limited number of tests on a compute node with up to shared memory cores).

The implemented code is designed such that it can run concurrently with our 3+1 evolution code Llama [31, 32]. This is important for future application in on-the-fly Cauchy-characteristic extraction, where the metric data will be transported to during Cauchy evolution without the need of an additional post-processing step (which is currently the case for the algorithm presented in [2, 3]). Furthermore, this is necessary for a future implementation of Cauchy characteristic matching[33, 7] in which the characteristic evolution is used to provide on-the-fly boundary data for a Cauchy evolution.

## 5 Results

### 5.1 Linearized solutions

To test the convergence of our numerical scheme on a dynamical spacetime, we use solutions to the linearized Einstein equations in Bondi-Sachs form on a Minkowski background (Section 4.3 of  [34]). We write

 Jlin = √(ℓ−1)ℓ(ℓ+1)(ℓ+2)2ZℓmR(Jℓ(r)eiνu), Ulin = √ℓ(ℓ+1)1ZℓmR(Uℓ(r)eiνu), βlin = ZℓmR(βℓeiνu), ^Wlin = ZℓmR(^Wℓ(r)eiνu), (62)

where , , , are in general complex, and taking the real part leads to and terms. The quantities and are real; while and are complex due to the terms and , representing different terms in the angular part of the metric. We require a solution that is well-behaved at future null infinity. We find [15], in the case ,

 β2 = β0, J2(r) = 24β0+3iνC1−iν3C236+C14r−C212r3, U2(r) = −24iνβ0+3ν2C1−ν4C236+2β0r+C12r2+iνC23r3+C24r4, ^W2(r) = 24iνβ0−3ν2C1+ν4C26+3iνC1−6β0−iν3C23r (63) −ν2C2r2+iνC2r3+C22r4,

with the (complex) constants , and freely specifiable; and in the case

 β3 = β0, J3(r) = 60β0+3iνC1+ν4C2180+C110r−iνC26r3−C24r4, U3(r) = −60iνβ0+3ν2C1−iν5C2180 +2β0r+C12r2−2ν2C23r3+5iνC24r4+C2r5, ^W3(r) = 60iνβ0−3ν2C1+iν5C215+iνC1−2β0+ν4C23r (64) −i2ν3C2r2−4iν2C2r3+5νC2r4+3C2r5.

We establish convergence by testing the evolution quantities against linearized solutions listed in Section 5.1. The linearized solution provides initial data for on a null cone, as well as boundary data for , , , , and at the worldtube . During evolution, we compute the error in all evolved quantities by comparing with the linearized solution. Since the code solves the general nonlinear case whereas the exact solution satisfies the linearized Einstein equations, we expect to converge towards zero at the order of accuracy of the numerical scheme only in a regime in which is much larger than any nonlinear contribution.

We have performed a number of test cases using linearized solutions, linearized solutions, and a superposition of both. In the latter case, we compute via

 Jlin=∑ℓ,mκℓ2ZℓmR(Jℓ(r)\rmeiνu), (65)

where . The remaining superposed linearized solutions for all other quantities are constructed in the same way, using the appropriate spin weight for the and the appropriate -dependent coefficients , respectively (compare Eq. (62)).

The linearized solutions depend on free parameters , , and which we have tested for a range of different values. Note that the amplitudes , , and must be of linear order ().

In all cases considered, we find better than 4th-order convergence until the error roughly reaches the square of the amplitude of the linearized solution, beyond which convergence deteriorates, as expected, due to the emergence of nonlinear behavior.

As a particular example, we show convergence of a superposed solution with parameters

 C1 =3×10−6,C2 =1×10−6, (66) β0 =i×10−6,ν =1.0. (67)

Fig. 1 plots the -norm of the error on two resolutions and (see Table 1) scaled for 4th-order convergence. We define the -norm in terms of the sum over all modes and radial points by

 ∥f∥=√∑i,ℓ,m(fℓm(xi))2. (68)

The appropriate convergence scaling can be determined from the convergence rate defined in terms of the grid spacing ,

 C=(Δxr0Δxr1)p, (69)

where is the expected order of convergence. By doubling the resolution, we expect the higher resolution error, , to be smaller by a factor of given our 4th-order accurate algorithm, i.e., we should get

 C=(∥ϵ(J)∥r0∥ϵ(J)∥r1)p=4=16. (70)

As shown in Fig. 1, this is indeed the case. We measure better than 4th-order convergence (see further below for a discussion). Furthermore, the evolution is still stable after (corresponding to cycles of the solution).

The grid settings and parameters for this test are given in Table 1. In all cases, we apply radial dissipation of amplitude . The inner boundary is located at