Calibration and Extrapolation Domains

Scrollable formula. Use Left and Right arrow keys to reveal hidden content.

Scrollable table. Use Left and Right arrow keys to reveal additional columns.

Scrollable figure. Use Left and Right arrow keys to inspect hidden content.

A calibrated waveform model has at least three different domains: the mathematical domain on which its equations evaluate, the support of the waveforms used to fit its adjustable terms, and the validation domain on which errors were independently tested. Only the first is an implementation fact. The other two determine what an accuracy statement can mean. This page records those distinctions for a dated SEOBNRv5HM calibration set and develops reproducible coverage diagnostics. It does not recalibrate the model, reproduce NR evolutions, or turn validation results into inference claims.

Helpful background. Einstein Equations supplies the underlying theory; Waves, Green Functions, and Geometric Optics supplies waveform comparison norms; Relativistic Two-Body Dynamics supplies the analytic EOB structure; Computational Relativity supplies numerical calibration products and their uncertainties; and Gravitational-Wave Astronomy treats detector-specific inference verdicts.

What a calibration record must contain

A usable record names the model and source version, fitted coefficients, calibration waveform identifiers, intrinsic coordinates and reference epochs, modes and frames, fit objective and weights, waveform alignment, train/validation split, numerical uncertainties, and domain boundaries. It must also distinguish interpolation from extrapolation and state what happens when a request lies beyond tested support.

For aligned-spin SEOBNRv5HM, the intrinsic calibration coordinates are $\Lambda=(q,\chi_1,\chi_2)$ . Here $q=m_1/m_2\ge1$ and body 1 is the heavier object, while

\chi_i=\frac{\boldsymbol S_i\cdot\widehat{\boldsymbol L}}{Gm_i^2}

is positive for spin aligned with the orbital angular momentum and negative for antialignment. Each value must be attached to the reference epoch stated by its source record; for SXS entries this is the relaxed post-junk reference data rather than an unspecified initial slice. Pompili and collaborators describe 442 NR waveforms: 441 SXS evolutions and one Einstein Toolkit evolution, supplemented by 13 Teukolsky waveforms at $q=10^3$ for test-mass information Pompili et al. 2023, §IV and Figure 2, pp. 11–12, PDF. The SXS records retain simulation identifiers, resolutions, extraction products, and version metadata needed to interpret those coordinates Boyle et al. 2019, §§2–4, pp. 7–30, PDF. The NR set spans $1\le q\le20$ , but its support is highly nonuniform: for $10\lt q\le20$ the available spins are only on the heavier hole and reach about $\lvert\chi_1\rvert\le0.5$ , or are nonspinning. A rectangular declaration such as “ $q\le20$ , $\lvert\chi_i\rvert\le1$ ” would therefore invent coverage in sparse corners.

The dated paper’s dominant-mode objective uses

(h_1|h_2)=4\operatorname{Re} \int_{f_l}^{f_h} \frac{\widetilde h_1(f)\widetilde h_2(f)^*}{S_n(f)}df

with a one-sided Advanced-LIGO design PSD, $f_h=2048\,\mathrm{Hz}$ , a start-frequency prescription, Planck taper, and maximization over relative time and phase. Calibration seeks maximum $(2,2)$ -mode unfaithfulness below $10^{-3}$ over detector-frame (redshifted) total masses $10M_\odot$ – $200M_\odot$ and merger-time error below $5GM$ Pompili et al. 2023, §IV.A–B, Eqs. (65)–(70), pp. 12–13, PDF. Those settings are part of the loss; changing PSD, mass range, modes, or alignment changes the objective.

Calibration support is not a box

Convex-hull membership is a first geometric warning: a point outside the hull of the calibration locations is certainly extrapolative under those coordinates. But the converse fails. A hull joins its extreme points with flat facets and can bridge a large empty region. A nearest-neighbor distance detects local sparsity that hull membership misses, while neither quantity accounts for anisotropic waveform variation.

This page normalizes locations as

X=\left( \frac{\ln q}{\ln20}, \frac{\chi_1+1}{2}, \frac{\chi_2+1}{2} \right).

