3.6 Characteristic Treatment of Binary 3 Characteristic Evolution Codes3.4 The Bondi Mass

3.5 3D Characteristic Evolution 

There has been rapid progress in 3D characteristic evolution. There are now two independent 3D codes, one developed at Canberra and the other at Pittsburgh (the PITT code), which have the capability to study gravitational waves in single black hole spacetimes, at a level still not mastered by Cauchy codes. Three years ago the Pittsburgh group established robust stability and second order accuracy of a fully nonlinear 3D code which calculates waveforms at null infinity [31Jump To The Next Citation Point In The Article, 30Jump To The Next Citation Point In The Article] and also tracks a dynamical black hole and excises its internal singularity from the computational grid [73Jump To The Next Citation Point In The Article, 70Jump To The Next Citation Point In The Article]. The Canberra group has implemented an independent nonlinear 3D code which accurately evolves the exterior region of a Schwarzschild black hole. Both codes pose data on an initial null hypersurface and on a worldtube boundary, and evolve the exterior spacetime out to a compactified version of null infinity, where the waveform is computed. However, there are essential differences in the underlying geometrical formalisms and numerical techniques used in the two codes and their success in evolving generic black hole spacetimes.

3.5.1 Geometrical formalism

The PITT code uses a standard Bondi-Sachs null coordinate system whereas the Canberra code employs a null quasi-spherical (NQS) gauge (not to be confused with the quasi-spherical approximation in which quadratically aspherical terms are ignored [31Jump To The Next Citation Point In The Article]). The NQS gauge takes advantage of the possibility of mapping the angular part of the Bondi metric conformally onto a unit sphere metric, so that tex2html_wrap_inline2075 . The required transformation tex2html_wrap_inline2077 is in general dependent upon u and r so that the quasi-spherical angular coordinates tex2html_wrap_inline2083 are not constant along the outgoing null rays, unlike the Bondi-Sachs angular coordinates. Instead the coordinates tex2html_wrap_inline2083 display the analogue of a shift on the null hypersurfaces tex2html_wrap_inline1695 . The radiation content of the metric is contained in a shear vector describing this shift. This results in the description of the radiation in terms of a spin-weight 1 field, rather than the spin-weight 2 field associated with tex2html_wrap_inline2061 in the Bondi-Sachs formalism. In both the Bondi-Sachs and NQS gauges, the independent gravitational data on a null hypersurface is the conformal part of its degenerate 3-metric. The Bondi-Sachs null data consists of tex2html_wrap_inline2061, which determines the intrinsic conformal metric of the null hypersurface. In the quasi-spherical case, tex2html_wrap_inline2093 and the shear vector comprises the only non-trivial part of the conformal 3-metric. Both the Bondi-Sachs and NQS gauges can be arranged to coincide in the special case of shear-free Robinson-Trautman metrics [49, 16Jump To The Next Citation Point In The Article].

The formulation of Einstein's equations in the NQS gauge is presented in Ref. [15] and the associated gauge freedom arising from (u, r) dependent rotation and boosts of the unit sphere is discussed in Ref. [16]. As in the PITT code, the main equations involve integrating a hierarchy of hypersurface equations along the radial null geodesics extending from the inner boundary to null infinity. In the NQS gauge the source terms for these radial ODE's are rather simple when the unknowns are chosen to be the connection coefficients. However, as a price to pay for this simplicity, after the radial integrations are performed on each null hypersurface a first order elliptic equation must be solved on each tex2html_wrap_inline2097 cross-section to reconstruct the underlying metric.

3.5.2 Numerical Methodology

The PITT code is an explicit second order finite difference evolution algorithm based upon retarded time steps on a uniform three-dimensional null coordinate grid. The straightforward numerical approach and the second order convergence of the finite difference equations has facilitated code development. The Canberra code uses an assortment of novel and elegant numerical methods. Most of these involve smoothing or filtering and have obvious advantage for removing short wavelength noise but would be unsuitable for modeling shocks.

