3.1 {1 + 1}-dimensional codes

It is often said that the solution of the general ordinary differential equation is essentially known, in light of the success of computational algorithms and present day computing power. Perhaps this is an overstatement because investigating singular behavior is still an art. But, in this spirit, it is fair to say that the general system of hyperbolic partial differential equations in one spatial dimension seems to be a solved problem in general relativity. Codes have been successful in revealing important new phenomena underlying singularity formation in cosmology [36] and in dealing with unstable spacetimes to discover critical phenomena [127Jump To The Next Citation Point]. As described below, characteristic evolution has contributed to a rich variety of such results.

One of the earliest characteristic evolution codes, constructed by Corkill and Stewart [79229], treated spacetimes with two Killing vectors using a grid based upon double null coordinates, with the null hypersurfaces intersecting in the surfaces spanned by the Killing vectors. They simulated colliding plane waves and evolved the Khan–Penrose [157] collision of impulsive (δ-function curvature) plane waves to within a few numerical zones from the final singularity, with extremely close agreement with the analytic results. Their simulations of collisions with more general waveforms, for which exact solutions are not known, provided input to the understanding of singularity formation which was unforeseen in the analytic treatments of this problem.

Many {1 + 1}-dimensional characteristic codes have been developed for spherically symmetric systems. Here matter must be included in order to make the system non-Schwarzschild. Initially the characteristic evolution of matter was restricted to simple cases, such as massless Klein–Gordon fields, which allowed simulation of gravitational collapse and radiation effects in the simple context of spherical symmetry. Now, characteristic evolution of matter is progressing to more complicated systems. Its application to hydrodynamics has made significant contributions to general relativistic astrophysics, as reviewed in Section 7.

The synergy between analytic and computational approaches has led to dramatic results in the massless Klein–Gordon case. On the analytic side, working in a characteristic initial value formulation based upon outgoing null cones, Christodoulou made a penetrating study of the spherically symmetric problem [697071727374]. In a suitable function space, he showed the existence of an open ball about Minkowski space data whose evolution is a complete regular spacetime; he showed that an evolution with a nonzero final Bondi mass forms a black hole; he proved a version of cosmic censorship for generic data; and he established the existence of naked singularities for non-generic data. What this analytic tour-de-force did not reveal was the remarkable critical behavior in the transition to the black hole regime, which was discovered by Choptuik [6768Jump To The Next Citation Point] in simulations using Cauchy evolution. This phenomenon has now been understood in terms of the methods of renormalization group theory and intermediate asymptotics, and has spawned a new subfield in general relativity, which is covered in the Living Review in Relativity on “Critical Phenomena in Gravitational Collapse” by Gundlach [127].

The characteristic evolution algorithm for the spherically symmetric Einstein–Klein–Gordon problem provides a simple illustration of the techniques used in the general case. It centers about the evolution scheme for the scalar field, which constitutes the only dynamical field. Given the scalar field, all gravitational quantities can be determined by integration along the characteristics of the null foliation. This is a coupled problem, since the scalar wave equation involves the curved space metric. It illustrates how null algorithms lead to a hierarchy of equations which can be integrated along the characteristics to effectively decouple the hypersurface and dynamical variables.

In a Bondi coordinate system based upon outgoing null hypersurfaces u = const. and a surface area coordinate r, the metric is

( ) ( ) ds2 = − e2βdu V-du + 2 dr + r2 d𝜃2 + sin2 𝜃dϕ2 . (3 ) r
Smoothness at r = 0 allows imposition of the coordinate conditions
V (u,r) = r + 𝒪 (r3) (4 ) β(u,r) = 𝒪 (r2).
The field equations consist of the curved space wave equation □ Φ = 0 for the scalar field and two hypersurface equations for the metric functions:
β,r = 2 πr(Φ,r)2, (5 ) V = e2β. (6 ) ,r
The wave equation can be expressed in the form
( V ) e−2βg □ (2)g − -- ------= 0, (7 ) r ,r r
where g = rΦ and □ (2) is the D’Alembertian associated with the two-dimensional submanifold spanned by the ingoing and outgoing null geodesics. Initial null data for evolution consists of Φ(u0,r ) at the initial retarded time u0.

