Numerical Relativity and Binary Black Hole Mergers

Numerical relativity is the art of solving Einstein's equations on a computer when nothing else will do. Its flagship problem is the coalescence of two black holes: the slow inspiral, the violent merger, and the ringdown of the single hole left behind. In that last handful of orbits the bodies move at half the speed of light, the field is as strong as gravity gets, and there is no symmetry and no small parameter. The post-Newtonian series in $v/c$ has diverged and perturbation theory does not yet apply, so the only way to learn what the spacetime does — and what it radiates — is to integrate the full nonlinear field equations directly. The strategy is to treat gravity as an initial-value problem, exactly as one treats the wave equation or fluid flow. The ADM 3+1 decomposition slices spacetime into a stack of spatial snapshots carrying the 3-metric $\gamma_{ij}$ and extrinsic curvature $K_{ij}$ , and Einstein's equations say how to compute the next snapshot from the current one. But the naive scheme is treacherous: the raw ADM evolution is only weakly hyperbolic, constraint-violating errors explode, and codes crash within milliseconds. The story of numerical relativity is the thirty-year hunt for the reformulations — BSSN, generalized harmonic gauge — the singularity treatments — excision, moving punctures — and the gauge conditions that finally tamed the equations. When they came together in 2005 they produced the merger–ringdown templates that let LIGO recognize GW150914 as two black holes, and the whole waveform-modelling industry that follows.

Key takeaways

Numerical relativity solves the full nonlinear Einstein equations $G_{\mu\nu}=8\pi T_{\mu\nu}$ as a Cauchy problem — foliate spacetime into spatial slices and evolve one forward to compute a merger waveform no perturbative method can reach.
A coalescence has three regimes: an adiabatic inspiral (post-Newtonian theory suffices), a strong-field merger (only full NR captures it), and a ringdown (quasinormal modes of the final Kerr hole).
The bare ADM system is only weakly hyperbolic and unstable; codes use BSSN with moving punctures, or generalized harmonic gauge with excision, to make the evolution strongly hyperbolic and well-posed.
The Hamiltonian and momentum constraints must hold on the initial slice and are preserved analytically; their numerical violation is a simulation's running health monitor.
An equal-mass nonspinning merger radiates about 4.8% of its rest mass and leaves a Kerr hole of spin $a_f\approx 0.69$ , which rings down at frequencies fixed by $M_f$ and $a_f$ alone.
NR waveforms calibrate the effective-one-body, phenomenological, and surrogate models LIGO–Virgo–KAGRA match against data.
Key idea: foliate spacetime into spatial slices, satisfy the constraints on the initial slice, and evolve the induced geometry with a hyperbolic, gauge-fixed reformulation of Einstein's equations — handling the singularities by excision or moving punctures — to compute the merger waveform that no perturbative method can reach.

The core insight

Foliate spacetime into spatial slices, satisfy the constraints on the initial slice, and evolve the induced geometry with a hyperbolic, gauge-fixed reformulation of Einstein's equations — handling the black-hole singularities by excision or moving punctures — to compute the merger waveform that no perturbative method can reach.

Where this lives today

Numerical relativity is now infrastructure. The SXS collaboration, the Einstein Toolkit community, and codes such as SpEC, BAM, and GRChombo have computed thousands of binary waveforms spanning mass ratio and spin — the ground truth for gravitational-wave astronomy. Because a single merger can cost weeks on a supercomputer, the catalogs are not searched directly; instead they calibrate the analytic models — effective-one-body Hamiltonians and fluxes, phenomenological models, and machine-learned surrogates — that generate a template in milliseconds and fill the parameter space LIGO–Virgo–KAGRA must cover. Every parameter estimate for every detected compact-binary event rests on this NR foundation. The frontier moves on: neutron-star mergers demand general-relativistic magnetohydrodynamics and neutrino transport for GW170817-style kilonovae; eccentric and precessing binaries stress the models; and next-generation detectors (Cosmic Explorer, the Einstein Telescope, LISA) will demand waveforms good to a part in $10^4$ across hundreds of cycles, pushing NR toward Cauchy-characteristic extraction and higher-order convergence.

The physical picture

