In this chapter we wish to construct initial data for quasiequilibrium binary black holes. That is, we seek solutions corresponding to two black holes in stable, nearly circular orbit about each other. In contrast to Newtonian theory, a stellar binary in general relativity can never be in strict equilibrium, with the companions moving in exactly circular orbits at constant separation for all time. Instead, gravitational radiation emission inevitably leads to loss of orbital energy and angular momentum, causing the orbit to decay. The resulting trajectory then traces out an inspiral rather than a perfect circle. For sufficiently large separations, the binary motion is nearly Newtonian, hence the orbit is nearly circular, decaying very little during one orbital period. If isolated from outside perturbations (e.g., gravitational encounters with other stars), it is expected that astrophysical binaries composed of compact stars (i.e., compact binaries) will ultimately evolve to a quasiequilibrium state following their formation at large separation. The reason is that gravitational radiation loss drives orbital circularization as well as decay, as we will discuss in the next section. Only when the orbits become very close and highly relativistic, just prior to radial plunge and binary merger, do the deviations from circular motion become large.

The construction of quasiequilibrium binary initial data poses a number of conceptual challenges. Getting started, however, is fairly straightforward. To find solutions we shall follow the approaches outlined in Chapter 3.

What is numerical relativity?

General relativity – Einstein's theory of relativistic gravitation – is the cornerstone of modern cosmology, the physics of neutron stars and black holes, the generation of gravitational radiation, and countless other cosmic phenomena in which strong-field gravitation plays a dominant role. Yet the theory remains largely untested, except in the weak-field, slow-velocity regime. Moreover, solutions to Einstein's equations, except for a few idealized cases characterized by high degrees of symmetry, have not been obtained as yet for many of the important dynamical scenarios thought to occur in nature. With the advent of supercomputers, it is now possible to tackle these complicated equations numerically and explore these scenarios in detail. That is the main goal of numerical relativity, the art and science of developing computer algorithms to solve Einstein's equations for astrophysically realistic, high-velocity, strong-field systems.

Numerical relativity has become one of the most powerful probes of relativistic spacetimes. It is the tool that allows us to recreate cataclysmic cosmic phenomena that are otherwise inaccessible in the conventional laboratory – like gravitational collapse to black holes and neutron stars, the inspiral and coalescence of binary black holes and neutron stars, and the generation and propagation of gravitational waves, to name a few. Numerical relativity picks up where post-Newtonian theory and general relativistic perturbation theory leave off. It enables us to follow the full nonlinear growth of relativistic instabilities and determine the final fate of unstable systems.

Binary black hole–neutron stars have received significantly less attention than binary black holes or binary neutron stars. No black hole–neutron star binary has been identified to date. However, stellar population synthesis models suggest that such systems represent a significant fraction of all compact binary mergers ultimately visible in gravitational waves by the LIGO detector. In addition, the study of black hole–neutron star mergers is important in light of the localizations of short-hard gamma-ray bursts. These GRB sources are found in galactic regions of low star-formation devoid of supernovae associations, ruling out massive stars as progenitors: massive stars have very short lifetimes and would need to be replenished more rapidly than is possible in low star-formation regions to account for these bursts. A more plausible progenitor for a short-hard GRB is a compact binary containing a neutron star, i.e., either a binary neutron star or binary black holeneutron star. The short-hard burst time scales and energetics are consistent with GRB models based on the coalescence of such compact binaries, and the evolution time scale of over 1 Gyr between formation and merger is consistent with the low star-formation rate.

Black hole-neutron star binaries can merge in two distinct ways. The neutron star may either be tidally disrupted by the black hole companion before being consumed, or it may be swallowed by the black hole more or less intact.

Spherically symmetric spacetimes, which we discussed in Chapter 8, do not admit gravitational radiation. Once we relax this symmetry restriction, as we shall do in the following chapters, we will encounter spacetimes that do contain gravitational radiation. In fact, simulating promising sources of gravitational radiation and predicting their gravitational wave signals are among the most important goals of numerical relativity. These goals are especially urgent in light of the new generation of gravitational wave laser interferometers which are now operational. A book on numerical relativity therefore would not be complete without a discussion of gravitational waves.

In this chapter we review several topics related to gravitational waves. We start in Section 9.1 with a discussion of linearized waves propagating in nearly Minkowski spacetimes and the role that these waves play even in the case of nonlinear sources of gravitational radiation. In Section 9.2 we survey plausible sources of gravitational waves, highlighting those that seem most promising from the perspective of gravitational wave detection. We briefly describe some of the existing and planned gravitational wave detectors in Section 9.3. Finally, in Section 9.4 we make contact with numerical relativity, and review different strategies that have been employed to extract gravitational radiation data from numerical relativity simulations.

Linearized waves

Most of this book deals with strong-field solutions of Einstein's equations, including black holes, neutron stars, and binaries containing these objects. As long as these solutions are dynamical and nonspherical, they will emit gravitational radiation.

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.