1 Introduction

The singularity theorems [244, 139, 140, 141] state that Einstein’s equations will not evolve generic, 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–Nordström 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 [95, 164] 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) [22, 187] 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 [97, 240].

Once the simplest, exactly solvable models are left behind, understanding of the singularity becomes more difficult. There has been significant analytic progress [245, 191, 219, 3]. However, until recently such methods have yielded 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 [100]. 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 [239], such configurations may yield a naked singularity. Analytic results suggest that one must go beyond the failure to observe an apparent horizon to conclude that a naked singularity has formed [245]. 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 [77]. In the initial study, the critical (Choptuik) solution is a zero mass naked singularity (visible from null infinity). It is a counterexample to the cosmic censorship conjecture [135]. However, it is a non-generic one since fine-tuning of the initial data is required to produce this critical solution. In a possibly related study, Christodoulou has shown [81] that for the spherically symmetric Einstein–scalar field equations, there always exists a perturbation that will convert a solution with a naked singularity (but of a different class from Choptuik’s) to one with a black hole. Reviews of critical phenomena in gravitational collapse can be found in [46, 126, 129, 131].

The second question which is now beginning to yield to numerical attack involves the stability of the Cauchy horizon in charged or rotating black holes. It has been conjectured [244, 73] 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 [240]. Numerical studies [56, 92, 63] show that a weak, null singularity forms first as had been predicted [212, 202].

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 [193, 223, 225, 25, 102, 62, 147, 216]. A coordinate invariant characterization of Mixmaster chaos has been formulated [85] which, while criticized in its details [194], has essentially resolved the question. In addition, a new extremely fast and accurate algorithm for Mixmaster simulations has been developed [39].

Belinskii, Khalatnikov, and Lifshitz (BKL) long ago claimed [17, 18, 19, 22, 21] that it is possible to formulate the generic cosmological solution to Einstein’s equations near the singularity as a Mixmaster universe at every spatial point. While others have questioned the validity of this claim [13], numerical evidence has been obtained for oscillatory behavior in the approach to the singularity of spatially inhomogeneous cosmologies [250, 43, 37, 41]. We shall discuss results from a numerical program to address this issue [42, 37, 31]. The key claim by BKL is that sufficiently close to the singularity, each spatial point evolves as a separate universe – most generally of the Mixmaster type. For this to be correct, the dynamical influence of spatial derivatives (embodying communication between spatial points) must be overwhelmed by the time dependence of the local dynamics. In the past few years, numerical simulations of collapsing, spatially inhomogeneous cosmological spacetimes have provided strong support for the BKL picture [42, 36, 44, 250, 43, 37, 41]. In addition, the Method of Consistent Potentials (MCP) [123, 37] has been developed to explain how the BKL asymptotic state arises during collapse. New asymptotic methods have been used to prove that open sets exist with BKL’s local behavior (although these are AVTD rather than of the Mixmaster type) [163, 173, 3]. Recently, van Elst, Uggla, and Wainwright developed a dynamical systems approach to cosmologies (i.e. those with 2 spatial symmetries) [242].