Picture spacetime as a deck of cards, each a three-dimensional snapshot of space at one instant. Numerical relativity holds one card — the geometry "now", stored as $\gamma_{ij}$ and its rate of change $K_{ij}$ — and uses Einstein's equations as a rule for drawing the next card a slice of time later. Two freedoms decide how: the lapse $\alpha$ , how much proper time elapses along the normal between cards, and the shift $\beta^i$ , how the spatial grid slides sideways from one card to the next. Neither changes the physics; both change how the coordinates crawl across the geometry, and choosing them badly is enough to crash a run. The first figure makes this foliation and its defining line element one picture.

Figure

The 3+1 foliation: a spatial slice Σ_t of geometry γ_ij is advanced to Σ_(t+dt) by stepping a proper time α dt along the unit normal n and sliding the coordinate grid by the shift β^i dt, so the same coordinate label x^i traces the blue worldline while the physics is carried by the dashed normal

The black holes are the hard part: at the centre of each the geometry runs to a singularity no coordinate patch can hold. Two tricks make it harmless. In excision one cuts a region out of the grid inside the horizon and evolves only the exterior — legitimate, since nothing inside a horizon can reach the outside physics. In the puncture method one factors the singular $1/r$ piece out analytically and evolves the smooth remainder; upgraded to moving punctures, the singular points are advected bodily through the grid, and with the right slicing the code never touches them. The second figure shows the payoff: the waveform of a real merger, the number the whole apparatus exists to produce.

Figure

Gravitational-wave strain of a binary black hole coalescence in the style of GW150914: a chirping inspiral of rising frequency and amplitude, a peak at merger, then a rapidly damped ringdown — source masses about 36 M☉ and 29 M☉ merge to a 62 M☉ Kerr hole, radiating roughly 3 M☉c² as the wave frequency sweeps from about 35 Hz to about 250 Hz

The wave itself is a faint ripple on a vastly stronger background curvature, so extraction is delicate: go far from the source, project the Weyl curvature onto an outgoing null tetrad to isolate its radiative component $\psi_4$ , and integrate twice in retarded time to recover $h_+ - i\,h_\times$ . That number in time is the chirp a LIGO sonification plays — rising in pitch and volume through the inspiral, peaking at merger, dying in a couple of ringdown oscillations. Why did it take three decades? Because the same equations that carry gentle waves also describe horizons forming and coordinates threatening to collapse; a careless discretization amplifies constraint-violating noise exponentially.

The core equations

Everything begins with the ADM split of the line element, trading the single four-metric for a lapse, a shift, and a spatial metric on each slice — the change of variables that turns "solve Einstein's equations" into "evolve a geometry forward in time".

\mathrm{d}s^2 = -\alpha^2\,\mathrm{d}t^2 + \gamma_{ij}\left(\mathrm{d}x^i + \beta^i\,\mathrm{d}t\right)\left(\mathrm{d}x^j + \beta^j\,\mathrm{d}t\right) \qquad (1)

where: $\alpha$ — lapse, the proper time per unit coordinate time along the normal; $\beta^i$ — shift, how the spatial coordinates are dragged slice to slice; $\gamma_{ij}$ — induced 3-metric on $\Sigma_t$ ; $x^i,\,t$ — spatial coordinates and the foliation time label. Reading it: lapse and shift are pure gauge — only the labelling of events — while all geometry lives in $\gamma_{ij}$ . Projecting the ten Einstein equations onto and across a slice separates them into four carrying no time derivative of $K_{ij}$ and six that do. The four are the constraints — elliptic conditions any valid slice must satisfy.

\mathcal{H} \equiv {}^{(3)}R + K^2 - K_{ij}K^{ij} - 16\pi\rho = 0,\qquad \mathcal{M}^i \equiv D_j\left(K^{ij} - \gamma^{ij}K\right) - 8\pi S^i = 0 \qquad (2)

where: ${}^{(3)}R$ — spatial Ricci scalar; $K_{ij},\,K=\gamma^{ij}K_{ij}$ — extrinsic curvature and its trace; $D_j$ — covariant derivative of $\gamma_{ij}$ ; $\rho,\,S^i$ — energy and momentum density seen by the normal observer, both zero in vacuum. Reading it: because these hold on every slice, the initial data are not free — you must solve $\mathcal{H}=\mathcal{M}^i=0$ before starting. The remaining six equations advance the metric in time. The first pair says $K_{ij}$ is essentially the velocity of the 3-metric, dragged by the shift.

