1 Introduction

The singularity theorems [178, 104, 105, 106] state that Einstein’s equations will not evolve regular initial data arbitrarily far into the future or the past. An obstruction such as infinite curvature or the termination of geodesics will always arise to stop the evolution somewhere. The simplest, physically relevant solutions representing for example a homogeneous, isotropic universe (Friedmann-Robertson-Walker (FRW)) or a spherically symmetric black hole (Schwarzschild) contain space-like infinite curvature singularities. Although, in principle, the presence of a singularity could lead to unpredictable measurements for a physically realistic observer, this does not happen for these two solutions. The surface of last scattering of the cosmic microwave background in the cosmological case and the event horizon in the black hole (BH) case effectively hide the singularity from present day, external observers. The extent to which this “hidden” singularity is generic, and the types of singularities that appear in generic spacetimes, remain major open questions in general relativity. The questions arise quickly, since other exact solutions to Einstein’s equations have singularities which are quite different from those described above. For example, the charged BH (Reissner-Nordstrom solution) has a time-like singularity. It also contains a Cauchy horizon (CH) marking the boundary of predictability of space-like initial data supplied outside the BH. A test observer can pass through the CH to another region of the extended spacetime. More general cosmologies can exhibit singularity behavior different from that in FRW. The Big Bang in FRW is classified as an asymptotically velocity term dominated (AVTD) singularity [69, 119], since any spatial curvature term in the Hamiltonian constraint becomes negligible compared to the square of the expansion rate as the singularity is approached. However, some anisotropic, homogeneous models exhibit Mixmaster dynamics (MD) [15, 139] and are not AVTD — the influence of the spatial scalar curvature can never be neglected. For more rigorous discussions of the classification and properties of the types of singularities see [71, 176].

Once the simplest, exactly solvable models are left behind, understanding of the singularity becomes more difficult. There has been significant analytic progress [179, 143, 163]. However, such methods yield either detailed knowledge of unrealistic, simplified (usually by symmetries) spacetimes or powerful, general results that do not contain details. To overcome these limitations, one might consider numerical methods to evolve realistic spacetimes to the point where the properties of the singularity may be identified. Of course, most of the effort in numerical relativity applied to BH collisions has addressed the avoidance of singularities [74]. One wishes to keep the computational grid in the observable region outside the horizon. Much less computational effort has focused on the nature of the singularity itself. Numerical calculations, even more than analytic ones, require finite values for all quantities. Ideally then, one must describe the singularity by the asymptotic non-singular approach to it. A numerical method which can follow the evolution into this asymptotic regime will then yield information about the singularity. Since the numerical study must begin with a particular set of initial data, the results can never have the force of mathematical proof. One may hope, however, that such studies will provide an understanding of the “phenomenology” of singularities that will eventually guide and motivate rigorous results. Some examples of the interplay between analytic and numerical results and methods will be given here.

In the following, we shall consider examples of numerical study of singularities both for asymptotically flat (AF) spacetimes and for cosmological models. These examples have been chosen to illustrate primarily numerical studies whose focus is the nature of the singularity itself. In the AF context, we shall consider two questions. The first is whether or not naked singularities exist for realistic matter sources.

One approach has been to explore highly non-spherical collapse looking for spindle or pancake singularities. If the formation of an event horizon requires a limit on the aspect ratio of the matter [175], such configurations may yield a naked singularity. Recent analytic results suggest that one must go beyond the failure to observe an apparent horizon to conclude that a naked singularity has formed [179].

Another approach is to probe the limits between initial configurations which lead to black holes and those which yield no singularity at all (i.e. flat spacetime plus radiation) to explore the singularity as the BH mass goes to zero. This quest led naturally to the discovery of critical behavior in the collapse of a scalar field [59]. The critical (Choptuik) solution is a zero mass naked singularity (visible from null infinity). It is a counterexample to the cosmic censorship conjecture [102]. However, it is a non-generic one since (in addition to the fine-tuning required for this critical solution) Christodoulou has shown [62] that for the spherically symmetric Einstein-scalar field equations, there always exists a perturbation that will convert a solution with a naked singularity to one with a black hole. Reviews of critical phenomena in gravitational collapse can be found in [33, 93, 96].

The other question which is now beginning to yield to numerical attack involves the stability of the Cauchy horizon (CH) in charged or rotating black holes. It has been conjectured [178, 56] that a real observer, as opposed to a test mass, cannot pass through the CH since realistic perturbed spacetimes will convert the CH to a strong spacelike singularity [176]. Numerical studies [40, 67, 47] show that a weak, null singularity forms first, as had been predicted [157, 151].

In cosmology, we shall consider both the behavior of the Mixmaster model and the issue of whether or not its properties are applicable to generic cosmological singularities. Although numerical evolution of the Mixmaster equations has a long history, developments in the past decade were motivated by inconsistencies between the known sensitivity to initial conditions and standard measures of the chaos usually associated with such behavior [145, 164, 166, 19, 76, 45, 111, 160]. Most recently, a coordinate invariant characterization of Mixmaster chaos has been formulated [64] and a new extremely fast and accurate algorithm for Mixmaster simulations developed [28].

Belinskii, Khalatnikov, and Lifshitz (BKL) long ago claimed [11, 12, 13, 15, 14] that it is possible to formulate the generic cosmological solution to Einstein’s equations near the singularity as a Mixmaster universe at every point. While others have questioned the validity of this claim [8], it is only very recently that evidence for oscillatory behavior in the approach to the singularity of spatially inhomogeneous cosmologies has been obtained [181, 31]. We shall discuss a numerical program to address this issue [27].

