C Reisswig

N T Bishop

D Pollney

B Szilágyi

Characteristic extraction in numerical relativity: binary black hole merger waveforms at null infinity

Keywords: numbers: 04.25.dg, 04.30.Db, 04.20.Ha, 04.30.Nk 1

published or not. The documents may come

I. INTRODUCTION

The computation of gravitational radiation from black hole merger events has attracted considerable attention, since the pioneering work by Smarr and collaborators [START_REF] Smarr | [END_REF][START_REF] Smarr | [END_REF][3]. With the advent of ground-based laser interferometric gravitational wave detectors, as well as the prospect of the Laser Interferometer Space Antenna (LISA), interest in the problem has increased considerably in recent years. The measurement of gravitational waves is expected to provide an important probe of strong-field nonlinear gravity, thereby enabling empirical testing of this regime, as well as to provide direct observations of processes of fundamental interest in astrophysics. The sensitivities of LISA, and of the upcoming advanced ground-based detectors AdLIGO and AdVirgo, are high enough that even an error in aspects of the waveform calculation of 0.1% could lead to an incorrect interpretation of the astrophysical properties of a source, or of a test of general relativity [4].

Nowadays, there are several codes that can produce a stable and convergent simulation of a black hole spacetime [5][6][7][8][9][10][11][12][13][14], and of particular topical relevance is the construction of long and accurate waveforms which to be used to analyze distinguishable features within the signals [15], and also to construct a family of templates with the goal of improving gravitational wave detection algorithms. Here the requirements are particularly challenging for numerical simulations, requiring waveforms which are accurate in phase and amplitude over many cycles to allow for an unambiguous matching to post-Newtonian waveforms at large separation. Some recent studies have shown very promising results in this direction for particular binary black hole models [16][17][18][19][20][21][22][23][24][25][26][27][28].

However, a particular difficulty arises from the fact that gravitational energy cannot be defined locally in general relativity, is well-defined only at future null infinity, J + . Since standard numerical evolutions are carried out on finite domains, there is a systematic error caused by estimating the gravitational radiation from fields on a world-tube at finite radius and the uncertainty in how it relates to measurement at J + [29]. Even though this error may be small, the expected sensitivity of upcoming detectors means that it is important to obtain an accurate result.

A rigorous formalism for the global measurement of gravitational energy at J + has been in place since the pioneering work of Bondi, Penrose and collaborators in the 1960s [30,31]; and subsequently, techniques for calculating gravitational radiation at J + have been developed. The idea which we pursue here is to combine a standard Cauchy or "3 + 1" numerical relativity code with a characteristic code [32] in order to transport the wave information to J + using the full Einstein equations. This method involves a Cauchy computation within a domain bounded by a world-tube, Γ, together with a characteristic computation with inner boundary Γ. If the data is consistently passed back and forth between the evolution schemes at Γ during the evolution, the method is known as Cauchy-characteristic matching (CCM). Given astrophysical initial data, such a method has only discretisation error [33], and a complete mathematical specification has been developed [START_REF] Bishop | Black Holes, Gravitational Radiation and the Universe[END_REF]. However, efforts to implement CCM encountered stability problems at the interface, and in the general case have not been successful. An alternative approach is Cauchycharacteristic extraction (CCE), or characteristic extraction. In this case, an independent Cauchy computation is carried initially, with its usual timelike artificial outer boundary and knows nothing about a characteristic code. Data on the world-tube Γ is stored, and then used to provide inner boundary data for a characteristic code -but not vice-versa. Initial implementations of CCE [START_REF] Babiuc | [END_REF]36] have shown the effectiveness of the method for a number of test problems.

In this paper we describe the implementation of a CCE code, and present results obtained from the code for the real astrophysical problem of the inspiral and merger waveform of two black holes.

The gravitational waveforms are calculated at J + , and are thus, for this case, the first waveforms that are unambiguous in the sense of being free of gauge and finite-radius effects. Further, the code is general purpose, in that it is independent of the details of the Cauchy code: it runs as a post-processing operation, taking as input only the required geometrical data on a world-tube.

Thus it can be used in combination with other Cauchy codes, and applied to other astrophysical problems.

The next section summarizes notation and results needed from other work. Then, Sec. III describes in some detail our implementation of characteristic extraction. Sec. IV presents results of testing the code against analytic solutions. The computation, including results, of the inspiral and merger of two black holes is described in Sec. V. Some technical details of the characteristic and extraction code parameters are described in Sec. VI.

II. REVIEW OF RESULTS NEEDED FROM OTHER WORK

A. The Bondi-Sachs metric

The formalism for the numerical evolution of Einstein's equations, in null cone coordinates, is well known [30,33,[37][38][39][40]. For the sake of completeness, we give a summary of those aspects of the formalism that will be used here. We start with coordinates based upon a family of outgoing null hypersurfaces. We let u label these hypersurfaces, x A (A = 2, 3), label the null rays, and r is a surface area coordinate. In the resulting x α = (u, r, x A ) coordinates, the metric takes the Bondi-Sachs form [30,41]

ds 2 = -e 2β (1 + W c r) -r 2 h AB U A U B du 2 -2e 2β dudr -2r 2 h AB U B dudx A + r 2 h AB dx A dx B , ( 1
)
where h AB h BC = δ A C and det(h AB ) = det(q AB ), with q AB a metric representing a unit 2-sphere embedded in flat Euclidean 3-space; W c is a normalized variable used in the code, related to the usual Bondi-Sachs variable V by V = r + W c r 2 . As discussed in more detail below, we represent q AB by means of a complex dyad q A . Then, for an arbitrary Bondi-Sachs metric, h AB can be represented by its dyad component

J = h AB q A q B /2, ( 2
)
with the spherically symmetric case characterized by J = 0. We also introduce the fields

K = 1 + J J, U = U A q A , (3)
as well as the (complex differential) eth operators ð and ð [42].

The 10 Einstein equations R αβ = 8π(T αβ -1 2 g αβ T ) are classified as: (i) hypersurface equations -R 11 , q A R 1A , h AB R AB -forming a hierarchical set for β, U and W c ; (ii) evolution equation q A q B R AB for J; and (iii) constraints R 0α . An evolution problem is normally formulated in the region of spacetime between a timelike or null world-tube, Γ, and future null infinity (J + ), with (free) initial data, J, given on u = 0, and with boundary data for β, U, W c , J satisfying the constraints given on the inner world-tube. The inclusion of J + in the computational grid is made possible by compactifying the radial coordinate, r, for instance by means of a transformation of the form

r → x = r r + r wt . ( 4
)
In characteristic coordinates (but not in general), the Einstein equations remain regular at J + under such a transformation.

B. The spin-weighted formalism and the ð operator

A complex dyad q A is a 2-vector whose real and imaginary parts are unit vectors that are orthogonal to each other. Further, q A represents the metric, and has the properties

q A q A = 0, q A qA = 2, q AB = 1 2 (q A qB + qA q B ). (5)
Note that q A is not unique, up to a unitary factor: if q A represents a given 2-metric, then so does q A = e iα q A . Thus, considerations of simplicity are used in deciding the precise form of dyad to represent a particular 2-metric. The dyad commonly used to represent the unit sphere metric in the stereographic coordinates used here, is

ds 2 = 4(dq 2 + dp 2 ) (1 + q 2 + p 2 ) 2 , q A = 2 1 + q 2 + p 2 (1, i). ( 6
)
Having defined a dyad, we may construct complex quantities representing all manner of tensorial objects, for example

X 1 = T A q A , X 2 = T AB q A qB , X 3 = T AB C qA qB qC . Each object has
no free indices, and has associated with it a spin-weight s defined as the number of q factors less the number of q factors in its definition. For example, s(X 1 ) = 1, s(X 2 ) = 0, s(X 3 ) = -3, and, in general, s(X) = -s( X). We define derivative operators ð and ð acting on a quantity V with spin-weight s

ðV = q A ∂ A V + sΓV, ðV = qA ∂ A V -s ΓV (7)
where the spin-weights of ðV and ðV are s + 1 and s -1, respectively, and where

Γ = - 1 2 q A qB ∇ A q B . ( 8
)
In the case of stereographic coordinates, Γ = q + ip.

The spin-weights of the quantities used in the Bondi-Sachs metric are

s(W c ) = s(β) = 0, s(J) = 2, s( J) = -2, s(K) = 0, s(U ) = 1, s( Ū ) = -1. ( 9
)
We will be using spin-weighted spherical harmonics [43,44] s Y m (the suffix s denotes the spin-weight), using the formalism described in [45], and note that this reference gives explicit formulas for the s Y m in stereographic coordinates. In the case s = 0, the s will be omitted, i.e.

