Gravitational wave extraction based on Cauchy-characteristic extraction and characteristic evolution

We implement a code to find the gravitational news at future null infinity by using data from a Cauchy code as boundary data for a characteristic code. This technique of Cauchy-characteristic Extraction (CCE) allows for the unambiguous extraction of gravitational waves from numerical simulations.


I. INTRODUCTION
The importance of gravitational waveform templates for gravitational wave detectors implies a need for accurate 3D numerical simulations of isolated sources such as binary black hole mergers. These simulations are often done with Cauchy codes based on a "3+1" slicing of spacetime. With the slicing conditions most commonly used in numerical simulations, two problems arise: artificial boundary conditions must be placed on the computational domain, and information such as the gravitational news cannot be extracted at future null infinity I + .
To avoid these problems two possible approaches have been suggested. One way is to evolve the entire spacetime. For example, a hyperboloidal slicing of the spacetime [1,2,3] allows information to propagate to I + in finite time. Conformal compactifications, such as those suggested by the conformal field equations [4,5] or employed in [6] allow I + to be located at a finite computational coordinate in a regular way. Another example is characteristic evolution based on a Bondi-Sachs line element , which gives a natural description of the radiation-region of spacetime extending to I + . Characteristic numerical codes have been used to study tail decay [7], critical phenomena [8,9,10,11,12], singularity structure [13,14,15,16] and fluid collapse [17,18,19], to list just a few examples.
Unfortunately, characteristic methods suffer from the appearance of caustics in the inner strong field region. The problem of caustics can be avoided by evolving the strong field region with a standard Cauchy slicing whilst using the characteristic approach for the exterior. This technique of Cauchy-characteristic matching (CCE) has proved successful in numerical evolution of the spherically symmetric Klein-Gordon-Einstein field equations [20], for 3D non-linear wave equations [21] and for the linearized harmonic Einstein system [22].
The second way of avoiding problems with standard Cauchy codes is a perturbative approach. The standard way of extracting gravitational waves from a numerical simulation is based on perturbation theory, either using the quadrupole formula or first order gauge invariant formalisms based on the work of Zerilli and Moncrief [23,24]. The vast majority of waveform templates are currently based upon these approaches. Numerical codes able to solve a generalization of the Zerilli equation to a three dimensional Cartesian coordinate system and to extract the gravitational signal are reported in [25,26,27,28]. Extraction of gravitational waves based on the Zerilli-Moncrief formalism in fully three-dimensional simulations are presented also in [29,30].
Perturbative methods have been also used to provide boundary conditions at the outer boundary of the Cauchy grid. This approach of Cauchy-perturbative matching has been implemented numerically in [31,32,33]. All those works are impressive, but discrepancies between the results of a perturbative approach and the full non-linear theory cannot be determined without solving the fully non-linear problem, although error estimates have been made [34].
Of these approaches, CCM has many appealing properties. The characteristic description of the exterior, particularly at I + , allows for a natural extraction of gravitational information such as the news. Matching to a standard Cauchy code in the interior implies that the methods employed in current 3D numerical simulations may immediately be used. As both Cauchy and characteristic approaches have been well tested and commonly used, CCM is a natural way of avoiding the potential problems of caustics in the characteristic region and artificial boundaries in the Cauchy region.
In this work we take a step towards the full CCM method by using a Cauchy code to provide boundary data for a characteristic code which propagates the solution to I + to extract the waveform. This procedure of Cauchy-characteristic extraction (CCE) allows the computation of gravitational waves in an unambiguous fashion. CCE was succesfully implemented in the quasispherical approximation in [35]; here we demonstrate it in the fully non-linear case.
In this work the outer region extending to I + is numerically evolved by the Pitt Null Code [20,35,36,37,38,39,40,41,42]. The link between the Cauchy and the characteristic modules is done by a non-linear 3D CCE algorithm [35,41,42,43,44]. At the outer edge of the characteristic grid the Bondi news is computed [37,39,45,46]. We have imported the CCE code into the Cactus computational infrastructure [47,48,49,50]. Within this infrastructure we have been able to use two separate Cauchy codes implementing the BSSN and harmonic 3 + 1 formulations of the Einstein equations. This allows us to test the robustness of the CCE approach.
We compare CCE with the Zerilli extraction technique [23,24,51] for both non-radiative and radiative spacetimes. The Zerilli formalism is based on the paper of Nagar and Rezzolla [52] which reviews and collects the relevant expressions related to the Regge-Wheeler and Zerilli equations for the odd and even-parity perturbations of a Schwarzschild spacetime. The conventions presented in their review are implemented in the Wave Extract code, included in the Cactus computational toolkit open source infrastructure [49].
For non-radiative spacetimes, the comparison gives results consistent with reference [34]: the CCE algorithm is O(∆ 2 ) accurate (where ∆ is the computational grid-step), while the accuracy of the perturbative (Zerilli) approach is O(r −2 ) where r is the radius of the wave extraction sphere. Since the signal is O(r −1 ), this implies an O(r −1 ) relative error in the Zerilli approach. As a result, it is shown in [34] that CCM is more efficient computationally in the sense that, as the desired error goes to zero, the amount of computation required for CCM becomes negligible compared to a pure Cauchy computation.
We further apply our algorithms to the study of the propagation of a linearized Teukolsky gravitational wave [27,53,54,55] and we compare the CCE and Zerilli wavesignal. For the case of a large extraction radius, we find that both methods give very good results and we demonstrate convergence to the analytical waveform of the Teukolsky solution for the CCE news. We show that the CCE waveform does not depend upon the extraction radius, which is a major advantage of the CCE method. A small extraction radius has no effect on the CCE waveform but it introduces errors in the Zerilli waveform.
In Sec. II we outline some general background and notation. In Sec. III we detail the transformation from the Cauchy slice to the characteristic slice. In Sec. IV A we give a brief summary of the characteristic evolution code. In Sec. IV B we describe how the news is computed at I + . Finally, in Sec. V we give two sets of tests in full 3D numerical relativity to validate the accuracy, convergence and robustness of the CCE algorithm. The first set contains four non-radiative tests, with no gravitational wave content. The second one is a radiative three-dimensional test. The results show that CCE is valid and accurate in both non-radiative and truly radiative situations.

