Numerical-Relativity Waveform Products

A reusable numerical-relativity (NR) waveform is not merely a plotted strain curve. It is a versioned product: simulation identity, physical parameters, resolutions, extraction method, asymptotic frame, waveform quantity, harmonic convention, junk-radiation treatment, and checksums must travel together. This page defines that product and tests one public example. Evolution and extraction algorithms are treated in Computational Relativity; here the concern is whether their output is sufficiently specified and tested to enter approximants, hybrids, and surrogates. Detector-specific adequacy remains a question for Gravitational-Wave Astronomy.

Helpful background. Einstein Equations supplies the evolved field equations; Waves, Green Functions, and Geometric Optics supplies outgoing-wave and Fourier concepts; Relativistic Two-Body Dynamics supplies source parameters and balance quantities; Computational Relativity supplies production evolution and extraction algorithms; and Gravitational-Wave Data Analysis and Inference supplies application-specific comparison criteria.

What the product must identify

For each simulation, record the catalog and release, simulation identifier, numerical resolution labels, initial and reference masses and spins, eccentricity, reference time, and remnant metadata. For each waveform, additionally record

• whether the quantity is strain $h$, news, or a convention-defined $\Psi_4$;
• the retarded-time and radius normalization, included $(\ell,m)$ modes, and harmonic phases;
• finite-radius, extrapolated, or characteristic extraction and all associated orders or worldtubes;
• Bondi, asymptotically inertial, center-of-mass, source, or coprecessing frame information;
• treatment of junk radiation, offsets, memory, drift, and missing early or late time;
• file version, checksum, access date, and software environment.

The 2019 SXS catalog describes its metadata, finite-radius products, polynomial extrapolation, and center-of-mass correction in Boyle et al. 2019, §§2.2–2.4, pp. 9–19, PDF. Its preferred strain file stores nonuniform samples of $u/M$, $\operatorname{Re}(rh_{\ell m}/M)$, and $\operatorname{Im}(rh_{\ell m}/M)$ in groups for extrapolation orders $N=2,3,4$ Boyle et al. 2019, Appendix A.3.1 and Table A2, pp. 37–39, PDF.

This page uses $h=h_+-ih_\times$ and unit-normalized spin-weight $-2$ harmonics. The SXS strain product uses that complex polarization, the $(-,+,+,+)$ signature, and an inertial angular basis. Its curvature file, however, uses $r\Psi_4=-r\ddot h$ Boyle et al. 2019, Appendix C, Eq. (C.24), pp. 41–44, PDF, whereas this volume’s inherited curvature convention has $\Psi_4^0=\ddot H$. A curvature product must therefore be translated as a complete convention package; the benchmark below deliberately uses strain only.

Extrapolation, characteristic extraction, and frames

Finite-radius modes contain near-zone, gauge, and retarded-time effects. Polynomial extrapolation fits the waveform’s slowly varying amplitude and phase in inverse radius and then returns to an inertial frame. Increasing the polynomial order need not improve the result: early inspiral may favor higher order while merger can amplify short-wavelength numerical noise. The SXS catalog consequently supplies several orders and recommends varying them Boyle et al. 2019, §2.4.1, pp. 13–15, PDF.

Characteristic extraction instead propagates data from a finite worldtube to future null infinity using a characteristic evolution. Its independent controls include worldtube radius, characteristic resolution, initial null data, and Bondi-frame fixing. “At null infinity” alone does not make two products equivalent: residual BMS transformations can alter cuts, modes, memory offsets, and charges.

The SXS filename ending in CoM denotes a fitted translation and boost intended to suppress coordinate center-of-mass drift Boyle et al. 2019, §2.4.2, pp. 15–19, PDF. It is not a claim that every proper-supertranslation ambiguity has been fixed. A hybrid or surrogate must name the exact frame it consumes rather than silently aligning away one visible discrepancy.

A pinned public-product test

