In Chapter 2 we performed a 3 + 1 decomposition of Einstein's field equations and have seen that these can be split into two distinct sets: constraint equations and evolution equations. The constraint equations contain no time derivatives and relate field quantities on a given t = constant spacelike hypersurface. The evolution equations contain first-order time derivatives that tell us how the field quantities change from one hypersurface to the next. In Chapter 3 we have brought the constraint equations into a form that is suitable for numerical implementation, that is, we cast the equations in terms of spatial differential operators that can be inverted with standard numerical techniques. We will provide a brief introduction to some common numerical algorithms for solving these (elliptic) equations in Chapter 6. The 3 + 1 evolution equations that we derived, e.g., equation (2.134) for γij, and equation (2.135) for Kij, are not quite ready for numerical integration. For one thing, we have yet to impose coordinate conditions by specifying the lapse function α and the shift vector βithat appear in these equations. The lapse and shift are freely specifiable gauge variables that need to be chosen in order to advance the field data from one time slice to the next. As it turns out, finding kinematical conditions for the coordinates that allow for a well-behaved, long time evolution is nontrivial in general. However, geometric insight and numerical experimentation can be combined to produce good gauge choices for treating many of the most important physical and astrophysical problems requiring numerical relativity for solution, as we shall see.

This brings us to the end of our introduction to numerical relativity. A quick glance at the table of contents shows that we have covered a wide range of subjects, starting with the foundations of numerical relativity and continuing with applications to different areas of gravitational physics and astrophysics. Despite the breadth of our survey, we have had to be selective in our choice of topics. Our focus has been on solving the Cauchy problem in general relativity for dynamical, asymptotically flat spacetimes, with applications to compact objects and compact binaries. There are a number of alternative approaches for solving Einstein's equations that we did not touch on at all, such as the characteristic approach and the Regge calculus. We also did not discuss in any detail applications involving strictly stationary spacetimes, such as gas accretion onto Kerr black holes, although some of the same schemes we described for matter evolution have been used successfully to treat problems with fixed background metrics.We trust that interested readers will find discussions of the subjects we omitted elsewhere in the literature.

We hope that our treatment laying out the foundations of numerical relativity will remain relevant for the foreseeable future. However, we suspect that some of the large-scale simulations we have chosen to illustrate different implementations will be superseded by more sophisticated calculations.

Black holes are characterized by the horizons surrounding them. Clearly, then, the numerical simulation of black holes requires the ability to locate and analyze black hole horizons in numerically generated spacetimes. In this chapter we first review different concepts of horizons in asymptotically flat spacetimes, and then discuss how these horizons can be probed numerically.

Concepts

Several different notions of horizons exist in general relativity. The defining property of a black hole is the presence of an event horizon (Section 7.2), but, as we will see, apparent horizons (Section 7.3) also play an extremely important role in the context of numerical relativity. In addition, the concepts of isolated and dynamical horizons (Section 7.4) serve as useful diagnostics in numerical spacetimes containing black holes.

A black hole is defined as a region of spacetime from which no null geodesic can escape to infinity. The surface of a black hole, the event horizon, acts as a one-way membrane through which light and matter can enter the black hole, but once inside, can never escape. It is the boundary in spacetime separating those events that can emit light rays that can propagate to infinity and those which cannot. More precisely, the event horizon is defined as the boundary of the causal past of future null infinity. It is a 2 + 1 dimensional hypersurface in spacetime formed by those outward-going, future-directed null geodesics that neither escape to infinity nor fall toward the center of the black hole.

In this chapter we assemble some of the elements of Einstein's theory of general relativity that we will be working with in later chapters. We assume that the geometric objects and equations that we list, as well as their interpretation, are already very familiar to readers. The discussion below should serve simply as a checklist of a few of the basics that we need to pack with us before embarking on our voyage into numerical spacetime.