Distances are Euclidean only in this declared diagnostic coordinate. They are not waveform mismatches or physical priors. An alternative transform, spin combination, or local metric must be revalidated rather than silently substituted.

The source is arXiv:2303.18039v3 together with its published v1 ancillary file NR_simulations.json, accessed 2026-08-18, MD5 23038ccb1a01d3b2d4d78dd5c1685b66. The 442 actual locations divide as follows:

Mass-ratio region	Calibration-waveform count
$1\le q\le3$	278
$3\lt q\le10$	143
$10\lt q\le20$	21

Only $21/442=4.75\%$ lie above $q=10$ . A model may still extrapolate well there because its analytic backbone constrains behavior, but density alone provides little empirical protection.

Held-out domain test with controlled response

The ancillary file provides locations and physical metadata, not the calibration-posterior coefficients needed to reproduce the SEOBNRv5 fit. The benchmark therefore keeps the real 442-point geometry but assigns a clearly synthetic smooth scalar response plus heteroscedastic deterministic noise. Its purpose is to test splitting, weights, fit order, hull membership, and sparse-boundary reporting—not to estimate SEOBNRv5 error.

With random seed 1729, 88 locations are held out and the remaining 354 form a nested training pool. A quadratic polynomial is fit by inverse-variance weighted ridge regression. The validation target is the known noiseless response. For each subset, the convex hull is rebuilt; a held-out point is labeled boundary-risk if it is outside that hull or farther than $0.16$ from the nearest training point in $X$ . The baseline has 82 interior and 6 boundary-risk validation points.

The exact code was run with Python 3.9.6, NumPy==2.0.2, and SciPy==1.13.1:

Training grid count	Validation RMSE error	Interior error	Boundary error
88	2.633372e-3	1.477441e-3	3.484109e-3
176	2.622911e-3	2.549494e-3	2.954275e-3
354	2.584828e-3	2.582578e-3	2.615382e-3

The overall error decreases only modestly, and the composition of the interior/boundary subsets changes as the hull grows. These rows show training-subset sensitivity for this split-specific scalar fit; three decreasing values do not establish convergence to the true response or a universal learning curve. In particular, the lower interior error of the smallest subset is not evidence that it is the best model: its interior subset contains different validation points.

Fit and boundary sensitivities

Keeping the same 88 held-out locations gives:

Variation	Validation error	Interior error	Boundary error
Quadratic weighted baseline	2.584828e-3	2.582578e-3	2.615382e-3
Half training subset	2.622911e-3	2.549494e-3	2.954275e-3
Unweighted quadratic	2.506417e-3	2.468827e-3	2.972897e-3
Linear fit	2.827491e-2	2.685670e-2	4.322087e-2
Regularized cubic	6.669098e-4	6.556735e-4	8.049027e-4

Unweighted fitting slightly lowers the global error while worsening the boundary error; a single aggregate score would hide that tradeoff. The linear fit fails by more than an order of magnitude. The cubic succeeds on this controlled response, but a higher-order fit is not automatically preferable for noisy physical calibration data—it must survive independent splits and regularization changes.

Four query diagnostics further separate “inside” from “supported”:

Query $(q,\chi_1,\chi_2)$	Inside hull	Nearest normalized distance	Minimum hull margin
Central $(2,0,0)$	yes	0.000071	0.231428
Exact ET:AEI:0004 metadata	yes	0	approximately 0
Sparse high- $q$ spin $(15,0.8,0.8)$	no	0.202118	-0.170974
Outside mass ratio $(25,0,0)$	no	0.074549	-0.066840

The margin is the minimum signed hull half-space margin: positive inside and negative outside. ET:AEI:0004 is itself an extreme calibration point, so zero nearest distance does not provide an interior buffer. Conversely, $q=25$ is closer to a boundary point than the sparse $q=15$ spin query, even though both are outside. Hull flag and local distance answer different questions.