II. NOTATION, GEOMETRY AND METRICS
The implementation of CCE described here follows previous descriptions of Cauchy-characteristic extraction and matching in the literature. Much of the work has been presented earlier [35,41,42,43,44]. Here we briefly outline the notation, geometry and metrics used.
The geometry is described by two separate foliations, neither of which cover the entire spacetime. The Cauchy foliation is described using a standard 3 + 1 ADM type metric [56], (2.1) Many different formulations can be used, given initial and boundary data, to evolve the 3-metric γ ij . In what follows we shall either consider the BSSN formulation [57,58] as implemented in [59] or the generalized harmonic formulation [60] as implemented in the Abigel code [61]. In both cases all interaction between the Cauchy and characteristic foliations will be performed in terms of the ADM metric, Eq. (2.1). The Cauchy slice does not extend to asymptotic infinity. Instead an artificial boundary is placed at |x i | = L i . Within this artificial boundary a world-tube Γ is constructed such that its intersection with any t = const. Cauchy slice is defined as a Cartesian sphere x 2 + y 2 + z 2 = R 2 , with angular coordinates labeled byỹÃ,Ã = (2, 3). The world-tube Γ is then used as the inner boundary of a characteristic foliation which uses the standard Bondi-Sachs metric [62,63] Here u labels the outgoing null hypersurfaces, y A =ỹÃ the angular coordinates (the null rays emanating from the world-tube), and r is a radial surface area distance. The angular metric h AB obeys the condition where q AB is the unit sphere metric. This coordinate system consistently covers the world-tube Γ and the exterior spacetime as long as it is free of caustics.
The free variables in the Bondi-Sachs metric are then V, β, U A and h AB . The physical interpretation of these variables is that h AB contains the 2 radiative degrees of freedom, e 2β measures the expansion of the nullcone, and V /r is the counterpart of the Newtonian gravitational potential. The 2+1 decomposition of the intrinsic metric on the r = const. world-tube identifies r 2 h AB as the intrinsic 2-metric of the u foliation, −U A as the shift vector and e 2β V /r as the square of the lapse.