Y m = 0 Y m . Note that the effect of the ð operator acting on Y m is

ðY m = ( + 1) 1 Y m , ð 2 Y m = ( -1) ( + 1)( + 2) 2 Y m . ( 10
)
The mathematical theory relating metric quantities at future null infinity, J + , to gravitational radiation, and using the present formalism, is given in [START_REF] Babiuc | [END_REF]38,46]. In the original work of Bondi et al. [30], it was possible to assume that the coordinates at J + were such that β = J = 0 there, and in that case the gravitational news takes the very simple form

N = 1 2 ∂ u ∂ ρ J (11)
(where we have written the radial derivative in terms of ρ := r -1 rather than the usual notation, , to avoid confusion with the spherical harmonic index ). However, the coordinates used in the characteristic code are fixed at the inner boundary Γ, and in general the Bondi gauge conditions are not satisfied. Previous work has presented the formalism for calculating, in a general gauge, the gravitational news [38], as well as ψ 4 [START_REF] Babiuc | [END_REF].

Geometrically, the Bondi gauge condition J = 0 means that, for large r, a 2-surface of constant (u, r) is spherical rather than being, for example, an ellipsoid; and this condition can be expressed algebraically by saying that the 2-surface has constant curvature. An expression for the news, in a general gauge, must take (implicit) account of the transformation to Bondi gauge coordinates,

r → r [B] = ω(u, x A )r, x A → x A [B] = x A [B] (u, x A ). (12)
The constant curvature condition leads to the factor ω satisfying

2K -ð ðK + 1 2 ð2 J + ð 2 J + 1 4K ð JðJ -ðJð J = 2 ω 2 + h AB D A D B log ω ( 13
)
with all metric quantities evaluated at J + . Alternatively, we can initialize J = 0 at u = 0, so that ω = 1 at u = 0, and then evolve ω by

∂ u ω = - 1 2 ( Ū ðω + U ðω) - ω 4 (ð Ū + ðU ). (14)
The Bondi gauge condition β = 0 at J + means that, for large r, coordinate time is the same as proper time, and the implementation of this condition is straightforward.

While the transformation to Bondi gauge coordinates can be done explicitly [46], in the present code it is done implicitly, leading to an expression for the news in terms of the code metric variables and coordinates. The formula is long and complicated [38] and is not repeated here. However, in the linearized case, the news function is much simpler [47]. The results presented later are approximately in the linearized regime, and so it can be instructive to examine this case (even though the code is fully nonlinear). The linearized formula for the news is

N = 1 2 ∂ u ∂ ρ J + 1 2 ð 2 (ω + 2β). (15)
In this case ω is simply related to J at J + , if the metric is decomposed into spin-weighted spherical harmonics,

J = m J m 2 Y m , β = m β m Y m , ( 16
)
Then the news is

N = m N m 2 Y m = m 1 2 ∂ u ∂ ρ J m - ( + 1) 4 J m + β m ( -1) ( + 1)( + 2) 2 Y m . ( 17
)
In characteristic work, it is conventional to work with a quantity Ψ 4 rather than the usual ψ 4 used in current Cauchy codes (see e.g. [9]). The relationship is

Ψ 4 = -(1/2) ψ4 , (18)
and we will usually transform the characteristic Ψ 4 to the standard ψ 4 . The formulas for Ψ 4 are long and complicated and are not given here, but in the linearized case, Ψ 4 can be expressed in a simple way in terms of metric quantities at J + [START_REF] Babiuc | [END_REF] (in which the notation Ψ rather than Ψ 4 is used)

Ψ 4 = 1 2 ∂ 2 u ∂ ρ J - 1 2 ∂ u J - 1 2 ðU - 1 8 ð 2 (ð Ū + ðU ) + ∂ u ð 2 β. ( 19
)
It can be shown that the linearized Einstein equations lead to the relation Ψ 4 = ∂ u N .

D. Characteristic extraction

We refer to [START_REF] Bishop | Black Holes, Gravitational Radiation and the Universe[END_REF] (see also [36,48]) for a full description of the characteristic extraction process.

While characteristic extraction appears to be only a coordinate transformation, it is actually rather more complicated. The difficulty is that Bondi-Sachs coordinates use a surface area radial coordinate, and this coordinate can be constructed only once the angular coefficients of the metric have been found. Thus the construction proceeds in two stages.

In the first stage, we use an affine coordinate λ in the radial direction, and find the transformed metric and its first ,λ derivatives at angular grid-points of the extraction world-tube Γ. The process is illustrated in Fig. 1, and can be summarized as follows: • Define a world-tube Γ by x 2 +y 2 +z 2 = R 2 with R constant, and induce angular coordinates φ A on Γ as though in Euclidean space.

s α x α = x α ( ; , φ λ u ) A N(u) S : 2 + y 2 + z 2 = R x α t φ A Γ H(t)
• Let H be a hypersurface of constant t, and define the 2-surface S = H ∩ Γ.

• Let t α be a unit normal to H, and let s α be normal to S in H.

• Construct outgoing null vectors α = t α + s α , and then outgoing null geodesics in direction α with affine parameter λ.

• The union of such geodesics is the null cone N , labeled by u being the Cauchy time t where N meets Γ.

• Construct the Jacobian for the coordinate transformation (t, x, y, z) → (u, λ, φ A ).

• Find the transformed metric at the angular grid-points of Γ.

• Find the first ,λ derivatives of the transformed metric at the angular grid-points of Γ.

In the second stage, we make the transformation to a surface area coordinate r. The difficulty here is that, in general, r is not constant on Γ. Thus in order to set data on an inner world-tube of the characteristic grid, we need the metric quantities on the world-tube, as well as their first derivatives off the world-tube. The process can be summarized as follows:

• Make the coordinate transformation λ → r = r(u, λ, φ A ) with r defined by the condition that it is a surface area coordinate.

• Find r and r ,λ at the grid points of Γ.

• Find the metric and its first ,λ derivatives at the grid points of Γ.

• Find J, U , W c and β at the grid points of Γ.

• Using ∂ r = ∂ λ /r ,λ , find J ,r , U ,r , W c,r and β ,r at the grid points of Γ.

• Set J, U , U ,r , W c and β on an inner world-tube of the characteristic grid.

The value used for U ,r is second-order accurate on the world-tube Γ, and it is this value which is used at the inner world-tube of the characteristic grid, so, in general, it is only first-order accurate.

However, in practice, the overall performance of the code remains second-order accurate.

E. The Cauchy evolution code and finite radius measurements

A fundamental problem with early attempts to implement CCM was the sensitivity of the 3+1 (Cauchy) implementation to boundary perturbations. In recent years, however, a number of strongly hyperbolic formulations have been developed and proven to be numerically robust for binary black hole evolutions using artificial outer boundary conditions We evolve the spacetime in the neighborhood of the binary black holes using a variant of the "BSSNOK" evolution system [49][50][51], according to the implementation described in [9]. This formulation uses a conformally transformed metric, γab , and (traceless) extrinsic curvature, Ãab , as evolution variables, where the conformal transformation is defined by

φ = 1 12 ln det γ ab , γab := e -4φ γ ab . ( 20
)
Instead of evolving φ directly, we use the concomitant

φ6 := (det γ ab ) -1/6 = φ -1/2 , ( 21
)
following the suggestion of [52].

Initial data for binary black holes is determined using the puncture method [53], and corresponds to a conformally flat solution of the Hamiltonian constraint with the extrinsic curvature specified according to Bowen and York [54]. Parameters for quasi-circular orbits starting from a given separation are determined by a post-Newtonian estimate [14], resulting in a low eccentricity (e 0.04) inspiral for the equal-mass non-spinning binary considered below.

The lapse and shift are evolved according to the "1+log" [55] and " Γ-driver" [56], respectively, including shift-advection terms which have proven to be important for maintaining phase accuracy in the puncture motion [7,57]. This combination of gauges are key ingredients of the "moving puncture" method [5,7] and the control of evolution of the coordinate singularity within the black hole throat [58].

The fields are discretize on regular locally Cartesian grids which cover the neighborhood of the binary to some finite coordinate radius. Finite differences at 8th-order are used to numerically approximate derivatives. Boundaries between coordinate patches are determined by 5th-order Lagrange polynomial interpolation. Fields are evolved in time according to the BSSNOK prescription, using a 4th-order Runge-Kutta integrator. A detailed description of the numerical scheme can be found in [9].