\partial_t \gamma_{ij} = -2\alpha K_{ij} + D_i\beta_j + D_j\beta_i \qquad (3)

where: $\partial_t$ — coordinate time derivative at fixed $x^i$ ; $-2\alpha K_{ij}$ — geometric change as the slice bends; $D_{(i}\beta_{j)}$ — Lie drag of the metric by the shift. Reading it: with zero shift the metric changes purely because the slice is extrinsically curved. The second pair is the dynamical heart — the analogue of "acceleration equals force" for geometry, sourced even in vacuum by the spatial Ricci tensor and the lapse gradient.

\partial_t K_{ij} = -D_iD_j\alpha + \alpha\!\left({}^{(3)}R_{ij} - 2K_{ik}K^k{}_j + K K_{ij}\right) - 8\pi\alpha\!\left(S_{ij} - \tfrac{1}{2}\gamma_{ij}(S-\rho)\right) + \beta^k D_k K_{ij} + K_{ik}D_j\beta^k + K_{kj}D_i\beta^k \qquad (4)

where: ${}^{(3)}R_{ij}$ — spatial Ricci tensor; $S_{ij},\,S$ — spatial stress and trace (zero in vacuum); $D_iD_j\alpha$ — acceleration from the lapse gradient, driving the dynamics even with no matter. Reading it: this is where curvature propagates, and its principal part in ${}^{(3)}R_{ij}$ is exactly what makes bare ADM only weakly hyperbolic. The cure behind moving-puncture binaries is the BSSN change of variables: peel off a conformal factor, split the trace of $K_{ij}$ , and promote the contracted connection to an independent field.

\tilde\gamma_{ij} = e^{-4\phi}\gamma_{ij},\quad e^{4\phi}=\det(\gamma_{ij})^{1/3};\qquad \tilde A_{ij} = e^{-4\phi}\!\left(K_{ij}-\tfrac{1}{3}\gamma_{ij}K\right),\qquad \tilde\Gamma^i = -\partial_j \tilde\gamma^{ij} \qquad (5)

where: $\phi$ — conformal factor; $\tilde\gamma_{ij}$ — conformal metric of unit determinant; $\tilde A_{ij}$ — conformal trace-free extrinsic curvature; $K$ — trace of $K_{ij}$ , evolved separately; $\tilde\Gamma^i$ — conformal connection functions, promoted to independent variables. Reading it: promoting $\tilde\Gamma^i$ is the decisive step — it turns the Ricci principal part into a clean wave operator (Derivation 2). An independent route to well-posedness is the generalized harmonic formulation Pretorius used in 2005: impose wave equations on the coordinates, so every metric component obeys a manifestly hyperbolic equation.

\Box x^\mu = \gamma^{\alpha\beta}\,{}^{(4)}\Gamma^\mu{}_{\alpha\beta} = -H^\mu,\qquad g^{\alpha\beta}\partial_\alpha\partial_\beta g_{\mu\nu} = \text{(lower-order)} - 2\partial_{(\mu}H_{\nu)} \qquad (6)

where: $H^\mu$ — gauge source functions ( $H^\mu=0$ is pure harmonic gauge); $\Box$ — curved-space d'Alembertian on the coordinate functions; $g_{\mu\nu}$ — full 4-metric. Reading it: $H^\mu$ steers the gauge while the principal part is guaranteed to be the wave operator $g^{\alpha\beta}\partial_\alpha\partial_\beta$ — a symmetric-hyperbolic system. Whatever the formulation, the gauge must keep the coordinates alive through merger. The moving-puncture recipe pairs $1+\log$ slicing with a $\Gamma$ -driver shift.

\partial_t\alpha = -2\alpha K,\qquad \partial_t\beta^i = \tfrac{3}{4}B^i,\quad \partial_t B^i = \partial_t\tilde\Gamma^i - \eta B^i \qquad (7)

where: $\alpha$ — lapse, driven to collapse near a puncture so the slicing avoids the singularity; $B^i$ — auxiliary shift; $\eta$ — damping of order $1/M$ that stops shift drift. Reading it: the collapsing lapse freezes evolution deep inside the horizon while the $\Gamma$ -driver lets punctures coast across the grid. Finally, the observable: the strain is read from the Weyl scalar $\psi_4$ far from the source, where the radiative part of the curvature falls off most slowly and dominates.

