Data Science Analytics and method calculating black-hole properties from gravitational-wave data

This abstract delves into the realm of Data Science Analytics applied to the computation of black-hole properties using gravitational-wave data. The research revolves around harnessing cutting-edge analytical methods to derive valuable insights from gravitational-wave signals emitted during black-hole mergers. These astrophysical events hold great significance, acting as crucial sources for gravitational-wave detectors. Through the integration of numerical simulations and innovative data analysis techniques, researchers aim to unravel essential black-hole properties, leading to a profound understanding of the dynamics inherent in these cosmic phenomena. This research carries immense importance in advancing our knowledge of astrophysics and gravitational-wave astronomy.

Dr Francesco Dergano
64 min readFeb 1, 2020


The longstanding, unsolved problem in theoretical physics lies in comprehending the predictions of general relativity concerning the interactions between two black holes. Black-hole mergers are profound events in astrophysics, releasing tremendous energy as gravitational radiation and playing a pivotal role as sources for gravitational-wave detectors. The dynamics and gravitational waveforms resulting from these mergers can only be calculated through numerical simulations of Einstein’s equations. For years, numerical relativists struggled with their codes crashing after simulating only a fraction of a binary orbit. However, recent breakthroughs in numerical relativity have enabled stable and robust black-hole merger simulations, leading to a profound understanding of black-hole binary dynamics. This remarkable progress has far-reaching implications for astrophysics, gravitational-wave astronomy, and other scientific domains.


In the preamble, we delve into the tremendous energy released during the final merger of two black holes in a binary system, emitting gravitational waves that traverse the Universe at the speed of light, carrying the signature waveform of the merger. The existing ground-based gravitational-wave detectors are well-equipped to capture stellar black-hole binary mergers, remnants of massive stars. Moreover, plans are underway for a space-based detector to observe mergers of colossal black holes at the centers of galaxies, boasting masses ranging from approximately 10^4 to 10^9 times that of the Sun.

The calculation of these waveforms requires solving the Einstein equations of general relativity on a computer, which is an immensely challenging task. Previous attempts at numerical relativity encountered various puzzling instabilities, leading to code crashes before successfully simulating a significant portion of a binary orbit. However, recent breakthroughs in numerical relativity have revolutionized the field, enabling robust and accurate simulations of black-hole mergers.

The article reviews these groundbreaking achievements and the valuable insights gained about black-hole mergers, emphasizing their applications in astrophysics and gravitational-wave data analysis. Our focus lies in comparable-mass black-hole binaries, where the component mass ratios range from 1 to 10. We refer to the symmetric mass ratio, represented as q = (Mi * M2) / (Mi + M2), which plays a significant role in our discussions.

For simplicity, we adopt the convention of setting the speed of light (c) and gravitational constant (G) to unity. With this normalization, we can scale the dynamics and waveforms for black-hole binaries based on the total system mass (M). This allows us to express length and time scales in terms of the mass, facilitating our analysis with M ~ 5 × 10^(-5) M_⊙ / M and M ~ 1.5M_⊙ / km.

We commence by providing the scientific and historical contexts. In Section II, we offer a brief overview of astrophysical black-hole binaries as crucial sources for gravitational-wave detectors. Moving on to Section III, we trace the historical progression of evolving black-hole mergers on computers over more than four decades, culminating in recent breakthroughs.

With the stage set, we delve deeper into the key components essential for successful black-hole merger simulations in Section IV. This includes an exploration of computational methodologies, such as numerical-relativity techniques and black-hole binary initial data. The heart of this review lies in Section V, where we extensively discuss the significant findings from numerical relativity simulations of black-hole mergers. We follow a historical development, focusing on the dynamics of mergers and the resulting gravitational waveforms, which have opened up exciting applications in general relativity, gravitational waves, and astrophysics.

In Section VI, we explore the synergistic interactions between numerical relativity and analytic approaches to model gravitational dynamics and waveforms. Section VII delves into the applications of these results to gravitational-wave data analysis. Moreover, in Section VIII, we discuss the impact of merger simulations on astrophysics, touching upon recoiling black holes and potential electromagnetic signatures of the final merger. Our journey concludes in Section IX, where we look at the frontiers and future directions of this field.

Prior to delving into the main content, we provide references to additional resources that may interest our readers. These include review articles by Lehner (2001) and Baumgarte and Shapiro (2003) on numerical relativity before the recent breakthroughs in black-hole merger simulations. Pretorius (2009) offers an early review of the recent successes, covering topics similar to those discussed here. Hannam (2009) reviews the status of black-hole simulations producing long waveforms and their relevance to gravitational-wave data analysis. For further mathematical and computational aspects of numerical relativity, readers can refer to the textbooks by Bona and Palenzuela (2005) and Alcubierre (2008), which complement our discussions.


Black holes and gravitational waves are fascinating predictions in physics, and they come together in black-hole binaries, which are among the strongest sources of gravitational radiation.

A Basic Properties

In the section titled “Basic Properties,” we present fundamental characteristics of black holes and gravitational waves. For more in-depth discussions and details, readers can refer to Miner et al. (1973) and Schutz (2009).

  1. Black-Hole Basics: A black hole is formed when matter collapses to infinite density, creating a singularity of infinite curvature in spacetime. Each black hole is surrounded by an event horizon, where the escape velocity equals the speed of light. The event horizon is defined by the paths of photons marking the boundary between trajectories that fall inward and those that escape to infinity. These photons hover at a finite radius at the black hole’s surface. Knowing the entire future development of the system is necessary to locate the event horizon, as mass or energy falling into it can change the hovering location of these photons. When black holes merge, a single event horizon forms with an area equal to or larger than the sum of the individual horizons. To determine when this occurs during calculations, numerical relativists use an apparent horizon, which relies solely on spacetime properties at a given time (Poisson, 2004). For quiescent black holes, the apparent and event horizons coincide, but for more general cases, the apparent horizon is always inside the event horizon, ensuring no physical phenomenon inside the apparent horizon affects the spacetime outside. The simplest black hole, a non-rotating one, is described by the spherically symmetric Schwarzschild solution. It is fully characterized by its mass M, and its horizon is located at r = 2M (in Schwarzschild coordinates), with an area of 4π(2M)^2. More general black holes can have charge and spin. In astrophysics, charged black holes are rapidly neutralized by surrounding plasma, so rotating, uncharged black holes are considered. Stationary black holes are described by the axisymmetric Kerr solution, characterized by mass M and angular momentum per unit mass a. The event horizon is at the Boyer-Lindquist radius r+, given by r+ = M + (M^2 – a^2)^1/2. When a = M, the black hole is “extremal” or maximally rotating. As the black-hole spin increases, the horizon moves deeper into the potential well. Near a single black hole, photons and test particles can experience either stable or unstable orbits. The innermost stable circular orbit (ISCO) for massive test particles occurs at r = 6M for a Schwarzschild black hole of mass M. For a Kerr black hole, the ISCO is closer in for co-rotating test particles and farther out for counter-rotating particles. While strictly defined for massive test particles, the concept of ISCO has proven useful in studying spacetime around two black holes spiraling together on quasicircular orbits. It helps identify the separation at which angular momentum is minimized, indicating the black holes’ motion toward each other, even without gravitational radiation.
  2. The Gravitational Wave Primer section introduces gravitational waves as ripples in spacetime, carrying energy and momentum while traveling at the speed of light, indicating disturbances in the gravitational field. Like electromagnetic waves, gravitational waves can be broken down into multipolar contributions that reflect the source generating them. Just as electromagnetic radiation lacks a monopole contribution due to charge conservation, gravitational waves lack monopole radiation due to mass-energy conservation. While electromagnetic waves can have dipole character from variations in charge and currents, dipolar gravitational waves are impossible due to linear and angular momentum conservation. Consequently, the leading-order contribution to gravitational radiation is quadrupolar. Gravitational waves are produced by systems with time-varying mass quadrupole moments (Flanagan and Hughes, 2005; Misner et al., 1973). In the wave zone, a gravitational wave is described as a perturbation (h) to a smooth underlying spacetime, and its amplitude is proportional to the source’s quadrupole moment (Q), mass (M), and distance from the source (r). The strongest gravitational waves originate from large masses moving at high velocities, such as binary systems of compact stars and black holes. Gravitational waves have purely transverse effects, acting tidally in directions perpendicular to their propagation direction. When gravitational waves interact with a detector of length L, they produce a length change SL/L ~ h/2. Typical astrophysical sources yield wave amplitudes of h < 10^-21 at Earth, making precise measurements essential for detections. Gravitational waves have two polarization components, h+ and hx, for a linearly polarized wave. Figure I illustrates the corresponding lines of force for sinusoidal gravitational waves propagating along the z-axis. In the purely + polarization, the wave stretches along one axis and squeezes along the other, alternating sinusoidally as it passes. The × polarization acts similarly but stretches and squeezes along axes rotated by 45°. In general, a gravitational wave is a superposition of these two states, expressed as a complex waveform strain h, where h = h+ +ihx.
FIG. 1 Lines of force for plane gravitational waves propagat-
ing along the z axis. The wave on the left is purely in the
+ polarization state, and the one on the right is purely in
the x polarization state. The gravitational waves produce
tidal forces in planes transverse to the propagation direction.
From (Abramovici et al., 1992). Reprinted with permission from AAAS.

B. Astrophysical Black Holes

Astronomers have made remarkable discoveries of black holes at various scales throughout the Universe. Among these, stellar black holes stand out as the smallest, with masses ranging from approximately 3 to 30 times that of our Sun (Mo). They result from the evolution of massive stars and are supported by dynamical measurements of compact objects in transient systems experiencing X-ray outbursts. These observations validate that objects with masses greater than about 3 Mo must be black holes since neutron stars cannot reach such high masses (Remillard and McClintock, 2006).

Moving up the scale, intermediate-mass black holes (IMBHs) weigh approximately 102 to 101 times that of our Sun. IMBHs may emerge from multiple mergers of smaller objects within dense stellar clusters in the present Universe or from the evolution of exceptionally massive stars, forming black-hole “seeds” at the centers of massive halos. Currently, the most compelling observational evidence for IMBHs arises from models of ultra-luminous X-ray sources (Colbert and Miller, 2005).

At the largest scale, massive black holes (MBHs) weigh approximately 104 to 109 Mo and are found at the centers of galaxies, including our Milky Way. The existence of MBHs is strongly supported by dynamical models of stars and gas moving in the gravitational potential well of the central MBH.

This article concentrates on black-hole binaries, where each component is a black hole. We focus on comparable-mass binaries, as they are expected to produce robust gravitational-wave signals. Stellar black-hole binaries may form through binaries composed of two massive stars or from dynamic processes within dense stellar environments. Similarly, IMBH binaries can arise through dynamical processes in stellar clusters and through mergers of massive halos at high redshifts. While observing these black-hole binaries through electromagnetic radiation presents challenges, the detection of gravitational radiation from these systems will greatly enhance our understanding of their dynamics (Bulik and Belezynski, 2009; Miller, 2009).

MBH binaries emerge when galaxies merge, considering almost all galaxies contain an MBH at their center and have likely experienced at least one merger throughout the Universe’s history. However, the detection of MBH binaries through electromagnetic observations is currently limited due to vast cosmic distances and small angular separations on the sky. Initially, MBH binaries have wide separations, and the gravitational radiation they emit is weak, insufficient to cause coalescence within the age of the Universe. Nevertheless, various processes, such as gaseous dissipation and interactions with stars, can gradually reduce their orbital energy, eventually leading to their coalescence and the emission of powerful gravitational waves (Sesana et al., 2009).

C. Gravitational Waves from Black-Hole Binaries

Gravitational waves emitted by merging black-hole binaries, which go through three stages – inspiral, merger, and ringdown, are expected to be powerful sources of these waves. In the inspiral stage, the black holes gradually approach each other on nearly circular orbits, emitting gravitational waves in a chirp-like pattern with increasing frequency and amplitude. In the merger stage, the black holes plunge together, forming a single, highly distorted remnant black hole, and this phase requires numerical relativity simulations. Finally, the remnant black hole settles down into a quiescent state through gravitational-wave emission, a process known as ringdown. These gravitational waves cover different frequency regimes and will be detected by ground-based and space-based laser interferometer detectors, as well as pulsar timing arrays.


Over the course of more than four decades, researchers have been engaged in a quest to compute the gravitational-wave signals resulting from the merger of black holes. This journey has witnessed significant advancements in theoretical and experimental general relativity, astrophysics, and computational science. In this section, we will provide an overview of these developments, followed by a closer examination of specific milestones that have contributed to the successful simulation of black-hole mergers.