Both codes require the ability to handle tensor fields and their derivatives on the sphere. Spherical coordinates and spherical harmonics are natural analytic tools for the description of radiation but their implementation in computational work requires dealing with the impossibility of smoothly covering the sphere with a single coordinate grid. Polar coordinate singularities in axisymmetric systems can be regularized by standard tricks. In the absence of symmetry, these techniques do not generalize and would be especially prohibitive to develop for tensor fields.

A crucial ingredient of the PITT code is the eth -module [72] which incorporates a computational version of the Newman-Penrose eth-formalism [107]. The eth-module covers the sphere with two overlapping stereographic coordinate grids (North and South). It provides everywhere regular, second order accurate, finite difference expressions for tensor fields on the sphere and their covariant derivatives. The eth-calculus simplifies the underlying equations, avoids spurious coordinate singularities and allows accurate differentiation of tensor fields on the sphere in a computationally efficient and clean way. Its main weakness is the numerical noise introduced by interpolations (4th order accurate) between the North and South patches. For parabolic or elliptic equations on the sphere, the finite difference approach of the eth-calculus would be less efficient than a spectral approach, but no parabolic or elliptic equations appear in the Bondi-Sachs evolution scheme.

The Canberra code handles fields on the sphere by means of a 3-fold representation [20]: (i) discretized functions on a spherical grid uniformly spaced in standard tex2html_wrap_inline2099 coordinates, (ii) fast-Fourier transforms with respect to tex2html_wrap_inline2099 (based upon a smooth map of the torus onto the sphere), and (iii) spectral decomposition of scalar, vector and tensor fields in terms of spin-weighted spherical harmonics. The grid values are used in carrying out nonlinear algebraic operations; the Fourier representation is used to calculate tex2html_wrap_inline2099 -derivatives; and the spherical harmonic representation is used to solve global problems, such as the solution of the first order elliptic equation for the reconstruction of the metric, whose unique solution requires pinning down the tex2html_wrap_inline2105 gauge freedom. The sizes of the grid and of the Fourier and spherical harmonic representations are coordinated. In practice, the spherical harmonic expansion is carried out to 15th order in tex2html_wrap_inline2107 but the resulting coefficients must then be projected into the tex2html_wrap_inline2109 subspace in order to avoid inconsistencies between the spherical harmonic and the grid and Fourier representations.

The Canberra code solves the null hypersurface equations by combining an 8th order Runge-Kutta integration with a convolution spline to interpolate field values. The radial grid points are dynamically positioned to approximate ingoing null geodesics, a technique originally due to Goldwirth and Piran [67] to avoid inaccuracies near a horizon resulting from a uniform r -grid due to the coordinate singularity in the case of a stationary horizon. The time evolution uses the method of lines with a 4th order Runge-Kutta integrator, which introduces further high frequency filtering.

3.5.3 Stability

PITT Code : Analytic stability analysis of the finite difference equations has been crucial in the development of a stable evolution algorithm, subject to the standard Courant-Friedrichs-Lewy (CFL) condition for an explicit code. Linear stability analysis on Minkowski and Schwarzschild backgrounds showed that certain field variables must be represented on the half-grid [74, 31Jump To The Next Citation Point In The Article]. Nonlinear stability analysis was essential in revealing and curing a mode coupling instability that was not present in the original axisymmetric version of the code [30Jump To The Next Citation Point In The Article, 98]. Stability persists even in the regime that the u -direction, along which the grid flows, becomes spacelike, such as outside the velocity of light cone in a rotating coordinate system. Severe tests were used to verify stability. In the linear regime, robust stability was established by imposing random initial data on the initial characteristic hypersurface and random constraint violating boundary data on an inner worldtube. The code ran stably for 10,000 grid crossing times under these conditions [31Jump To The Next Citation Point In The Article, 138Jump To The Next Citation Point In The Article]. Except for a linearized version of an Arnowitt-Deser-Misner code [138Jump To The Next Citation Point In The Article] (see also Sec.  4), the PITT code is the only 3D general relativistic code which has passed this test of robust stability. The use of random data is only possible in sufficiently weak cases where effective energy terms quadratic in the field gradients are not dominant. Stability in the highly nonlinear regime was tested in two ways. Runs for a time of 60,000 M were carried out for a moving, distorted Schwarzschild black hole (of mass M), with the marginally trapped surface at the inner boundary tracked and its interior excised from the computational grid [70Jump To The Next Citation Point In The Article, 143]. This is the longest black hole simulation carried out to date. Furthermore, the scattering of a gravitational wave off a Schwarzschild black hole was successfully carried out in the extreme nonlinear regime where the backscattered Bondi news was as large as N =400 (in dimensionless geometric units) [30Jump To The Next Citation Point In The Article], showing that the code can cope with the enormous power output tex2html_wrap_inline2121 in conventional units. This exceeds the power that would be produced if, in 1 second, the entire galaxy were converted to gravitational radiation.

