Spinning Effective-to-Backwards One Body ($\texttt{SEBOB}$): combining Effective One-Body inspirals and Backwards One-Body merger-ringdowns for aligned spin black hole binaries

TL;DR

The paper introduces SEBOB, a hybrid gravitational waveform model combining Effective-One-Body and Backwards-One-Body approaches, aiming for high accuracy, physical insight, and reduced reliance on numerical relativity data.

Contribution

SEBOB is a novel hybrid framework that integrates EOB and BOB models, enhancing physical interpretability and reducing calibration dependence in gravitational waveform modeling.

Findings

SEBOB achieves accuracy comparable to state-of-the-art models.

The BOB-informed NQC corrections improve waveform smoothness.

SEBOB is computationally efficient and extensible.

Abstract

High-fidelity gravitational waveform models are essential for realizing the scientific potential of next-generation gravitational-wave observatories. While highly accurate, state-of-the-art models often rely on extensive phenomenological calibrations to numerical relativity (NR) simulations for the late-inspiral and merger phases, which can limit physical insight and extrapolation to regions where NR data is sparse. To address this, we introduce the Spinning Effective-to-Backwards One Body (SEBOB) formalism, a hybrid approach that combines the well-established Effective-One-Body (EOB) framework with the analytically-driven Backwards-One-Body (BOB) model, which describes the merger-ringdown from first principles as a perturbation of the final remnant black hole. We present two variants building on the state-of-the-art $\texttt{SEOBNRv5HM}$ model: $\texttt{seobnrv5_nrnqc_bob}$, which retains standard NR-calibrated non-quasi-circular (NQC) corrections and attaches a BOB-based merger-ringdown; and a more ambitious variant, $\texttt{seobnrv5_bob}$, which uses BOB to also inform the NQC corrections, thereby reducing reliance on NR fitting and enabling higher-order ($\mathcal{C}^2$) continuity by construction. Implemented in the open-source $\texttt{NRPy}$ framework for optimized C-code generation, the SEBOB model is transparent, extensible, and computationally efficient. By comparing our waveforms to a large catalog of NR simulations, we demonstrate that SEBOB yields accuracies comparable to the highly-calibrated $\texttt{SEOBNRv5HM}$ model, providing a viable pathway towards more physically motivated and robust waveform models

Tables

Table 1. Table 1 : Summary of NR inputs for the (2,2) mode of the SEOBNRv5 waveform model. The first three parameters are calibration parameters, tuned via Bayesian inference against NR waveforms. The remaining parameters are derived from NR fits and are specified, prior to calibration, for NQC corrections or the merger-ringdown model. Nominal values for these parameters are found in [ 15 ] , except for M f M_{f} and a f a_{f} , which are given in the listed citations in the Method column. The abbreviation NBI stands for Nested Bayesian Inference.

Parameter	Description	Pre Calibration	Method	In seobnrv5_bob
$a_6$	Effective 5PN non-spinning inspiral coefficient	-	NBI	✓
$d_{\mathrm{SO}}$	Effective 4.5PN inspiral spin-orbit coefficient	-	NBI	✓
$\Delta t^{\mathrm{ISCO}}$	Time shift between EOB ISCO passage and (2,2) peak strain	-	NBI	✓
${\left|h_{22}^{\mathrm{NR}}\right|}^{\mathrm{peak}}$	(2,2) mode peak strain amplitude	✓	NR fit	✓
${\ddot{\left|h\right|}}_{22}^{\mathrm{peak}}$	(2,2) mode strain amplitude 2nd deriv. at peak strain time	✓	NR fit	-
${\omega }_{22}^{\mathrm{peak}}$	(2,2) mode frequency at peak strain time	✓	NR fit	✓
${\dot{\omega }}_{22}^{\mathrm{peak}}$	(2,2) mode frequency deriv. at peak strain time	✓	NR fit	-
$M_f$	Remnant Mass	✓	NR fit [26]	✓
$a_f$	Remnant Spin	✓	NR fit [27]	✓
$c_{1,f}$	Merger-ringdown amplitude coefficient	✓	NR fit	-
$c_{2,f}$	Merger-ringdown amplitude coefficient	✓	NR fit	-
$d_{1,f}$	Merger-ringdown phase coefficient	✓	NR fit	-
$d_{2,f}$	Merger-ringdown phase coefficient	✓	NR fit	-

Table 2. Table 2 : Frequency of ODE stop conditions for the SEOBNRv5 waveform model. The conditions are calculated over 5000 randomly sampled points with 1 ≤ q ≤ 20 1\leq q\leq 20 , − 0.999 ≤ χ i ≤ .999 -0.999\leq\chi_{i}\leq.999 , and 5 × 10 − 3 ≤ Ω start ≤ 1.5 × 10 − 2 5\times 10^{-3}\leq\Omega_{\mathrm{start}}\leq 1.5\times 10^{-2} . The frequency for the 607 quasi-circular, spin-aligned, black hole binaries simulated by the SXS catalog is given in brackets. For each stop condition, we specify how often the ODE was terminated before remnant ISCO and/or t 22 p ​ e ​ a ​ k t_{22}^{peak} crossing. The 2 non ISCO crossings when r < r stop < r ISCO r<r_{\mathrm{stop}}<r_{\mathrm{ISCO}} , happens because d ​ t ODE ≈ .1 ≡ d ​ t sample dt_{\mathrm{ODE}}\approx.1\equiv dt_{\mathrm{sample}} at last time leading to exclusion of r ISCO r_{\mathrm{ISCO}} in the fine sampling.

Condition	Total	ISCO	Peak
$\Omega$ peak	3220(369)	78(0)	530(20)
$p_{r_{\ast }}$ peak	1080(103)	0(0)	0(0)
$\dot{r}>0$	6(0)	0(0)	0(0)
	684(132)	2(0)	0(0)
and $\Omega >1$	10(3)	0(0)	0(0)

Equations

$${\mathcal{H}}_{\mathrm{real}}=\frac{1}{\eta }\sqrt{1+2\eta \left({\mathcal{H}}_{\mathrm{eff}}-1\right)}.$$

$$g^{00}{\mathcal{H}}_{\mathrm{eff}}^2+g^{0i}{\mathcal{H}}_{\mathrm{eff}}{p}_i+g^{ij}{p}_i{p}_j+\mathcal{Q}(x^j,p_j)=-1.$$

$${\mathcal{H}}_{\mathrm{eff}}=H_{\mathrm{odd}}+H_{\mathrm{even}},$$

$$H_{\mathrm{even}}{|}_{a_i\to 0}=\sqrt{A_{\mathrm{noS}}\left(1+A_{\mathrm{noS}}{D}_{\mathrm{noS}}{p}_{r_{\ast }}^2+\frac{p_{\varphi }^2}{r^2}+Q_{\mathrm{noS}}\right)}.$$

$$\xi =A_{\mathrm{noS}}\sqrt{D_{\mathrm{noS}}},p_{r_{\ast }}=\xi p_r,\frac{d{r}_{\ast }}{dr}=\frac{1}{\xi }.$$

$${\mathcal{F}}_{\varphi }$$

$${\mathcal{F}}_r$$

$$\frac{dr}{dt}=\xi \frac{\partial {\mathcal{H}}_{\mathrm{real}}}{\partial p_{r_{\ast }}},$$

$$\frac{d\varphi }{dt}=\frac{\partial {\mathcal{H}}_{\mathrm{real}}}{\partial p_{\varphi }}\equiv \Omega ,$$

$$\frac{d{p}_{r_{\ast }}}{dt}=-\xi \frac{\partial {\mathcal{H}}_{\mathrm{real}}}{\partial r}+{\mathcal{F}}_r,$$

$$\frac{d{p}_{\varphi }}{dt}={\mathcal{F}}_{\varphi }.$$

