Analytic–Numerical Comparisons for Strong Scattering

Testing an analytic strong-scattering prediction against numerical relativity requires reconstructing common asymptotic invariants and separating numerical resolution, extraction, alignment, and analytic-truncation errors. A finite-separation coordinate deflection becomes an asymptotic scattering estimate only after a frame and extrapolation procedure are specified. The examples below progress from exact test-mass motion and a published finite-mass comparison to an independently reproducible extraction from two released numerical-relativity trajectories.

Helpful background. The Einstein Equation fixes the theory, Weak Metric Perturbations and Gauge Freedom supplies extraction variables, Presymplectic Potentials and Currents supplies ADM and radiative charges, Computational Relativity supplies convergence methods, Relativistic Two-Body Dynamics supplies analytic scattering predictions, and Black-Hole Geometry and Mechanics supplies the Schwarzschild benchmark.

Scattering data and reference calculations

The geodesic calculation can be run immediately with Python 3 and mpmath==1.3.0:

python3 -m pip install mpmath==1.3.0
python3 schwarzschild-scattering.py

Download the linked script before running the second command. It prints its environment and checks an exact test-mass sequence by quadrature, elliptic integrals, and independent orbit integration. The published-table comparison uses only the Python standard library and recomputes the arithmetic in the worked example below. Neither program evolves Einstein’s equations.

The SXS extraction calculation goes further: it downloads original horizon trajectories at three numerical resolutions, verifies their checksums, varies the extrapolation model and fitting windows, and compares the result with successive PM orders. These are separate checks. A stable trajectory fit does not establish the accuracy of a truncated analytic model, and a small resolution difference does not remove finite-separation extrapolation bias.

For any comparison, preserve the source’s definitions of masses, spins, energy and angular momentum, including its treatment of initial junk radiation. Initial ADM charges and charges measured after the junk pulse are different inputs. Radiative losses or a covariance matrix should be propagated when supplied; they cannot be reconstructed from a scattering angle alone.

Reconstructing the scattering angle

At large separation, fit incoming and outgoing motion separately in a center-of-mass asymptotic frame. If $\phi_{\mathrm{in}}^\infty$ and $\phi_{\mathrm{out}}^\infty$ are the unwrapped polar angles of the relative separation vector on the incoming and outgoing branches,

$$\chi =\phi_{\mathrm{out}}^\infty -\phi_{\mathrm{in}}^\infty-\pi,$$

with orientation and branch fixed continuously across the sequence. The incoming relative velocity points opposite the incoming separation vector; this explains the subtraction of $\pi$. If velocities or momenta are used instead, their directed deflection already excludes that straight-path angle. Their endpoint directions alone determine the angle only modulo $2\pi$; retain the orbital winding when following a zoom–whirl sequence. Do not fit through the strong-field encounter.

At each resolution $h$ and extraction radius $R_{\mathrm{ext}}$, construct

$$\chi(h,R_{\mathrm{ext}}) =\chi_\infty +a h^p +\sum_{k=1}^{K} \frac{b_k}{R_{\mathrm{ext}}^k} +\cdots.$$

Refine $h$ and $R_{\mathrm{ext}}$ independently. Vary convergence order $p$, polynomial order $K$, and fit windows; treat the spread as model-selection uncertainty. Propagate the covariance of incoming $(E,J)$ because near threshold $\partial\chi/\partial J$ can be very large.

Radiated four-momentum provides an independent balance check:

$$P_{\mathrm{out}}^\mu =P_{\mathrm{in}}^\mu-P_{\mathrm{rad}}^\mu.$$

Here $P^\mu$ is the sum of the bodies’ momenta. Its spatial direction measures the recoil of their center of mass, not their relative deflection. For example, equal-mass nonspinning scattering has zero total spatial momentum by symmetry while its scattering angle can be large. Reconstructing the angle from momenta therefore requires the individual outgoing momenta, boosted to the declared comparison frame, as well as an unwrapped orbital branch. Total balance alone cannot determine it.

Exact and finite-ratio benchmarks

For Schwarzschild geodesics, evaluate

$$\chi_{\mathrm{geo}}(E,J) =2\int_{r_{\min}}^\infty \frac{d\phi}{dr}\,dr-\pi$$

with high-precision mpmath quadrature. Locate the turning point independently and compare with the analytic weak-deflection limit and logarithmic near-critical slope. Tighten precision until both integral and root are stable beyond the target tolerance.

For the finite-ratio benchmark, use at least one published numerical-relativity sequence with multiple resolutions and extraction radii. Preserve the source’s definition of incoming invariants and quote its DOI and data version. Parameters used to calibrate the analytic model are excluded from validation metrics.