\psi_4 = -C_{\alpha\beta\gamma\delta}\,n^\alpha \bar m^\beta n^\gamma \bar m^\delta,\qquad \lim_{r\to\infty}\psi_4 = \ddot h_+ - i\,\ddot h_\times \qquad (8)

where: $C_{\alpha\beta\gamma\delta}$ — Weyl conformal curvature tensor; $n^\alpha,\,\bar m^\beta$ — legs of a null tetrad adapted to the outgoing radiation; $h_+,\,h_\times$ — the two polarization strains. Reading it: expand $\psi_4$ in spin-weight $-2$ harmonics, keep the dominant $(2,2)$ mode, integrate twice in retarded time for the template.

Derivation ASplitting Einstein's equations into constraints and evolution

The goal is to see why ten covariant equations become four constraints plus six evolution equations, and why the four contain no second time derivative — the fact that makes them conditions on initial data rather than dynamical laws. We start from $G_{\mu\nu}=8\pi T_{\mu\nu}$ , introduce the unit normal $n^\mu$ to a slice, and project the equations into normal and tangential parts using the Gauss and Codazzi relations, tracking exactly where time derivatives survive.

Step 1. Take the future-pointing unit normal $n^\mu$ to $\Sigma_t$ and the projector $\gamma^\mu{}_\nu = \delta^\mu{}_\nu + n^\mu n_\nu$ ; the extrinsic curvature $K_{ij}=-\gamma^\mu{}_i\gamma^\nu{}_j\nabla_\mu n_\nu$ measures how the normal tips along the slice. Why: this machinery expresses four-dimensional curvature in terms of a slice's intrinsic curvature and how it sits in spacetime.

Step 2. Contract the Einstein equation twice with $n^\mu$ . The Gauss relation gives ${}^{(3)}R + K^2 - K_{ij}K^{ij} = 16\pi\rho$ — the Hamiltonian constraint, with $\rho = n^\mu n^\nu T_{\mu\nu}$ . Why: the twice-normal projection kills every second time derivative, so it constrains data on one slice rather than evolving it.

Step 3. Contract once normal and once tangentially. The Codazzi relation gives $D_j(K^{ij}-\gamma^{ij}K)=8\pi S^i$ — the momentum constraint, with $S^i = -\gamma^{i}{}_\mu n_\nu T^{\mu\nu}$ . Why: again no second time derivative survives, so together these four equations pin down which $(\gamma_{ij},K_{ij})$ are admissible initial data.

Step 4. The purely tangential projection keeps the second time derivatives and yields the six evolution equations (3)–(4), lapse and shift free. The contracted Bianchi identity then makes $\partial_t\mathcal{H}$ and $\partial_t\mathcal{M}^i$ vanish whenever constraints and evolution hold — the constraints propagate. Why: this is the licence to enforce the constraints once on the initial slice and thereafter read their violation as pure discretization error. *Verified: the projections reproduce (2)–(4); constraint propagation follows from $\nabla_\mu G^{\mu\nu}=0$ . MTW ch. 21, Wald ch. 10.*

Go deeper · Derivation B · 1 minWhy BSSN is well-posed where ADM crashes

Derivation BWhy BSSN is well-posed where ADM crashes

The aim is to expose, at the level of the principal part, why the bare ADM evolution is only weakly hyperbolic — so small errors need not stay bounded — and how the BSSN variables repair it into a strongly hyperbolic system. Well-posedness is no nicety: a weakly hyperbolic system can admit solutions that grow without bound as the grid is refined, so no resolution saves it. The character of the equations is fixed by the second-derivative terms in the Ricci tensor.

Step 1. Write the spatial Ricci tensor schematically as ${}^{(3)}R_{ij} = -\tfrac{1}{2}\gamma^{kl}\partial_k\partial_l\gamma_{ij} + \partial_{(i}\Gamma_{j)} + \dots$ , with $\Gamma^i=\gamma^{kl}\Gamma^i{}_{kl}$ . The first term is a clean wave operator; the mixed second derivatives in $\partial\Gamma$ spoil it. Why: the type of a PDE system is fixed by its highest derivatives, so these mixed terms decide stability.

Step 2. In ADM those terms stay in the principal part and the system is only weakly hyperbolic: its characteristic matrix cannot be diagonalized, admitting Jordan blocks whose modes grow polynomially and, coupled nonlinearly, without bound. Why: a non-diagonalizable principal symbol is the algebraic signature of an ill-behaved initial-value problem.