2 Singularities in AF Spacetimes

While I have divided this topic into three subsections, there is considerable overlap. The primary questions can be formulated as the cosmic censorship conjecture. The weak cosmic censorship conjecture [154] requires a singularity formed from regular, asymptotically flat initial data to be hidden from an external observer by an event horizon. An excellent review of the meaning and status of weak cosmic censorship has been given by Wald [179]. Counter examples have been known for a long time but tend to be dismissed as unrealistic in some way. The strong form of the cosmic censorship conjecture [155] forbids timelike singularities, even within black holes.

2.1 Naked Singularities and the Hoop Conjecture

2.1.1 Overview

Perhaps the first numerical approach to study the cosmic censorship conjecture consisted of attempts to create naked singularities. Many of these studies were motivated by Thorne’s “hoop conjecture” [175] that collapse will yield a black hole only if a mass M is compressed to a region with circumference C ≤ 4πM in all directions. (As is discussed by Wald [179], one runs into difficulties in any attempt to formulate the conjecture precisely. For example, how does one define C and M, especially if the initial data are not at least axially symmetric?) If the hoop conjecture is true, naked singularities may form if collapse can yield C ≥ 4πM in some direction. The existence of a naked singularity is inferred from the absence of an apparent horizon (AH) which can be identified locally by following null geodesics. Although a definitive identification of a naked singularity requires the event horizon (EH) to be proven to be absent, to identify an EH requires knowledge of the entire spacetime. While one finds an AH within an EH [120, 121], it is possible to construct a spacetime slicing which avoids the AH even though an EH is present [180]. Methods to find an EH in a numerically determined spacetime have only recently become available and have not been applied to this issue [132, 136].

2.1.2 Naked Spindle Singularities?

In the best known attempt to produce naked singularities, Shapiro and Teukolsky (ST) [169] considered collapse of prolate spheroids of collisionless gas. (Nakamura and Sato [146] had previously studied the collapse of non-rotating deformed stars with an initial large reduction of internal energy and apparently found spindle or pancake singularities in extreme cases.) ST solved the general relativistic Vlasov equation for the particles along with Einstein’s equations for the gravitational field. Null geodesics were followed to identify an AH if present. The curvature invariant I = Rμνρσ Rμνρσ was also computed. They found that an AH (and presumably a BH) formed if C ≤ 4πM < 1 everywhere but no AH (and presumably a naked singularity) in the opposite case. In the latter case, the evolution (not surprisingly) could not proceed past the moment of formation of the singularity. In a subsequent study, ST [170] also showed that a small amount of rotation (counter rotating particles with no net angular momentum) does not prevent the formation of a naked spindle singularity. However, Wald and Iyer [180] have shown that the Schwarzschild solution has a time slicing whose evolution approaches arbitrarily close to the singularity with no AH in any slice, (but of course, with an EH in the spacetime). This may mean that there is a chance that the increasing prolateness found by ST in effect changes the slicing to one with no apparent horizon just at the point required by the hoop conjecture. While, on the face of it, this seems unlikely, Tod gives an example where the AH does not form on a chosen constant time slice — but rather different portions form at different times. He argues that a numerical simulation might be forced by the singularity to end before the formation of the AH is complete. Such an AH would not be found by the simulations [177]. In response to such a possibility, Shapiro and Teukolsky considered equilibrium sequences of prolate relativistic star clusters [171]. The idea is to counter the possibility that an EH might form after the simulation must stop. If an equilibrium configuration is non-singular, it cannot contain an EH, since singularity theorems say that an EH implies a singularity. However, a sequence of nonsingular equilibria with rising I ever closer to the spindle singularity would lend support to the existence of a naked spindle singularity since one can approach the singular state without formation of an EH. They constructed this sequence and found that the singular end points were very similar to their dynamical spindle singularity. Wald believes, however, that it is likely that the ST slicing is such that their singularities are not naked — the AH is present but has not yet appeared in their time slices [179].

Another numerical study of the hoop conjecture was made by Chiba et al [58]. Rather than a dynamical collapse model, they searched for AH’s in analytic initial data for discs, annuli, and rings. Previous studies of this type were done by Nakamura et al [147] with oblate and prolate spheroids and by Wojtkiewicz [182] with axisymmetric singular lines and rings. The summary of their results is that an AH forms if C ≤ 4πM ≤ 1.26. (Analytic results due to Barrabes et al [5, 4] and Tod [177] give similar quantitative results with different initial data classes and (possibly) definition of C.)

There is strong analytical evidence against the development of spindle singularities. It has been shown by Chrusciel and Moncrief that strong cosmic censorship holds in AF electrovac solutions which admit a G2 symmetric Cauchy surface [24]. The evolutions of these highly nonlinear equations are in fact non-singular.

2.1.3 Going Further

Motivated by ST’s results [169], Echeverria [70] numerically studied the properties of the naked singularity that is known to form in the collapse of an infinite, cylindrical dust shell [175]. While the asymptotic state can be found analytically, the approach to it must be followed numerically. The analytic asymptotic solution can be matched to the numerical one (which cannot be followed all the way to the collapse) to show that the singularity is strong (an observer experiences infinite stretching parallel to the symmetry axis and squeezing perpendicular to the symmetry axis). A burst of gravitational radiation emitted just prior to the formation of the singularity stretches and squeezes in opposite directions to the singularity. This result for dust conflicts with rigorously nonsingular solutions for the electrovac case [24]. One wonders then if dust collapse gives any information about singularities of the gravitational field.