$$\begin{gathered} h_{22}^{\mathrm{insp}-\mathrm{plunge}}=h_{22}^{\mathrm{F}}\left[1+\frac{p_{r_{\ast }}^2}{r^2{\Omega }^2}\left(a_1+\frac{a_2}{r}+\frac{a_3}{r^{3/2}}\right)\right]\times \\ {e}^{i\left(b_1\frac{p_{r_{\ast }}}{r\Omega }+b_2\frac{p_{r_{\ast }}^3}{r\Omega }\right)}. \end{gathered}$$

$$t_{22}^{\mathrm{peak}}=t^{\mathrm{ISCO}}-\Delta t^{\mathrm{ISCO}}.$$

$$\left(\begin{array}{ccc}{\mathcal{Q}}_1& {\mathcal{Q}}_2& {\mathcal{Q}}_3\\ {\dot{\mathcal{Q}}}_1& {\dot{\mathcal{Q}}}_2& {\dot{\mathcal{Q}}}_3\\ {\ddot{\mathcal{Q}}}_1& {\ddot{\mathcal{Q}}}_2& {\ddot{\mathcal{Q}}}_3\end{array}\right)\left(\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right)={\left(\begin{array}{c}|h_{22}^{\mathrm{NR}}|-|h_{22}^{\mathrm{F}}|\\ \frac{d}{dt}|h_{22}^{\mathrm{NR}}|-\frac{d}{dt}|h_{22}^{\mathrm{F}}|\\ \frac{d^2}{d{t}^2}|h_{22}^{\mathrm{NR}}|-\frac{d^2}{d{t}^2}|h_{22}^{\mathrm{F}}|\end{array}\right)}_{t_{22}^{\mathrm{peak}}},$$

$$\mathcal{Q}={\left[|h_{22}^{\mathrm{F}}|\left\{\frac{p_{r_{\ast }}^2}{r^2{\Omega }^2},\frac{p_{r_{\ast }}^2}{r^3{\Omega }^2},\frac{p_{r_{\ast }}^2}{r^{7/2}{\Omega }^2}\right\}\right]}_{t=t_{22}^{\mathrm{peak}}}.$$

$$\left(\begin{array}{cc}{\mathcal{P}}_1& {\mathcal{P}}_2\\ {\dot{\mathcal{P}}}_1& {\dot{\mathcal{P}}}_2\end{array}\right)\left(\begin{array}{c}{b}_1\\ b_2\end{array}\right)={\left(\begin{array}{c}{\omega }_{22}^{\mathrm{NR}}-{\omega }_{22}^{\mathrm{F}}\\ {\dot{\omega }}_{22}^{\mathrm{NR}}-{\dot{\omega }}_{22}^{\mathrm{F}}\end{array}\right)}_{t_{22}^{\mathrm{peak}}},$$

$$\mathcal{P}={\left[\left\{\frac{p_{r_{\ast }}}{r\Omega },\frac{p_{r_{\ast }}^3}{r\Omega }\right\}\right]}_{t=t_{22}^{\mathrm{peak}}}.$$

$$h_{22}^{\mathrm{merger}-\mathrm{ringdown}}(t)={\tilde{A}}_{22}(t)e^{i{\tilde{\varphi }}_{22}(t)}{e}^{-i{\sigma }_{22}^{\mathrm{QNM}}(t-t_{22}^{\mathrm{peak}})},$$

$${\tilde{A}}_{22}(t)=\eta c_{1,c}\tanh \left[c_{1,f}\left(t-t_{22}^{\mathrm{peak}}\right)+c_{2,f}\right]+c_{2,c},$$

$${\tilde{\varphi }}_{22}={\varphi }_0-d_{1,c}\log \left[\frac{1+d_{2,f}{e}^{-d_{1,f}\left(t-t_{22}^{\mathrm{peak}}\right)}}{1+d_{2,f}}\right],$$

$${\sigma }_{22}^{\mathrm{QNM}}=-\frac{i}{{\tau }_{22}^{\mathrm{QNM}}}+{\omega }_{22}^{\mathrm{QNM}}.$$

$$A_{\mathrm{noS}}=1-2u+2\eta u^{3}\\ +\eta\left(\frac{94}{3}-\frac{41\pi^{2}}{32}\right)u^{4}\\ +\Bigg{[}\eta\left(\frac{2275\pi^{2}}{512}-\frac{4237}{60}+\frac{128\gamma_{\mathrm{E}}}{5}+\frac{256\ln 2}{5}\right)\\ +\eta^{2}\left(\frac{41\pi^{2}}{32}-\frac{221}{6}\right)\\ +\frac{64}{5}\eta\ln u\Bigg{]}u^{5}\\ +\left[\eta a_{6}+\left(-\frac{144\eta^{2}}{5}-\frac{7004\eta}{105}\right)\ln u\right]u^{6},$$

$$H_{\mathrm{odd}}=\frac{p_{\varphi }\left(g_{a_{+}}{a}_{+}+g_{a_{-}}{a}_{-}\right)+\eta d_{\mathrm{SO}}{u}^3{p}_{\varphi }{a}_{+}+G_{a^3}^{\mathrm{align}}}{r^3+a_{+}^2\left(r+2\right)},$$

$$A_{lm}(t)\equiv |{\psi }_{4,lm}|=A_{p,lm}\mathrm{sech}\left(\frac{t-t_{p,lm}}{{\tau }_{lm}}\right),$$

$$|{\mathcal{N}}_{lm}{|}^2\propto m^2\frac{d{\Omega }^2}{dt}.$$

$$\omega (t)={\left\{{\omega }_0^4+k\left(\tanh \left(\frac{t-t_p}{\tau }\right)-\tanh \left(\frac{t_0-t_p}{\tau }\right)\right)\right\}}^{1/4},$$

$$k=\frac{{\omega }_{\mathrm{QNM}}^4-{\omega }_0^4}{1-\tanh \left(\frac{t_0-t_p}{\tau }\right)},$$

$$\begin{gathered} \varphi (t)={\varphi }_0+{\omega }_{+\infty }\tau \left({\arctan }_{+}+{\mathrm{arctanh}}_{+}\right) \\ -{\omega }_{-\infty }\tau \left({\arctan }_{-}+{\mathrm{arctanh}}_{-}\right), \end{gathered}$$

$${\arctan }_{\pm }$$

$${\mathrm{arctanh}}_{\pm }$$

$$h(t)=\frac{|h_0|{\omega }_0^2}{\omega (t{)}^2}\frac{\mathrm{sech}\left(\frac{t-t_p}{\tau }\right)}{\mathrm{sech}\left(\frac{t_0-t_p}{\tau }\right)}\exp \left(i\varphi (t)\right),$$

$$t_p=t_0-2\tau \ln \left(\frac{{\omega }_0}{{\omega }_{\mathrm{QNM}}}\right),$$

$${\tanh }^{-1}(x)=\frac{1}{2}\ln \left(\frac{1+x}{1-x}\right)$$

$$\frac{d{p}_r}{dt}\equiv \frac{\partial {\mathcal{H}}_{\mathrm{real}}(r,p_{r_{\ast }}=0,p_{\varphi })}{\partial r}$$

$$\frac{d\varphi }{dt}\equiv \frac{\partial {\mathcal{H}}_{\mathrm{real}}(r,p_{r_{\ast }}=0,p_{\varphi })}{\partial p_{\varphi }}$$

$${\dot{r}}_{\mathrm{dyn}}$$

Full text

I Introduction

The last decade of gravitational-wave (GW) astronomy, driven by observations from the LIGO and Virgo detectors, has produced a catalog of $\mathcal{O}(100)$ compact binary coalescences (CBCs) [1]. This has inaugurated an era of gravitational-wave astrophysics, yielding unprecedented insights into the populations and properties of stellar remnants, enabling novel tests of General Relativity [2], new measurements of the cosmic expansion rate [3], and a deeper understanding of compact object populations [4]. Looking ahead, next-generation observatories such as the space-based Laser Interferometer Space Antenna (LISA) [5], and the planned third-generation ground-based detectors Einstein Telescope (ET) [6] and Cosmic Explorer (CE) [7], promise order-of-magnitude improvements in sensitivity. These instruments will regularly produce signals with signal-to-noise ratios (SNRs) in the hundreds or thousands, demanding GW waveform models of correspondingly higher fidelity. Indeed, even at current detector sensitivities, systematic errors in waveform modeling can appreciably bias inferred source parameters [8], and reducing such modeling errors will be absolutely critical for future high-SNR science.

