Waveform Construction and Cross-Regime Consistency

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 beyond-GR waveform must carry the same theory parameters through compact-body matching, conservative motion, radiation generation, merger modeling, and propagation. Adding separately derived “generation” and “propagation” corrections is safe only after checking that they do not encode the same physical effect twice. The bounded calculation below makes that bookkeeping test executable: one parameter modifies both the inspiral chirp and the propagation phase, and a deliberately duplicated model exposes the resulting bias.

Helpful background. Einstein Equations supplies the GR benchmark; Gravitational Actions and Boundary Terms fixes theory parameters and field normalization; Formulations of General Relativity supplies initial-data and branch distinctions; Relativistic Two-Body Dynamics supplies the adiabatic chirp; Gravitational Radiation and Waveforms supplies frequency-domain phase conventions; and Gravity Beyond General Relativity supplies the shared theory record and cutoff.

One parameter, two physical stages

For a circular nonspinning inspiral, use the Newtonian GR chirp rate with a synthetic $-1$ PN generation correction,

\dot f=K f^{11/3}\left[1+\epsilon \left(\frac{f}{f_{\rm ref}}\right)^{-2/3}\right], \qquad K=\frac{96}{5}\pi^{8/3} \left(\frac{G\mathcal M}{c^3}\right)^{5/3}.

The exponent mimics the relative scaling of dipole flux, but the benchmark is not a waveform claimed by any specific theory. The accumulated time and phase obey

\frac{dt}{df}=\frac1{\dot f}, \qquad \frac{d\Phi_{\rm gen}}{df}=\frac{2\pi f}{\dot f}.

These equations give time and the time-domain accumulated wave phase as functions of instantaneous frequency. A propagation phase multiplies a Fourier amplitude, so it must be combined with the Fourier-domain generation phase. For a positive-frequency analytic strain $h(t)=A(t)e^{-i\Phi_{\rm gen}(t)}$ and $\widetilde h(f)=\int h(t)e^{2\pi ift}\,dt$ , stationary phase gives

\Psi_{\rm gen}(f)=2\pi f\,t(f)-\Phi_{\rm gen}(f), \qquad \frac{d\Psi_{\rm gen}}{df}=2\pi t(f).

We suppress the common constant $-\pi/4$ . The initial conditions $t(f_{\min})=\Phi_{\rm gen}(f_{\min})=0$ set both the initial phase and time origin; they give $\Psi_{\rm gen}(f_{\min})=\Psi_{\rm gen}'(f_{\min})=0$ . Using $\Phi_{\rm gen}(f)$ directly as a Fourier phase would miss the essential $2\pi ft$ term.

The same dimensionless coupling $\epsilon$ enters a synthetic dispersive propagation map,

\delta\Psi_{\rm prop}(f) =\epsilon\kappa\left(\frac{f_{\rm ref}}{f} -\frac{f_{\rm ref}}{f_{\min}}\right).

Both Fourier phases are anchored at $f_{\min}$ . We retain the propagation group delay: its first derivative need not vanish at that boundary. The consistent and deliberately flawed maps, with the same time and phase origins, are

\Psi_{\rm correct} =\Psi_{\rm gen}(\epsilon)+\delta\Psi_{\rm prop}(\epsilon),

\Psi_{\rm double} =\Psi_{\rm correct} +\left[\Psi_{\rm gen}(\epsilon)-\Psi_{\rm GR}\right].

The bracket in the second line is the same generation correction added a second time. Berti, Yagi, and Yunes distinguish generation from propagation and describe the theory-to-parameter mapping in Berti, Yagi, and Yunes 2018, §§ 3.1–3.2, eqs. (19)–(22), pp. 14–18. A theory-specific map must additionally import charges and PN coefficients from a calculation such as Mirshekari and Will 2013, §§ V–VII, pp. 15–31.

Reproducible benchmark

The retained environment is Python 3.9.6, NumPy 2.0.2, and SciPy 1.13.1. Frequencies are detector-frame hertz at zero synthetic redshift; masses and constants are SI inside $K$ . Inputs are $\mathcal M=26M_\odot$ , $f\in[20,200]\,{\rm Hz}$ , $f_{\rm ref}=100\,{\rm Hz}$ , $\epsilon=2.0\times10^{-3}$ , and $\kappa=50\,{\rm rad}$ . Initial data are $t(f_{\min})=\Phi_{\rm gen}(f_{\min})=0$ . DOP853 uses the displayed relative tolerances, componentwise absolute tolerance $10^{-3}$ times the relative tolerance, and maximum step $(f_{\max}-f_{\min})/(N-1)$ .

