Waveform Hybridization

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.

Waveform hybridization constructs one long signal from an early-time analytic approximant and a late-time numerical-relativity (NR) waveform. The operation is not a splice at a visually convenient point: the inputs must describe the same intrinsic source, use compatible modes and frames, be aligned by declared freedoms, overlap where both are credible, and be blended with controlled differentiability. This chapter consumes dynamics and NR products supplied elsewhere. Its responsibility is the construction and its error budget, not a new binary evolution or a detector-specific inference verdict.

Helpful background. Einstein Equations supplies the common spacetime theory; Waves, Green Functions, and Geometric Optics supplies phase, Fourier, and windowing methods; Relativistic Two-Body Dynamics supplies analytic inspiral dynamics; Computational Relativity supplies the NR waveform product; and Gravitational-Wave Data Analysis and Inference supplies downstream application-specific norms.

Compatible inputs before alignment

The analytic and NR inputs must agree on masses, spin vectors at a declared reference time, eccentricity assumptions, retarded-time units, complex strain $h=h_+-ih_\times$ , unit-normalized spin-weight $-2$ harmonics, and the $e^{-i\omega t}$ convention. They must also represent the same amplitude normalization. One convenient source-frame choice is the distance-rescaled mode $\widehat H_{\ell m}=H_{\ell m}/(GM_{\rm src})$ tabulated against time in units of $GM_{\rm src}$ ; an analytic detector strain cannot be aligned directly to an NR product tabulated as $rh_{\ell m}/M$ until both are converted to a common object. Record whether modes are in an asymptotically inertial, source, coprecessing, or Bondi frame and translate any residual BMS transformation before fitting time and phase. A rigid time and phase shift cannot repair a boost, supertranslation, normalization mismatch, different spin-weighted-harmonic phase, or unremoved center-of-mass drift.

The overlap interval must begin after junk radiation and end before the analytic approximant loses validity. It must also contain enough chirp to break the near-degeneracy between time and phase. MacDonald, Nissanke, and Pfeiffer formulate this as a frequency interval mapped to a time interval and show why a window shorter than about a cycle makes the recovered time shift sensitive to small phase oscillations MacDonald, Nissanke, and Pfeiffer 2011, §§3.3 and 4.1–4.2, pp. 13–17, PDF.

For precessing or multimode data, alignment can require a three-dimensional frame rotation as well as time and orbital phase. The rotation must act coherently on all $m$ at fixed $\ell$ ; independently phasing every mode can erase physical relations. Varma and collaborators give a concrete multimode frame alignment and blending construction for hybridized NR waveforms in Varma et al. 2019, §§IV–V, PDF.

Alignment and blend

For the nonprecessing $(2,2)$ example, let $h_A$ and $h_N$ denote the analytic and NR representations of that same normalized mode. Determine $(\Delta t,\Delta\varphi)$ by minimizing a declared residual over $[t_1,t_2]$ ,

\mathcal R^2(\Delta t,\Delta\varphi) =\frac{1}{t_2-t_1} \int_{t_1}^{t_2} \left\lvert h_A(t+\Delta t)e^{i\Delta\varphi}-h_N(t) \right\rvert^2dt.

Here $\Delta\varphi$ is the waveform-mode phase; an orbital phase rotation would enter an $(\ell,m)$ mode as $e^{-im\Delta\Phi}$ . The fit norm, amplitude weighting, optimizer bounds, interpolation, modes, and frame are part of the answer. MacDonald et al. instead minimize a phase residual and then use a cosine blend MacDonald, Nissanke, and Pfeiffer 2011, §3.3, Eqs. (20)–(23), pp. 13–14, PDF; neither choice is universally privileged.

With $x=(t-t_1)/(t_2-t_1)$ and the quintic smoothstep

w(x)=6x^5-15x^4+10x^3,

the hybrid is

