Merger and Peak Radiation

Near coalescence, an adiabatic orbital description gives way to a short, strongly nonlinear burst and then to damped radiation from the remnant. There is no source-independent instant called “the merger”: the peak of strain amplitude, the peak of luminosity, common-horizon formation, and the onset of a useful ringdown fit are distinct, convention-dependent landmarks. This page defines those landmarks operationally and measures their ordering in one public numerical-relativity waveform.

Helpful background. Einstein Equations supplies the nonlinear field equations; Waves, Green Functions, and Geometric Optics supplies retarded radiation and polarization; Relativistic Two-Body Dynamics supplies the accepted binary parameters and inspiral history; Gravitational Radiation and Waveforms supplies strain modes and flux normalization; and Computational Relativity supplies the numerical-evolution and extraction context.

Operational landmarks near peak radiation

In a source inertial frame, decompose the distance-rescaled strain using unit-normalized spin-weight $-2$ harmonics,

$$H(u,\theta,\phi) =r\bigl(h_+-ih_\times\bigr) =\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell} H_{\ell m}(u)\,{}_{-2}Y_{\ell m}(\theta,\phi).$$

The time convention is $e^{-i\omega u}$. Define $M$ as the initial total mass, $\hat u=u/(GM)$, and $\hat H_{\ell m}=H_{\ell m}/(GM)$. For a declared mode cutoff $\ell_{\max}$, two rotation-invariant scalars are

$$\mathcal A(\hat u) =\left(\sum_{\ell=2}^{\ell_{\max}}\sum_m \lvert\hat H_{\ell m}\rvert^2\right)^{1/2}, \qquad G\mathcal P(\hat u) =\frac{1}{16\pi} \sum_{\ell=2}^{\ell_{\max}}\sum_m \left\lvert\frac{d\hat H_{\ell m}}{d\hat u}\right\rvert^2.$$

Here $\mathcal P$ is the luminosity at future null infinity. We call $u_{\mathcal A}$ the maximum of $\mathcal A$ and $u_{\mathcal P}$ the maximum of $\mathcal P$. These definitions avoid choosing a line of sight, but they still depend on the mode set, retarded-time construction, Bondi frame, interpolation, and numerical waveform. A peak of $\lvert H_{22}\rvert$ would be another legitimate landmark and need not coincide with either one.

Common-horizon formation, denoted $t_{\rm CH}$, is found on a Cauchy slice rather than at null infinity. Comparing its coordinate time with an extrapolated retarded time is therefore a useful product-level diagnostic, not an invariant causal delay. Likewise, ringdown is a regime in which damped remnant modes provide a controlled description. Because $H_{\ell m}$ above is a spherical mode while Kerr quasinormal modes are naturally spheroidal, the projection has the form

