The numerical relativity breakthrough for binary black holes

The interaction of two point masses is probably the simplest and most fundamental dynamical problem one can conceive of in a theory of gravity. This problem is sufficiently simple in Newton's theory such that it can be solved analytically. In spite of the simplifications, the Newtonian two-body problem describes with high accuracy a wide class of physical systems, ranging from the planetary orbits in the solar system to the motion of the spacecraft that carried humans to the moon in 1969. Observed deviations in the motion of Uranus from the Newtonian predictions led to the prediction of a further planet, Neptune, that was indeed identified in 1846. For a while, a similar explanation was considered for anomalies observed in the perihelion precession of Mercury. The conjectured planet 'Vulcan', however, has never been found and in this case the explanation came in the form of a modified theory of gravity, namely Einstein's general relativity (GR).

GR differs from Newtonian gravity not only in terms of quantitative predictions but also presents a conceptually totally different description of gravity. Acceleration of objects due to gravitational interaction is no longer the result of a force but due to the curvature of spacetime itself. Mathematically, the spacetime is described in terms of a manifold $\mathcal{M}$ equipped with a metric ${{g}_{\alpha \beta }}$ which is determined through Einstein's field equations

$$\begin{eqnarray}&&{{R}_{\alpha \beta }}-\frac{1}{2}{{g}_{\alpha \beta }}R+\Lambda {{g}_{\alpha \beta }}=\frac{8\pi G}{{{c}^{4}}}{{T}_{\alpha \beta }}.\end{eqnarray} \tag{ 1 }$$

Here, Greek indices range from 0 to 3, ${{R}_{\alpha \beta }}$ is the Ricci tensor, R the Ricci scalar, Λ the cosmological constant and ${{T}_{\alpha \beta }}$ the energy momentum tensor. Unless specified otherwise, we will work in units where the gravitational constant and speed of light are unity, $G=1=c$. Much of this work will focus on the case of vacuum and asymptotically flat spacetimes where $\Lambda =0$ and ${{T}_{\alpha \beta }}=0$ and the Einstein equations become ${{R}_{\alpha \beta }}=0$.

The Ricci tensor ${{R}_{\alpha \beta }}$ is a nonlinear function of the spacetime metric components ${{g}_{\alpha \beta }}$ and their first and second derivatives. The Einstein equations couple space and time in a complex manner to the gravitational sources which is encapsulated in Wheeler's popular phrase 'Matter tells spacetime how to curve, spacetime tells matter how to move'. As we shall discuss in more detail in section 3, this nonlinear coupling of geometry and sources¹ makes the two-body problem much more complicated in GR but also lends a vast richness of new physics to these seemingly simple systems. Most importantly, we have the following conceptual differences between the Newtonian and the general relativistic case. (i) Sources of finite mass-energy cannot be point-like in GR but inevitably represent extended regions of non-vanishing curvature. The closest approximation to a point-like source in GR is a black hole (BH) and the two-body problem in GR therefore is a BH binary. (ii) The interaction of the two BHs is dissipative as energy and momentum can be radiated away from the binary in the form of gravitational waves (GW) which are subject of large-scale efforts for direct detection. Bound systems therefore eventually result in the merger of the two constituents.

In view of these special features of GR, the BH binary problem is often regarded as a three-stage process: (i) an extended phase of the interaction of two separate BHs, often referred to as the inspiral for bound systems, (ii) the merger, and (iii) the ringdown, a process of damped sinusoidal oscillations as the post-merger remnant sheds all structure beyond mass and angular momentum and settles down into a stationary Kerr BH. Unbound systems do not undergo stage (ii) and (iii) of this process but may still interact in a highly nonlinear manner during stage (i). The challenge to accurately model all possible stages of the dynamics of a BH binary then consists in solving Einstein's vacuum equations ${{R}_{\alpha \beta }}=0$, a system of 10 coupled, nonlinear second-order partial differential equations (PDEs). This challenge has often been referred to as the Holy Grail of numerical relativity (NR) and how it has eventually been met is the subject of this review.

For understanding the magnitude of this challenge and the particular issues arising in GR, it will be instructive to contrast it with its Newtonian counterpart and we will therefore start in section 2 with a brief review of the Newtonian two-body problem. The GR case is then summarized in section 3 including a 'to do list' of items that specifically arise in solving Einstein's rather than Newton's equations. The methodology to address these items is discussed in section 4. We continue in section 5 with an overview of the historical progress of the community which culminated in the 2005 breakthroughs by Pretorius [1] as well as the moving puncture method of the Brownsville (now Rochester) [2] and the NASA Goddard [3] groups who, quite remarkably, presented two rather different methods to solve the BH binary problem within some months. In section 6 we summarize the physical features of the dynamics of BH binary systems and briefly discuss the importance of the GR two-body problem in contemporary physics. We conclude in section 7 where we also provide references for further reading.

In the Newtonian two-body problem, we consider two point masses m₁ and m₂ moving in a background space and time. The two masses are separated by a distance vector $\vec{r}\equiv {{\vec{r}}_{1}}-{{\vec{r}}_{2}}$, and we denote by $\vec{F}$ the gravitational force exerted by m₂ on m₁. By Newton's laws m₁ acts on m₂ with $-\vec{F}$ and we have

$$\begin{eqnarray}&&\vec{F}=-\frac{G{{m}_{1}}{{m}_{2}}}{{{r}^{2}}}\widehat{{\vec{r}}},\end{eqnarray} \tag{ 2 }$$

where $r\equiv |\vec{r}|$ and $\widehat{{\vec{r}}}\equiv \vec{r}/r$ is the unit vector pointing from m₂ to m₁; see figure 1. The equations of motion for the two particles are given by

$$\begin{eqnarray}&&{{m}_{1}}\frac{{{{\rm d}}^{2}}{{{\vec{r}}}_{1}}}{{\rm d}{{t}^{2}}}=\vec{F}=-G\frac{{{m}_{1}}{{m}_{2}}}{{{r}^{2}}}\widehat{{\vec{r}}}=-{{m}_{2}}\frac{{{{\rm d}}^{2}}\overrightarrow{{{r}_{2}}}}{{\rm d}{{t}^{2}}}\;.\end{eqnarray} \tag{ 3 }$$

These equations are most conveniently solved by introducing the reduced mass $\mu ={{m}_{1}}{{m}_{2}}/({{m}_{1}}+{{m}_{2}})$ and rewriting (3) as the equation of motion for a single particle of mass μ

$$\begin{eqnarray}&&\mu \frac{{{{\rm d}}^{2}}\vec{r}}{{\rm d}{{t}^{2}}}=\vec{F}\;.\end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

Without loss of generality, we choose Cartesian coordinates $x,\;y,\;z$ such that the particles' motion takes place in the plane z = 0 and we furthermore introduce polar coordinates $r,\;\theta$ with $x=r{\rm cos} \theta$, $y=r{\rm sin} \theta$. We thus obtain two constants of motion, the energy E and the angular momentum L given by

$$\begin{eqnarray}&&E=\frac{1}{2}\mu \left[ {{\left( \frac{{\rm d}r}{{\rm d}t} \right)}^{2}}+{{r}^{2}}{{\left( \frac{{\rm d}\theta }{{\rm d}t} \right)}^{2}} \right]-G\frac{{{m}_{1}}{{m}_{2}}}{r},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&L=\mu {{r}^{2}}\frac{{\rm d}\theta }{{\rm d}t}\;.\end{eqnarray} \tag{ 6 }$$

By substituting ${\rm d}\theta /{\rm d}t$ in (5) in terms of L through (6) and solving the two equations for ${\rm d}r/{\rm d}t$ and ${\rm d}\theta /{\rm d}t$, respectively, we obtain a differential equation for r regarded now as a function of θ:

$$\begin{eqnarray}&&\frac{{\rm d}r}{{\rm d}\theta }=\frac{{\dot{r}}}{{\dot{\theta }}}=\frac{\sqrt{2\mu }{{r}^{2}}}{L}\sqrt{E-\frac{{{L}^{2}}}{2\mu {{r}^{2}}}+G\frac{{{m}_{1}}{{m}_{2}}}{r}},\end{eqnarray} \tag{ 7 }$$

where a 'dot' denotes a time derivative ${\rm d}/{\rm d}t$. A solution to equation (7) is given in closed analytic form by

$$\begin{eqnarray}&&r=\frac{{{r}_{0}}}{1+\epsilon {\rm cos} \theta },{{r}_{0}}=\frac{{{L}^{2}}}{\mu G{{m}_{1}}{{m}_{2}}},\epsilon =\sqrt{1+\frac{2E{{L}^{2}}}{\mu {{(G{{m}_{1}}{{m}_{2}})}^{2}}}},\end{eqnarray} \tag{ 8 }$$

where the semilatus rectum r₀ and the eccentricity are determined completely in terms of the constants of motion E, L.

The solutions to equation (7) can be classified into the following four types.

• (1)  
  Circular orbits given by $\epsilon =0$, where E takes on its minimal possible value ${{E}_{{\rm min} }}$ and the solutions are circles $r(\theta )={{r}_{0}}$.
• (2)  
  Kepler ellipses given by $0\lt \epsilon \lt 1$ or, equivalently ${{E}_{{\rm min} }}\lt E\lt 0$. In this case, the solution (8) can be written as

  $$\begin{eqnarray}&&\frac{{{(1-{{\epsilon }^{2}})}^{2}}}{r_{0}^{2}}{{\left( x+\frac{\epsilon {{r}_{0}}}{1-{{\epsilon }^{2}}} \right)}^{2}}+{{y}^{2}}\frac{1-{{\epsilon }^{2}}}{r_{0}^{2}}=1,\end{eqnarray} \tag{ 9 }$$

  which is of the general form ${{(x-{{x}_{0}})}^{2}}/{{a}^{2}}+{{(y-{{y}_{0}})}^{2}}/{{b}^{2}}=1$ for an ellipse centred on $({{x}_{0}},{{y}_{0}})$.
• (3)  
  Parabola given by $\epsilon =1\Leftrightarrow E=0$ in which case (8) takes on the form $2{{r}_{0}}x+{{y}^{2}}=r_{0}^{2}$.
• (4)  
  Hyperbolic orbits given by $\epsilon \gt 1\Leftrightarrow E\gt 0$ where the solution (8) can be written as $-({{\epsilon }^{2}}-1){{x}^{2}}+2{{r}_{0}}\epsilon x+{{y}^{2}}=r_{0}^{2}$.

The elliptic type of solution is often extended to include the circular case (1) and we have the three classic cone cross sections of elliptic, parabolic and hyperbolic particle curves.

We have seen that the Newtonian two-body problem can be formulated as one ordinary differential equation (7) for which the initial data at t = 0 need to be specified in the form of the initial position $(r,\theta )$ and the velocity components $\dot{r}$ and $\dot{\theta }$ or, alternatively, as ${\rm d}r/{\rm d}\theta$. The free parameters of the system are given by the masses m₁ and m₂ as well as the constants of motion E and L.