Step 3. Promote $\tilde\Gamma^i$ to an independent field with its own evolution equation, substituting it wherever $\partial_j\tilde\gamma^{ij}$ appears. The mixed derivatives are absorbed, and the conformal metric obeys $\partial_t^2\tilde\gamma_{ij}\sim \alpha^2\,\tilde\gamma^{kl}\partial_k\partial_l\tilde\gamma_{ij}$ — a clean wave operator. Why: trading a second metric derivative for a first derivative of an auxiliary field is the standard route from weak to strong hyperbolicity.

Step 4. Close with $1+\log$ slicing and the $\Gamma$ -driver shift (7); with the conformal-traceless split of $K_{ij}$ this is the moving-puncture BSSN system that runs stably through merger and ringdown on ordinary grids. Why: a well-posed continuum problem is the necessary condition for a convergent numerical scheme. *Verified: promoting $\tilde\Gamma^i$ yields a wave-type principal part; moving-puncture BSSN is strongly hyperbolic. Baumgarte–Shapiro; Alcubierre.*

Go deeper · Derivation C · 1 minQuasinormal ringdown and reading off the remnant

Derivation CQuasinormal ringdown and reading off the remnant

The purpose is to justify why the post-merger signal is a sum of exponentially damped sinusoids whose frequencies scale as $1/M_f$ and depend only on the remnant spin, and why measuring them tests the no-hair theorem. After merger the spacetime is a single, slightly distorted Kerr hole, so linear perturbation theory applies, and its natural oscillations are the quasinormal modes of the Teukolsky equation with outgoing conditions at infinity and ingoing conditions at the horizon.

Step 1. Linearize about Kerr. The perturbation separates into spin-weight $-2$ spheroidal harmonics $(\ell,m)$ , and the dissipative boundary conditions permit only a discrete set of complex frequencies $\omega_{\ell m n} = \omega_R - i/\tau$ , indexed by overtone $n$ . Why: those conditions make the problem non-Hermitian, so eigenfrequencies are complex and the modes decay.

Step 2. The only scale is $M_f$ , so every frequency is $\omega_{\ell m n} = \hat\omega_{\ell m n}(a_f)/M_f$ with $a_f=J/M_f^2$ a pure number; thus $\omega_R M_f$ and $Q=\tfrac12\omega_R\tau$ depend on spin alone. Why: in geometric units mass is a Kerr hole's only dimensional parameter, so it sets the scale, never the dimensionless shape.

Step 3. The strain is $h_{\text{ring}}(t)\propto \sum_{\ell m n} A_{\ell m n}\,e^{-t/\tau_{\ell m n}}\cos(\omega_{R,\ell m n}\,t+\phi_{\ell m n})$ , dominated by $(2,2,0)$ . Perturbation theory fixes every $\omega$ and $\tau$ but not the amplitudes $A$ and phases $\phi$ — those depend on how the merger excited the hole, and only NR supplies them. Why: this division of labour is why NR is indispensable even where perturbation theory applies.

Step 4. Since two parameters $(M_f,a_f)$ fix all the frequencies, one measured mode determines them and a second over-determines the system — a consistency test of Kerr-ness. Why: this is the operational content of the no-hair theorem, which NR both predicts and validates. *Verified: Schwarzschild $(2,2,0)$ $M\omega \approx 0.37367 - 0.08896\,i$ (Leaver 1985), giving $Q\approx 2.10$ ; the $1/M_f$ scaling is exact.*

Limiting cases

The three regimes of a coalescence, and the two mass-ratio extremes, show precisely where full numerical relativity is the only tool.

Case	Result	Meaning
Inspiral (wide separation)	$v/c \ll 1,\ r \gg M$	Adiabatic, weak field: post-Newtonian theory gives the waveform analytically. NR cannot afford the thousands of early orbits, so PN and NR are hybridized at the boundary.
Merger (strong field)	$v/c \sim 0.5,\ r \sim M$	Highly nonlinear, no symmetry, no small parameter. Only full NR captures the last orbits, plunge, and horizon formation.
Ringdown (post-merger)	$h \sim \sum_{\ell m n} A_{\ell m n}\, e^{-t/\tau_{\ell m n}}\cos(\omega_{\ell m n} t + \phi_{\ell m n})$	The remnant settles to Kerr by radiating quasinormal modes; perturbation theory fixes the frequencies, NR the amplitudes, confirming $M_f,a_f$ .
Extreme mass ratio	$q = m_2/m_1 \to 0$	Perturbation theory and gravitational self-force take over; 3D NR chokes on the disparate length scales.
Test-field / linear limit	$\|h\|\to 0$	Curvature perturbation on a fixed background — linearized gravity and the Teukolsky equation, recovered as NR's weak-field limit.

