#
Phys. Rev. D 64, 024007 (2001)

Gravitational waveforms with controlled accuracy

###### Abstract

A partially first-order form of the characteristic formulation is introduced to control the accuracy in the computation of gravitational waveforms produced by highly distorted single black hole spacetimes. Our approach is to reduce the system of equations to first-order differential form on the angular derivatives, while retaining the proven radial and time integration schemes of the standard characteristic formulation. This results in significantly improved accuracy over the standard mixed-order approach in the extremely nonlinear post-merger regime of binary black hole collisions.

###### pacs:

04.20.Ex, 04.25.Dm, 04.25.Nx, 04.70.Bw## I Introduction

One of the most impressive achievements of the characteristic approach is its demonstrated ability to carry out long-term stable evolutions [1] of generic single black hole spacetimes. Work is under way to extend these numerical simulations to the post-merger regime of binary black hole coalescence, by first computing the gravitational radiation emitted during a white hole fission [2], and then extending these results to compute the gravitational waveforms emitted during the post-merger phase of a binary black hole collision [3].

The present work is an outgrowth of one of these projects, having been motivated by numerical difficulties that we have encountered with the present implementation of the characteristic code [4]. In the course of our simulations of a white-hole fission [2], we have noticed an angular oscillation mode in some of the metric variables of the characteristic formulation, similar in appearance to a red-black [5] decoupling. Numerical experiments reveal that although the area that it affects can be somewhat reduced, its amplitude does not converge to zero with increasing grid resolution for practical grid sizes in the rapidly changing environment near the time of merger of a binary black hole. We had expected to see convergence in this regime, given the measured second order convergence of the code on less rapidly changing solutions, even at fairly low resolutions [4, 6].

This effect is clearly of numerical nature, and would not have been observed in earlier work, as it is only triggered when, in addition to extreme nonlinearities in the metric functions (such as those found in the post-merger phase of the head-on collision of two black holes, or in the time-reversed picture, during the fission of a single black-hole), a fairly large degree of asymmetry and of time dependence in the metric functions is present. Moreover, it appears that even when some of these conditions are satisfied, the dissipation built into the algorithm used for the evolution the conformal metric of the spheres [7] can, in certain circumstances, be sufficient to suppress the spurious oscillations, provided the boundary data is specified analytically [8].

In the systems that we have previously considered, such as Robinson-Trautman spacetimes, a perturbed Kerr black hole [8] and the collapse of matter into a black hole [6], our choice of the gauge conditions at the boundary effectively precluded fulfilling the second condition above. In all cases, the boundary data was given analytically, so no extraneous effects could be observed.

This oscillatory mode is not an unstable mode of the finite difference approximation in the usual sense of exponential growth of the error. This is confirmed by tests where random initial and/or boundary data [9] is prescribed, and the code run for many crossing times. But it is worth pointing out that while those tests are useful to spot exponentially growing modes, they would not detect a non-exponentially growing oscillatory mode such as the one we address here.

We observe, in fact, an inaccuracy in the computation of the metric variable , which does not converge to zero, and gets triggered in situations where there are large angular variations and rapid time changes in the metric variables, specially in , Eq. (1). Because is related to the time slicing of the boundary (and in the black hole case, this boundary is the event horizon ), will indeed change substantially near the time of coalescence of a binary black hole horizon (or, in the time-reversed description actually used by our numerical code, near the time of fission of a white hole).

The PITT null code runs stably in the usual sense (without exponentially growing modes) in the configurations that we have explored. The effect of the angular noise is nevertheless important for our present purposes, since the angular error it introduces in the metric functions is most pronounced near null infinity, and this in turn spoils the accuracy of the waveforms computed.