III. THE CCE ALGORITHM
The crucial task of the CCE algorithm is to take Cauchy data given in the ADM form Eq. (2.1) in a neighborhood of the world-tube Γ and to transform it into boundary data for the Bondi-Sachs metric Eq. (2.2). Then the Bondi code can use the hypersurface equations to evolve the appropriate quantities out to I + so that the gravitational news may be extracted there. Much of the present version of the CCE algorithm has been presented in earlier work [35,41,42,43,44], so in this section we will give a brief description, highlighting a few new features.
In section II the world-tube was defined as a Cartesian sphere x 2 + y 2 + z 2 = R 2 , with angular coordinates labeled bỹ yÃ,Ã = (2, 3). In addition to the angular coordinates, we set u = t on the world-tube and choose the fourth coordinate, λ to be the affine parameter along the radial direction, with λ |Γ = 0. The characteristic cones are constructed such that λ is the future oriented, outgoing null direction normal to the foliation of Γ. (In order to avoid a singular Jacobian for the Cauchy-characteristic coordinate transformation, we require that the world-tube Γ be timelike.) The affine parameter λ is used because the world-tube is not constructed to be a surface of constant Bondi r.
The choice of the angular coordinates is determined following [64] by the use of two stereographic patches. The use of two patches avoids numerical difficulties in taking derivatives near the poles when standard spherical coordinates are used. The use of only two patches may not give the most accurate numerical results; various ways of discretizing the 2-sphere on multiple patches are discussed in [65]. The two patches are centered around the North and South poles with the stereographic coordinates related to the usual spherical coordinates (θ, φ) by and We also introduce the complex vector on the sphere q A = (P/2)(δ A 2 + iδ A 3 ) and its co-vector q A = (2/P )(δ 2 A + iδ 3 A ), with P = 1 + ξξ = 1 + q 2 + p 2 . The orthogonality conditionq A q A = 2 is satisfied by construction. The unit sphere metric corresponding to these coordinates is The angular subspace in the Bondi code is treated by use of the eth formalism and spin-weighted quantities. (See [39,64] for details.) The CCE algorithm starts by interpolating the Cauchy metric γ ij , lapse α and shift β i and their spatial derivatives onto the world-tube. Time derivatives are computed via backwards finite-differencing done along Γ, e.g., at t = t N we write Knowledge of the Cauchy metric and its 4-derivative is enough to compute the affine metricηαβ as a Taylor series expansion around the world-tubeηαβ Next a second null coordinate system (the Bondi-Sachs system) is introduced The fourth coordinate y 4 = r is a surface area coordinate, defined by The Bondi-Sachs metric on the extraction world-tube can be computed as Note that the metric on the sphere is unchanged by this coordinate transformation, i.e., η AB =ηÃB. Therefore one only needs to work with the Jacobian components that correspond to derivatives of r.
In terms of q A ,q A the two dimensional metric η AB can be encoded into the metric functions The determinant condition Eq. (2.3) translates into With h AB symmetric and of fixed determinant, there are only two degrees of freedom in the angular metric that are encoded into the complex function J.
This quantity is a measure of the expansion of the light rays as they propagate outwards. The CCE and the Bondi codes have been implemented with the assumption that r ,λ > 0 (or that β is real). The radial-angular components η rA can be represented by In addition to these quantities, in [40] the auxiliary variables have been introduced to eliminate the need to explicitly use second angular derivatives in the Bondi evolution code. The required boundary data is J, β, U, ∂ r U, W, ν, k, and B. (See Sec. IV A.) Notice that once the Bondi-Sachs metric is known on the worldtube one can only obtain J, β, U, W . In order to provide the rest of the necessary boundary data we need the radial derivative of the Bondi-Sachs metric Withηαβ ,λ already known, the only non-trivial parts in Eq.(3.16) are the Jacobian terms These depend on the second derivatives of the Cauchy metric. In order to avoid possible numerical problems caused by interpolating second derivatives onto the world-tube, we calculate r ,λÃ by taking centered derivatives of r ,λ on the world-tube and we calculate r ,λu by backwards differencing in time along the world-tube. The remaining term r ,λλ is calculated using the identity 18) and the characteristic equation (The right hand side of Eq. (3.19) can be computed in terms of J ,r = J ,λ /r ,λ , which, in turn, can be computed in terms of already known quantities.) Knowledge of the radial derivative of the Bondi metric is important not only for obtaining (∂ r U ) |Γ but also because the gridstructure of the Bondi code is based on the radial coordinate r, and the extraction world-tube will not, in general, coincide with any of these radial grid-points. We need to use, therefore, Taylor series expansions to fill the Bondi gridpoints surrounding Γ with the necessary boundary data. We write, e.g., (3.20) Another problem is the need to provide the auxiliary angular variables, B, k and ν on the world-tube. These have been defined as the ð-derivatives of Bondi fields in the y α frame, while taking angular derivatives on the world-tube amounts to computing ð derivatives in theỹα frame. The correction term between the two frames is The quantities J,ðJ and J ,λ are known from the interpolated Cauchy data, while J ,λλ is not computed in CCE. Thus one can compute ðJ but not ðJ ,λ . As a consequence we cannot use a simple Taylor series expansion to place ν on the Bondi grid to O(λ 2 ) accuracy. The solution to this problem is to first place J on Bondi grid-points surrounding the world-tube, then compute ðJ on those points (i.e., compute angular derivatives in the Bondi-Sachs frame), and then calculate ∂ λ ðJ on the world-tube by use of the neighboring Bondi grid values of ðJ. With ν = ðJ and its radial derivative known on the world-tube, we can then use the standard O(λ 2 ) expansion to provide ν at Bondi gridpoints surrounding the world-tube. This way we make maximal use of the Cauchy data as interpolated onto the world-tube and minimal use of the finite difference ð algorithm that has its own discontinuous O(∆ 2 ) error at the edges of the stereographic patches. A similar approach is used for k = ðK and B = ðβ.