Because any two-dimensional geometry is conformally flat, the surface integral of □(2)g over a null parallelogram Σ gives exactly the same result as in a flat 2-space, and leads to an integral identity upon which a simple evolution algorithm can be based [122Jump To The Next Citation Point]. Let the vertices of the null parallelogram be labeled by N, E, S, and W corresponding, respectively, to their relative locations (North, East, South, and West) in the 2-space, as shown in Figure 2View Image. Upon integration of Equation (7View Equation), curvature introduces an integral correction to the flat space null parallelogram relation between the values of g at the vertices:

∫ ( ) 1- V- g- gN − gW − gE + gS = − 2 du dr r r . (8 ) Σ ,r
View Image

Figure 2: The null parallelogram. After computing the field at point N, the algorithm marches the computation to ℐ+ by shifting the corners by N → n, E → e, S → E, W → N.

This identity, in one form or another, lies behind all of the null evolution algorithms that have been applied to this system. The prime distinction between the different algorithms is whether they are based upon double null coordinates, or upon Bondi coordinates as in Equation (3View Equation). When a double null coordinate system is adopted, the points N, E, S, and W can be located in each computational cell at grid points, so that evaluation of the left hand side of Equation (8View Equation) requires no interpolation. As a result, in flat space, where the right hand side of Equation (8View Equation) vanishes, it is possible to formulate an exact evolution algorithm. In curved space, of course, there is a truncation error arising from the approximation of the integral, e.g., by evaluating the integrand at the center of Σ.

The identity (8View Equation) gives rise to the following explicit marching algorithm, indicated in Figure 2View Image. Let the null parallelogram lie at some fixed 𝜃 and ϕ and span adjacent retarded time levels u0 and u0 + Δu. Imagine for now that the points N, E, S, and W lie on the spatial grid, with rN − rW = rE − rS = Δr. If g has been determined on the entire initial cone u0, which contains the points E and S, and g has been determined radially outward from the origin to the point W on the next cone u0 + Δu, then Equation (8View Equation) determines g at the next radial grid point N in terms of an integral over Σ. The integrand can be approximated to second order, i.e. to 𝒪 (Δr Δu ), by evaluating it at the center of Σ. To this same accuracy, the value of g at the center equals its average between the points E and W, at which g has already been determined. Similarly, the value of (V ∕r),r at the center of Σ can be approximated to second order in terms of values of V at points where it can be determined by integrating the hypersurface equations (5View Equation, 6View Equation) radially outward from r = 0.

After carrying out this procedure to evaluate g at the point N, the procedure can then be iterated to determine g at the next radially outward grid point on the u0 + Δu level, i.e. point n in Figure 2View Image. Upon completing this radial march to null infinity, in terms of a compactified radial coordinate such as x = r∕(1 + r), the field g is then evaluated on the next null cone at u0 + 2Δu, beginning at the vertex where smoothness gives the startup condition that g(u, 0) = 0.

In the compactified Bondi formalism, the vertices N, E, S, and W of the null parallelogram Σ cannot be chosen to lie exactly on the grid because, even in Minkowski space, the velocity of light in terms of a compactified radial coordinate x is not constant. As a consequence, the fields g, β, and V at the vertices of Σ are approximated to second order accuracy by interpolating between grid points. However, cancellations arise between these four interpolations so that Equation (8View Equation) is satisfied to fourth order accuracy. The net result is that the finite difference version of Equation (8View Equation) steps g radially outward one zone with an error of fourth order in grid size, 𝒪 ((Δu )2(Δx )2). In addition, the smoothness conditions (4View Equation) can be incorporated into the startup for the numerical integrations for V and β to insure no loss of accuracy in starting up the march at r = 0. The resulting global error in g, after evolving a finite retarded time, is then 𝒪 (Δu Δx ), after compounding errors from 1∕(Δu Δx ) number of zones.