The increasing sensitivity and anticipated event rates necessitate highly accurate and computationally efficient gravitational waveform models. These models are the cornerstone of GW astronomy, essential for detecting signals, accurately estimating source parameters, and performing robust astrophysical inference [9]. Current data analysis pipelines, particularly those based on Bayesian inference, often require millions of waveform evaluations for a single event, demanding models that are both precise and fast.

Several approaches model CBC waveforms. Numerical Relativity (NR) simulations provide the most accurate solutions to Einstein’s equations for the strong-field merger phase [10], serving as crucial benchmarks. Post-Newtonian (PN) theory offers analytical approximations for the early inspiral [11] but breaks down near merger. Phenomenological models [12] and surrogate models [13] combine analytical insights with fits to NR data. The Effective-One-Body (EOB) formalism [14] maps the two-body problem to an effective particle in a deformed metric, incorporating PN information, resummation techniques, and NR calibration. State-of-the-art EOB models like SEOBNRv5HM [15] (which includes higher-order modes and gravitational self-force information) and TEOBResumS [16] are widely used. While the ultimate goal includes precessing spins (e.g., via SEOBNRv5PHM development), accurate aligned-spin models are a crucial foundation.

This paper focuses on advancing aligned-spin EOB models. While SEOBNRv5HM [15] is highly accurate, its merger-ringdown and non-quasi-circular (NQC) corrections often rely on phenomenological fits to NR data (see Table 1). This empirical tuning can obscure underlying physics and limit extrapolation, especially to regions like extreme mass ratios or nearly extremal spins where NR calibration is sparse.

To address this, we combine EOB with the Backwards-One-Body (BOB) formalism [17]. BOB describes the late inspiral, merger, and ringdown by working backwards from the final remnant black hole, effectively treating the late-stage binary as a linear perturbation of the stationary remnant. It leverages the idea that outgoing GWs are driven by perturbations near the remnant’s light ring [18], with an underlying orbital frequency smoothly transitioning to the remnant’s quasinormal mode (QNM) frequencies [19, 20, 21]. Notably, the original BOB model reproduced NR merger-ringdown waveforms to within NR uncertainty without phenomenological tuning parameters, suggesting an accurate analytic description of merger physics from first principles.

Our primary contribution is the Spinning Effective-to-Backwards One Body (SEBOB) model. We present two distinct strategies for this integration. The first, termed seobnrv5_nrnqc_bob, uses the standard NR-calibrated Non-Quasi-Circular (NQC) corrections from SEOBNRv5HM for the inspiral-plunge phase, then smoothly attaches a BOB model to analytically model the plunge, merger, and ringdown. The second, more ambitious strategy, seobnrv5_bob, utilizes BOB not only for the merger-ringdown but also to inform the NQC corrections themselves, thereby reducing reliance on direct NR fits and facilitating, by construction, higher-order continuity (e.g., $\mathcal{C}^{2}$) at the transition point. SEBOB aims to be an analytically driven model for the most nonlinear stage, enhancing robustness where NR calibrations are unreliable. In comprehensive comparisons to SXS NR waveforms (focusing on the dominant $(\ell,m)=(2,2)$ mode), both SEBOB variants achieve median noise-free mismatches of $\mathcal{M}\simeq 2\times 10^{-4}$, competitive with SEOBNRv5 and TEOBResumS, even as the BOB-informed variant reduces direct NR calibration (see Sec. V).

We implement SEBOB within NRPy [22], an open-source Python package for symbolic manipulation and optimized C-code generation. NRPy translates EOB Hamiltonians, fluxes, and waveform quantities from a human-readable symbolic form into highly efficient, automatically generated C code, using advanced common subexpression elimination (CSE) on Hamiltonian, flux, and strain expressions. This offers the transparency and flexibility advantages of developing in Python (making it an ‘experimental playground’), while retaining the performance benefits of C. In wall-time benchmarks (Sec. V), our NRPy-based C implementation yields $\sim 3\times$ speedups over pySEOBNR [23] for equal and moderate mass ratios, while remaining competitive with other state-of-the-art generators, including lalsuite’s C-based SEOBNRv5_ROMs [24](when taking into account the ROM is already in frequency domain). Our analysis also identifies the main sources of residual error near the peak amplitude, pointing to improvements in the treatment of amplitude curvature that we defer to future work.

In summary, SEBOB addresses several needs for next-generation gravitational-wave modeling. It aims for improved accuracy in the late merger-ringdown through analytic insight, which reduces dependence on empirical calibration to NR. Furthermore, it is designed for computational efficiency and scalability suitable for the data rates and parameter-space demands of future observatories. Finally, it promotes enhanced reproducibility and transparency through its open-source implementation. The remainder of this paper is organized as follows: Sec. II details the SEOBNRv5HM formalism. Sec. III introduces BOB and our SEBOB integration prescriptions. Sec. IV discusses the NRPy implementation. Sec. V presents accuracy and performance comparisons. Sec. VI concludes with a discussion of the results and future directions.

Throughout this paper, we employ geometrized units ($c=G=M=1$). Time derivatives are denoted by an overdot (e.g., $\dot{X}\equiv dX/dt$). The symmetric mass ratio of a binary is given by $\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$, where $m_{1}$ and $m_{2}$ are the masses of the component black holes. The dimensionless spin angular momenta of the individual black holes are represented by $\chi_{i}\equiv\vec{S}_{i}/m_{i}^{2}$. Waveform angular frequency i.e, derivative of the waveform phase, is denoted by $\omega$ and orbital frequency i.e, derivatives of the azimuthal angle, is denoted by $\Omega$.

II Spinning EOB Waveforms

The original EOB formalism [14] canonically mapped the non-spinning, second-order PN two-body problem to the geodesic motion of an effective test particle in a modified Schwarzschild spacetime. Subsequent advances, including higher PN orders, spins, resummation techniques, and calibration against NR data, have yielded state-of-the-art models such as SEOBNRv5HM[15] and TEOBResumS[25]. While these models share foundational EOB principles, this work focuses specifically on SEOBNRv5HM. This section details the analytical foundations of SEOBNRv5HM, discusses the incorporation of NR data into the formalism, and outlines the NR parameters essential for constructing inspiral-merger-ringdown (IMR) waveforms.

Within the EOB formalism, the conservative Hamiltonian $\mathcal{H}_{\mathrm{real}}$ of the PN dynamics is canonically mapped to an effective Hamiltonian $\mathcal{H}_{\mathrm{eff}}$ via the relation:

[TABLE]

The effective Hamiltonian $\mathcal{H}_{\mathrm{eff}}$ describes the geodesic motion of an effective test particle in a modified Kerr metric $g_{\mu\nu}(x^{j})$, with canonical coordinates and momenta denoted by $\{x^{j},p_{j}\}$. It is derived from the relativistic mass-shell constraint, $g^{\mu\nu}p_{\mu}p_{\nu}=-1$, supplemented by a potential $\mathcal{Q}(x^{j},p_{j})$ that accounts for quartic ($\mathcal{O}(p_{j}^{4})$) and higher-order corrections to geodesic motion. The effective Hamiltonian is given by:

[TABLE]

Expressed in spherical polar coordinates $(r,\phi)$ with associated momenta $(p_{r},p_{\phi})$, the effective Hamiltonian also takes the form:

[TABLE]