In order to illustrate the many fundamental differences that arise in the general relativistic two-body problem, it is helpful to first consider the Einstein equations (1) in a time–space split form. This is conveniently achieved with the canonical '3+1' split of the Einstein equations originally developed by Arnowitt, Deser and Misner (ADM) [4] and later reformulated by York [5, 6]; for a detailed review see [7]. We consider for this purpose a manifold $\mathcal{M}$ equipped with a spacetime metric ${{g}_{\alpha \beta }}$ and assume that there exists a function $t:\mathcal{M}\to \mathbb{R}$ that satisfies the following two properties. (i) The one-form ${\bf d}t$ is timelike everywhere, and (ii) the hypersurfaces ${{\Sigma }_{t}}$ defined by $t={\rm const}$ are non-intersecting and ${{\cup }_{t\in \mathbb{R}}}{{\Sigma }_{t}}=\mathcal{M}$. The resulting sequence of hypersurfaces ${{\Sigma }_{t}}$ is often referred to as a foliation of the spacetime and we denote by $n\equiv -{\bf d}t/||{\bf d}t||$ the future pointing unit normal field of the ${{\Sigma }_{t}}$. We furthermore define coordinates ${{x}^{\alpha }}$ to be adapted to the foliation if ${{x}^{0}}=t$ and the xⁱ, $i=1,\ldots ,3$, form a coordinate system in each hypersurface ${{\Sigma }_{t}}$.

It turns out to be convenient to define the lapse function and shift vector by

$$\begin{eqnarray}&&\alpha \equiv \frac{1}{||{\bf d}t||},{{\beta }^{\mu }}\equiv {{({{\partial }_{t}})}^{\mu }}-\alpha {{n}^{\mu }}\;.\end{eqnarray} \tag{ 10 }$$

Lapse and shift relate the unit normal direction $n$ to the direction ${{\partial }_{t}}$ of the coordinate time t which is illustrated in figure 2. One straightforwardly shows that $\langle {\bf d}t,\alpha n\rangle =1$ and, together with $\langle {\bf d}t,{{\partial }_{t}}\rangle =1$, it follows that the shift vector $\beta$ is tangent to ${{\Sigma }_{t}}$. Finally, the lapse function relates proper time τ as measured by an observer with four-velocity ${{n}^{\alpha }}$ to the coordinate time t: $\Delta \tau =\alpha \;\Delta t$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Having decomposed the spacetime into a one-parameter family of spatial hypersurfaces, we next consider projections of tensors. For this purpose, we define the projection operator ${{\bot }^{\alpha }}_{\mu }\equiv {{\delta }^{\alpha }}_{\mu }+{{n}^{\alpha }}{{n}_{\mu }}$ and the projection of an arbitrary tensor ${{T}^{{{\mu }_{1}}{{\mu }_{2}}\ldots }}_{{{\nu }_{1}}{{\nu }_{2}}\ldots }$ by

$$\begin{eqnarray}&&{{(\bot T)}^{{{\alpha }_{1}}{{\alpha }_{2}}\ldots }}_{{{\beta }_{1}}{{\beta }_{2}}\ldots }\equiv {{\bot }^{{{\alpha }_{1}}}}_{{{\mu }_{1}}}{{\bot }^{{{\alpha }_{2}}}}_{{{\mu }_{2}}}\ldots {{\bot }^{{{\nu }_{1}}}}_{{{\beta }_{1}}}{{\bot }^{{{\nu }_{2}}}}_{{{\beta }_{2}}}\ldots {{T}^{{{\mu }_{1}}{{\mu }_{2}}\ldots }}_{{{\nu }_{1}}{{\nu }_{2}}\ldots }\;.\end{eqnarray} \tag{ 11 }$$

In particular, the spatial projection of the metric gives us the first fundamental form or spatial metric

$$\begin{eqnarray}&&{{\gamma }_{\alpha \beta }}\equiv {{\bot }^{\mu }}_{\alpha }{{\bot }^{\nu }}_{\beta }{{g}_{\mu \nu }}={{g}_{\alpha \beta }}+{{n}_{\alpha }}{{n}_{\beta }}={{\bot }_{\alpha \beta }}\;.\end{eqnarray} \tag{ 12 }$$

$\gamma$ and $\bot$ thus represent the same tensor and we shall use both symbols depending on whether the emphasis is on the projection operation or the geometry of the spatial slices. It is straightforward to show that the components of the spacetime metric in adapted coordinates are related to the spatial metric, lapse and shift according to

$$\begin{eqnarray}{{g}_{\alpha \beta }}=\left( \begin{array}{ccccccccccccccc} -{{\alpha }^{2}}+{{\beta }_{m}}{{\beta }^{m}} & {{\beta }^{j}} \\ {{\beta }^{i}} & {{\gamma }_{ij}} \\ \end{array} \right)\Leftrightarrow {{g}^{\alpha \beta }}=\left( \begin{array}{ccccccccccccccc} -{{\alpha }^{-2}} & {{\alpha }^{-2}}{{\beta }^{j}} \\ {{\alpha }^{-2}}{{\beta }^{i}} & {{\gamma }^{ij}}-{{\alpha }^{-2}}{{\beta }^{i}}{{\beta }^{j}} \\ \end{array} \right).\end{eqnarray} \tag{ 13 }$$

Here, Latin indices $i,\;j,\;\ldots$ extend from 1 to 3 and spatial indices are raised and lowered with the spatial metric ${{\gamma }_{ij}}$ and its inverse ${{\gamma }^{ij}}$. The spatial metric furthermore defines a unique torsion-free and metric-compatible connection $\Gamma _{jk}^{i}=\frac{1}{2}{{\gamma }^{im}}({{\partial }_{j}}{{\gamma }_{km}}+{{\partial }_{k}}{{\gamma }_{mj}}-{{\partial }_{m}}{{\gamma }_{jk}})$ and an associated covariant derivative for arbitrary spatial tensors given by

$$\begin{eqnarray}&&{{D}_{\gamma }}{{S}^{{{\alpha }_{1}}{{\alpha }_{2}}\ldots }}_{{{\beta }_{1}}{{\beta }_{2}}\ldots }={{\bot }^{\lambda }}_{\gamma }{{\bot }^{{{\alpha }_{1}}}}_{{{\mu }_{1}}}{{\bot }^{{{\alpha }_{2}}}}_{{{\mu }_{2}}}\ldots {{\bot }^{{{\nu }_{1}}}}_{{{\beta }_{1}}}{{\bot }^{{{\nu }_{2}}}}_{{{\beta }_{2}}}\ldots {{\nabla }_{\lambda }}{{S}^{{{\mu }_{1}}{{\mu }_{2}}\ldots }}_{{{\nu }_{1}}{{\nu }_{2}}\ldots }\;.\end{eqnarray} \tag{ 14 }$$

The final ingredient we shall need in the space-time split of the Einstein equations is the second fundamental form or extrinsic curvature

$$\begin{eqnarray}&&{{K}_{\alpha \beta }}=-\bot {{\nabla }_{\beta }}{{n}_{\alpha }}\;.\end{eqnarray} \tag{ 15 }$$

Here, the minus sign is a common convention in NR but the extrinsic curvature is sometimes also defined with a plus sign in the literature. Furthermore, the definition (15) implies the relation ${{K}_{\alpha \beta }}=-\frac{1}{2}{{\mathcal{L}}_{n}}{{\gamma }_{\alpha \beta }}$, where $\mathcal{L}$ denotes the Lie derivative. It turns out to be convenient to also introduce the following projections of the energy momentum tensor

$$\begin{eqnarray}&&\rho ={{T}_{\mu \nu }}{{n}^{\mu }}{{n}^{\nu }},\;\;\;\;{{j}_{\alpha }}=-{{\bot }^{\nu }}_{\alpha }{{T}_{\mu \nu }}{{n}^{\mu }},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{S}_{\alpha \beta }}={{\bot }^{\mu }}_{\alpha }{{\bot }^{\nu }}_{\beta }{{T}_{\mu \nu }},\;\;\;\;S={{\gamma }^{\mu \nu }}{{S}_{\mu \nu }}.\end{eqnarray} \tag{ 17 }$$

The space–time split of the Einstein equations ${{G}_{\alpha \beta }}=8\pi {{T}_{\alpha \beta }}$ is then obtained through a lengthy calculation whose details can be found for example in [7]. This calculation gives six second-order in time evolution equations as well as the Hamiltonian and three momentum constraints

$$\begin{eqnarray}&&{{\partial }_{t}}{{\gamma }_{ij}}={{\beta }^{m}}{{\partial }_{m}}{{\gamma }_{ij}}+{{\gamma }_{mj}}{{\partial }_{i}}{{\beta }^{m}}+{{\gamma }_{im}}{{\partial }_{j}}{{\beta }^{m}}-2\alpha {{K}_{ij}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\partial }_{t}}{{K}_{ij}} & = & {{\beta }^{m}}{{\partial }_{m}}{{K}_{ij}}+{{K}_{mj}}{{\partial }_{i}}{{\beta }^{m}}+{{K}_{im}}{{\partial }_{j}}{{\beta }^{m}}-{{D}_{i}}{{D}_{j}}\alpha \\ {} & {} & +\alpha ({{\mathcal{R}}_{ij}}+K{{K}_{ij}}-2{{K}_{im}}{{K}^{m}}_{j})+4\pi \alpha [(S-\rho ){{\gamma }_{ij}}-2{{S}_{ij}}], \\ \end{array}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&0=\mathcal{R}+{{K}^{2}}-{{K}^{mn}}{{K}_{mn}}-16\pi \rho ,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&0={{D}_{i}}K-{{D}_{m}}{{K}^{m}}_{i}+8\pi {{j}_{i}}\;.\end{eqnarray} \tag{ 21 }$$

Here, ${{\mathcal{R}}_{ij}}$ and $\mathcal{R}$ denote the Ricci tensor and scalar associated with ${{\gamma }_{ij}}$. Note that we assume here coordinates $(t,\;{{x}^{i}})$ adapted to the foliation and therefore have replaced spacetime indices $\alpha ,\;\beta ,\;\ldots$ with spatial indices $i,\;j,\;\ldots \$. We also see that the extrinsic curvature K_ij allows us to write the second-order in time evolution equations as a first-order system. Finally, the constraint equations (20) and (21) are preserved under the evolution equations because of the Bianchi identities.

At this point it is worth taking a break to consider our situation in comparison with the Newtonian two-body problem. In place of one ordinary differential equation, we now have a system of coupled PDEs which forms an initial-boundary-value problem (IBVP). This system consists of six second-order in time evolution equations written in equations (18), (19) in first-order form as well as four constraints (20), (21). Even though the constraints are preserved under the evolution equations in the continuum limit, care needs to be taken that constraint violations due to numerical inaccuracies do not grow out of bounds. Furthermore, the initial data need to satisfy the constraints which requires solving a set of elliptic differential equations. Note that the Einstein equations in ADM form (18)–(21) make no predictions about the lapse function α and the shift vector β. Instead, these functions represent the diffeomorphism invariance or gauge freedom of GR; they can be freely specified to fix the coordinates but, as it turns out, it is highly non-trivial to find gauge conditions that ensure numerically stable evolutions.