A. The Bondi code
The inner workings of the Bondi code had been described in detail elsewhere (see [20,35,36,37,38,39,40,41,42]) so here we give only a brief overview of the algorithm. As already stated in Sec. III, the variables of the code are J, β, U, W as well as ν, k, B. Out of these the equation for J is the only one to contain a time derivative. For this reason J is updated via an evolution stencil that involves two time-levels. The rest of the variables are integrated radially from the world-tube out to I + . The integration constants are set by CCE. All of these radial integration equations are first differential order in r except for the U equation which contains U ,rr . For this reason integrating U requires two constants, which explains the need to provide, at the world-tube, U as well as U ,r .

B. The News algorithm
The calculation of the Bondi news function on I + is based on the algorithm developed in [39] with the modifications introduced in [45,46] (an alternative calculation for the news was recently introduced in [66]). Here we present an overview of the algorithm.
The Bondi-Sachs metric Eq. (2.2) in a neighborhood of I + in (u, x A , l = 1/r) coordinates (after multiplying by a conformal factor l 2 ) has the form where the tilde denotes a quantity defined with respect to inertial observers. Note thatm (AmB) = q AB as required. We define a complex vector F a , analogous tom a (i.e. F (AF B) = H AB ), adapted to the coordinates used in the characteristic evolution code via where where q A is the dyad defined in Sec. III, J 0 = q A q B H AB /2, and K 0 = q AqB H AB /2. F a andm a are related bỹ The γñ a term will not enter the news calculation.
To calculate the Bondi news one needs to evolve two scalar quantities δ, the phase factor in Eq. (4.7), and the conformal factor ω, as well the relations ξ(ξ B ) between the angular coordinates used in the characteristic evolution and inertial angular coordinates, and u B (u, ξ B ) between the inertial time slicing and the time slicing of the characteristic evolution code. Then δ, ξ(ξ B ), and u B (u, ξ B ) are evolved using the following ODE's along the null generators of I + where U 0 = q A L A . Also ω is evolved using the PDE The Bondi news function up to a phase factor of e −2iδ is given by where D A is the covariant derivative with respect to H AB . The Bondi news is calculated in three steps. First Eq. (4.12) is evaluated, ignoring the e −2iδ phase factor, as a function of the evolution coordinates (u, ξ). Then using the relation ξ(ξ B ), the news is interpolated onto a fixed inertial angular grid (i.e. N (u, ξ B )) and multiplied by the phase factor e −2iδ (which is only known on the inertial grid). Finally, in post-processing, the news is interpolated in time onto fixed inertial time slices (i.e .  N (u B , ξ B )). Once the news is obtained in inertial coordinates it can be decomposed into spin-weighted spherical harmonics.