A Setting the Stage

In the late 18th century, Michell and Laplace speculated that a star could become so compact that its escape velocity would surpass the speed of light, using Newtonian gravity. Later in the 20th century, scientists realized that such objects, known as black holes, could form as a result of total gravitational collapse in the framework of general relativity. John Wheeler popularized the term “black hole” for such entities.

From the 1960s, various highly energetic astrophysical phenomena were discovered, indicating strong gravitational fields as the underlying mechanisms, including quasars and X-ray binaries like Cygnus X-1, which was the first credible black-hole candidate. Astrophysical black holes are now believed to exist on different scales throughout the Universe, with black-hole binaries being considered significant sources of gravitational waves.

Gravitational waves were first recognized as solutions to the linearized, weak-field Einstein equations in the early 20th century. Over time, gravitational waves were acknowledged as real physical phenomena, carrying energy and capable of producing responses in detectors. The detection of two neutron stars in a binary system provided an indirect confirmation of gravitational radiation, as their observed orbit matched the predictions of general relativity. This led to the Nobel Prize being awarded to Hulse and Taylor in 1993. The development of numerical relativity was facilitated by advancements in computer power and computational methods, allowing researchers to explore physics beyond the scope of perturbation theory.

Numerical-relativity simulations generally involve decomposing four-dimensional spacetime into three-dimensional spacelike slices connected by timelike curves, using the “3+1" approach pioneered by Arnowit, Deser, and Misner (ADM). This Hamiltonian formulation involves configuration variables represented by the three-metric on spatial slices and conjugate momenta based on the extrinsic curvature. Solving the Cauchy problem with appropriate gauge conditions is essential for the successful simulation of black-hole mergers.

B. Numerical Relativity Milestones

Numerical relativity milestones related to simulating black-hole mergers have been achieved over several decades. In 1964, Hahn and Lindquist made the first attempt at simulating the head-on collision of two equal-mass black holes using a two-dimensional axisymmetric approach. Smarr and Eppley revisited the problem in the mid-1970s with improvements in coordinate conditions, but instability and numerical errors posed challenges.

In the 1990s, with the LIGO project advancing and interest in observing gravitational waves from black-hole mergers, numerical relativists reexamined the head-on collision problem with modern techniques and more powerful computers. In parallel, the National Science Foundation funded a large collaboration aimed at evolving black-hole mergers in three dimensions and calculating gravitational waveforms. Around the same time, a significant numerical relativity group emerged at the Albert-Einstein Institut in Germany.

Throughout the late 1990s and early 2000s, advancements were made in various aspects of numerical relativity. These milestones include developing initial data for binary black holes, new methods to represent black holes on computational grids, recognizing the importance of hyperbolicity in formulating Einstein’s equations for numerical solution, improved formulations of the Einstein equations, and the use of fully three-dimensional evolution codes to study various scenarios like distorted black holes and collisions.

Over time, the length of black-hole evolutions gradually increased, and novel approaches, such as the Lazarus project, combined relativistic short-duration binary simulations with perturbation techniques for late-time behavior. By late 2003, the first complete orbit of two equal-mass, non-spinning black holes was achieved. Progress was challenging and incremental but set the stage for significant advancements in the field of numerical relativity.

C. Breakthroughs and the Gold Rush

In early 2005, Frans Pretorius achieved a groundbreaking feat by evolving an equal-mass black-hole binary through its final orbit, merger, and ringdown using innovative techniques different from those previously employed (Pretorius, 2005a). Shortly after, independent breakthroughs were made by the groups at the University of Texas at Brownsville (UTB) and NASA’s Goddard Space Flight Center, who discovered the “moving punctures” method, leading to successful black-hole mergers (Baker et al., 2006b; Campanelli et al., 2006a).

The moving-puncture approach was rapidly adopted by most researchers due to its underlying techniques and led to a flurry of results in 2006. Several firsts were achieved, including simulations of unequal-mass black holes and the resulting recoil of the remnant hole, mergers of spinning black holes, long waveforms (~7 orbits), comparisons with post-Newtonian results, and systematic parameter studies in numerical relativity.

The following years saw further advancements, including the discovery of “superkicks” with recoil velocities exceeding 1000 km/s in 2007, and successful black-hole binary mergers with a mass ratio of q = 10 in 2008. Long-term evolutions of generic spinning and precessing black-hole binaries were accomplished in 2009. The state of the art continued to progress with the first simulations using spectral numerical techniques and initial steps towards modeling gas flows around merging black holes.

These breakthroughs in numerical relativity have found various applications, such as comparisons with post-Newtonian expressions for waveforms, astrophysical computations of black-hole merger rates, and the development of templates for gravitational-wave data analysis. The study of black-hole mergers using numerical relativity has seen significant advancements and remains a thriving area of research.


In this section, we delve into the mathematical and numerical principles that form the basis of present-day black-hole merger simulations. We emphasize the critical challenges in achieving accurate evolutions. For more in-depth information, interested readers can refer to Alcubierre (2008) or Gourgoulhon (2007).

A. Einstein’s Equations

In numerical relativity, the central task revolves around solving Einstein’s field equations, represented as Gun = 87 Tuv. Here, Einstein’s tensor Ga describes spacetime curvature, while the energy-momentum tensor Tu contains matter sources, with indices u = 0 for “time” components and 4 = 1, 2, 3 for “space” components. The dependence of Guv on the first and second derivatives of the metric tensor 9uv is crucial.

The metric tensor characterizes spacetime geometry and is defined by the infinitesimal spacetime interval ds through the equation ds? = gunda”dx”. The implications of the metric are evident, as it determines the paths light rays must follow in spacetime when ds? = 0.

In the context of vacuum black-hole spacetimes, Einstein’s equations simplify to •guv – tuv = 0, with D being the flat-space wave-operator and tuv encompassing all nonlinear terms in the metric. This form resembles a simple wave equation, enabling its Cauchy problem to be solved by specifying the metric and its first derivative on an initial spatial surface and then integrating in time.

Numerical simulations of black-hole spacetimes involve various formulations and coordinate conditions tailored to the specific problem. Nonetheless, common practices entail starting with an initial three-dimensional slice of spacetime and evolving it forward in time. These developments have enabled the generation of gravitational waveforms from black-hole binary sources over multiple cycles, marking a significant achievement in this field.

B. The Cauchy Problem

The Cauchy problem in numerical relativity involves finding solutions to the field equations based on initial data specified on an initial spatial hypersurface. To simplify the investigation, a “3+1" formulation is commonly used, which divides the spacetime into three-dimensional spatial slices parametrized by a time coordinate. This 3+1 formulation, inspired by ADM (Arnowitt et al., 1962), expresses the line element using the lapse function (a), shift vector (B’), and spatial three-metric (V = 9ij).

The lapse and shift represent coordinate freedom in the metric, allowing us to choose these quantities arbitrarily. However, the three-metric V; and its spatial derivatives determine the intrinsic curvature of the slice, carrying crucial information about the gravitational field, and thus they are constrained by the physics.

To understand the meaning of the lapse and shift, consider two successive spatial slices separated by a small time interval dt. An observer along a normal vector to the first slice will measure an elapsed proper time of dr = adt in evolving to the second slice, accompanied by a change in spatial coordinate of da’ = -B’dt.

To deal with Einstein’s second-order time equations, we need to specify the initial time derivative of the three-metric. Instead of directly specifying this derivative, a new quantity called the extrinsic curvature (Kij) is introduced. It measures the change of the normal vector along the slice and provides an extrinsic measure of the curvature of the three-dimensional spatial slice with respect to its embedding in four-dimensional spacetime.

Of the ten component equations from Eq. (6), six determine the time-evolution of the metric, while the remaining four are constraint equations that must be satisfied on a spatial slice at any given time. These constraint equations, under appropriate choices of time coordinate and assuming vacuum spacetime, are equivalent to the condition Gov = 0. The Hamiltonian constraint (Goo = 0) and the momentum constraint (Goi = 0) involve conditions on the extrinsic curvature Kij, which play a critical role in calculating initial data.

In summary, solving these constraint equations is essential for obtaining an initial spatial slice consistent with Einstein’s equations in numerical simulations.

FIG. 2 The 3+1 split into space and time. Two spatial slices at t = 0 and t = dt are depicted. @ is the lapse, and adt represents the proper time lapse between slices. B° is the shift, and B° dt represents the amount by which the spatial coordinates shift between slices. n° is normal to the slice at t = 0. If a ray parallel to n° intersects the t = 0 slice at a point a”, then it will intersect the t = dt slice at a” – B° dt. t° = an° + Ba is a coordinate time vector. If a ray parallel to t° intersects the t = 0 slice at point x”, then it will also intersect the t = dt slice at y°

C. Representing Black Holes in Numerical Spacetimes

In numerical simulations, representing a complex object like a black hole poses challenges due to its physical and coordinate singularities. However, two successful strategies have emerged to address this issue.

One approach takes advantage of the unique topology of black holes. As originally shown by Einstein and Rosen (1935), a black hole can be seen as a “wormhole” connecting two separate worldsheets. By constructing a spatial slice of Schwarzschild spacetime while excluding the interior of the event horizon, we can identify the spherical boundaries of two identical copies of this space. These identified boundaries form the “throat” of the wormhole. Using Brill and Lindquist coordinates, we can continuously cover both worldsheets, allowing us to represent the black hole without encountering the physical singularity. Although this method introduces a coordinate singularity called a “puncture” at r’ = 0, it can be confined to a single scalar variable, making it manageable in numerical simulations.

Another strategy, proposed by Unruh and referred to as excision, involves excluding the interior of the black hole from the computational domain. Since no physical information can escape the event horizon, and all physical information propagates inwards toward the physical singularity, the excised region does not affect the exterior. Extrapolation can be used for boundary conditions, and any non-physical numerical error that escapes the horizon is typically negligible. This approach allows more flexibility in the choice of coordinates but requires precise positioning of the excision boundary and accurate extrapolation.

Both the “puncture” and excision methods prove useful in representing black holes on the initial data slice. Remarkably, these representations can also be made robust enough to persist as the black holes evolve during the numerical simulation.

FIG. 3 In the wormhole representation of a black hole, the initial slice typically just touches the horizon. The upper “sheet” represents the exterior space, while the lower sheet is a duplicate, joined to the upper sheet by a “throat”

D. Initial Data

To successfully simulate the merger of black holes, the first step is to find initial data for astrophysically realistic inspiraling black holes. In Newtonian gravity, simulating orbits of stars is straightforward, where we specify masses, spins, positions, and velocities of stars as point particles and numerically evolve the system until it settles into orbits appropriate for finite-sized bodies. However, in general relativity, the initial data must satisfy the constraint equations (11) – (12), which involve the 3-metric (vi) and the extrinsic curvature (Ki). Since there is no apparent connection between these field variables and the astrophysical properties of inspiraling black holes, obtaining suitable initial data presents a significant challenge.

Building on Lichnerowicz’s earlier work (1944), York developed a general procedure in the 1970s to solve the constraint equations and generate initial data for the Cauchy problem (see York Jr. (1979) for a review). This involves solving a coupled elliptic system of four nonlinear field equations. Although some simplifying assumptions are made, resulting in a loss of generality and astrophysical realism for two black holes, it has been observed that the final orbits and waveform signatures from the black-hole evolution are largely insensitive to this level of detail, especially for equal-mass non-spinning black holes.

The first simplification is to choose traceless extrinsic curvature (K = 0), which separates the Hamiltonian constraint from the momentum constraints, enabling separate solutions. To handle multiple black holes, the initial slice is generally assumed to be conformally flat, where the three metric is the product of a scalar conformal factor with a flat metric (Vij = W* Sig). This simplifies the problem into solving a typically nonlinear equation for the scalar field v.

Brill and Lindquist (1963) found a simple solution for N black holes momentarily at rest, but these black holes lack momentum and spin, making it less useful for general applications. Bowen and York Jr. (1980) provided an explicit solution for the momentum constraint, paving the way for calculating initial data for multiple black holes with specified linear and angular momenta. The widely used Brandt-Brügmann puncture data allows arbitrary momenta and spin due to its ease of implementation. Another approach, the Conformal-Thin-Sandwich (CTS) approach developed by York Jr. (1999), offers additional advantages, such as the ability to achieve higher spin parameters.