Nakamura et al (NSN) [148] conjectured that even if naked spindle singularities could exist, they would either disappear or become black holes. This demise of the naked singularity would be caused by the back reaction of the gravitational waves emitted by it. While NSN proposed a numerical test of their conjecture, they believed it to be beyond the current generation of computer technology.

The results of all these searches for naked spindle singularities are controversial, but they could be resolved if the presence or absence of the EH could be determined. One could demonstrate numerically whether or not Wald’s view of ST’s results is correct by using existing EH finders [132, 136] in a relevant simulation.

2.2 Critical Behavior in Collapse

2.2.1 Gravitational collapse simulations

We now consider an effect that was originally found by Choptuik [59] in a numerical study of the collapse of a spherically symmetric massless scalar field. For recent reviews see [93, 96]. We note that this is the first completely new phenomenon in general relativity to be discovered by numerical simulation. In collapse of a scalar field, essentially two things can happen: Either a black hole (BH) forms, or the scalar waves pass through each other and disperse. Choptuik discovered that for any 1-parameter set of initial data labeled by p, there is a critical value p* such that p > p* yields a BH. He found

$${M_{BH}} \approx {C_F}{(p - p*)^\gamma}$$
(1)

where MBH is the mass of the eventual BH. The constant CF depends on the parameter of the initial data that is selected, but γ ≈ .37 is the same for all choices. Furthermore, in terms of logarithmic variables ρ = ln r + κ = ln(T0* − T0) + κ (T0 is the proper time of an observer at r = 0 with T0* the finite proper time at which the critical evolution concludes and κ is a constant which scales r), the waveform X repeats (echoes) at intervals of Δ in τ if ρ is rescaled to ρ − Δ, i.e. X(ρ − Δ, τ − Δ) ≈ X(ρ, τ). The scaling behavior (1) demonstrates that the minimum BH mass (for bosons) is zero. The critical solution itself is a counter-example to cosmic censorship, (since the formation of the zero mass BH causes high curvature regions to become visible at r = ∞). (See, e.g., the discussion in Hirschmann and Eardley [109].) The numerical demonstration of this feature of the critical solution was provided by Hamadé and Stewart [102]. This result caused Hawking to pay off a bet to Preskill and Thorne [44, 123].

Soon after this discovery, scaling and critical phenomena were found in a variety of contexts. Abrahams and Evans [1] discovered the same phenomenon in axisymmetric gravitational wave collapse with a different value of Δ and, to within numerical error, the same value of γ. (Note that the rescaling of r with eΔ ≈ 30 required Choptuik to use adaptive mesh refinement (AMR) to distinguish subsequent echoes. Abrahams and Evans’ smaller Δ (eΔ ≈ 1.8) allowed them to see echoing with their 2 + 1 code without AMR.) Garfinkle [81] confirmed Choptuik’s results with a completely different algorithm that does not require AMR. He used Goldwirth and Piran’s [86] method of simulating Christodoulou’s [61] formulation of the spherically symmetric scalar field in null coordinates. This formulation allowed the grid to be automatically rescaled by choosing the edge of the grid to be the null ray that just hits the central observer at the end of the critical evolution. (Missing points of null rays that cross the central observer’s world line are replaced by interpolation between those that remain.) Hamadé and Stewart [102] have also repeated Choptuik’s calculation using null coordinates and AMR. They are able to achieve greater accuracy and find γ = .374.

2.2.2 Critical Solutions as an Eigenvalue Problem

Evans and Coleman [72] realized that self-similar rather than self-periodic collapse might be more tractable both numerically (since ODE’s rather than PDE’s are involved) and analytically. They discovered that a collapsing radiation fluid had that desirable property. (Note that self-similarity (homothetic motion) is incompatible with AF. However, most of the action occurs in the center so that a match of the self-similar inner region to an outer AF one should always be possible.) In a series of papers, Hirschmann and Eardley [108, 109] developed a (numerical) self-similar solution to the spherically symmetric complex scalar field equations. These are ODE’s with too many boundary conditions causing a solution to exist only for certain fixed values of Δ. Numerical solution of this eigenvalue problem allows very accurate determination of Δ. The self-similarity also allows accurate calculation of γ as follows: The critical p = p* solution is unstable to a small change in p. At any time t (where t < 0 is increasing toward zero), the amplitude a of the perturbation exhibits power law growth:

$$a \propto (p - p*){(- t)^{- \kappa}}$$
(2)

where κ > 0. At any fixed t, larger a implies larger MBH. Equivalently, any fixed amplitude a = δ will be reached faster for larger eventual MBH. Scaling arguments give the dependence of MBH on the time at which any fixed amplitude is reached:

$${M_{BH}} \propto (- {t_1})$$
(3)

where

$$(p - p*){(- {t_1})^\kappa} \propto \delta.$$
(4)

Thus

$${M_{BH}} \propto {(p - {p^*})^{1/\kappa}}.$$
(5)