h_H(t)= \begin{cases} h_A^{\mathrm{aligned}}(t),&t\lt t_1,\\ [1-w(x)]h_A^{\mathrm{aligned}}(t)+w(x)h_N(t),&t_1\le t\le t_2,\\ h_N(t),&t\gt t_2. \end{cases}

Because $w' = w''=0$ at both endpoints, this construction preserves two derivatives when both inputs are $C^2$ . That formal smoothness does not imply accuracy: if the inputs disagree coherently, the blend merely distributes the disagreement.

It also changes derivative-based observables inside the overlap. Direct differentiation gives

\dot h_H=(1-w)\dot h_A+w\dot h_N +\dot w(h_N-h_A),

where $h_A$ now denotes the aligned input. The last term can create an artificial contribution to news and energy flux when the inputs disagree, even though the strain looks smooth. Compare derivative quantities as well as strain, and never silently fill shifted samples beyond an input’s support with a constant endpoint value.

Controlled synthetic benchmark

The deterministic test uses one $(2,2)$ mode in an inertial frame; numerical $t$ below denotes the dimensionless coordinate $u/(GM)$ . The exact family is

h_*(t)=\left(0.12+0.88x^2\right) \exp\left\{-2i\left[0.018(t+1000)+1.1\times10^{-5}(t+1000)^2\right]\right\},

with $x=(t+1000)/1100$ clipped to $[0,1]$ . The analytic input receives a growing $0.10x^5$ phase perturbation and $0.003x^3$ amplitude perturbation. The NR input receives early-decaying perturbations, then an injected coordinate offset $\tau=3.25GM$ and mode phase $0.37$ rad. The baseline window is $[-500,-300]GM$ ; 2,000 uniformly spaced shift trials cover $[-4.5,-2]GM$ at 4,001 samples. The phase is minimized analytically for each trial. The norm is the time-domain complex $L^2$ norm, not a detector match.

The code below was run with Python 3.12.14 and NumPy==2.3.5, with no random state or external inputs. Shifted samples outside the analytic record are marked unavailable, and every relative-error comparison uses samples within the same supported interval $-995.5\le t\le100$ , valid for all trial shifts. The trapezoidal integral starts at the first retained sample, so its lower endpoint can differ from $-995.5$ by less than one time step; this is part of the discretization error tested by refinement. The sampling test holds the physics, window, blend, and shift bounds fixed:

Grid samples	Boundary derivative residual
2001	3.672277e-2
4001	1.835653e-2
8001	9.178899e-3

The first-order decrease is the expected finite-difference truncation behavior. It checks implementation continuity; it does not estimate modeling error.

The program also emits this sensitivity table:

Variation	Recovered shift $\Delta t/(GM)$	Recovered phase (rad)	Alignment RMS	Hybrid relative error
Baseline	-3.89594797	0.33393807	$7.154006\times10^{-4}$	$1.433480\times10^{-3}$
Wider window	-4.15857929	0.31813641	$9.047054\times10^{-4}$	$2.203663\times10^{-3}$
Linear blend	-3.89594797	0.33393807	$7.154006\times10^{-4}$	$1.434619\times10^{-3}$
Coarse sampling	-3.85185185	0.33670966	$7.166770\times10^{-4}$	$1.279681\times10^{-3}$
Doubled analytic error	-4.50000000	0.30163926	$1.209132\times10^{-3}$	$2.478658\times10^{-3}$
Doubled NR error	-3.99974987	0.32624040	$9.366536\times10^{-4}$	$2.119996\times10^{-3}$

The baseline does not recover the injected $-3.25GM$ shift exactly: analytic and NR phase errors bias the optimizer by $0.64595GM$ . The doubled-analytic case reaches the search boundary, a visible optimizer failure. Window and approximant variation are therefore physical-systematics tests, not optional plotting choices.

Executable rerun

import numpy as np

TAU_INJECTED = 3.25
PHASE_INJECTED = 0.37