Both the CTS and Brandt-Brügmann methods provide ansatzes for constructing three-dimensional binary black-hole initial field data with specified particle-like parameters, masses, positions, momenta, and spins. However, some spurious radiation may be generated during initial transient dynamics as the system relaxes, due to unavoidable consequences of conformal flatness. Although these transient effects are generally negligible for simulations lasting several orbits, they can affect the quality of the simulation numerically.

Most black-hole binary simulations aim to represent astrophysical systems that have circularized before significant gravitational radiation. Therefore, initial data configurations for quasicircular orbits are used. The CTS approach naturally allows for the imposition of initially circular motion. For more realistic inspiraling trajectories, post-Newtonian (PN) approximations or iterative procedures involving short numerical evolutions can be used to achieve higher accuracy in the initial data.

E. Numerically Friendly Formulations of the Evolution Equations

The Einstein evolution equations, responsible for the time-development of the initial data, form a set of at least six coupled, nonlinear propagation equations. The specific formulation of these equations depends on the chosen evolved variables and the incorporation of constraint identities. While numerous formulations are possible, not all of them are equivalent from a numerical standpoint.

The choice of formulation can significantly impact how numerical errors evolve over time and whether these errors can converge to zero as the computational grid resolution increases.

  1. Hyperbolicity and Well-Posedness: An ideal formulation is considered well-posed if it satisfies two conditions: (i) for a given set of initial data, a unique solution exists, and (ii) solutions depend continuously on perturbations of the initial data. To assess whether a formulation is well-posed, it can be represented as a matrix equation (Eq. 17). Strongly hyperbolic formulations possess real eigenvalues and complete eigenvectors, ensuring well-posedness and bounded solutions. Conversely, weakly hyperbolic formulations have real eigenvalues but incomplete eigenvectors, leading to unbounded solutions that may grow at rates depending on the initial data and resolution. Strong hyperbolicity is crucial for stable black-hole merger simulations.
  2. Harmonic Formulations: One line of formulation development involves “harmonic coordinates,” satisfying the wave equation, leading to manifestly hyperbolic Einstein’s equations. Generalized harmonic coordinates were later introduced to retain hyperbolicity while relaxing restrictions. Additionally, constraint-damping terms were incorporated to enhance stability. The formulation can be second-order in both time and space or first-order-in-time for more efficient numerical integration.
  3. ADM-based Formulations: Another line of development started with the ADM formulation of Einstein’s equations. Refinements led to strongly hyperbolic versions by promoting certain derived quantities to independently evolved variables. Among these formulations, the BSSN (Baumgarte-Shapiro-Shibata-Nakamura) formulation emerged as successful and robustly stable for numerical purposes. Analytic justification for its success was found, and subsequent refinements, such as specific terms to damp constraints, further improved stability. For evolution of punctures, adjustments were made in the form of the conformal factor to avoid singularities.

In other words, finding numerically friendly formulations is essential for accurate and stable simulations of black-hole mergers. Strong hyperbolicity, stability, and avoidance of singularities are key factors in determining the success of a particular formulation.

F. Gauge Conditions

The choice of gauge or coordinate conditions, similar to the formulation choice, significantly impacts the numerical stability of simulations, especially in the presence of extreme black hole conditions like physical singularities, coordinate singularities, strong-field gradients, and dynamic surrounding space-time. The selected coordinates must effectively handle these features in a numerically tractable manner.

  1. Choosing the Slicing and Shift: An early concern was avoiding the physical singularity of a black hole through a suitable slicing condition. Maximal slicing, which required solving a numerically expensive elliptic equation, was one such choice. Later, the Bona-Massó family, a generalization of harmonic slicing, was introduced as a hyperbolic evolution system, with the 1+log slicing being more computationally efficient while still avoiding singularities. To counter large field gradients near black holes, a hyperbolic shift condition called the “Gamma driver” was employed, which effectively stabilized the coordinates and accommodated black hole mergers.
  1. Moving Punctures: Initially designed for punctures at fixed coordinates (comoving coordinates), these conditions encountered problems with large gradients and pathologies during black hole mergers. Modifications were made to allow punctures to move with robust stability. These changes, along with refinements to eliminate zero-speed modes, led to the successful “moving puncture” method widely used in numerical relativity simulations. The gauge conditions, coupled with the BSSN equations, also ensured strong hyperbolicity.
  1. Generalized Harmonic Coordinates: Generalized harmonic coordinates allowed source terms in the wave equation for each coordinate, providing more flexibility in gauge driving. A successful example involved keeping the lapse close to unity while the spatial coordinates remained harmonic. Another gauge driver was developed to dampen extraneous gauge dynamics during binary black hole inspiral and merger.
  1. Other Coordinate Techniques: Incorporating multiple coordinate patches allowed for local physical geometry considerations to optimize accuracy and construct physical boundary conditions between patches. “Dual coordinates” were introduced to combine the advantages of comoving and non-comoving coordinates, mitigating steep field gradients and avoiding interpolation errors while allowing stationary excision boundaries.

Overall, selecting appropriate gauge conditions is crucial for achieving numerical stability and accurate simulations of black hole mergers. Different gauge choices offer various advantages and trade-offs, and the moving puncture method, in particular, has proven to be very successful in the numerical-relativity community.

G. Numerical Approximation Methods

After making analytic choices for initial data, formulation, and gauge conditions, the focus shifts to finite, discrete computations that numerically approximate these analytical specifications. The simulated domain must have a finite extent, and various conformal compactification schemes have been tested to map spatial or null infinity to a finite boundary. However, it is more common to impose artificial boundaries at finite spatial coordinates and use either radiative boundary conditions or constraint-preserving boundary conditions.

To handle inward-propagating errors from the outer boundaries, mesh refinement is often employed to push the outer boundaries away from regions of interest, such as wave sources and extraction regions. Grid points’ density is strategically selected to be highest near strong field gradients of black holes, moderate in wave propagation regions, and coarser beyond that. For highly dynamic simulated black holes, an algorithm for automatic mesh refinement is necessary, and interpolation is needed at interfaces between refinement regions.

Spatial derivatives are computed using either finite differencing stencils across uniformly spaced grid points, achieving up to eighth-order accuracy, or the spectral approach, which involves computing coefficients of an expansion in basis functions on collocation points to obtain analytical derivatives. Dissipative terms are often added to the evolution equations in finite differencing methods to reduce noise, while spectral methods may include a smoothing step.

Time advancement occurs in discrete steps, and the widely used fourth-order Runge-Kutta algorithm is preferred due to its superior accuracy and efficiency. In some codes, the timestep size varies with spatial resolution for greater efficiency, albeit with the complexity of time interpolation. Implicit timestepping schemes are also being explored to significantly increase the timestep size.

H. Extracting the Physics

One of the most significant outcomes of simulating merging black holes is computing the emitted gravitational radiation. For this purpose, the “Weyl tensor,” Cabed, is calculated. Unlike the Einstein tensor, the Weyl tensor has degrees of freedom that can exist without a massive source, making it a useful measure for radiation. In a special “transverse-traceless” spherical coordinate system, the Weyl scalar V4 can be expressed as a complex quantity. This scalar is associated with outgoing radiation at spatial infinity and approximates the radiative behavior at large but finite radii in numerical simulations.

To extract specific modes of the field, V4 is integrated over a sphere multiplied by spherical harmonic functions, particularly spin-weight -2 spherical harmonics. The dominant contribution is usually the l = 2 quadrupole moment.

Another essential aspect of analyzing spacetimes is obtaining information about the black holes. This is commonly done through apparent horizons, surfaces where all light rays must go inward. During the simulation, surfaces with minimized light ray expansion are found to locate the apparent horizons. From these horizons, properties like mass, spin magnitude, and linear momentum can be defined, which characterize the physical properties of the black holes and are used to label spacetimes when comparing with PN theory.


The final merger of a binary black hole occurs in the strong-field regime of general relativity. Numerical relativists have long been curious about what they would find when exploring this area through black-hole merger simulations. Questions arise about the appearance of merger waveforms, the impact of artifacts in the initial data on the mergers, and how unequal masses and spins might alter the outcomes compared to the simplest case of nonspinning, equal-mass mergers. Moreover, researchers are interested in uncovering any unexpected phenomena in the process. Recent breakthroughs in black-hole merger simulations have provided valuable tools to address these inquiries.

In this section, we delve into the dynamics of black-hole mergers and the resulting gravitational waveforms. We start with the Lazarus approach, which combines analytic methods with short numerical simulations to produce hybrid waveforms. Then, we embark on a comprehensive discussion of black-hole merger physics, beginning with nonspinning, equal-mass mergers. We subsequently examine unequal-mass mergers and mergers involving spins. Throughout our exploration, we provide a historical context while emphasizing the key physical findings.

A. First Glimpses of the Merger: The Lazarus Approach

In the late 1990s, numerical relativity had progressed to a point where it was possible to perform brief simulations of two black holes in three dimensions. These simulations lasted approximately 10M, a small fraction of the estimated binary orbital period near the Innermost Stable Circular Orbit (ISCO). While most researchers focused on extending the duration of these simulations, a small group of collaborators adopted a different strategy (Baker et al., 2000).

Their innovative idea was to use full numerical-relativity simulations to calculate the strong-field approach to the merger during the initial 10M of simulation time. Then, they employed perturbative techniques to calculate the remaining evolution (Baker et al., 2002a). They started with traditional puncture black holes on nearly circular orbits near the ISCO and evolved them towards merger. Just before the simulation became unstable, they stopped the calculation and extracted the physical data for the merging black holes and the gravitational waves produced. They interpreted this data as initial data for a highly distorted single black hole and evolved it using black-hole perturbation theory. This method of reviving a nearly completed calculation was named the Lazarus approach.

The hybrid waveforms produced by the Lazarus collaboration provided the first insights into what might be expected from the final merger of black-hole binaries (Baker et al., 2001, 2002b). They conducted a series of simulations, varying the initial black-hole separations and the time at which they transitioned from fully numerical evolution to perturbed black-hole evolution. The dominant quadrupole component of the gravitational waveform, corresponding to the + polarization, for equal-mass, nonspinning black holes is shown in Figure 1. The waveform exhibited a simple shape, smoothly connecting the end of an inspiral chirp with a damped ringdown sinusoid.

By adjusting the mass and spin of the background hole for the perturbative evolution, Baker et al. (2002b) determined that the merger remnant was a spinning Kerr black hole with mass MKerr ~ 0.97 Minitial and spin aKen ~ 0.7 MKerr. The Lazarus method was also applied to mergers of black holes with spins either aligned or anti-aligned with the orbital angular momentum, and in all cases, a similar waveform shape was observed. The researchers wondered if this simple shape was generic and whether it would arise in situations where the black holes completed one or more orbits before merging.

(Figure 4 shows Lazarus waveforms for equal-mass, nonspinning black-hole mergers, from Baker et al. (2002b). The real, l = 2, m = 2 part of the Weyl scalar, corresponding to the + polarization measured on the system’s axis, is depicted for ten simulations with different initial black-hole separations, each with two transition times to perturbative evolution.)

B. Mergers of Equal-Mass, Nonspinning Black Holes

  1. The First Merger Waveforms: In the early 2000s, Frans Pretorius achieved a groundbreaking simulation of two equal-mass black holes undergoing a plunge orbit, merger, and ringdown, providing the first robust gravitational waveforms (Pretorius, 2005a). The resulting waveforms showed a simple shape, and the merger yielded a spinning black hole with a spin parameter af ~ 0.70M, where M is the mass of the final black hole. Pretorius used a generalized harmonic formulation and excised the black holes, evolving their motion across the grid using adaptive mesh refinement.
FIG. 5 The first gravitational waveform for black holes evolving through a single plunge orbit, merger, and ringdown, achieved by Pretorius (2005a). The waves were extracted at four radii, and shifted in time to overlap for comparison. Note the amplitudes decrease for larger r, due to lower resolution in the outer regions

Around the same time, two independent research groups at UTB and Goddard developed the “moving puncture” approach, successfully evolving mergers of puncture black holes within the traditional BSSN approach (Baker et al., 2006b; Campanelli et al., 2006a). Their simulations produced similar waveforms, confirming the robustness of the approach.

FIG. 6 The trajectories of the black-hole punctures from a run by Campanelli et al. (2006a). The individual apparent horizons are shown at times t = 0,10M and 18.8M. The solid circles denote the centroids of the apparent horizons every 2.5M during the evolution. The first common horizon, marking the time of merger, forms at 18.8M, just before the punctures complete 1/2 orbit.