Common trapA stable-looking run is automatically a correct one.

It is not. First, hyperbolicity: the bare ADM equations are only weakly hyperbolic, so constraint-violating modes need not stay bounded and can blow up in milliseconds — a reformulation (BSSN or generalized harmonic) is mandatory, not cosmetic. Second, subtler: a run can stay numerically stable while quietly drifting off the constraint surface or the correct gauge, producing a smooth but wrong waveform. So a real NR result is never just "the code finished" — one must show the constraints stay small, demonstrate convergence at the expected finite-difference order, and check that $\psi_4$ is independent of extraction radius after extrapolation to infinity. Well-posedness, convergence, and constraint monitoring are as much physics as the field equations.

Worked examples

1 · Energy radiated and final spin of an equal-mass nonspinning merger

Two nonspinning black holes of equal mass spiral together from a quasicircular orbit. Estimate the fraction of the total mass radiated as gravitational waves and the dimensionless spin of the Kerr remnant, working in geometric units $G=c=1$ .

Step 1. Fix the bookkeeping. With $M=m_1+m_2$ and $m_1=m_2$ , the symmetric mass ratio is $\nu = m_1 m_2/M^2 = 1/4$ , its maximum. The remnant mass is $M_f = M - E_{\text{rad}}$ , with $E_{\text{rad}}$ read from the flux $\mathrm{d}E/\mathrm{d}t = \tfrac{1}{16\pi}\oint |\dot h|^2\,\mathrm{d}\Omega$ integrated over the coalescence.

Step 2. Read the radiated energy. Long convergent simulations (SXS, Einstein Toolkit) give $E_{\text{rad}}/M \approx 0.048$ — about 4.8% of the rest mass leaves as gravitational waves, most of it from the final orbits and merger where the $(2,2)$ amplitude peaks.

Step 3. Read the spin from the final horizon. The holes bring in orbital angular momentum, but the merger cannot retain it all — the excess is radiated, leaving a single Kerr hole with $a_f = J_f/M_f^2 \approx 0.69$ , below the extremal bound $a_f=1$ .

Step 4. Sanity-check the efficiency. Hydrogen fusion liberates $\sim 0.007$ of rest mass; here $\approx 0.048$ leaves as radiation in a fraction of a second, at a peak luminosity of order $c^5/G$ times a dimensionless factor — briefly the most luminous event in the universe. The headline numbers are $E_{\text{rad}}/M\approx 4.8\%$ and $a_f\approx 0.69$ , reproduced by independent codes.

2 · Ringdown frequency and damping time of the GW150914 remnant

The GW150914 coalescence left a Kerr remnant of mass $M_f\approx 62\,M_\odot$ and spin $a_f\approx 0.67$ . Predict the frequency and damping time of its dominant $(2,2,0)$ quasinormal mode, and compare with the observed post-merger signal.

Step 1. Convert the mass to a time. Since $GM_\odot/c^3 = 4.925\times10^{-6}\,\text{s}$ , the remnant timescale is $M_f = 62\times 4.925\times10^{-6}\,\text{s} = 3.05\times10^{-4}\,\text{s}$ ; every frequency comes out as a pure number over this.

Step 2. Use the fitted dimensionless frequency. A standard fit to the Teukolsky spectrum gives $M_f\,\omega_R = 1.5251 - 1.1568\,(1-a_f)^{0.1292}$ ; with $a_f=0.67$ this is $\approx 0.523$ , so $\omega_R \approx 0.523/(3.05\times10^{-4}\,\text{s}) \approx 1.71\times10^{3}\,\text{rad/s}$ .

Step 3. Convert to a frequency: $f = \omega_R/2\pi \approx \mathbf{2.7\times10^{2}\ Hz}$ — about 270 Hz, squarely in the band of the observed GW150914 post-merger signal near 250 Hz.

