Binary Black Hole Waveform Extraction at Null Infinity
Abstract
In this work, we present a work in progress towards an efficient and economical computational module which interfaces between Cauchy and characteristic evolution codes. Our goal is to provide a standardized waveform extraction tool for the numerical relativity community which will allow CCE to be readily applied to a generic Cauchy code. The tool provides a means of unambiguous comparison between the waveforms generated by evolution codes based upon different formulations of the Einstein equations and different numerical approximation.
pacs:
04.20Ex, 04.25Dm, 04.25Nx, 04.70BwI Introduction
The study of gravitational waves will brighten unexplored features of our universe that are otherwise invisible to conventional astronomy and will increase our knowledge about the very nature of time and space (1). Gravitational wave detectors are already operating, and results from the first LIGO and Virgo collaboration were recently published in Nature (2). The signal predicted by numerical relativity will provide a template bank used for filtering the noise, indispensable to the success of gravitational wave detectors such as LIGO, Virgo, and LISA. The current sensitivity levels of the detectors will be improved substantially in nextgeneration detection estimated by 2015. Although existing simulations are sufficiently accurate for populating the parameter space in current searches of groundbased detectors, the new generation of advanced detectors will be 1015 times more sensitive by 2015.
Ideally, the emitted gravitational wave signature should be extracted at spatial or null infinity. However, most present codes impose artificial, finite outer boundaries and are performing the waveform extraction at finite radius. This method introduces systematic errors, especially for higher modes, which is the main obstacle in reaching the desired accuracy. With the exception of Pretorius (3), who uses coordinates that compactify spatial infinity, all the other codes use a computational domain with a finite outer boundary and Sommerfeldlike approximate outer boundary conditions must be imposed, which introduce errors in the computation of gravitational waves. The choice of proper boundary conditions is complicated by gauge freedom and constraint preservation (4). The emitted gravitational wave signature is calculated at finite distance, using either the NewmanPenrose Weyl scalar (5), or the odd and even parity functions , in the ZerilliMoncrief formalism (6). The strain of the wave used in detection is obtained performing one time integration from the ZerilliMoncrief multipoles, and two time integrations from the NewmanPenrose curvature. The waveform is affected by gauge ambiguities which are magnified by the integration (7).
Cauchycharacteristic extraction (CCE) (8), which is one of the pieces of the CCM strategy (9), offers a means to avoid the error introduced by extraction at a finite worldtube. In CCE, the inner worldtube data supplied by the Cauchy evolution is used as boundary data for a characteristic evolution to future null infinity , where the waveform can be unambiguously computed by geometric methods. This characteristic initialboundary value problem based upon a timelike worldtube (10) has been implemented as a mature evolution code, the PITT null code (11); (12), which incorporates a Penrose compactification of the spacetime. By itself, CCE does not use the characteristic evolution to inject outer boundary data for the Cauchy evolution, which can be a source of instability in full CCM.
The PITT code has been tested to be second order convergent in a wide range of testbeds extending from the perturbative regime (13) to highly nonlinear single black hole spacetimes (12). However, in cases which require high resolution, such as the inspiral of matter into a black hole, the error in CCE has been a troublesome factor in the postprocessing phase (14). This has motivated a recent project (15) to increase the accuracy of the PITT code. Other results achieved with previous versions of the PITT have been recently reported (16); (17). Recently, the code underwent major improvements and corrections to previous versions, to improve accuracy and convergence (18).
Here we test this improved version of CCE on a realistic application involving a Cauchy evolution of the inspiral and merger of two equal mass nonspinning black holes. We use the same code specifications described in (15) except that the accuracy of angular derivatives has been increased to a 4th order finite difference approximation. The results presented here are a work in progress towards our goal to develop CCE as a reliable and accurate waveform extraction tool for the numerical relativity community. This paper addresses the first two objectives:

To create a robust and flexible interface between a binary blackhole Cauchy evolution code and a characteristic code for wave extraction at infinity.