Therefore, one need only identify the growth rate of the unstable mode to obtain an accurate value of γ = 1/κ. It is not necessary to undertake the entire dynamical evolution or probe the space of initial data. Hirschmann and Eardley obtain γ = 0.387106 for the complex scalar field solution, while Koike et al [129] obtain γ = 0.35580192 for the Evans-Coleman solution. Although the similarities among the critical exponents γ in the collapse computations suggested a universal value, Maison [135] used these same scaling-perturbation methods to show that γ depends on the equation of state p = kρ of the fluid in the Evans-Coleman solution. Gundlach [92] used a similar approach to locate Choptuik’s critical solution accurately. This is much harder, due to its discrete self-similarity. Gundlach reformulates the model as a nonlinear hyperbolic boundary value problem with eigenvalue Δ and finds Δ = 3.4439. As with the self-similar solutions described above, the critical solution is found directly without the need to perform a dynamical evolution or explore the space of initial data. Hara et al extended the renormalization group approach of [129] to the discretely-self-similar case [103]. (For a recent application of renormalization group methods to cosmology see [118].)

2.2.3 Recent Results

Recently, Gundlach [95] completed his eigenvalue analysis of the Choptuik solution to find the growth rate of the unstable mode to be γ = .374 ± .001. He also predicted a periodic “wiggle” in the Choptuik mass scaling relation. This was later observed numerically by Hod and Piran [114]. Self-similar critical behavior has been seen in string theory related axion-dilaton models [68, 101] and in the nonlinear σ-model [110]. Garfinkle and Duncan have shown that subcritical collapse of a spherically symmetric scalar field yields a scaling relation for the maximum curvature observed by the central observer with critical parameters that would be expected on the basis of those found for supercritical collapse [83].

Choptuik et al [60] have generalized the original Einstein-scalar field calculations to the Einstein-Yang-Mills (EYM) (for SU(2)) case. Here something new was found. Two types of behavior appeared depending on the initial data. In Type I, BH formation had a non-zero mass threshold. The critical solution is a regular, unstable solution to the EYM equations found previously by Bartnik and McKinnon [9]. In Type II collapse, the minimum BH mass is zero with the critical solution similar to that of Choptuik (with a different γ ≈ 0.20, Δ ≈ 0. 74). Gundlach has also looked at this case with the same results [94]. The Type I behavior arises when the collapsing system has a metastable static solution in addition to the Choptuik critical one [98].

Brady, Chambers and Goncalves [54, 37] conjectured that addition of a mass to the scalar field of the original model would break scale invariance and might yield a distinct critical behavior. They found numerically the same Type I and II “phases” seen in the EYM case. The Type II solution can be understood as perturbations of Choptuik’s original model with a small scalar field mass μ. Here small means that λμ ≪ 1 where λ is the spatial extent of the original nonzero field region. (The scalar field is well within the Compton wavelength corresponding to μ.) On the other hand, λμ ≫ 1 yields Type I behavior. The minimum mass critical solution is an unstable soliton of the type found by Seidel and Suen [168]. The massive scalar field can be treated as collapsing dust to yield a criterion for BH formation [87].

The Choptuik solution has also been found to be the critical solution for charged scalar fields [98, 113]. As p → p*, Q/M → 0 for the black hole. Q obeys a power law scaling. Numerical study of the critical collapse of collisionless matter (Einstein-Vlasov equations) has yielded a non-zero minimum BH mass [161]. Bizoh and Chmaj [34] have considered the critical collapse of skyrmions.

An astrophysical application of BH critical phenomena has been considered by Nimeyer and Jedamzik [149] and Yokoyama [183]. They consider its implications for primordial BH formation and suggest that it could be important.

2.2.4 Going Further

The question is then why these critical phenomena should appear in so many collapsing gravitational systems. The discrete self-similarity of Choptuik’s solution may be understood as scaling of a limit cycle [103]. (The self-similarity of other systems may be understood as scaling of a limit point.) Garfinkle [82] has conjectured that the scale invariance of Einstein’s equations might provide an underlying explanation for the self-similarity and discrete-self-similarity found in collapse.

Until recently, only Abrahams and Evans [1] had ventured beyond spherical symmetry. The first additional departure has been made by Gundlach [97]. He considered spherical and non-spherical perturbations of \(P = {1 \over 3}\rho\) perfect fluid collapse. Only the original (spherical) growing mode survived.

2.3 Nature of the Singularity in Charged or Rotating Black Holes

2.3.1 Overview

Unlike the simple singularity structure of the Schwarzschild solution, where the event horizon encloses a spacelike singularity at r = 0, charged and/or rotating BH’s have a much richer singularity structure. The extended spacetimes have an inner Cauchy horizon (CH) which is the boundary of predictability. To the future of the CH lies a timelike (ring) singularity [178]. Poisson and Israel [157, 158] began an analytic study of the effect of perturbations on the CH. Their goal was to check conjectures that the blue-shifted infalling radiation during collapse would convert the CH into a true singularity and thus prevent an observer’s passage into the rest of the extended regions. By including both ingoing and back-scattered outgoing radiation, they find for the Reissner-Nordström (RN) solution that the mass function (qualitatively Rαβγδ ∝ M/r3) diverges at the CH (mass inflation). However, Ori showed both for RN and Kerr [151, 152] that the metric perturbations are finite (even though Rμνρσ Rμνρσ diverges) so that an observer would not be destroyed by tidal forces (the tidal distortion would be finite) and could survive passage through the CH. A numerical solution of the Einstein-Maxwell-scalar field equations could test these perturbative results.

2.3.2 Numerical Studies