An objective of this paper is to compare the gravitational waveforms measured at J + with the locally determined measures which are in common use in numerical relativity. These local measurements generally evaluate the geometrical fields (the metric and associated derivatives) on a coordinate sphere of constant radius. The measurement spheres are assumed to be in the "wavezone", which is commonly taken to be from r = 18M [59] to r = 225M [8] from the source. The particular numerical grid structure employed in this paper follows that of companion studies in which we have demonstrated accurate and convergent wave measurements to r = 1000M [9,60].

Polynomial extrapolations to spatial r → ∞ are applied in order to reduce the systematic error associated with the finite radius measurements. Below, we demonstrate that the error associated with this procedure in comparison with the CCE result is, however, still larger than the numerical error for the models we have studied.

Two methods are used in order to estimate the gravitational radiation at a given radius. The first evaluates the Weyl tensor in a particular null frame (n, l, m, m) [31], which is chosen to be oriented along the coordinate radial direction, and orthonormalised using the local metric. According to the peeling property of asymptotically flat spacetimes, the component

ψ 4 := C n mn m ( 22
)
is associated with the outgoing gravitational radiation. As an alternative to ψ 4 , we also measure waves by evaluating the Zerilli-Moncrief variables, assuming the spacetime metric at large radius from the source to be representable as perturbations of a fixed Schwarzschild background [61][62][63]. We evaluate 1st-order gauge invariant odd-parity and even parity multipoles, (Q × lm and Q + lm , respectively) [64,65]. The implementation of the ψ 4 and Zerilli-Moncrief measurements follow the outline described in [9,10].

Our evolution code is built using the Cactus computational framework [START_REF] Goodale | Vector and Parallel Processing -VECPAR'2002, 5th International Conference[END_REF][START_REF]Cactus Computational Toolkit home page[END_REF]. The grid structures are implemented using extensions to the Carpet mesh refinement driver [START_REF] Schnetter | [END_REF][69][START_REF]Mesh Refinement with Carpet[END_REF], allowing the evolution domain to be covered by multiple patches with generalized local coordinates [9]. Initial data is computed using the TwoPunctures code [START_REF] Ansorg | [END_REF]), which have been used

for previous Cartesian evolutions.

III. IMPLEMENTATION OF CHARACTERISTIC EXTRACTION AS POST-PROCESSING

CODE

The characteristic extraction is independently performed as a post-processing step after the full Cauchy evolution. During the Cauchy evolution, we output the Cauchy 4-metric and its derivatives. The variables are decomposed in terms of spherical harmonic amplitudes on a set of extraction spheres with fixed radii around the origin of the black hole spacetime. On the characteristic side, this data can then be read and reconstructed to generate inner boundary data on the worldtube.

Since we obtain the fields at J + in terms of coordinate time u, we finally have to transform the fields to constant Bondi time u B . Furthermore, due to the rather complicated procedure of calculating the conformal factor ω for the calculation of the news, we have implemented a linearized approximation that turns out to be numerically more accurate than the full non-linear computation.

In the next subsections, we will describe these procedures in more detail.

A. World-tube boundary data

Boundary data for the characteristic evolution are constructed during Cauchy evolution in terms of the ADM variables on coordinate spheres with fixed radius R. The time series of the coordinate spheres at successive times generates the discretize world-tube.

In order to be able to transform the ADM variables to a representation of the full 4-metric in terms of Bondi coordinates and in terms of the characteristic evolution variables in a neighborhood of the world-tube, we have to compute the following quantities: the ADM metric components γ ij , as well as the time-derivatives ∂ t γ ij and Cartesian derivatives ∂ k γ ij , the lapse α, as well as time and Cartesian derivatives ∂ t α and ∂ i α, the shift vector β i , and again time and Cartesian derivatives

∂ t β i and ∂ i β k .
In the specific formulation of the evolution system that is used in this paper, we can compute these quantities as follows. The BSSNOK evolution variables are transformed to the ADM variables after each evolution step and are therefore known everywhere. What remains is the calculation of time and spatial derivatives.

The time derivatives of the lapse and shift vector are known directly from the RHS of the "1+log" slicing condition and the hyperbolic Γ-driver condition (see [9] for details of the implementation of these gauge conditions). The RHSs of these equations have to be computed at each time-step in order to evolve the gauge variables with a method of lines (MoL) scheme, similar to the evolution equations.

The time derivatives of the metric components can be obtained via the ADM relation

∂ t γ ij = -2αK ij + D i β j + D j β i , ( 23
)
where D i is the covariant derivative operator compatible with the spatial 3-metric γ ij and K ij the extrinsic curvature. In order to compute the covariant derivative, it is necessary to invert the 3-metric and to calculate the Christoffel symbols on each grid point. Because the shift vector is given in contravariant form, we have to transform it to its covariant representation with index down. After looping over all points on the computational grid, we store the outcome in grid functions, i.e., one per metric component.

Next, obtain the Cartesian derivatives of the metric, lapse and shift. This is straight-forwardly done by applying finite difference operators to the variables. Note that in case of non-Cartesian local coordinates, the derivatives are calculated via the global derivatives as described in [9]. For practical reasons (particularly storage and memory requirements), it is sufficient to calculate the radial instead of the Cartesian derivatives. The Cartesian derivatives can later be reconstructed in terms of angular derivatives of the spherical harmonics which are known analytically. This saves two additional grid functions per ADM variable. We obtain the radial derivatives in terms of the Cartesian derivatives via

∂ r = 1 r (x ∂ x +y ∂ y +z ∂ z ) . ( 24
)
This completes the first step of calculating all variables that are necessary for the construction of characteristic boundary data. What remains is the projection of these variables on to a coordinate sphere with fixed radius R and a subsequent decomposition in terms of scalar spherical harmonics.

To project the variables on to the sphere, we use fourth-order Lagrange interpolation. Afterwards, the variables are decomposed as

f m = S 2 dΩ Ȳ m f (Ω) , (25)
for all variables f whether scalar, vector or tensor quantities. The projection and decomposition are both done on polar spherical coordinate spheres where the surface integration is done with a Gauss quadrature scheme as described in [9], and the resulting array of spherical harmonic modes is stored in a file for further processing. The reconstruction from the data in that file will be described in the next subsection.

B. Reconstruction from harmonic modes

Once the boundary data in terms of the ADM variables and their time and radial derivatives are known on a coordinate sphere with radius R and stored in a file, the characteristic boundary module reads and interprets the data in that file in order to reconstruct the variables on S2 . That is, the characteristic code is run as a post-processing tool to obtain the gravitational radiation at J + .

The variables are reconstructed via

f = ,m f m Y m . ( 26
)
for any variable f .

The advantage of reconstructing the variables from harmonic modes is the independence in terms of angular coordinates, 1 as well as resolution. Additionally, the necessary finite -mode cutoff can act as as a filter that factors out angular high-frequency noise. In practice for the binary black hole models considered here, the amplitudes of higher modes falls off rapidly so that they can be considered essentially irrelevant numerically for, say, ≥ 10.

Note that for practical reasons, we have preferred to calculate the radial instead of the Cartesian derivatives of the ADM variables. As a consequence, we have the harmonic modes of the radial derivatives. In order to obtain the Cartesian derivatives instead, one can take angular derivatives of the spherical harmonics and then apply a Jacobian transformation from stereographic to Cartesian coordinates

∂ ∂ x i f = ∂ r ∂ x i ∂ r f m Y m (q, p) + f m ∂ q ∂ x i ∂ Y m ∂ q + f m ∂ p ∂ x i ∂ Y m ∂ p , ( 27
)
where q, p denote stereographic angular coordinates. By using the relation

A ,q = ðA + ðA -2ipsA 1 + q 2 + p 2 , A ,p = i( ðA -ðA + 2qsA) 1 + q 2 + p 2 , ( 28
)
for a quantity A with spin-weight s, we can express partial derivatives with respect to angular coordinates in terms of the eth-derivative

∂ ∂ x i f = ∂ r ∂ x i ∂ r f m Y m + f m ∂ q ∂ x i ðY m + ðY m 1 + q 2 + p 2 + f m ∂ p ∂ x i i ðY m -ðY m 1 + q 2 + p 2 , ( 29
)
where we have used the fact that the scalar spherical harmonics have spin-weight s = 0. This can be further simplified by replacing the eth-derivatives of the spherical harmonics by spin-weight s = ±1 spherical harmonics. Finally, we can write

