It is my pleasure to review progress in numerical relativity based upon characteristic evolution. In the spirit of Living Reviews in Relativity, I invite my colleagues to continue to send me contributions and comments at winicour@pitt.edu.

We are now in an era in which Einstein’s equations can effectively be considered solved at the local level. Several groups, as reported here and in other Living Reviews in Relativity, have developed 3D codes which are stable and accurate in some sufficiently local setting. Since my last review, the pioneering works [195] (based upon a harmonic formulation) and [66, 19] (based upon BSSN formulations [218, 32]) have initiated dramatic progress in the ability of Cauchy codes to simulate the inspiral and merger of binary black holes, the premier problem in classical relativity. These codes evolve the spacetime inside an artificially constructed outer boundary. Global solutions of binary black holes are another matter. In particular, none of the existing codes are capable of computing the waveform of the gravitational radiation emanating from a binary inspiral with only the error introduced by the numerical approximation. The best that these Cauchy codes by themselves can provide is an extraction of the waveform based upon the analytical approximation that the outer boundary is at infinity. Just as several coordinate patches are necessary to describe a spacetime with nontrivial topology, the most effective attack on the binary black hole waveform might involve a global solution patched together from pieces of spacetime handled by a combination of different codes and techniques.

Most of the effort in numerical relativity has centered about the Cauchy {3 + 1} formalism [253], with the gravitational radiation extracted by perturbative methods based upon introducing an artificial Schwarzschild background in the exterior region [1, 4, 2, 3, 208, 204, 177]. These wave extraction methods have not been thoroughly tested in a nonlinear 3D setting. A different approach which is specifically tailored to study radiation is based upon the characteristic initial value problem. In the 1960s, Bondi [53, 54] and Penrose [187] pioneered the use of null hypersurfaces to describe gravitational waves. Their new approach has flourished in general relativity. It led to the first unambiguous description of gravitational radiation in a fully nonlinear context. By formulating asymptotic flatness in terms of characteristic hypersurfaces extending to infinity, they were able to reconstruct, in a nonlinear geometric setting, the basic properties of gravitational waves which had been developed in linearized theory on a Minkowski background. The major new nonlinear features were the Bondi mass and news function, and the mass loss formula relating them. The Bondi news function is an invariantly defined complex radiation amplitude , whose real and imaginary parts correspond to the time derivatives and of the “plus” and “cross” polarization modes of the strain incident on a gravitational wave antenna. The corresponding waveforms are important both for the design of detection templates for a binary black hole inspiral and merger and for the determination of the resulting recoil velocity.

The recent success of Cauchy evolutions in simulating binary black holes emphasizes the need of applying these global techniques to accurate waveform extraction. This has stimulated several attempts to increase the accuracy of characteristic evolution. The successful Cauchy simulation of binary black holes has been based upon adaptive mesh techniques and upon fourth order finite difference approximations. The initial characteristic codes were developed with unigrid second order accuracy. One of the prime factors affecting the accuracy of any characteristic code is the introduction of a smooth coordinate system covering the sphere, which labels the null directions on the outgoing light cones. This is also an underlying problem in meteorology and oceanography. In a pioneering paper on large-scale numerical weather prediction, Phillips [189] put forward a list of desirable features for a mapping of the sphere to be useful for global forecasting. The first requirement was the freedom from singularities. This led to two distinct choices which had been developed earlier in purely geometrical studies: stereographic coordinates (two coordinate patches) and cubed-sphere coordinates (six patches). Both coordinate systems have been rediscovered in the context of numerical relativity (see Section 4.1). The cubed-sphere method has stimulated two new attempts at improved codes for characteristic evolution (see Section 4.2.3).

Another issue affecting code accuracy is the choice between a second vs first differential order reduction of the evolution system. Historically, the predominant importance of computational fluid dynamics has favored first order systems, in particular the reduction to symmetric hyperbolic form. However, in acoustics and elasticity theory, where the natural treatment is in terms of second order wave equations, an effective argument for the second order form has been made [159]. In general relativity, the question of whether first or second order formulations are more natural depends on how Einstein’s equations are reduced to a hyperbolic system by some choice of coordinates and variables. The second order form is more natural in the harmonic formulation, where the Einstein equations reduce to quasilinear wave equations. The first order form is more natural in the Friedrich-Nagy formulation, which includes the Weyl tensor among the evolution variables, and was used in the first demonstration of a well-posed initial-boundary value problem for Einstein’s equations. Recent investigations of first order formulations of the characteristic initial value problem are discussed in Section 4.2.2.

The major drawback of a stand-alone characteristic approach arises from the formation of caustics in the light rays generating the null hypersurfaces. In the most ambitious scheme proposed at the theoretical level such caustics would be treated “head-on” as part of the evolution problem [233]. This is a profoundly attractive idea. Only a few structural stable caustics can arise in numerical evolution, and their geometrical properties are well enough understood to model their singular behavior numerically [99], although a computational implementation has not yet been attempted.