where $H_{\mathrm{odd}}$ encodes spin-orbit interactions, $H_{\mathrm{even}}$ describes orbital motion, and $p_{\phi}$ represents the system’s canonical angular momentum. For a detailed expansion of the Hamiltonian, readers are referred to Section II and Appendix A of Ref. [15]. In the non-spinning limit ($a_{i}\rightarrow 0$, denoted by the subscript $\mathrm{noS}$), the orbital component $H_{\mathrm{even}}$ simplifies to:

[TABLE]

The functions $A_{\mathrm{noS}}$, $D_{\mathrm{noS}}$, and $Q_{\mathrm{noS}}$ incorporate the canonically mapped PN information in this limit. To improve computational efficiency and numerical stability closer to the effective horizon, the coordinates $(r,p_{r})$ are transformed to tortoise coordinates $(r_{*},p_{r_{*}})$ via:

[TABLE]

Inspiral gravitational waveform modes $h_{lm}$ are constructed using a factorized representation, $h^{\mathrm{F}}_{lm}$, derived from PN waveform modes. These modes depend on the binary parameters ($\eta,\chi_{1},\chi_{2}$), the canonical phase space variables $(\vec{r},\vec{p})$, the EOB Hamiltonian $\mathcal{H}_{\mathrm{real}}$, the instantaneous orbital frequency $\Omega=d\phi/dt$, and $\Omega_{\mathrm{circ}}$, which is the corresponding frequency for a circular orbit ($p_{r_{*}}\rightarrow 0$) at radius $r$. Gravitational radiation reaction is incorporated by resumming factorized modes $h^{\mathrm{F}}_{lm}$ into a dissipative flux term that modify Hamilton’s equations. The angular momentum flux $\mathcal{F}_{\phi}$ and the corresponding radial flux $\mathcal{F}_{r}$ are given by:

[TABLE]

The canonical coordinates and momenta are evolved numerically using Hamilton’s equations of motion, with the fluxes $\mathcal{F}_{\phi}$ and $\mathcal{F}_{r}$ included as external force terms:

[TABLE]

The equations of motion are typically integrated until the orbital frequency $\Omega$ or the time derivative of the radial momentum $dp_{r_{*}}/dt$ reaches an extremum. At each integration timestep, auxiliary quantities ($\mathcal{H}_{\mathrm{real}}$, $\Omega$, and $\Omega_{\mathrm{circ}}$) that are necessary for computing the factorized waveform modes $h^{\mathrm{F}}_{lm}$ are evaluated alongside the canonical variables.

II.1 Numerical relativity in the EOB formalism

The EOB formalism, as presented thus far, relies primarily on analytical PN approximations of general relativity. To enhance the accuracy of inspiral waveforms for gravitational-wave signal parameter estimation, corrections from NR are essential. This subsection details the three primary NR inputs used to construct the full IMR waveform: non-quasi-circular (NQC) corrections, the merger-ringdown (MR) waveform model, and inspiral calibration parameters.

First, NQC corrections are applied to improve waveform accuracy in the plunge regime, where the quasi-circular orbit approximation used for $h^{\mathrm{F}}_{lm}$ becomes invalid. After obtaining the EOB trajectory and factorized waveform, the NQC-corrected inspiral-plunge waveform $h^{\mathrm{insp-plunge}}_{22}$ is computed as:

[TABLE]

Coefficients $\{a_{1},a_{2},a_{3}\}$ and $\{b_{1},b_{2}\}$ are determined by matching $h^{\mathrm{insp-plunge}}_{22}$’s amplitude and frequency (defined as the phase derivative) to corresponding NR values at a specific time, $t_{22}^{\mathrm{peak}}$. This $t_{22}^{\mathrm{peak}}$ corresponds to the peak amplitude of the NR $(2,2)$ mode strain, related to EOB dynamics via the tunable parameter $\Delta t^{\mathrm{ISCO}}$:

[TABLE]

Here, $t^{\mathrm{ISCO}}$ denotes the time the EOB trajectory crosses the innermost stable circular orbit (ISCO) radius of the remnant black hole. The remnant’s mass $M_{f}$ and dimensionless spin $a_{f}$, essential for ISCO radius calculation, are obtained from NR fitting formulae [26, 27]. To determine the amplitude coefficients $\{a_{1},a_{2},a_{3}\}$, the factorized EOB waveform and trajectory are interpolated near $t_{22}^{\mathrm{peak}}$. These coefficients are then found by enforcing agreement between the NQC-corrected EOB amplitude $|h^{\mathrm{insp-plunge}}_{22}|$ (and its first two time derivatives) and the corresponding NR amplitude $|h_{22}^{\mathrm{NR}}|$ (and its first two time derivatives, with the first derivative being zero at peak amplitude), all evaluated at $t_{22}^{\mathrm{peak}}$:

[TABLE]

Similarly, phase correction coefficients $\{b_{1},b_{2}\}$ are determined by matching the NQC-corrected frequency $\omega_{22}=d\phi_{22}/dt$ and its time derivative $\dot{\omega}_{22}$ of the NQC-corrected EOB waveform to corresponding NR values at $t_{22}^{\mathrm{peak}}$:

[TABLE]

Given the trajectory and the factorized waveform, the quantities $\mathcal{Q}$ and $\mathcal{P}$ can be interpolated, spline-differentiated, and evaluated at $t_{22}^{\mathrm{peak}}$. Equations 14 and 16 are then solved numerically to determine the NQC coefficients $\{a_{1},a_{2},a_{3}\}$ and $\{b_{1},b_{2}\}$ respectively. Second, the EOB formalism incorporates a separate MR model to specify the full gravitational strain of a black hole binary merger. In SEOBNRv5, MR waveform modes are modeled as the emission of quasi-normal modes (QNMs)[19] from the remnant black hole, with amplitudes and phases modulated by phenomenological functions fitted to NR data:

[TABLE]

Here, $\omega_{22}^{\mathrm{QNM}}$ and $\tau_{22}^{\mathrm{QNM}}$ are the frequency and damping time of the dominant $(2,2)$ QNM. These are calculated as functions of the remnant mass $M_{f}$ and spin $a_{f}$ (obtained from NR fits) using tools such as the qnm package [28]. Coefficients $\{c_{1,f},c_{2,f},d_{1,f},d_{2,f}\}$ are determined by fitting to a catalog of NR waveforms. The remaining coefficients $\{c_{1,c},c_{2,c},d_{1,c},\phi_{0}\}$ are fixed by requiring $C^{1}$ continuity (continuous value and first derivative) when attaching the MR waveform to the NQC-corrected inspiral-plunge waveform at $t_{22}^{\mathrm{peak}}$.

Third, after NQC corrections and the MR model structure are incorporated, the overall accuracy of the resulting IMR waveform is further enhanced by tuning three parameters $\left\{a_{6},d_{\mathrm{SO}},\Delta t^{\mathrm{ISCO}}\right\}$ within the EOB dynamics against a large set of NR waveforms. The parameter $a_{6}$ is a pseudo-PN coefficient entering the non-spinning inspiral potential $A_{\mathrm{noS}}$ at 5PN order:

[TABLE]

where $u\equiv 1/r$ is the inverse of the radial separation. The parameter $d_{\mathrm{SO}}$ is a pseudo-PN parameter that enters the spin-orbit Hamiltonian at the 4.5PN order:

[TABLE]