2. Universal Waveform: Astrophysical black-hole binaries generally undergo many orbits before merging, leading to quasicircular orbits at the final inspiral. Surprisingly, simulations of equal-mass, nonspinning black-hole mergers starting from different initial separations produced a universal waveform shape, subject only to rescaling with the total system mass M. Longer waveforms starting many cycles before merger were also generated, allowing for more advanced applications in gravitational-wave data analysis and comparisons with PN analytical treatments.

FIG. 7 Gravitational waveforms from the black-hole merger simulations by Baker et al. (2006b). The real part of the l= 2, m = 2 mode of rV was extracted from the numerical simulation on spheres of radii r = 20M, and 40M for the medium- and high-resolution runs. The waveforms extracted at different radii have been rescaled by r, and shifted in time to account for the wave propagation time between the two extraction spheres. Curves are shown for medium (M /24 in the finest grid) and high (M/32 resolutions. A waveform from a Lazarus calculation starting with the same initial data is shown for comparison (Baker et al., 2002b). The Lazarus waveform shown in this figure is scaled differently from those in Fig. 14.

3. Longer Waveforms: Research groups extended the simulations to cover more orbits before merger, resulting in longer waveforms. The UTB group achieved waveforms starting from ~ 7 orbits before merger, while the Jena group simulated a binary inspiralling for nine orbits before merger. The Caltech-Cornell group achieved the longest evolution, starting 16 orbits before merger with very small initial eccentricity. Comparisons between different groups’ waveforms showed good consistency, indicating sufficient accuracy for detection with current and future gravitational-wave detectors. The Samurai project further explored waveform consistency among different numerical codes.

FIG. 8 Puncture tracks from equal-mass, nonspinning black-hole merger simulations starting at increasingly larger separations by Baker et al. (2006a). For clarity, the trajectory of only one black hole in each run is shown. The tracks lock on to a universal trajectory ~ 1 orbit before merger (denoted by *).
FIG. 9 The universal waveform (l = 2, m = 2 mode) produced by the four simulations whose puncture tracks are shown in Fig. 8 (Baker et al., 2006a). The signals have been shifted in time so that the peak radiation amplitude occurs at t = 0. The agreement among the waveforms is excellent, with differences ~ 1% for the final burst of radiation starting at t ~ -50Mp.
FIG. 10 Comparison of three waveforms from simulations of equal-mass black-hole mergers computed by Pretorius and the UTB and Goddard groups (Baker et al., 2007b). The UTB and Goddard simulations used nonspinning black holes, whereas Pretorius’ initial conditions produced co-rotating black holes, each with a small spin (a/M)1,2 = 0.08. The differences between Pretorius’ waveform and the others during the ringdown, t/My ~ 25, are consistent with his simulation producing a slightly faster rotating final black hole due to these initial spins. Reproduced by permission of the Institute of Physics.
FIG. 11 Puncture tracks for equal-mass black-hole mergers with initial orbital eccentricity e = 0.1 (left) and e = 0.3 (right) from Hinder et al. (2008). In these cases, the initial eccentricity is radiated away and the binaries circularize be fore merger.
FIG. 12 Real part of (2, 2) strain mode (labeled h+ here) for black-hole mergers shown in Fig. 11|from Hinder et al. (2008). These runs have initial orbital eccentricity e = 0.1 (top) and e = 0.3 (bottom) and enter a universal plunge by t ~ 50M g before the gravitational radiation reaches its peak value.

3. Extending Waveforms to Longer Durations

The previous black-hole simulations covered only the last few orbits before the merger, marking a significant achievement in solving the two-body problem for the final merger of equal-mass Schwarzschild black holes driven by gravitational radiation reaction. These simulations offered valuable insights into the dynamics and waveforms of black holes merging in the strong-field regime. They also paved the way for detecting the final gravitational-wave burst from black-hole mergers (Baumgarte et al., 2008). However, for more advanced applications in gravitational-wave data analysis and comparisons with PN analytical treatments, longer waveforms that start many cycles before the merger become crucial.

The Goddard group made a breakthrough by producing the first long waveforms for equal-mass, nonspinning black holes starting approximately seven orbits or 14 gravitational-wave cycles before the merger (Baker et al., 2007c,d). They used improved initial conditions with low eccentricity to enhance accuracy during the long inspiral phase. These long waveforms were successfully compared with PN results, with a particular focus on waveform phases (Baker et al., 2007d).

FIG. 13 Trajectories for the merger of equal-mass, nonspin-ning black holes computed by the Caltech-Cornell group using their spectral evolution code (Scheel et al., 2009). The red and black circles/ellipses are the initial/final coordinate locations of the apparent-horizon surfaces, and the blue thick contour is the common apparent horizon just after it appeared. The black holes complete 16 orbits before merging. Figure kindly provided by H. Pfeiffer.
FIG. 14 Gravitational waveforms from the Caltech-Cornell merger simulation seen in Fig. 13| showing the l = 2, m = 2 component of Re(rV4) (Scheel et al., 2009). The left panel shows a zoom of the inspiral waveform, and the right panel a zoom of the merger and ringdown.

Subsequently, the Jena group simulated a binary inspiralling for nine orbits (18 gravitational-wave cycles) before the merger, allowing for quantitative comparisons with PN phase and amplitude. They used higher-order finite differencing and initial binary parameters calculated using the PN approximation to reduce the initial eccentricity significantly, enabling precise measurements of the waveform phase (Hannam et al., 2008c).

FIG. 15 The gravitational-wave phase • as a function of the dimensionless frequency Mw for the five codes from the Samurai comparison (Hannam et al., 2009a). The scale along the top of the panel labels the frequency in kHz scaled with respect to the total binary mass in solar units.
FIG. 16 Amplitudes of the l = 2, m = 2 component of the gravitational waves produced by the five codes in the Samurai comparison (Hannam et al., 2009a) are shown as a function of frequency.

The Caltech-Cornell group holds the current record for the longest and most accurate black-hole binary evolution, starting 16 orbits and 32 gravitational-wave cycles before the merger. They achieved this using their spectral code and very small initial orbital eccentricity, resulting in highly accurate simulations throughout the long inspiral, merger, and ringdown (Scheel et al., 2009).

Comparing waveforms from different groups remains important to ensure higher accuracy and their applicability in gravitational-wave data analysis. The Samurai project has been at the forefront of studying the consistency of black-hole binary waveforms. They focus on equal-mass, nonspinning binaries, starting with at least six orbits (or 12 gravitational-wave cycles) before the merger and continuing through the ringdown. The amplitude and phase of the l = 2, m = 2 mode of the gravitational waveforms are compared among five independent numerical codes, indicating good consistency and sufficient accuracy for detection with current and planned ground-based detectors.

C. Black Hole Mergers with Different Masses

In 2006, numerical relativists advanced their simulations by studying nonspinning black-hole binaries with unequal masses, introducing the mass ratio g as a new parameter. This required adaptive mesh refinement to achieve sufficient resolution around the smaller black hole, which moves faster and demands smaller timesteps for its evolution, making the simulations more challenging. As of now, numerical relativists can simulate mass ratios up to q = 10 (González et al., 2009).

FIG. 17 Puncture tracks for an unequal-mass nonspinning binary with q = 2 produced by the Caltech-Cornell spectral code; figure kindly provided by H. Pfeiffer. The solid (red) curve gives the trajectory of the smaller black hole, and the dotted (blue) curve the path of the larger one. The red and blue contours are the final coordinate locations of the apparent-horizon surfaces, and the black contour is the common apparent horizon just after it appeared.

1. Mode Analysis and Gravitational Waveforms: Decomposing the gravitational waveforms into spin-weighted spherical harmonic modes offers a comprehensive study of black-hole mergers. For the equal-mass case (q = 1), the dominant mode is the l = 2, m = ‡2 quadrupole mode, with the odd-m modes suppressed by symmetry. As the mass ratio increases, the sub-dominant modes become more significant, affecting both the source’s evolution and the emitted radiation. Studies by Berti et al. (2007b) on unequal-mass mergers with mass ratios in the range 1 < q < 4 revealed that the total energy emitted during merger scales ~ n? while the final black hole’s spin scales ~ n, where n is the symmetric mass ratio. They also observed that higher modes carry larger fractions of the total energy as q increases, with the l = 3 mode generally carrying ~ 10% of the emitted energy for q > 2. The Goddard group (Baker et al., 2008b) complemented this study with nonspinning unequal-mass mergers for mass ratios 1 < q < 6. They found that the simple waveform shape observed in equal-mass mergers extends to cases with g > 1 when gravitational waves are scaled by n. The waveforms show modest amplitude and phase modulations when viewed off-axis while preserving the overall simplicity in shape. Interestingly, they found that each individual spherical harmonic component is circularly polarized during the inspiral, merger, and ringdown, and the relationship between phase and frequency of the l = m modes remains consistent during merger and ringdown. González et al. (2009) conducted the highest mass ratio merger simulation to date with q = 10, finding that the gravitational waveforms’ simple shape extends through l = m = 5, the highest mode they studied.

FIG. 18 The amplitudes of several (l, m) modes of Mr¢4 for mass ratio q = 2 from simulations by Berti et al. (2007b). The initial burst of radiation is an artifact of the initial data, and the wiggles at late times result from numerical noise.
FIG. 19 Strain waveforms scaled by n for different mass ratios, from Baker et al. (2008b). h+ is calculated for different mass ratios using the l = 2 and l = 3 modes. The observer is located at distance r along the axis of the system.
FIG. 20 The l = m = 2 and l = m = 3 component of Re(rM W4) from the 10:1 mass ratio simulation by González et al. (2009)

2. A New Phenomenon: Recoils or Kicks: Mergers of unequal-mass black-hole binaries introduce a new phenomenon called recoils or kicks. As discussed in Sec. TI.A, gravitational waves emitted during a merger carry away linear momentum. If the pattern of gravitational-wave emission is asymmetrical, the center of mass and the remnant black hole recoil in the opposite direction to conserve momentum. Unequal-mass mergers, where the smaller black hole moves faster, result in more forward beaming of its wave emission than the larger hole, leading to net momentum flux parallel to the smaller hole’s velocity and an opposing kick at the center of mass. Since the black holes spiral together, instantaneous kicks accumulate, causing the final merged black hole to have a non-zero net recoil in the orbital plane. Accurate calculations of the kicks require full numerical relativity simulations of the final stages of binary inspiral, merger, and ringdown. Studies by the Goddard group (Baker et al., 2006c) and others revealed that the kick velocity depends on the initial separation of the black holes and can reach significant values, impacting astrophysical scenarios of black-hole growth and retention.

FIG. 21 Schematic drawing of the physical basis of the recoil or kick produced by a merging black-hole binary with unequal masses. A net flux of momentum Peiected is emitted parallel to the smaller hole’s velocity. Momentum conservation requires that the system center of mass recoil in the opposite direc-tion, Precoil. Reproduced with kind permission from Springer Science + Business Media: Fig. 1 of Hughes et al. (2005).
FIG. 22 The magnitude of the kick velocity for a nonspinning, unequal mass merger with g = 1.5 from Baker et al. (2006c). Results are shown from three simulations, with increasingly wider initial black-hole coordinate separations dinit. The small-dotted (purple) curve shows a 2 PN calculation starting with initial conditions commensurate with the dinit = 7.0M case (blue dotted line). Reproduced by permission of the AAS.

D. Mergers of Spinning Black Holes

Astrophysical black holes are typically spinning, which adds six more dimensions to the parameter space when considering the spin vectors of both black holes in a binary, making it a total of seven dimensions when including the mass ratio q. Initial explorations of this large parameter space began in 2006, starting with simple cases of equal masses and spins, and gradually expanding to more generic binaries.

  1. Gravitational Waveforms: The UTB group (Campanelli et al., 2006d) performed the first simulations of black-hole binary mergers with spin. They studied three cases: aligned spins, anti-aligned spins, and a nonspinning case for comparison. The aligned and anti-aligned cases showed longer merger times and more orbits due to a spin-orbit interaction, causing an effective repulsive or attractive force between the black holes. However, all three cases produced similar waveforms with a clean, simple shape. Moving beyond simple binary dynamics, the RIT group (Campanelli et al., 2009) pioneered the study of precessing black-hole binary mergers with unequal masses and non-aligned precessing spins. They found that the precession of the total system spin leads to the precession of the orbital plane, resulting in amplitude modulations in the waveforms’ sub-dominant modes.

2. Spinning Binary Mergers and Spin Flips: When nonspinning black holes merge, the final black hole’s spin direction aligns with the orbital angular momentum. However, when the individual black holes have spins not aligned with the orbital angular momentum, the final black hole’s spin direction may differ significantly from the initial spins, known as a “spin flip.” The RIT group (Campanelli et al., 2007a, b) explored spin flips in generic binaries, demonstrating how the final spin direction can deviate from the initial spins.

FIG. 24 The difference in the black-hole trajectories 21 – ⅝2 for the generic binary evolution from Campanelli et al. (2009). The precession of the system spin induces precession of the orbital plane.
FIG. 25 Precession-induced amplitude modulations in the waveforms from the generic black-hole binary merger from Campanelli et al. (2009). The l = 2, m = 1 mode of the strain (h+) is shown for the numerical-relativity simulation (solid) and a 3.5 PN calculation (dotted). The amplitude oscillations are induced by precession.
FIG. 26 Black hole trajectories for the spin flip case studied by Campanelli et al. (2007a). The black holes have equal masses, and equal spins initially pointing 45° above the initial orbital plane. The spin directions are shown as arrows along each trajectory every 4M until merger.

3. Kicks from Mergers of Spinning Black Holes: The mergers of asymmetrical spinning binaries can produce recoiling remnant black holes, known as kicks. Early simulations by Herrmann et al. (2007a), Koppitz et al. (2007), and Pollney et al. (2007) showed that spinning black-hole mergers could produce significantly larger kick velocities than nonspinning mergers. The phenomenon of “superkicks” was introduced by the Jena group (González et al., 2007a), where specific configurations resulted in extremely large kick velocities. Further studies by the RIT group (Campanelli et al., 2007c) and Schnittman et al. (2008) demonstrated the role of spin in generating superkicks and the importance of certain multipole moments in predicting recoil velocities accurately.

FIG. 27 The spin flip produced by the evolution in Fig. 26, from Campanelli et al. (2007a). The spin direction of the component holes, shown by the red arrows plotted every 4M until merger, varies continuously due to precession. The blue arrow shows the final black-hole spin, with direction flipped discontinuously from that of the component holes just prior to merger.
FIG. 28 Puncture trajectories from a superkick simulation carried out by González et al. (2007a). The black holes have equal masses and spins (a/M)1,2 ~ 0.8, initially oppositely directed in the orbital plane. During the evolution, the black holes move out of the original plane, and the final black hole recoils with velocity kick ~ 2500 km sI in the negative z-direction.


Einstein’s general relativity lays out predictions for the gravitational dynamics and generation of gravitational waves in the context of black-hole binary systems, encompassing all phases of their merging. As discussed, numerical relativity has become a valuable tool for obtaining concrete predictions, commencing from the final orbits and extending through the coalescence process until the remnant black hole reaches the ringdown phase. While there isn’t a strict cutoff for when numerical relativity can be applied in simulations, the substantial computational demands tend to restrict the simulations to a few tens of orbital periods. Post-Newtonian (PN) expansion, though an approximation of general relativity, is adept at producing precise and efficient forecasts over extended time spans, provided the black holes remain well-separated. However, as the black-hole velocities escalate during the later stages of the merger, PN theory’s accuracy diminishes.

A comprehensive comprehension of general relativity’s forecasts for the merger procedure necessitates the integration of both methodologies. The maturation of numerical relativity has driven endeavors to synthesize a comprehensive comprehension of black-hole mergers by incorporating findings and techniques from both fields. A comprehensive overview of PN formalism and outcomes exceeds the scope of this discussion, as decades of exploration have already covered this territory (Blanchet, 2006, 2010; Blanchet et al., 2008; Damour and Nagar, 2010; Schäfer, 2010). Our focus here lies solely on pivotal areas where symbiotic interactions between numerical relativity and PN findings are enriching the general-relativistic comprehension of black-hole binary systems.

The extensive array of research at the nexus of these two approaches can be categorized into four groups, each representing distinct manners in which PN theory influences numerical-relativity research: 1) furnishing independent results for verification and comparison with numerical relativity; 2) supplying models for the initial parameters essential to initiate astrophysically realistic numerical-relativity simulations; 3) offering insights for interpreting numerical-relativity outcomes; and 4) potentially forming the foundation for empirical models that encapsulate the amalgamated insights drawn from both PN and numerical-relativity investigations.