Save this source exactly as waveform_cross_regime.py:

import sys
import numpy as np
import scipy
from scipy.integrate import quad, solve_ivp

EXPECTED = ("3.9.6", "2.0.2", "1.13.1")
FOUND = (sys.version.split()[0], np.__version__, scipy.__version__)
if FOUND != EXPECTED:
    raise SystemExit(f"version mismatch: expected {EXPECTED}, found {FOUND}")

G_SI = 6.67430e-11
C_SI = 299792458.0
M_SUN_SI = 1.98847e30
MCHIRP = 26.0 * M_SUN_SI
F_MIN, F_MAX, F_REF = 20.0, 200.0, 100.0  # Hz
EPSILON = 2.0e-3
KAPPA_PROP = 50.0  # radians; normalization of the synthetic propagation map
K_CHIRP = (96.0 / 5.0) * np.pi ** (8.0 / 3.0) * (
    G_SI * MCHIRP / C_SI**3
) ** (5.0 / 3.0)


def fdot(frequency, epsilon):
    """SI Hz/s, with a -1PN synthetic generation correction."""
    return K_CHIRP * frequency ** (11.0 / 3.0) * (
        1.0 + epsilon * (frequency / F_REF) ** (-2.0 / 3.0)
    )


def integrate_phase(n_points, rtol, epsilon):
    grid = np.linspace(F_MIN, F_MAX, n_points)

    def rhs(frequency, state):
        rate = fdot(frequency, epsilon)
        return (1.0 / rate, 2.0 * np.pi * frequency / rate)

    solution = solve_ivp(
        rhs,
        (F_MIN, F_MAX),
        (0.0, 0.0),
        method="DOP853",
        t_eval=grid,
        rtol=rtol,
        atol=(rtol * 1.0e-3, rtol * 1.0e-3),
        max_step=(F_MAX - F_MIN) / (n_points - 1),
    )
    if not solution.success:
        raise SystemExit(solution.message)
    return grid, solution.y[0], solution.y[1]


def propagation_phase(grid, epsilon):
    # The phase origin is fixed at F_MIN, so generation and propagation
    # use the same boundary convention.
    return epsilon * KAPPA_PROP * (
        F_REF / grid - F_REF / F_MIN
    )


def gr_phase_exact(grid):
    time = (3.0 / (8.0 * K_CHIRP)) * (
        F_MIN ** (-8.0 / 3.0) - grid ** (-8.0 / 3.0)
    )
    orbital = (6.0 * np.pi / (5.0 * K_CHIRP)) * (
        F_MIN ** (-5.0 / 3.0) - grid ** (-5.0 / 3.0)
    )
    return 2.0 * np.pi * grid * time - orbital


LEVELS = (
    (257, 2.0e-7),
    (513, 5.0e-8),
    (1025, 1.25e-8),
    (2049, 3.125e-9),
)
rows = []
previous_endpoint = None
finest = None
for n_points, rtol in LEVELS:
    f, time_gen, orbital_gen = integrate_phase(n_points, rtol, EPSILON)
    phase_gen = 2.0 * np.pi * f * time_gen - orbital_gen
    _, time_gr, orbital_gr = integrate_phase(n_points, rtol, 0.0)
    phase_gr = 2.0 * np.pi * f * time_gr - orbital_gr
    prop = propagation_phase(f, EPSILON)
    phase_correct = phase_gen + prop
    phase_double = phase_correct + (phase_gen - phase_gr)
    endpoint = phase_correct[-1] - phase_gr[-1]
    change = np.nan if previous_endpoint is None else abs(endpoint - previous_endpoint)
    rms_double_bias = np.sqrt(
        np.trapezoid((phase_double - phase_correct) ** 2, f) / (F_MAX - F_MIN)
    )
    rhs_phase_level = 2.0 * np.pi * f / fdot(f, EPSILON)
    ode_residual_level = np.max(
        np.abs(np.gradient(orbital_gen, f, edge_order=2) - rhs_phase_level)
    ) / np.max(np.abs(rhs_phase_level))
    rows.append((n_points, rtol, endpoint, change, rms_double_bias, ode_residual_level))
    previous_endpoint = endpoint
    finest = (f, phase_gen, phase_gr, phase_correct, phase_double)