For a waveform application, add physically targeted holdouts: a high-spin slice, a high-mass-ratio slice, or an entire simulation family. A random split can leave nearly identical physical configurations on both sides even when identifiers differ. Such group holdouts probe transfer across the intended parameter regime; the random split below is adequate for its declared synthetic-regression example but does not establish that stronger physical generalization.

Executable rerun

import hashlib
import json
from pathlib import Path

import numpy as np
from scipy.spatial import ConvexHull, Delaunay, distance_matrix

SEED = 1729
path = Path("NR_simulations.json")
expected_md5 = "23038ccb1a01d3b2d4d78dd5c1685b66"
assert hashlib.md5(path.read_bytes()).hexdigest() == expected_md5
records = json.loads(path.read_text())

q = np.array([item["reference_mass_ratio"] for item in records])
chi1 = np.array([item["reference_dimensionless_spin1"] for item in records])
chi2 = np.array([item["reference_dimensionless_spin2"] for item in records])
coordinates = np.column_stack((np.log(q)/np.log(20), (chi1+1)/2, (chi2+1)/2))
print("q_region,count")
for label, use in (
    ("1_to_3", q <= 3),
    ("3_to_10", (q > 3) & (q <= 10)),
    ("10_to_20", q > 10),
):
    print(f"{label},{np.count_nonzero(use)}")

rng = np.random.default_rng(SEED)
permutation = rng.permutation(len(records))
validation_indices = permutation[:88]
pool_indices = permutation[88:]

def response(indices, with_noise=True):
    x = coordinates[indices, 0]
    c1, c2 = chi1[indices], chi2[indices]
    exact = 0.15*x + 0.20*(q[indices]*c1+c2)/(q[indices]+1)
    exact += -0.08*((c1-c2)/2)**2 + 0.05*x*c1 + 0.02*np.sin(np.pi*x)*c2
    sigma = 0.001*(1 + 2*x + abs(c1) + abs(c2))
    if not with_noise:
        return exact, sigma
    noise = np.random.default_rng(SEED+1).normal(size=len(records))
    return exact + sigma*noise[indices], sigma

def design(points, degree):
    x, y, z = points.T
    columns = [np.ones(len(points)), x, y, z]
    if degree >= 2:
        columns += [x*x, y*y, z*z, x*y, x*z, y*z]
    if degree >= 3:
        columns += [x**3, y**3, z**3, x*x*y, x*x*z, y*y*x,
                    y*y*z, z*z*x, z*z*y, x*y*z]
    return np.column_stack(columns)

def fit_and_score(train_indices, degree=2, ridge=1e-8, weighted=True):
    matrix = design(coordinates[train_indices], degree)
    values, sigma = response(train_indices)
    weights = 1/sigma**2 if weighted else np.ones_like(sigma)
    normal = matrix.T @ (weights[:, None]*matrix) + ridge*np.eye(matrix.shape[1])
    coefficients = np.linalg.solve(normal, matrix.T @ (weights*values))
    validation_matrix = design(coordinates[validation_indices], degree)
    exact, _ = response(validation_indices, with_noise=False)
    residual = validation_matrix @ coefficients - exact

    hull = ConvexHull(coordinates[train_indices])
    triangulation = Delaunay(coordinates[train_indices][hull.vertices])
    inside = triangulation.find_simplex(coordinates[validation_indices]) >= 0
    nearest = distance_matrix(coordinates[validation_indices],
                              coordinates[train_indices]).min(axis=1)
    boundary = (~inside) | (nearest > 0.16)
    interior_rmse = np.sqrt(np.mean(residual[~boundary]**2))
    boundary_rmse = np.sqrt(np.mean(residual[boundary]**2))
    return np.sqrt(np.mean(residual**2)), interior_rmse, boundary_rmse, coefficients

print("training_count,validation_rmse_error,interior_error,boundary_error")
for count in (88, 176, 354):
    score = fit_and_score(pool_indices[:count])
    print(f"{count},{score[0]:.6e},{score[1]:.6e},{score[2]:.6e}")