A. Autonomous Post-Newtonian Dynamics and Waveforms

As previously indicated, Post-Newtonian (PN) theory builds upon an approximate expansion of Einstein’s equations using velocity powers, where E” ~ (v/c)^n. This provides corrections from general relativity to the Newtonian framework, particularly in the context of small velocities. This approach has a lengthy history of success, being a primary framework for deriving specific predictions from general relativity in various physical and astrophysical observations. Given the generally modest velocities observed, low-order PN theory has been effective in testing general relativity across contexts such as Earth-orbit experiments, solar-system dynamics, and even precise pulsar timing in binary neutron star systems (Will, 2006).

In the domain of gravitational waves, researchers have employed PN theory to represent signal expectations based on general relativity for gravitational-wave searches, both for anticipated observations of black-hole and neutron-star binary inspirals, and in the design of current and future gravitational-wave instruments. For near-circular inspirals, where the binary separation R scales with the total system mass M as R ~ M, the PN expansion offers reliable predictions, particularly at larger separations.

FIG. 29 Comparison of NR and PN waveforms from Baker et al. (2007d) provided mutual affirmation of the results from each approach and showed that in combination NR and PN results could treat the complete signal.

However, as the black holes lose energy and draw closer, their velocities increase, necessitating higher-order PN terms for accurate predictions. Presently, PN predictions are available up to the 3.5 PN order (2.5 PN for spin terms). While these expansions provide satisfactory predictions for data analysis from current instruments for most of the merger process, their accuracy diminishes as the separation nears a critical point (R/M ~ 10). Consequently, PN theory becomes less reliable during the final orbits and the merger phase.

The most substantial gravitational radiation occurs in the late stages of inspiral or merger, where the internal consistency of the PN approach remains challenging to assess. Numerical simulations can handle late stages, but computational constraints limit the simulation duration and the initial black-hole separation. Thus, there’s a need to determine if numerical simulations can overlap with the aspects handled successfully by PN methods or if an intermediate approach is necessary.

FIG. 30 Phase differences between a numerical simulation and various PN models from Boyle et al. (2007). The left panel shows phase differences from an initial (waveform) angular frequency of Mw = 0.04 up to Mw = 0.063, while the comparison in right panel extends this range to Mw = 0.10.

Early comparisons between numerical and PN results yielded promising signs that the gap between Numerical Relativity (NR) and PN could be bridged, although PN results were not uniquely determined for the final phase of inspiral. These comparisons verified the numerical results’ internal consistency and explored the convergence towards Einstein’s equations as resolution increased. External verification would strengthen the validity of these results.

The first comprehensive comparison emerged around 2006, examining equal-mass non-spinning black-hole mergers. The 3.5 PN waveform aligned with numerical outcomes over the last seven orbits before merger. This demonstrated the accuracy of PN waveforms approaching the final orbits, suggesting that PN and NR combined could analyze the complete waveform signal.

Different PN waveform variants (approximants) don’t always agree as well with numerical results. PN approximation involves an expansion in the merging black holes’ speed, v/c. PN waveforms express orbital phase and polarization amplitudes separately, given as expansions in the time to merger v/c or orbital frequency v/c. The TaylorT4 approximant showed strong agreement with numerical results, indicating that it could be accurate for late orbits approaching merger.

PN spin effects are derived only up to 2.5 PN order, which might limit waveform accuracy in certain scenarios. Some studies explored other configurations, such as precessing systems or eccentric non-spinning mergers. In comparing PN Fourier-domain approximants with numerical outcomes for different mass ratios, truncating Fourier-domain waveforms produced Gibbs oscillations resembling the physical quasinormal ring-down that ends numerical waveforms. This led to the development of phenomenological full-waveform models.

In other words, the interplay between PN theory and numerical relativity has unveiled the potential to combine their strengths for a more comprehensive understanding of black-hole binary systems, leading to accurate predictions across various stages of the merger process.

B. Comprehensive Analytical Full-Waveform Models

During the initial stages of waveforms, multiple Post-Newtonian (PN) approximations find consensus up to the 3.5 PN order – the most advanced level currently available – offering remarkably precise estimations. However, as time advances, discrepancies arise among various approximations, necessitating reliance on numerical-simulation outcomes for accurate predictions derived from general relativity. This prompts the question: Can an analytical description of waveforms be developed to effectively capture insights from both of these complementary theoretical approaches? The relatively straightforward features observed in numerically simulated merger waveforms offer promise for creating reasonably accurate approximate analytical waveform descriptions. These descriptions might have simple dependencies on binary system parameters, making them efficient tools for interpolating from sparsely sampled parameter spaces that have been explored through numerical simulations. Such waveform models would offer computational efficiency, compared to time-consuming numerical simulations, and could prove valuable across a diverse range of observational data analysis applications.

Notably, Ajith and collaborators have devised a Fourier-domain full-waveform model using a phenomenological approach akin to the method used by Pan et al. (2008). They place greater emphasis on aligning with numerical waveforms (Ajith, 2008; Ajith et al., 2008, 2007). Their process begins by combining information from PN and numerical-relativity (NR) based waveforms, stitching together time-series data for the dominant (l = 2, m = 2) component waveforms, which involves joining a lengthy PN precursor to a numerical merger waveform. This allows them to fit Fourier transforms of waveforms to a parameterized model that resembles waveform families used in phenomenological treatments of purely PN waveforms (Arun et al., 2006; Buonanno et al., 2003). Their waveform model for non-spinning binary mergers has the structure:

\[ h(f) = \dfrac{C}{d} \times

. \begin{cases}

. \left( \dfrac{f}{f_m} \right)^{-7/6} & \text{if } f < f_m \\

. \left( \dfrac{f}{f_m} \right)^{-2/3} & \text{if } f \geq f_m

. \end{cases}


The waveform amplitude A_eff is modeled in two separate power-law segments before and after the merger frequency \(f_m\), and a Lorentzian decay beyond the transition frequency \(f_r\) that marks the ringdown phase of the final black hole. The Fourier phase is represented by a single power-law expansion:

\[ \Phi(f) = \Phi_0 + \sum_{k=0}^{\infty} \dfrac{c_k}{(k – 5/3)} \left( \dfrac{f}{f_m} \right)^{k-5/3} \]

The free amplitude and phase coefficients are determined by fitting each waveform to stitched NR-PN hybrid waveforms. It has been observed that the dependence on the mass ratio parameter can be modeled by simple fits to quadratic functions of the symmetric mass ratio.

These waveforms have found application in gravitational-wave data analysis studies (refer to Section IVII.A; Ajith and Bose, 2009). A generalized version of this model for non-precessing spinning black-hole systems has also been developed by Ajith et al. (2009).

Alternatively, a more direct alignment between PN and NR can be achieved in the full waveform model, bypassing the need to stitch together PN and NR waveforms. The Effective-One-Body (EOB) family of PN Hamiltonian models shows promise as a tunable means of integrating both PN and NR results. Adjustments to the EOB model allow it to match waveforms derived from numerical-relativity simulations. While consistent with PN theory, this approach can match numerical-relativity waveforms without violating PN consistency.

In summary, advancements in waveform modeling seek to integrate insights from both PN and numerical relativity, offering efficient and accurate approximations of black-hole merger waveforms, and paving the way for enhanced observational data analysis in the field of gravitational waves.

FIG. 31 Comparison of EOB-model and numerical-relativity waveforms from Buonanno et al. (2007b) for the h+ component of the gravitational-wave strain from the merger of binary with mass ratio 4:1. The waveform is for an observer located an at inclination angle 0 = /3 from the axis of rotation.
FIG. 32 Agreement of refined EOB-model l = m = 2 waveform (Damour and Nagar, 2009) (red) with numerical result from Ref. (Scheel et al., 2009) for equal-mass non-spinning merger (black). Slight differences near the merger time are difficult to perceive without color.

C. Incorporating Post-Newtonian Insights into Numerical Initial Data Models

Prior to the initiation of a numerical simulation, a numerical relativist must define the initial configuration of the two black holes and the geometric fields that represent them. This process involves not only solving the constraint equations of general relativity but also providing a suitable ansatz for the free data, as discussed in Section IV.D.

In some scenarios, elements of information from Post-Newtonian (PN) theory can be employed to create an initial data model that more accurately captures the desired physical setup. To grasp how PN information can be effectively applied, it’s helpful to revisit the discussion surrounding numerical initial data. Typical simulations commence with a modeling ansatz, like the Brandt-Brügmann model, that generates initial field data based on a set of particle-like parameters tied to the black holes’ masses, spins, initial positions, and momenta. Usually, these initial data models don’t aim to account for the gravitational radiation fields previously produced by the black holes’ motion. Furthermore, they don’t utilize PN information to streamline the reduction of field degrees of freedom to particle parameters. However, PN theory’s insights are frequently utilized in selecting specific particle parameters aligned with the intended simulation, especially for circular inspiral setups.

Efforts have also been directed towards leveraging PN information more extensively to enhance the field ansatz. The majority of discussed simulations aimed to depict black holes in circularized orbits. Before enduring numerical simulation outcomes were accessible, comparative assessments between numerical initial data sets and PN-derived information provided a benchmark for results (e.g., Baker et al., 2002b; Caudill et al., 2006; Cook and Pfeiffer, 2004). These studies, without actual evolutions, concentrated on theoretical aspects of the black-hole configuration space, such as the Innermost Stable Circular Orbit (ISCO) (refer to Section II.A.1).

Nonetheless, subtle imbalances in the initial data, like excess angular momentum for the chosen separation, can induce eccentricity in the simulation, impacting the simulated radiation waveforms (Boyle et al., 2008b; Pfeiffer et al., 2007). While residual eccentricity can be mitigated through iterative processes (Pfeiffer et al., 2007), this necessitates multiple simulations spanning several orbits, a resource-intensive approach. To minimize eccentricity without resorting to iteration, several groups employ trajectory information from PN theory when setting up parameters for numerical initial data. This strategy facilitates simulations with low eccentricities (e < 0.002) for equal-mass non-spinning mergers with initial separations \(d_{\text{init}} > 10M\) (Husa et al., 2008b). This approach proves beneficial for simulations involving unequal masses and nonvanishing spins, although residual eccentricity is typically more pronounced (Walther et al., 2009).

Ongoing work involves developing techniques to more comprehensively integrate PN insights into numerical initial data modeling. The prevalent assumption of a conformally flat spatial metric in numerical initial data contradicts PN results at the 2 PN order (Damour et al., 2000), likely introducing modeling errors that generate spurious initial transients in simulations.

While not as mature as widely used models, an alternative approach utilizes metric information from PN theory for a nontrivial initial conformal metric (Jaranowski and Schäfer, 1998; Ohta et al., 1974; Schäfer, 1985). These techniques also enable the incorporation of information about prior radiation generated by the system before the “initial” time of the simulation (Kelly et al., 2007; Tichy et al., 2003).

D. Leveraging Post-Newtonian Theory to Interpret Numerical Findings

We’ve highlighted numerous specific domains where research in numerical simulations intersects with Post-Newtonian (PN) theory. However, these instances alone don’t present a complete understanding of the interplay between these two theoretical approaches. In a broader sense, PN theory serves as a foundational framework for interpreting numerical outcomes.

Numerical relativists widely draw upon the PN-derived background knowledge of black-hole binaries. Many pivotal phenomena, such as the Innermost Stable Circular Orbit (ISCO), spin-orbit interaction, spin precession, and the concept of ‘gravitational rocket’ kicks, were initially explored through the PN approximation. This foundational understanding laid the groundwork for subsequent numerical investigations. Throughout this review, numerous smaller examples underline this interrelation. Notably, we emphasize the significance of PN outcomes in comprehending the mergers of spinning black holes (refer to Section V.D) and in understanding kicks (Section V.C). The development of analytic formulas to approximately describe the final state of black holes in terms of mass, spin, and momentum has also heavily relied on insights derived from the PN treatment (see Section VIII.A.11).

Utilization in Gravitational Wave Data Analysis

A range of detectors, either active or in planning, aims to capture gravitational-wave signals originating from astrophysical phenomena. Notably, the ground-based interferometric detectors like LIGO (U.S.A.), GEO (Germany), Virgo (Italy), TAMA (Japan), and AIGO (Australia), with frequency sensitivity spanning 10¹ – 10⁴ Hz, as well as upcoming instruments like the Einstein Telescope (Freise et al., 2009), and the Laser Interferometer Space Antenna (LISA) with sensitivity in the range 10⁻⁴ – 10⁻¹ Hz. The latter operates in a complementary frequency domain. Noise interference poses a significant challenge for these instruments. The design sensitivity curve, taking into account various noise sources, reflects this challenge as shown in Figure 33 for LIGO and Virgo detectors.

FIG. 33 Design sensitivity curves for the LIGO & Virgo detectors (dashed lines), as well as the approximate curves used for the NINJA project (solid lines). LIGO sensitivity is effectively zero below ~ 30Hz, while Virgo does a little better (note that achieved sensitivities may not perfectly mirror the design curves). Taken from Aylott et al. (2009). Reproduced by permission of the Institute of Physics.

Efficient signal extraction demands matched filtering, where the noisy input signal is optimally correlated with template waveforms (Flanagan and Hughes, 1998; Wainstein and Zubakov, 1962). The degree of overlap between a physical signal hi(t) and a template waveform h2(t) is gauged through a frequency-space inner product (l•) defined as (Cutler and Flanagan, 1994):

= 4Re

Sn (f)

12 hi (f) h2* (f) df

Sn (f)


For effective application, a set of template waveforms hm(t; θ) is required, representing analytic model waveforms. Even in the case of a seemingly straightforward system like a black-hole binary, the combined potential configurations for source and detector yield a 17-dimensional parameter space. Certain parameters are intrinsic to the source: total binary mass M, symmetric mass ratio n, spins χ₁ and χ₂, coalescence time tc, eccentricity e, and eccentric phase φ. The remaining parameters pertain to the source-detector relation: luminosity distance DL, inclination angle i, orbital phase φ, waveform polarization ψ, and the binary’s sky position (θ, ϕ). Circularized binaries can disregard eccentricity parameters, simplifying the parameter space. Additionally, all observables display straightforward dependence on the total mass and redshifted mass, making (1+2)M and te sometimes treated as extrinsic parameters. Numerical simulations are left with a seven-dimensional parameter space (n, χ₁, χ₂), presenting a considerable challenge to cover thoroughly.

With a parametric set of template waveforms hm(t; θ), we evaluate the signal-to-noise ratio (SNR) for an incoming detector data stream st(t) = he(t) + n(t), incorporating gravitational wave he(t) and detector noise n(t):

P(θ) =

(hm (θ) | s)

√( hm (θ) | hm (θ)) (41)

In cases where the signal he(t) matches the model waveform hm(t; θ) for a parameter set θ = θ₀, the optimal SNR is achieved: Popt = P(θ₀).

Real model waveforms may be imprecise due to limited understanding of the underlying physics. In such cases, the achieved SNR cannot be optimal. The fitting factor (FF), defined as the reduction in signal strength due to imperfect templates, helps quantify this effect. If FF = 1, the signal is perfectly modeled in the parameter space of the templates. Achieving an accurate parameter estimation demands that the template waveforms closely correspond to the underlying source parameters. This is termed the fitting factor (FF) requirement, which is vital for identifying the source of radiation.

FIG. 34 Relationship between exact waveform he, model waveform with the same physical parameters hm, actual “best fit model waveform hm, and template bank waveforms hb, hw. Figure taken from Lindblom et al. (2008)]

Traditionally, templates for the detection of comparable-mass black-hole binaries have heavily relied on PN theory. While these templates provide accurate phase information during the inspiral phase, their accuracy diminishes as the merger phase approaches. Numerical insights into the final moments of the merger provide a novel perspective on waveform accuracy in these crucial regions.

A. Direct Influence of Merger Waveforms in Data Analysis

With the availability of full numerical merger waveforms, research teams scrutinized older predictions regarding the impact of merger segments on observability. The numerical outcomes for the equal-mass nonspinning scenario were consistent with the predictions made by Flanagan and Hughes (1998) for initial LIGO, though the merger SNR was smaller (Baker et al., 2007c; Buonanno et al., 2007a). Advanced LIGO and LISA experienced a significantly more substantial SNR boost from the merger, depicted in Figure 35 from Baker et al. (2007c).