The reproducible benchmark uses SXS:BBH:0305, the fifth Zenodo deposit (v5), record 3301877, accessed 2026-08-18 SXS Collaboration, versioned dataset. Inputs are metadata for Lev2–Lev6 and the Lev6 file rhOverM_Asymptotic_GeometricUnits_CoM.h5, MD5 7e7cc4b5a929d92ca49d2a60af19f4de. The file contains 77 modes through $\ell=8$ on $-107.8674\le u/M\le3844.4222$ and the $N=2,3,4$ extrapolated groups. Its top-level VersionHist.ver dataset has the single row commit ef51849550a1d8a5bbdd810c7debf0dd839e86dd and message Overall sign change in complex strain h. The rerun asserts that exact post-sign-change product instead of inferring a strain sign from the filename. The environment is Python 3.9.6, NumPy==2.0.2, and h5py==3.11.0.

The first check compares remnant masses carried by independently evolved resolutions. Lev4’s small reversal is retained rather than hidden; adaptive simulations do not promise a monotone error sequence.

Resolution level	Remnant mass $M_f/M$	Absolute difference from Lev6
2	0.951970728685	$6.2211019\times10^{-5}$
3	0.952002973090	$2.9966614\times10^{-5}$
4	0.952000206769	$3.2732935\times10^{-5}$
5	0.952020101220	$1.2838484\times10^{-5}$
6	0.952032939704	$0$

The waveform check removes the first $500M$, as in the catalog quality study, interpolates the $(2,2)$ mode to 16,384 uniform samples, applies a declared $200M$ sine-squared taper at each end, and compares both real quadratures of this complex mode with a flat one-sided PSD. They are its real part and minus its imaginary part, not the full plus and cross strain along an observer’s line of sight; that strain requires a harmonic sum. The check uses every rFFT bin from DC through Nyquist, $0\le f\le[2\Delta u]^{-1}$, with trapezoidal endpoint weights $1/2$ at DC and Nyquist and unit weights on interior bins. The common $4\Delta f/S_n$ factor cancels in the normalized overlap. The check holds the catalog retarded-time and phase origins fixed; there is no time or phase maximization. Thus $1-\mathcal O_{22}$ below is a reproducible product-consistency residual, not the 20-direction, all-mode, optimized mismatch distribution in the paper.

Extrapolation order comparison	One-sided-flat-PSD residual $1-\mathcal O_{22}$
2–3	$5.2369960\times10^{-7}$
2–4	$4.1078488\times10^{-7}$
3–4	$5.4616699\times10^{-8}$

This scale is consistent with the SXS catalog’s finding that its $N=2$ versus $N=3$ distribution peaks below $10^{-6}$, but the comparison is conditional on this simulation, mode, taper, time interval, flat PSD, fixed alignment, and product version Boyle et al. 2019, §4 and Figure 11, pp. 24–30, PDF.

Executable rerun

The program emits the two tables above. It also verifies every checksum and performs a deliberately naive energy-balance check so that differentiation and extraction failure are visible.

import hashlib
import json
import re
from pathlib import Path

import h5py
import numpy as np

root = Path("data")
wave_path = root / "SXS_BBH_0305_Lev6_CoM.h5"
expected_wave_md5 = "7e7cc4b5a929d92ca49d2a60af19f4de"
assert hashlib.md5(wave_path.read_bytes()).hexdigest() == expected_wave_md5

expected_meta_md5 = {
    2: "9743f16509e3792daf3acb453edba87c",
    3: "e17321064d863adc8cb3d45cca1ec28f",
    4: "8ac4060359a22bb086a723e4853968e6",
    5: "ffd14ae47788556b0c69f776a4a89c0f",
    6: "88a25dd6158f2dc9a267ab16c9aa8d3e",
}
metadata = {}
for lev, expected in expected_meta_md5.items():
    path = root / f"Lev{lev}.json"
    assert hashlib.md5(path.read_bytes()).hexdigest() == expected
    metadata[lev] = json.loads(path.read_text())

reference_mass = metadata[6]["remnant_mass"]
print("resolution,remnant_mass,absolute_difference_from_Lev6")
for lev, record in metadata.items():
    mass = record["remnant_mass"]
    print(f"{lev},{mass:.12f},{abs(mass-reference_mass):.9e}")