def exact_waveform(t):
    x = np.clip((t + 1000.0) / 1100.0, 0.0, 1.0)
    amplitude = 0.12 + 0.88 * x**2
    orbital_phase = 0.018*(t + 1000.0) + 1.1e-5*(t + 1000.0)**2
    return amplitude * np.exp(-2j * orbital_phase)

def blend_weight(x, kind):
    if kind == "linear":
        return x
    if kind == "cubic":
        return 3*x**2 - 2*x**3
    return 6*x**5 - 15*x**4 + 10*x**3

def align(pn_time, pn, nr_time, nr, window, trials):
    use = (nr_time >= window[0]) & (nr_time <= window[1])
    target = nr[use]
    shifts = np.linspace(-4.5, -2.0, trials)
    best = (np.inf, 0.0, 0.0)
    for shift in shifts:
        query = nr_time[use] + shift
        if query.min() < pn_time[0] or query.max() > pn_time[-1]:
            raise ValueError("Alignment query lies outside the analytic record")
        candidate = np.interp(query, pn_time, pn.real)
        candidate = candidate + 1j*np.interp(query, pn_time, pn.imag)
        phase = np.angle(np.vdot(candidate, target))
        residual = np.mean(abs(candidate*np.exp(1j*phase) - target)**2)
        if residual < best[0]:
            best = (residual, shift, phase)
    return best[1], best[2], np.sqrt(best[0])