It is instructive to count the number of physical degrees of freedom contained in the system (18)–(21). We have ten components of the Einstein metric ${{g}_{\alpha \beta }}$ corresponding to the ten functions ${{\gamma }_{ij}}$, ${{\beta }^{i}}$ and α in the ADM formulation. Four of these, the lapse and shift, are freely specifiable and do not contain physical information. The constraints impose four further conditions on the remaining functions ${{\gamma }_{ij}}$ that must be satisfied on each hypersurface ${{\Sigma }_{t}}$ and we are left with two gravitational degrees of freedom which correspond to the + and × GW polarization modes; see e.g. [8]. The two gravitational degrees of freedom are recovered even more elegantly in the characteristic formulation of the Einstein equations developed by Bondi, Sachs and collaborators [9, 10]. Here, one chooses at least one coordinate to be null and thus foliates spacetime in terms of light cones. The Einstein equations assume a natural hierarchy of two evolution equations, four hypersurface equations relating variables inside the hypersurfaces, three supplementary and one trivial equation; for details see [11] and references therein. Codes based on the characteristic formulation have been applied with great success in the presence of special spacetime symmetries and indeed been the first to model single BH spacetimes with long-term stability [12, 13]. In spite of the formalism's appealing properties, however, characteristic codes have as yet not been successfully generalized to BH binaries because the formation of caustics causes a breakdown of the coordinate system. It remains to be seen whether this obstacle can be overcome in future investigations; for a recent study see [14].

In the case of a non-vanishing energy momentum tensor ${{T}_{\alpha \beta }}$, there may be additional matter degrees of freedom. We also note that BH spacetimes with a BH mass M contain various different length scales, the BH horizon which has a size of² $\mathcal{O}(M)$, the wave length of GW signals, typically of the order $\mathcal{O}({{10}^{2}}\;M)$, and the wave zone of $\mathcal{O}({{10}^{3}}\;M)$ where perturbation theory permits a precise definition of GWs. Finally, we need to specify outer boundary conditions such that there is no ingoing gravitational radiation from infinity.

Bearing in mind all these features of the Einstein equations, we face the following list of tasks to obtain stable, accurate numerical simulations of the binary BH problem in GR.

• Formulate the Einstein equations in a manner that admits a well-posed IBVP, i.e. ensures a continuous dependence of the spacetime solution on the initial data.
• Choose numerically suitable gauge conditions.
• Discretize the resulting PDEs for a computer based treatment.
• Specify physically correct boundary conditions that also satisfy the constraints.
• Find a numerical treatment of the singularities inherent in BH spacetimes that avoids the appearance of non-assigned numbers.
• Calculate initial data which satisfy the constraints and represent a realistic snapshot of the initial state of the physical system under consideration.
• Implementation of mesh refinement or similar methods using multiple domains to accurately handle the different length scales and parallelize the resulting algorithms for multi-processor computation.
• Extract physical results in a gauge-invariant manner from the numerical data.

In the next section we will discuss the most important methods which have been developed for handling these tasks and have made possible the NR breakthroughs in solving the binary BH problem in GR.

The techniques for addressing the above list of tasks have been developed over several decades through an interplay of numerical experiments and theoretical studies carried out by numerous groups and researchers. Many numerical investigations, especially the earlier ones, were performed without a comprehensive understanding of all the difficulties associated with this list of tasks. In hindsight, it is therefore not too surprising that they met with limited success. And yet, as is the nature of scientific exploration, all these attempts taught us valuable lessons and contributed to the gradual assembly of the complete picture we are going to describe in this section. We shall present this more technical description not in chronological order but, for reasons of clarity, topic by topic. A brief historical review of the applications will be given in section 5 below.

4.1. Formulations of the Einstein equations

Any successful attempt at numerically solving the Einstein equations must be based on a well-posed IBVP, i.e. a computational algorithm that results in a time evolution which depends continuously on the initial data. Given the inevitability of numerical noise present in the form of round-off error in any numerical initial data, algorithms that do not meet this criterion are evidently unsuitable for obtaining reliable results. The well-posedness of a numerical implementation depends on many aspects including the specific formulation of the differential equations, gauge and boundary conditions and the discretization schemes. Here we are concerned with the conditions a formulation of the Einstein equations must satisfy to admit a well-posed IBVP.

The suitability of a formulation is commonly studied in the form of the hyperbolicity properties of the PDEs which are related to the symmetrizeability of the principal symbol (in simple words, the coefficient matrix of the terms containing the highest derivatives) of the system of PDEs. A system is called strongly hyperbolic if the principal symbol has only imaginary eigenvalues and a complete set of linearly independent eigenvectors [15]. If the eigenvectors are not linearly independent, the system is called weakly hyperbolic. A system is called symmetric hyperbolic if there exists a conserved, positive energy norm. We skip the technical details here, but the interested reader will find extended discussions in [16, 17] and references therein. For us, the most important conclusions are the following. (i) Of the three notions of weak, strong and symmetric hyperbolicity, each is a stronger condition than the previous one; see section 3.1.4 in [17]. (ii) Strong hyperbolicity is a necessary condition for a well-posed IBVP [18, 19]. (iii) The ADM evolution equations (18), (19) have been shown to be weakly but not strongly hyperbolic for fixed gauge [15] and a first-order version of the ADM equations has been shown to be only weakly hyperbolic in [20]. While these studies do not constitute a rigorous proof of the unsuitability of the ADM equations for numerical evolutions, they strongly suggest a search for alternative formulations for which strong hyperbolicity can be established.

Explorations of modifications of the ADM equations or alternative formulations for use in NR already began in the late 1980s, before the full impact of the hyperbolicity properties of the different formulations had been realized. Over the course of the ensuing 25 years, a great variety of different formulations has been developed and implemented in numerical codes; for an overview see for example [21] and in particular figure 4 therein. Here, we shall focus on those two formulations that underlie the NR breakthroughs of 2005.

Zoom In Zoom Out Reset image size

The Baumgarte–Shapiro–Shibata–Nakamura (BSSN) formulation [22–24] is directly derived from the ADM equations but works with conformally rescaled variables, a trace split of the extrinsic curvature and promotes the contracted Christoffel symbols to the status of independent variables. Specifically, the BSSN variables χ, ${{\tilde{\gamma }}_{ij}}$, K, ${{\tilde{A}}_{ij}}$ and ${{\tilde{\Gamma }}^{i}}$ are defined by

$$\begin{eqnarray}&&\chi ={{\gamma }^{-1/3}},{{\tilde{\gamma }}_{ij}}=\chi {{\gamma }_{ij}},K={{\gamma }^{mn}}{{K}_{mn}},{{\tilde{A}}_{ij}}=\chi \left( {{K}_{ij}}-\frac{1}{3}{{\gamma }_{ij}}K \right),{{\tilde{\Gamma }}^{i}}={{\tilde{\gamma }}^{mn}}\tilde{\Gamma }_{mn}^{i},\end{eqnarray} \tag{ 22 }$$

where $\gamma \equiv {\rm det} \;{{\gamma }_{ij}}$ and $\tilde{\Gamma }_{jk}^{i}$ are the Christoffel symbols associated with the conformal metric ${{\tilde{\gamma }}_{ij}}$. The BSSN system has also been used with χ replaced by the variables $\phi =-({\rm ln} \;\chi )/4$ or $W=\sqrt{\chi }$. Inserting the definition (22) into the ADM equations (18), (19) and using the Hamiltonian and momentum constraints respectively in the resulting evolution equations for K and ${{\Gamma }^{i}}$ yields

$$\begin{eqnarray}&&{{\partial }_{t}}\chi ={{\beta }^{m}}{{\partial }_{m}}\chi +\frac{2}{3}\chi (\alpha K-{{\partial }_{m}}{{\beta }^{m}}),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\tilde{\gamma }}_{ij}}={{\beta }^{m}}{{\partial }_{m}}{{\tilde{\gamma }}_{ij}}+2{{\tilde{\gamma }}_{m(i}}{{\partial }_{j)}}{{\beta }^{m}}-\frac{2}{3}{{\tilde{\gamma }}_{ij}}{{\partial }_{m}}{{\beta }^{m}}-2\alpha {{\tilde{A}}_{ij}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\partial }_{t}}K & = & {{\beta }^{m}}{{\partial }_{m}}K-\chi {{\tilde{\gamma }}^{mn}}{{D}_{m}}{{D}_{n}}\alpha +\alpha {{\tilde{A}}^{mn}}{{\tilde{A}}_{mn}}+\frac{1}{3}\alpha {{K}^{2}} \\ {} & {} & +4\pi \alpha [S+\rho ], \\ \end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\partial }_{t}}{{\tilde{A}}_{ij}} & = & {{\beta }^{m}}{{\partial }_{m}}{{\tilde{A}}_{ij}}+2{{\tilde{A}}_{m(i}}{{\partial }_{j)}}{{\beta }^{m}}-\frac{2}{3}{{\tilde{A}}_{ij}}{{\partial }_{m}}{{\beta }^{m}}+\alpha K{{\tilde{A}}_{ij}} \\ {} & {} & -2\alpha {{\tilde{A}}_{im}}{{\tilde{A}}^{m}}_{j}+\chi {{\left( \alpha {{\mathcal{R}}_{ij}}-{{D}_{i}}{{D}_{j}}\alpha -8\pi \alpha {{S}_{ij}} \right)}^{{\rm TF}}}, \\ \end{array}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\partial }_{t}}{{\tilde{\Gamma }}^{i}} & = & {{\beta }^{m}}{{\partial }_{m}}{{\tilde{\Gamma }}^{i}}+\frac{2}{3}{{\tilde{\Gamma }}^{i}}{{\partial }_{m}}{{\beta }^{m}}-{{\tilde{\Gamma }}^{m}}{{\partial }_{m}}{{\beta }^{i}}+{{\tilde{\gamma }}^{mn}}{{\partial }_{m}}{{\partial }_{n}}{{\beta }^{i}} \\ {} & {} & +\frac{1}{3}{{\tilde{\gamma }}^{im}}{{\partial }_{m}}{{\partial }_{n}}{{\beta }^{n}}-{{\tilde{A}}^{im}}\left[ 3\alpha \frac{{{\partial }_{m}}\chi }{\chi }+2{{\partial }_{m}}\alpha \right]+2\alpha \tilde{\Gamma }_{mn}^{i}{{\tilde{A}}^{mn}} \\ {} & {} & -\frac{4}{3}\alpha {{\tilde{\gamma }}^{im}}{{\partial }_{m}}K-16\pi \alpha {{j}^{i}}. \\ \end{array}\end{eqnarray} \tag{ 27 }$$

Here, 'TF' denotes the tracefree part and ${{\mathcal{R}}_{ij}}$ the Ricci tensor associated with the physical three-metric ${{\gamma }_{ij}}$. The promotion of auxiliary variables to independent status introduces three additional constraints to the BSSN system given by

$$\begin{eqnarray}&&{\rm det} \;{{\tilde{\gamma }}_{ij}}=1,{{\tilde{\gamma }}^{mn}}{{\tilde{A}}_{mn}}=0,{{\mathcal{G}}^{i}}\equiv {{\tilde{\Gamma }}^{i}}-{{\tilde{\gamma }}^{mn}}\tilde{\Gamma }_{mn}^{i}=0.\end{eqnarray} \tag{ 28 }$$