Canberra code : Direct stability analysis of the underlying finite difference equations is impractical because of the extensive mix of spectral techniques, higher order methods and splines. Similarly, there is no clear-cut CFL limit on the code, although stability tests show that there is a limit on the time step. The damping of high frequency modes due to the implicit filtering would be expected to suppress numerical instability, but the stability of the Canberra code is nevertheless subject to two qualifications [18, 19, 21]: (i) At late times (less than 100 M), the evolution terminates as it approaches an event horizon, apparently because of a breakdown of the NQS gauge condition, although an analysis of how and why this should occur has not yet been given. (ii) Numerical instabilities arise from dynamic inner boundary conditions and restrict the inner boundary to a fixed Schwarzschild horizon. Tests in the extreme nonlinear regime were not reported.

3.5.4 Accuracy

PITT Code : Second order accuracy was verified in an extensive number of testbeds [31Jump To The Next Citation Point In The Article, 30Jump To The Next Citation Point In The Article, 70Jump To The Next Citation Point In The Article], including new exact solutions specifically constructed in null coordinates for the purpose of convergence tests:

Canberra code : The complexity of the algorithm and NQS gauge makes it problematic to establish accuracy by direct means. Exact solutions do not provide an effective convergence check because the Schwarzschild solution is trivial in the NQS gauge and other known solutions in this gauge require dynamic inner boundary conditions which destabilize the present version of the code. Convergence to linearized solutions is a possible check but has not yet been performed. Instead indirect tests by means of geometric consistency and partial convergence tests are used to calibrate accuracy. The consistency tests are based on the constraint equations, which are not enforced during null evolution except at the inner boundary. The balance between mass loss and radiation flux through tex2html_wrap_inline1655 is a global version of these constraints. No appreciable growth of the constraints is noticeable until within 5 M of the final breakdown of the code. In weak field tests where angular resolution does not dominate the error, partial convergence tests based upon varying the radial grid size show the 8th order convergence in the shear expected from the Runge-Kutta integration and splines. When the radial source of error is small, reduced error with smaller time step can also be discerned.

In practical runs, the major source of inaccuracy is the spherical harmonic resolution, currently restricted to tex2html_wrap_inline2129 by hardware limitations. Truncation of the spherical harmonic expansion has the effect of modifying the equations to a system for which the constraints are no longer satisfied. The relative error in the constraints is tex2html_wrap_inline2131 for strong field simulations [17Jump To The Next Citation Point In The Article].

3.5.5 Nonlinear Scattering Off a Schwarzschild Black Hole

A natural physical application of a characteristic evolution code is the nonlinear version of the classic problem of scattering off a Schwarzschild black hole, first solved perturbatively by Price [116]. Here the inner worldtube for the characteristic initial value problem consists of the ingoing r =2 m surface (the past horizon), where Schwarzschild data are prescribed. The nonlinear problem of a gravitational wave scattering off a Schwarzschild black hole is then posed in terms of data on an outgoing null cone consisting of an incoming pulse with compact support. Part of the energy of this pulse falls into the black hole and part is backscattered out to tex2html_wrap_inline1655 . This problem has been investigated using both the PITT and Canberra codes.