print("variation,validation_error,interior_error,boundary_error")
for name, count, degree, ridge, weighted in (
    ("baseline", 354, 2, 1e-8, True),
    ("half_subset", 176, 2, 1e-8, True),
    ("unweighted", 354, 2, 1e-8, False),
    ("linear_fit", 354, 1, 1e-8, True),
    ("cubic_regularized", 354, 3, 1e-4, True),
):
    score = fit_and_score(pool_indices[:count], degree, ridge, weighted)
    print(f"{name},{score[0]:.6e},{score[1]:.6e},{score[2]:.6e}")

hull = ConvexHull(coordinates)
triangulation = Delaunay(coordinates[hull.vertices])
queries = {
    "central": (2.0, 0.0, 0.0),
    "ET_AEI_0004": (7.992026232813275, 0.851333746356202, 0.8503915330008012),
    "sparse_high_q_spin": (15.0, 0.8, 0.8),
    "outside_mass_ratio": (25.0, 0.0, 0.0),
}
print("query,inside_hull,nearest_distance,hull_margin")
for name, (mass_ratio, spin1, spin2) in queries.items():
    point = np.array([np.log(mass_ratio)/np.log(20), (spin1+1)/2, (spin2+1)/2])
    inside = triangulation.find_simplex(point) >= 0
    nearest = np.linalg.norm(coordinates-point, axis=1).min()
    margin = np.min(-(hull.equations[:, :3]@point + hull.equations[:, 3]))
    print(f"{name},{inside},{nearest:.6f},{margin:.6f}")

Save the Python program above as calibration_domain_check.py. The program emits the coverage, convergence, fit-sensitivity, and query tables. The exact rerun recipe is:

python3 -m venv .venv
. .venv/bin/activate
python3 -m pip install numpy==2.0.2 scipy==1.13.1
curl -L "https://arxiv.org/src/2303.18039v1/anc/NR_simulations.json" -o NR_simulations.json
python3 calibration_domain_check.py

Failure modes and interpretation limits

Calibration evidence fails if identifiers or source versions are missing, private inputs cannot be distinguished from public ones, reference spins are mixed across epochs, NR resolution/extraction uncertainty is omitted, fitted and held-out cases leak across near-duplicate configurations, or the loss is reported without its PSD, band, modes, and alignment. A fit also fails review when subset, weight, order, regularization, or boundary results change materially without entering the error budget.

Pompili et al. estimate NR uncertainty by comparing highest and next-highest resolutions and extrapolation orders, and show that worst cases concentrate near sparse parameter-space corners Pompili et al. 2023, §V.A–C, Figures 5–12 and Tables I–II, PDF. Their published mismatch fractions are conditional on the 2023 model version, listed modes, NR products, total-mass intervals, PSD, inclination averaging, and alignment. They validate model behavior—including comparisons to NRHybSur3dq8—but do not turn a convex hull into an accuracy theorem.

Exercise

Compute the fraction of calibration waveforms with $q\gt 10$ . Then compare the linear-fit boundary error with its interior error and with the baseline boundary error. What distinct failures do the ratios diagnose?

Solution

The high-mass-ratio fraction is

\frac{21}{442}=0.04751=4.751\%.

For the linear fit, the boundary-to-interior ratio is

\frac{4.322087\times10^{-2}}{2.685670\times10^{-2}} =1.609.

The linear boundary error relative to the quadratic baseline boundary error is

\frac{4.322087\times10^{-2}}{2.615382\times10^{-3}} =16.53.

The first ratio isolates added boundary risk within one underfit model. The second shows a much larger global model-complexity failure: a linear response is inadequate even before extrapolating. Sparse support and fit bias are separate error sources and should not be represented by one domain flag.

References

Boyle, Michael, et al. “The SXS Collaboration Catalog of Binary Black Hole Simulations.” Classical and Quantum Gravity 36 (2019): 195006. DOI. Open PDF.
Pompili, Lorenzo, et al. “Laying the Foundation of the Effective-One-Body Waveform Models SEOBNRv5: Improved Accuracy and Efficiency for Spinning Nonprecessing Binary Black Holes.” Physical Review D 108 (2023): 124035. DOI. Open PDF.