The analytic PM/EOB comparison and invariant framing are discussed in Damour 2016, §§III–VIII. The expected near-critical behavior and its finite-mass numerical setting are developed in Pretorius and Khurana 2007, §§2–4, pp. S87–S101.

Worked comparison with published numerical results

An instructive historical test uses equal-mass, nonspinning black holes. The incoming energies are near $E_{\mathrm{in}}/M=1.0226$, with $G=c=1$ in the source. Its parameter $b_{\mathrm{NR}}$ labels the initial-data family; use the reported $E_{\mathrm{in}},J_{\mathrm{in}}$ for invariant comparisons rather than identifying that label with an exact asymptotic impact parameter. Selected entries from Damour et al. 2014, Tables I–II, pp. 2–3, Open PDF are:

$b_{\mathrm{NR}}/M$	Numerical angle	2014 EOB angle	3PN angle
10	$222.9\pm1.7$	$227\pm5$	$126\pm1$
12	$120.7\pm1.5$	$120.77\pm0.06$	$104.34\pm0.04$
16	$70.7\pm1.9$	$70.666\pm0.002$	$69.122\pm0.002$

Angles are in degrees; the uncertainties are the authors’ estimates. “2014 EOB” denotes their $\mathrm{EOBNR}_{5\mathrm{PNlog}}$ model, not a complete 5PN expansion or a current EOB implementation. At $b_{\mathrm{NR}}/M=10$, the EOB and PN central differences from the numerical result are about $1.84\%$ and $-43.47\%$. The contrast is a result for this sequence, not evidence that every resummation improves every encounter.

The paper compares conservative predictions at averaged charges,

$$\bar E=E_{\mathrm{in}}-\frac{\Delta E}{2}, \qquad \bar J=J_{\mathrm{in}}-\frac{\Delta J}{2}, \qquad \chi_{\mathrm{ret}}\simeq\chi_{\mathrm{cons}}(\bar E,\bar J).$$

This relation retains the first-order effect of radiation reaction in the stated perturbative treatment; it is not an exact rule for large losses. The first selected row gives $\bar E/M=1.0172121$ and $\bar J/M^2=1.087323$. The comparison script reproduces these conversions and the angle differences directly from the versioned source-table entries.

The source estimated numerical error from two resolutions with an assumed fourth-order convergence law plus extraction variation. That is weaker than independently measuring convergence with three or more resolutions. These small programs do not reconstruct the source’s continuum extrapolation, rerun its EOB model, or create a new blind test. They supply a concrete comparison and make its reproducibility boundary explicit.

The associated radiated charges, useful for checking the averaged-charge conversion, are retained here from the same source’s Table I. Parentheses denote the authors’ uncertainty in the final quoted digits.

$b_{\mathrm{NR}}/M$	$E_{\mathrm{in}}/M$	$J_{\mathrm{in}}/M^2$	$\Delta E/M$	$\Delta J/M^2$
10	1.0225791(50)	1.145523(38)	0.010734(75)	0.1164(14)
12	1.0225907(50)	1.374658(45)	0.001980(13)	0.04022(53)
16	1.0225932(50)	1.832883(58)	0.0003338(77)	0.0162(11)

Extracting angles from released SXS trajectories

The releases SXS:BBH:4292v3.0 and SXS:BBH:3999v3.0 contain horizon-center trajectories at three adaptive numerical resolutions, Lev1 through Lev3. Both systems are nominally equal mass and nonspinning. We use the supplied inertial coordinates in the nominal initial center-of-mass frame and measure each hole’s unwrapped deflection, then average:

$$\theta_i=\phi_{i,\mathrm{out}}^\infty -\phi_{i,\mathrm{in}}^\infty-\pi, \qquad \theta_{\mathrm{avg}}=\frac{\theta_1+\theta_2}{2}.$$

For unequal masses this average is generally distinct from the relative scattering angle defined using the final center-of-mass frame. The distinction is explained in Long et al. 2025, Eqs. (14), (26)–(27), PDF. Here small residual initial momenta and out-of-plane motion are recorded, not assumed to vanish exactly.

Let $D$ be the distance between horizon centers in units of the sum of their initial Christodoulou masses. Fit each distant leg with both a cubic polynomial in $1/D$ and a Keplerian relation $1/D=A\cos\phi+B\sin\phi+C$. We use incoming windows $[d_{\mathrm{in}},200]$, with $d_{\mathrm{in}}=50,55,\ldots,80$, and outgoing windows $[d_{\mathrm{out}},350]$, with $d_{\mathrm{out}}=80,90,\ldots,230$. Each model is used on both legs; the $2\times7\times16=224$ model/window pairs receive equal weight. The models and window ranges follow Long et al. 2025, §II.D, Eqs. (8)–(12), PDF; the grid increments and averaging rule are specified here for reproduction. The Keplerian asymptote is solved directly from its fitted equation, with the branch checked against a rotated exact hyperbola.