When implemented on a grid based upon the (u, r) coordinates, the stability of this algorithm is subject to a Courant–Friedrichs–Lewy (CFL) condition requiring that the physical domain of dependence be contained in the numerical domain of dependence. In the spherically symmetric case, this condition requires that the ratio of the time step to radial step be limited by (V ∕r)Δu ≤ 2Δr, where Δr = Δ [x∕(1 − x)]. This condition can be built into the code using the value 2H V ∕r = e, corresponding to the maximum of V ∕r at + ℐ. The strongest restriction on the time step then arises just before the formation of a horizon, where V∕r → ∞ at ℐ+. This infinite redshift provides a mechanism for locating the true event horizon “on the fly” and restricting the evolution to the exterior spacetime. Points near ℐ+ must be dropped in order to evolve across the horizon due to the lack of a nonsingular compactified version of future time infinity + I.

The situation is quite different in a double null coordinate system, in which the vertices of the null parallelogram can be placed exactly on grid points so that the CFL condition is automatically satisfied. A characteristic code based upon double null coordinates was developed by Goldwirth and Piran in a study of cosmic censorship [109Jump To The Next Citation Point] based upon the spherically symmetric gravitational collapse of a massless scalar field. Their early study lacked the sensitivity of adaptive mesh refinement (AMR) which later enabled Choptuik to discover the critical phenomena appearing in this problem. Subsequent work by Marsa and Choptuik [171Jump To The Next Citation Point] combined the use of the null related ingoing Eddington–Finklestein coordinates with Unruh’s strategy of singularity excision to construct a 1D code that “runs forever”. Later, Garfinkle [104Jump To The Next Citation Point] constructed an improved version of the Goldwirth–Piran double null code which was able to simulate critical phenomena without using adaptive mesh refinement. In this treatment, as the evolution proceeds on one outgoing null cone to the next, the grid points follow the ingoing null cones and must be dropped as they cross the origin at r = 0. However, after half the grid points are lost they are then “recycled” at new positions midway between the remaining grid points. This technique is crucial for resolving the critical phenomena associated with an r → 0 size horizon. An extension of the code [105] was later used to verify that scalar field collapse in six dimensions continues to display critical phenomena.

Hamadé and Stewart [135Jump To The Next Citation Point] also applied a double null code to study critical phenomena. In order to obtain the accuracy necessary to confirm Choptuik’s results they developed the first example of a characteristic grid with AMR. They did this with both the standard Berger and Oliger algorithm and their own simplified version, with both versions giving indistinguishable results. Their simulations of critical collapse of a massless scalar field agreed with Choptuik’s values for the universal parameters governing mass scaling and displayed the echoing associated with discrete self-similarity. Hamadé, Horne, and Stewart [134] extended this study to the spherical collapse of an axion/dilaton system and found in this case that self-similarity was a continuous symmetry of the critical solution.

Brady, Chambers, and Gonçalves [55] used Garfinkle’s [104Jump To The Next Citation Point] double null algorithm to investigate the effect of a massive scalar field on critical phenomena. The introduction of a mass term in the scalar wave equation introduces a scale to the problem, which suggests that the critical point behavior might differ from the massless case. They found that there are two different regimes depending on the ratio of the Compton wavelength 1∕m of the scalar mass to the radial size λ of the scalar pulse used to induce collapse. When λm < < 1, the critical solution is the one found by Choptuik in the m = 0 case, corresponding to a type II phase transition. However, when λm >> 1, the critical solution is an unstable soliton star (see [216Jump To The Next Citation Point]), corresponding to a type I phase transition where black hole formation turns on at a finite mass.