Throughout this book we adopt the (−+++) metric signature together with all the sign conventions of Misner et al. (1973). Following that book, but in this chapter only, we will display a tensor in spacetime by a symbol in boldface when emphasizing its coordinatefree character, or by its components when the tensor has been expanded in a particular set of basis tensors. However, unlike that book, we will use Latin indices a, b, … instead of Greek letters to denote the spacetime indices of the tensor components, with the values of the indices running from 0 to 3. This choice anticipates a switch we will make to abstract index notation in all subsequent chapters of this book. We will introduce this switch in Section 2.1. We adopt the usual Einstein convention of summing over repeated indices. Finally, here and throughout we will use geometrized units in which both the gravitational constant and the speed of light are assigned the values of one, G = c = 1.

Rotating stars in general relativity are of long-standing theoretical interest. With the discovery of radio pulsars in 1967 and their identification as magnetized, rotating neutron stars, rotating relativistic stars have also become objects of intense observational scrutiny. Radio pulsars have become particularly useful for relativity, since they provide excellent cosmic clocks. The original idea of Baade and Zwicky (1934) that neutron stars could be formed during supernovae events, coupled with the discovery of radio pulsars in some supernova remnants like the Crab, firmly linked rotating neutron stars with stellar collapse and supernovae explosions. Following the detection of radio pulsars, X-ray pulsars were discovered in 1971. X-ray pulsars proved the existence of rotating neutron stars that are accreting gas supplied by a companion star in a binary system. The discovery of the first binary pulsar by Hulse and Taylor in 1974 proved that rotating, relativistic stars can reside in binary systems, making neutron stars even more interesting to relativists. General relativity is required to describe the gravitational field of neutron stars, since their compaction, M/R, where M is the mass and R is the characteristic radius of the star, is large (˜0.1–0.2). Numerical relativity is important for treating many dynamical processes involving neutron stars accurately. It provides a valuable tool for probing neutron star formation from stellar core collapse in a supernova and for tracking the dynamical evolution and assessing the final fate of neutron stars subject to instabilities.

Turn now from the most general spacetime in full 3 + 1 dimensions to the special case of spherical symmetry. Why should we do this? Actually, spherical systems provide very useful computational, physical and astrophysical insight and working with them serves multiple purposes. The field equations reduce to 1 + 1 dimensions – variables may be written as functions of only two parameters, a time coordinate t and a suitable radial coordinate r – and are much simpler to solve in spherical symmetry. Solving them is a very cost-effective way of probing dynamical spacetimes with strong gravitational fields, including spacetimes with black holes. After all, nonrotating stars and black holes are themselves spherical, so many important aspects of gravitational collapse, including black hole formation and growth, can be studied in spherical symmetry. For example, the numerical study of spherically symmetric collapse to black holes led to the discovery of critical phenomena in black hole formation. The simplification in the equations, together with the reduction in the number of spatial dimensions, means that the system of spherical equations can be solved more quickly, in terms of both human input and computer time, and with much higher accuracy, than the set required for more general spacetimes. As a result, tackling problems in spherical symmetry provides an excellent starting point for learning how to do numerical relativity. It also serves as a convenient laboratory for experimenting with different gauge choices (coordinates) and for generating high precision, test-bed solutions for numerical codes designed to work in higher dimensions.

As we have seen, Einstein's field equations in 3 + 1 form consist of a set of nonlinear, multidimensional, coupled partial differential equations in space and time. The equations of motion of the matter fields that may be present are typically of a similar nature. Except for very idealized problems with special symmetries, such equations must be solved by numerical means, often on supercomputers. Just as there is no unique analytic formulation of the 3 + 1 field equations, there is no unique prescription by which a partial differential equation may be cast into a form suitable for numerical integration. Standard numerical algorithms for treating such equations may be found in many textbooks on numerical methods, as well as in textbooks, monographs and review articles on compuational physics. This branch of applied mathematics is a rich area of ongoing investigation; it progresses with each advance in computer technology. It would take us too far a field to review the subject in any depth here. Instead, we shall present a brief introduction to some of the basic numerical concepts and associated techniques, focusing on those most often employed to solve the partial differential equations that arise in numerical relativity. Although our treatment is rudimentary, we hope that it is sufficient to convey the flavor of the subject, especially to readers unfamiliar with the basic ideas. Throughout our discussion we shall refer the reader to some of the literature where further details and other references can be found.