where the spin combinations $a_{\pm}=\chi_{1}m_{1}\pm\chi_{2}m_{2}$, $g_{a_{\pm}}$ are the corresponding gyrogravitomagnetic factors, and $G^{\mathrm{align}}_{a^{3}}$ contains cubic-in-spin couplings (See Equations(26)-(30) of  [29]). The third $\Delta t^{\mathrm{ISCO}}$ parameter, defined earlier, encodes the time between ISCO crossing and the peak of the NR waveform. These three parameters are calibrated via Nested Bayesian inference, minimizing mismatch against a set of 446 SXS [10] (and one Einstein Toolkit) waveforms. The inference proceeds in stages: first, non-spinning waveforms constrain $a_{6}$ and the spin-independent components of $\Delta t^{\mathrm{ISCO}}$; subsequently, $d_{\mathrm{SO}}$ and the spin-dependent components of $\Delta t^{\mathrm{ISCO}}$ are constrained. Crucially, this calibration is performed on NQC-corrected, merger-attached waveforms. Consequently, NR-derived parameters for NQC corrections ($|h_{22}^{\mathrm{NR}}|^{\mathrm{peak}},\ddot{|h|}_{22}^{\mathrm{peak}},\omega_{22}^{\mathrm{peak}},\dot{\omega}_{22}^{\mathrm{peak}}$) and MR model fits ($c_{i,f},d_{i,f},M_{f},a_{f}$) are considered “pre-calibration” inputs, fixed prior to the final tuning of $\{a_{6},d_{\mathrm{SO}},\Delta t^{\mathrm{ISCO}}\}$. The complete set of NR-derived parameters for the SEOBNRv5 model is summarized in Table 1.

The SEOBNRv5 waveform model, which combines novel PN approximation approaches with advances in accurate NR simulations, provides highly accurate templates for analyzing gravitational wave signals from black hole binary mergers observed by LIGO, Virgo, and KAGRA. The model [15] is implemented in Cython within the open-source package pySEOBNR[23]. For faster waveform evaluation, a reduced-order surrogate model, SEOBNRv5_ROM, is also implemented in C as part of lalsuite[24].

III Backwards One Body Waveforms

The Backwards One Body (BOB) formalism [17] models the late inspiral and merger-ringdown stages of black hole binary coalescence. It describes these phases as dynamical perturbations propagating outwards through null congruences near the remnant’s light ring, superposed on the stationary spacetime of the merged remnant black hole. This section first reviews BOB’s mathematical formulation, then presents and justifies two distinct prescriptions for matching EOB inspiral waveforms. Throughout this section, $\omega(t)$ and related quantities like $\omega_{0}$ and $\omega_{\mathrm{QNM}}$ refer to the gravitational waveform frequency for a given $(l,m)$ mode.

The BOB formalism models gravitational wave perturbation amplitudes using the Newman-Penrose scalar [18] $\psi_{4}$, which represents the second derivative of the gravitational wave strain, $\ddot{h}$. Based on the divergence of null congruences near the light ring, the amplitude of each $(l,m)$ mode is given by (cf. Eq. 5 in Ref. [17]):

[TABLE]

where $(t_{p,lm},A_{p,lm})$ correspond to the time and value of the peak $\psi_{4,lm}$ amplitude for that mode, and $\tau_{lm}$ is the damping time of the $(l,m)$ quasinormal mode (QNM) [19] of the remnant black hole.

The waveform phase derivation assumes a single underlying orbital frequency $\Omega(t)$, which governs all multipole frequencies via $\omega_{lm}=m\Omega$. This frequency evolution is inferred from the identity [30, 31] relating the news function, $\mathcal{N}_{lm}\equiv\dot{h}_{lm}$, to the frequency evolution:

[TABLE]

Combined with the adiabatic relation $|\psi_{4,lm}|^{2}\approx m^{2}\Omega^{2}|\mathcal{N}_{lm}|^{2}$ (more precisely, $|\psi_{4,lm}|^{2}\propto m^{4}\Omega^{3}\frac{d\Omega}{dt}$ which leads to the form below) that assumes a negligible contribution from the news amplitude derivative to $\psi_{4}$, we obtain an expression for $\omega_{lm}=m\Omega$. Dropping mode labels, the frequency $\omega(t)$ is given by:

[TABLE]

with $k$ defined as:

[TABLE]

where $t_{0}$ and $\omega_{0}$ are the time and waveform frequency at a reference point, and $\omega_{\mathrm{QNM}}$ is the fundamental QNM frequency of the remnant black hole. Integrating $\omega(t)$ then yields the waveform phase $\phi(t)$:

[TABLE]

where $\phi_{0}$ is the phase at the reference time $t_{0}$, $\omega_{\pm\infty}$ are the frequencies evaluated in the limit $t\to\pm\infty$, and the terms $\mathrm{arctan}_{\pm}$/$\mathrm{arctanh}_{\pm}$ are defined as:

[TABLE]

where $\mathrm{tan}^{-1}$ and $\mathrm{tanh}^{-1}$ are the inverse tangent and inverse hyperbolic tangent respectively.

To relate $\psi_{4}$ to the gravitational wave strain $h$, we assume the phase evolves linearly at leading order near merger, yielding $\psi_{4}=\ddot{h}\approx-\omega(t)^{2}h$. Under this approximation, the complex strain is given by:

[TABLE]

where $|h_{0}|$ is the strain amplitude at the reference time $t_{0}$.

A condition relating $t_{p}$ to $t_{0}$ is required. For matching with SEOBNRv5, $t_{0}$ is chosen as the time of peak strain amplitude, $t^{\mathrm{peak}}_{22}$, the natural end point of the EOB inspiral. Requiring $d|h|/dt=0$ at $t=t_{0}$ leads to the following expression:

[TABLE]

where the identity

[TABLE]

is used to invert the hyperbolic tangent function. Equation (32) differs from the $t_{p}$ condition derived in Eq. 9 of Ref. [17], which imposes $\mathcal{C}^{2}$ phase continuity and depends on $\dot{\omega}_{0}$ for an arbitrary matching time $t_{0}$. However, within the EOB framework, non-quasi-circular (NQC) corrections ensure a peak in the strain at the matching time. This choice for $t_{0}$ therefore ensures the matching point corresponds to the peak strain amplitude in both the EOB and BOB models. This approach, implemented in the seobnrv5_nrnqc_bob model, utilizes BOB to model the merger-ringdown, attached to the NQC-corrected EOB inspiral.

Alternatively, BOB can directly provide NQC corrections. The conventional SEOBNRv5 NQC corrections [cf. Eq. 12] rely on NR fits to the peak amplitude, frequency, and their derivatives at the peak strain time. With the BOB approach, these derivatives are instead computed analytically from the BOB model (Eqs. 31 and 26) evaluated at the peak strain time. The functional form of the NQC-corrected inspiral-plunge waveform from Eq. (12) is maintained, but the coefficients $\{a_{1},a_{2},a_{3},b_{1},b_{2}\}$ are instead determined by requiring agreement with the BOB-predicted values for the strain amplitude $|h_{\mathrm{BOB}}|$, its first time derivative (which is zero at $t_{0}$ by construction using Eq. (32)), its second time derivative $d^{2}|h_{\mathrm{BOB}}|/dt^{2}$, as well as the frequency $\omega_{\mathrm{BOB}}$ and its first time derivative $d\omega_{\mathrm{BOB}}/dt|$, all evaluated at $t_{0}$. This model is implemented in NRPy as seobnrv5_bob.

The approach of using BOB to inform NQCs is motivated by the observation that the peak strain occurs close to the perturber crossing the remnant innermost stable circular orbit (ISCO), a regime particularly apt for BOB’s description. It reduces reliance on direct NR fitting for NQC parameters by leveraging BOB’s analytical structure. Primary NR inputs are thus limited to the remnant properties (mass and spin, which determine $\omega_{\mathrm{QNM}}$ and $\tau$) and the strain amplitude $|h_{0}|$ and frequency $\omega_{0}$ at the matching time $t_{0}$. Furthermore, while the SEOBNRv5 merger-ringdown attachment guarantees only $\mathcal{C}^{1}$ continuity (continuous complex $h$ and $\dot{h}$) at the matching point $t_{0}$, employing BOB for NQC computation inherently links inspiral corrections to the merger-ringdown model, thus enabling $\mathcal{C}^{2}$ continuity (continuous $h,\dot{h},\ddot{h}$) by construction.

Conversely, NQC corrections are applied to the entire inspiral waveform via their dependence on EOB dynamical variables ($p_{r_{*}},r,\Omega$). For generic non-circular orbits, the physically correct dependence may involve post-Newtonian (PN) waveform predictions. In such cases, BOB may not provide a complete description of the inspiral’s NQC corrections. Therefore, employing BOB solely as a merger-ringdown attachment (seobnrv5_nrnqc_bob) retains the NR-calibrated NQCs for the inspiral and could be advantageous.