A code based upon Bondi coordinates, developed by Husa and his collaborators [149], has been successfully applied to spherically symmetric critical collapse of a nonlinear σ-model coupled to gravity. Critical phenomena cannot be resolved on a static grid based upon the Bondi r-coordinate. Instead, the numerical techniques of Garfinkle were adopted by using a dynamic grid following the ingoing null rays and by recycling radial grid points. They studied how coupling to gravity affects the critical behavior previously observed by Bizoń [51] and others in the Minkowski space version of the model. For a wide range of the coupling constant, they observe discrete self-similarity and typical mass scaling near the critical solution. The code is shown to be second order accurate and to give second order convergence for the value of the critical parameter.

The first characteristic code in Bondi coordinates for the self-gravitating scalar wave problem was constructed by Gómez and Winicour [122Jump To The Next Citation Point]. They introduced a numerical compactification of ℐ+ for the purpose of studying effects of self-gravity on the scalar radiation, particularly in the high amplitude limit of the rescaling Φ → aΦ. As a → ∞, the red shift creates an effective boundary layer at + ℐ which causes the Bondi mass MB and the scalar field monopole moment Q to be related by √ -- MB ∼ π|Q |∕ 2, rather than the quadratic relation of the weak field limit [122]. This could also be established analytically so that the high amplitude limit provided a check on the code’s ability to handle strongly nonlinear fields. In the small amplitude case, this work incorrectly reported that the radiation tails from black hole formation had an exponential decay characteristic of quasinormal modes rather than the polynomial 1∕t or 1 ∕t2 falloff expected from Price’s [198Jump To The Next Citation Point] work on perturbations of Schwarzschild black holes. In hindsight, the error here was not having confidence to run the code sufficiently long to see the proper late time behavior.

Gundlach, Price, and Pullin [129130] subsequently reexamined the issue of power law tails using a double null code similar to that developed by Goldwirth and Piran. Their numerical simulations verified the existence of power law tails in the full nonlinear case, thus establishing consistency with analytic perturbative theory. They also found normal mode ringing at intermediate time, which provided reassuring consistency with perturbation theory and showed that there is a region of spacetime where the results of linearized theory are remarkably reliable even though highly nonlinear behavior is taking place elsewhere. These results have led to a methodology that has application beyond the confines of spherically symmetric problems, most notably in the “close approximation” for the binary black hole problem [199Jump To The Next Citation Point]. Power law tails and quasinormal ringing have also been confirmed using Cauchy evolution [171Jump To The Next Citation Point].

The study of the radiation tail decay of a scalar field was subsequently extended by Gómez, Schmidt, and Winicour [125Jump To The Next Citation Point] using a characteristic code. They showed that the Newman–Penrose constant [181] for the scalar field determines the exponent of the power law (and not the static monopole moment as often stated). When this constant is non-zero, the tail decays as 1∕t on ℐ+, as opposed to the 1∕t2 decay for the vanishing case. (They also found −n t log t corrections, in addition to the exponentially decaying contributions of the quasinormal modes.) This code was also used to study the instability of a topological kink in the configuration of the scalar field [23]. The kink instability provides the simplest example of the turning point instability [152226] which underlies gravitational collapse of static equilibria.

Brady and Smith [57] have demonstrated that characteristic evolution is especially well adapted to explore properties of Cauchy horizons. They examined the stability of the Reissner–Nordström Cauchy horizon using an Einstein–Klein–Gordon code based upon advanced Bondi coordinates (v,r) (where the hypersurfaces v = const are ingoing null hypersurfaces). They studied the effect of a spherically symmetric scalar pulse on the spacetime structure as it propagates across the event horizon. Their numerical methods were patterned after the work of Goldwirth and Piran [109Jump To The Next Citation Point], with modifications of the radial grid structure that allow deep penetration inside the black hole. In accord with expectations from analytic studies, they found that the pulse first induces a weak null singularity on the Cauchy horizon, which then leads to a crushing spacelike singularity as r → 0. The null singularity is weak in the sense that an infalling observer experiences a finite tidal force, although the Newman–Penrose Weyl component Ψ2 diverges, a phenomenon known as mass inflation [193]. These results confirm the earlier result of Gnedin and Gnedin [108] that a central spacelike singularity would be created by the interaction of a charged black hole with a scalar field, in accord with a physical argument by Penrose [188] that a small perturbation undergoes an infinite redshift as it approaches the Cauchy horizon.