As we learned in Chapter 8, where we studied spherical systems, collisionless clusters provide a simple relativistic source for exploring the nature of Einstein's equations and experimenting with numerical techniques to solve them. Once we relax the restriction to spherical symmetry, the spacetimes can exhibit two new dynamical features: rotation and gravitational waves. Not much is known about nonspherical collisionless systems in general relativity, even for stationary equilibria. Some interesting results have emerged by exploiting numerical relativity to investigate the equilibrium structure and collapse of nonspherical rotating and nonrotating clusters in axisymmetry.1 To highlight the power of the technique, we shall summarize a few of the simulations and their key findings in this chapter.

The examples discussed below are chosen to demonstrate how numerical relativity, quite apart from providing accurate quantitative solutions to dynamical scenarios involving strong gravitational fields, can provide qualitative insight into Einstein's equations in those cases where uncertainty still prevails. It can even be helpful as a guide to proving (or disproving) theorems about strong-field spacetimes in those instances where analytic means alone have not been adequate.

Collapse of prolate spheroids to spindle singularities

It is well-known that classical general relativity admits solutions with singularities, and that such solutions can be produced by the gravitational collapse of nonsingular, asymptotically flat initial data. The Cosmic Censorship Conjecture of Penrose states that such singularities will always be clothed by event horizons and hence can never be visible from the outside (no naked singularities).

The dynamical simulation of the head-on collision of two black holes in axisymmetry was an early success of numerical relativity (see Chapter 10.2). Based on this success one might surmise that the subsequent simulation of the inspiral and merger of binary black holes initially in circular orbit represented a straightforward generalization of the head-on case. It turned out, however, that relaxing the assumption of axisymmetry to treat a binary in circular orbit, and then tracking the resulting evolution, presented several nontrivial challenges. As a consequence, dynamical simulations of these binaries were stalled for many years until these challenges were finally overcome. Today, the inspiral and merger of binary black holes is essentially a solved problem, constituting one of the major triumphs of numerial relativity.

One obvious complication that arises in numerical simulations when moving from two to three spatial dimensions is the burden of increased computational resources required to cover the added dimension. Even though computers have become considerably faster and can handle far more memory than the machines available when the first head-on black hole simulations were performed, the resources required to resolve inspiraling black holes in the strong-field, near-zone while simultaneously tracking and ultimately extracting gravitational waves in the weak-field, far-zone remain formidable. Different investigators have adopted different approaches to address this problem of “dynamic range”, including the use of fixed or adaptive mesh refinement or the construction of novel coordinate systems that allocate gridpoints where they are most needed.

Our book is intended both as a general reference for researchers and as a textbook for use in a formal course that treats numerical relativity. We envision that there are at least two different ways in which the book can be used in the classroom: as the main text for a one-semester course on numerical relativity for students who have already taken an introductory course in general relativity, or as supplementary reading in numerical relativity at the end of an introductory course in general relativity. There may be more material in the book than can be covered comfortably in a single semester devoted entirely to numerical relativity. There certainly is more material than can be integrated into a supplementary unit on numerical relativity in an introductory course on general relativity.

The latter may be true even when such a course is taught as a two-semester sequence, if the course is already broad and comprehensive without numerical relativity. There are several ways to design a shortened presentation of the material in our book without sacrificing the core concepts or interfering with the logical flow. The amount of material that must be cut out from any course depends, naturally, on the amount of time that is available to devote to the subject. One means of reducing the content while retaining the fundamental ideas in a self-contained format is to restrict the discussion to pure vacuum spacetimes, i.e., spacetimes with no matter sources.

