Numerical Waveform Error Budgets

A numerical waveform error budget must distinguish evolution discretization, finite-radius extraction, extrapolation or characteristic transport, initial-data transients, frame correction, interpolation, and comparison preprocessing. A single highest-minus-next-resolution mismatch mixes these effects and can even become nonmonotonic. This page gives a reproducible convergence control and a defensible combination workflow. Current detector curves and event-specific parameter biases are outside its scope.

Required background. Probability and Statistical Inference supplies covariance, correlation, and uncertainty propagation.

Helpful background. Einstein Equations supplies the continuum target. Waves, Green Functions, and Geometric Optics supplies extraction and spectral comparison. Gravitational Radiation and Waveforms supplies modes, frames, and mismatch. Computational Relativity supplies discretization, constraints, and numerical validation.

Decompose before combining

Write the numerical waveform schematically as

$$h_{N,r,k}(u)=h_\infty(u) +e_{\rm disc}(N) +e_{\rm radius}(r) +e_{\rm extrap}(k) +e_{\rm junk} +e_{\rm frame} +e_{\rm interp}.$$

$N$ denotes an evolution resolution, $r$ an extraction radius, and $k$ an extrapolation order. Here each $e_i$ is locally defined as a numerical-minus-continuum contribution; it is distinct from the chapter’s inference residual $\delta h=h_*-h_m$. This is bookkeeping, not an assumption of independence. For example, finite resolution changes coordinate trajectories used in extrapolation, while center-of-mass drift can differ across resolutions. Measure signed waveform differences before reducing them to scalar norms.

Declare the simulation identifier and version, resolution labels, extraction product, extrapolation orders, retarded time, BMS/center-of-mass correction, modes, junk cut, interpolation grid, alignment window, reference waveform, PSD or time norm, and frequency band. Without those inputs an error number cannot be reproduced.

Three-resolution convergence

For a scalar phase diagnostic with grid scale $h_g$ and leading error

$$\phi(h_g)=\phi_0+Ch_g^p+O(h_g^{p+1}),$$

three resolutions with constant refinement ratio $r=h_{g,c}/h_{g,m}=h_{g,m}/h_{g,f}$ give

$$p_{\rm obs} =\frac{\ln\left\lvert (\phi_c-\phi_m)/(\phi_m-\phi_f) \right\rvert}{\ln r}.$$

Only intervals with consistent signs and stable $p_{\rm obs}$ support Richardson extrapolation. Apply the check separately to unwrapped phase, relative amplitude, constraints, and waveform norms. A mismatch can converge while phase does not because time/phase alignment projects away leading error.

Reproducible fourth-order control

The synthetic control uses an exact phase $\phi_0=1.2$ rad and injects $0.8(64/N)^4$. It tests the analysis program, not a numerical-relativity code. Use CPython 3.12.7 and NumPy 2.1.3. Save the following program as waveform_error_budget.py. The script writes waveform_error_budget.csv, emits the rows, and estimates the observed order.

import numpy as np

resolution = np.array([64.0, 96.0, 144.0, 216.0])
phi_exact = 1.2
p_true = 4.0
phase = phi_exact + 0.8 * (64.0 / resolution) ** p_true
error = np.abs(phase - phi_exact)
ratio = resolution[1:] / resolution[:-1]
orders = []
for index in range(len(phase) - 2):
    numerator = abs(phase[index] - phase[index + 1])
    denominator = abs(phase[index + 1] - phase[index + 2])
    orders.append(np.log(numerator / denominator) / np.log(ratio[index]))
rows = np.column_stack((resolution, phase, error))
np.savetxt("waveform_error_budget.csv", rows, delimiter=",",
           header="resolution_n,phase_result,error", comments="")
for row in rows:
    print(f"{int(row[0]):4d} {row[1]:.9f} {row[2]:.9f}")
print("observed orders:", *(f"{value:.6f}" for value in orders))

Create an isolated environment and rerun exactly:

python3.12 -m venv .venv-wave-error
.venv-wave-error/bin/python -m pip install numpy==2.1.3
.venv-wave-error/bin/python waveform_error_budget.py

The expected refinement table is:

resolution $N$	phase result (rad)	absolute error (rad)
64	2.000000000	0.800000000
96	1.358024691	0.158024691
144	1.231214754	0.031214754
216	1.206165877	0.006165877

Both triplets recover $p_{\rm obs}=4$. Changing the alignment independently at each resolution would contaminate this known benchmark and is therefore a useful failure test.

A published catalog-scale benchmark