Step 4. Get the damping time from the quality factor $Q = \tfrac12\omega_R\tau \approx 0.700 + 1.419\,(1-a_f)^{-0.499} \approx 3.2$ , so $\tau = 2Q/\omega_R \approx \mathbf{3.7\ ms}$ . The mode completes only $Q/\pi \approx 1$ e-folding cycle before fading — a ~270 Hz ringdown lasting a few milliseconds, the brief final "note" of the detection. Measuring $f$ and $\tau$ returns $M_f$ and $a_f$ , closing the loop with the inspiral masses.

Physics of the portrait

Numerical relativity is a fifty-year story of an equation that resisted being solved — from the founders of the initial-value formulation to the person who first watched two black holes merge inside a computer.

Origins · 1 minFrans Pretorius and how this was discovered

Frans Pretorius

b. 1973 · first stable binary black hole merger, 2005

Frans Pretorius, a South African–Canadian physicist trained under Matthew Choptuik, ended a thirty-year impasse in 2005 with the first stable, long-term evolution of two black holes through inspiral, merger, and ringdown. The problem was old: Arnowitt, Deser and Misner cast general relativity as a Hamiltonian initial-value problem in 1962, and Hahn and Lindquist attempted a two-hole evolution in 1964, but every direct attack on the ADM equations crashed — the system is weakly hyperbolic, and constraint-violating noise grew explosively however fine the grid. Smarr and Eppley computed the first head-on collision waveform in the late 1970s; Choptuik found critical phenomena in collapse in 1993; and through the 1990s the Shibata–Nakamura and Baumgarte–Shapiro conformal reformulation (BSSN) tamed single-star evolutions. Pretorius took a different road: a first-order generalized harmonic formulation, coordinates satisfying wave equations, combined with excision, constraint damping, and adaptive mesh refinement. Within a year the Campanelli–Lousto–Marronetti–Zlochower and Baker–Centrella groups independently found the moving-puncture method, achieving the same inside the simpler BSSN framework, which became the workhorse. Together these 2005–2006 breakthroughs turned general relativity's hardest dynamical problem into a routine computation — and, a decade later, supplied the merger–ringdown templates without which LIGO could not have identified GW150914 as two black holes.

History

1962

— Arnowitt, Deser and Misner cast general relativity as a Hamiltonian initial-value problem (the ADM 3+1 split)

1964

— Hahn and Lindquist attempt the first numerical two-black-hole evolution

1977

— Smarr and Eppley compute the first head-on black-hole collision waveform

2005

— Pretorius achieves the first stable binary black hole merger via generalized harmonic gauge

2006

— Campanelli et al. and Baker et al. independently discover the moving-puncture method

2015

— GW150914: numerical-relativity waveforms match the first detected black-hole merger

Going further

Texts. Baumgarte and Shapiro, Numerical Relativity (2010) — the standard text · Alcubierre, Introduction to 3+1 Numerical Relativity (2008) · MTW, Gravitation, ch. 21 · Wald, General Relativity, ch. 10 · Carroll, Spacetime and Geometry, ch. 4 · Schutz, A First Course in General Relativity · Hartle, Gravity.

Papers. Pretorius, Phys. Rev. Lett. 95, 121101 (2005) — first stable merger · Campanelli et al., PRL 96, 111101 (2006) and Baker et al., PRL 96, 111102 (2006) — moving punctures · Leaver, Proc. R. Soc. Lond. A 402, 285 (1985) — quasinormal modes · Abbott et al., PRL 116, 061102 (2016) — GW150914.

Provenance. The ADM split and constraints follow MTW ch. 21 and Wald ch. 10 (rederive from Gauss–Codazzi); BSSN and hyperbolicity follow Baumgarte–Shapiro and Alcubierre; generalized harmonic gauge follows Pretorius (2005). The equal-mass nonspinning $E_{\text{rad}}/M\approx 0.048$ and $a_f\approx 0.69$ are standard (Campanelli et al. 2006; Baker et al. 2006; SXS); the Schwarzschild $(2,2,0)$ QNM $M\omega\approx 0.37367 - 0.08896\,i$ is Leaver (1985); the GW150914 remnant values ( $M_f\approx 62\,M_\odot$ , $a_f\approx 0.67$ , $f_{\text{ring}}\approx 270\,\text{Hz}$ , $\tau\approx 3.7\,\text{ms}$ ) use the Berti–Cardoso–Will fits with $GM_\odot/c^3=4.925\times10^{-6}\,\text{s}$ . Geometric units $G=c=1$ .