∂ ∂ x i f = f m ( + 1) 1 + q 2 + p 2 [ 1 Y m (q ,i -ip ,i ) --1 Y m (q ,i + ip ,i )] + r ,i f m ,r Y m , ( 30
)
where (32) p ,i = ∂ i ±y 2 ∓xy, ±r + z ∓ y 2 /r, -y ∓ yz/r , (33) where the upper sign is valid for the north patch and the lower sign is valid for the south patch.

r ,i = ∂ i x 2 + y 2 + z 2 = x i r , (31)
q ,i = ∂ i x x 2 + y 2 + z 2 ± z = 1 (r ± z) 2 r ± z -x 2 /r, -xy/r, -xz/r ∓ x ,
x 2 + y 2 + z 2 ± z = 1 (r ± z)
This means that we do not have to calculate any extra derivatives numerically and by virtue of Eq. ( 30), we immediately obtain the Cartesian derivatives from the radial derivatives.

Once the ADM variables and their derivatives are reconstructed on the coordinate sphere, the code executes the computations discussed in Sec. II D to obtain the boundary data in terms of the characteristic evolution variables and in terms of the Bondi coordinates. Once this step is complete, the code can use the boundary data to evolve the characteristic evolution variables further in time.

Results for the application of this method to a binary black hole merger are presented in Sec. V. However, the fields at J + first have to be transformed to constant Bondi time u B . This is described in the next subsection.

After the evolution of the characteristic variables to J + , the news function N and the Weyl scalar Ψ 4 are computed in terms of inertial spatial Bondi coordinates (r B , p B , q B ), but are not yet given in terms of the inertial Bondi time coordinate u B , i.e., we have N = N (r B , p B , q B , u). Also, one deals with the complete 2D data at J + , but for analysis purposes, it is more convenient to have the data decomposed in terms of spherical harmonic modes.

In a previous implementation, the transformation to constant u B was achieved by storing the news function as a function on S 2 at each time-step to a file. In a post-processing step, these timeslices were then read and interpolated to constant Bondi time. This, however, is rather inconvenient to use in practice. Rather, one would like to have the variables automatically transformed to the correct time coordinate and decomposed in terms of spherical harmonic modes, avoiding the need for additional post-processing tools. For this reason, we have extended the module for calculating the news function and the Weyl scalar Ψ 4 to take care of these issues.

During characteristic evolution, we know the inertial Bondi time in terms of coordinate time and angular coordinates u B = u B (u, p, q) at each time-step at J + . We have found that in practice the difference between u and u B is very small, so that by keeping five time levels of u B , N and Ψ 4 in memory at each time-step, it is possible to interpolate N and Ψ 4 to u B = const. for each point on S 2 . This is done by means of fourth-order Lagrange interpolation on the fly during evolution.

In practice, we average the Bondi time over S 2 on the past-past time-level to find the target time u B = const. We use the past-past time-level to have the interpolation stencil centered around the target interpolation time to maintain maximal accuracy in the interpolation. Once the target time is known, N and Ψ 4 are interpolated for each point on S 2 to u B = const. so that we finally have

N = N (r B , p B , q B , u B ) and Ψ 4 = Ψ 4 (r B , p B , q B , u B ).
Afterwards and if requested by the user, the code can automatically decompose the two functions in terms of spin-weighted spherical harmonic modes. Therefore, the data is available in a format that is convenient for analysis without any further post-processing. A set of post-processing tools can read in the harmonic modes of Ψ 4 and calculate radiated energy, linear and angular momentum as well as the wave-strain h.

However, there is one word of caution required. As mentioned earlier, the characteristic code computes a quantity Ψ 4 which is related to the commonly used ψ 4 via Eq. (18). It follows that the harmonic modes are related by

ψ m 4 = -2 Ψ ,-m 4 (-1) m . ( 34
)

D. The linearized conformal factor

In order to compute the news N and Ψ 4 , it is first necessary to calculate the conformal factor ω of the conformal compactification of the Bondi metric at J + . However, the full non-linear computation of this factor can display a high level of noise when it is implemented numerically, reducing its accuracy. Since in all cases that we consider the fields at J + are in a linear regime, we can make use of a relation between the conformal factor ω and J at J + in the linearized approximation.

Given J at J + in terms of s = 2 spin-weighted spherical harmonics, i.e.,

J| J + (u, x A ) = ≥2,m J ,m (u) 2 Y ,m , (35)
we can write ω in terms of scalar spherical harmonics as

ω(u, x A ) = 1 + ≥2,m ω ,m (u)Y ,m , (36)
where

ω ,m (u) = - 1 4 ( + 1) ( -1)( + 2) J ,m (u) + (-1) m J ,-m (u) . ( 37
)
The linearized computation has been implemented as an alternative to the full non-linear computation of ω and either method can be selected at run-time.

IV. TESTS OF THE CHARACTERISTIC EXTRACTION CODE

In this section, we report on tests we have performed in order to calibrate the characteristic extraction code. For this, we put an analytical solution to the ADM metric variables on the Cauchy side, generate the associated spherical harmonic modes on each timestep and afterward run the characteristic extraction code.

A. Linearized solutions

We have tested convergence of the entire extraction process against an analytical (linearized) solution [47] by transforming this solution to the ADM variables, i.e., lapse α, shift β i and the 3-metric components γ ij , on the Cauchy side. According to the description in Sec. III, we decompose the ADM variables on a world-tube after each timestep in order to produce spherical harmonic modes that can be used as boundary data for the subsequent characteristic evolution. On the characteristic side, we measure the error, that is, the difference between evolved variables and exact fields. We do so by performing a set of two resolutions, and we can assume convergence against the linearized solution given that the amplitudes of the fields are small, i.e., they are in the linearized regime.

The solutions to be used are those given in Sec. 4.3 of ref. [47] for the case of a dynamic spacetime on a Minkowski background with = 2, m = 0. We write

J = ( -1) ( + 1)( + 2) 2 Y m (J (r)e iνu ), U = ( + 1) 1 Y m (U (r)e iνu ), β = Y m (β e iνu ),
W c = Y m (W c (r)e iνu ), (38)
where J (r), U (r), β , W c (r) are in general complex, and taking the real part leads to cos(νu)

and sin(νu) terms. The quantities β and W c are real; while J and U are complex due to the terms 2 Y m and 1 Y m , representing different terms in the angular part of the metric. In the case = 2

β 2 = β 0 , J 2 (r) = 24β 0 + 3iνC 1 -iν 3 C 2 36 + C 1 4r - C 2 12r 3 , U 2 (r) = -24iνβ 0 + 3ν 2 C 1 -ν 4 C 2 36 + 2β 0 r + C 1 2r 2 + iνC 2 3r 3 + C 2 4r 4 , W c2 (r) = 24iνβ 0 -3ν 2 C 1 + ν 4 C 2 6 + 3iνC 1 -6β 0 -iν 3 C 2 3r - ν 2 C 2 r 2 + iνC 2 r 3 + C 2 2r 4 , ( 39
)
with the (complex) constants β 0 , C 1 and C 2 freely specifiable.

We present results for the case = 2, m = 0, ν = 1, with

C 1 = 10 -6 , C 2 = 3 × 10 -6 , β 0 = 0 (40)
in Eq. (39).

The numerical simulations use a compactified radial coordinate, Eq. ( 4), with r wt = 9. Data is prescribed at time u = 0 and extracted at a number of different world-tube locations, specifically at r = {30, 60, 90, 120}. The stereographic grids (with ghost zones excluded) are

Coarse: n x = n q = n p = 41, Fine:

n x = n q = n p = 81; (41)
J| h=0.20M J + × 2 J| h=0.40M J +
FIG. 2: First-order convergence in the L 2 -norm of J at J + against a linearized solution as given by (39). and there is no overlap between the two patches, i.e., we set the code parameter q size = 1 which means that on the nominal grid, the holomorphic coordinate function ζ = q + ip takes values in q, p ∈ [-1, 1]. The compactified coordinate ranges from, x ∈ [x in , 1] such that the world-tube is located on the characteristic grid and close to the inner boundary to avoid accuracy losses. For example, for a world-tube located at r = 30M, we set x in = 0.74. For different world-tube locations, x in may be changed accordingly.

The grids on the Cauchy side have a radial resolution of h = 0.4M and h = 0.2M with the number of angular points per direction and per patch set to n = 21 and n = 41, respectively. The timestep size is taken to be ∆t = 0.1M and ∆t = 0.05M.