In practical applications, it turns out necessary for numerical stability to control these auxiliary constraints in the following manner. (i) Enforce ${{\tilde{\gamma }}^{mn}}{{\tilde{A}}_{mn}}=0$ and (ii) either add the constraint ${{\mathcal{G}}^{i}}$ to the right-hand side of equation (27) [25] or, alternatively, substitute on the right-hand side of equation (27) all ${{\tilde{\Gamma }}^{i}}$ that appear in undifferentiated form by their definition in terms of the metric ${{\tilde{\gamma }}_{ij}}$ [26].

Empirical studies quickly demonstrated that the BSSN system provides superior numerical stability when compared with the ADM equations (e.g. [24]) and mathematical studies demonstrated BSSN to provide a strongly hyperbolic formulation of the Einstein equations [27]. The BSSN formulation is employed in the binary BH breakthroughs by the Brownsville/Rochester and the NASA Goddard groups [2, 3].

The other formulation instrumental for the breakthroughs is based on the Einstein equations in harmonic gauge [28] defined by the spacetime coordinates satisfying the condition $\square {{x}^{\alpha }}=-{{g}^{\mu \nu }}\Gamma _{\mu \nu }^{\alpha }=0$. In this form, the Ricci tensor takes on the form

$$\begin{eqnarray}&&{{R}_{\alpha \beta }}=\;-\frac{1}{2}{{g}^{\mu \nu }}{{\partial }_{\mu }}{{\partial }_{\nu }}{{g}_{\alpha \beta }}+\ldots ,\end{eqnarray} \tag{ 29 }$$

where the dots denote terms containing at most first derivatives of the spacetime metric. In this form, the principal part of the Einstein equations is that of the scalar wave operator and the equations are symmetric hyperbolic. Harmonic coordinates have been used in the first proofs of the local uniqueness of the Cauchy problem in GR [29–31]. A generalization of this particularly appealing form of the Einstein equations to arbitrary gauge is realized by promoting the functions

$$\begin{eqnarray}&&{{H}^{\alpha }}\equiv \square {{x}^{\alpha }}=\;-{{g}^{\mu \nu }}\Gamma _{\mu \nu }^{\alpha },\end{eqnarray} \tag{ 30 }$$

to the status of independently evolved variables [32, 33]. The resulting system is often referred to as the generalized harmonic gauge (GHG) formulation and considers the generalized set of equations

$$\begin{eqnarray}&&{{R}_{\alpha \beta }}-{{\nabla }_{(\alpha }}{{\mathcal{C}}_{\beta )}}=8\pi \left( {{T}_{\alpha \beta }}-\frac{1}{2}T{{g}_{\alpha \beta }} \right),\end{eqnarray} \tag{ 31 }$$

with the auxiliary constraints ${{\mathcal{C}}^{\alpha }}\equiv {{H}^{\alpha }}-\square {{x}^{\alpha }}$. A solution to the Einstein equations is obtained by solving equation (31) subject to the condition ${{\mathcal{C}}^{\alpha }}=0$. In practice, this is conveniently achieved by prescribing initial data for ${{g}_{\alpha \beta }}$ and ${{\partial }_{t}}{{g}_{\alpha \beta }}$ and initializing the ${{H}^{\alpha }}$ through equation (30). If the initial data furthermore satisfy the Hamiltonian and momentum constraints (20), (21), this can be shown to imply ${{\partial }_{t}}{{\mathcal{C}}^{\alpha }}=0$. The Bianchi identities then ensure that the auxiliary constraint is preserved under time evolution so that ${{\mathcal{C}}^{\alpha }}=0$ at all times in the continuum limit. For controlling violations of these constraints at the discretized level in numerical evolutions, Gundlach et al [34] suggested the addition of constraint damping terms which turned out crucial in achieving the numerical stability required for binary BH simulations [1]. With these terms, the generalized Einstein equations (31) can be written in the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{g}^{\mu \nu }}{{\partial }_{\mu }}{{\partial }_{\nu }}{{g}_{\alpha \beta }} & = & -2{{\partial }_{\nu }}{{g}_{\mu (\alpha \,{{\partial }_{\beta )}}}}{{g}^{\mu \nu }}-2{{\partial }_{(\alpha }}{{H}_{\beta )}}+2{{H}_{\mu }}\Gamma _{\alpha \beta }^{\mu }-2\Gamma _{\nu \alpha }^{\mu }\Gamma _{\mu \beta }^{\nu } \\ {} & {} & -8\pi {{T}_{\alpha \beta }}+4\pi T{{g}_{\alpha \beta }}-2\kappa \left[ 2{{n}_{(\alpha }}{{\mathcal{C}}_{\beta )}}-\lambda {{g}_{\alpha \beta }}{{n}^{\mu }}{{\mathcal{C}}_{\mu }} \right]. \\ \end{array}\end{eqnarray} \tag{ 32 }$$

Here, κ and λ are user-specified constant parameters which control the constraint damping. The GHG formulation has been used in Pretorius' breakthrough simulations [1]; see also [35].

4.2. Gauge conditions

In the previous section, we have seen that some of the evolution variables are not determined by the Einstein equations. The lapse function α and the shift vector ${{\beta }^{i}}$ are freely specifiable in the ADM system (18)–(21) or the BSSN equations (23)–(27) and the ${{H}^{\alpha }}$ are undetermined in the GHG formulation (32). Instead, these functions represent the coordinate or gauge freedom of GR, and their choice leaves the physical properties of the spacetime invariant. As one might expect from this shared property, the two sets of gauge functions are related; the normal component and spatial projection of the ${{H}^{\alpha }}$ can be expressed in terms of the lapse α and shift ${{\beta }^{i}}$ respectively as given in equation (18), (19) of [35]. The geometrical meaning of the gauge functions is more intuitively encoded in lapse and shift and most investigations into the numerical properties of different gauge choices have been carried out in terms of these variables.

The simplest choice would appear to be given by $\alpha =1$ and ${{\beta }^{i}}=0$, referred to as geodesic slicing with vanishing shift. The problems of using this gauge in numerical simulations, however, have been demonstrated as early as 1978 by Smarr and York's [36] time evolutions of the time symmetric slice of the Kruskal–Schwarzschild spacetime. Setting $\alpha =1$ and $\beta =0$, the numerically constructed hypersurfaces encounter the BH singularity after an evolution time $t=\pi M$; see the upper panel of their figure 2. This behaviour highlights one key requirement to be met by any numerically suitable set of gauge conditions: The evolution in proper time should be slowed down in regions where the hypersurfaces approach a singularity. This feature is commonly referred to as singularity avoiding slicing and has been suggested first in the form of maximal slicing K = 0 [37]. Singularity avoidance is achieved by letting the lapse function vary in space and time such that it drops towards zero in the vicinity of spacetime singularities. For an illustration of this effect, we refer again to figure 2 in [36] which contrasts maximal with geodesic slicing. A wider class of singularity avoiding slicings has been studied in the form of the Bona–Massó family [38] which includes maximal slicing as a special case; see also [39] and, for the case of harmonic coordinates, [33]. It has been noticed, however, that due to the different advance in proper time in different regions of the spacetime during the numerical evolution, neighbouring grid points of a computational domain may correspond to increasingly distant points in the spacetime. This phenomenon, often referred to as grid or slice stretching, needs to be cured by a 'suitable shifting of grid points' through the use of a non-zero shift vector; see e.g. [40].

Geometrically motivated shift conditions were used in the already mentioned work by Smarr and York [36] in the form of the minimal distortion gauge. Considering for example a small sphere on a given hypersurface ${{\Sigma }_{t}}$, it can be shown that the minimal distortion gauge preserves the spherical shape at leading order whereas in general the shape will be sheared into an ellipse. The maximal slicing condition furthermore preserves the volume of the sphere. The numerical implementation of this shift is complicated by the necessity to solve elliptic equations for the ${{\beta }^{i}}$; for details see section 4 and, in particular equation (4.11) in [36]. In practice, it is much simpler to evolve the shift in time according to parabolic or hyperbolic differential equations which can be achieved with so-called 'driver' conditions [41]. Alcubierre et al [42] have obtained such equations for the shift by relating ${{\partial }_{t}}{{\beta }^{i}}$ or $\partial _{t}^{2}{{\beta }^{i}}$ to the elliptic operator obtained from the 'Gamma freezing' condition ${{\partial }_{t}}{{\tilde{\Gamma }}^{i}}=0$, where $\tilde{\Gamma }$ is the BSSN variable defined in equation (22). They thus arrive at the hyperbolic 'Γ-driver' condition $\partial _{t}^{2}{{\beta }^{i}}=\zeta {{\partial }_{t}}{{\tilde{\Gamma }}^{i}}-\xi {{\partial }_{t}}{{\beta }^{i}}$, where ζ and ξ are specifiable positive functions. In a similar way, the Bona–Massó family replaces the elliptic maximal slicing condition K = 0 with an 'easier to implement' hyperbolic condition ${{\partial }_{t}}\alpha -{{\alpha }^{2}}f(\alpha )[K-K(t=0)]$ where $f(\alpha )$ is a positive function. By using a specific version of this slicing and Γ-driver shift, Alcubierre et al [42] managed to extract GWs from the evolution of a distorted BH and drive the coordinates to a frame where the system remains almost static at late times. A specific version of the Bona–Massó family is the so-called '1+log' slicing which sets $f(\alpha )=2/\alpha$ and enabled Alcubierre et al [40] to evolve BH data for long times above $1000\;M$.

The breakthrough simulations of [2, 3] obtained with the BSSN formulation employ variants of the 1+log slicing and the Γ-driver shift condition given by

$$\begin{eqnarray}&&{{\partial }_{t}}\alpha ={{\beta }^{m}}{{\partial }_{m}}-2\alpha K,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\beta }^{i}}={{\beta }^{m}}{{\partial }_{m}}{{\beta }^{i}}+\frac{3}{4}{{B}^{i}},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{B}^{i}}={{\beta }^{m}}{{\partial }_{m}}{{B}^{i}}+{{\partial }_{t}}{{\tilde{\Gamma }}^{i}}-\eta {{B}^{i}},\end{eqnarray} \tag{ 35 }$$

or some minor modification of these equations; see [43]. Here, η is a user specified constant or function of the coordinates.

The GHG formulation is motivated by the beneficial properties of the Einstein equations in harmonic gauge, but stable numerical evolutions of binary BH spacetimes have so far required at least some component of the ${{H}^{\alpha }}$ to be non-zero. Pretorius [1] sets H_i = 0 and evolves the t component according to

$$\begin{eqnarray}&&\square {{H}_{t}}=\;-{{\xi }_{1}}\frac{\alpha -1}{{{\alpha }^{\eta }}}+{{\xi }_{2}}{{n}^{\mu }}{{\partial }_{\mu }}{{H}_{t}},\end{eqnarray} \tag{ 36 }$$

where ${{\xi }_{1}}=19/m$, ${{\xi }_{2}}=2.5/m$, $\eta =5$ and m denotes the mass of one of the two (equal-mass) BHs. This choice prevents the lapse from deviating too much from unity which may have caused instabilities in earlier GHG evolutions [44]. For some further discussion of gauge choices in the GHG formulation, see e.g. [45, 46].