Equation (32) implies that the merger-ringdown waveform is determined solely by the strain amplitude $|h_{0}|$ and frequency $\omega_{0}$ at the matching time $t_{0}$, along with the remnant mass and spin used to calculate $\omega_{\mathrm{QNM}}$ and $\tau$. Thus, BOB waveforms eliminate all NR-informed degrees of freedom in modeling the merger-ringdown evolution, provided that the remnant mass and spin themselves are predicted from inspiral parameters without direct NR fitting (e.g., using methods like those in Ref. [32]). This can also aid in calibrating the EOB Hamiltonian to accurately predict remnant mass and spin directly from the inspiral dynamics (see [32]). Although the BOB formalism can model any gravitational wave mode $(l,m)$, our analysis in this paper is restricted to the dominant $(2,2)$ mode. The combination of BOB’s analytical model for the merger-ringdown with the EOB inspiral’s analytical foundations into an Effective-to-Backwards-One-Body (EBOB) framework aims to describe the complete black hole binary coalescence waveform analytically and from physical motivations. In Section V, we will assess the fidelity of BOB predictions against numerical relativity data and examine the impact of using BOB-derived NQCs on overall waveform accuracy.

IV Numerical Implementation

The SEBOB waveform model, comprising the seobnrv5_nrnqc_bob and seobnrv5_bob variants, is implemented in the open-source NRPy framework [22, 33]. NRPy, built upon the SymPy library [34] for symbolic manipulation, translates all required EOB and BOB mathematical expressions into highly optimized, modular C code. This symbolic-to-numerical workflow, which extensively applies Common Subexpression Elimination (CSE) to reduce runtime costs, facilitates transparency, extensibility, and high performance.

The generated C code is self-contained and can be compiled as an executable or linked as a library. It relies on the GNU Scientific Library (GSL) [35] for core numerical tasks, including root-finding, ordinary differential equation (ODE) integration, interpolation, and linear algebra. This structure enables seamless integration into gravitational-wave data analysis packages such as LALSuite [24] by avoiding Python-to-C overhead, and supports deployment within infrastructures such as BlackHoles@Home [33] for large-scale model calibration and waveform generation.

IV.1 EOB Dynamics: Initial Conditions and Evolution

IV.1.1 Initial Conditions

EOB evolution begins by determining initial conditions for the canonical phase-space variables $\{r,\phi,p_{r_{*}},p_{\phi}\}$. First, for a user-specified initial orbital frequency $\Omega_{0}$, the initial orbital separation $r_{0}$ and azimuthal momentum $p_{\phi,0}$ for quasi-circular orbits are obtained by solving the two-dimensional non-linear algebraic system:

[TABLE]

This system is solved using GSL’s hybrid Powell multi-dimensional root-finding algorithm. The first condition is equivalent to no radial motion since $p_{r_{*}}=0$ guarantees $dr/dt=0$ since the Hamiltonian is quadratic in $p_{r_{*}}$.

Second, the initial radial momentum $p_{r_{*,0}}$, typically small but non-zero for dissipative inspirals, is determined by solving a one-dimensional root-finding problem using GSL’s Brent-Dekker solver. The radial momentum is derived by equating the initial conservative radial velocity $\dot{r}_{\mathrm{dyn}}=\xi\partial\mathcal{H}_{\mathrm{real}}/\partial p_{r_{*}}$ to the radial velocity due to radiation reaction $\dot{r}_{\mathrm{rad}}=\dot{r}_{\mathrm{dyn}}$ in the adiabatic limit, with $\dot{r}_{\mathrm{rad}}$ given by [32]

[TABLE]

The initial phase $\phi_{0}$ is set to zero.

IV.1.2 Orbital Evolution

With initial conditions established, Hamilton’s equations of motion are integrated forward in time using a GSL adaptive 8th-order Runge-Kutta Prince-Dormand method (or similar) with error control. Tortoise coordinates $(r_{*},p_{r_{*}})$ are used to enhance numerical stability, particularly near the EOB horizon analogue. Integration terminates after a radial separation of $6M$ is reached and one of the following conditions are met:

• •

The orbital frequency $\Omega$ reaches a maximum and begins to decrease, indicating an unphysical peak in frequency evolution.

• •

$dp_{r_{*}}/dt$ becomes positive or peaks, indicating an unphysical outspiral.

• •

$dr/dt$ becomes positive, indicating an unphysical outspiral.

• •

The orbital separation $r$ drops below a predefined minimum (e.g., related to $r_{\mathrm{ISCO}}$). This is only realized in cases where the parameter $\Delta t^{\mathrm{ISCO}}$ is positive, indicating the NR waveform strain peaks before the EOB trajectory reaches the ISCO of the remnant black hole.

• •

The circular orbit frequency $\Omega_{\mathrm{circ}}$ exceeds unity at radius $r<3M$. This is indicative of unphysically high rotation speeds and is used to terminate the integration.

Throughout the integration, quantities required for radiation reaction fluxes ($\mathcal{F}_{r},\mathcal{F}_{\phi}$) and waveform modes ($h^{\mathrm{F}}_{lm}$), such as $\mathcal{H}_{\mathrm{real}}$, $\Omega$, and $\Omega_{\mathrm{circ}}$, are computed from the phase-space state and binary parameters. These quantities are stored at variable resolution, with coarser sampling during the early, slower inspiral and finer sampling of the last $250M$ in time for performing NQC corrections.

IV.1.3 Factorized Waveform Generation

Factorized inspiral waveform modes $h^{\mathrm{F}}_{lm}$ are computed at each step of the finely sampled EOB trajectory. These modes are constructed symbolically within NRPy and generated into optimized C code. This symbolic generation significantly facilitates faster computation of waveform modes, particularly by enabling pre-computation of coefficients that depend only on constant binary parameters. This benefit is particularly important for precessing binaries, where these coefficients become time-dependent due to spin vector evolution and thus cannot be precomputed during integration.

IV.2 Non-Quasi-Circular Corrections and Merger Attachment

To accurately model the late inspiral and plunge, Non-Quasi-Circular (NQC) corrections are applied to the dominant $(l=2,m=2)$ factorized EOB mode $h^{\mathrm{F}}_{22}$, as described by Eq. (12). This involves determining the NQC amplitude coefficients $\{a_{1},a_{2},a_{3}\}$ and phase coefficients $\{b_{1},b_{2}\}$ by solving $3\times 3$ and $2\times 2$ linear systems (Eqs. (14) and (16)), respectively, using GSL routines for LU decomposition and forward/backward substitution.

These systems require quantities computed at the predicted matching time, $t_{0}$, defined as:

[TABLE]

where $t^{\mathrm{ISCO}}$ is the time the EOB trajectory crosses the innermost stable circular orbit (ISCO) radius of the remnant black hole, and $\Delta t^{\mathrm{ISCO}}$ is a calibration parameter that sets the time at which the NR waveform strain peaks relative to the ISCO time. Both $t^{\mathrm{ISCO}}$ and $t_{0}$ are determined from the finely sampled portion of the trajectory by identifying the nearest index in the radial or time arrays, followed by interpolation between the 10 nearest sample points at higher resolution.

In a minority of cases, if ODE integration terminates before the perturber crosses the remnant ISCO or before $t_{0}$, then $t_{0}$ is set to the last sample point of the trajectory. Consistency in NQC coefficient calculation is crucial, as GSL cubic spline interpolation (used for time derivatives at the last sample point) applies natural boundary conditions, while pySEOBNR (also used) employs not-a-knot boundary conditions, leading to significant differences. To ensure consistency, we modified pySEOBNR to use natural boundary conditions for all its cubic spline interpolations.

IV.2.1 seobnrv5_nrnqc_bob: NR-informed NQCs and BOB Merger-Ringdown