In Fig. 2, we show convergence in the L 2 -norm of J at J + to the linearized solution, starting off from a world-tube at r = 30M. We observe clear first-order convergence. This is consistent with the fact that the boundary data for U ,r is only known to first-order off the world-tube. However, for other models presented in subsequent sections, it appears that the coefficient of numerical error arising from U ,r is small enough that a second-order convergence exponent can be measured, corresponding to the interior finite difference accuracy of the characteristic code.

B. Gauge invariance of a shifted Schwarzschild black hole

An important aspect of CCE as a wave extraction method is the fact that it should produce invariant results independent of the gauge at the world tube. As a simple test problem which demonstrates this property, we place a single static Schwarzschild black hole at the origin of our spacetime, and apply a shift to the solution, according to

β x = Aω cos(ωt) , β y = β z = 0 , ( 42
)
so that the x-coordinate oscillates as

x = A sin(ωt) , ( 43
)
where the amplitude is set to A = 1.0 and the frequency is set to ω = 0.1. As the spacetime is spherically symmetric and static, the resulting radiation content must be zero, and hence, in numerical simulations the residual in Ψ 4 and the news N must converge to zero. This is indeed what we find. Fig. 3 shows the ( , m) = (2, 0) mode of Ψ 4 and N at two resolutions h = 0.4M and h = 0.2M scaled for second-order convergence. All other modes converge to zero at the same order.

In contrast to this, finite radius extraction methods based on ψ 4 and on gauge-invariant perturbations of Schwarzschild fail to be convergent. This is shown in Fig. 4, where we plot the perturbative Zerilli-Moncrief (even) master function Q + =2,m=0 [64,72,73] and ψ 4 for the two resolutions h = 0.4M and h = 0.2M without any rescaling. The two waves are on top of each other indicating that the results do not converge to zero.

Although the effect is exaggerated by this artificial test, the point remains that the commonly used finite radius extraction methods are susceptible to gauge variations which are difficult to disentangle from the genuine wave signal. The CCE method as implemented removes this systematic error completely.

V. RESULTS FOR BINARY BLACK HOLE INSPIRALS

We have carried out full three-dimensional, non-linear evolutions of two binary black hole systems for the whole spacetime out to J + using the Cauchy and characteristic formulations as described in the previous sections. The first configuration represents an equal-mass non-spinning binary, with modes up to = 8 of the CCE waveform illustrated in Fig. 5; while the second configuration is an equal-mass binary with individual spins that are aligned with the orbital angular momentum and hence the system is not subject to precession effects. In this section, we compare the differences in the waveform measures at finite radius and at J + , and compute various features of the two binaries. and ψ 4 computed at a finite radius (right) extracted at r = 100M . For each case, we plot two resolutions h = 0.20M and h = 0.40M of the ( , m) = (2, 0) mode without any rescaling. We find that the waveforms do not converge to zero, as should be expected for this test case.

Q + ×10 -4 Q +,h=0.40M Q +,h=0.20M
We demonstrate that for these tests, the wave-extraction at J + contains only numerical error.

Additionally, we measure the energy, linear and angular momenta radiated through gravitational waves, both at finite radii and at J + and test for conservation of energy and angular momentum.

As the emission of gravitational radiation is symmetric in the first case, but anti-symmetric in the second case, the net-linear momentum carried by the gravitational waves is non-zero for the second binary. Hence, this binary gives rise to a gravitational recoil, or kick, of the merger remnant. We measure this effect to high precision using gravitational waveform data at finite radii and at J + .

A. Initial data

For the first equal-mass non-spinning binary, the numerical evolution starts from an initial separation d/M = 11.0, and goes through approximately 8 orbits (a physical time of around 1360M), merger and ring-down. The masses of the punctures are set to m = 0.4872 and are initially placed on the x-axis with momenta p = (∓0.000728, ±0.0903, 0), giving the initial slice an ADM mass M ADM = 0.99051968 ± 2 × 10 -8 . The z-axis is chosen to coincide with the initial orbital angular momentum and is given in Table I. The y-axis is then chosen perpendicular to the other axis.

The second binary has been investigated mainly so as to be able to measure the kick-velocity that the merger remnant acquires due to the asymmetric emission of gravitational radiation. As this effect depends largely on the very last orbit (compare e.g. [10]), it is sufficient to consider a binary with a closer initial separation which in our case has been chosen to be d/M = 4.0.

The initial masses of the punctures are set to m = 0.2999 and are initially placed on the xaxis with momenta p = (∓0.00211428, ±0.112539, 0), giving the initial slice an ADM mass

M ADM = 0.98826246 ± 2 × 10 -8 .
The coordinate system is chosen such that the initial orbital angular momentum is aligned with the z-axis of the coordinate system, and we represent the initial spin vectors by

S 1 = (0, 0, a 1 ) , S 2 = (0, 0, a 2 ) , ( 44
)
where we have chosen (a 1 , a 2 ) = (0, 0) for the first binary and (a 1 , a 2 ) = (0.8, -0.8) for the second one.

The initial data parameters were determined using a post-Newtonian evolution from large initial separation, following the procedure outlined in [14], with the conservative part of the Hamiltonian accurate to 3PN, and radiation-reaction to 3.5PN. The measured eccentricity of the resulting orbits is 0.004 ± 0.0005

The most convenient way to set the initial data on the initial null-hypersurface of the characteristic grid is to set J(u = 0) = 0 everywhere, since the wave-extraction algorithm requires J(u = 0) |J + = 0. Fortunately, it turns out that J = 0 on the extraction world-tube Γ at u = 0 is consistent with conformal flatness, although in general there will be an inconsistency at points on the null cone off the initial Cauchy hypersurface.

B. Grid setup

The binary black hole evolution was carried out on a 7-patch grid structure, as described in [9,60], incorporating a Cartesian mesh-refined region which covers the near-zone, and six radially oriented patches covering the wave-zone. The inner boundary of the radial grids was placed at r t = 35.2M relative to the centre of the Cartesian grid. The outer boundary for the spherical grids was chosen based on the expected time duration of the measurement and radius of the furthest extraction world-tube, in order to remove any influence of the artificial outer boundary condition.

In particular, given that the evolution takes a time T m for the entire inspiral, merger and ringdown, and a world-tube location at a finite radius r Γ , we would like to ensure that a disturbance traveling at the speed of light c = 1 [74,75] from the outer boundary does not reach the world-tube radius (see Fig. 6) 2 . Thus we place the boundary at

r b > T m + 2r Γ . ( 45
)
For the particular evolutions considered here, for the non-spinning binary and for the spinning The near-zone grids incorporate 5 levels of 2:1 mesh refinement, covering regions centred around each of the black holes. For the highest resolution we have considered here, the finest grid (covering the BH horizon) has a grid spacing of dx = 0.02M. The wave-zone grids have an inner radial resolution which is commensurate with the coarse Cartesian grid resolution, dr = 0.64M in this case. This resolution is maintained essentially constant to the outermost measurement radius (r = 1000M), at which point we apply a gradual decrease in resolution (see [9]) over a distance of r = 500M. From r = 1500M to the outer boundary, we maintain a resolution of dx = 2.56M, sufficient to resolve the inspiral frequencies of the dominant ( , m) = (2, 2) mode of the gravitational wave signal. The angular coordinates have 31 points (30 cells) in ρ and σ on each of the 6 patches. The time-step of the wave-zone grids is dt = 0.144, and we take wave measurements at each iteration.

We have carried out evolutions at three resolutions, h = 0.64M, h = 0.80M and h = 0.96M, in order to estimate the convergence of our numerical methods (compare [9]).

The characteristic extraction world-tubes Γ have been placed at r Γ = 100M and r Γ = 200M

for the first binary, while for the second binary we have chosen r Γ = 100M and r Γ = 250M. We note that over the course of the evolution, these world-tubes are causally disconnected from the artificial outer boundary.

C. Invariance with respect to world-tube location

Invariance with respect to the world tube location is demonstrated in Fig. 7. We have considered the differences between waveforms at J + resulting from two independent characteristic evolutions using boundary data at r Γ = 100M and r Γ = 200M, respectively, and for two resolutions, h = 0.96M and h = 0.64M. The difference between the results should be entirely due to the discretisation error, and indeed this is what we find. The differences converge to zero with approximately second-order accuracy, as expected for the null evolution code. The figure displays the differences in the amplitude and phase of the wave mode ψ 4 ( = 2, m = 2) for resolution h = 0.96M and h = 0.64M scaled for second-order convergence. The same order of convergence is also obtained for higher order modes such as ψ 4 ( = 4, m = 4) and ψ 4 ( = 6, m = 6) during inspiral. However, the higher modes are more sensitive to the order of convergence on the Cauchy side during merger, and hence, the error goes down with a higher convergence order (note that we use 8th-order finite difference operators). The phase of ( , m) = (6, 6) mode converges at sixth-order while the amplitude converges at third order.