Burko [60Jump To The Next Citation Point] has confirmed and extended these results, using a code based upon double null coordinates which was developed with Ori [61] in a study of tail decay. He found that in the early stages the perturbation of the Cauchy horizon is weak and in agreement with the behavior calculated by perturbation theory.

Brady, Chambers, Krivan, and Laguna [56] have found interesting effects of a non-zero cosmological constant Λ on tail decay by using a characteristic Einstein–Klein–Gordon code to study the effect of a massless scalar pulse on Schwarzschild–de Sitter and Reissner–Nordström–de Sitter spacetimes. First, by constructing a linearized scalar evolution code, they show that scalar test fields with ℓ ⁄= 0 have exponentially decaying tails, in contrast to the standard power law tails for asymptotically flat spacetimes. Rather than decaying, the monopole mode asymptotes at late time to a constant, which scales linearly with Λ, in contrast to the standard no-hair result. This unusual behavior for the ℓ = 0 case was then independently confirmed with a nonlinear spherical characteristic code.

Using a combination of numerical and analytic techniques based upon null coordinates, Hod and Piran have made an extensive series of investigations of the spherically symmetric charged Einstein–Klein–Gordon system dealing with the effect of charge on critical gravitational collapse [142] and the late time tail decay of a charged scalar field on a Reissner–Nordström black hole [143146144145]. These studies culminated in a full nonlinear investigation of horizon formation by the collapse of a charged massless scalar pulse [147Jump To The Next Citation Point]. They track the formation of an apparent horizon which is followed by a weakly singular Cauchy horizon which then develops a strong spacelike singularity at r = 0. This is in complete accord with prior perturbative results and nonlinear simulations involving a pre-existing black hole. Oren and Piran [182] increased the late time accuracy of this study by incorporating an adaptive grid for the retarded time coordinate u, with a refinement criterion to maintain Δr ∕r = const. The accuracy of this scheme is confirmed through convergence tests as well as charge and constraint conservation. They were able to observe the physical mechanism which prohibits black hole formation with charge to mass ration Q ∕M > 1. Electrostatic repulsion of the outer parts of the scalar pulse increases relative to the gravitational attraction and causes the outer portion of the charge to disperse to larger radii before the black hole is formed. Inside the black hole, they confirm the formation of a weakly singular Cauchy horizon which turns into a strong spacelike singularity, in accord with other studies.

Hod extended this combined numerical-analytical double null approach to investigate higher order corrections to the dominant power law tail [140], as well as corrections due to a general spherically symmetric scattering potential [139] and due to a time dependent potential [141]. He found (logt)∕t modifications to the leading order tail behavior for a Schwarzschild black hole, in accord with earlier results of Gómez et al. [125]. These modifications fall off at a slow rate so that a very long numerical evolution (t ≈ 3000 M)is necessary to cleanly identify the leading order power law decay.

The foregoing numerical-analytical work based upon characteristic evolution has contributed to a very comprehensive classical treatment of spherically symmetric gravitational collapse. Sorkin and Piran [225] have investigated the question of quantum corrections due to pair creation on the gravitational collapse of a charged scalar field. For observers outside the black hole, several analytic studies have indicated that such pair-production can rapidly diminish the charge of the black hole. Sorkin and Piran apply the same double-null characteristic code used in studying the classical problem [147] to evolve across the event horizon and observe the quantum effects on the Cauchy horizon. The quantum electrodynamic effects are modeled in a rudimentary way by a nonlinear dielectric 𝜖 constant that limits the electric field to the critical value necessary for pair creation. The back-reaction of the pairs on the stress-energy and the electric current are ignored. They found that quantum effects leave the classical picture of the Cauchy horizon qualitatively intact but that they shorten its “lifetime” by hastening the conversion of the weak null singularity into a strong spacelike singularity.