This variant mirrors the SEOBNRv5HM approach for NQC corrections. The target values for the right-hand side of Eqs. (14) and (16)—specifically, the amplitude $|h_{22}^{\mathrm{NR}}|$, its first and second time derivatives, the frequency $\omega^{\mathrm{NR}}_{22}$, and its first time derivative—are derived from fits to NR simulations, evaluated at the peak of the (2,2) mode of the strain. After constructing the NQC-corrected inspiral-plunge waveform $h^{\mathrm{insp-plunge}}_{22}$, the BOB model (Sec. III) is attached at $t_{0}$ to describe the merger and ringdown. This attachment procedure ensures $\mathcal{C}^{1}$ continuity of the complex strain $h$ and its first time derivative $\dot{h}$ at $t_{0}$.

IV.2.2 seobnrv5_bob: BOB-informed NQCs and Merger-Ringdown

This model integrates BOB more deeply into the inspiral by leveraging it to inform the NQC corrections. The functional form of the NQC corrections (Eq. (12)) is maintained, but the target values for the linear systems are sourced analytically from the BOB model itself (Eqs. (31), (26), and (32)), evaluated at $t_{0}=t_{22}^{\mathrm{peak}}$. Specifically, the required amplitude and frequency derivatives at the matching time, ($d^{2}|h_{\mathrm{BOB}}|/dt^{2}|_{t_{0}}$ and $d\omega_{\mathrm{BOB}}/dt|_{t_{0}}$) are computed analytically from the BOB equations, with NRPy performing symbolic differentiation and generating C code for their evaluation. Remnant parameters $M_{f}$ and $a_{f}$ (essential for computing $\omega_{\mathrm{QNM}}$ and $\tau$ within the BOB model) are obtained from NR fits [26, 27]. The QNM frequency $\omega_{\mathrm{QNM}}$ and damping time $\tau$ are computed by interpolating over a set of 107 values of the final spin and the corresponding QNM values computed in Python using the qnm package [28].

Once NQC coefficients are found using either BOB-derived or NR-informed targets, $h^{\mathrm{insp-plunge}}_{22}$ is constructed by multiplying the factorized waveform $h^{\mathrm{F}}_{22}$ by the NQC correction given in Eq. (12). Since the factorized waveform is computed at the sample points of the ODE trajectory, no additional interpolation is required. For $t\geq t_{0}$, the waveform is subsequently described by the same BOB model used for NQC targets. This construction, leveraging BOB-derived NQC targets and the BOB model for merger-ringdown, inherently ensures $\mathcal{C}^{2}$ continuity of $h,\dot{h},\text{and }\ddot{h}$ at $t_{0}$.

In order to perform consistency checks against the SEOBNRv5HM model, we also implement the native SEOBNRv5HM merger ringdown model given in Eq. (18). This model is implemented in NRPy as seobnrv5_nrpy.

The SEBOB waveform code used in this paper can be generated with the following commands (flag -h provides instructions to generate specific variants):

1  pip install git+https://github.com/nrpy/nrpy.git@e16a1192  python3 -m nrpy.examples.seobnrv5_aligned_spin_inspiral -h The C code can be built and run with the following commands (flag -h provides instructions to pass input parameters):

1  cd project/seobnrv5_aligned_spin_inspiral2  make3  ./seobnrv5_aligned_spin_inspiral -hTo access up-to-date versions of developmental features discussed in this paper, we recommend directly cloning the repository – git clone git@github.com:nrpy/nrpy.git.

V Results

We present a systematic evaluation of the SEBOB waveform model. After verifying numerical consistency between our NRPy implementation and the public pySEOBNR code, we compare SEBOB and other state-of-the-art effective-one-body (EOB) models to all 607 quasi-circular, aligned-spin black hole binary, numerical-relativity (NR) waveforms from the SXS catalog (v2025.0.10 release [36]). We then test the key assumptions underlying the BOB merger–ringdown approximation and, finally, benchmark the computational efficiency of our implementation.

Throughout, waveform accuracy is quantified by the noise-free mismatch of the dominant $(2,2)$ strain mode,

[TABLE]

with inner product

[TABLE]