The Pittsburgh group studied the backscattered waveform (described by the Bondi news function) as a function of incoming pulse amplitude. The computational eth module smoothly handled the complicated time dependent transformation between the non-inertial computational frame at tex2html_wrap_inline1655 and the inertial (Bondi) frame necessary to obtain the standard ``plus'' and ``cross'' polarization modes. In the perturbative regime, the news corresponds to the backscattering of the incoming pulse off the effective Schwarzschild potential. When the energy of the pulse is no larger than the central Schwarzschild mass, the backscattered waveform still depends roughly linearly on the amplitude of the incoming pulse. However, for very high amplitudes the waveform behaves quite differently. Its amplitude is greater than that predicted by linear scaling and its shape drastically changes and exhibits extra oscillations. In this very high amplitude case, the mass of the system is completely dominated by the incoming pulse, which essentially backscatters off itself in a nonlinear way.

The Canberra code was used to study the change in Bondi mass due to radiation [17]. The Hawking mass tex2html_wrap_inline2139 is first calculated as a function of radius and retarded time, with the Bondi mass tex2html_wrap_inline2141 then obtained by taking the limit tex2html_wrap_inline2051 . The limit has good numerical behavior. For a strong initial pulse with tex2html_wrap_inline2145 angular dependence, in a run from u =0 to u =70 (in units where the interior Schwarzschild mass is 1), the Bondi mass drops from 1.8 to 1.00002, showing that almost half of the initial energy of the system was backscattered and that a surprisingly negligible amount falls into the black hole. The Bondi mass decreases monotonically in time, as theoretically necessary, but its rate of change exhibits an interesting pulsing behavior whose time scale cannot be obviously explained in terms of quasinormal oscillations. The Bondi mass loss formula is confirmed with relative error of less than tex2html_wrap_inline2131 . This is impressive accuracy considering the potential sources of numerical error introduced by taking the limit of the Hawking mass. The Canberra group also studied the effect of initial pulse amplitude on the waveform of the backscattered radiation but did not extend their study to the very high amplitude regime in which qualitatively interesting nonlinear effects occur.

3.5.6 Black Hole in a Box

The PITT code has also been implemented to evolve along a family of ingoing null cones, with data given on a worldtube at their outer boundary and on an initial ingoing null cone. The code was used to evolve a black hole in the region interior to the worldtube by locating the marginally trapped surface (MTS) on the ingoing cones and excising its singular interior [73Jump To The Next Citation Point In The Article]. The code tracks the motion of the MTS and measures its area during the evolution. It was used to simulate a distorted ``black hole in a box'' [70Jump To The Next Citation Point In The Article]. Data at the outer worldtube was induced from a Schwarzschild or Kerr spacetime but the worldtube was allowed to move relative to the stationary trajectories; i.e. with respect to the grid the worldtube is fixed but the black hole moves inside it. The initial null data consisted of a pulse of radiation which subsequently travels outward to the worldtube where it reflects back toward the black hole. The approach of the system to equilibrium was monitored by the area of the MTS, which also equals its Hawking mass. When the worldtube is stationary (static or rotating in place), the distorted black hole inside evolves to equilibrium with the boundary. A boost or other motion of the worldtube with respect to the black hole does not affect this result. The marginally trapped surface always reaches equilibrium with the outer boundary, confirming that the motion of the boundary is ``pure gauge''.

The code essentially runs ``forever'' even when the worldtube wobbles with respect to the black hole to produce artificial periodic time dependence. An initially distorted, wobbling black hole was evolved for a time of 60,000 M, longer by orders of magnitude than stability permits for any other existing 3D black hole code. This exceptional performance opens a promising new approach to handle the inner boundary condition for Cauchy evolution of black holes by the matching methods reviewed in Sec.  4 .

3.6 Characteristic Treatment of Binary 3 Characteristic Evolution Codes3.4 The Bondi Mass

image Characteristic Evolution and Matching
Jeffrey Winicour
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de