The resulting average angles, in degrees, are:

SXS release	Lev1	Lev2	Lev3
4292v3.0	179.793754	179.798525	179.798517
3999v3.0	88.823840	88.827834	88.827840

At Lev3, the paired model/window estimates range from $179.790707^\circ$ to $179.805446^\circ$ for 4292, and from $88.817092^\circ$ to $88.837097^\circ$ for 3999. These are extraction sensitivities, not rigorous bounds or Gaussian confidence intervals. The final resolution changes are only about $8.0\times10^{-6}$ and $6.2\times10^{-6}$ degrees. Those small changes do not establish the total accuracy: the extrapolation models still disagree at a much larger scale. The calculation also repeats the fits on a coarser window grid and at half and quarter sampling density. Lev labels do not supply uniform mesh spacings, so no Richardson convergence order is inferred.

Analytic order at the same initial parameters

For this numerical example set $G=c=1$. With $M=m_1+m_2$ and $\nu=m_1m_2/M^2$, define

$$\Gamma=\frac{E_{\mathrm{ADM}}}{M}, \qquad \gamma=1+\frac{\Gamma^2-1}{2\nu}, \qquad j=\frac{J_{\mathrm{ADM},z}-m_1^2\chi_{1z}-m_2^2\chi_{2z}} {m_1m_2}.$$

The tiny residual spins enter the conversion from total to orbital angular momentum; the analytic angle itself is nonspinning. At the actual Lev3 inputs, $\nu\simeq1/4$ and $\Gamma\simeq1.022642$ for both releases, while

Release	$\gamma$	$j$
4292v3.0	1.0915935213	4.8000782422
3999v3.0	1.0915936633	6.4001055573

The full angle through successive PM orders is

$$\begin{aligned} \theta_{1\mathrm{PM}} &=\frac{2(2\gamma^2-1)}{j\sqrt{\gamma^2-1}},\\ \theta_{2\mathrm{PM}} &=\theta_{1\mathrm{PM}}+ \frac{3\pi(5\gamma^2-1)}{4\Gamma j^2}. \end{aligned}$$

The second line includes both orders; its last term is the 2PM correction. These formulas are in radians; they double the half-angle coefficients in Damour 2020, Eqs. (1.4), (1.9)–(1.11), (2.45)–(2.46), PDF. Converting to degrees gives:

Release	1PM	Through 2PM	Extracted Lev3 mean
4292v3.0	75.440309	103.846283	179.798517
3999v3.0	56.580201	72.558559	88.827840

Increasing the analytic order moves both predictions toward the numerical angles, but a large discrepancy remains, especially for the more strongly deflected encounter. Numerical refinement cannot repair that analytic truncation. The difference between two PM orders is not a proven bound on the omitted terms. Radiation changes scattering from 3PM Damour 2020, §V, Eqs. (25)–(31), PDF, so the remaining discrepancy cannot be assigned wholly to conservative dynamics. No radiation flux at infinity is recomputed from these horizon files. The historical radiated-charge table above is a separate source, not a substitute for missing fluxes of these releases.

Repeat the extraction

Save the extraction script and input manifest in one folder, then run:

python3 -m venv .venv
.venv/bin/python -m pip install numpy==2.5.3 h5py==3.16.0
.venv/bin/python extract-scattering.py --fetch --cache sxs-cache --output my-results.json

Use Python 3.12; on Windows the environment’s interpreter is .venv\Scripts\python.exe. The first run downloads about 3.1 MB of original files and verifies their checksums; later runs can omit --fetch and use the cache offline. The recorded results retain every model/window angle and both refinement studies. The method and reproduction notes give exact release provenance, full-fit options and conventions. The release inputs do not coincide exactly with rows of the article’s ancillary table; this calculation neither substitutes a nearby row nor claims to reproduce its published central angle.

Error model and pass criteria

The following criteria apply to a claim of quantitative analytic–numerical agreement. The low-order PM examples above expose a discrepancy; they do not pass such an accuracy test.

For prediction $\chi_{\mathrm A}$ and numerical estimate $\chi_{\mathrm N}$,

$$\Delta\chi=\chi_{\mathrm A}-\chi_{\mathrm N},$$ $$\sigma_\Delta^2 =\sigma_{\mathrm{num}}^2 +\nabla\chi_{\mathrm A}^{\mathsf T} \Sigma_{EJ} \nabla\chi_{\mathrm A} +\sigma_{\mathrm{analytic}}^2,$$