f, phase_gen, phase_gr, phase_correct, phase_double = finest
exact_gr = gr_phase_exact(f)
gr_error = np.max(np.abs(phase_gr - exact_gr))
boundary_mismatch = max(abs(phase_correct[0]), abs(phase_gr[0]))
rhs_phase = 2.0 * np.pi * f / fdot(f, EPSILON)
finite_difference = np.gradient(orbital_gen, f, edge_order=2)
ode_residual = np.max(np.abs(finite_difference - rhs_phase)) / np.max(
    np.abs(rhs_phase)
)
gr_limit = np.max(np.abs(phase_gr - exact_gr))
double_identity = np.max(
    np.abs((phase_double - phase_correct) - (phase_gen - phase_gr))
)

# Independent stationary-phase integral: integrate Psi''=2*pi/fdot twice
# with Psi(F_MIN)=Psi'(F_MIN)=0. This catches a time/Fourier phase swap.
psi_quad = quad(
    lambda frequency: 2.0 * np.pi * (F_MAX - frequency)
    / fdot(frequency, EPSILON), F_MIN, F_MAX,
    epsabs=1.0e-10, epsrel=1.0e-12,
)[0]
spa_error = abs(phase_gen[-1] - psi_quad)

print("versions:", *FOUND)
print("N     rtol       delta_psi_end[rad]   change[rad]   rms_double_bias[rad]   ODE_residual")
for n_points, rtol, endpoint, change, bias, residual in rows:
    change_text = "---" if np.isnan(change) else f"{change:.3e}"
    print(
        f"{n_points:4d}  {rtol:.3e}  {endpoint: .9f}          "
        f"{change_text:>9s}       {bias:.9f}              {residual:.3e}"
    )
print(f"analytic_GR_max_error_rad={gr_error:.3e}")
print(f"normalized_ODE_residual={ode_residual:.3e}")
print(f"boundary_mismatch_rad={boundary_mismatch:.3e}")
print(f"GR_limit_error_rad={gr_limit:.3e}")
print(f"double_count_identity_error_rad={double_identity:.3e}")
print(f"stationary_phase_integral_error_rad={spa_error:.3e}")

tests = {
    "endpoint convergence": rows[-1][3] < 2.0e-7,
    "analytic GR benchmark": gr_error < 2.0e-6,
    "ODE residual convergence": (
        all(rows[i + 1][5] < rows[i][5] for i in range(len(rows) - 1))
        and ode_residual < 2.0e-4
    ),
    "boundary convention": boundary_mismatch < 1.0e-12,
    "GR limit": gr_limit < 2.0e-6,
    "stationary-phase integral": spa_error < 2.0e-8,
    "double-count control": rows[-1][4] > 0.1 and double_identity < 1.0e-10,
}
failed = [name for name, passed in tests.items() if not passed]
if failed:
    raise SystemExit("FAILED: " + ", ".join(failed))
print("PASS:", ", ".join(tests))

The program must produce the refinement table, diagnostic scalars, and final PASS line shown below. Rerun from a clean directory containing that file:

set -eu
test -f waveform_cross_regime.py
python3 waveform_cross_regime.py

The independently rerun output is:

$N$	relative tolerance	endpoint theory–GR phase (rad)	successive change (rad)	RMS duplicated-correction bias (rad)	normalized ODE defect
257	$2.000\times10^{-7}$	$-5.275253694$	—	$2.721781725$	$3.566\times10^{-3}$
513	$5.000\times10^{-8}$	$-5.275253694$	$3.411\times10^{-13}$	$2.721773164$	$9.453\times10^{-4}$
1025	$1.250\times10^{-8}$	$-5.275253694$	$5.343\times10^{-12}$	$2.721771024$	$2.435\times10^{-4}$
2049	$3.125\times10^{-9}$	$-5.275253694$	$2.274\times10^{-13}$	$2.721770489$	$6.181\times10^{-5}$

The analytic GR time and Fourier phase in the same boundary convention are

\begin{aligned} t_{\rm GR}(f)&=\frac{3}{8K} \left(f_{\min}^{-8/3}-f^{-8/3}\right),\\ \Psi_{\rm GR}(f)&=2\pi f t_{\rm GR}(f) -\frac{6\pi}{5K}\left(f_{\min}^{-5/3}-f^{-5/3}\right). \end{aligned}