V. TESTS
We apply the algorithm described above for two sets of tests: non-radiative spacetimes and radiative spacetimes. In the first, non-radiative set of tests (section V A), we analyze Minkowsy and Schwazchild spacetimes. The Minkowski (Sec. V A 1) or small perturbation of Minkowski (Sec. V A 2) tests are used to show the stability of the code and the errors due to transforming between the different coordinate systems and sets of variables.
The Schwazchild tests describe a static spherically symmetric black hole in a centered frame (Sec. V A 3) and in an oscillating frame (Sec. V A 4). These tests indicate the accuracy of the code in non-trivial spacetimes.
The last test (Sec. V B) is a standard radiative test: namely a linearized Teukolsky wave, and indicates that CCE is a valid and accurate method in truly radiative situations.
In the non-radiative cases, the extracted gravitational wave signal should vanish identically. As this is true uniformly for all points at I + we simplify the computation for the norms of the gravitational wave signal. When the results of the CCE are shown the norm used is the simple L 2 norm over all grid points at I + , without any weighting by the area element. In the radiative case, we plot the real part of the extracted news waveform (∂ t h + ), versus time.
The Zerilli Moncrief formalism (see [23,24,51] for the implementation used here) approximates the spacetime as a linearized perturbation of a spherically symmetric spacetime. Such perturbations are decomposed in terms of spherical harmonics which are then expressed in terms of some basis set. An example using the Regge-Wheeler set is where the even parity metric perturbations are decomposed as Once the basis functions such as G, K, h 1 , H 2 and b are found from the decomposition, the even parity master function Q + can be computed from Equivalent definitions in terms of other basis functions, alternative conventions and expressions for the odd parity sector can be found in [52].
Using this notation where Q × gives the odd parity master function, the asymptotic value for the wave forms is For the non-radiative cases, when we compute the wavesignal from the Zerilli approximation, we use the norm where dΩ is the element of solid angle. The last equality Eq. (5.5) follows from the orthogonality of the spin-weighted tensor spherical harmonics. We expect h + − ih × to vanish for all angles so this norm should be identically zero. In what follows we shall look at the spherical harmonics up to ℓ = m = 7 although in practice there is little contribution from the higher modes. For the Teukolsky wave used here we expect all terms except for the Q + 2,0 term in Eq. (5.3) to be negligible asymptotically, which leads to a simple but tedious computation of the wavesignal.

Minkowski Flat Spacetime
Minkowski space in standard coordinates is evolved in the harmonic Abigel code, with the metric given analytically on all Cauchy slices. The domain extends between x i ∈ [−10, 10], the extraction world-tube is placed at r = 7, and the simulation is run until t = 10. The coarsest simulation used 50 3 points for the Cauchy grid and 35 2 × 31 points for the characteristic grid. Simulations to test convergence scaled all grids by factors of 2. For this test, extraction using the Zerilli method gives results that are identically zero, as expected.
This test indicates the level of numerical round off error in transforming from the Cauchy variables to the Bondi variables on the world-tube. As shown in Fig. 1, the extracted news is small in all cases but the level of truncation error increases linearly with the number of points on the extraction world-tube. Since the construction of the boundary data for the CCE code involves derivatives of the Cauchy metric, this sensitivity to noise in the Cauchy data is expected. However, given the extremely low amplitude of the noise, it is unlikely that this will cause problems in practical simulations.

Random Perturbations around Minkowski Flat Spacetime
This test is in the spirit of the robust stability tests of [22,67,68]. The stable evolution and extraction of white noise initial data with no frequency dependent growth of the wavesignal is a good indication that the combined evolution and extraction codes are stable. The domain extends between x i ∈ [−10, 10], the extraction world-tube is placed at r = 7, and the simulation is run until t = 100. The coarsest Cauchy grid has 51 3 points and the coarsest characteristic grid 35 2 × 31 points.
The results in Fig. 2 show that the error is independent of the resolution of both Cauchy and characteristic grids. This is strong evidence that the CCE code is stable against small perturbations that in a practical run would be induced by numerical error. The wavesignal from the Zerilli extraction is about an order of magnitude smaller. This is to be expected as the Zerilli extraction approach uses the field values. The error in the CCE news has a jump at the beginning. This arises because on the initial characteristic data is set to zero so that it takes some time for the noise to build up on the outgoing null cone.