We have found that the difficulty lies with the characteristic evolution equations themselves, and it is not a boundary-induced effect. Although the boundary (horizon) data is smooth and, by construction, satisfies the consistency conditions [10], the evolution and hypersurface equations themselves introduce the above mentioned oscillation mode, which is shown clearly in Fig. 1. We have found clear evidence that the origin of the difficulty lies with those terms in the equations that contain second angular derivatives. These type of terms are present both in the evolution equation for the conformal metric of the surfaces at constant luminosity distance , and in the hypersurface equation for the metric function itself. Most notably, when this terms are suppressed, so is the angular oscillation mode in . The most troublesome terms are those containing second angular derivatives of the metric variable , of the form and .

Second angular derivatives of the conformal metric , Eq. (1), also appear in the hypersurface equation for , in the form of and operators acting on quantities defined in terms of the conformal metric. However, these terms do not seem to lead to numerical difficulties as pronounced as those caused by the terms involving second angular derivatives of , although the reason behind this is not entirely clear at the moment. It is possible that in the present regime the contribution from these terms is not as important.

Our rationale for implementing the characteristic code in mixed first and second order form has been that doing so is entirely straightforward, and leads to an accurate and stable discretization with the least number of variables. The argument of using the least variables is compelling; however one must consider whether relaxing this requirement might not lead to better overall numerical behavior. We will show that this is indeed the case.

Rigorous stability arguments can be given for the linearized system [11] in the axisymmetric case, with the analysis extending unmodified to the full three-dimensional problem as shown in [9]. Later work [7] has shown that special care must be taken in the nonlinear case, as some unstable modes are present that are not revealed by the Von-Neuman analysis of the linear problem. One must keep in mind however that a stable discretization of a model problem as in [7] may not address other issues such as those raised here. The numerical analysis for mixed order systems is not yet at the level of sophistication achieved for first differential order systems. Long-term stable evolutions [1] have been achieved, and this is a consequence of the remarkably good choice of variables and coordinates inherent to Bondi’s system of equations, as well as of the care taken in implementing them numerically.

We have attempted to solve the problem at hand by the expeditious method of introducing further dissipation in the equations, in the spirit of [7]. Regrettably, these numerical remedies do not seem to have the desired effect, apparently because they do not go to the root of the problem.

We consider here a different solution, which is analytic in nature, but motivated by a numerical application, specifically the need to attain controlled accuracy in the computation of gravitational waveforms out of highly distorted single black hole systems. Our approach is to remove the offending second angular derivatives from our system of equations altogether, by introducing additional variables, in terms of which the enlarged system can be written in first differential order form in the angular coordinates. We thus implement a partial reduction to a first differential order system, leaving the radial part of the evolution equations in second order form, and allowing for mixed radial and angular derivatives, both in the original equations and in the auxiliary equations introduced. One goal is to preserve, where appropriate, the radial and time integration schemes developed for the PITT null code, whose long-term stability (in the usual sense of not having exponentially growing modes) has been amply demonstrated [1]. At the same time, we perform the necessary modification to the system of equations to make the PITT null code useful for the problem at hand.

The plan of the present paper is as follows: in Sec. II we briefly review the characteristic formulation in its standard second-order form. This material has appeared in Refs. [4, 9], and we include it here in a condensed form, as needed for the derivation of the partially reduced system. Our new results are contained in Sec. III and onward, where we present a partial reduction of the characteristic formulation, with all second order angular derivatives eliminated with the aid of a (minimal) set of additional variables. In Sec. IV, we review the discretization strategy for the system of equations introduced in Sec. III. Finally, we present an application of the partially reduced system in Sec. V to the computation of the spacetime exterior to a white hole horizon, comparing our results with those obtained with the standard null cone formulation.

## Ii Null cone formulation in second order form

The material in this section has appeared in Refs. [4, 9] which should be consulted for a detailed derivation. Here we present the main equations (including all the nonlinear terms), in a slightly more compact form than that of Ref. [4], for ease of reference and to provide a starting point for the approach developed in Sec. III.