In the typical setting for the characteristic initial value problem, the domain of dependence of a single smooth null hypersurface is empty. In order to obtain a nontrivial evolution problem, the null hypersurface must either be completed to a caustic-crossover region where it pinches off, or an additional inner boundary must be introduced. So far, the only caustics that have been successfully evolved numerically in general relativity are pure point caustics (the complete null cone problem). When spherical symmetry is not present, the stability conditions near the vertex of a light cone place a strong restriction on the allowed time step [124]. Nevertheless, point caustics in general relativity have been successfully handled for axisymmetric vacuum spacetimes [120]. Progress toward extending these results to realistic astrophysical sources has been made in two Ph.D. theses based upon characteristic evolution codes. Florian Siebel’s thesis work [219], at the Technische Universität München, integrated an axisymmetric characteristic gravitational-hydro code with a high resolution shock capturing code for the relativistic hydrodynamics. This coupled general relativistic code has been thoroughly tested and has yielded state-of-the-art results for the gravitational waves produced by the oscillation and collapse of a relativistic star (see Sections 7.1 and 7.2). Nevertheless, computational demands to extend these results to 3D evolution would be prohibitive using current generation supercomputers, due to the small timestep required at the vertex of the null cone (see Section 3.3). This is unfortunate because, away from the caustics, characteristic evolution offers myriad computational and geometrical advantages. Vacuum simulations of black hole spacetimes, where the inner boundary can be taken to be the white hole horizon, offer a scenario where both the timestep and caustic problems can be avoided and three-dimensional simulations are practical. In Yosef Zlochower’s thesis work [255], at the University of Pittsburgh, the gravitational waves generated from the post-merger phase of a binary black black hole were computed using a fully nonlinear three-dimensional characteristic code [256] (see Section 4.5). He showed how the characteristic code could be employed to investigate the nonlinear mode coupling in the response of a black hole to the infall of gravitational waves.

At least in the near future, fully three-dimensional computational applications of characteristic evolution are likely to be restricted to some mixed form, in which data is prescribed on a non-singular but incomplete initial null hypersurface and on a second inner boundary , which together with the initial null hypersurface determine a nontrivial domain of dependence. The hypersurface may be either (i) null, (ii) timelike or (iii) spacelike, as schematically depicted in Figure 1. The first two possibilities give rise to (i) the double null problem and (ii) the nullcone-worldtube problem. Possibility (iii) has more than one interpretation. It may be regarded as a Cauchy initial-boundary value problem where the outer boundary is null. An alternative interpretation is the Cauchy-characteristic matching (CCM) problem, in which the Cauchy and characteristic evolutions are matched transparently across a worldtube W, as indicated in Figure 1.

In CCM, it is possible to choose the matching interface between the Cauchy and characteristic regions to be a null hypersurface, but it is more practical to match across a timelike worldtube. CCM combines the advantages of characteristic evolution in treating the outer radiation zone in spherical coordinates which are naturally adapted to the topology of the worldtube with the advantages of Cauchy evolution in treating the inner region in Cartesian coordinates, where spherical coordinates would break down.

In this review, we trace the development of characteristic algorithms from model 1D problems to a 2D axisymmetric code which computes the gravitational radiation from the oscillation and gravitational collapse of a relativistic star and to a 3D code designed to calculate the waveform emitted in the merger to ringdown phase of a binary black hole. And we trace the development of CCM from early feasibility studies to successful implementation in the linear regime and through current attempts to treat the binary black hole problem.

A major achievement has been the successful application of CCM to the linearized matching problem between a 3D characteristic code and a 3D Cauchy code based upon harmonic coordinates [236] (see Section 5.8). Here the linearized Cauchy code satisfies a well-posed initial-boundary value problem, which seems to be a critical missing ingredient in previous attempts at CCM in general relativity. Recently a well-posed initial-boundary value problem has been established for fully nonlinear harmonic evolution [162] (see Section 5.3), which should facilitate the extension of CCM to the nonlinear case.

In previous reviews, I tried to include material on the treatment of boundaries in the computational mathematics and fluid dynamics literature because of its relevance to the CCM problem. The fertile growth of this subject warrants a separate Living Review in Relativity on boundary conditions, which is presently under consideration. In anticipation of this, I will not attempt to keep this subject up to date except for material of direct relevance to CCM. See [213] for an independent review of boundary conditions that have been used in numerical relativity.

The problem of computing the evolution of a neutron star, in close orbit about a black hole is of clear importance to the new gravitational wave detectors. The interaction with the black hole could be strong enough to produce a drastic change in the emitted waves, say by tidally disrupting the star, so that a perturbative calculation would be inadequate. The understanding of such nonlinear phenomena requires well behaved numerical simulations of hydrodynamic systems satisfying Einstein’s equations. Several numerical relativity codes for treating the problem of a neutron star near a black hole have been developed, as described in the Living Review in Relativity on “Numerical Hydrodynamics in General Relativity” by Font [91]. Although most of these efforts concentrate on Cauchy evolution, the characteristic approach has shown remarkable robustness in dealing with a single black hole or relativistic star. In this vein, state-of-the-art axisymmetric studies of the oscillation and gravitational collapse of relativistic stars have been achieved (see Section 7.2) and progress has been made in the 3D simulation of a body in close orbit about a Schwarzschild black hole (see Sections 4.6 and 7.3).

http://www.livingreviews.org/lrr-2009-3 |
This work is licensed under a Creative Commons License. Problems/comments to |