To prove the robustness of the interface by performing precise computations of gravitational waveforms at infinity from binary blackhole, using this Cauchycharacteristic extraction approach.
We construct an interface that takes the Cartesian data from a Cauchy evolution and converts it into boundary data on a spherical grid for the characteristic evolution. The data are evolved to future null infinity, where it is used to compute the gravitational waveform. The flexibility of the interface is due to implementation of a spectral decomposition of data. This implementation has been tested with a realistic application involving a binary blackhole inspiral. In Sec. II we review the formalism underlying CCE, including enough details of the patching, evolution and extraction, to make clear the difficulties underlying the calculation of an accurate waveform at . In Sec. III we briefly describe the initial data for Cauchy and characteristic evolution. In Sec. IV, we present the details of the CCE interface which allows the data from a Cauchy evolution to be used as boundary data on an inner worltube for a characteristic evolution to , where the waveform is extracted. In Sec. V, we test the CCE interface by extracting the waveform from a Cauchy evolution of a binary blackhole inspiral and merger, and by comparing it to the waveform obtained by other standard method in current practice.
Ii The CCE Formalism
ii.1 CauchyCharacteristic Patching
Characteristic data are provided by the Cauchy evolution on a worldtube , free initial data being given on the initial null hypersurface , which sets the metric on the entire initial cone (fig. 1).
The metric data from a Cauchy evolution are interpolated onto a timelike inner worldtube to extract the boundary data for the characteristic evolution. The characteristic evolution is embedded into the Cauchy evolution and is extending to future null infinity , where the waveform can be unambiguously computed using the geometric methods developed by Bondi et al (19), Sachs (20) and Penrose (21). The extraction process involves carrying out the complicated Jacobian transformation between the Cartesian coordinates used in the Cauchy evolution and the spherical null coordinates used in the characteristic evolution (the full details are given in (22).)
ii.2 Characteristic Evolution
The characteristic formalism is based upon a family of outgoing null hypersurfaces, emanating from some inner worldtube, which extend to infinity where they foliate into spherical slices. We let label these hypersurfaces, be angular coordinates which label the null rays and be a surface area coordinate. (fig. 2).
In the resulting coordinates, the metric takes the BondiSachs form (19); (20)
(1)  
where is the BondiSachs conformal 2metric with . The code introduces an auxiliary unit sphere metric , with associated complex dyad satisfying . For a general BondiSachs metric, the full nonlinear is uniquely determined by the dyad component , since the other dyad component is constrained by the determinant condition . The spherically symmetric case characterized by . We introduce the spinweighted fields and , where
(2) 
as well as the (complex differential) operators and . Refer to (23); (8) for further details regarding numerical implementation. The auxiliary variables
(3) 
are also introduced to eliminate all second angular derivatives. In certain applications this has been found to give rise to increased accuracy by suppressing short wavelength error (24).
In this formalism, the Einstein equations decompose into hypersurface equations, evolution equations and conservation conditions on the inner worldtube. As described in more detail in (25); (26), the hypersurface equations take the form
(4)  
(5)  
(6)  
(7) 
where (23)
(8) 
is the curvature scalar of the 2metric . Those equations have a hierarchical structure in such that the right hand sides, e..g. only depend upon previous variables and their derivatives intrinsic to the hypersurface.
The evolution equation takes the form
(9) 
where, , , , and are nonlinear terms which vanish for spherical symmetry. Expressions for these terms as complex spinweighted fields and a discussion of the conservation conditions are given in (8).
The characteristic Einstein equations are evolved in a domain between an inner radial boundary at the interior worldtube, and an outer boundary at future null infinity. The characteristic evolution code implements this formalism as an explicit finite difference scheme, based upon the compactified radial coordinate
(10) 
so that at . Here is a parameter based upon the extraction worldtube, which in the CCE module is chosen as the radius of the extraction worldtube, as determined by in terms of the Cartesian coordiates used in the Cauchy evolution code. The boundary data for , , , , and on the worldtube supply the integration constants for a radial numerical integration of the hypersurface Einstein equations. The finite difference scheme for integrating the hypersurface and evolution equations is based on the marching equation for a spherically symmetric scalar field :
(11) 
where the point N is the ”new” point in the evolution scheme, and is defined by the spherically symmetric version of the BondiSachs metric given above. The evolution scheme in the full gravitational case used to determine the metric at the next point on the null hypersurfaces is modeled after this example (see (12); (24) for details).
ii.3 Gravitational Radiation Calculation
The theoretical derivation of the waveform at infinity is carried out in terms of an inverse surfacearea coordinate , where at . In the resulting coordinates, the physical spacetime metric (1) has the conformal compactification , where is smooth at and takes the form (10)
(12) 
As described in (15), the Bondi news function and the NewmanPenrose Weyl tensor component which describe the waveform are both determined by the asymptotic limit at of the tensor field
(13) 
constructed from the leading coefficients in an expansion of the metric in powers of
(14)  
(15)  
(16)  
(17) 
where and are the 2dimensional curvature scalar and covariant derivative associated with .
The expansion coefficients , , and (all functions of and ) completely determine the radiation field. Before the gravitational radiation is calculated from the metric in the neighborhood of , it is necessary to determine the conformal factor relating to a unit sphere metric , i.e. to an inertial conformal Bondi frame (10) satisfying
(18) 
The news function is directly computed by the code in terms of the computational coordinates , as opposed to the inertial coordinates on corresponding to an idealized distant observatory. The transformation to inertial coordinates proceeds first by introducing the conformally rescaled metric in which the crosssections of have unit sphere geometry, in accord with (18). Then the rescaled null vector is the generator of time translations on , i.e. . The inertial coordinates thus satisfy the propagation equations
(19) 
where in terms of the computational coordinates. The inertial coordinates are obtained by integrating (19), thus establishing a second pair of stereographic grid patches corresponding to . Then the news function is transformed into .
The Bondi news function is given by (20),
(20) 
where is the Lie derivative with respect to . The NewmanPenrose Weyl tensor component is given by (21)
(21) 
In the inertial Bondi coordinates, the expression for the news function (20) reduces to the simple form
(22) 
and (21) reduces to the single term
(23) 
This is related to the expression for the news function in inertial Bondi coordinates by
(24) 
Equation (24) holds true in the linearized approximation of the Einstein equations. In the nonlinear case, the full expression for news and must be used in the code. This introduces additional challenges to numerical accuracy due to high order angular derivatives of and large number of terms.
Iii Initial Data
iii.1 Initial Cauchy Data
For the Cauchy evolution we used the LazEv code (27); (28) along with the Cactus framework cac () and Carpet (30) mesh refinement driver. LazEV is an eighthorderaccurate finitedifference code based upon the BaumgarteShapiroShibataNakamura (BSSN) formulation (31); (32) of Einstein’s equations, which deals with the internal singularities by the moving puncture approach (27); (33). Our simulation used 9 levels of refinement with finest resolution of , and outer Cauchy boundary at . The initial data consisted of a close quasicircular blackhole binary with orbital frequency , leading to more than a complete orbit before merger (see (34)). We output the metric data on the extraction worldtube every .
iii.2 Initial Characteristic data
The initial data for the characteristic evolution consist of the values of on the initial hypersurface . One way of supressing incoming radiation in the data would be to set the NewmanPenrose Weyl tensor component on the initial null hypersurface. For a perturbation of the Schwarzschild metric, this condition implies no incoming radiation in the linearized approximation. However, in order to avoid shocks arising from incompatibility with the Cauchy data on the extraction world tube (with given by 10), we also need to require that and are continuous. In the linearized approximation, the condition that implies that . The combination of those requirements leads to , which would imply that at . For technical simplicity we avoid this complication by initializing according to
(25) 
which matches the Cauchy data and the derivatives at and is consistent with asymptotic flatness. Since this choice of vanishes at infinity, the initial slice of has a unit sphere metric so that the conformal factor has the simple initialization .
Iv Computational interface
We have designed an interface that takes Cartesian grid data from a Cauchy evolution and converts it into boundary data for characteristic evolution on a spherical grid extending to . We treat each component of the Cauchy metric as a scalar function in the Cartesian coordinates which are used in the evolution.
In order to make the interface as flexible as possible for future development as a community tool for waveform extraction, we have based it upon a spectral decomposition of the Cauchy data in the region between two world tubes or radii and , where is the Cartesian coordinate radius. Then at a given time , we decompose in terms of Tchebychev polynomials of the second kind and spherical harmonics , where are related to in the standard way. The Tchebychev polynomials are conventionally defined as functions on the interval . Here we map them to the interval by the transformation
where the extraction shell thickness is determined by the number of Tchebychev polynomials used. (In tests of binary black holes with mass M we use a relatively small range , a larger value of would be needed if the range were expanded). Thus, for , we expand
(26) 
For the applications to waveform extraction given in this paper, it is sufficient to consider , where , and , where . The coefficients then allow us to reconstruct a spherical harmonic decomposition of each component of the Cauchy metric on the extraction worldtube , i.e.
(27) 
This decomposition is carried out at a sequence of Cauchy time steps , where is chosen to be much smaller than the characteristic time scale of the problem but, for purposes of economy, larger than the time step used for the Cauchy evolution. A fifthorder polynomial interpolation is carried out locally over the to provide characteristic boundary data at any time in analytic form.
The extraction module also requires the derivatives and at the extraction worldtube. The derivative is constructed by a fourthorderaccurate finitedifference stencil using the surrounding Cauchy times . The derivative is obtained analytically, at each time level , by differentiation of the Tchebychev polynomials.
The spherical harmonic interpolator from the Cartesian to the spherical coordinates is part of the extraction module, but its resolution is controlled by the Cauchy evolution.
The stereographic coordinates used to label the outgoing null rays in the Bondi metric are matched to the spherical coordinates induced by the Cartesian Cauchy coordinates on the extraction worldtube by a standard transformation, using the conventions in (23). The value of the surfacearea coordinate in the BondiSachs metric is obtained on the extraction worldtube from the 2determinant of the Cartesian metric on the surfaces . As a result the radius of the Bondi coordinate on the extraction worldtube. The metric has to be calculated at a common value of the surface coordinate , because the original Cauchy extraction was at constant R. In order to make this calculation possible, the transformation from Cartesian coordinates to BondiSachs coordinates is carried out via an intermediate Sachs coordinate system (20) where is an affine parameter along the outgoing null rays. The affine freedom allows us to set on the extraction worldtube . After carrying out the Jacobian transformation from to , the Cartesian metric and its first derivatives at the extraction worldtube provide a firstorder Taylor expansion in (about ) of the null metric in Sachs coordinates. The corresponding Taylor expansion of the metric in BondiSachs coordinates then follows from the computed value of and at , which are obtained from the 2determinant of the Cartesian metric. In order to obtain a firstorder Taylor expansion for the Bondi metric variable , the hypersurface equation (4) must be used to evaluate at the extraction worldtube. All other metric variables are then initialized consistent with second order accuracy. Taylor expansions are also needed to start up the radial integration equations for the auxiliary variables (3) used to convert angular derivatives to firstorder form. These expansions are obtained from applying the operator to the Taylor expansion of the underlying metric. This is a complicated process because the operator intrinsic to the extraction worldtube is not the same as the operator intrinsic to the Bondi spheres (see (18) for a discussion). The low order intermediate Taylor expansions limit the accuracy of the result. A new approach that avoids entirely the use of the Taylor expansion and gives better accuracy is presented in (18). In the original approach, used for this paper, the resulting Taylor expansion of the evolution variables is used to fill the points of the BondSachs grid to start the integration of the characteristic hypersurface and evolution equations (4)  (9). The integration proceeds from the extraction world tube to on a radial grid based upon the compactified coordinate (10).
Domain of dependence considerations place a constraint between the characteristic time step and the size of the characteristic grid analogous to the CFL condition for the Cauchy evolution. For an estimate, consider the Minkowski space case with the conformally rescaled metric
(28) 
where the unit sphere metric takes the form
(29) 
The past light cone is determined by
(30) 
For typical characteristic grid parameters, , the resulting restriction is
(31) 
For a Cauchy simulation of a binary blackhole system of total mass with timestep (sufficient to describe the typical frequencies of a binary system), (31) leads to
(32) 
for the choice of characteristic timestep . The corresponding number of radial gridpoints must roughly satisfy . This places no limit of practical concern on the resolution of the characteristic evolution even for the small extraction radius . Thus, for purposes of CCE, there are no demanding CFL restrictions.
The interface was debugged and calibrated using the analytic Schwarzschild metric in KerrSchild coordinates ,
(33) 
where .
V Results
We present results for the characteristic extracted waveform either in terms of , related to the Bondi news by in the linearized regime, or, when comparing to the perturbative waveform, in terms of the NewmanPenrose component . The relationship between the Cauchy and the characteristic waveforms is: . We decompose the signal in spherical harmonic modes but, for illustrative purposes, we concentrate on the dominant and subdominant modes. The Cauchy data were given at the extraction radii . The relationship between the Cauchy radius and the characteristic worldtube radius is: . The characteristic extraction module was run with the following specifications: angular gridpoints = radial gridpoints = , and timestep , where . The test was run until , using order accurate angular derivatives, on stereographic patches with circular boundaries and angular dissipation (see (15) for details on how the angular dissipation is added to the evolution equation 9). The results are shown for the highest resolution. Table 1 gives the convergence rates for the worldtube variables obtained with a small extraction radius at a time corresponding to the peak of the signal (). The rates are given for the real and the imaginary part. All quantities are very close to second order convergent, including , which is the term which determines the waveform.
We are reporting only firstorder convergence rates at future null infinity for the Bondi News and the Weyl complonent , but the error is relatively small ( during the late inspiral). From an extraction point of view, these errors are smaller than the error in the Cauchy code data and are of little concern. The data are not convergent at early time when high frequencies dominate the error. For a thorough analysis of the causes for the firstorder accurate results and major improvements to the code see (18).
Figure 4 compares the imaginary and real parts of the mode of the Cauchy with the complex conjugate of the mode of the characteristic . We obtain very good amplitude match. Also, we observe improved phase agreement as the extraction radius is increased (from to ), because the phase error in is reduced with increased extraction radius.
Figure 5 compares the imaginary and real parts of the mode of the Cauchy with the complex conjugate of the mode of the characteristic . Here we see two effects. First the improved phase agreement as , but also an attenuation of the amplitude due to dissipation of higherorder modes. Also, the noise is apparent for the mode.
Figure 6 compares the amplitudes and the phases between the absolute value of the mode of the Cauchy extracted at , and the absolute value of the mode of the characteristic for the same extraction radius, at the highest resolution (). The difference in amplitude is relatively small, maximum of the Cauchy amplitude in the wave zone.
Figure 7 compares the real part of the modes of the characteristic extracted at three different extraction radii: , and . The waveform extracted at has the biggest amplitude, and a very small attenuation of the signal with the radius is observed.
Vi Conclusion
We have presented here a method for interfacing outer boundary data from a Cauchy evolution with inner boundary data for a characteristic evolution so that the waveform can be accurately extracted at infinity. We have demonstrated how the PITT null code can be interfaced with the LazEv code, which is a finitedifference BSSN code, to produced calibrated waveforms from a binary blackhole inspiral. The extraction interface has been implemented as a thorn in the Einstein Computational Toolkit ein () In this paper we are reporting only preliminary results (see (18) improvements). Although we are aware of the deficiencies in the characteristic waveform extraction tool presented here, there is pressing interest from several numerical relativity groups to apply the tool to extract waveforms from binary blackhole inspirals.
References
 B.S. Sathyaprakash, and B. F. Shultz, Living Rev. Rel., 12 2 (2009), arXiv:0903.0338
 The LIGO Scientific Collaboration and the Virgo Collaboration, Nature, Vol. 460 990 (2009)
 F. Pretorius, Class. Quant. Grav, 22 425 (2005), grqc/0407110
 L. Lehner, and O. M. Moreschi, Phys. Rev. D, 76 124040 (2007), arXiv:0706.1319
 E. T. Newman, and R. Penrose, J. Math. Phys., 3 566 (1962)
 V. Moncrief, Ann. Phys, 88 323342 (1974)
 L. Lindblom, B. J. Owen, D. A. Brown,
 N.T. Bishop, R. Gómez, L. Lehner, and J. Winicour, Phys. Rev. D, 54 6153 (1996).
 Winicour, J. (2005), Living Rev. Relativity 8, 10.
 L.A. Tamburino and J. Winicour, Phys. Rev. 150 1039, 1966.
 R.A. Isaacson, J.S. Welling and J. Winicour, J. Math. Phys. 24 1824 (1983).
 N. T. Bishop, R. Gómez, L. Lehner, M. Maharaj and J. Winicour, Phys. Rev. D,56 6298 (1997).
 M. Babiuc, B. Szilágyi, I. Hawke, and Y. Zlochower, Class. Quantum Grav. 22 5089 (2005)
 N. T. Bishop, R. Gómez, S. Husa, L. Lehner, J. Winicour Phys.Rev. D 68, 084015 (2003)
 M.C. Babiuc, N.T. Bishop, B. Szilágyi and J. Winicour, Phys. Rev. D, 79 084011 (2009).
 C. Reisswig, N. T. Bishop, D. Pollney abd B. Szilágyi, Phys. Rev. Lett., 95, 221101 (2009).
 C. Reisswig, N. T. Bishop, D. Pollney abd B. Szilágyi, Class. Quantum Grav. 27, 075014 (2010).
 M.C. Babiuc, B. Szilágyi, J. Winicour and Y. Zlochower [arXiv:1011.4223]
 H. Bondi, M.J.G. van der Burg and A.W.K. Metzner, Proc. R. Soc. A 269 21 (1962).
 R.K. Sachs, Proc. R. Soc. A 270 103 (1962).
 R. Penrose, Phys. Rev. Letters, 10 66 (1963).
 N.T. Bishop, R. Gómez, L. Lehner, B. Szilágyi, J. Winicour and R. A. Isaacson, Black Holes, Gravitational Radiation and the Universe, Kluwer Academic Publishers, Dordrecht, 1998
 R. Gómez, L. Lehner, P. Papadopoulos and J. Winicour, Class. Quantum Grav. 14 977, 1997.
 R. Gómez, Phys. Rev. D 64, 1–8 (2001).
 J. Winicour, J. Math. Phys. 24 1193 (1983).
 J. Winicour, J. Math. Phys. 25 2506 (1984).
 M. Campanelli, C. O. Lousto, P. Marronetti and Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006) [arXiv:grqc/0511048].
 Y. Zlochower, J. G. Baker, M. Campanelli and C. O. Lousto, Phys. Rev. D 72, 024021 (2005) [arXiv:grqc/0505055].
 Cactus Computational Toolkit home page: http://www.cactuscode.org.
 E. Schnetter, S. H. Hawley and I. Hawke, Class. Quant. Grav. 21, 1465 (2004) [arXiv:grqc/0310042].
 M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995).
 T. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 024007 (1999).
 J. G. Baker, J. Centrella, D.I. Choi, M. Koppitz, and J. van Meter, Phys.Rev.Lett. 96, 111102 (2006).
 M. Campanelli, C. O. Lousto and Y. Zlochower, Phys. Rev. D 74, 041501 (2006) [arXiv:grqc/0604012].
 Einstein Computational Toolkit home page: http://einsteintoolkit.org.