The SXS catalog records exact simulation identifiers SXS:BBH:dddd, resolution labels LevN, multiple extrapolation orders, strain and $\Psi_4$ products, and center-of-mass corrections. Its catalog-wide comparison found highest-versus-second-resolution flat-weight mismatches broadly around $10^{-6}$–$10^{-3}$, with typical scale near $10^{-4}$, while a nontrivial subset was nonconvergent under the adjacent-resolution test. Comparing extrapolation orders $N=2$ and $N=3$ peaked below $10^{-6}$ but had a tail reaching $10^{-3}$. These are population-level scales, not an error assigned to every entry; see Boyle et al. 2019, §4, pp. 24–28 and Figs. 9–11, PDF.

The same study emphasizes that highest-minus-second differences estimate the second-highest error, mix truncation with extraction, and can be distorted when resolutions acquire slightly different physical parameters. A real rerun must name the chosen SXS:BBH:dddd, every LevN, extrapolation order, mode set, and preprocessing; this page does not copy or redistribute a catalog waveform.

Extraction and frame stress tests

For each selected waveform, repeat the comparison under the following variations:

• use at least three evolution resolutions and report nonmonotonic intervals;
• compare extrapolation orders $k=2,3,4$ and, when available, characteristic extraction;
• shift junk-removal and alignment windows rather than choosing the minimum error;
• compare dominant-mode and multimode norms in the same inertial/BMS frame;
• vary interpolation sampling and taper;
• repeat with different reference resolutions and both time- and frequency-domain diagnostics.

Finite-radius error should decrease with radius only in a controlled common gauge. High extrapolation order can amplify noise, so increasing $k$ is not automatically refinement. $\Psi_4$ and strain comparisons must control two integrations, low-frequency content, and BMS frame.

The distinction among coordinate-radius extraction, extrapolation, and Cauchy-characteristic extraction—including their separate assumptions and comparisons—is reviewed in Bishop and Rezzolla 2016, §3.3, pp. 22–28; §§6–7, pp. 51–76; and §8, pp. 76–79, PDF.

Combining the budget

If the total error waveform is $e=\sum_i e_i$, bilinearity gives the deterministic identity

$$\lVert e\rVert^2 =\sum_{ij}(e_i\mid e_j).$$

No probabilistic covariance assumption is needed for this identity. Cross terms measure whether the signed errors reinforce or cancel in the declared waveform norm. For a probabilistic budget, take an ensemble expectation and specify the means and cross-covariances separately.

If cross terms are unknown, the triangle bound

$$\lVert e\rVert \le\sum_i\lVert e_i\rVert$$

is conservative. Adding deterministic norms in quadrature is exact for orthogonal errors. Independent, zero-mean random errors justify quadrature of expected squared norms, not exact cancellation of cross terms in each realization. Coherent phase errors can exceed a quadrature estimate. Avoid counting extraction disagreement once as “finite radius” and again inside the resolution difference that already contains it.

Failure modes

Phase wrapping can create false jumps; align and unwrap with a recorded rule. Optimizing time and phase separately for every pair can cancel the leading discretization error. Spectral adaptive mesh labels are not universal grid spacings, so a formal order may not exist. Junk radiation changes initial masses and contaminates early alignment. Frame drift mixes modes. Near merger, refinement can be nonmonotonic. Interpolation and resampling can dominate a very accurate evolution. Insufficient resolution can make all three levels agree on the same biased branch.

These are reasons to report failed convergence, not to delete inconvenient intervals. Production algorithms and continuum validation remain part of Computational Relativity.

Exercises

1. Recover the order and continuum phase. Use the final three table rows with refinement ratio $r=1.5$ to estimate $p$ and Richardson-extrapolate from $N=144,216$.

Solution

The difference ratio is

$$\frac{1.358024691-1.231214754} {1.231214754-1.206165877} =5.0625=1.5^4,$$

so $p=4$. The fine-grid continuum estimate is

$$\phi_0\simeq\phi_f+ \frac{\phi_f-\phi_m}{r^p-1} =1.206165877+ \frac{1.206165877-1.231214754}{1.5^4-1} =1.200000000.$$

It recovers the injected exact value because the control contains a single pure fourth-order term.

2. Correlated components. Two phase-error waveforms are identical and each has norm $0.01$. Compare quadrature with the actual combined norm.

Solution

Quadrature would give $\sqrt{0.01^2+0.01^2}=0.0141$. Because the signed errors are identical, they add coherently and the true norm is $0.02$. Their cross term is positive and equal to each self term. The example shows why correlated extraction and calibration effects cannot be combined as independent scalars without evidence.

References

• Bishop, Nigel T., and Luciano Rezzolla. “Extraction of Gravitational Waves in Numerical Relativity.” Living Reviews in Relativity 19 (2016): 2. DOI. Open PDF.
• Boyle, Michael, et al. “The SXS Collaboration Catalog of Binary Black Hole Simulations.” Classical and Quantum Gravity 36 (2019): 195006. DOI. Open PDF.