Schwarzchild Black Hole in a "centered" frame
To test a simple black hole spacetime we use a Schwarzschild black hole in ingoing Eddington-Finklestein coordinates (t,x i ), with the line element This is manifestly static in these coordinates. The spacetime is evolved using the BSSN code and excision methods described in [69]. The coarsest Cauchy grid has 29 3 points with ∆x i = 0.4M , whilst the characteristic grid has 35 2 × 31 points. The world-tube is at r = 7M . In the Cauchy evolution domain octant symmetry is used. The evolution is only performed for a short time (to t = 100M ). Over this timescale we see second order convergence for the CCE news until t ≈ 7M and second order convergence in the waveform from Zerilli extraction until t ≈ 20M , as seen in Fig. 3. By varying the location of the world-tube or the Zerilli extraction sphere we can see that the errors come from a variety of locations.
The difference between the times at which the two extraction methods lose convergence is probably due to the greater differencing error seen in the Zerilli extraction. At early times (where finite differencing error dominates) the absolute value of the error is larger for Zerilli extraction by a factor ≈ 4, as seen in Fig. 3. At late times (where outer Cauchy boundary errors dominate) the error for Zerilli extraction is of the same order as for the CCE extraction.
The non-convergent errors are probably caused by the boundary conditions on the Cauchy grid which do not satisfy the constraints. As in [69] we are simply applying Sommerfeld type boundary conditions to all fields. This condition does not a The news displays very slow growth that is independent of the frequency of the initial data. This is a good indication that the CCE method is stable. The right hand figure shows that the norm of the waveform extracted from the Cauchy grid using the Zerilli method is about one order of magnitude smaller.
priori satisfy the constraints and is not known to be well-posed. Thus we might expect errors to be induced by the use of these boundary conditions. It is likely that constraint satisfying boundary conditions or Cauchy-characteristic matching would solve or at least greatly reduce this problem. Due to the complicated pattern of the reflected errors arising from a boundary condition on a cubical surface, we have been unable to determine exactly how this error depends on the location of the outer Cauchy boundary. FIG. 3: The scaled L2 norm of the CCE news for a static spherically symmetric black hole where the Cauchy slice is evolved using the BSSN formalism. At early times, when the world-tube is causally disconnected from the outer boundary of the Cauchy slice, the CCE news converges to zero at second order as it should. Later, when the outer boundary is causally connected to the world-tube, deviations from second order convergence appear, indicating the effect of the Cauchy boundary conditions on the extracted wavesignal. The left panel shows the CCE news extracted at I + , whilst the right panel shows the "norm" of the wavesignal extracted by the Zerilli method at r = 7.

Schwarzchild Black Hole in an oscillating frame
We use the line element for the standard Schwarzschild black hole in ingoing Eddington-Finklestein coordinates, given in Eq. (5.6). The moving coordinate frame (t, x i ) is given by where B i are parameters specifying the velocity and b(t) is a simple periodic function turned on after some time t > t 0 ; here we use The CCE news converges to zero as it should. However, Zerilli extraction does not give a signal that converges to zero as grid resolution is refined. Instead, an erroneous non-trivial signal is computed.
FIG. 4: The scaled L2 norm of the CCE news for a static spherically symmetric black hole in an oscillating frame. The CCE news converges to zero at second order as it should.

