3.2 {2 + 1} Codes3 Characteristic Evolution Codes3 Characteristic Evolution Codes

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 the same vein, it is fair to say that the general system of hyperbolic partial differential equations in one spatial dimension is a solved problem. At least, it seems to be true in general relativity.

One of the earliest characteristic evolution codes, constructed by Corkill and Stewart [47, 133], 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 [96] collision of impulsive (tex2html_wrap_inline1719 -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 unforseen 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, which is beginning to make significant contributions to general relativistic astrophysics, is reviewed in Sec.  5 .

The synergy between analytic and computational approaches has already 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 existence and uniqueness of solutions to this problem [41, 40, 43, 42]. He showed that weak initial data evolve to Minkowski space asymptotically in time, but that sufficiently strong data form a horizon, with nonzero Bondi mass. In the latter case, he showed that the geometry is asymptotically Schwarzschild in the approach to tex2html_wrap_inline1721 (future time infinity) from outside the horizon, thus establishing a rigorous version of the no-hair theorem. What this analytic tour-de-force did not reveal was the remarkable critical behavior in the transition between these two regimes, which was discovered by Choptuik [38, 39Jump To The Next Citation Point In The Article] by computational simulation based upon Cauchy evolution. This phenomena 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 a Living Review on Critical Phenomena in Gravitational Collapse by C. Gundlach [79].

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 tex2html_wrap_inline1695 and a surface area coordinate r, the metric is


Smoothness at r =0 allows the coordinate conditions


The field equations consist of the wave equation tex2html_wrap_inline1729 for the scalar field and two hypersurface equations for the metric functions:



The wave equation can be expressed in the form


where tex2html_wrap_inline1731 and tex2html_wrap_inline1733 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 tex2html_wrap_inline1735 at initial retarded time tex2html_wrap_inline1673 .

Because any two-dimensional geometry is conformally flat, the surface integral of tex2html_wrap_inline1739 over a null parallelogram tex2html_wrap_inline1741 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 [76Jump To The Next Citation Point In The Article]. 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 Fig.  2 . Upon integration of (7Popup Equation), curvature introduces an area integral correction to the flat space null parallelogram relation between the values of g at the vertices:



Click on thumbnail to view image

Figure 2: The null parallelogram. After computing the field at point N, the algorithm marches the computation to tex2html_wrap_inline1655 by shifting the corners by tex2html_wrap_inline1657, tex2html_wrap_inline1659, tex2html_wrap_inline1661, tex2html_wrap_inline1663 .

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 Bondi coordinates as in Eq. (3Popup 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 Eq. (8Popup Equation) requires no interpolation. As a result, in flat space, where the right hand side of Eq. (8Popup 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 by evaluating the integrand at the center of tex2html_wrap_inline1741 .

The identity (8Popup Equation) gives rise to the following explicit marching algorithm, indicated in Fig.  2 . Let the null parallelogram lie at some fixed tex2html_wrap_inline1775 and tex2html_wrap_inline1777 and span adjacent retarded time levels tex2html_wrap_inline1673 and tex2html_wrap_inline1781 . Imagine for now that the points N, E, S and W lie on the spatial grid, with tex2html_wrap_inline1791 . If g has been determined on the initial cone tex2html_wrap_inline1673, which contains the points E and S, and radially outward from the origin to the point W on the next cone tex2html_wrap_inline1781, then Eq. (8Popup Equation) determines g at the next radial grid point N in terms of an integral over tex2html_wrap_inline1741 . The integrand can be approximated to second order, i.e. to tex2html_wrap_inline1811, by evaluating it at the center of tex2html_wrap_inline1741 . 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 tex2html_wrap_inline1823 at the center of tex2html_wrap_inline1741 can be approximated to second order in terms of values of V at points where it can be determined by integrating the hypersurface equations (5Popup Equation) and (6Popup 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 tex2html_wrap_inline1781 level, i.e. point n in Fig.  2 . 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 tex2html_wrap_inline1845, 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 tex2html_wrap_inline1741 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, tex2html_wrap_inline1863 and V at the vertices of tex2html_wrap_inline1741 are approximated to second order accuracy by interpolating between grid points. However, cancellations arise between these four interpolations so that Eq. (8Popup Equation) is satisfied to fourth order accuracy. The net result is that the finite difference version of (8Popup Equation) steps g radially outward one zone with an error of fourth order in grid size, tex2html_wrap_inline1871 . In addition, the smoothness conditions (4Popup Equation) can be incorporated into the startup for the numerical integrations for V and tex2html_wrap_inline1863 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 tex2html_wrap_inline1881, after compounding errors from tex2html_wrap_inline1883 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 tex2html_wrap_inline1887, where tex2html_wrap_inline1889 . This condition can be built into the code using the value tex2html_wrap_inline1891, corresponding to the maximum of V / r at tex2html_wrap_inline1655 . The strongest restriction on the time step then arises just before the formation of a horizon, where tex2html_wrap_inline1897 at tex2html_wrap_inline1655 . 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 tex2html_wrap_inline1655 must be dropped in order to evolve across the horizon due to the lack of a nonsingular compactified version of future time infinity tex2html_wrap_inline1721 .

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 used by Goldwirth and Piran in a study of cosmic censorship [67Jump To The Next Citation Point In The Article] based upon the spherically symmetric gravitational collapse of a scalar field. Their early study lacked the sensitivity of adaptive mesh refinement which later enabled Choptuik to discover the critical phenomena appearing in this problem. Subsequent work by Marsa and Choptuik [103Jump To The Next Citation Point In The Article] combined the use of the null related ingoing Eddington-Finklestein coordinates with singularity excision to construct a 1D code that ``runs forever''. Later, Garfinkle  [63] 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 the remaining grid points. This technique is crucial to resolving the critical phenomena associated with an tex2html_wrap_inline1907 sized horizon. An extension of the code [64] was later used to verify that scalar field collapse in six dimensions continues to display critical phenomena.

Hamadé and Stewart [84] 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 adaptive mesh refinement (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 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 [83] extended this study to the spherical collapse of an axion/dilaton system and found in this case that the self-similarity was a continuous symmetry of the critical solution.

A code based upon Bondi coordinates, developed by S. Husa and his collaborators [88], has been successfully applied to spherically symmetric critical collapse of a nonlinear tex2html_wrap_inline1909 -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 P. Bizon [32] 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 [76Jump To The Next Citation Point In The Article]. They introduced a numerical compactification of tex2html_wrap_inline1655 for the purpose of studying effects of self-gravity on the scalar radiation, particularly in the high amplitude limit of the rescaling tex2html_wrap_inline1915 . As tex2html_wrap_inline1917, the red shift creates an effective boundary layer at tex2html_wrap_inline1655 which causes the Bondi mass tex2html_wrap_inline1921 and the scalar field monopole moment Q to be related by tex2html_wrap_inline1925, rather than the quadratic relation of the weak field limit [76]. 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 tex2html_wrap_inline1929 falloff expected from Price's [116Jump To The Next Citation Point In The Article] 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 [80, 81] 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 [117Jump To The Next Citation Point In The Article]. Power law tails and quasinormal ringing were also confirmed using Cauchy evolution [103Jump To The Next Citation Point In The Article].

The study of the radiation tail decay of a scalar field was subsequently extended by Gómez, Schmidt and Winicour [75] using a characteristic code. They showed that the Newman-Penrose constant [108] for the scalar field determines the exponent of the power law (not the static monopole moment as often stated). When this constant is non-zero, the tail decays as 1/ t, as opposed to the tex2html_wrap_inline1929 decay for the vanishing case. (There are also tex2html_wrap_inline1935 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 [13]. The kink instability provides the simplest example of the turning point instability [90, 129] which underlies gravitational collapse of static equilibria.

The Southampton group has constructed a {1 + 1}-dimensional characteristic code for spacetimes with cylindrical symmetry [46Jump To The Next Citation Point In The Article, 54Jump To The Next Citation Point In The Article]. The original motivation was to use it as the exterior characteristic code in a test case of Cauchy-characteristic matching (CCM). (See Sec.  4.4.1 of this review for the application to matching.) However, U. Sperhake, K. R. P. Sjödin and J. A. Vickers [130, 131] have recently 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 [65] with respect to the translational Killing vector, the global problem is reduced to a {2 + 1}-dimensional asymptotically flat spacetime, so that tex2html_wrap_inline1655 can 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 [152Jump To The Next Citation Point In The Article], Xanthopoulos' rotating solution [155], and in tests of wave scattering by a string. The results show damped ringing of the string after a incoming Weber-Wheeler pulse has excited it and scattered out to tex2html_wrap_inline1655 . The ringing frequencies are independent of the details of the pulse but are inversely proportional to the masses of the scalar and vector fields.

The group at the Universidad de Oriente in Venezuela has applied characteristic methods to study the self-similar collapse of spherical matter and charge distributions [11, 14, 12Jump To The Next Citation Point In The Article]. The assumption of self-similarity reduces the problem to a system of ODE's, subject to boundary conditions determined by matching to an exterior Reissner-Nordstrom-Vaidya solution. Heat flow in the internal fluid is balanced at the surface by the Vaidya radiation. Their solutions illustrate how a nonzero total charge can halt gravitational collapse and produce a final stable equilibrium [12]. It is interesting that the pressure vanishes in the final equilibrium state so that hydrostatic support is completely supplied by Coulomb repulsion.

3.2 {2 + 1} Codes3 Characteristic Evolution Codes3 Characteristic Evolution Codes

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