where $\Sigma_{EJ}$ is the incoming-invariant covariance. This sum assumes that the displayed contributions are uncorrelated. In general, $\operatorname{Var}(\chi_{\mathrm A}-\chi_{\mathrm N})= \operatorname{Var}(\chi_{\mathrm A})+\operatorname{Var}(\chi_{\mathrm N}) -2\operatorname{Cov}(\chi_{\mathrm A},\chi_{\mathrm N})$; shared extraction data can make the covariance consequential. Analytic order variation is a diagnostic rather than automatically a Gaussian standard deviation. Report residuals and the full covariance across a sequence; points sharing an extrapolation are not independent.

Pass criteria require observed numerical convergence compatible with the discretization order; extraction-window stability; independent test-mass agreement; analytic-order variation consistent with the truncation estimate; withheld finite-ratio residuals covered at the declared confidence level when a probabilistic uncertainty model is justified; and balance of incoming, outgoing, and radiated four-momentum.

Failure modes

Common failures are mismatched $E$ or $J$, a finite-coordinate deflection mistaken for an asymptotic angle, mixing charges defined before and after the initial junk pulse, correlated fit errors treated as independent, near-threshold loss of significant digits, and validation on calibrated points. Warning signs include analytic residuals that shrink only after refitting every point, a critical value that moves outside uncertainty under window changes, and nonmonotonic resolution behavior without a conservative error enlargement.

Keep machine-readable values alongside a comparison plot: original case identifiers, aligned invariants, extrapolated angles, analytic predictions, and the uncertainty information actually available. The SXS example retains individual fit values instead of inventing a covariance matrix. This allows another reader to change the averaging rule and assess its effect.

Exercises

Why must extraction radius and numerical resolution be refined separately?

Solution

They represent different limits. Increasing grid resolution at fixed finite radius does not remove extraction bias, and moving extraction outward on an underresolved grid does not establish continuum convergence. A combined fit without independent variation can trade one error against the other.

Why is a calibrated point unsuitable as a validation point?

Solution

Its value influenced the model’s free parameters, so a small residual partly measures fitting capacity. Predictive validation requires withheld data. Calibrated points remain useful for checking implementation of the fitting objective.

In the SXS:BBH:4292 example, the last resolution change is about $8\times10^{-6}$ degrees, while the Lev3 model/window range has width about $1.47\times10^{-2}$ degrees. Which controls the accuracy one may claim?

Solution

The resolution change checks only one source of numerical error. The much larger extraction sensitivity prevents claiming total accuracy of $8\times10^{-6}$ degrees. Even the model/window range is not a rigorous bound: all selected models can share a bias. Report both measures, the chosen windows and models, and the independent sampling check.

References

• Damour, Thibault, Federico Guercilena, Ian Hinder, Seth Hopper, Alessandro Nagar, and Luciano Rezzolla. “Strong-Field Scattering of Two Black Holes: Numerics versus Analytics.” Physical Review D 89 (2014): 081503(R). DOI. Open PDF, arXiv v1.
• Damour, Thibault. “Classical and Quantum Scattering in Post-Minkowskian Gravity.” Physical Review D 102 (2020): 024060. DOI. Open PDF, arXiv v3.
• Damour, Thibault. “Gravitational Scattering, Post-Minkowskian Approximation and Effective One-Body Theory.” Physical Review D 94 (2016): 104015. DOI.
• Damour, Thibault. “Radiative Contribution to Classical Gravitational Scattering at the Third Order in G.” Physical Review D 102 (2020): 124008. DOI. Open PDF.
• Long, Oliver, Harald P. Pfeiffer, Alessandra Buonanno, Gustav Uhre Jakobsen, Gustav Mogull, Antoni Ramos-Buades, Hannes R. Rüter, Lawrence E. Kidder, and Mark A. Scheel. “Highly Accurate Simulations of Asymmetric Black-Hole Scattering and Cross Validation of Effective-One-Body Models.” Physical Review D 112 (2025): 124039. DOI. Open PDF, arXiv v2.
• Pretorius, Frans, and Deepak Khurana. “Black Hole Mergers and Unstable Circular Orbits.” Classical and Quantum Gravity 24 (2007): S83–S108. DOI.
• SXS Collaboration. “Binary Black-Hole Simulation SXS:BBH:3999,” version 3.0 (2025). Dataset DOI.
• SXS Collaboration. “Binary Black-Hole Simulation SXS:BBH:4292,” version 3.0 (2025). Dataset DOI.