4.3. Boundary conditions and singularity treatment

Astrophysical BHs are commonly modeled as asymptotically flat spacetimes, i.e. the spacetime approaches the Minkowski limit far away from the BH regions of strong curvature. Strictly speaking, this is an approximation to the cosmological spacetimes that describe our universe, but for most practical applications, as for example the modeling of GW signals expected to be observed with laser interferometric detectors, it is sufficient to include cosmological effects in the form of a redshift factor $1+z$ multiplying the source mass. Asymptotically flat spacetimes are of infinite extent and the challenge in NR is to describe these inside compact computational domains. The most elegant way to achieve this goal is to compactify the spacetime coordinates and cover dimensions of infinite extent with a finite coordinate interval as for example using maps of the type

$$\begin{eqnarray*}&&r\in [0,\infty )\to x\equiv 1/(r+1)\in (0,1].\end{eqnarray*}$$

Applied to 3+1 splits of the spacetime, however, this often results in an asymptotically infinite blue shifting of gravitational radiation; the wavelength of the radiation asymptotically shrinks to zero as measured in the compactified coordinate and fails to be resolved numerically. Characteristic formulations of the Einstein equations [11], on the other hand, are ideally suited for such a treatment, as the coordinates are constructed in terms of light cones, GW signals have constant phase along the characteristic coordinate curves and no resolution problems arise. Here lies one of the attractive features of characteristic formulations mentioned above in section 3. Inside 3+1 formulations, such behaviour can be obtained by slicing the spacetime with hypersurfaces that are spacelike everywhere, but become asymptotically null at infinity. This type of slicing can be obtained for example using hyperbolic slicing $K={\rm const}\ne 0$ and plays an important role in the so-called conformal field equations (see [47] and references therein), but has, to our knowledge, not yet been applied successfully to simulate BH binaries.

In practice, most numerical applications model only a finite subset of the total spacetime and impose boundary conditions at large but finite distance from the BHs. Ideally, such boundary conditions satisfy the following three requirements. (i) Ensure well posedness of the IBVP, (ii) compatibility with the Einstein constraint equations, and (iii) a correct representation of the physical boundary conditions, typically minimization of the ingoing gravitational radiation [48]. Such conditions have been studied mostly for the GHG formulation; see [49–51] and references therein.

Boundary conditions meeting the above criteria have not yet been developed for the BSSN formulation³ and numerical applications of the BSSN system therefore resort to an approximation using outgoing radiation or Sommerfeld conditions. The assumption underlying this approach is that the evolution variables f approach a constant background value f₀ far away from the strong-field sources and deviations from this value at finite radius r can be written as $f={{f}_{0}}+u(t-r)/{{r}^{n}}$ with a positive, integer n. The outgoing radiation condition ${{\partial }_{t}}u+{{\partial }_{r}}u=0$ for the radiative deviations then translates into the boundary condition [40]

$$\begin{eqnarray}&&{{\partial }_{t}}f+n\frac{f-{{f}_{0}}}{r}+\frac{{{x}^{i}}}{r}{{\partial }_{i}}f=0,\end{eqnarray} \tag{ 37 }$$

where xⁱ denote Cartesian coordinates and ${{r}^{2}}={{\sum }_{i}}{{({{x}^{i}})}^{2}}$. These conditions are not without problems: (i) the system is over-determined because the number of conditions imposed exceeds that of the ingoing characteristics; (ii) Sommerfeld conditions are not constraint satisfying, and (iii) the non-exact outgoing nature of these conditions at finite radii may introduce spurious reflections. In spite of these caveats, Sommerfeld conditions turn out to work rather well and robustly in many numerical applications (see e.g. [50]) and are the method of choice for the moving puncture breakthroughs of the Brownsville/RIT and Goddard groups [2, 3]. Pretorius [1], instead, uses a compactified domain and overcomes the problem of under-resolving the blue-shifted radiation by damping the radiation through numerical viscosity and thus effectively emulates no-ingoing-radiation conditions.

A second type of boundary condition arises in NR applications through the presence of the spacetime singularities. These singularities typically manifest themselves in the form of diverging or vanishing metric components as for example the ${{g}_{rr}}={{(1-2\;M/r)}^{-1}}$ in the Schwarzschild metric in Schwarzschild coordinates. Computer simulations react with non-assigned numbers to the resulting infinities which rapidly swamp the entire computational domain and render the simulation practically useless. One elegant approach to handle this problem employs so-called 'puncture' initial data (see section 4.5) and factors out the singular part of the BH data throughout the time evolution; see e.g. [40, 53] for details of these fixed puncture evolutions. In this approach, the BHs remain at fixed coordinate location throughout the evolution and it appears to be difficult to construct long-term stable coordinate conditions for BH inspiral in this approach; see [54] for the most advanced application of this type leading to about one orbit of BH inspiral.

An alternative method to handle singularities which has become popular over the years is the BH or singularity excision technique attributed to Unruh [55]. By virtue of Penrose's cosmic censorship conjecture (see e.g. [56]), spacetime singularities should be cloaked inside an event horizon such that the spacetime exterior to the horizon is causally disconnected from events inside the horizon. In particular, the exterior spacetime should not be affected by completely removing a finite region around the singularity from the numerical evolution as long as the excised region remains inside the event horizon or, as usually done in practice, is located inside the apparent horizon (AH) [57] on each hypersurface ${{\Sigma }_{t}}$. This is illustrated in figure 3 where the large circle represents the AH and the region consisting of the white (empty) small circles is removed from the computational domain whereas black (filled) points are updated regularly. The update in time at a given grid point requires data from neighbouring points to evaluate spatial derivatives, and therefore the excision boundary (grey circles) may require some special treatment. This can be achieved either by sideways differencing operators [35], extrapolation from data on grid points further out [58] or calculating the function values from the spectral expansion in spectral codes [59]. The first of these schemes is the method employed in Pretorius' work [1]. Rather astonishingly, the moving puncture method [2, 3] does not implement an explicit excision scheme but instead uses finite differencing stencils right across the BH singularities. The surprising success of this method has been explored in more depth in [60–62] and references therein. The singularity of puncture type initial data is a coordinate singularity that contains spatial infinity of the far side of the wormhole geometry compactified into a single point. In moving puncture evolutions, however, these initial data rapidly change from a wormhole to a so-called 'trumpet' geometry which is only partially covered by the computational domain because of the discrete structure of the numerical grid; see figure 1 in [61]. The singularity, instead, 'falls through the grid' and the moving puncture technique can therefore be interpreted as an indirect excision method provided by the finite grid resolution.

4.4. Discretization and mesh refinement

4.5. Initial data

The task in generating initial data in NR is two-fold. (i) The data must satisfy the Einstein constraint equations (20), (21) and (ii) they need to represent a realistic snapshot of the physical system under study. A natural starting point for the construction of initial data is to apply modifications to existing analytic BH solutions. An analytic solution of particular relevance for this purpose is the Schwarzschild solution in isotropic coordinates

$$\begin{eqnarray}&&{\rm d}{{s}^{2}}=\;-{{\left( \frac{2r-M}{2r+M} \right)}^{2}}{\rm d}{{t}^{2}}+{{\left( 1+\frac{M}{2r} \right)}^{4}}\left( {\rm d}{{r}^{2}}+{{r}^{2}}{\rm d}{{\theta }^{2}}+{{r}^{2}}{{{\rm sin} }^{2}}\theta \;{\rm d}{{\phi }^{2}} \right).\end{eqnarray} \tag{ 38 }$$

Its scope for generalization becomes clear in a systematic approach to solving the Einstein constraints based on the York–Lichnerowicz split [85–87]. This method consists of a conformal transformation of the spatial metric ${{\gamma }_{ij}}={{\psi }^{4}}{{\bar{\gamma }}_{ij}}$ and a conformal traceless split of the extrinsic curvature

$$\begin{eqnarray}&&{{K}_{ij}}={{A}_{ij}}+\frac{1}{3}{{\gamma }_{ij}}K,\quad {{A}^{ij}}={{\psi }^{-10}}{{\bar{A}}^{ij}},\quad {{A}_{ij}}={{\psi }^{-2}}{{\bar{A}}_{ij}}\;.\end{eqnarray} \tag{ 39 }$$

By further decomposing ${{\bar{A}}_{ij}}$ into a longitudinal piece plus a transverse traceless part, the momentum constraints simplify considerably [88]. Similar simplifications are achieved by using instead of equation (39) a physical transverse traceless split or the so-called thin-sandwich decomposition; for details see [88–90] and references therein.

The Hamiltonian and momentum constraints take on a particularly simple form if one further requires that (i) the trace of the extrinsic curvature vanishes K = 0, (ii) the conformal metric is flat ${{\bar{\gamma }}_{ij}}={{f}_{ij}}$ where f_ij is the Euclidean metric (not necessarily in Cartesian coordinates), and (iii) the spatial metric is asymptotically flat, ${{{\rm lim} }_{r\to \infty }}\psi =1$. Then the momentum constraints decouple from the Hamiltonian constraint and form a set of equations for the ${{\bar{A}}_{ij}}$. The simplest solution for these equations is the trivial one ${{\bar{A}}_{ij}}=0$ in which case the Hamiltonian constraint becomes $\bar{\Delta }\psi =0$, where $\bar{\Delta }$ is the Laplace operator associated with the flat metric f_ij. This Laplace equation for the conformal factor ψ is solved by the spatial part of the Schwarzschild metric (38) ${{\gamma }_{ij}}={{\psi }^{4}}{{f}_{ij}}$ with $\psi =1+M/(2r)$. By linearity of the Laplace equation, we can obtain initial data containing multiple BHs by superposing N solutions of this type according to $\psi =1+\sum _{A=1}^{N}{{m}_{A}}/(2|\vec{r}-{{\vec{r}}_{A}}|)$, where m_A and ${{\vec{r}}_{A}}$ denote the mass and location of the BHs. This solution is known as Brill–Lindquist [91] data and represents a snapshot of N BHs at the moment of time symmetry. A similar type of data differing only in the symmetry conditions at the throat of the worm hole(s) has been constructed by Misner [92]. Both, Brill–Lindquist and Misner data formed the starting point of many BH evolutions over the previous decades.

Quite remarkably, under the assumption of conformal flatness and K = 0, the momentum constraints even admit non-vanishing analytic solutions for the extrinsic curvature, the Bowen–York [93] data

$$\begin{eqnarray}&&{{\bar{A}}_{ij}}=\frac{3}{2{{r}^{2}}}\left[ {{P}_{i}}{{n}_{j}}+{{P}_{j}}{{n}_{i}}-({{f}_{ij}}-{{n}_{i}}{{n}_{j}}){{P}^{k}}{{n}_{k}} \right]+\frac{3}{{{r}^{3}}}\left( {{\epsilon }_{kil}}{{S}^{l}}{{n}^{k}}{{n}_{j}}+{{\epsilon }_{kjl}}{{S}^{l}}{{n}^{k}}{{n}_{i}} \right),\end{eqnarray} \tag{ 40 }$$