The differences between the waveforms at J + for resolution h = 0.64M are of order of 0.03% in amplitude with a dephasing of 0.002 radians over the course of the evolution.

It should be noted that the differences in the waves between the two world-tube locations is one order of magnitude smaller than the numerical errors in amplitude and phase inherent in the Cauchy evolution [9]. We therefore expect that the dominant error in the waveforms is due to the discretisation error during Cauchy evolution. Indeed, using the three resolutions h = 0.96M, h = 0.80M and h = 0.64M, we find a very similar convergence order for (a) CCE with fixed worldtube location, and (b) for a pure Cauchy evolution with the same fixed finite radius extraction (compare [9]).

However, it should also be noted that the (non-convergent) differences between extrapolated waveforms and waveforms obtained at J + as found in [4] are of the same order as the (convergent) discretisation error inherent in the Cauchy evolution for our given resolutions as found in [60].

Further increasing the resolution of the Cauchy evolution will therefore lead to an error that is much smaller than the error made by finite radius extrapolation. This highlights the importance of the systematic error due to extrapolation from finite radius, which is removed by performing computations at J + via CCE.

D. Conservation of energy and angular momentum

We have used the measured waveforms to compute the initial mass and angular momentum of the modeled spacetimes. The initial mass is given by the ADM mass as computed at infinity on the initial slice using the methods described in [START_REF] Ansorg | [END_REF]. Since the initial slice is conformally flat, the initial angular momentum of the spacetime is simply given by (see for example [76][77][78])

J ADM = C 1 × p 1 + C 2 × p 2 + S 1 + S 2 . ( 46
)
Here C i , p i and S i are the position, the linear momentum and the spin of the i-th black hole.

After evolution, i.e., after the horizon of the final black hole has equilibriated, we compute the mass and spin of the remnant on its apparent horizon using the isolated horizon formalism [79][80][81],

or alternatively, by using the horizon circumference formula [82] (the two methods agree within where the latter is scaled for second-order convergence. The top right panels show the the same for the ( , m) = (4, 4) mode of ψ 4 . The lower panels show convergence of the ( , m) = (6, 6) mode. The phase error for the high resolution is rescaled for sixth-order convergence while the amplitude error is rescaled for third-order convergence.

the reported errors as given in Table I, where we use the values obtained from the circumference formula 3 ). In addition, we compute the radiated energy and angular momentum carried by the gravitational waves using the results of finite radius computations, from extraction at J + , and from an extrapolation of the results at several finite radii. The finite radius computations are obtained at an extraction world-tube r Γ = 100M, which is typical for current numerical relativity simulations (see e.g., [83,84]).

The radiated energy and angular momentum flux can be calculated in terms of spin-weighted spherical harmonic coefficients of ψ 4 as derived in [85]. The radiated energy flux is given by

dE dt = lim r→∞ r 2 16π ,m t -∞ d tψ m 4 2 (47)
and the components of the radiated angular momentum flux read

dJ x dt = lim r→∞ r 2 32π ,m t -∞ dt t -∞ dt ψ ,m 4 × t -∞ dt f ,m ψ ,m+1 4 + f l,-m ψ ,m-1 4 , ( 48
)
dJ y dt = lim r→∞ r 2 32π ,m t -∞ dt t -∞ dt ψ ,m 4 × t -∞ dt f ,m ψ ,m+1 4 -f l,-m ψ ,m-1 4 , (49)
dJ z dt = lim r→∞ r 2 16π ,m t -∞ dt t -∞ dt ψ ,m 4 × t -∞ dt ψ ,m 4 , (50)
where

f ,m = ( + 1) -m(m + 1). (51)
Note that we can also compute the expressions for radiated energy and angular momentum in terms of the news function N by using the fact that in a Bondi-frame, Ψ 4 is related to the news via

Ψ 4 = ∂ u N . ( 52
)
Hence, we can replace all first time integrals of ψ 4 in ( 47)-( 50) by N .

We find that modes are well resolved only up to ( , m) = (8, 8), and hence the sum in the expression above is to = 8. We have computed the sums for waveforms extracted at finite radii and at J + . The finite radius computations have been performed on extraction Model (a 1 , a 2 ) = (0, 0) (a 1 , a 2 ) = (0.8, -0.8)

total ADM mass M ADM 0.99051968 ± 2 × 10 -8 0.98826246 ± 2 × 10 -8 initial angular momentum J ADM 0.993300 ± 1.0 × 10 -16 0.900312 ± 1.0 × 10 -16

irr. Mass M f,irr 0.884355 ± 2.0 × 10 -5 0.883636

± 1.0 × 10 -4 Spin S f /M 2 f 0.686923 ± 1.0 × 10 -5 0.686300 ± 2.7 × 10 -4
Mass M f 0.951764 ± 2.0 × 10 -5 0.950829 ± 3.8 × 10 -5

Angular momentum S f 0.622252 ± 1.0 × 10 -5 0.620469 ± 1.9 × 10 at the same fixed resolution. The error in the initial angular momentum J ADM is at machine precision as it is known analytically in terms of the initial parameters of the black holes. All other errors are given by the Richardson expanded differences between a medium and a high resolution run as described in [9]. The last three rows display the differences between extrapolated results and results at J + , which we note is of the order of the discretisation error as indicated by the finite radius error-bars.

spheres r = {280M, 300M, 400M, 500M, 600M, 1000M} for the non-spinning binary, and r = {100M, 250M, 500M} for the spinning binary. Given the radiated energy and angular momentum on the three extraction spheres, we do a linear extrapolation by fitting a function of the form

Q(r) = Q 0 + Q 1 /r , ( 53
)
where Q is either E rad or J rad .

The values for the initial mass M ADM , initial angular momentum J ADM , final (Christodoulou) mass M f , final angular momentum S f , as well as radiated energy E rad and radiated angular momentum J rad are reported in Table I for both binary systems. The reported errors refer to Richardson expanded differences between a high and a medium resolution run as described in [9], except for the error in the extrapolation and the error in CCE, where for the former we report the error in the fit, and for the latter report the differences between two given extraction world-tubes for the same fixed resolution. Hence, they do not contain the discretisation error of the evolution, which, in all cases, is one order of magnitude larger than the error in the extrapolation and in CCE.

Table I lists the radiated quantities determined from the measured values of ψ 4 . However, since we also determine the news function N at J + , we can independently compute E rad,CCE and J rad,CCE based on N in order to cross check the result. We find that the differences in all cases are on the order of the reported error bars, hence further validating the results based on CCE.

Another important point to make is that the differences between extrapolated data and data obtained at J + are of the same order of magnitude as the discretisation error in the Cauchy evolution. This is indicated by the last three rows of Table I. The discretisation error of these quantities due to the the Cauchy evolution is given by the error-bars of the finite radius quantities.

Due to conservation of energy and angular momentum, the following relations must hold exactly

S f + J rad = J ADM , M f + E rad = M ADM . ( 54
)
The residuals of these relations are reported in Table II for the quantities at a finite radius, extrapolated to infinity, and computed at J + , for both binary systems.

Note that the residuals in the extrapolated values of the two binaries are of comparable accuracy to the values obtained via CCE. This is surprising, considering the error in the extrapolation procedure itself, and demonstrates that extrapolation, even in the most extreme spinning cases, delivers comparable performance. Note, however, that our extrapolation uses large extraction radii

|S f + J rad,finite -J ADM | 3.1 × 10 -2 1.3 × 10 -2 |S f + J rad,extrap -J ADM | 1.0 × 10 -3 3.9 × 10 -3 |S f + J rad,CCE -J ADM | 3.4 × 10 -4 7.1 × 10 -3 |M f + E rad,finite -M ADM | 1.4 × 10 -3 2.1 × 10 -3 |M f + E rad,extrap -M ADM | 4.1 × 10 -5 2.4 × 10 -4 |M f + E rad,CCE -M ADM | 7.2 × 10 -5 2.6 × 10 -4
TABLE II: Residuals in conservation of energy and angular momentum (54) for the non-spinning and spinning configuration and for radiated mass E rad and angular momentum J rad as calculate at a finite radius r = 100M , extrapolated from finite radii computations to infinity, and by the CCE computation at J + . Both, extrapolated quantities and those directly computed at J + satisfy the conservation conditions with about the same amount of accuracy. Not using any extrapolation but relying on a pure finite radius computation leads to an error that is of about 1-2 orders of magnitude larger than the one found by the other two methods.