Gnedin and Gnedin [85] have numerically evolved the spherically symmetric Einstein-Maxwell with massless scalar field equations in a 2 + 2 formulation. The initial conditions place a scalar field on part of the RN event horizon (with zero field on the rest). An asymptotically null or spacelike singularity whose shape depends on the strength of the initial perturbation replaces the CH. For a sufficiently strong perturbation, the singularity is Schwarzschild-like. Although they claim to have found that the CH evolved to become a spacelike singularity, the diagrams in their paper show at least part of the final singularity to be null or asymptotically null in most cases.

Brady and Smith [40] used the Goldwirth-Piran formulation [86] to study the same problem. They assume the spacetime is RN for v < v0. They follow the evolution of the CH into a null singularity, demonstrate mass inflation, and support (with observed exponential decay of the metric component g) the validity of previous analytic results [157, 158, 151, 152] including the “weak” nature of the singularity that forms. They find that the observer hits the null CH singularity before falling into the curvature singularity at r = 0. Whether or not these results are in conflict with Gnedin and Gnedin [85] is unclear [35]. However, it has become clear that Brady and Smith’s conclusions are correct. The collapse of a scalar field in a charged, spherically symmetric spacetime causes an initial RN CH to become a null singularity, except at r = 0, where it is spacelike. The observer falling into the BH experiences (and potentially survives) the weak, null singularity [151, 152, 36] before the spacelike singularity forms. This has been confirmed by Droz [67] using a plane wave model of the interior and by Burko [47] using a collapsing scalar field. See also [49, 51].

Recently, numerical studies of the interiors of non-Abelian black holes have been carried out by Breitenlohner et al [41, 42] and by Gal’tsov et al [66, 78, 80, 79] (see also [184]). Although there appear to be conflicts between the two groups, the differences can be attributed to misunderstandings of each other’s notation [43]. The main results include an interesting oscilliatory behavior of the metric.

2.3.3 Going Further

Replacing the AF boundary conditions with Schwarzschild-de Sitter and RN-de Sitter BH’s was long believed to yield a counterexample to strong cosmic censorship. (See [137, 138, 156, 53] and references therein for background and extended discussions.) The stability of the CH is related to the decay tails of the radiating scalar field. Numerical studies recently determined these to be exponential [38, 53, 55] rather than power law as in AF spacetimes [50]. The decay tails of the radiation are necessary initial data for numerical study of CH stability [40] and are crucial to the development of the null singularity. Analytic studies had indicated that the CH is stable in RN-de Sitter BH’s for a restricted range of parameters near extremality. However, Brady et al [39] have shown (using linear perturbation theory) that, if one includes the backscattering of outgoing modes which originate near the event horizon, the CH is always unstable for all ranges of parameters. Thus RN-de Sitter BH’s appear not to be a counterexample to strong cosmic censorship. Numerical studies are needed to demonstrate the existence of a null singularity at the CH in nonlinear evolution.

As a potentially useful approach to the numerical study of singularities, we consider Hübner’s [115, 116, 117] numerical scheme to evolve on a conformal compactified grid using Friedrich’s formalism [77]. He considers the spherically symmetric scalar field model in a 2 + 2 formulation. So far, this code has been used to locate singularities and to identify Choptuik’s scaling [59].

3 Singularities in Cosmological Models

3.1 Singularities in Spatially Homogeneous Cosmologies

The generic singularity in spatially homogeneous cosmologies is reasonably well understood. The approach to it asymptotically falls into two classes. The first, called asymptotically velocity term dominated (AVTD) [69, 119], refers to a cosmology that approaches the Kasner (vacuum, Bianchi I) solution [124] as τ → ∞. (Spatially homogeneous universes can be described as a sequence of homogeneous spaces labeled by τ. Here we shall choose τ so that τ = ∞ coincides with the singularity.) An example of such a solution is the vacuum Bianchi II model [174] which begins with a fixed set of Kasner-like anisotropic expansion rates, and, possibly, makes one change of the rates in a prescribed way (Mixmaster-like bounce) and then continues to τ = ∞ as a fixed Kasner solution. In contrast are the homogeneous cosmologies, which display Mixmaster dynamics such as vacuum Bianchi VIII and IX [15, 139, 100] and Bianchi VI0 and Bianchi I with a magnetic field [131, 20, 130]. Jantzen [122] has discussed other examples. Mixmaster dynamics describes an approach to the singularity which is a sequence of Kasner epochs with a prescription, originally due to Belinskii, Khalatnikov, and Lifshitz (BKL) [15], for relating one Kasner epoch to the next. Some of the Mixmaster bounces (era changes) display sensitivity to initial conditions one usually associates with chaos, and in fact Mixmaster dynamics is chaotic [64]. The vacuum Bianchi I (Kasner) solution is distinguished from the other Bianchi types in that the spatial scalar curvature 3R, (proportional to) the minisuperspace (MSS) potential [139, 167], vanishes identically. But 3R arises in other Bianchi types due to spatial dependence of the metric in a coordinate basis. Thus an AVTD singularity is also characterized as a regime in which terms containing or arising from spatial derivatives no longer influence the dynamics. This means that the Mixmaster models do not have an AVTD singularity, since the influence of the spatial derivatives (through the MSS potential) never disappears — there is no last bounce.

3.2 Symplectic Numerical Methods

Symplectic numerical methods have proven useful in studies of the approach to the singularity in cosmological models [27]. Symplectic ODE and PDE [75, 141] methods comprise a type of operator splitting. An outline of the method (for one degree of freedom) follows. Details of the application to the Gowdy model (PDE’s in one space and one time direction) are given elsewhere [30].

For a field q(t) and its conjugate momentum p(t) split the Hamiltonian operator into kinetic and potential energy subhamiltonians. Thus,