$$H_{\ell m}(u)\sim \sum_{\ell',n} \mu_{\ell\ell'm}(a_f\omega_{\ell'mn}) C_{\ell'mn}e^{-i\omega_{\ell'mn}(u-u_0)}, \qquad \operatorname{Im}\omega_{\ell'mn}\lt0,$$

where $\mu_{\ell\ell'm}$ is the spherical–spheroidal overlap; a one-family formula is a weak-mixing truncation, and any mirror branch must be declared separately. Ringdown is not a discontinuous phase beginning at a unique $u_0$. The original effective-one-body construction makes the inspiral-to-plunge transition explicit Buonanno and Damour 1999, §VII; it does not itself supply a complete inspiral–merger–ringdown waveform. Numerical fits show that adding overtones can move a useful ringdown fit start toward the peak of $\lvert H_{22}\rvert$ Giesler et al. 2019, §§II–IV. Neither statement makes a fitted decomposition unique or identifies that single-mode peak with $u_{\mathcal A}$.

Pinned waveform, modes, and frame

The benchmark uses SXS:BBH:0314, catalog tag 3.0.0, accessed 2026-08-18. It is a mildly unequal-mass, nearly circular, nonprecessing binary. The catalog lists the initial mass ratio $q_{\rm init}=m_1/m_2=1.22785$; the highest-resolution reference values are $q_{\rm ref}=1.2278017358$, eccentricity $1.19\times10^{-4}$, and approximately aligned spins $\chi_{1z}=0.3101$, $\chi_{2z}=-0.4599$. The public SXS:BBH:0314 v3.0 data record supplies four Cauchy resolutions, Lev1 through Lev4.

For every resolution the input is Strain_N2.h5: asymptotically extrapolated strain with extrapolation order $N=2$, modes $2\leq\ell\leq8$, and $H/(GM)$ and $u/(GM)$ already dimensionless. The catalog MD5 checksums are d07039260929f361b70473b8b571f0fc (Lev1), 75578609550c88b3a2857d1f522ef8f4 (Lev2), ab4315caa36ff5d7e77b6cd4d82522c3 (Lev3), and cc7244bb8b08c2192f1ec8ed4f8f59d2 (Lev4). The central calculation retains complete multiplets through $\ell=4$.

The loader returns inertial-frame modes. The v3 product metadata record a center-of-mass space translation and boost before the analysis and include the nonlinear memory contribution; those operations reduce center-of-mass drift but do not select a unique frame under the full Bondi–Metzner–Sachs group. The benchmark therefore declares a frame treatment rather than calling the modes frame invariant. SXS catalog waveforms, extraction choices, and resolution comparisons are documented in Boyle et al. 2019, §§2–3 and §4, printed pp. 24–28, especially eqs. (22)–(24).

The SXS files use $G=c=1$ internally. Every tabulated time below is therefore the dimensionless source-frame coordinate $\hat u=u/(GM_{\rm src})$; the symbol $M$ is not used as a time unit. At cosmological redshift $z$, an observed interval is $\Delta t_{\rm det}=(1+z)GM_{\rm src}\Delta\hat u=GM_z\Delta\hat u$, where $M_z=(1+z)M_{\rm src}$. The corresponding observed strain has the scale $GM_z/D_L$; mass ratios and dimensionless spins are unchanged.

Reproducible landmark calculation

The calculation discards all samples with $\hat u\lt1500$, later than every resolution’s catalog relaxation time, to exclude the initial junk burst. It forms $\mathcal A$ and $G\mathcal P$ from all modes through $\ell=4$, differentiates a cubic spline of each complex mode, and locates each maximum on a 2001-point interpolation grid spanning five native samples. Save the following as merger_peak_sxs.py. The script writes the comparison table and diagnostics to standard output after downloading the pinned public inputs.

import numpy as np
import sxs
from scipy.interpolate import CubicSpline

CATALOG, SIMULATION = "3.0.0", "SXS:BBH:0314v3.0"
LEVS, LMAX, T_START = (1, 2, 3, 4), 4, 1500.0
sxs.load("dataframe", tag=CATALOG, download=True, cache=True)

def peak(t, y):
    i = int(np.argmax(y))
    if i < 2 or i + 2 >= len(t):
        raise ValueError("Peak lacks two neighboring samples on each side")
    j = slice(i - 2, i + 3)
    spline = CubicSpline(t[j], y[j])
    fine = np.linspace(t[j][0], t[j][-1], 2001)
    k = int(np.argmax(spline(fine)))
    return float(fine[k]), float(spline(fine[k])), float(t[i]), float(y[i])

def landmarks(lev, ell_max=LMAX):
    sim = sxs.load(f"{SIMULATION}/Lev{lev}", download=True, cache=True)
    h = sim.h
    t, lm, data = np.asarray(h.t), np.asarray(h.LM), np.asarray(h)
    cut, modes = t >= T_START, lm[:, 0] <= ell_max
    t, data = t[cut], data[cut][:, modes]
    amplitude = np.sqrt(np.sum(np.abs(data) ** 2, axis=1))
    dh_dt = CubicSpline(t, data, axis=0).derivative()(t)
    luminosity = np.sum(np.abs(dh_dt) ** 2, axis=1) / (16.0 * np.pi)
    return peak(t, amplitude), peak(t, luminosity), float(sim.metadata.common_horizon_time)

rows = [(lev, *landmarks(lev)) for lev in LEVS]
t_ref, a_ref = rows[-1][1][0:2]
print("Lev  hat_u_A    A_peak       |dhat_u_A|  |dA|/A_ref")
for lev, a, luminosity, horizon in rows:
    print(f"{lev:>3}  {a[0]:9.4f}  {a[1]:.9f}  "
          f"{abs(a[0]-t_ref):10.4f}  {abs(a[1]-a_ref)/a_ref:.3e}")

a4, l4, ch4 = rows[-1][1:]
a8, l8, _ = landmarks(4, ell_max=8)
print(f"Lev4: dhat_u_LA={l4[0]-a4[0]:.4f}, "
      f"G*P_peak={l4[1]:.9f}, dhat_u_CHA={ch4-a4[0]:.4f}")
print(f"raw-grid shifts: dhat_u_A={a4[2]-a4[0]:+.4f}, "
      f"dA/A={(a4[3]-a4[1])/a4[1]:+.3e}")
print(f"ell<=8 shifts: dhat_u_A={a8[0]-a4[0]:+.4f}, "
      f"dA/A={(a8[1]-a4[1])/a4[1]:+.3e}")

The simulation identifier includes v3.0 as well as the catalog tag: catalog metadata and a simulation’s data release are separate version choices. The endpoint check rejects a junk cut or truncated record that puts the apparent maximum outside a valid interpolation neighborhood.

The tested environment is Python 3.12.9, sxs==2026.0.2, numpy==2.2.6, and scipy==1.15.3. From an empty directory containing the saved script, the exact rerun is

python3.12 -m venv sxs-peak-env
sxs-peak-env/bin/python -m pip install sxs==2026.0.2 numpy==2.2.6 scipy==1.15.3
. sxs-peak-env/bin/activate
export SXSCACHEDIR="$PWD/sxs-cache"
python3 merger_peak_sxs.py

Resolution benchmark and landmark ordering

The Lev4 result is the comparison reference, not an assumed continuum value. The absolute peak time retains the simulation’s catalog time origin.

Resolution	Peak-time result $\hat u_{\mathcal A}$	Peak-amplitude result $\mathcal A_{\max}$	$\lvert\Delta\hat u_{\mathcal A}\rvert$ versus Lev4	Fractional amplitude error versus Lev4
Lev1	3677.5052	0.553345462	6.9305	$1.331\times10^{-5}$
Lev2	3686.2951	0.553356568	1.8594	$6.762\times10^{-6}$
Lev3	3681.6531	0.553395014	2.7826	$7.624\times10^{-5}$
Lev4	3684.4357	0.553352827	0	0

The sequence is not monotone. “Lev” labels adaptive numerical-error settings, not a known uniform grid-spacing sequence, so a Richardson order inferred from these four numbers would be spurious. The most resolved pair differs by $\lvert\Delta\hat u_{\mathcal A}\rvert=2.7826$ in accumulated peak time and $7.624\times10^{-5}$ fractionally in peak amplitude; the full dimensionless time spread is $8.7899$. This is exactly why a catalog product must not be treated as an exact continuum waveform.

At Lev4, the three operational landmarks have the following ordering after setting $u_{\mathcal A}=0$:

Landmark	Relative dimensionless time $\Delta\hat u$
Common-horizon coordinate time	$-7.4543$
Strain-norm peak, $\ell\leq4$	$0$
Luminosity peak, $\ell\leq4$	$+6.6519$

The computed peak luminosity is $G\mathcal P_{\max}=1.010493\times10^{-3}$. Its numerical value is conditional on spline differentiation and the mode cutoff; the table establishes an ordering for this product, not a universal chronology for all binaries.

Error budget and physical checks

Resolution spread dominates the measured peak-time variations reported in this diagnostic. Using the nearest native sample instead of the local spline shifts the Lev4 strain peak by $\Delta\hat u=-0.0426$ and its amplitude by $-2.314\times10^{-6}$ fractionally. Extending the sum from $\ell\leq4$ to every available mode through $\ell=8$ shifts the peak by $\Delta\hat u=+0.0016$ and the amplitude by $+5.768\times10^{-5}$. These variations are separate from, and should not be added blindly to, the correlated resolution error.

Several unmeasured systematics remain. Only the $N=2$ extrapolation product is held fixed, so finite-radius and extrapolation-order error is not estimated. A different center-of-mass or BMS transformation can mix modes and alter peaks; residual drift is not a proof of a bad evolution, but it must be diagnosed. Moving the junk cut too early can select the unphysical initial burst. Time alignment can hide accumulated phase error. Omitting weak modes can matter much more for an asymmetric or precessing binary. Differentiation amplifies high-frequency numerical noise, and inconsistent remnant mass or spin would move the candidate quasinormal frequencies. These are limitations of the benchmark, not production-algorithm prescriptions.

Two limiting checks are especially useful. In the early adiabatic regime, $\mathcal A$ and $\mathcal P$ grow slowly compared with an orbit; a large isolated early maximum signals junk or corrupted data. In the stationary-remnant limit, oscillatory ringdown decays and $\mathcal P\to0$. The v3 strain may retain a nonzero memory offset, so demanding $H\to0$ would be the wrong stationarity test.

Exercise: convert the measured delay to detector time

Take the Lev4 luminosity-peak delay $\Delta\hat u=6.6519$ for a source with $M_{\rm src}=60M_\odot$ at $z=0.20$. Find the detector-frame delay. Then convert the Lev3–Lev4 strain-peak difference $\lvert\Delta\hat u\rvert=2.7826$ to the same units and interpret it.

Solution

The redshifted mass is $M_z=(1+z)M_{\rm src}=72M_\odot$. Restoring $c$ for the conversion and using $GM_\odot/c^3=4.92549\,\mu\mathrm{s}$ gives

$$\Delta t_{\rm det} =6.6519\frac{GM_z}{c^3} =6.6519(72)(4.92549\,\mu\mathrm{s}) =2.359\,\mathrm{ms}.$$

The resolution difference becomes

$$2.7826\frac{GM_z}{c^3}=0.987\,\mathrm{ms}.$$

The first number is the detector-frame image of this benchmark’s source-frame landmark separation. The second is not an observational error bar: it is a two-resolution discrepancy in the absolute peak time for one numerical simulation, correlated with its accumulated phase error and time convention.

References

• Boyle, Michael, et al. “The SXS Collaboration Catalog of Binary Black Hole Simulations.” Classical and Quantum Gravity 36 (2019) 195006. DOI. Open PDF
• Buonanno, Alessandra, and Thibault Damour. “Effective One-Body Approach to General Relativistic Two-Body Dynamics.” Physical Review D 59 (1999) 084006. DOI. Open PDF
• Giesler, Matthew, Maximiliano Isi, Mark A. Scheel, and Saul A. Teukolsky. “Black Hole Ringdown: The Importance of Overtones.” Physical Review X 9 (2019) 041060. DOI. Open PDF
• SXS Collaboration, “SXS:BBH:0314,” catalog version 3.0 waveform and metadata record, CaltechDATA (released 2025; accessed 2026-08-18). Data record