Prerequisites. ADM 3+1 Decomposition · Hamiltonian and Momentum Constraints · Initial-Value (Cauchy) Problem · BSSN Formulation · Generalized Harmonic Gauge. Leads to. Black Hole Excision and Puncture Methods · Apparent Horizons and Horizon Finders · Gauge Conditions in Numerical Relativity · Conformal Method for Constraint Solving.

Ask the tutor: Ask me to walk through the Gauss–Codazzi projection that turns Einstein's equations into constraints plus evolution, or to estimate the ringdown note of any remnant mass and spin.

Why can't post-Newtonian theory just be pushed to cover the merger?

The post-Newtonian expansion is a series in $v/c$ , valid while the bodies move slowly and stay far apart. In the final orbits the orbital speed passes half the speed of light and the separation shrinks to a few times the total mass, so the expansion parameter is no longer small and the series stops converging usefully. Worse, the merger involves horizons forming and joining — a genuinely non-perturbative event with no counterpart in the slow-motion expansion. PN theory stays excellent for the long early inspiral and is stitched to NR at the boundary, but the plunge and merger can only be obtained by integrating the full equations.

If the lapse and shift are pure gauge, why do they matter so much?

Because gauge choices, though physically irrelevant, are numerically decisive. The lapse and shift do not change the spacetime, but they change how the coordinate grid crawls across it — and a grid that piles up, or marches into a singularity, destroys the computation regardless of the physics. The $1+\log$ slicing drives the lapse toward zero near a puncture, freezing evolution where the geometry is pathological, while the $\Gamma$ -driver shift lets punctures glide through the grid instead of being smeared across it. Good gauge is what keeps the coordinates alive long enough to reach the answer.

What exactly is a puncture, and how can the code evolve a singularity?

A puncture is a point where the conformal factor blows up, encoding a black hole's interior singularity in a controlled analytic form. Because that singular piece is known in closed form it can be factored out: the computer evolves only the regular remainder, never the divergent part. In the moving-puncture method the singular points are advected through the grid like tracked particles, and the collapsing lapse ensures the region near each puncture is barely evolved. Combined with the fact that nothing inside a horizon can influence the exterior, this lets a finite-difference code carry black holes for hundreds of orbits without ever resolving the true singularity.

How is the strain extracted, and why far from the binary?

One computes the Weyl tensor on a large coordinate sphere, projects it onto a null tetrad adapted to outgoing light to isolate $\psi_4$ , expands in spin-weight $-2$ harmonics, and integrates the dominant modes twice in retarded time to get $h_+ - i\,h_\times$ . It must be done far out because only in the wave zone does the radiative $1/r$ part of the curvature cleanly dominate the near-zone Coulombic and gauge pieces that fall faster — the peeling property. In practice one extracts at several radii and extrapolates to infinity, or uses Cauchy-characteristic extraction to reach null infinity, which matters for the part-in- $10^4$ accuracy next-generation detectors will demand.

What does constraint violation tell you during a run?

On the continuum the constraints, once satisfied on the initial slice, hold for all time — the evolution preserves them exactly. On a discrete grid they do not: truncation error seeds a small violation every step. Monitoring $\mathcal{H}$ and $\mathcal{M}^i$ therefore gives a direct measure of numerical error, one that should shrink at the expected order as the grid is refined. If instead it grows exponentially, the formulation is ill-posed or the gauge has failed; constraint-damping terms are added precisely to drive small violations back toward zero rather than let them amplify.

Why does the ringdown last only a couple of oscillations?

Because the fundamental quasinormal mode of a moderately spinning hole has a low quality factor — for the GW150914 remnant, $Q\approx 3$ . The quality factor counts oscillations per e-folding of amplitude, so $Q/\pi\approx 1$ means the mode loses most of its energy within about one cycle. A perturbed black hole is a very leaky resonator: it sheds its deformation into gravitational waves so efficiently that the ringdown is a brief few-millisecond chirp, not a long tone. Higher-spin remnants ring a little longer and overtones decay even faster, which is why detecting a second mode — black-hole spectroscopy — is so demanding.