with h5py.File(wave_path, "r") as waveform:
    history = waveform["VersionHist.ver"].asstr()[()]
    assert history.tolist() == [[
        "ef51849550a1d8a5bbdd810c7debf0dd839e86dd",
        "Overall sign change in complex strain h.",
    ]]

    modes = {}
    for order in (2, 3, 4):
        data = waveform[f"Extrapolated_N{order}.dir/Y_l2_m2.dat"][:]
        modes[order] = data[:, 1] + 1j * data[:, 2]
        if order == 2:
            time = data[:, 0]
        else:
            assert np.array_equal(time, data[:, 0]), "Extrapolation grids differ"

    mask = time >= time[0] + 500.0
    time_cut = time[mask]
    uniform_time = np.linspace(time_cut[0], time_cut[-1], 16384)
    taper = np.ones_like(uniform_time)
    width = 200.0
    left = uniform_time < uniform_time[0] + width
    right = uniform_time > uniform_time[-1] - width
    taper[left] = np.sin(np.pi*(uniform_time[left]-uniform_time[0])/(2*width))**2
    taper[right] = np.sin(np.pi*(uniform_time[-1]-uniform_time[right])/(2*width))**2

    spectra = {}
    for order, mode in modes.items():
        plus = np.interp(uniform_time, time_cut, mode[mask].real) * taper
        cross = -np.interp(uniform_time, time_cut, mode[mask].imag) * taper
        spectra[order] = (np.fft.rfft(plus), np.fft.rfft(cross))

    def overlap(order_a, order_b):
        ap, ax = spectra[order_a]
        bp, bx = spectra[order_b]
        weights = np.ones(ap.size)
        weights[[0, -1]] = 0.5
        product = np.vdot(weights*bp, ap).real + np.vdot(weights*bx, ax).real
        norm_a = np.vdot(weights*ap, ap).real + np.vdot(weights*ax, ax).real
        norm_b = np.vdot(weights*bp, bp).real + np.vdot(weights*bx, bx).real
        return product / np.sqrt(norm_a * norm_b)

    print("extrapolation_orders,one_sided_flat_PSD_residual")
    for first, second in ((2, 3), (2, 4), (3, 4)):
        print(f"{first}-{second},{1-overlap(first, second):.9e}")

    initial_adm = metadata[6]["initial_ADM_energy"]
    final_mass = metadata[6]["remnant_mass"]
    print("extrapolation_order,radiated_energy,balance_residual")
    for order in (2, 3, 4):
        group = waveform[f"Extrapolated_N{order}.dir"]
        radiated = 0.0
        for key in group:
            if not re.fullmatch(r"Y_l\d+_m-?\d+\.dat", key):
                continue
            data = group[key][:]
            mode_time = data[:, 0]
            mode = data[:, 1] + 1j * data[:, 2]
            news = np.gradient(mode, mode_time, edge_order=2)
            radiated += np.trapezoid(abs(news)**2, mode_time) / (16*np.pi)
        residual = initial_adm - final_mass - radiated
        print(f"{order},{radiated:.12f},{residual:.9e}")

Save the Python program above as sxs_product_check.py. The exact retrieval and rerun recipe is:

python3 -m venv .venv
. .venv/bin/activate
python3 -m pip install numpy==2.0.2 h5py==3.11.0
mkdir -p data
for lev in 2 3 4 5 6; do curl -L "https://zenodo.org/api/records/3301877/files/SXS:BBH:0305/Lev${lev}/metadata.json/content" -o "data/Lev${lev}.json"; done
curl -L "https://zenodo.org/api/records/3301877/files/SXS:BBH:0305/Lev6/rhOverM_Asymptotic_GeometricUnits_CoM.h5/content" -o data/SXS_BBH_0305_Lev6_CoM.h5
python3 sxs_product_check.py

Balance law as a failure diagnostic

Restoring $G$, the positive outward energy flux represented by the strain modes is