where r is the areal radius associated with the flat metric f_ij, nⁱ the unit, outgoing radial vector and Pⁱ, Sⁱ are free parameters that correspond to the total linear and angular momentum of the initial hypersurface [94]. The momentum constraints are linear in ${{\bar{A}}_{ij}}$, so that multiple solutions of the type (40) can be superposed. Equation (40) gives the generalization of Brill–Lindquist data for non-zero BH momenta. For generalization of Misner data, one needs to construct inversion-symmetric data of the type (40) using the method of images (see section 3.2.1 in [88]).

In the conformal-flatness approximation with K = 0 and non-vanishing Bowen–York extrinsic curvature (40), the Hamiltonian constraint becomes

$$\begin{eqnarray}&&\bar{\Delta }\psi +\frac{1}{8}{{K}^{mn}}{{K}_{mn}}{{\psi }^{-7}}=0\;.\end{eqnarray} \tag{ 41 }$$

This elliptic equation is often solved by decomposing the conformal factor into a Brill–Lindquist piece ${{\psi }_{{\rm BL}}}=1+M/(2r)$ plus a regular contribution u. Under such decomposition, the existence and uniqueness of C² regular solutions u to equation (41) has been proven by Brandt and Brügmann [95] and the data thus constructed are commonly referred to as puncture data. The Schwarzschild solution (38) in isotropic coordinates is the simplest non-trivial solution of this type, a single BH with zero linear and angular momentum. Alternative to the puncture method, the initial data formalism summarized here has also been applied, in some flavour or other, to the construction of BH excision data; see e.g. [67, 83, 96].

In spite of their popularity in time evolutions, the conformal flatness nature of puncture data results in some restrictions. In particular, the Kerr [97] spacetime describing a single rotating BH does not contain a maximal, conformally flat hypersurface [98, 99]. Puncture data with non-zero Bowen–York angular momentum Sⁱ therefore not only contain a rotating BH but some further gravitational fields which manifest themselves as pulses of spurious radiation colloquially referred to as 'junk radiation'. While this spurious GW pulse is often small, it increases nonlinearly with the Bowen–York parameters Pⁱ and Sⁱ. In particular, this imposes a practical limit of the initial dimensionless spin parameter of BH configurations of $\approx 0.93$ [100, 101] and has motivated the construction of initial data without the conformal-flatness assumption by either applying most of the puncture formalism to a non-conformally flat background metric [102, 103] or applying the conformal-thin-sandwich method to a background of superposed Kerr–Schild data [104, 105].

Puncture data have been the starting point for the majority of BH binary simulations in the past decade and, as suggested by the name, were also used in the moving puncture breakthroughs [2, 3]. A conceptually rather different approach was used in Pretorius' simulations. Instead of using initial data containing BHs, the simulations start with matter in the form of scalar field clouds concentrated and boosted such that they rapidly collapse into a BH with velocity corresponding approximately to a binary configuration in quasi-circular orbit [1].

4.6. Diagnostics

The extraction of physical information from a BH binary simulation is typically a diagnostic process that uses the numerically constructed fields but has at most minor impact⁴ on the actual time integration of the fields. Their significance for BH binaries therefore mostly consists in understanding the results rather than solving the two-body problem itself and we shall only briefly summarize here the most important diagnostic quantities but provide references for more details. The diagnostic quantities can be loosely classified into three groups: properties of the global spacetimes, the GW signal and the BH horizons.

Global quantities.

For asymptotically flat spacetimes, the total mass-energy and the linear momentum of a spacetime is given by the ADM mass and momentum [4]. If the coordinate system is chosen such that in the limit of infinite separation r from the strong-field sources the spacetime metric deviates from the Minkowski metric ${{\eta }_{\mu \nu }}$ according to ${{g}_{\mu \nu }}={{\eta }_{\mu \nu }}+\mathcal{O}(1/r)$, the ADM mass and momentum is given in terms of the ADM variables by the integrals

$$\begin{eqnarray}&&M=\frac{1}{16\pi }\mathop{{\rm lim} }\limits_{r\to \infty }{{\oint }_{{{S}_{r}}}}{{\delta }^{mn}}({{\partial }_{n}}{{\gamma }_{mk}}-{{\partial }_{k}}{{\gamma }_{mn}}){{\hat{r}}^{k}}\;{\rm d}S,\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{{P}_{i}}=\frac{1}{8\pi }\mathop{{\rm lim} }\limits_{r\to \infty }{{\oint }_{{{S}_{r}}}}({{K}_{mi}}-{{\delta }_{mi}}K){{\hat{r}}^{m}}\;{\rm d}S,\end{eqnarray} \tag{ 43 }$$

where the components ${{\gamma }_{mn}}$, K_mn are in Cartesian coordinates, ${{\hat{r}}^{i}}={{x}^{i}}/r$ is the outgoing unit normal to the surface of integration, S_r denotes the 2 sphere of coordinate radius r and dS is the standard surface element of the 2 sphere. Under a more restricted class of gauge conditions (see [7] for details), one can also calculate the angular momentum of the spacetime from

$$\begin{eqnarray}&&{{J}_{i}}=\frac{1}{8\pi }\mathop{{\rm lim} }\limits_{r\to \infty }{{\oint }_{{{S}_{r}}}}({{K}_{jk}}-K{{\gamma }_{jk}})\xi _{(i)}^{j}\hat{r}k\;{\rm d}S,\end{eqnarray} \tag{ 44 }$$

where the ${{\xi }_{(i)}}$ are the Killing vectors associated with the asymptotic rotational symmetry. For an extended discussion of the ADM variables, the reader is referred to section 7 of [7] and references therein.

Horizons.

Event horizons are a defining criterion of BH spacetimes and mark the boundary between points from which null geodesics can reach infinity and points from which they cannot. Even though numerical algorithms have been developed for the calculation of event horizons in BH spacetimes [106–108], it is often more convenient to instead calculate the AH [57]. An AH is defined as the outermost marginally trapped surface on a spatial slice ${{\Sigma }_{t}}$. This condition can be shown to result in an elliptic equation for the outgoing normal direction sⁱ of the two-dimensional AH (see e.g. [109])

$$\begin{eqnarray}&&{{q}^{mn}}{{D}_{m}}{{s}_{n}}-K+{{K}_{mn}}{{s}^{m}}{{s}^{n}}=0,\end{eqnarray} \tag{ 45 }$$

where q_mn is the induced two-metric on the horizon surface. Unlike an event horizon, the AH can be calculated independently from the data on each hypersurface without further knowledge of the spacetime. Under the assumption of cosmic censorship and certain energy conditions, it can furthermore be shown that if a hypersurface ${{\Sigma }_{t}}$ contains an AH, it will coincide or lie within the event horizon's cross section with ${{\Sigma }_{t}}$ [110, 111]. The irreducible mass of a BH is directly encoded in the AH surface area by $M_{{\rm irr}}^{2}=16\pi {{A}_{{\rm AH}}}$ and, combined with the BH spin S gives the total BH mass through ${{M}^{2}}=M_{{\rm irr}}^{2}+{{S}^{2}}/(4M_{{\rm irr}}^{2})$ [112]. This formula also provides an estimate for the dimensionless BH spin $j=S/{{M}^{2}}$ in terms of the AH area and equatorial circumference $2\pi {{A}_{{\rm AH}}}/C_{e}^{2}=1+\sqrt{1-{{j}^{2}}}$; see e.g. [113]. The importance of horizons in the analysis of BH spacetimes follows to a large extent from the isolated and dynamic horizon framework developed by Ashtekar and coworkers [114].

Gravitational waves.

Arguably the most fundamental difference between the Newtonian and the general relativistic two-body problem is the dissipative character of the latter; energy and momenta of the binary are not conserved but partly radiated away in the form of GWs. Large-scale international efforts are dedicated to directly detect GWs with ground-based laser interferometric detectors LIGO, VIRGO, GEO600, KAGRA, future space missions of LISA type or pulsar timing arrays [115–120] and the theoretical prediction of the expected GW signals from astrophysical sources has been one of the main motivations of NR.

The most common approach to calculate GWs in BH simulations is based on the Newman–Penrose formalism [121] where the ten independent components of the Weyl tensor are projected onto a tetrad consisting of one outgoing and one ingoing null vector $\ell$ and $k$ as well as two complex spatial null vectors $m$ and $\bar{m}$. The interpretation of the resulting five complex scalars ${{\Psi }_{n}}$, $n=0,\;\ldots ,\;4$ in terms of GWs is based on the work of Bondi, Sachs and Penrose and coworkers [9, 10, 122, 123] and application of this formalism in NR requires a careful choice of the tetrad. In particular, the tetrad must correspond to a Bondi frame which can be realized, for example, by choosing a so-called quasi-Kinnersley tetrad [124, 125], i.e. a tetrad that converges to the Kinnersley tetrad [126] as the spacetime approaches Petrov type D. A particularly convenient choice is realized in the so-called transverse frame where the outgoing gravitational radiation is encoded in one complex scalar (see e.g. [127])

$$\begin{eqnarray}&&{{\Psi }_{4}}=-{{C}_{\alpha \beta \gamma \delta }}{{\ell }^{\alpha }}{{\bar{m}}^{\beta }}{{\ell }^{\gamma }}{{\bar{m}}^{\delta }}\;.\end{eqnarray} \tag{ 46 }$$

Even though ${{\Psi }_{4}}$ is well defined at infinity only, it is in practice often extracted at large but finite radii and this procedure generates various potential systematic errors which are discussed in [128]. The effect of these ambiguities is sometimes mitigated by extrapolating results at different finite radii to infinity [129] which suggests that the errors thus obtained are of the order $\mathcal{O}(\%)$. The Newman–Penrose scalar ${{\Psi }_{4}}$ is commonly decomposed into multipoles ${{\psi }_{\operatorname{lm}}}$ according to

$$\begin{eqnarray}&&{{\Psi }_{4}}(t,\theta ,\phi )=\mathop{\sum }\limits_{l,m}{{\psi }_{\operatorname{lm}}}(t)Y_{\operatorname{lm}}^{-2}(\theta ,\phi ),\end{eqnarray} \tag{ 47 }$$

where $Y_{\operatorname{lm}}^{-2}$ are spherical harmonics of spin weight −2. Often, the radiation is dominated by one or a few multipoles which can then be displayed in the form of functions of time. The Newman–Penrose scalar furthermore provides the energy, linear and angular momentum carried by the GWs which are obtained from straightforward integrals of ${{\Psi }_{4}}$ and its projections onto asymptotic Killing vectors [130].

Other methods for estimating the gravitational radiation have been applied in NR. (i) Perturbative wave extraction is based on the Regge–Wheeler–Zerilli–Moncrief formalism [131–133] and constructs master functions from the deviation of the spacetime metric from a Schwarzschild background; see e.g. [134, 135]. (ii) The Landau–Lifshitz pseudo tensor [136] is constructed by mapping the curved, physical spacetime onto an auxiliary flat spacetime with metric ${{\eta }_{\mu \nu }}$. This leads to expressions for the radiated energy and momenta; for applications see e.g. [137]. (iii) In characteristic formulations of the Einstein equations, the Bondi news function [9] provides a direct measure of the GW signal which has also been used in 3+1 NR through Cauchy-characteristic extraction [138, 139]. This method provides a particularly accurate extraction since it is performed by construction at infinity; for a comparison with other methods see [140].