where $\tilde{h}(f)$ is the Fourier transform of the real part of the time-domain strain $h(t)$. The mismatches are calculated using the open source python package PyART(https://github.com/RoxGamba/PyART).

V.1 Numerical Consistency

A known implementation difference between NRPy and pySEOBNR is the choice of cubic-spline boundary conditions, which in turn affects non-quasi-circular (NQC) coefficients when the peak strain time is not attained by the trajectory. To isolate purely numerical effects we consider two modified pySEOBNR variants:

pyseobnr_pert, which applies random $10^{-15}$-level perturbations to the mass ratio and spins, and 2. 2.

pyseobnr_nat, which replaces all interpolations with natural splines 111Natural splines set the second derivative to zero at each endpoint, whereas the public code employs “not-a-knot” conditions..

Figure 1 shows that altering the interpolation boundary conditions induces slightly larger mismatches than round-off perturbations.

By contrast, Figure 2 demonstrates round-off-level agreement between NRPy and pyseobnr_nat, confirming that our code faithfully reproduces the core SEOBNRv5 dynamics.

In Table 2, we state the typical frequency of cases when the estimated peak time is not encountered. We profile stop conditions for a random combination of 5000 mass ratios, spins, and start frequencies to estimate how often this can be encountered during parameter estimation. We note that NQC conditions evaluated at the endpoint of the trajectory in 10% of the random sample and 3% of SXS simulations.

V.2 Accuracy of SEBOB

Figure 3 compares mismatches against the SXS catalog for all publicly available EOB models. Both seobnrv5_bob and seobnrv5_nrnqc_bob are competitive with existing models, with a small loss in median accuracy due to the replacement of NR-fitted information in the merger-ringdown regime with BOB.

To highlight how BOB-modelled merger-ringdown and NQC corrections manifest in the waveform, we show the waveform strain amplitude (top) and frequency (bottom) of the $(2,2)$ mode for the the case where BOB is least accurate against SXS compared to SEOBNRv5 (Figure 4) and the case where BOB is most accurate against SXS compared to SEOBNRv5 (Figure 5). Comparing both cases, we find that BOB provides a more accurate description of the frequency evolution compared to the amplitude, even when used as a replacement for NR-fitted NQC corrections. Particularly, in Figure 5, we note that seobnrv5_bob poorly models the curvature (second derivative) of the amplitude close to the peak. Since BOB uses the assumption of adiabaticity ($|h|\approx|\psi_{4}|\omega^{2}$, i.e $|\psi_{4}|$ is slowly varying), the variation of the strain amplitude is not described accurately if the variation in the $\psi_{4}$ amplitude is non-negligible. In the subsequent discussion, we will assess the validity and impact of the BOB postulates and assumptions.

V.3 Validity of BOB

The BOB derivation assumes (i) a common phase $\phi$ for $\psi_{4}$, the news $\mathcal{N}$, and the strain $h$,

[TABLE]

and (ii) an adiabatically varying amplitude,

[TABLE]

where $\omega\equiv\dot{\phi}$. To quantify departures from these postulates we compute four root-mean-squared errors (RMSEs) over the interval $[t_{0},t_{0}+20M]$ following peak strain $t_{0}$:

[TABLE]

Figure 6 shows that BOB’s core relations hold to good accuracy for quasi-circular, aligned-spin binaries; the largest mismatches coincide with elevated $\epsilon_{1}$ and $\epsilon_{2}$. Improvements will therefore focus on the news–frequency linkage and $\psi_{4}$ amplitude model [38].

V.4 NQC Derivatives from BOB

Unlike standard SEOBNRv5, which uses NR-fitted derivatives at peak strain, seobnrv5_bob computes the second derivative of the amplitude and the first derivative of the frequency analytically from BOB. Their accuracy is assessed in Figure 7 222The highest error case in the amplitude and frequency fits corresponds to SXS:BBH:1110 ($q=7,\chi_{1}=3\times 10^{-6},\chi_{2}=-2\times 10^{-7}$), which was deprecated during final development of the paper.. The quality of the NR fits and the BOB-derived quantities are measured by relative error,

[TABLE]

where $y_{\mathrm{NR}}$ is the value of the quantity obtained at the time of the peak strain directly from waveforms in the SXS catalog and $y_{\mathrm{pred}}$ is the predicted quantity either from the fits given in SEOBNRv5 or derived from BOB.

The frequency derivative is reproduced to $\sim 4\%$, close to the $\sim 1\%$ NR-fit uncertainty, whereas the amplitude’s second derivative shows $\mathcal{O}(1)$ relative error. This deficiency stems from the assumption of a slowly varying $\psi_{4}$ amplitude. However, even these sizeable parameter errors translate into only mild increases in mismatch; this is evidence that BOB captures the dominant merger-ringdown physics despite its simplified assumptions. We note that these errors can be decreased significantly by formulating BOB to describe the news, rather than $\psi_{4}$, and by more accurately calculating the amplitude without assuming adiabaticity, as will be shown in an upcoming paper [38]. Additionally, we note that the single outlying case in the amplitude and frequency fits, as seen in the bottom panels of Figure 6 and three of the panels in Figure 7, corresponds to SXS:BBH:1110 ($q=7,\chi_{1}=3\times 10^{-6},\chi_{2}=-2\times 10^{-7}$), where the SXS waveform has particularly noisy behavior in the vicinity of merger, as demonstrated in Figure 8 [39].

V.5 Efficiency

Because NRPy generates highly optimised C code, we anticipate that our implementation leads to faster waveform generation than the native pySEOBNR implementation. Since pySEOBNR only implements a CSE optimized Hamiltonian, the equivalently optimized flux terms in our implementation should significantly reduce computational time in the ODE right-hand-sides. Figure 9 shows that seobnrv5_bob is $\sim 3$ times faster than pySEOBNR for $q=[1,3]$ cases, and is competitive with other state-of-the-art waveform models. The cases shown in Figure 9 are a representative comparison of walltimes, since we found that changing the spin configurations did not affect the relative performance of the waveform models. All timings exclude I/O and post-processing overhead. seobnrv4_opt [40], NRHybSur3dq8 [13] and SEOBNRv5_ROM are timed using their LALSuite; TEOBResumS-GIOTTO uses its standalone C code without post-adiabatic speed-ups. One-time Cython compilation costs for pySEOBNR are omitted. In all cases, the walltime is computed by averaging over the walltime of 100 waveform generations. The computations were done using the Thorny Flat cluster at West Virginia University equipped with the Intel® Xeon® Gold 6138 Processor.

VI Conclusion

We have introduced SEBOB, a hybrid waveform framework that combines an Effective-One-Body (EOB) inspiral with a Backwards-One-Body (BOB) description of the merger–ringdown for aligned-spin binary black holes. Two realizations were constructed atop SEOBNRv5HM: (i) seobnrv5_nrnqc_bob, which retains NR-informed non–quasi-circular (NQC) corrections and attaches a BOB-based merger–ringdown, and (ii) seobnrv5_bob, which uses BOB to supply the NQC targets themselves, reducing reliance on phenomenological fits and enabling $\mathcal{C}^{2}$ continuity at the attachment time. Both are implemented in the open-source NRPy ecosystem, combining a user-friendly, transparent Python environment and a modularized, efficient C code.

Across a large set of SXS simulations for quasi-circular, aligned-spin binaries, both realizations attain median mismatches (with a flat psd) of $\mathcal{M}\simeq 2\times 10^{-4}$ for the dominant $(2,2)$ mode (Fig. 3), demonstrating that an analytically motivated merger–ringdown can maintain state-of-the-art fidelity. Diagnostics identify a localized imperfection in the curvature of the strain amplitude near its peak when NQC corrections are informed by the current $\psi_{4}$-based version of BOB; by contrast, the frequency evolution is captured to within a few percent, with typical peak $|\dot{\omega}|$ errors of $\sim 4\%$ (Fig. 7). Catalog-wide tests indicate that the core BOB postulates are satisfied to useful accuracy, with only weak correlations between their measured deviations and the overall mismatch (Fig. 6). Even with the localized amplitude shortfall, total mismatches remain at the few-$10^{-4}$ level.

A key practical outcome is computational performance. With symbolic optimization and comprehensive common–subexpression elimination applied to both the Hamiltonian and radiation–reaction fluxes, the NRPy-generated C implementation achieves $\sim 3\times$ shorter wall times than pySEOBNR for equal and moderate mass ratios while remaining competitive with other time-domain generators (Fig. 9). Since we do not benchmark the time taken to output and process waveforms, a caveat to Fig. 9 is that we do not take into account the time taken to produce a template bank or to perform parameter inference as part of our walltime analysis.

A key limitation of the BOB-informed NQC corrections is a lack of overlap between the EOB and BOB physics in the window where the corrections are necessary. EOB uses PN factorized waveform and fluxes, assumes circularity, and requires late-inspiral corrections to account for highly non-circular orbital motion and the general breakdown of PN close to merger. Upcoming work will explore constructing BOB directly from the news function without the assumption of adiabaticity, which should significantly reduce errors in the amplitude derivatives close to the peak [38]. Within the current paradigm of NQC corrections, this version of BOB can mitigate the need for higher derivatives of NR data to accurately model the amplitude corrections in the late inspiral.

A long-term goal for the SEBOB formalism is to eliminate the remaining direct reliance on NR fittings entering the EOB inspiral and dynamics. EOB models rely on extensive calibration of inspiral waveforms to numerical relativity (NR) waveforms to enhance accuracy. In SEOBNRv5, the underlying waveforms already incorporate NR-derived through NQC corrections and the phenomenological merger-ringdown (see Table 1 for a breakdown of pre-calibration NR parameters). A key advantage of SEBOB is that it limits pre-calibration NR information to 4 parameters ($|h|_{0},\omega_{0},M_{f},a_{f}$), significantly reducing the need for NR to model the merger-ringdown within calibration accuracy. It also reduces the need for extensive calibrations by providing a substantially accurate, physically motivated map from inspiral to merger-ringdown physics. Furthermore, self-consistently modeling non-circularities into the EOB fluxes and inferring remnant properties directly from the EOB inspiral (as suggested in [41]) can potentially reduce the number of NR parameters needed in the inspiral sector. Such improvements would significantly limit the need for NR calibration to a narrow window (between the EOB plunge and the BOB merger) and fewer parameters.

The present study is limited to spin-aligned, quasi-circular binaries and to the $(2,2)$ mode, and it continues to use NR fits for remnant properties. Several targeted extensions are natural: (i) improving the amplitude model by formulating BOB in terms of the news rather than $\psi_{4}$ and relaxing adiabaticity to sharpen NQC derivatives [38]; (ii) incorporating higher harmonics and, ultimately, generic precessing spins; (iii) tightening the inspiral-to-remnant mapping to further curtail NR-derived inputs; and (iv) developing an open-source pipeline to perform targeted calibrations to NR. By substantially shrinking the role of phenomenological calibration in the most nonlinear phase while retaining competitive accuracy and speed, SEBOB offers a robust, interpretable, and extensible path toward waveform modeling in the high-SNR era.

Acknowledgements.

We would like to acknowledge Anuj Kankani, Matthew Cerep, Lorenzo Pompili, Rossella Gamba, and Aasim Jan for useful comments and discussion. Scientific libraries numpy [42] and scipy [43] were used for part of the accuracy analysis. SM and STM were supported in part by National Science Foundation CAREER grant 1945130 and NASA grants 22-LPS22-0022 and 24-2024EPSCoR-0010. We acknowledge the computational resources provided by the WVU Research Computing Spruce Knob HPC cluster, which is funded in part by NSF grant EPS-1003907, and the Thorny Flat HPC cluster, which is funded in part by NSF grant OAC-1726534. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452.