We use coordinates based upon a family of outgoing null hypersurfaces, letting label these hypersurfaces, , label the null rays, and be a surface area coordinate, such that in the coordinates, the metric takes the Bondi-Sachs form [12, 13]

(1) | |||||

where is related to the Bondi-Sachs variable by ; and where and , with a unit sphere metric, given in terms of a complex dyad satisfying , , , with and .

We represent tensors on the sphere by spin-weighted variables [14]. The conformal metric , which satisfies the condition , is represented by the complex function , and by the real function , where . The metric functions are similarly encoded in the complex function [4, 9]. Angular derivatives are expressed in terms of and operators, for details, see Ref. [14].

The equations for the null cone formulation, in second order form, are obtained directly from the appropriate components of the Ricci tensor [4, 9]. These equations form a hierarchy, splitting into hypersurface equations, which involve only derivatives on the null cone, i.e. provides an equation for in terms of , while gives in terms of and , and the trace yields in terms of , and , respectively.

(2) | |||||

(3) | |||||

(4) | |||||

(5) | |||||

and the evolution equation for , obtained from ,

(8) |

where

(9) |

In Refs. [4, 9] we have made a partial reduction to a first order system by introducing the auxiliary variable , which plays the role of . We sometimes refer to this system as first order (in time) because only the first (retarded) time derivative of appears. See Refs. [4, 9] for details of the discretization of the above system of equations.

## Iii Null cone formulation in partially reduced form

A reduction to first order form of the system of equations in Ref. [9], in the quasi-spherical approximation, was considered by Frittelli and Lehner [15]. In their work, the main variables are , , and , which measure deviation from Schwarzschild. The auxiliary variables they introduce to reduce the system to first order form are the following radial and angular derivatives of the metric functions, , , and .

Here we do not work with tensor components directly, but instead, following Refs. [9, 4], we represent tensors in terms of spin-weighted variables. For example, instead of writing the equations in terms of , we work with , and express its angular derivatives in terms of and operators. Second angular derivatives of can in turn be expressed in terms of and .

Eliminating all second angular derivatives in the hypersurface and evolution equations necessitates the introduction of exactly three complex variables. Since by symmetry there are only three independent components of the conformal metric , there are then only six possible angular derivatives , which are encoded in the three complex quantities

(10) | |||||

(11) | |||||

(12) |

with spin weight , and , respectively. The first order system for the quasi-spherical case considered in [15] introduces six new variables just for the angular derivatives of the conformal metric because no use is made there of the determinant condition. The conformal metric satisfies the condition , i.e. it contains only two independent pieces of information, which are given in our notation by the complex function . The remaining dyad component, given by the real function , is completely determined by the relation , which is a consequence of the determinant condition.

Therefore only two of the three complex quantities in Eqs. (11),(12) can be given independently, with the remaining one fixed by the determinant condition. Our choice is to introduce the additional variables and . By inspection of Eqs. (2)–(8), we note that to put the system in first order form in the angular derivatives, we need only introduce one more auxiliary variable, .

The new variables are initialized at the boundary with , and respectively, and the consistency conditions

(13) | |||||

(14) | |||||

(15) |

propagate them to the interior, being the additional hypersurface equations for the new variables.

The complete system of equations consists of Eqs. (7)–(12) of Ref. [4], together with Eq. (15), i.e., the evolution equation for

(16) | |||||

and the hierarchy of hypersurface equations

(17) | |||||

(18) | |||||

(19) | |||||

(20) | |||||

(21) | |||||

(22) | |||||

(23) | |||||

where

(24) |

In Eqs. (16)–(24) we have expressed the second angular derivatives in terms of operators acting on , and , replacing also and operators acting upon , and with the corresponding auxiliary variables whenever possible.

As in [9], we have regularized the equation by setting , and as in Appendix A of [4], we split the terms which vanish in the quasi-spherical approximation [9] into two parts, one which contains only hypersurface derivatives,

(25) | |||||

and another, where we have explicitly isolated the only non-linear term where (retarded) time derivatives of appear,