$$H = \int\nolimits {dx} \{{1 \over 2}{p^2} + V(q)\} = {H_1}(p) + {H_2}(q).$$
(6)

if the vector X = (p, q) defines the variables at time t, then the time evolution is given by

$${{dX} \over {dt}} = {\{H,X\} _{PB}} \equiv AX$$
(7)

where { }PB is the Poisson bracket. The usual exponentiation yields an evolution operator

$${e^{A\Delta t}} = {e^{{A_1}}}^{(\Delta t/2)}{e^{{A_2}\Delta t}}{e^{{A_1}}}^{(\Delta t/2)} + O(\Delta {t^3})$$
(8)

for A = A1 + A2 the generator of the time evolution. Higher order accuracy may be obtained by a better approximation to the evolution operator [172, 173]. This method is useful when exact solutions for the subhamiltonians are known. For the given H, variation of H1 yields the solution

$$q = {q_0} + {p_0}\,\Delta t\quad, \quad p = {p_0},$$
(9)

while that of H2 yields

$$q = {q_0}\quad, \quad p = {p_0} - {\left. {{{dV} \over {dq}}} \right|_{{q_0}}}\Delta t \, \, .$$
(10)

Note that H2 is exactly solvable for any potential V no matter how complicated, although the required differenced form of the potential gradient may be non-trivial. One evolves from t to t + Δt using the exact solutions to the sub-hamiltonians according to the prescription given by the approximate evolution operator (8). Extension to more degrees of freedom and to fields is straightforward [30, 21].

3.3 Mixmaster Dynamics

3.3.1 Overview

Belinskii, Khalatnikov, and Lifshitz [15] (BKL) described the singularity approach of vacuum Bianchi IX cosmologies as an infinite sequence of Kasner [124] epochs whose indices change when the scalar curvature terms in Einstein’s equations become important. They were able to describe the dynamics approximately by a map evolving a discrete set of parameters from one Kasner epoch to the next [15, 57]. For example, the Kasner indices for the power law dependence of the anisotropic scale factors can be parametrized by a single variable u ≥ 1. BKL determined that