FIG. 35 The importance of including NR merger waveforms in data analysis can be seen here for Advanced LIGO (top) and LISA (bottom), taken from Baker et al. (2007c). In both panels, the dashed curve is the achieved SNR (for a given redshifted source mass) when only the early inspiral signal (up to ~ 1000M before merger) is used; the dotted curve uses the signal between this time and when final plunge commences (1000M to 50M before merger; the thin solid curve uses just the merger waveform (starting 50M before merger). The thick solid curve is the combined result of using the entire waveform. The SNR, and hence distance reach, of the detector is clearly greatly enhanced by the final merger portion.

The significance of the final plunge and merger phases to the overall SNR depends on the binary’s total mass and its position within the detector’s sensitivity range. As illustrated in Figure 36, adapted from Baker et al. (2007c), LISA’s “characteristic signal strain” hchar(f) = 2fh(f) (dominant quadrupole radiation) reveals that low-mass binaries possess higher-frequency waveforms across all dynamic stages. For a binary with total mass M < 101Mo, the late-merger and ringdown signals extend beyond LISA’s sensitivity range. Hence, for lower masses, data analysis purposes are well served by inspiral-only waveforms.

FIG. 36 The importance of binary mass for the LISA observation window, demonstrated by the characteristic amplitudes (Baker et al., 2007c) of three different sources, relative to the rms noise amplitude of the LISA detector. On each hchar curve, we mark times before the peak amplitude (cir-cle) – one hour (square), one day (diamond) and one month (triangle).

The LIGO and Virgo detectors are sensitive to frequencies above ~30Hz, and the merger of binary systems generates a chirp gravitational-wave signal with a frequency that saturates at the dominant quasinormal mode. The frequency depends on the post-merger hole’s mass and spin, for instance, foNM ~ 1.75 × 101(M/Mo)-1 Hz for a nonspinning binary of total mass M. This implies that LIGO can’t observe merging binaries with a mass greater than ~600M. Post-Newtonian templates generally stop at fisco, the frequency at the innermost stable circular orbit, making their applicability limited near the last stages of merger.

FIG. 37 The power distribution in 10% segments of the waveform before merger, for a range of LIGO-appropriate equal-mass binaries. Note that as the total mass M increases, the final merger and ringdown account for a greater fraction of the total power. Adapted from Pan et al. (2008).

Research by Buonanno et al. (2009a) refined these estimates, revealing that inspiral-only PN templates suffice for all LIGO configurations for systems with M ~ 12M-. Beyond this mass, numerical results are essential to bridge the gap in reliable information. Figure 37, adapted from Pan et al. (2008), illustrates how the last stages of merger become increasingly crucial in initial LIGO for larger mass ranges.

FIG. 38 Maximum signal-to-noise ratio (SNR) p below which different “long” numerical waveforms are indistinguishable for the Advanced LIGO detector. Even the earliest, lowest-accuracy published long waveforms are indistinguishable for SN levels below 14, as demonstrated in this figure from the Samurai paper (Hannam et al., 2009a).

The “Samurai” project (Hannam et al., 2009a) showcased the consistency of long equal-mass nonspinning waveforms, validating their numerical accuracy. Direct phase and amplitude comparisons along with mismatch tests demonstrated the close match of quadrupole waveforms for different numerical codes. These results are significant for current ground-based detectors, displaying only minor differences for SNRs below 14.

Incorporating fuller harmonics for spinning and precessing binary systems enhances parameter estimation. Although these improvements vary among parameters, the introduction of full numerical merger waveforms is expected to further enhance parameter estimation, especially for portions of the detector’s window where the binary merger is observable. Ajith and Bose (2009) demonstrated improved parameter-estimate accuracy for Advanced LIGO and Advanced Virgo using phenomenological merger-inspiral-ringdown templates.

Estimates regarding LISA’s parameter estimation improvements are diverse (Babak et al., 2008; McWilliams et al., 2010; Thorpe et al., 2009). Discrepancies arise from factors such as signal-to-noise ratio (SNR) and the utilization of higher-order multipoles. While differing outcomes exist, this issue remains unresolved.

B. Development of Analytic Templates for Inspiral-Merger-Ringdown Waveforms

The long-term goal of numerical-relativity simulations of binary mergers is to create a comprehensive collection of gravitational-wave “templates” that encompasses the multidimensional space of astrophysical parameters. However, numerical simulations remain far too slow to run in real-time to filter incoming detector data streams. Therefore, there’s a pressing need to formulate analytic waveform expressions that encode the insights from these numerical simulations.

For nonspinning binary mergers, it has been demonstrated (Buonanno et al., 2009a) that inspiral-only Post-Newtonian (PN)-based templates become inconsistent with each other (and full waveforms) when the total binary mass exceeds 12M for both initial and advanced LIGO configurations. Beyond this mass threshold, it’s essential to employ full templates that incorporate numerical-simulation-based insights into the merger and ringdown phases.

FIG. 39 Accuracy in parameter estimation for EOBNR templates in the NINJA project. The top panel shows the fractional error in estimated total mass, (Minjected – Mdetected)/Minjected, for all detected signals, while the lower shows the error in merger time. Adapted from Aylott et al. (2009). Reproduced by permission of the Institute of Physics.

Pan et al. (2008) developed a set of templates that integrate numerical-relativity data. These templates were extensions of the Stationary-Phase Approximation (SPA) templates utilized in LIGO/Virgo data analysis. They incorporated input from equal-mass numerical waveforms provided by Frans Pretorius and also equal- and unequal-mass numerical waveforms from the Goddard group. These templates were effective, yielding fitting factors around 0.96, yet they extended source parameters into unphysical regions.

Currently, only a few faithful full-waveform template banks have been created, focusing on nonspinning binary systems with various mass ratios. The “phenomenological” templates by Ajith et al. (2008, 2007) use three-segment curves in frequency space, adjusting matching parameters based on numerical data from codes like BAM and CCATIE. The EOBNR templates by Boyle et al. (2008b) and Buonanno et al. (2007b) extend effective-one-body waveforms with a single “pseudo-4PN” parameter tuned by numerical data. Both template banks are faithful and have been used in data-analysis injection tests.

Figure 139 from Aylott et al. (2009) presents the accuracy of total mass M and time of merger extraction from NINJA detected injections when using EOBNR templates. These templates also excel in injection and filtering for the high-mass region of LIGO test analysis and the ringdown band. The NINJA project’s matched-filter analyses revealed that these templates performed comparably to inspiral-only templates for detection and significantly better for parameter estimation.

As significant spin is anticipated in astrophysical black holes at various scales, it’s crucial to incorporate spin into templates. Vaishnav et al. (2007) showed that matched filtering with nonspinning merger waveforms is effective in detecting binaries with substantial spins. However, their configurations lack certain spin-related effects. To address this, phenomenological templates by Ajith et al. have been extended to include non-precessing spinning systems (Ajith, 2009). These templates are shown to be effective and faithful in LIGO searches over non-precessing binary hybrid waveforms, a marked improvement over earlier zero-spin templates. These new templates are ultimately characterized by three parameters: total mass M, symmetric mass ratio n (1), and a single effective spin.

C. Utilizing Numerical Waveforms in Data Analysis Scenarios

An evident application of numerical waveforms in data analysis is their direct utilization for testing analysis algorithms. Given that current templates predominantly rely on Post-Newtonian (PN) information, which either extend heuristically to encompass merger and ringdown or cease before the merger phase, the existing detection and parameter estimation algorithms also operate on incomplete data.

A recent collaboration known as the NINJA project has taken a crucial step in this direction, directly inserting explicit numerical waveforms into a simulated LIGO dataset to evaluate detection efficiency and systematic errors (Aylott et al., 2009). This endeavor enlisted contributions from numerous active numerical-relativity groups, scaling waveforms to represent different-mass sources and producing 126 distinct injections. Nine data-analysis groups then assessed the post-injection data stream using a range of techniques, including matched filtering, burst analysis, and Bayesian methods.

In matched filtering, groups employed various templates such as the “TaylorF2" inspiral-only templates generated from PN theory, ringdown-only templates, and full inspiral-merger-ringdown templates like EOBNR and phenomenological waveforms. TaylorF2 and full templates achieved high detection rates, detecting approximately 80 out of the 126 injected signals in triple coincidence between LIGO detectors, whereas ringdown templates performed less favorably, detecting only 45 signals. For parameter estimation of total mass M, TaylorF2 estimates often exhibited errors of 40% or more, notably more than EOBNR and phenomenological estimates. While mass errors remained substantial with EOBNR templates, ranging between approximately -30% and +50% for most detected signals, they were considerably less biased compared to TaylorF2 templates.

Burst analysis, which involves extracting partial information from signals such as dominant frequency and approximate duration, was also employed. NINJA used two burst techniques – the “Q pipeline” and the Hilbert-Huang transform (HHT). The Q pipeline fits signals to a sine-Gaussian model, while HHT decomposes input data into intrinsic mode functions to generate a high-resolution time-frequency map. Both methods yielded competitive results with matched-filter searches.

Bayesian analysis aims to reconstruct the posterior probability density function based on the observed signal. NINJA employed two Bayesian approaches: Markov Chain Monte Carlo (MCMC) and Nested Sampling. Detection and parameter estimation were explored using various parametrized waveform models. While encouraging, current results are constrained by factors like the limited length of waveforms and simplified noise profiles, requiring further work to develop longer numerical waveforms spanning the LIGO sensitivity window and blending PN inspiral waveforms with late-inspiral merger numerical waveforms. Such efforts will likely extend to other gravitational-wave sources as well.


The recent achievements in simulating binary black-hole mergers have garnered significant interest from astronomers and astrophysicists alike. Although observations and theoretical investigations of black-hole binaries have been ongoing for an extended period, they have mostly revolved around the context of Newtonian gravity. The recent advancements in numerical relativity have finally provided the missing element: a comprehensive understanding of the impact of the final merger in the strong-field, general relativistic domain. In this segment, we outline three crucial realms of contemporary focus: the phenomenon of recoiling black holes, the spin characteristics of the merged black hole, and the potential electromagnetic signatures originating from the ultimate merger event.

A. Displacement of Black Holes due to Recoil

One of the profound implications arising from gravitational wave computations in relation to astrophysical observations is the potential for gravitational radiation-induced recoil to expel black holes from their resident galaxies. Numerical relativity has predicted recoil velocities that exceed the escape velocities of numerous galaxies, reaching around 1000 km/s. Evaluating the likelihood of such an event demands a comprehensive grasp of recoil’s impact on variables like mass ratios, spin magnitudes, and spin orientations, coupled with assumptions about the distribution patterns of these parameters.

1. Anticipating the Recoil

It would be beneficial to possess a simplified, analytical expression that encapsulates the relevant factors influencing recoil velocity. The intricate nonlinear interactions within merger dynamics govern the ultimate recoil velocity of the amalgamated remnant. Precise results necessitate numerical simulations. Nevertheless, when constructing a heuristic for a phenomenological equation, an initial assumption might mirror the form expected from Post-Newtonian (PN) predictions. This could offer insight into the mass and spin dependencies along with their primary exponents. For instance, PN analysis (Fitchett, 1983) suggests that recoil due to dissimilar masses could be proportional to v = An? 1 – 4n(1 + Bn) (referred to as the “Fitchett formula”), where n signifies the symmetric mass ratio (1).

Furthermore, the spin contributions, at leading PN orders, should stem from components like 4 = 52/M^2 – S1/M and § = S2 + S1 (Campanelli et al., 2007b; Kidder, 1995; Racine et al., 2009).

Such an assumption could be augmented and confined by principles of symmetry. A pragmatic approach involves representing the recoil velocity as a Taylor series encompassing all six spin components (three for each black hole), with coefficients functioning as mass and initial separation-dependent functions. The binary system’s rotation, parity, and exchange symmetries then establish connections among these coefficients (Boyle and Kesden, 2008; Boyle et al., 2008a). For instance, the component of recoil velocity parallel to the orbital angular momentum, through these arguments, can be tied, in spin’s leading order, to the dot product between the black holes’ coordinate separation vector and a specific linear combination of spins (Baker et al., 2008a; Boyle et al., 2008a).

Unknown coefficients within this heuristic can be fitted to numerical outcomes, and corrective terms can be introduced to achieve numerical alignment. Consequently, a tentative equation has been taking shape through collective efforts of the numerical relativity community, assuming the following form (Baker et al., 2007a, 2008a; Brügmann et al., 2008a; Campanelli et al., 2007c; González et al., 2007b; Lousto et al., 2010; Lousto and Zlochower, 2009; van Meter et al., 2010a):

Here, Um stands for the contribution attributed to mass asymmetry, vI for the contribution due to spin that imparts a kick perpendicular to the orbital angular momentum, vi for the contribution due to spin that causes a kick parallel to the orbital angular momentum, a; signifies the projection of black hole i’s dimensionless spin vector a; = S;/M^2 onto the orbital angular momentum, of refers to the angle formed by ∞- with respect to a reference angle in the orbital plane, and 1 and 2 represent constants corresponding to a given mass ratio and initial separation. The spins are assumed to have been measured before or ideally very close to the merger event. 1 and Ф2 account for each spin’s precession degree before the merger.

This formula can then be evaluated across ranges of anticipated mass ratios, spin magnitudes, and various spin orientations to determine the probability of a specific recoil velocity. For example, for mass ratios ranging from 1 to 10 and spin magnitudes of 0.9, the likelihood of exceeding 1000 km/s is projected to be approximately 10%. Studies on speed probability under various model-dependent scenarios are detailed by Schnittman (2007) and Baker et al. (2008a).

2. Ramifications of Black Hole Recoil

Gravitational radiation-induced recoil bears significant astrophysical consequences, notably impacting the growth rates of black holes. For instance, the recoil of massive black holes whose imparted velocity remains below their host galaxies’ escape velocity tends to influence the amount of accreted mass (Blecha and Lob, 2008). Recoil velocities surpassing galactic escape velocities can also affect the growth of massive black holes by diminishing the likelihood of subsequent mergers (Sesana et al., 2007; Volonteri, 2007). This, in turn, may modestly reduce the rate of observable coalescence events for the Laser Interferometer Space Antenna (LISA) (Sesana et al., 2007).

Recoil might also have observable direct consequences. Quasars, often the result of massive black hole coalescence during galactic mergers, could potentially be ejected from their host galaxies. A study by Bonning et al. (2007) found limited evidence of such occurrences, suggesting their rarity. Nonetheless, recent observations of two fast-moving extragalactic quasars (Shields et al., 2009) stand as strong candidates for recoiled black holes (Komossa et al., 2008).

B. Final Black Hole’s Spin

The ultimate spin of the merged outcome of a binary system also holds significance in astrophysics. Understanding the likelihood of different spin strengths aids in creating templates for matched filtering and assessing gravitational signal detectability in projects like LIGO and LISA (Berti et al., 2007a). Being able to predict the end spin based on initial binary properties implies that observing a black hole’s spin could provide insights into its origin. Particularly, grasping the connection between a black hole’s final spin and the orbital angular momentum of its binary predecessor might clarify the formation of X-shaped jets (Barausse and Rezzolla, 2009).

Similarly to the recoil effect, having a simple analytic formula for spin estimation would be beneficial. Insights from numerical simulations hint at the feasibility of this. For instance, in cases of nonspinning, equally massive binaries undergoing circular inspiral, the ultimate spin isn’t influenced by initial distance; instead, it’s governed by the final merger dynamics, leading to a consistent final spin of around 0.7 (Baker et al., 2006a; Campanelli et al., 2006a,b,d; Pretorius, 2005a). Additionally, in the case of dissimilar mass, non-spinning binaries, the end spin roughly scales linearly with the symmetric mass ratio (1) (González et al., 2007b).

Various strategies for constructing a presumption for the final spin have been put forth. In the scenario of aligned-spin, equal-mass mergers, Campanelli et al. (2006d) devised a straightforward quadratic formula involving the shared initial spin that adhered to the “cosmic censorship” idea: af/M < 1. For nonspinning mergers, Berti et al. (2007b) adopted af/Ms = an + br?, where a and b are fitting parameters, aligning well with data from non-spinning binaries.

Lousto et al. (2010) arrived at a more universal assumption for initial black holes of varying spins using the PN approximation. Barausse and Rezzolla (2009) presented an alternative expression for generic binaries, based on assumptions about the near-conservation of magnitudes and relative angles of specific angular momentum vectors. Various other formulas have been proposed (Buonanno et al., 2008; Rezzolla, 2009; Rezzolla et al., 2008a,b; Tichy and Marronetti, 2008), each in good agreement with certain subsets of available numerical data. However, a widely accepted analytical formula for the final spin in the context of generically precessing binaries is still pending.

C. Electromagnetic Signatures of Black Hole Mergers

We’ve previously established that black hole mergers generate strong gravitational-wave signals, but the question remains whether they will also emit detectable photons across the electromagnetic spectrum. In other words, will they be visually observable in addition to being “audible”?

  1. Astrophysical Factors

The presence of observable electromagnetic counterparts hinges on the distribution of gas and magnetic fields around the merging binary. In the case of stellar black hole binaries and intermediate-mass black hole binaries within star clusters, any material in an accretion disk would likely be rapidly consumed by the black holes, making the occurrence of electromagnetic counterparts improbable.

However, the situation differs for massive black hole binaries formed from galaxy mergers. In gas-rich mergers, there could be enough gas available to feed accretion disks around the black holes, leading to the possibility of detectable electromagnetic emissions. Even in gas-poor mergers, the presence of thin, hot gas in elliptical galaxies could still produce electromagnetic signatures.

Such signatures would hold great value for astrophysics. Locating the source in the sky would confirm and characterize the merger, while also offering insights into accretion processes. By measuring the redshift using electromagnetic radiation alongside the luminosity distance from gravitational waves observed by LISA, it could provide an independent calibration of distances and contribute to cosmological studies including the nature of dark energy.

FIG. 40 Magnetic and electric fields lines around inspiralling equal-mass, nonspinning black holes at ~ 40M before the merger, from Palenzuela et al. (2009). The electric field lines are twisted around the black holes, while the magnetic field lines are mainly aligned with the z axis
FIG. 41 Same as Fig. 40, except at a later time, ~ 20M before merger.

2. Simulations with Magnetic Fields or Gas

Efforts to understand electromagnetic signals from black hole mergers have centered around scenarios involving accretion disks. Some approaches explore emission generated when a merged black hole interacts with a surrounding disk. There’s also interest in signals emerging in the dynamic spacetime close to the binary during its final orbits and merger.

Simulations by Palenzuela et al. (2009) have studied the effects of merging black hole binaries on surrounding magnetic fields. For gas-rich mergers, accretion disks around each black hole might evolve into a circumbinary disk, potentially generating detectable electromagnetic emissions. Other scenarios involve solving the equations of general relativistic magnetohydrodynamics in a dynamic spacetime governed by Einstein’s equations. These simulations are ongoing and explore the potential for electromagnetic signatures during the last orbits and merger.

In summary, the presence of electromagnetic counterparts to black hole mergers depends on factors such as the amount of surrounding gas and magnetic fields. Ongoing simulations and studies aim to shed light on the feasibility and characteristics of these electromagnetic signals.

New Horizons and Future Avenues

A mere half-decade ago, uncertainties abounded regarding potential hurdles that might hinder the practical utilization of numerical-relativity simulations for deriving forecasts rooted in general relativity, aimed at tackling inquiries concerning astrophysical black hole binaries. Until then, numerical relativity researchers predominantly directed their efforts toward matters of theoretical and computational physics.

However, quite abruptly, those obstacles dissipated into thin air. As documented, numerical relativity has now gained the potential to decode general relativity’s anticipations concerning the intricate gravitational physics in strong-field environments, and to embark on delving into the realms of astrophysics and gravitational-wave data analysis. Anticipated progress in this arena will not stem from internal issues within gravitational theory, but rather from the realms of astrophysics and other domains where strong-field gravitational theory finds application. These queries will serve as the driving force behind the sustained, more intricate exploration of already unveiled phenomena, as well as the innovation of fresh capabilities to apply simulation findings to novel questions.

A. Gravitational-Wave Astronomical Investigations

While initial strides have been taken to incorporate numerical-relativity findings into the realm of gravitational-wave observations, a significant amount of further effort is necessary. It has become evident that, in numerous potential observations, the signal-to-noise ratio (SNR) from the final stages of the merger will overshadow the inspiral signal. Although numerical-relativity simulations have provided substantial insights into the fundamental characteristics of the burst of gravitational radiation that culminates in the predicted signals, applying these insights to data analysis demands a comprehensive quantitative understanding of the entire spectrum of merger signals across the parameter space.

Although we’ve discussed examples of mergers that explore the binary black hole parameter space along potentially primary axes, obtaining a complete grasp of the signal spectrum necessitates an extensive and systematic investigative endeavor. Even with the support of empirical models designed to encapsulate numerical outcomes, numerous simulations must be conducted to qualify plausible models. Initial models suitable for detecting and potentially estimating parameters of low-SNR signals require substantial refinement for their application to the high-SNR observations anticipated with future instrument upgrades like Advanced LIGO and future systems like LISA. The achievement of these goals will demand a significantly larger number of high-quality simulations.

Beyond gaining a deeper comprehension of signals generated by mergers throughout the majority of the parameter space, researchers in relativity must also expand the scope of parameter space covered by simulations. Current simulations are constrained in terms of the practical range of parameters they can explore. Presently favored initial data models, especially those assuming conformal flatness, are limited to (a/M)1,2 ≤ 0.93. Approaches more appropriate for rapidly spinning Kerr black holes are less developed in various aspects. Although astrophysical constraints on spin magnitudes remain uncertain, treating larger spins may necessitate novel initial data models. Challenges and limitations within simulations may also exist.

Existing simulations are also confined to comparable mass ratios. Simulations involving q ~ 5 and beyond demand considerable computational resources using current techniques. However, astrophysical mergers might occur across a broad range of mass ratios. In situations of extremely large mass ratios, approximations become feasible, treating the motion of the smaller object as a perturbed geodesic within the larger black hole’s spacetime. Analytic methods hold promise under such conditions where numerical simulations are less practical. Nevertheless, a significant range, say 10 < q < 100, might require innovative techniques.

An issue arises from the small spatial scale of the smaller black hole, which mandates exceedingly small timesteps for stable explicit time integrations, despite uneventful occurrences on these timescales. Alternatives like implicit schemes may pave the way for broader applications. Most of the simulations discussed have concentrated on circular or near-circular orbits. Although this is unquestionably a significant segment of the black hole binary population, scenarios have been proposed involving one black hole capturing another after nearly parabolic initial interactions. Unraveling the potential signals from such systems may necessitate new approaches suitable for extended periods of gradual, effectively Newtonian evolution, punctuated by brief episodes of intense gravitational interaction.

B. Diverse Astrophysical Implications

As previously mentioned, the numerical revelation of potent gravitational-wave recoils in asymmetric mergers of spinning black holes has reverberated across various realms of astrophysics, extending beyond the direct observation of gravitational waves.

These findings have ignited a sense of enthusiasm within the black-hole astronomy community, sparking contemplation about other potential applications of numerical relativity. Given that numerous astronomical domains rely mainly on electromagnetic observations, it’s only natural that such intersections stimulate inquiries into the interplay between black holes and matter visible via electromagnetism. A plethora of pertinent phenomena piques interest, including accretion disks encircling black holes that could be perturbed by gravitational recoil, binary systems involving black holes and neutron stars, neutron star-neutron star binaries, the emissions from active galactic nuclei, quasars, and the enigmatic genesis of gamma-ray bursts.

FIG. 42 Numerical simulations are applied to study the critical transition separating mergers from scattering events in high-velocity black-hole encounters (Sperhake et al., 2009). The curves show the path of one black hole in each of three simulations begun with different initial conditions near the critical impact parameter. The trajectories track closely to-get her coming in from the right until the black holes encounter each other near the origin. Figure kindly provided by U. Sper-hake.

A pivotal avenue for forthcoming research involves expanding the scope of physics beyond the realm of purely gravitational simulations discussed thus far. As previously explored, ongoing initiatives are directed towards implementing and advancing magnetohydrodynamics within dynamic spacetime. As these simulations are applied to increasingly realistic astrophysical scenarios, there will be a need to enhance and fine-tune the techniques employed. Effectively resolving shocks and turbulence in conjunction with magnetic fields, ensuring dependable precision amid powerful magnetic fields, and enforcing the magnetic fields’ divergence constraint within a curved spacetime context all pose challenges that necessitate solutions.

C. Exploring Alternative Physics for Data Science

Astrophysical mergers typically revolve around nearly circular orbits or, in more extreme cases, initially parabolic encounters. However, numerical simulations can be harnessed to investigate other forms of interactions involving black holes.

In instances of high-velocity black-hole encounters, an anticipation of a phase transition emerges, demarcating spacetimes in which the black holes merge from those in which they pass each other in a hyperbolic trajectory, forever avoiding close proximity – refer to Fig. 42. Numerical studies are delving into the critical behavior near this phase transition, offering insights into these captivating phenomena existing at the boundaries of powerful gravitational dynamics (Pretorius, 2006; Pretorius and Khurana, 2007; Sperhake et al., 2009). These phenomena have prompted discussions in the context of impending experimental observations.

Within the realm of trans-Planckian energy levels, certain particle collisions might be explicable through classical black-hole dynamics (Banks and Fischler, 1999; Eardley and Giddings, 2002). Speculative theories involving extensive extra spatial dimensions entertain the notion that sub-TeV scale physics could lead to black-hole-like collisions (Giddings, 2001; Giddings and Thomas, 2002). These possibilities serve as a driving force behind the initial numerical-simulation inquiries into ultra-relativistic black-hole collisions (Shibata et al., 2008; Sperhake et al., 2008b, 2009).

D. Strong Gravity and the Science of Observation

Almost a century ago, Einstein bestowed upon us the gift of general relativity. Throughout the subsequent hundred years, this theory, serving as our foundational model of gravitational physics, has successfully withstood experimental and observational trials. These assessments span an array of scales – from controlled physics experiments at the laboratory level and tests within the solar system to examinations of compact celestial entities and cosmological phenomena (Will, 2006). Thus far, these measurements have not necessitated any modifications to Einstein’s gravity theory, although certain models aimed at elucidating cosmic acceleration do venture into alternative gravitational theories (Silvestri and Trodden, 2009).

Today, numerical simulations are divulging the intricate forecasts of Einstein’s theory concerning the conclusive merger dynamics of binary black holes, as well as the signature it leaves behind in emitted gravitational radiation. The advent of gravitational-wave observations capturing such events will thrust these gravitational occurrences into measurements within an intensely strong-gravity environment, far surpassing any previous assessments. These encompass the pioneering evaluations of higher-order terms in the post-Newtonian expansion, along with the actual physical predictions that numerical relativity is presently uncovering.

The prospect that general relativity might encounter its boundaries within these observations ignites a need for an enhanced comprehension of the potential waveforms that could be anticipated by alternative gravitational theories (Yunes and Pretorius, 2009). It’s conceivable that models involving alternative gravity in the context of black hole mergers would also necessitate numerical simulation to formulate waveform predictions, although limited progress has been made in this direction thus far, with an example from Salgado et al. (2008).

Should Einstein’s theory stand as a valid predictor of the intricate physics unveiled by numerical relativity, it would stand as an astounding feat of scientific inference, unraveling the intricacies of phenomena that lie significantly beyond the original empirical observations that formed the basis of the theory. Conversely, if deviations emerge, these observations, paired with a robust comprehension of Einstein’s predictions as facilitated by numerical relativity simulations, might serve as the stepping stones for the subsequent theory of gravity.

In Conclusion

In closing, this review has been crafted through the collaborative efforts of a diverse research community, and its fruition wouldn’t have been possible without the invaluable contributions and assistance of our fellow colleagues. We wish to extend our profound gratitude to a select group of individuals who have made exceptionally significant inputs. William D. Boggs, Bernd Brügmann, Alessandra Buonanno, Mark Hannam, Richard Matzner, Cole Miller, and Harald Pfeiffer have all provided insightful and constructive feedback on our article. Special mention goes to Manuela Campanelli and Harald Pfeiffer, who graciously provided us with figures from their simulations that weren’t accessible in the published literature. Additionally, our discussions with Sean McWilliams have been instrumental in shaping the content of this review.

We must also acknowledge the financial support we’ve received, notably from NASA Grant No. 06-BEFS06–19. Furthermore, B.J.K. benefited from support through an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center. This program is administered by Oak Ridge Associated Universities under contract with NASA.

As we conclude this review, we extend our heartfelt appreciation to all those who have played a role in making this exploration of numerical relativity’s implications on astrophysical black holes a reality.



Dr Francesco Dergano

Founder and Chief Executive Officer (CEO) of SkyDataSol