(26) |

The procedure just described eliminates all second angular derivatives from the hypersurface and evolution equations. In computing the Bondi News, we must also evaluate the conformal factor , which relates the metric on the coordinates to the metric on a Bondi frame at [4]. (The relevant expressions are given in Appendix B of Ref. [4] and we will not repeat them here.) Inspection of Eqns. (B1)-(B6) of [4] reveals that the second angular derivatives of the conformal factor enter in the calculation as well. To remove these derivatives, we introduce the auxiliary variable . Since is defined solely on , we use the consistency relation to propagate along the generators of , initializing it with , where .

## Iv Numerical integration

We use the compactified coordinate , such that at future null infinity, . Thus, the relevant radial derivatives can be written as , and , all of which are regular at . As in [9], we carry out the radial integration with the right-hand side of Eqs. (21) evaluated at midpoints of the radial grid,

(27) |

and the same treatment is applied to Eq. (23). Note that the hypersurface and evolution equations are manifestly regular at when expressed in the variable , except for the terms in Eq. (21) and in Eq. (23) respectively, which have an extra factor of . This factor, which would make them singular at , is effectively canceled by the corresponding factor of which appears in the right-hand-side of Eq. (27), hence the equations are regular for . Note also the limiting form of Eq. (21) at , , which is useful in the discretization of . For details on the numerical implementation of the eth operators, see Ref. [14]. The radial integration algorithm is explained in detail in Refs. [11, 9, 4, 7] and we will not cover it here. A departure from [4] is that we use a 3-step iterative Crank-Nicholson scheme [16] to ensure stability of the time evolution [17]. For the time integration of the conformal factor we use a combination of second-order Runge-Kutta and mid-point rule [5] integration schemes.

## V Waveforms from a fissioning white hole

We have implemented the system of equations presented in Sec. III, and used it in the calculation of the waveforms emitted by a fissioning white hole. The details of that calculation are beyond the scope of this work and will appear elsewhere [2, 3]. The spacetime exterior to a white hole horizon is computed, by solving a double-null initial-boundary value problem. Boundary data is given on a white hole horizon, and initial data is specified on an outgoing null surface emanating from the horizon.

A procedure for the construction of the complete boundary data, i.e. the intrinsic and extrinsic geometry of the white hole horizon, and a description of how to obtain the boundary data necessary to compute the exterior space-time via characteristic evolution can be found in Ref. [10]. Some regularization of the equations is needed to ensure we deal only with fields which are regular in the entire exterior spacetime [2].

By reducing the system of equations presented in Refs. [10, 2] into the form presented in Sec. III, with the addition of the auxiliary variables , and , we are able to remove the angular oscillation clearly present in the metric function when this is computed with the standard approach.

We take as our benchmark case the most nonlinear (highest eccentricity) case considered in Ref. [10], with an eccentricity parameter of . This corresponds to a case where the white hole horizon pinches off before the expansion of the outward null rays goes to zero, a necessary condition for a Bondi evolution forward in time throughout the pre-fission period.

Figures 1 and 2 show the metric function at on the north stereographic patch for both approaches. There is a clearly visible high frequency angular mode present in the edge of the patch in the calculation performed with the standard approach, as seen in Fig. 1. This oscillatory mode is clearly absent from Fig. 2.

Using both approaches (standard and reduced), we have computed the waveforms, which correspond to the radiation emitted by a highly distorted white hole. In the time-reversed scenario, this would correspond to the radiation incident from , and the white-hole horizon, to the post-merger phase of a black hole horizon formed during a binary black hole collision. How to use this information to control the amount of incoming radiation and to gain insight into the radiation emitted during the post-merger phase of a black hole collision will be discussed elsewhere [18, 2, 3].