In this section we will briefly review the main developments in NR leading to the breakthroughs of 2005. It is beyond the scope of this work to present a comprehensive history of the enormous amount of work and publications generated in this field over the last 50 years. Our review should therefore be understood as a potentially biased precis intended to give the reader a rough guideline of NR's history. The articles quoted in this section contain many further references the reader will find a valuable source for a more thorough account.

The earliest documented effort to generate BH spacetimes by numerically solving Einstein's equations was made half a century ago by Hahn and Lindquist [141]. It is worth noting that at the time the notion of BHs, horizons and the area theorems were not yet understood. Furthermore, virtually nothing about the delicacies of all the ingredients discussed in the previous section was known at the time. It is thus not too surprising that they could evolve their data for very short times only. And yet, their first steps into uncharted territory demonstrated a genuinely new alternative for the modeling of BHs even if few, at the time, would have predicted what kind of avalanche of numerical explorations had been kicked loose.

Starting in the late 1960s, the problem was reinvestigated in an effort initiated by DeWitt which led to PhD theses by Čadež, Smarr and Eppley [142–144]. These works implemented the ADM equations in axisymmetry using a specific type of coordinates often referred to as Čadež coordinates and thus evolved head-on collisions starting with Misner data, testing several gauge conditions including maximal slicing, vanishing shift and minimal distortion shift. Their equal-mass head-on collisions predict a GW energy of about $0.1\%$ of the total mass, albeit with uncertainties of a factor a few; for details see [36, 145, 146]. This value turned out to be correct within a factor of about 2.

The next burst of efforts took place in the 1990s, much of it as part of the 'Binary Black Hole Grand Challenge Project' [147]. Earlier studies of this decade still employed axisymmetry and similar techniques as above, but on significantly improved computational architecture. Anninos et al [148–150] extracted the l = 2 and l = 4 multipoles of the emitted GW signal, calculated AHs and compared results obtained with different wave extraction methods; they confirmed the earlier estimates for the emitted GW energy within error bars and found good agreement with close-limit [151] predictions for small initial separations of the BHs; see also [152] for collisions of boosted BHs. One of the most memorable results of these axisymmetric BH simulations is the 'pair-of-pants' picture obtained when calculating the BH horizons in a merger process; see figure 10 in [153]. Using a special class of 'body fitting' coordinates, Anninos and Brandt [154] performed the first evolutions of unequal-mass BH head-on collisions. By extracting GW modes up to l = 4, they validated perturbative results in the close and large separation limit. Further axisymmetric studies include initially distorted, rotating BHs and the resulting GW emission [155] and accretion onto rotating BHs [156].

The first fully 3+1 dimensional BH simulations were presented in 1995 by Anninos et al [53] with the so-called 'G-code'. This code is based on the ADM formulation, uses Schwarzschild initial data in isotropic coordinates, different types of singularity avoiding slicings and a shift that locks the coordinate radius to the AH location. It produced numerically stable solutions of a Schwarzschild BH for up to $t\approx 50\;M$. The GW signal including BH ringdown obtained with the G-code were found to agree well with those from axisymmetric codes [157, 158]. Further development of their 3+1 code enabled the Grand Challenge Alliance to evolve a single BH that moves across the computational domain with 0.1 times the speed of light for a total time of about $60\;M$ [159]. Around the same time, mesh refinement was first used in 3+1 simulations of a BH by Brügmann [72]. The year 1997 saw the release of Cactus 1.0 [160], a freely available environment for the development of parallel, scalable, high-performance multidimensional component-based code for NR and other numerical applications.

The first binary BH merger in 3+1 NR was simulated by Brügmann [161] in grazing collisions using the ADM equations and the fixed puncture technique. The first grazing collisions of BHs using excision were performed with Agave, a revised version of the Grand Challenge code [162].

As the second millennium drew to a close, however, it was still a general feature of the space–-time-split based codes, in axisymmetry or full 3+1, to be limited by numerical instabilities to life times of the order of $\mathcal{O}(100\;M)$. This is in sharp contrast to the remarkable stability properties achieved at the same time using characteristic methods which facilitated long-term stable simulations of single distorted, moving or rotating BHs with lifetimes up to $60\;000\;M$ [12, 13]. As mentioned in section 3, characteristic codes have not yet been generalized successfully to binary BH spacetimes where the formation of caustics and the ensuing breakdown of the null coordinates have so far represented an insurmountable obstacle. A different approach named 'Lazarus' [75, 163] was developed around the turn of the millennium in order to maximize the scientific output obtained from 3+1 evolutions as available at the time. In this eclectic approach the relatively short numerical simulations are matched to perturbative calculations once a merger into a single BH has occurred and the spacetime is perturbatively close to a Kerr BH. The preceding inspiral phase, instead, is to be described by PN methods which provide initial data for the numerical computation; for applications of this method see [164, 165].

Progress in the stability properties of 3+1 codes accelerated considerably in the early 2000s as a wider range of formulations of the Einstein equations and gauge conditions were used in BH simulations. The first applications of the BSSN formulation focused on GW pulses, including collapse to BHs, and demonstrated significantly better stability properties than the ADM system [23, 24, 166]. Soon afterwards, this observation was confirmed for evolutions of boson stars, neutron stars and BHs [167]. Using the BSSN formulation combined with a 'K-freezing' slicing and Gamma-freezing shift (see their section 3 for details) and a 'simple excision' of a cubic region on whose boundary time derivatives are copied from neighbouring grid points, Alcubierre and Brügmann [26] were able to evolve a Schwarzschild BH in ingoing Eddington Finkelstein coordinates over many thousands of M with no signs of instability. These simulations were obtained in 3+1 dimensions with octant symmetry and the stability properties could be generalized to 3+1 grids with no symmetry by using the constraint ${{\mathcal{G}}^{i}}$ of equation (28) in the evolution equation (27) [25, 40]. Evolutions of distorted BHs with the BSSN system were pushed to a few $100\;M$, about twice as long as axisymmetric ADM simulations [42]. By evolving Brill–Lindquist data with BSSN, 1+log slicing and variants of the Gamma-freezing shift, combined with the fixed puncture technique (see section 4.3), Alcubierre et al [40] extracted GWs from BH collisions in good agreement with earlier axisymmetric ADM simulations but over much extended lifetimes of $\sim 1000\;M$. BSSN simulations of head-on collisions were extended to include mesh refinement in [135, 168] using the Carpet [79, 80] and Paramesh [74] packages, respectively. The first evolution of a quasi-circular BH binary extending over more than one orbital time scale was performed by Brügmann et al [54] in 2003 and explored in more detail in [169]. Around the same time BSSN simulations of inspiraling and merging neutron-star binaries were obtained by various groups [170–172]. While neutron star spacetimes contain complex matter sources, the spacetime curvature is significantly weaker than for BH binaries and there are no singularities (other than nearly stationary post-merger BHs). Possibly, therein lie the reasons why neutron star inspiral and merger simulations were achieved before their binary BH counterparts. Early 3+1 applications of the GHG formulation focused on the collapse of scalar fields and the approach to the formation of singularities using unigrids as well as mesh refinement [33, 35].

In 2005, the jigsaw was finally assembled. The first simulations of BH binaries through inspiral, merger and ringdown were obtained by Pretorius [1] and a few months later by the Brownsville/Rochester and NASA Goddard groups [2, 3]. The new ingredients which finally pushed the BH simulations 'over the cliff' can probably be summarized as follows. The GHG formulation was adjusted by the addition of constraint damping terms suggested by Gundlach et al [34]; see equation (32). Furthermore, Pretorius managed to specify conditions for the gauge functions ${{H}^{\alpha }}$ that avoided instabilities which had troubled earlier simulations using the GHG system; see equation (36). His compactification of the spacetime and the damping of signals near spatial infinity through numerical viscosity prevented outer boundary effects to affect the BH region in a significant manner. The main new feature of the moving puncture simulations is encapsulated in the word 'moving'. Up to that point, puncture simulations had kept the location of the singularity fixed and factored out the irregular part of the metric analytically. By using gauge conditions as given in equations (33)–(35), the Brownsville/Rochester and NASA Goddard groups obtained a description where the punctures could freely move across the grid. They furthermore managed to apply standard finite differencing across the puncture which effectively provided an exceptionally robust type of singularity excision; see the discussion in section 4.3.

These breakthroughs in simulating BH binaries triggered a veritable phase transition in the field of NR as the community gained unprecedented insight into the dynamics of BH binaries. This is the subject of our next section.

In spite of its dissipative character, the two-body problem in GR leads to classes of orbits similar, though not identical, to their Newtonian counterparts discussed at the end of section 2. The simplest configuration, and focus of the breakthroughs as well as early follow-up studies, is the quasi-circular inspiral of two non-spinning BHs of equal mass. The quasi-circular case is also of particular relevance for the modeling of GW sources for ground based detectors because GW emission has long since been known to efficiently decrease the orbital eccentricity even at large binary separation [173], so that BH binaries are expected to be circularized to high precision by the time they reach the frequency window of the detectors. The dynamics of the late stages of the quasi-circular inspiral and merger of an equal-mass, non-spinning BH binary are graphically illustrated in figure 5. The left panel shows the trajectory of one BH with that of the other obtained through reflection at the origin. The binary separation decreases through the emission of GWs which is dominated by the quadrupole shown in the form of the GW strain h in the right panel. The strain is a complex function containing the two GW polarization modes and related to the Newman–Penrose scalar ${{\Psi }_{4}}$ from which it is calculated by two integrations in time, and it is the quantity of choice used in GW data analysis; see e.g. [174, 175]. As the binary inspirals, the wave signal increases in amplitude and frequency which reach a maximum around the merger stage followed by the exponentially damped ringdown signal; for details see e.g. [176, 177]. Eventually, the spacetime settles down into a stationary Kerr configuration and GW emission ceases. Realistic astrophysical BH binaries start at a much larger separation than shown in figure 5 corresponding to a much longer inspiral phase of many thousands of cycles [129]. The breakthrough simulations studied short inspirals of about one orbit and found that about $3\%$ of the total mass is radiated in GWs leaving behind a spinning BH with dimensionless spin parameter of $\sim 0.7$. These results have been confirmed with good precision by many follow-up studies, but the earlier inspiral phase not covered in the rather short simulations provides additional energy release in GWs. The most accurate simulation to date is that of Scheel et al [178] and gives ${{E}_{{\rm GW}}}/M=0.048\;38\pm 0.000\;02$ for a 16 orbit inspiral and a dimensionless spin 0.686 46 ± 0.000 04 for the resulting Kerr BH.

Zoom In Zoom Out Reset image size