Simulating the late inspiral, merger and ringdown of a binary black hole was for many years the “holy grail” of numerical relativity. Dozens of researchers have spent many years attempting to formulate a stable algorithm capable of solving this important problem. As discussed in Chapter 13, several different techniques and corresponding codes now exist that can evolve coalescing binary black holes successfully. In this appendix we summarize two simulations in sufficient detail to enable the reader to assess the requirements and caliber of a typical computation. We employ one widely adopted method for the simulations, based on the BSSN evolution scheme and binary puncture initial data. The examples chosen here, and the quoted results, can be used as test problems for the construction of new algorithms and new modules.

We emphasize that performing a simulation of binary black hole coalescence provides a highly nontrivial and comprehensive laboratory for testing a 3 + 1 code capable of solving Einstein's field equations in vacuum in the strong-field regime of general relativity. It tests the caliber of not only the basic evolution scheme, but also the all-important diagnostic routines. These routines measure globally conserved quantities like the total mass and angular momentum of the system, the masses and spins of the individual black holes, the location and area of the black hole horizons, the asymptotic gravitational waveforms, the recoil kick imparted to the black hole remnant, and other “vital statistics” characterizing a dynamical spacetime containing black holes.

At this point we might suspect that we are all set to carry out dynamical simulations involving relativistic gravitational fields in three spatial dimensions, i.e., dynamical simulations in full 3 + 1 dimensional spacetimes. After all, we have derived the 3 + 1 evolution equations for the gravitational fields in Chapter 2, developed techniques for the construction of initial field data in Chapter 3, and discussed strategies for imposing suitable coordinate conditions in Chapter 4. Should the spacetime in which we are interested contain the most common matter sources, we can consult Chapter 5 for the matter source terms and the matter equations of motion. Beyond assembling all of the relevant equations, we have also sketched algorithms for solving them numerically in Chapter 6. If black holes are present, we have derived methods to locate and measure their horizons in Chapter 7; if gravitational waves are generated we have discussed how to extract radiation wave forms numerically in Chapter 9. Finally, we have evolved relativistic systems in spherical symmetry (1 + 1 dimensions) in Chapter 8 and in axisymmetry (2 + 1 dimensions) in Chapter 10. What else is left to do before going on to perform simulations in full 3 + 1 dimensions?

Suppose then we were to plunge ahead with the tools currently at our disposal to build a dynamical spacetime in 3 + 1 dimensions. We could gain some confidence and computational experience by choosing to explore as our first case a simple dynamical system with a known analytic solution.

Here we list the key equations describing a mean-field, particle simulation scheme that can treat the evolution of collisionless matter in axisymmetry according to general relativity. The scheme is a generalization of the one described in Chapter 8.2 for spherical systems and employs the standard ADM form of the field equations as listed in Box 2.1. We adopt spherical polar spacetime coordinates (t, r, θ,ϕ), assume axisymmetry and specialize to the case where there is no net angular momentum. In axisymmetry all quantities are functions only of (t, r, θ). We also impose maximal slicing and quasi-isotropic spatial coordinates as our gauge conditions. This spatial gauge condition reduces to isotropic coordinates for Schwarzschild geometry. The field equations listed below constitute a fully constrained approach to solving the Einstein field equations for this problem, i.e., one which solves all of the constraint equations in lieu of integrating evolution equations for some of the variables. Fully constraint schemes have the advantage over unconstrained schemes that the constraints are guaranteed to be satisfied at all times, which may in some cases also eliminate some instabilities associated with the evolution equations. Their disadvantage is that the constraints constitute elliptic equations, which typically require more computational resources to solve than explicit time evolution equations. This disadvantage is not so severe, however, in 1 + 1 or 2 + 1 spacetimes. A similar set of variables and field equations to the ones summarized below has been used to simulate the gravitational collapse of hydrodynamic fluids and vacuum gravitational waves6 in nonrotating, axisymmetric spacetimes, as well as the head-on collision of neutron stars.