The Southampton group has constructed a {1 + 1}-dimensional characteristic code for spacetimes with cylindrical symmetry [77Jump To The Next Citation Point87Jump To The Next Citation Point]. The original motivation was to use it as the exterior characteristic code in a test case of CCM (see Section 5.5.1 for the application to matching). Subsequently, Sperhake, Sjödin, and Vickers [223227] modified the code into a global characteristic version for the purpose of studying cosmic strings, represented by massive scalar and vector fields coupled to gravity. Using a Geroch decomposition [106] with respect to the translational Killing vector, they reduced the global problem to a {2 + 1}-dimensional asymptotically flat spacetime, so that ℐ+ could be compactified and included in the numerical grid. Rather than the explicit scheme used in CCM, the new version employs an implicit, second order in space and time, Crank–Nicholson evolution scheme. The code showed long term stability and second order convergence in vacuum tests based upon exact Weber–Wheeler waves [247Jump To The Next Citation Point] and Xanthopoulos’ rotating solution [252], and in tests of wave scattering by a string. The results show damped ringing of the string after an incoming Weber–Wheeler pulse has excited it and then scattered to ℐ+. The ringing frequencies are independent of the details of the pulse but are inversely proportional to the masses of the scalar and vector fields.

Frittelli and Gómez [101] have cast the spherically symmetric Einstein-Klein-Gordon problem in symmetric hyperbolic form, where in a Bondi-Sachs gauge the fundamental variables are the scalar field, lapse and shift. The Bondi-Sachs gauge conditions relate the usual ADM variables (the 3-metric and extrinsic curvature) to the lapse and shift, which obey simpler evolution equations. The resulting Cauchy problem is well-posed and the outer boundary condition is constraint preserving (although whether the resulting IBVP is well-posed is not addressed, i.e. whether the boundary condition is dissipative). A numerical evolution algorithm based upon the system produces a stable simulation of a scalar pulse Φ scattering off a black hole. The initial data for the pulse satisfies ∂tΦ = 0 so, as expected, it contains an ingoing part, which crosses the horizon, and an outgoing part, which leaves the grid at the outer boundary with a small amount of back reflection.

3.1.1 Adaptive mesh refinement

The goal of computing waveforms from relativistic binaries, such as a neutron star or stellar mass back hole spiraling into a supermassive black hole, requires more than a stable convergent code. It is a delicate task to extract a waveform in a spacetime in which there are multiple length scales: the size of the supermassive black hole, the size of the star, the wavelength of the radiation. It is commonly agreed that some form of mesh refinement is essential to attack this problem. Mesh refinement was first applied in characteristic evolution to solve specific spherically symmetric problems regarding critical phenomena and singularity structure [10413560].

Pretorius and Lehner [197Jump To The Next Citation Point] have presented a general approach for AMR to a generic characteristic code. Although the method is designed to treat 3D simulations, the implementation has so far been restricted to the Einstein–Klein–Gordon system in spherical symmetry. The 3D approach is modeled after the Berger and Oliger AMR algorithm for hyperbolic Cauchy problems, which is reformulated in terms of null coordinates. The resulting characteristic AMR algorithm can be applied to any unigrid characteristic code and is amenable to parallelization. They applied it to the problem of a spherically symmetric massive Klein–Gordon field propagating outward from a black hole. The non-zero rest mass restricts the Klein–Gordon field from propagating to infinity. Instead it diffuses into higher frequency components which Pretorius and Lehner show can be resolved using AMR but not with a comparison unigrid code.

  Go to previous page Go up Go to next page