compared to common practice in numerical relativity. By using a set of smaller extraction radii in the extrapolation, we expect that the errors increase up to one order of magnitude (compare [60]). In addition, the systematic error made in the extrapolation will not converge with numerical resolution, but only linearly with extraction radius O(r). This however, puts high computational demands on future high accuracy simulations.

In all cases, we observe that a pure finite radius computation based on an extraction sphere at r = 100M clearly performs worse and has an error that is larger by 1-2 orders of magnitude, which, in the case of angular momentum, corresponds to an error of the order of 1.0% -3.0%.

E. Recoil velocity

In addition to energy and angular momentum, gravitational waves also carry linear momentum.

Depending on the binary configuration, there might be a prefered direction at which radiation, and hence linear momentum, is beamed. Due to conservation of linear momentum, this gives rise to a gravitational recoil, or "kick" acquired by the merger remnant. An accurate measurement of this is effect is crucial in various fields of astrophysics as this determines whether the final black hole is kicked out of its host environment. Clearly, the existence or non-existence of a black hole can have a dramatic impact on the host and determines its further evolution.

Previous studies have already considered this effect in detail, and particularly, the binary configuration that we consider here belongs to a sequence that has already been examined in [86,87].

However, these studies estimated the recoil velocity using waveforms extracted at a finite radius.

In this subsection, we measure this effect where it is defined unambiguously, namely at J + .

The radiated linear momentum flux can be calculated in terms of spin-weighted spherical harmonic coefficients as derived in [85] and given by

dP x dt + i dP y dt = lim r→∞ r 2 8π ,m t -∞ dt ψ ,m 4 × t -∞ dt a ,m ψ ,m+1 4 + b l,-m ψ -1,m+1 4 -b l+1,m+1 ψ +1,m+1 4 , ( 55
)
dP z dt = lim r→∞ r 2 16π ,m t -∞ dt ψ ,m 4 × t -∞ dt c ,m ψ ,m 4 + d ,m ψ -1,m 4 + d +1,m ψ +1,m 4 , (56)
where

a ,m = ( -m)( + m + 1) ( + 1) , ( 57
) b ,m = 1 2 ( -2)( + 2)( + m)( + m -1) (2 -1)(2 + 1) , ( 58
) c ,m = 2m ( + 1) , ( 59
)
d ,m = 1 ( -2)( + 2)( -m)( + m) (2 -1)(2 + 1) . ( 60
)
Again, we have computed the expression above by including modes up to ≤ 8. In addition, we calculate this effect using a finite radius computation and a computation at J + . Given the three extraction spheres at r = {100M, 250M, 500M} for the spinning configuration, we have extrapolated the recoil velocity to infinity by means of Eq. ( 53).

It should be noted that it is crucial to include the appropriate vector integration constant when integrating the linear momentum flux in time. As analyzed in [10], there are two important contributions. The first mainly comes from the initial burst of junk radiation that adds significantly to the measured recoil at the beginning of the simulation. The second is due to the non-zero net-linear momentum that the binary had acquired if it had inspiralled from an infinite separation. We have measured the correct vector integration constant by applying the method described in [10] and corrected the computation accordingly. We find good agreement between finite radius computations, extrapolation to infinity, and computation at J + . The results are included in Table I.

F. Constraint preservation

We check the preservation of the constraints at the world-tube and at J + for different radii and resolutions. We use the constraint equations as implemented and first tested with linearized solutions in [88]. As a consequence of the Bondi system, there are three constraint equations

R 00 = 0 , R 01 = 0 , q A R 0A = 0 . ( 61
)
Analytically, if these constraints are satisfied initially and at the world-tube, they should be satisfied exactly everywhere and at all times. However, in computer simulations, the solution always contains numerical error so that the constraints can only be satisfied approximately. Since this error is only due to discretisation, it should converge away in the limit of infinite resolution. Indeed,

this is what we find.

Fig. 8 shows the convergence of R 00 at J + and q A R 0A and R 01 at the world-tube Γ for the nonspinning binary. While R 00 appears to be only first-order convergent at J + , we note that the error is already very small and at the order of 10 -10 and hence close to machine precision. The same applies to the constraint R 01 which exhibits a convergence exponent between first and secondorder. The constraint q A R 0A , on the other hand, exhibits an overall second-order convergence exponent at the world-tube. The origin of the spurious spikes occuring at the high-resolution run between -500 < t < -100 and reducing the convergence order, is currently unkown but may be due to the first-order knowledge of U ,r at the world-tube. This is subject to future investigation and improvements of the algorithm and code. We also note that during merger, the convergence of the q A R 0A constraint is reduced at the world-tube. This, again, may be due to the first-order knowledge of U ,r at the world-tube. However, the convergence of the constraints appeared to be problematic already in [88] and we emphasize that the differences in the waveforms at J + as extracted from different world-tube locations converge to zero at approximately second-order and independent of the constraints.

-1400 -1200 -1000 -800 -600 -400 -200

0 t/M 10 -10 10 -9 R 00 | J + R 00 | h=0.64M J + × 1.5 R 00 | h=0.96M J +
-1400 -1200 -1000 -800 -600 -400 -200 In all cases, the constraints exhibit a large spike at the time when the first non-trivial data, particularly the initial burst of junk radiation, reaches the world-tube location. Although the junk radiation should satisfy the constraint equations since it arises as part of the solution to the constraint preserving initial data at the initial Cauchy hypersurface, we note that there is principally an incompatibility between the conformal flatness assumption of the 3-metric on the initial Cauchy hypersurface and the initial data of the 4-metric, J(u = 0) = 0, on the initial null hypersurface.

0 t/M 10 -7 10 -6 10 -5 10 -4 10 -3 q A R 0A | Γ q A R 0A | h=0.64M Γ × 2.25 q A R 0A | h=0.96M Γ - 1400
This may lead to the observed initial spikes in the constraints. Further investigation is necessary to clearly identify the source of error and to improve the construction of characteristic initial data.

VI. NOTES ON CHARACTERISTIC CODE PARAMETERS

There are some subtleties regarding the choice of some code and grid parameters that have been found empirically during tests of the extraction method. Here, we collect some of these findings for future reference.

• We use conformally flat initial data on the initial null hypersurface, i.e.J(u = 0) = 0.

This assumption is in principle incompatible with the proper radiation content of a binary black hole spacetime, and it is subject to the same error as on the Cauchy side, where we use conformally flat initial data. In practice, the spurious burst of junk radiation at the beginning of the simulation quickly leaves the system and does not affect the radiation at the accuracy given by the numerics.

• The radial resolution of the physical coordinate on the characteristic grid should be at least as fine at the world-tube as the radial resolution on the Cauchy side. We can estimate the compactified radial resolution in terms of the physical resolution via

dr dx = r wt 1 1 -x + x (1 -x) 2 , ( 62
)
which asymptotes to r wt for x → 0 at the world-tube, i.e.

x → 0, dr dx → r wt . (63) It follows that

∆r| Γ ∼ r wt ∆x (64)
at the world-tube. We can then select an appropriate ∆x, so that ∆r| Γ at the world-tube is at least as small as the radial resolution ∆r on the Cauchy side.

It turns out that at least N x ≈ 200 radial points are required to resolve the wave properly at J + . Since ∆R = 0.64M at the extraction spheres for the highest resolution, the compactification parameter has to be r wt ≈ 256 in order to meet the condition ∆r ≤ ∆R. Further increasing the characteristic radial resolution while leaving the Cauchy ∆R = const. has no noticeable influence on the wave at J + .

• The linearized approximation for the conformal factor ω is used for the computation of the news N . The full non-linear computation results in spurious drifts and does not converge when considering the Schwarzschild test. This is potentially a bug in the existing implementation and affects only the news. Ψ 4 is only affected at higher-order which does not play a role in practice. However, as the fields at J + are in the linear regime, the linear approximation is a valid simplification.

• We use the first-order evolution system in the characteristic code, i.e., we use only first ethderivatives. The implementation with second eth-derivatives does not converge and gives inconsistent results for the binary waveform when considering different world-tube locations. This might also be related to the problem formulated in [39,89]. In principle, both methods should be equivalent and the non-convergence of the second-order system can be attributed to a bug in the code.