The maximum GR discrepancy is $5.798\times10^{-12}\,{\rm rad}$ in this rerun, the initial phase mismatch is zero, and the duplicated-correction identity closes in floating-point arithmetic. The independent quadrature of $\Psi_{\rm gen}''=2\pi/\dot f$ checks the non-GR Fourier phase as well. Last-bit roundoff may vary across machines. The finite-difference defect of the integrated time-domain phase falls by nearly four on each grid doubling, as expected for its second-order derivative check. It measures the sampled derivative, not the adaptive solver’s continuum error. The endpoint is already at roundoff; its irregular successive changes do not establish an ODE convergence order.

Cross-regime assembly contract

A production map begins with one theory record: action and physical matter metric; background and strong-field branch; compact-body masses, charges, and spins; PN and EFT cutoffs; and the cosmology used in propagation. A schematic detector-frame waveform is

\widetilde h(f)=\frac{\mathcal A_{\rm src}(f;\theta)}{d_L^{\rm GW}} \exp\!\left\{i\left[\Psi_{\rm GR}+\delta\Psi_{\rm gen} +\delta\Psi_{\rm strong}+\delta\Psi_{\rm prop}\right]\right\},

where every correction is a function of the same parameter vector $\theta$ . If PN and strong-field models overlap, blending weights must sum to one and the phase and its frequency derivative must agree across the matching window. A merger correction calibrated on a different solution branch, matter coupling, or asymptotic scalar value cannot be inserted as though it shared $\theta$ . When no controlled merger calculation exists, the honest product is an inspiral-only waveform with an explicit termination frequency, not an unqualified inspiral–merger model.

Scientific interpretation

At $200\,{\rm Hz}$ the correct Fourier-phase theory–GR shift is $-5.275253694\,{\rm rad}$ . Of this, $-0.45\,{\rm rad}$ is the anchored propagation term and $-4.825253694\,{\rm rad}$ is the generation term. Duplicating generation changes the endpoint to $-10.100507388\,{\rm rad}$ and creates a $2.721770489\,{\rm rad}$ band-averaged RMS discrepancy at the finest displayed grid. Refinement cannot cure that error because it is a formulation error. These values use the stated initial-frequency alignment; they are not detector-weighted mismatches or parameter biases after fitting time and phase offsets.

This benchmark tests bookkeeping, convergence, a boundary convention, and the GR limit. It does not model spins, eccentricity, strong-field charge evolution, merger, ringdown, cosmological redshift, detector calibration, or parameter estimation. A physical waveform must match the isolated-body branch to PN charges, use a common source/detector-frame map, transition to a merger model without overlap or gaps, and keep every frequency below the PN and field-theory cutoffs. Berti, Yagi, and Yunes emphasize that generic deformation bases are useful but do not replace a theory-specific map in Berti, Yagi, and Yunes 2018, §§ 3 and 5, pp. 12–19 and 30–32.

Production calculations must diagnose more failure modes than this smooth scalar ODE: wrong strong-field branches, boundary contamination, stiffness, loss of hyperbolicity, unresolved gradients, cancellation between sectors, solver stagnation, spurious eigenvalues in mode calculations, cutoff violation, and sensitivity to initial guesses. A small phase change under resolution refinement does not detect those failures by itself.

Exercises

Choose the consistent formulation. At the upper boundary, the rerun finds a Fourier generation shift $-4.825253694\,{\rm rad}$ and a propagation shift $-0.45\,{\rm rad}$ from the same $\epsilon$ . Compare (A) one shared parameter applied once in each sector, (B) the same parameter with generation added twice, and (C) two independent generation and propagation parameters.

Solution

Formulation A gives $-4.825253694-0.45=-5.275253694\,{\rm rad}$ and is the common-parameter map declared here. Formulation B gives $2(-4.825253694)-0.45=-10.100507388\,{\rm rad}$ , an endpoint error of $-4.825253694\,{\rm rad}$ . Formulation C can be a legitimate agnostic null-test basis, but it is not this prediction unless a theory permits the two coefficients to vary independently. When the action fixes both from one $\epsilon$ , comparison with data must preserve that relation rather than fit two couplings and interpret them as independent theory parameters.

References

Berti, Emanuele, Kent Yagi, and Nicolás Yunes. “Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences: (I) Inspiral–Merger.” General Relativity and Gravitation 50 (2018): 46. DOI; Open PDF.
Mirshekari, Saeed, and Clifford M. Will. “Compact Binary Systems in Scalar-Tensor Gravity: Equations of Motion to 2.5 Post-Newtonian Order.” Physical Review D 87 (2013): 084070. DOI; Open PDF.