Theoretical part
The Teukolsky solution [27,53,55] to the linearized Einstein equation is a weak gravitational wave propagating through space, and represents one of the most valuable standard test cases for numerical codes that implement wave extraction techniques.
The general form of the spacetime metric is ds 2 = −dt 2 + (1 + Af rr )dr 2 + 2Bf rθ rdrdθ + 2Bf rφ r sin θdrdφ + (1 + Cf θθ )r 2 dθ 2 + 2(A − 2C)f θφ r 2 sinθdθdφ + (1 + Cf (1) φφ + Af where the angular functions f ij , corresponding to an l = 2, m = 0, spin − weight = 2 spherical harmonic, are (5.10) The functions A, B and C are given in terms of a free generating function F (x) = Axe −x 2 /λ 2 /λ 2 , where A is the amplitude and λ determines the width of the wave, by Here, x = t − r corresponds to an outgoing wave. To obtain the ingoing wave solution coresponding to x = t + r, we change the sign in front of all the terms with odd numbers of F . We consider a superposition of ingoing t + r and outgoing t − r waves [44], centered at the origin of the coordinate system at t = 0, which provides a moment of time symmetry. This solution, after is linearized in amplitude, is implemented in the harmonic Abigel code.
We further give the analytical form of the Teukolsky wave signal. The real part of the Bondi news function is proportional to the time derivative of the "plus" polarization mode of the gravitational wave signal. The time derivative of the CCE news satisfiesṄ = lim r→∞ rΨ 4 (5.13) where Ψ 4 is the Newman-Penrose component of the Weyl tensor Ψ 4 = −C µρντ n µmρ n νmτ (5.14) with n ν an ingoing null vector, and m ρ a complex unit vector, oriented in the angular directions. The connection between m ρ and the dyad vector q A on the unit sphere is Moreover, n µ ∂ µ = ∂ x . A straightforward computation starting from Eq. (5.13), giveṡ At infinity, in the linearized regime, only the O(1/r) terms in the Teukolsky functions A, B and C contribute to the news. The term that determines the gravitational wave signal is d 4 x F . We finḋ Alternatingly, computation of the linearized expression for the Bondi news from Eq. (4.12), gives

C. Numerical part
We give the metric specified by the Teukolsky solution analytically on the Cauchy grid. Then data is extracted at the worldtube and evolved by the CCE code. We carry out a series of simulations, varying both the location of the characteristic world-tube and the radius of the Zerilli sphere. We use 80 3 , 120 3 , and 160 3 gridpoints for the Cartesian Cauchy grid and 60 2 ×80, 90 2 ×120, and 120 2 ×160 gridpoints for the null grid. The domain extends between x i ∈ [− 15,15] and the simulations are run until t = 30.
We study the dependence of the signal with the amplitude and conclude that CCE code resolves correctly amplitudes of A ≥ 10 −8 . At smaller amplitudes, clean convergence behavior is contaminated by round-off error. We show results only for A = 10 −5 , λ = 1. We study also the dependence of the wave signal with the world-tube radius and conclude that the accuracy of the computed CCE news is preserved even for radii as small as r = 5. Fig. 6 shows the convergence of the CCE news to the analytical solution, for a world-tube radius of r = 10. The convergence ratec r of the CCE news to the analytical data, is given bỹ The convergence rate of the computed CCE news to the analytic value Eq. (5.18) at t − r = 0, corresponding to the peak of the radiated signal, isc Hence, the CCE wavesignal is second order convergent and independent of the world-tube radius, as expected. Fig. 7 shows that for small extraction radii, the Zerilli waveform has a slight asymmetry, which is caused by the dependence of the Zerilli formalism upon the extraction radius. We have to increase the extraction radius in order to decrease this error. Because the error does not fall off sufficiently fast with radius in the near zone, where we can realistically carry out the simulation, we instead analytically compute the Zerilli news at the extraction radius r = 300. Fig. 8 demonstrates that the agreement between the computed CCE news at small radii and the analytical Zerilli news at big radii, is very good.

VI. CONCLUSIONS
We have demonstrated the accuracy of Cauchy-characteristic extraction in 3D numerical relativity. The interface at the extraction world-tube does not introduce significant error with either the BSSN [59] or harmonic [22,61] formulations.  The CCE news is stable against small perturbations on the Cauchy grid, as shown in Sec. V A 2, and the level of truncation error caused by transforming from the Cauchy to the characteristic variables is small, as shown in Sec. V A 1. In a non-trivial black hole spacetime, CCE performs as expected. Second order convergence is found in the moving Schwarzhild black hole test in Sec. V A 4 and in the early time behavior of the BSSN test in Sec. V A 3. However, as seen in the late time behavior in Sec. V A 3, the effect of improper BSSN outer boundary conditions on the Cauchy grid is clearly visible in the extracted CCE news.
In the linearized Teukolsky wave test in Sec. V B, we demonstrate that the computed CCE news converges to the analytic Teukolsky waveform and does not depend on the world-tube radius. The accuracy of the CCE news is preserved even for small radii, while the Zerilli news is affected by near zone error. When the extraction radius is sufficiently large, both Zerilli and CCE give excellent results. The advantage of CCE over the Zerilli method is that the extracted CCE news does not depend on the extraction radius. Cauchy-characteristic extraction can be seen either as a first step towards a full Cauchy-characteristic matching code or as an improved gravitational wave extraction method in its own right. The results of this paper show that as a stand alone wave extraction method, it produces the correct results in situations where other methods, such as Zerilli extraction, may fail. Such a situation is seen in Sec. V A 4, where the Zerilli method fails to converge as a result of the pure gauge motion of the Schwarzchild metric.
Finally, these tests also show the need for improvement of the current implementation of the code. Any angular dependence in the solution whether through genuine physics, gauge effects, or due to the imposition of boundary conditions on the cubical Cauchy boundary, leads to short wavelenght error that is poorly resolved by the CCE code. This is particularly noticeable for features in the region where the stereographic patches overlap. Improvements in this area of the implementation will enhance the accuracy of the code in extracting the gravitational news from astrophysical simulations requiring high resolution. To this end, we are investigating the use of different multiple patch implementations such as [65].