def construct(samples=4001, window=(-500.0, -300.0), kind="quintic",
              pn_error=1.0, nr_error=1.0):
    pn_time = np.linspace(-1000.0, 100.0, samples)
    nr_time = pn_time.copy()
    x_pn = np.clip((pn_time + 1000.0) / 1100.0, 0.0, 1.0)
    pn = exact_waveform(pn_time)
    pn *= (1 + pn_error*0.003*x_pn**3)*np.exp(-1j*pn_error*0.10*x_pn**5)

    physical_nr_time = nr_time - TAU_INJECTED
    target = exact_waveform(physical_nr_time)*np.exp(1j*PHASE_INJECTED)
    x_nr = np.clip((physical_nr_time + 650.0) / 750.0, 0.0, 1.0)
    nr = target*(1 + nr_error*0.001*(1-x_nr)**3)
    nr *= np.exp(-1j*nr_error*0.006*(1-x_nr)**4)

    shift, phase, rms = align(pn_time, pn, nr_time, nr, window,
                              max(301, samples//2))
    pn_aligned = np.interp(nr_time+shift, pn_time, pn.real, left=np.nan, right=np.nan)
    pn_aligned = pn_aligned + 1j*np.interp(
        nr_time+shift, pn_time, pn.imag, left=np.nan, right=np.nan)
    pn_aligned *= np.exp(1j*phase)
    scaled = np.clip((nr_time-window[0])/(window[1]-window[0]), 0.0, 1.0)
    weight = blend_weight(scaled, kind)
    hybrid = np.where(weight == 1, nr, (1-weight)*pn_aligned + weight*nr)
    comparison = nr_time >= -995.5
    assert np.all(np.isfinite(hybrid[comparison]))
    error = np.sqrt(np.trapezoid(abs(hybrid[comparison]-target[comparison])**2,
                                nr_time[comparison])
                    / np.trapezoid(abs(target[comparison])**2, nr_time[comparison]))

    derivative = np.gradient(hybrid, nr_time, edge_order=2)
    jumps = []
    for edge in window:
        index = np.argmin(abs(nr_time-edge))
        left = (hybrid[index]-hybrid[index-1])/(nr_time[index]-nr_time[index-1])
        right = (hybrid[index+1]-hybrid[index])/(nr_time[index+1]-nr_time[index])
        jumps.append(abs(right-left)/max(abs(derivative[index]), 1e-15))
    return shift, phase, rms, error, max(jumps)

print("samples,boundary_derivative_residual")
for samples in (2001, 4001, 8001):
    shift, phase, rms, error, jump = construct(samples=samples)
    print(f"{samples},{jump:.6e}")

print("case,shift,phase,alignment_rms,hybrid_relative_error,boundary_residual")
cases = [
    ("baseline", 4001, (-500., -300.), "quintic", 1., 1.),
    ("wide_window", 4001, (-550., -250.), "quintic", 1., 1.),
    ("linear_blend", 4001, (-500., -300.), "linear", 1., 1.),
    ("coarse_sampling", 2001, (-500., -300.), "quintic", 1., 1.),
    ("larger_PN_error", 4001, (-500., -300.), "quintic", 2., 1.),
    ("larger_NR_error", 4001, (-500., -300.), "quintic", 1., 2.),
]
for name, samples, window, kind, pn_error, nr_error in cases:
    values = construct(samples, window, kind, pn_error, nr_error)
    shift, phase, rms, error, jump = values
    print(f"{name},{shift:.8f},{phase:.8f},{rms:.6e},{error:.6e},{jump:.6e}")

Save the Python program above as hybrid_benchmark.py. The program emits both tables. The exact rerun is:

python3.12 -m venv .venv
. .venv/bin/activate
python3 -m pip install numpy==2.5.3
python3 hybrid_benchmark.py

Error budget and stopping conditions

Vary, at minimum, analytic approximant and formal order, NR resolution and extraction product, matching-window location and width, blend, sample rate and interpolation, mode set, and alignment bounds. Keep analytic truncation, NR error, alignment bias, blend error, mode omission, and frame uncertainty as separate terms before any combined norm. MacDonald et al. find in their equal-mass nonspinning study that PN-approximant differences dominate several numerical effects, but that conclusion is conditional on their sources, PSDs, mass scaling, alignment, and mode content MacDonald, Nissanke, and Pfeiffer 2011, §§4.3–5, pp. 17–23, PDF.

Stop if the optimum reaches a bound, multiple minima exchange under small perturbations, the recovered shift changes substantially with window, the inputs disagree outside their estimated errors, an endpoint derivative fails to converge, or important modes cannot be put in one frame. If a detector mismatch is later quoted, declare the one-sided PSD, frequency band, total-mass scaling, modes, frame, taper, and maximized shifts. Such a number is conditional evidence, not an inference verdict.

Exercise

Verify the endpoint differentiability of the quintic blend. Then use the baseline row to compute the recovered time- and phase-shift biases relative to the injected values $-3.25GM$ and $0.37$ rad.

Solution

Direct differentiation gives

w'(x)=30x^2(x-1)^2, \qquad w''(x)=60x(2x^2-3x+1).

Hence $w(0)=0$ , $w(1)=1$ , and both $w'$ and $w''$ vanish at $0$ and $1$ . The piecewise hybrid is therefore $C^2$ at the endpoints when both inputs are $C^2$ . A linear blend has nonzero endpoint slope and does not provide this guarantee.

The recovered biases are

\left\lvert-3.89594797-(-3.25)\right\rvert =0.64594797\,GM,

\left\lvert0.33393807-0.37\right\rvert =0.03606193\ \mathrm{rad}.

These are much larger than the shift-grid spacing. They arise from the deliberately inconsistent analytic and NR perturbations, showing that an optimizer can be numerically converged yet physically biased.

References

MacDonald, Ilana, Samaya Nissanke, and Harald P. Pfeiffer. “Suitability of Post-Newtonian/Numerical-Relativity Hybrid Waveforms for Gravitational Wave Detectors.” Classical and Quantum Gravity 28 (2011): 134002. DOI. Open PDF.
Varma, Vijay, Scott E. Field, Mark A. Scheel, Jonathan Blackman, Lawrence E. Kidder, and Harald P. Pfeiffer. “Surrogate Model of Hybridized Numerical Relativity Binary Black Hole Waveforms.” Physical Review D 99 (2019): 064045. DOI. Open PDF.