Our motivation in introducing the partially reduced system was to control the numerical difficulties encountered in the computation of the radiation, as encoded in the Bondi news. The quality of the numerical waveforms is greatly improved in the partially reduced system. This is readily apparent comparing Figs. 3 and 4, which display the real part of the Bondi News, at a fixed time, for the standard and the reduced form of the equations. (In the plots, the region of the patch below the equator has been masked out.) While the evolution of the exterior spacetime is stable with the standard approach, the numerical noise swamps the signal. The noise is more pronounced near the edges of the stereographic patch, and it propagates to the interior. Figure 4 shows that, with the partially reduced system, the Bondi News at is perfectly smooth.

## Conclusions

We have introduced a system of equations which takes the traditional null cone formulation and casts it as a system of equations which is of first differential order in the angular variables.

The effectiveness of this approach is clearly demonstrated by comparing the metric functions at null infinity, and specially the waveforms at when computed via this approach versus those obtained with the traditional characteristic approach, in which second-order angular derivatives are present.

We have thus shown that the partially reduced system of equations is highly effective in suppressing numerical errors, which otherwise would adversely affect the calculations of waveforms in highly nonlinear scenarios.

We are currently using this approach to compute gravitational waveforms emitted during a white hole fission [2] and during the post-merger phase of a binary black hole collision [3], and we are studying their application to the characteristic problem with matter sources.

###### Acknowledgements.

We thank J. Winicour for a careful reading of the manuscript. We have benefited from conversations with S. Frittelli, and with N.T. Bishop and L. Lehner. This work has been supported by NSF PHY 9800731 to the University of Pittsburgh. Computer time for this project was provided by the Department of Physics and Astronomy, by the Pittsburgh Supercomputing Center and by NPACI.## References

- [1] R. Gómez et. al. Phys. Rev. Lett. 80, 3915 (1998).
- [2] R. Gómez, S. Husa, L. Lehner, and J. Winicour, “White hole fission”, (in preparation).
- [3] N. T. Bishop, R. Gómez, L. Lehner, and J. Winicour, “Post merger waveforms”, (in preparation).
- [4] N. T. Bishop, R. Gómez, L. Lehner, M. Maharaj, and J. Winicour, Phys. Rev. D 56, 6298 (1997).
- [5] B. P. Flannery, W. H. Press, S. A. Teukolski, and W. T. Vetterling, Numerical Recipes in Fortran 77, second ed. (Cambridge University Press, Cambridge, England, 1992).
- [6] N. T. Bishop, R. Gómez, L. Lehner, M. Maharaj, and J. Winicour, Phys. Rev. D 60, 024005 (1999).
- [7] L. Lehner, J. Comput. Phys. 149, 59 (1999).
- [8] R. Gómez, L. Lehner, R.L. Marsa, and J. Winicour, Phys. Rev. D 57, 4778 (1998)
- [9] N. T. Bishop, R. Gómez, L. Lehner, and J. Winicour, Phys. Rev. D 54, 6153 (1996)
- [10] R. Gómez, S. Husa, and J. Winicour, Phys. Rev. D (to be published), gr-qc/0009092.
- [11] R. Gómez, P. Papadopoulos, and J. Winicour, J. Math. Phys. 35, 4184 (1994)
- [12] H. Bondi, M.J.G. van der Burg and A.W.K. Metzner, Proc. R. Soc. London A269, 21 (1962).
- [13] R. K. Sachs, Proc. R. Soc. London A270, 103 (1962).
- [14] R. Gómez, L. Lehner, P. Papadopoulos, and J. Winicour, Class. Quantum Grav. 14, 977 (1997).
- [15] S. Frittelli and L. Lehner, Phys. Rev. D 59, 084012 (1999).
- [16] S. Frittelli and R. Gómez, J. Math. Phys. 41, 5535 (2000).
- [17] S. A. Teukolsky, Phys. Rev. D 61, 087501 (2000).
- [18] M. Campanelli, R. Gómez, S. Husa, J. Winicour, and Y. Zlochower, Phys. Rev. D 63, 124013 (2001).