$$\mathcal F_E^\infty =\frac{1}{16\pi G} \sum_{\ell,m}\left\lvert\dot H_{\ell m}\right\rvert^2.$$

In the file’s $G=c=M=1$ normalization, the program differentiates every available mode and compares $E_{\mathrm{ADM}}-M_f-E_{\mathrm{rad}}$. It produces:

Extrapolation order	Radiated energy $E_{\mathrm{rad}}/M$	Balance residual $E_{\mathrm{ADM}}-M_f-E_{\mathrm{rad}}$
2	0.042072358310	$-1.7511125\times10^{-3}$
3	0.043201548367	$-2.8803026\times10^{-3}$
4	0.044317661189	$-3.9964154\times10^{-3}$

The worsening residual is a failed convergence test, not evidence that higher extrapolation order radiates more physical energy. Direct finite differences amplify early-time and high-order extrapolation noise, and the full time series includes junk. A reliable balance computation should use a controlled news product or documented filtering, vary junk cuts and differentiation, include mode truncation, compare extrapolation order or CCE worldtubes, and also test angular and linear momentum against remnant spin and recoil. Computational Relativity develops methods for calculating those quantities. The comparison here shows why a balance check that fails to close cannot validate the waveform.

The displayed subtraction is deliberately incomplete as a physical conservation test. In one common asymptotic inertial frame, the final energy is $E_f=\sqrt{M_f^2+\lvert\boldsymbol P_f\rvert^2}$, not the remnant rest mass alone when recoil is nonzero. Initial energy, radiative flux, and remnant four-momentum must be expressed in that same frame. A finite waveform also omits any radiation outside its time interval; after removing junk one must use the corresponding early Bondi energy or account for the removed flux. Thus the residual does not, by itself, isolate a discretization error or establish the continuum radiated energy.

Failure modes and error budget

The main separable contributions are evolution truncation, finite-radius or characteristic evolution error, extrapolation order, angular-mode truncation, interpolation and differentiation, junk excision, finite waveform duration, BMS-frame ambiguity, and horizon/remnant measurement error. Resolution differences estimate a combination of these, not solely the highest-resolution truncation error. The SXS study explicitly notes nonconvergent adaptive cases and that its highest-versus-next-highest mismatch measures the lower resolution’s error Boyle et al. 2019, §4, pp. 24–29, PDF.

Checksum failure or incompatible strain versions must be resolved before reuse. Missing metadata, unexplained resolution reversals, extraction sensitivity beyond the claimed tolerance, or a failed balance diagnostic prevent the corresponding accuracy claim until investigated. A deliberately inadequate energy estimator does not by itself invalidate every use of the strain product. No mismatch threshold is universal: any detector match must state its one-sided PSD, band, total-mass scaling, modes, frame, taper, and maximized shifts.

Exercise

For the $N=2$ row, compute the balance residual both as a percentage of the initial ADM energy $0.992354185511M$ and as a percentage of the computed radiated energy. Explain why quoting only the smaller percentage would be misleading.

Solution

The two relative residuals are

$$100\frac{1.7511125\times10^{-3}}{0.992354185511} =0.1765\%,$$ $$100\frac{1.7511125\times10^{-3}}{0.042072358310} =4.162\%.$$

The first denominator includes the nearly unchanged total mass and makes the closure failure look small. The second ratio exposes a discrepancy equal to roughly four percent of the reconstructed radiated energy; it is not a measured four-percent continuum error, because the missing interval, frame, and recoil corrections have not been isolated. Together with the worsening $N=2,3,4$ trend, it rules out this naive procedure as a precision balance calculation even though the $(2,2)$ product-consistency residual is below $10^{-6}$.

References

• Boyle, Michael, et al. “The SXS Collaboration Catalog of Binary Black Hole Simulations.” Classical and Quantum Gravity 36 (2019): 195006. DOI. Open PDF.
• SXS Collaboration. “Binary Black-Hole Simulation SXS:BBH:0305.” Zenodo, record 3301877, versioned data deposit. Accessed August 18, 2026. DOI.