• The news and Ψ 4 are computed using the second-order implementation of the extraction module, i.e., using second eth-derivatives.

• All quantities are decomposed in terms of the complex-valued spin-weighted spherical harmonics s Y m .

VII. CONCLUSIONS

We have demonstrated results from an implementation of characteristic extraction in numerical relativity, applied to the astrophysically relevant test cases of binary black holes. Our implementation carries out the characteristic extraction as a post-processing step. As such, the code is general purpose and can be applied to boundary data generated by any Cauchy simulation, independent of numerical scheme or formulation of the Einstein equations. The code has been tested against analytic solutions, and it was found that it exhibits the expected convergence properties, as well as gauge-independence, a property lacking in finite radius calculations. Applied to two test cases involving the inspiral and merger of binary black holes, we have determined the emitted gravitational radiation, and found that the CCE method does indeed generate results which are independent of the world-tube on which the data has been extracted. We have also determined the residuals in the conservation of mass and angular momentum, as well as the recoil velocity (when the black hole spins are anti-aligned). We find that the residuals of the CCE method are very similar to those obtained by extrapolation to infinity of results obtained using a pure Cauchy evolution, providing an important confirmation of the usefulness of such methods. However, the (non-convergent) differences between both measures is already on the order of the discretisation error of the Cauchy evolution, therefore indicating that as the resolution is further increased, the finite radius extrapolation will be the dominant source of error in the wave-extraction. The error in CCE, as measured here for a particular set of resolutions resolutions is one order of magnitude smaller, and will converge away if the resolution is further increased.

Characteristic extraction yields results that are free of a number of systematic errors in the modeling of the physical problem. Placing the outer boundary at a sufficiently large distance, we are able to causally disconnect the artificial boundary from the extraction world-tube. However, there is still the possibility of a systematic error in the extent to which there is an unphysical radiation content in the initial data. Further, it should be noted that there is, in general, an inconsistency between the conformally flat initial data used in the Cauchy computation, and that of J = 0 used on the initial characteristic null cone. It is usually assumed that the spurious initial radiation will quickly be flushed out of the system in a burst of initial junk radiation. The observed convergence of the differences between the radiation signals at different extraction radii, indicates that, after the burst of junk radiation, the effect of spurious, unphysical initial data is very small and below the discretisation error.

Future developments will clearly include applying characteristic extraction to a variety of other physical systems. It will also be interesting to revisit the problem of Cauchy-characteristic matching, in which the extraction world-tube and the Cauchy outer boundary coincide and boundary information is passed in both directions during the course of the evolution. In this case evolution in the Cauchy domain between, say r = 100M and r = 3000M, is avoided so the method is significantly more efficient; and further one can avoid inconsistencies in the initial data between the two domains. Finally, it would be useful to construct initial data that correctly models the physical problem; this is, in principle, possible at least for the characteristic domain.

Geodesic 2 αFIG. 1 :

 21 FIG. 1: Schematic illustration of the (first stage) construction of characteristic coordinates and metric.

FIG. 3 :

 3 FIG. 3: Second-order convergence of the gravitational wave signal to zero for the shifted-Schwarzschild test case. The left panel shows the = 2, m = 0 mode of the Weyl scalar Ψ 4 for the two resolutions h = 0.40M and h = 0.20M where the latter is scaled for second-order convergence. The right panel shows the ( , m) = (2, 0) mode of the news N for the same resolutions and the same rescaling. Similar plots for all other modes exhibit the same order of convergence.

4 FIG. 4 :

 44 FIG.4: Non-convergent behavior of finite-radius measures of the gravitational radiation for the shifted-Schwarzschild test case. Shown are the perturbative Zerilli-Moncrief (even) master function Q + (left)

4 ×10 - 1 l = 8 , m = 8 l = 6 , m = 6 l = 4 , m = 4 l = 2 , m = 2 FIG. 5 :

 41886644225 FIG.5: Dominant spherical harmonic modes of ψ 4 for ≤ 8 extracted at J + for the model a 1 = a 2 = 0.

FIG. 6 :

 6 FIG. 6: Portions of the Kruskal diagram which are determined by the numerical evolution. The green horizontal lines indicate the region spacetime that is determined by the Cauchy evolution and which has finite spatial extent with artificial outer boundary at r b . The red diagonal lines indicate the region that isdetermined by characteristic evolution, and which has an inner boundary at a world-tube Γ located at r Γ , consisting of data from the intersecting slices (green horizontal lines) of the Cauchy evolution. Note that due to the final spatial extent of the latter, the future Cauchy horizon of the Cauchy initial data, indicated by the dotted diagonal line parallel to J + , principally prohibits the determination of the spacetime out to timelike infinity i + . However, as long as the world-tube Γ is located within the future Cauchy horizon, the numerically evolved subset of the spacetime is determined consistently. Particularly, the artificial outer boundary r b is sufficiently removed such that the evolved binary black hole spacetime forms a merged black hole before the Cauchy horizon reaches the world-tube.

  binary respectively, we have evolution times T m 1350M and T m 550M, and outermost extraction world-tubes r Γ = 1000M and r Γ = 500M. We have placed the outer boundary of the evolution domains at r b = 3600M and r b = 2000M, respectively.

FIG. 7 :

 7 FIG. 7: Convergence of the differences in wave-signal to zero for various test cases. The top left panels show the amplitude and phase differences between world-tube data at r Γ = 100M and r Γ = 200M of the ( , m) = (2, 2) mode of the Weyl scalar ψ 4 for the two resolutions h = 0.96M and h = 0.64M

25 RFIG. 8 :

 258 FIG.8:The left panel shows the L 2 -norm of the constraint R 00 at J + for the two resolutions h = 0.96M and h = 0.64M scaled for first-order convergence. The right figure shows the L 2 -norm of the constraint q A R 0A at the world-tube Γ for the two resolutions h = 0.96M and h = 0.64M scaled for second-order convergence. The bottom panel shows the same for the constraint R 01 at the world-tube Γ. All figures refer to the non-spinning binary.

  Radiated energy E rad,finite,r=100 0.037354 ± 3.9 × 10 -5 0.035369 ± 3.7 × 10 -4

			-4
	Recoil velocity v kick,CCE	-	364.30 ± 0.32 km/s
	Radiated angular momentum J rad,CCE	0.370704 ± 6.2 × 10 -5	0.286908 ± 2.6 × 10 -4
	Radiated energy E rad,CCE	0.038828 ± 2.5 × 10 -6	0.037689 ± 4.3 × 10 -5
	Recoil velocity v kick,finite,r=100	-	350.37 ± 1.46 km/s
	Radiated angular momentum J rad,finite,r=100	0.339725 ± 4.0 × 10 -4	0.266654 ± 5.7 × 10 -3
	Residual |v kick,CCE -v kick,extrap |	-	1.43km/s
	Residual |J rad,CCE -J rad,extrap |	7.0 × 10 -4	3.1 × 10 -3
	Residual |E rad,CCE -E rad,extrap |	1.1 × 10 -4	5.0 × 10 -4

Recoil velocity v kick,extrap -362.87 ± 1.19 km/s Radiated angular momentum J rad,extrap 0.370007 ± 6.8 × 10 -5 0.283738 ± 2.1 × 10 -3 Radiated energy E rad,extrap 0.038715 ± 2.2 × 10 -6 0.037189 ± 2.2 × 10 -4

TABLE I :

 I Properties of the initial spacetime, of the merger remnant, and of the emitted radiation calculated from the appropriate waveform. The reported errors of the extrapolated values refer to the error in the extrapolation and the errors of the values at J + refer to the error from two different extraction world-tubes

On the Cauchy side, the ADM quantities are represented in polar-spherical coordinates on S

. The characteristic code, on the other hand, works with a stereographic coordinate mapping of S 2 .

The 1 + log slicing condition which we use propagates at √ 2c[55], however this is a gauge mode and empirically we find it to have negligible effect on finite radius measurements. Particularly, since CCE is gauge invariant, this does not affect the waveforms at J + .

The close agreement between both methods further supports the idea that the final black hole is indeed a Kerr black hole.

Acknowledgments

We thank Stanislav Babak, Ian Hinder, Luciano Rezzolla for useful discussions. CR and DP thank Rhodes University, and NTB thanks Max-Planck-Institut für Gravitationsphysik, for hospitality. This work was supported by international collaboration grants funded by National Research Foundation, South Africa, Bundesministerium für Bildung und Forschung, Germany, and DFG grant SFB/Transregio 7 "Gravitational Wave Astronomy". DP was supported by a grant from the