For the reasons mentioned above, eccentric binaries have received less attention in the NR studies so far. In fact, a good deal of effort has been spent on measuring the eccentricity in binaries and reducing the residual eccentricity in the initial data as much as possible in order to generate high-precision GW templates for quasi-circular BH systems [179–182]. Still, binaries with sizable eccentricity have been studied in their own right [183, 184] and also compared with PN predictions [185]. Due to the GW emission, the eccentricity is a time dependent function in GR; see equation (5.13) in [173]. For high values of , the BHs do not complete a single orbit but rather plunge into each other. For mild eccentricities, the periastron advance obtained numerically shows some deviations from post-Newtonian predictions [181] but excellent agreement with perturbative calculations [186]. An interesting feature has been found for eccentric binaries when the initial momenta of the BHs are fine-tuned, a threshold of immediate merger where the binary spends some time in near-circular orbits before eventually merging or separating [187]; see also [188, 189]. In this regime, also identified in geodesic calculations [190, 191], the number of zoom-whirl orbits exhibits a logarithmic dependence on the deviation of the impact parameter b = L/P (L is the initial orbital angular momentum of the binary and P the momentum of either BH in the centre-of-mass frame) from the threshold of immediate merger b* [113, 187]. For a fixed value of the momentum parameter P, numerical studies have revealed three possible regimes separated by two special values of the impact parameter, the threshold of immediate merger b* and the scattering threshold ${{b}_{{\rm scat}}}$: (i) prompt mergers resulting from a plunge or inspiral for $b\lt b*$, (ii) non-prompt mergers where the BHs experience a close encounter, then separate but loose enough energy in GWs to eventually form a bound system and merge for $b*\lt b\lt {{b}_{{\rm scat}}}$, and (iii) scattering configurations for $b\gt {{v}_{{\rm scat}}}$, where the BHs separate to infinity [113, 192]. Here, the latter two regimes can only be obtained for sufficiently large P. We note the similarity between configurations with ${{b}_{{\rm scat}}}$ and parabolic orbits in Newtonian gravity which forms the boundary between bound and unbound configurations. Furthermore, BH collisions at velocities close to the speed of light can generate enormous amounts of gravitational radiation of up to $\sim 50\%$ of the total centre-of-mass energy [70] making them ideal scenarios to probe GR in its most violent regime.

Unequal mass-ratios $q={{m}_{2}}/{{m}_{1}}$ (here defined such that $q\leqslant 1$) and/or non-zero BH spins naturally affect the dynamics of the BHs interaction, but as yet there are no indications that the above picture of the different types of orbits is changed in a qualitative manner. For many practical applications, however, spins and unequal masses have vital implications. Unequal mass ratios lead to asymmetric emission of GWs which imparts a recoil or kick on the post-merger remnant BH of up to 175 km s⁻¹ [193–195] and specific spin-orientations can lead to so-called superkicks of several 1000 km s⁻¹ [196–199]. Spins aligned with the orbital angular momentum cause a so-called hang-up effect extending the late inspiral stage and leading to an increase in GW emission by about a factor of two compared with the non-spinning case [64, 200]. Signatures of spins and mass ratio furthermore leave specific imprints such as amplitude modulation or relatively enhanced higher-order multipoles, which are important effects in the analysis of observational data of GW detectors and, hence, require extensive modeling combining NR methods with post-Newtonian and/or perturbative techniques; see [129, 201–204] and references therein.

It is well beyond the scope of this article to discuss all these features and the wide range of applications of NR simulations of BH binaries in detail. Instead, we will refer the interested reader in the next section to various reviews that have appeared over the past decade and provide extended overviews of the many, exciting and sometimes unexpected developments that have been triggered by the 2005 breakthroughs in the modeling of BH binaries in GR.

The 2005 breakthroughs in the modeling of binary BHs mark a phase transition in the field of NR. Even though some very specific BH systems (e.g. head-on collisions) had been modeled successfully before, these were few isolated and idealized examples. In 2005, virtually the entire space of BH binary systems was opened up for accurate, quantitative modeling. The years following 2005 have sometimes been referred to as the 'goldrush years of NR' as the community suddenly had available the tools to study a wealth of phenomena previously subject to speculations rather than precision studies. BH kicks, for example, had been known for decades to result from GW emission with potentially dramatic astrophysical consequences but the magnitude of this effect remained largely shrouded in mystery. The superkicks of thousands of km s⁻¹ have been one of the most remarkable and unexpected results of post-2005 NR, a fact still in the process of digestion in the interpretation of astrophysical observations [205].

A lot of motivation for the long-term effort of NR came from the modeling of GW sources in support of the ongoing search for direct detection with detectors such as LIGO, VIRGO, GEO600 or KAGRA. Here the impact of the NR breakthroughs is both short and long-term. A wealth of numerical studies has been devoted to exploring the parameter space of BH binaries in the attempt to construct waveform catalogues for use in GW data analysis. The parameter space, however, has at least seven dimensions (mass ratio and three spin parameters for each BH) making a dense coverage in every dimension prohibitively costly from a computational point of view. Instead, numerical as well as analytic studies have looked for systematic dependencies of the GW signals on some of the parameters and thus effectively reduce the parameter space to be explored numerically; see e.g. [206–208]. At the same time, the codes have matured considerably in accuracy and efficiency and are now capable of generating large numbers of waveforms [129, 209, 210]. And yet, a comprehensive modeling of the BH parameter space certainly cannot be done exclusively with NR because the inspiral signal of a BH binary in the sensitivity band of ground (or space based) detectors typically contains many thousands of orbits rather than the $\mathcal{O}(10)$ orbits covered in numerical simulations. Waveform templates are therefore constructed by stitching together post-Newtonian with numerical waveforms as for example in [201] or using effective one body models [211, 212] and calibrating free parameters by comparison with NR predictions; see e.g. [203, 213]. Furthermore, an important class of sources for space based laser interferometers of LISA type are so-called extreme-mass-ratio inspirals with mass ratios down to $\mathcal{O}({{10}^{-6}})$ which the current NR codes cannot handle; the smallest mass ratios achieved to date is $q=1/100$ for a small number of orbits [214]. These as well as intermediate mass ratio inspirals with $q=\mathcal{O}({{10}^{-3}})$ require modeling through perturbative techniques; see [215, 216] and references therein.

In recent years, BH binary codes have been extended to include various forms of matter with particular focus on the identification of electromagnetic counterparts to GW signals. As an example, we note the variation of the Blandford–Znajek effect in the inspiral of two BHs in the presence of a circumbinary disk whose magnetic field is capable of extracting rotational energy from the orbital motion of the binary giving rise to single or dual jets [217, 218]. Other studies address accretion of matter onto BHs, the impact of a recoiling BH on the surrounding disk material, periodic oscillations in the matter due to the binary motion or the tidal disruption of white dwarfs in the field of a BH [219–222]; for an overview see also [223]. Non-vacuum spacetimes also include some classes of alternative theories of gravity. NR simulations of BH binaries have so far focused on scalar-tensor theories of gravity which are conveniently implemented in existing GR codes by merely adding a minimally coupled scalar field in the so-called Einstein frame [224, 225]. Extension of these studies to a wider class of alternative theories of gravity may be a highly non-trivial task as the hyperbolicity properties of the underlying evolution equations and, thus, their suitability for numerical treatment remain unclear for most theories; for such an investigation into dynamical Chern–Simons theory see [226].

One of the most remarkable developments of NR in the last ∼10 years is the wide range of applications in areas of physics far outside the more traditional regime of GW and astrophysics. High-energy collisions of BHs may provide valuable inside into the cross section of proton–proton collisions performed at the Large Hadron Collider to probe the possibility of BH formation as conjectured [227–229] in the so-called TeV gravity scenarios [230–233]. The gauge-gravity duality, also often referred to as the anti-de Sitter (AdS)/conformal field theory correspondence [234–236] relates BH spacetimes in asymptotically AdS spacetimes to physical systems in the strongly coupled regime of certain gauge theories. To name but a few examples, this duality has been used to model the thermalization of balls of deconfined plasma in heavy-ion collisions (see [237–241] and references therein), the optical conductivity of so-called strange metals [242, 243] or condensed matter systems (see [244] for an overview). Numerical studies of BHs and their formation still teaches us unexpected lessons about the fundamental properties of GR as for example in Choptuik's critical collapse study [71] or the surprising instability observed in the extension of Choptuik's work to asymptotically AdS spacetimes by Bizoń and Rostworowski [245]. Tests of the cosmic censorship conjecture are now possible over a wide range of physical scenarios including high-energy collisions of BHs [246], BH strings in higher dimensions [247] or BH collisions in asymptotically de Sitter spacetimes [248]. The numerical simulation of lattices consisting of multiple BHs are being used to model large scale structures in cosmological spacetimes [249, 250] and NR tools are developed for the modeling of the early universe [251].

This short (and incomplete) summary demonstrates that the breakthroughs of 2005 have opened the door onto a vast garden of opportunities probably well beyond the wildest expectations the NR community has harbored during the 40 year path towards the holy grail. We conclude here with a list of suggestions for further reading on the various topics whose surfaces have been scratched on the preceding pages.

Extended books on the methodology of NR have been written by Alcubierre [252], Bona et al [253] as well as Baumgarte and Shapiro [254]; we also note Gourgoulhon's review of the '3+1' formulation of GR [7] and the comprehensive review by Centrella et al [255]. An earlier review also covering the techniques for simulating matter is given by Lehner [256]. The field of GW physics and the use of BH binary simulations therein is discussed in various articles [257–260]. Comparisons of GW signals obtained numerically with those from various semi-analytic calculations are reviewed in [261]. The reader will find overviews of several applications including astrophysics and high-energy physics in [44, 262]. NR applications outside the more traditional areas of astrophysics and GW physics are the focus of [263] and in particular the review [264]. A description of more technical details of BH simulations in asymptotically AdS spacetimes is given in [240]; see also [265].

If the past ∼10 years of NR have shown anything it is an astounding potential to surprise us with unexpected results and new areas of applications. Aside from the above mentioned articles in the present literature, the reader will undoubtedly be richly rewarded by following online and print journals to remain up to date on the avalanche of results and applications of BH simulations triggered by the milestone of the 2005 breakthroughs.

The author thanks Emanuele Berti, Bernd Brügmann, Vitor Cardoso, Pau Figueras, Jose González, Leonardo Gualtieri, Mark Hannam, Carlos Herdeiro, David Hilditch, Sascha Husa, Bernard Kelly, Pablo Laguna, Luis Lehner, Christian Ott, Frans Pretorius, Harvey Reall, Christian Reisswig, Erik Schnetter, Deirdre Shoemaker, Ken Smith, Carlos Sopuerta, Helvi Witek, and Miguel Zilhão for many fruitful discussions. This work was supported by the FP7-PEOPLE-2011-CIG CBHEO grant no. 293412, the FP7-PEOPLE-2011-IRSES NRHEP grant no. 295189, the STFC grant nos. ST/I002006/1 and ST/L000636/1, the Trestles system of the San Diego Supercomputing Centre (SDSC), Stampede of the Texas Advanced Computing Center (TACC) and Kraken of the National Centre for Supercomputing Applications through XSEDE grant no. PHY-090003 by the National Science Foundation, the COSMOS Shared Memory system at DAMTP, University of Cambridge, operated on behalf of the DiRAC HPC Facility and funded by BIS National E-infrastructure capital grant no. ST/J005673/1 and STFC grant nos. ST/H008586/1 and ST/K00333X/1, the Centro de Supercomputación de Galicia (CESGA) under grant no. ICTS-2013-249, and the European Union's FP7 ERC Starting Grant DyBHo-256667.