$${u_{n + 1}} = \left\{{\begin{array}{*{20}c} {{u_n} - 1\quad, \quad 2 \le {u_n}}\quad \quad \quad \\ {{{({u_n} - 1)}^{- 1}}\quad, \quad 1 \le {u_n} \le 2} \\ \end{array} \quad.} \right.$$
(11)

The subtraction in the denominator for 1 ≤ un ≤ 2 yields the sensitivity to initial conditions associated with Mixmaster dynamics (MD). Misner [139] described the same behavior in terms of the model’s volume and anisotropic shears. A multiple of the scalar curvature acts as an outward moving potential in the anisotropy plane. Kasner epochs become straight line trajectories moving outward along a potential corner while bouncing from one side to the other. A change of corner ends a BKL era when u → (u − 1)−1. Numerical evolution of Einstein’s equations was used to explore the accuracy of the BKL map as a descriptor of the dynamics as well as the implications of the map [145, 164, 166, 19]. Rendall has studied analytically the validity of the BKL map as an approximation to the true trajectories [162].

Later, the BKL sensitivity to initial conditions was discussed in the language of chaos [6, 125]. An extended application of Bernoulli shifts and Farey trees was given by Rugh [165] and repeated by Cornish and Levin [63]. However, the chaotic nature of Mixmaster dynamics was questioned when numerical evolution of the Mixmaster equations yielded zero Lyapunov exponents (LE’s) [76, 45, 111]. (The LE measures the divergence of initially nearby trajectories. Only an exponential divergence, characteristic of a chaotic system, will yield positive exponent.) Other numerical studies yielded positive LE [160]. This issue was resolved when the LE was shown numerically and analytically to depend on the choice of time variable [164, 18, 73]. Although MD itself is well-understood, its characterization as chaotic or not had been quite controversial [112].

LeBlanc et al [131] have shown (analytically and numerically) that MD can arise in Bianchi VI0 models with magnetic fields (see also [133]). In essence, the magnetic field provides the wall needed to close the potential in a way that yields the BKL map for u [20]. A similar study of magnetic Bianchi I has been given by LeBlanc [130]. Jantzen has discussed which vacuum and electomagnetic cosmologies could display MD [122].

3.3.2 Recent Developments

Recently, Cornish and Levin (CL) [64] identified a coordinate invarient way to characterize MD. Sensitivity to initial conditions can lead to qualitatively distinct outcomes from initially nearby trajectories. While the LE measures the exponential divergence of the trajectories, one can also “color code” the regions of initial data space corresponding to particular outcomes. A chaotic system will exhibit a fractal pattern in the colors. CL defined the following set of discrete outcomes: During a numerical evolution, the BKL parameter u is evaluated from the trajectories. The first time u > 7 appears in an approximately Kasner epoch, the trajectory is examined to see which metric scale factor has the largest time derivative. This defines three outcomes and thus three colors for initial data space.

To study the CL fractal and ergodic properties of Mixmaster evolution [64], one could take advantage of a new numerical algorithm due to Berger, Garfinkle, and Strasser [28]. Symplectic methods are used to allow the known exact solution for a single Mixmaster bounce [174] to be used in the ODE solver. Standard ODE solvers [159] can take large time steps in the Kasner segments but must slow down at the bounce. The new algorithm patches together exactly solved bounces. Tens of orders of magnitude improvement in speed are obtained while the accuracy of machine precision is maintained. In [28], the new algorithm was used to distinguish Bianchi IX and magnetic Bianchi VI0 bounces. This required an improvement of the BKL map (for parameters other than u) to take into account details of the exponential potential.

3.3.3 Going Further

One can easily invent prescriptions other than that given by Cornish and Levin [64] which would yield discrete outcomes. The fractal nature of initial data space should be common to all of them. It is not clear how the value of the fractal dimension as measured by Cornish and Levin would be affected.

There are also recent numerical studies of Mixmaster dynamics in other theories of gravity. For example, Carretero-Gonzalez et al [52] find evidence of chaotic behavior in Bianchi IX-Brans-Dicke solutions while Cotsakis et al [65] have shown that Bianchi IX models in 4th order gravity theories have stable non-chaotic solutions. Barrow and Levin find evidence of chaos in Bianchi IX Einstein-Yang-Mills (EYM) cosmologies [7]. Their analysis may be applicable to the corresponding EYM black hole interior solutions.

Finally, we remark on a successful application of numerical Regge calculus in 3 + 1 dimensions. Gentle and Miller have been able to evolve the Kasner solution [84].

3.4 Inhomogeneous Cosmologies

3.4.1 Overview

BKL have conjectured that one should expect a generic singularity in spatially inhomogeneous cosmologies to be locally of the Mixmaster type [15]. For a review of homogeneous cosmologies, inhomogeneous cosmologies, and the relation between them, see [134]. The main difficulty with the acceptance of this conjecture has been the controversy over whether the required time slicing can be constructed globally [8, 89]. Montani [144], Belinskii [10], and Kirillov and Kochnev [128, 127] have pointed out that if the BKL conjecture is correct, the spatial structure of the singularity could become extremely complicated as bounces occur at different locations at different times. A class of cosmological models which might have local MD are vacuum universes on T3 × R with a U(1) symmetry [142]. The non-commuting Killing vectors of local MD can be constructed since only one Killing vector is already present. The two commuting Killing vectors of the even simpler plane symmetric Gowdy cosmologies [88, 17] preclude their use to test the conjecture. However, these models are interesting in their own right since they have been conjectured to possess an AVTD singularity [90].

Polarized plane symmetric cosmologies have been evolved numerically using standard techniques by Anninos, Centrella, and Matzner [2, 3]. A long-term project involving Berger, Garfinkle, and Moncrief and our students to study the generic cosmological singularity numerically has been reviewed in detail elsewhere [27] but will be discussed briefly here.

3.4.2 Gowdy Cosmologies

The Gowdy model on T3 × R is described by gravitational wave amplitudes P(θ, τ) and Q(θ, τ) which propagate in a spatially inhomogeneous background universe described by λ(θ, τ). (We note that the physical behavior of a Gowdy spacetime can be computed from the effect of the metric evolution on a test cylinder [29].) We impose 0 ≤ 6 ≤ 2π and periodic boundary conditions. The time variable τ measures the area in the symmetry plane with τ = ∞ a curvature singularity.

Einstein’s equations split into two groups. The first is nonlinearly coupled wave equations for dynamical variables P and Q (where, a = ∂/∂a) obtained from the variation of [140]

$$\begin{array}{*{20}c} {H \quad = \quad {1 \over 2}\int\limits_0^{2\pi} {d\theta} [\pi _P^2 + {e^{- 2P}}\pi _Q^2]\quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ \quad {1 \over 2}\int\limits_0^{2\pi} {d\theta} [{e^{- 2\tau}}(P,_\theta ^2 + {e^{2P}}Q,_{\theta}^2)] = {H_1} + {H_2}.} \\ \end{array}$$
(12)

This Hamiltonian has the form required by the symplectic scheme. If the model is, in fact, AVTD, the approximation in the symplectic numerical scheme should become more accurate as the singularity is approached. The second group of Einstein equations contains the Hamiltonian and θ-momentum constraints respectively. These can be expressed as first order equations for λ in terms of P and Q. This break into dynamical and constraint equations removes two of the most problematical areas of numerical relativity from this model — the initial value problem and numerical preservation of the constraints.

For the special case of the polarized Gowdy model (Q = 0), P satisfies a linear wave equation whose exact solution is well-known [17]. For this case, it has been proven that the singularity is AVTD [119]. This has also been conjectured to be true for generic Gowdy models [90].

AVTD behavior is defined in [119] as follows: Solve the Gowdy wave equations neglecting all terms containing spatial derivatives. This yields the AVTD solution [30]. If the approach to the singularity is AVTD, the full solution comes arbitrarily close to an AVTD solution at each spatial point as τ → ∞. As τ → ∞, the AVTD solution becomes

$$P(\theta, \, \tau) \to v (\theta)\tau \quad, \quad Q(\theta, \tau) \to {Q_0}(\theta)$$
(13)

where υ > 0. Substitution in the wave equations shows that this behavior is consistent with asymptotic exponential decay of all terms containing spatial derivatives only if 0 ≤ υ ≤ 1 [90]. We have shown that, except at isolated spatial points, the nonlinear terms in the wave equation for P drive υ into this range [25, 27]. The exceptional points occur when coefficients of the nonlinear terms vanish and are responsible for the growth of spiky features seen in the wave forms [30, 25]. We conclude that generic Gowdy cosmologies have an AVTD singularity except at isolated spatial points [25, 27]. A claim to the contrary [107] is incorrect [26]. Recently, it has been proven analytically that Gowdy solutions with 0 < υ < 1 and AVTD behavior almost everywhere are generic [126].

Addition of a magnetic field to the vacuum Gowdy models (plus a topology change) which yields the inhomogeneous generalization of magnetic Bianchi VI0 models provides an additional potential which grows exponentially if 0 < υ < 1. Local Mixmaster behavior has recently been observed in these magnetic Gowdy models [181].

3.4.3 U(1) Symmetric Cosmologies

Given the success of the symplectic method in studying the singularity behavior of the Gowdy model, we can consider its extension to the case of U(1) symmetric cosmologies. Moncrief has shown [142] that cosmological models on T3 × R with a spatial U(1) symmetry can be described by five degrees of freedom {x, z, Λ, φ, ω} and their respective conjugate momenta {px, pz, pΛ, p, r} All variables are functions of spatial variables u, v and time, τ. Einstein’s equations can be obtained by variation of

$$\begin{array}{*{20}c} {H\quad = \quad \oint {\oint {du\,dv\,{\mathcal H}}} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \,\,} \\ {\quad = \quad \oint {\oint {du\,dv\,\left({{1 \over 8}p_z^2 + {1 \over 2}{e^{4z}}p_x^2 + {1 \over 8}{p^2} + {1 \over 2}{e^{4\varphi}}{r^2} - {1 \over 2}p_\Lambda ^2 + 2{p_\Lambda}} \right)}}} \\ {+ {e^{- 2\tau}} = \oint {\oint {du\,dv} \left\{{\left({{e^\Lambda}{e^{ab}}} \right){,_{ab}} - \left({{e^\Lambda}{e^{ab}}} \right){,_a}\Lambda {,_b}\quad \quad} \right.} \quad} \\ {+ {e^\Lambda}\left[ {\left({{e^{- 2z}}} \right){,_u}x{,_v} - \left({{e^{- 2z}}} \right){,_v}x{,_u}} \right]\quad \quad \quad \quad \quad \quad \quad \quad} \\ {\left. {+ 2{e^\Lambda}{e^{ab}}\varphi {,_a}\varphi {,_b} + {1 \over 2}{e^\Lambda}{e^{- 4\varphi}}{e^{ab}}\omega {,_a}\omega {,_b}} \right\}\quad \quad \quad \quad \quad \quad} \\ {= \quad {H_1} + {H_2}.\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$
(14)

Here φ and ω are analogous to P and Q while eΛ is a conformal factor for the metric eab(x, z) in the u-v plane perpendicular to the symmetry direction. Note particularly that H1 contains two copies of the Gowdy H1 plus a free particle term and is thus exactly solvable. The potential term H2 is very complicated. However, it still contains no momenta, so its equations are trivially exactly solvable. However, issues of spatial differencing are problematic. (Currently, a scheme due to Norton [150] is used. A spectral evaluation of derivatives [74] which has been shown to work in Gowdy simulations [16] does not appear to be helpful in the U(1) case.) A particular solution to the initial value problem is used, since the general solution is not available [27]. Currently, except as discussed below, the constraints are monitored but not explicitly enforced.

Despite the current limitations of the U(1) code, conclusions can be confidently drawn for polarized U(1) models in our restricted class of initial data. No numerical difficulties arise in this case. Polarized models have r = ω = 0. Grubišić and Moncrief [91] have conjectured that these polarized models are AVTD. The numerical simulations provide strong support for this conjecture [27, 32].

3.4.4 Going Further

Methods similar to those used in [91], indicate that the term containing gradients of ω in (14) acts as a Mixmaster-like potential to drive the system away from AVTD behavior in generic U(1) models [23]. Numerical simulations provide support for this suggestion [27, 31]. Whether this potential term grows or decays depends on a function of the field momenta. This in turn is restricted by the Hamiltonian constraint. However, failure to enforce the constraints can cause an erroneous relationship among the momenta to yield qualitatively wrong behavior. There is numerical evidence that this error tends to suppress Mixmaster-like behavior leading to apparent AVTD behavior in extended spatial regions [21, 22]. In fact, it has been found very recently [31] that when the Hamiltonian constraint is enforced at every time step, the predicted local oscillatory behavior of the approach to the singularity is observed.

Mixmaster simulations with the new algorithm [28] can easily evolve more than 250 bounces reaching |Ω| ≈ 1062. This compares to earlier simulations yielding 30 or so bounces with |Ω| ≈ 108. The larger number of bounces quickly reveals that it is necessary to enforce the Hamiltonian constraint. A similar conclusion for gravitational wave equations was obtained by Gundlach and Pullin [99]. Although this particular analysis may be incorrect [46], it is still likely that constraint enforcement will be essential for sufficiently long simulations. An explicitly constraint enforcing U(1) code was developed some years ago by Ove (see [153] and references therein).

4 Discussion

Numerical investigation of singularities provides an arena for the close coupling of analytic and numerical techniques. The references provided here contain many examples of analytic results guided by numerical results and numerical studies to demonstrate whether or not analytic conjectures are valid.

Even more striking is the convergence of the separate topics of this review. While the search for naked singularities in the collapse of highly prolate systems has yielded controversial results, a naked singularity was discovered in the collapse of spherically symmetric scalar fields. The numerical exploration of cosmological singularities has yielded strong evidence that the asymptotic behavior is local — each spatial point evolves toward the singularity as a separate universe. This means that conclusions from these studies should be relevant in any generic collapse. This area of research then should begin to overlap with the studies of black hole interiors (see for example [48]).