1 Introduction↩︎

Since the first detection of gravitational waves from a coalescing binary black hole (BBH)[1], the LIGO-Virgo-KAGRA collaboration has confirmed over 390 gravitational wave (GW) signals from merging compact objects such as stellar mass black holes and neutron stars [2]–[7]. Waveform modeling has become a central pillar of GW astronomy: accurate waveform models are needed not only for detection, but also for precise inference of source properties and for tests of general relativity. Detector sensitivity will continue to improve with future ground-based detectors like the Einstein Telescope [8] and Cosmic Explorer [9] as well as space-based detectors such as the Laser Interferometer Space Antenna (LISA) [10], [11], TianQin [12], and Taiji [13]. As such, previously subdominant waveform systematics will become increasingly important [14], [15]. This is particularly true for long-lived signals which accumulate many inspiral cycles in band.

Including orbital eccentricity in waveform models has become increasingly important as evidence for eccentricity in detected signals has mounted, though not to the point of a definitive detection of eccentricity [16]–[26]. While gravitational radiation tends to circularize binaries during their inspiral [27]–[30], a range of astrophysical formation channels in globular clusters and active galactic nuclei can produce systems that retain measurable eccentricity when they enter the observational band of current and future detectors [31]–[36]. This is especially relevant for next-generation observatories which will be more sensitive to the inspiral where eccentricity effects are most pronounced [11].

As such, many waveform modeling frameworks have been extended to include orbital eccentricity [37]–[50]. In this work, we focus on the effective-one-body (EOB) framework [51]–[55]. It provides a powerful description of compact-binary dynamics by resumming information from post-Newtonian (PN) [56]–[62], gravitational self-force (GSF) [63]–[65], and recently post-Minkowskian (PM) theory [66]–[71], into an effective system of a test mass moving in a perturbed black hole spacetime. The model is calibrated against numerical relativity (NR) simulations to ensure agreement in the late inspiral and merger-ringdown [72]–[74]. The EOB paradigm has proven highly accurate and forms the basis of several state-of-the-art waveform models used in GW data analysis [74]–[86] with recent extensions including orbital eccentricity [87]–[105] and spin precession [106]–[108].

This accuracy comes at a cost: the inspiral dynamics must be obtained by numerically solving the EOB Hamiltonian equations of motion. In the quasicircular case, this can be bypassed with the “post-adiabatic approximation” [109]: one perturbs around a quasicircular orbit in order to recast the EOB equations of motion into a set of algebraic equations that can be solved analytically [110]–[112]. This approach is ill-suited to eccentric systems, so one must integrate the inspiral dynamics numerically. This is particularly costly for low-mass, large-mass-ratio systems whose inspirals contribute many waveform cycles within the detector band.

To address this problem, we adapt techniques that have proven highly successful in the modeling of eccentric extreme-mass-ratio inspirals (EMRIs) [113]–[116]. We first recast the eccentric EOB dynamics in terms of osculating orbital elements (OOE) [117], [118] (also known as the method of variation of constants [119], [120]) using a quasi-Keplerian parametrization à la Ref. [121]. This makes explicit the separation between the slowly evolving orbital elements and the rapidly varying orbital phase. We then apply near-identity (averaging) transformations (NITs) [122] to recast the short-timescale oscillatory structure from the inspiral dynamics as higher-order secular information needed for faithful waveform generation [123]–[131]. The result is a new set of averaged equations of motion that can be solved on the radiation-reaction timescale rather than the orbital timescale. Unlike in the EMRI case, we do not expand in the mass ratio, but rather in the small parameter \(\epsilon\) that characterizes the separation of timescales between the orbital and radiation-reaction timescales. This is suppressed not only by \(\nu\) but also by the fact that radiation reaction begins at 2.5PN order [27], [28]. We also extend the calculation beyond the next-to-leading-order terms needed for EMRIs [132], [133], to include next-to-next-to-leading-order terms. The resulting averaged dynamics are then mapped back to the full EOB dynamics for the final transition to plunge.

We validate the method by testing it on an inspiral-only, nonspinning EOB model. We compare waveforms generated from the averaged dynamics against those from the full EOB evolution across a range of eccentric configurations. We find that next-to-next-to-leading-order is needed to model comparable-mass systems accurately with the resulting model reproducing the full dynamics to mismatches \(\leq 8.05 \times 10^{-5}\) while substantially reducing computational cost by a factor of \(1.5 - 8 \times\). NITs thus provide a practical route to accelerating eccentric EOB waveform generation without sacrificing the accuracy required for data-analysis applications.

The remainder of this paper is organized as follows. In Sec. 2, we review the nonspinning eccentric EOB inspiral dynamics used in this work. In Sec. 3, we rewrite these dynamics using osculating orbital elements. In Sec. 4, we derive the near-identity averaged equations of motion through second post-adiabatic order. In Sec. 5, we describe waveform generation from both the full and averaged dynamics, and in Sec. 6 we summarize the numerical implementation. We present accuracy and runtime comparisons in Sec. 7, and conclude in Sec. 8 with a discussion of limitations and future extensions. In this work, we use geometric units where \(G = c = 1\).

2 Nonspinning effective one body dynamics↩︎

We now describe a two-body system consisting of a primary of mass \(m_1\) and a secondary of mass \(m_2\). For convenience, we define the total mass \(M = m_1 + m_2\) and the reduced mass \(\mu = m_1 m_2 / M\). In the effective one-body framework, we model the binary as a test mass moving in a perturbed black hole spacetime where the perturbation parameter is the symmetric mass ratio \(\nu = \mu / M = m_1 m_2 / (m_1 + m_2)^2\). We focus on the nonspinning sector, but our methods extend straightforwardly to (anti-)aligned spin effects.

The metric of this perturbed spacetime in EOB coordinates \((t, r,\phi)\) is given by: \[ds^2 = -A(\nu;r) dt^2 + \frac{dr^2}{A(\nu;r) \bar{D}(\nu;r)} + r^2 d\phi^2\] Note that in the test mass limit \(\nu \rightarrow 0\), \(A \rightarrow 1 - 2 M/r\) and \(\bar{D} \rightarrow 1\), recovering the Schwarzschild metric. Many choices for \(A\) and \(\bar{D}\) exist in the literature, but here we adopt those of SEOBNRv5 [83], [84].

Note that the test mass does not follow a geodesic in the perturbed spacetime as it is slowly inspiraling due to radiation reaction. Instead, the test mass follows a perturbed orbit satisfying: \[u^\alpha p_\alpha + Q(\nu;r,p_{r_*}) = g^{\beta \alpha} p_\beta p_\alpha + Q(\nu;r,p_{r_*}) = -1\] Here, \(u^\alpha\) and \(p_\alpha\) represent the four-velocity and four-momentum, respectively, of the test mass (all normalized by a factor of \(\mu\)), \(p_{r_*} = \xi(r) p_r\) is the tortoise radial momentum, where \(\xi(r) = A \sqrt{\bar{D}}\). This \(Q\) potential captures the effect of the slow inspiral on the conservative dynamics. As such, \(Q \rightarrow 0\) in the test mass limit. It is given as a PN series that appears at 3PN order and is expressed as an expansion in \(p_{r_*}\) (or \(p_r\)) and has the form \[\label{eq:Q95Pot} Q = Q_4(\nu,r) p_{r_*}^4 + Q_6(\nu,r) p_{r_*}^6 + Q_8(\nu,r) p_{r_*}^8 + \mathcal{O}(p_{r_*}^{10}).\tag{1}\] Again, we use the same expressions as SEOBNRv5 where the full expressions for \(Q_4,Q_6,\) and \(Q_8\) for nonspinning binaries can be found in appendix A of Ref. [83].

From the above equation, one derives the effective Hamiltonian by noting that \(p_t = - H_{\text{eff}}\), converting from \(p_r\) to \(p_{r_*}\) and rearranging: \[\label{eq:Heff} H_{\text{eff}}^2 = p_{r_*}^2 + A(\nu,r)\left[1 + \frac{p_\phi^{2}}{r^2} + Q(\nu;r,p_{r_*})\right]\tag{2}\]

The effective Hamiltonian is related to the real (EOB) Hamiltonian by the EOB energy map: \[H_\text{EOB} = M \sqrt{1 + 2 \nu (H_\text{eff} -1)}.\] From this Hamiltonian, one obtains the equations of motion via Hamilton’s equations: \[\begin{align} \label{eq:EOB95EOM} \begin{aligned} \frac{d r}{dt} &= \xi(r) \frac{\partial H_{\text{EOB} }}{\partial p_{r_*}}, \quad \frac{d \phi}{dt} = \frac{\partial H_{\text{EOB}}}{\partial p_\phi} = \Omega_\phi, \end{aligned}\\ \begin{align} \frac{d p_{r_*}}{dt} &= - \xi(r) \frac{\partial H_{\text{EOB}}}{\partial r} + \xi(r) F_r, \quad \frac{d p_\phi}{dt} = F_\phi, \end{align} \end{align}\tag{3}\] where \(F_r\) and \(F_\phi\) are the radiation-reaction force components. We use the nonspinning limit of the 3PN radiation-reaction force components from Ref. [62] expressed in EOB coordinates \((r,p_r,p_\phi)\); we provide them in Appendix 11.

3 Osculating Orbital Elements↩︎

3.1 Recasting into quasi-Keplerian parametrization↩︎

We now recast these equations into a set of osculating orbital element equations using the following quasi-Keplerian parametrization: \[\label{eq:r95of95p95e95xi} r = \frac{p M }{1 + e \cos \zeta},\tag{4}\] where \(p\) is the semilatus rectum, \(e\) is the eccentricity, and \(\zeta\) is the relativistic anomaly. This is commonly decomposed into \(\zeta =\chi - w\), where \(\chi\) is a radial phase parameter and \(w\) is the argument of periapsis. The semilatus rectum and eccentricity can be related to the radial roots, \(r_1\) and \(r_2\), via: \[p = \frac{2 r_1 r_2}{M(r_1 + r_2)}, \quad e = \frac{r_1 - r_2}{r_1 + r_2}.\] These can be found by taking the roots of the radial potential, i.e.: \[\frac{d r}{dt} = \xi(r) \frac{\partial H_\text{EOB}}{\partial p_{r_*}} \Bigg \vert _{r_1,r_2} = 0.\]

In this parametrization, one can find expressions for the other constants of motion, such as the effective energy and angular momentum [121]: \[\begin{equation} H_{\text{eff}} =-p_t = \frac{2 \sqrt{e} \sqrt{A(\nu;r_1) A(\nu;r_2)}}{\sqrt{(1+e)^2 A(\nu;r_2) - (1 - e)^2 A(\nu;r_1)}}, \end{equation} \begin{equation} \label{eq:pphi95of95p95e} p_\phi = \sqrt{\frac{p^2 M^2 (A(\nu;r_2) - A(\nu;r_1))}{(1-e)^2 A(\nu;r_1) - (1+e)^2 A(\nu;r_2)}}. \end{equation}\tag{5}\] Inverting the original expression for \(H_{\text{eff}}(\nu;r,p_{r_*},p_\phi)\) yields an expression for \(p_{r_*}(\nu;p,e,\zeta)\). This can be done exactly in the test mass limit or when the Q potential consists only of its leading order term [121]. However, including more terms in the Q potential generally precludes an exact analytical inverse, so one must numerically root find for \(p_{r_*}(\nu;p,e,\zeta)\). We avoid this by exploiting that \(Q\) is a power series expansion in \(p_{r_*}\) and perturbatively invert the series to obtain an approximate expression for \(p_{r_*}(\nu;p,e,\zeta)\). Since Q is known to \(\mathcal{O}(p_{r_*}^8)\), this approximation remains very accurate even at large eccentricities. We give an explicit expression for this expansion in Appendix 9.

3.2 The method of osculating orbital elements (OOE)↩︎

To uniquely specify an eccentric, planar orbit, one need specify \(p\), \(e\), and \(w\). We denote these as the orbital elements \(I^A = (p,e,w)\). In the absence of radiation reaction, they are constant.

We now recast the EOB equations of motion into evolution equations for the orbital elements via the method of osculating orbital elements (OOE) [117], [118]. At each instant in time (\(t\)) one maps the worldline of the inspiraling test mass \(x^\alpha\) onto a tangent conservative orbit sharing the same orbital elements \(I^A\) via the osculating conditions: \[\begin{equation} x^\alpha(t) = x_\text{cons}^\alpha(I^A(t), t) \end{equation} \begin{equation} \frac{d x^\alpha}{d t} (t) = \frac{\partial x_\text{cons}^\alpha}{\partial t}(I^A(t), t). \end{equation}\] From these, one can derive the osculating orbital equations: \[\begin{equation} \frac{\partial x_{\text{cons}}^\alpha}{\partial I^A} \frac{d I^A}{d t} = 0, \end{equation} \begin{equation} \frac{\partial p^{\text{cons}}_\alpha}{\partial I^A} \frac{d I^A}{d t} = F_{\alpha}, \end{equation}\] For simplicity, we drop the “cons” subscript henceforth, but it is important to remember that all quantities in the above equations are evaluated on the tangent conservative orbit. We also make use of the relationship: \[F_t = - \frac{1}{\mu}\frac{d H_{\text{EOB}}}{d t} = \left( \frac{d r}{dt}F_r + \frac{d \phi}{d t} F_\phi \right).\] With this in hand, using the second osculating condition gives us the following simultaneous equations: \[\begin{equation} \frac{\partial p_t}{\partial p} \frac{d p}{d t} + \frac{\partial p_t}{\partial e} \frac{d e}{d t} = - \frac{d H_{\text{eff}}}{d H_{\text{EOB}}} \frac{d H_{\text{EOB}}}{d t} = \frac{H_{\text{EOB}}}{M} F_t \end{equation} \begin{equation} \frac{\partial p_\phi}{\partial p} \frac{d p}{d t} + \frac{\partial p_\phi}{\partial e} \frac{d e}{d t} = F_\phi, \end{equation}\] which has the solution: \[\begin{align} \begin{aligned} \frac{d p}{d t} & = \frac{1}{D} \left( \frac{\partial p_t}{\partial e} F_\phi - \frac{\partial p_\phi}{\partial e} \frac{H_{\text{EOB}}}{M} F_t \right), \end{aligned} \\ \begin{align} \frac{d e}{d t} & = \frac{1}{D} \left( \frac{\partial p_t}{\partial p} F_\phi - \frac{\partial p_\phi}{\partial p} \frac{H_{\text{EOB}}}{M} F_t \right), \end{align} \\ \begin{align} D & = \frac{\partial p_t}{\partial p} \frac{\partial p_\phi}{\partial e} - \frac{\partial p_t}{\partial e} \frac{\partial p_\phi}{\partial p}. \end{align} \end{align}\] Here \(D\) corresponds to the Jacobian determinant of the transformation from \((H_\text{eff},p_\phi)\) to \((p,e)\) which is singular as \(e \rightarrow 0\).

Using the radial component of the first osculating condition gives us: \[\frac{\partial r}{\partial p} \frac{d p}{d t} + \frac{\partial r}{\partial e} \frac{d e}{d t} + \frac{\partial r}{\partial w} \frac{d w}{d t} = 0\] Rearranging the above gives us: \[\label{eq:dwdt} \frac{d w}{d t} = - \frac{1}{ \partial r / \partial w} \left( \frac{\partial r}{\partial p} \frac{d p}{d t} + \frac{\partial r}{\partial e} \frac{d e}{d t} \right).\tag{6}\] Note that \(\partial r / \partial w \rightarrow 0\) in the circular-orbit limit \(e \rightarrow 0\), making the above expression singular. ¹

For planar orbits that we wish to orbit average, it is far more convenient to parametrize the system with the radial phase parameter \(\chi\) rather than with the time coordinate \(t\) [123]. Using the chain rule, we find that \[\frac{d t}{d \chi} = \frac{d r / d \chi}{d r / dt } = \frac{d r / d \chi}{ \xi(r) \partial H_\text{EOB} / \partial p_{r_*} }.\] Note that this expression yields \(0/0\) at turning points where \(d r / d \chi = d r / dt =0\). In the Schwarzschild limit, one can simplify this expression to remove these unphysical singularities, but this is not always possible in the EOB case.

Moreover, rather than evolving \(w\) on its own, we use \(\zeta = \chi - w\) to absorb its evolution into that of the phase \(\zeta\), thus reducing the orbital element parameter space by one. Overall, we end up with equations of motion that take the following form: \[\begin{align}\label{eq:OOE95EoM} \begin{aligned} \frac{d p}{d \chi} &= \frac{d p }{d t} \frac{d t}{d \chi} = \epsilon F_p (\nu; p,e,\zeta), \end{aligned} \\ \begin{align} \frac{d e}{d \chi} &= \frac{d e }{d t} \frac{d t}{d \chi} = \epsilon F_e (\nu; p,e,\zeta) , \end{align}\\ \begin{align} \frac{d \zeta}{d \chi} &= \frac{d \chi}{d \chi} - \frac{d w }{d t} \frac{d t}{d \chi} = 1 + \epsilon f_{w} (\nu; p,e,\zeta) , \end{align} \\ \begin{aligned} \frac{d t}{d \chi} & = f_t (\nu; p,e,\zeta) \end{aligned}\\ \begin{align} \frac{d \phi}{d \chi} &= \frac{d \phi }{d t} \frac{d t}{d \chi} = f_\phi (\nu; p,e,\zeta) \end{align}, \end{align}\tag{7}\] where the bookkeeping parameter \(\epsilon\) marks which terms depend on radiation reaction. We summarize osculating-orbital-element equations schematically as: \[\label{eq:OOE95EoM95Eccentric95Param} \begin{align} \begin{aligned} & \dot{P_j} = \epsilon F_j (\nu,\mathbf{P}, \zeta) , \end{aligned}\\ \begin{align} & \dot{\zeta} = 1 + \epsilon f_{w}(\nu; \mathbf{P}, \zeta), \end{align}\\ \begin{align} & \dot{S}_k = s_k(\nu;\mathbf{P}, \zeta), \end{align} \end{align}\tag{8}\] where we denote the orbital elements \(P_j = (p,e)\), we treat \(\zeta\) as the radial phase, \(S_k = (t,\phi)\) are the “extrinsic" quantities, and an overdot denotes the rate of change with respect to \(\chi\). We stress that at this stage, the only approximation that has been made is for \(p_{r*}(\nu;p,e,\zeta)\); otherwise these equations are exactly equivalent to the original EOB equations of motion.

3.3 Regularization of forcing terms via Fourier decomposition↩︎

The forcing terms \(F_j\), \(f_{w}\), and \(s_k\) are all periodic functions of \(\zeta\) with period \(2 \pi\). However, they also have unphysical numerical singularities at the radial turning points where \(\zeta = n \pi\) and \(n \in \mathbb{Z}\). These singularities must be regularized before the equations of motion can be integrated numerically. We decompose functions into Fourier series using the convention: \[\label{eq:Fourier} A(\mathbf{P},\zeta) = \sum_{\kappa} A_{\kappa}(\mathbf{P}) e^{i \kappa \zeta},\tag{9}\] where \(A\) represents any of the functions \(F_j\), \(f_{w}\), or \(s_k\). Based on this, we can split the function into an averaged piece \[\label{eq:average} \left\langle A \right\rangle (\mathbf{P}) = A_{0}(\mathbf{P}) = \frac{1}{(2\pi)} \int_0^{2\pi} A(\mathbf{P},\zeta) d \zeta,\tag{10}\] and an oscillating piece \[\label{eq:oscillating} \Breve{A}(\mathbf{P},\zeta) = A(\mathbf{P},\zeta) - \left\langle A \right\rangle(\mathbf{P}) = \sum_{\kappa \neq 0 } A_{\kappa}(\mathbf{P}) e^{i \kappa \zeta}.\tag{11}\]

One way to compute the Fourier coefficients \(A_\kappa(\mathbf{P})\) is to sample \(A(\mathbf{P},\zeta)\) at equal steps in \(\zeta\) over \([0,2 \pi ]\) and apply a fast Fourier transform (FFT). However, this requires evaluating the function at the turning points where it is singular. We therefore sample at the same steps but with a constant offset \(\delta = 0.1\) to avoid the singular points. Using an FFT on this list gives us the coefficients \(\hat{A}_\kappa(\mathbf{P})\), related to the original coefficients via \[A_\kappa(\mathbf{P}) = \exp(i \kappa \delta) \hat{A}_\kappa(\mathbf{P}).\] We can reconstruct the original function \(A(\mathbf{P},\zeta)\) via Eq. 9 without ever evaluating it at the singular points.

Moreover, the Fourier coefficients decay exponentially with \(|\kappa|\) so only a finite number of modes are needed to accurately reconstruct the function. For our work, we keep only \(\kappa_\text{max} = 30\). Functions such as \(f_t\) have Fourier modes that decay slowly with \(\kappa\) but their reciprocal \(f_t^{-1}\) decays much faster. We exploit this fact to accurately reconstruct \(f_t\) using the method outlined in Appendix 10.

To obtain smooth functions over the parameter space, we numerically compute the Fourier coefficients for each of the forcing functions \(F_j\), \(f_{w}\), and \(s_k\) on a grid in \(\nu\), \(p\) and \(e\). This yields regularized forcing functions evaluable anywhere in the parameter space, allowing us to integrate the osculating orbital element equations numerically.

The resulting trajectories are much slower to compute than the EOB trajectories and do not match them perfectly owing to numerical errors from interpolation, Fourier-mode truncation, and the numerical ODE solver. Since we use the same interpolation grid, number of Fourier modes, and ODE solver throughout, comparing the OOE trajectories with the full EOB trajectories acts as a control for these errors and isolates the errors introduced by the near-identity averaging transformations of the next section.

4 Near-Identity Averaging Transformations↩︎

Near-identity (averaging) transformations (NITs) are a well-known technique in applied mathematics and celestial mechanics [122]. This technique involves making small transformations to the equations of motion, so that the short timescale physics is averaged out, while the long term evolution is retained. This is well suited to EMRIs, where we do not need the trajectory on the orbital timescale, so long as we can accurately track its evolution on the radiation-reaction timescale [123]–[131]. We will show that NITs are also well suited to accelerating eccentric EOB inspiral dynamics.

We first define NIT variables \(\tilde{P}_j\), \(\tilde{\zeta}\) and \(\tilde{S}_k\), which are related to the OOE variables \(P_j\), \(\zeta\) and \(S_k\) via the following near-identity transformations: \[\tag{12} \begin{align} \begin{aligned}\tag{13} \tilde{P}_j &= P_j + \epsilon Y_j^{(1)}(\mathbf{P}, \zeta) + \epsilon^2 Y_j^{(2)}(\mathbf{P}, \zeta) + \mathcal{O}( \epsilon^{ 3 } ), \end{aligned}\\ \begin{align} \tilde{\zeta} &= \zeta + \epsilon X^{(1)}(\mathbf{P}, \zeta) + \mathcal{O}( \epsilon^{ 2 } ), \end{align}\\ \begin{align} \tilde{S}_k &= S_k + Z_k^{(0)}(\mathbf{P}, \zeta) +\epsilon Z_k^{(1)}(\mathbf{P}, \zeta) + \mathcal{O}( \epsilon^{ 2 } ). \end{align} \end{align}\] Note that the above functions also depend on the symmetric mass ratio \(\nu\), which we suppress for notational convenience. Here, the transformation functions \(Y_j^{(n)}\), \(X_i^{(n)}\), and \(Z_k^{(n)}\) are required to be smooth, periodic functions of the radial phase \(\zeta\). At leading order, Eqs. 12 are identity transformations for \(P_k\) and \(\zeta\) but not for \(S_k\) owing to the zeroth order transformation term \(Z_k^{(0)}\). The inverse transformations can be found for \(P_j\) and \(\zeta\) by requiring that their composition with the transformations in Eqs. 12 must recover the identity transformation. Expanding order by order in \(\epsilon\), this gives us \[\tag{14} \begin{align} \begin{aligned} P_j & = \tilde{P_j} - \epsilon Y_j^{(1)}(\mathbf{\tilde{P}},\tilde{ \zeta}) - \epsilon^2 \tilde{Y}_j^{(2)}(\mathbf{\tilde{P}},\tilde{ \zeta}) + \mathcal{O}(\epsilon^3), \end{aligned}\\ \begin{align} \zeta & = \tilde{ \zeta} - \epsilon X^{(1)}(\mathbf{\tilde{P}},\tilde{ \zeta}) + \mathcal{O}(\epsilon^2), \end{align}\\ \begin{align} \tag{15} S_k & = \tilde{S}_k -Z_k^{(0)}(\mathbf{\tilde{P}},\tilde{ \zeta}) - \epsilon \tilde{Z}_k^{(1)}(\mathbf{\tilde{P}},\tilde{ \zeta}) + \mathcal{O}(\epsilon^2), \end{align} \end{align}\] where \[\label{eq:inverse95tramsformation95definitions} \begin{equation} \begin{align} \tilde{Y}_j^{(2)} &= \left( Y_j^{(2)} - \frac{\partial Y_j^{(1)}}{\partial \tilde{P_k}} Y_k^{(1)} - \frac{\partial Y_j^{(1)}}{\partial \tilde{ \zeta}} X^{(1)} \right), \end{align} \end{equation} \begin{equation} \begin{align} \tilde{Z}_k^{(1)} &= \left( Z_k^{(1)} - \frac{\partial Z_k^{(0)}}{\partial \tilde{P_k}} Y_k^{(1)} - \frac{\partial Z_k^{(0)}}{\partial \tilde{ \zeta}} X^{(1)} \right). \end{align} \end{equation}\tag{16}\]

To derive the averaged equations of motion, one takes the \(\chi\) derivative of Eqs. 12 , substitutes the OOE equations of motion 8 , converts from OOE variables to NIT variables using Eqs. 14 , and expands order by order in \(\epsilon\). At each order, one then fixes the form of the transformation functions \(Y_j^{(n)}, X^{(n)}\) and \(Z_k^{(n)}\) such that the resulting equations of motion for \(\tilde{P}_j, \tilde{\zeta}\) and \(\tilde{S}_k\) take the following form \[\label{eq:transformed95EoM} \begin{align} \begin{aligned} \dot{\tilde{P}}_j & = \epsilon \tilde{F}_j^{(0)}(\mathbf{\tilde{P}}) + \epsilon^2 \tilde{F}_j^{(1)}(\mathbf{\tilde{P}}) + \epsilon^3 \tilde{F}_j^{(2)}(\mathbf{\tilde{P}}) + \mathcal{O}( \epsilon^{ 4 } ), \end{aligned} \\ \begin{align} \dot{\tilde{\zeta}} & = 1 +\epsilon \tilde{f}_{w}^{(1)}(\mathbf{\tilde{P}}) +\epsilon^2 \tilde{f}_{w}^{(2)}(\mathbf{\tilde{P}})+ \mathcal{O}( \epsilon^{ 3 } ), \end{align}\\ \begin{align} \dot{\tilde{S}}_k & = \tilde{s}_k^{(0)}(\mathbf{\tilde{P}}) + \epsilon\tilde{s}_k^{(1)}(\mathbf{\tilde{P}}) + \epsilon^2 \tilde{s}_k^{(2)}(\mathbf{\tilde{P}}) + \mathcal{O}( \epsilon^{ 3 } ). \end{align} \end{align}\tag{17}\] where the superscripts in parentheses indicate the post-adiabatic (PA) order of the terms.

Here we use "post-adiabatic order" to refer to corrections to the adiabatic evolution of the orbital elements of an eccentric inspiral that arise from the NIT expansion. Leading order in the expansion corresponds to the adiabatic order (0PA), next-to-leading order corresponds to the first post-adiabatic order (1PA), and next-to-next-to-leading order corresponds to the second post-adiabatic order (2PA).

Note that we do not perform an expansion in the mass-ratio. As such, our usage differs slightly from that of gravitational self-force, where the post-adiabatic order is also tied to the corresponding order of the self-force in the mass-ratio expansion. It also differs from the quasi-circular post-adiabatic expansion of Refs. [109]–[112], which expands about the adiabatic evolution of a quasicircular inspiral.

Crucially, the above equations of motion are now independent of the orbital phase \(\zeta\). In this procedure, one fixes only the oscillating parts of the transformation functions and is free to choose their orbit averaged parts. Many such choices are explored in Ref. [123], but for simplicity we set all the orbit averaged parts of the transformation functions to be zero. For brevity, we provide an ancillary Mathematica file that derives the relationships between (\(\tilde{F}_j^{(n)},\tilde{f}_w^{(n)}\) and \(\tilde{s}_k^{(n)}\)) to the terms in the OOE equations of motion and summarize the results below.

4.1 Adiabatic terms↩︎

The adiabatic (0PA) order averaged terms are simply given by: \[\begin{equation} \tilde{F}_j^{(0)} = \left<F_{j}\right>, \end{equation} \begin{equation} \label{eq:0PA95f95t95f95phi} \tilde{s}_k^{(0)} = \left<s_{k}\right>. \end{equation}\tag{18}\] At this order, the model reduces to the standard adiabatic model where one simply orbit averages the forcing terms in the OOE equations of motion. In deriving these relations we have constrained the oscillating pieces of the NIT transformation functions to be \[\begin{equation}\tag{19} \Breve{Y}_j^{(1)} = \sum_{\kappa \neq 0} \frac{i}{\kappa} F_{j,\kappa} e^{i \kappa \zeta}, \end{equation} \begin{equation}\tag{20} \Breve{Z}_k^{(0)} = \sum_{\kappa \neq 0} \frac{i}{\kappa} s_{k,\kappa}e^{i \kappa \zeta}. \end{equation}\] This leading-order transformation between OOE and NIT variables is crucial for setting consistent initial conditions, ensuring we are comparing the same physical systems. Thus, we refer to the model that uses only the above 0PA equations of motion and leading-order transformations for the initial conditions as the “0PA model."

4.2 First post-adiabatic terms↩︎

At first postadiabatic order (1PA), one obtains: \[\label{eq:1PA95Averaged95Terms95Ecc} \begin{equation} \tilde{F}_j^{(1)} = \left<\frac{\partial \Breve{Y}_j^{(1)}}{\partial \tilde{P}_k} \Breve{F}_k \right> + \left<\frac{\partial \Breve{Y}_j^{(1)}}{\partial \tilde{\zeta}} \Breve{f}_{w} \right>, \end{equation} \begin{equation} \tilde{f}_{w}^{(1)}= \left<f_{w}\right>, \end{equation} \begin{equation} \tilde{s}_k^{(1)} = \left<\frac{\partial \Breve{Z}_k^{(0)}}{\partial \tilde{P}_j} \Breve{F}_j \right> + \left<\frac{\partial \Breve{Z}_k^{(0)}}{\partial \tilde{\zeta}} \Breve{f}_{w} \right> , \end{equation}\tag{21}\] where we have constrained \[\begin{align} \begin{aligned}\tag{22} \Breve{Y}_j^{(2)} = & \sum_{\kappa \neq 0} \frac{ i e^{ i \kappa \zeta}}{\kappa} \Biggl[ - \frac{\partial \tilde{F}_j^{(0)}}{\partial \tilde{P}_k} \Breve{Y}_{k, \kappa}^{(1)} + \\ & \sum_{\kappa' \neq 0} \left( \frac{\partial \Breve{Y}_{j, \kappa}^{(1)}}{\partial \tilde{P}_k} \Breve{F}_{k,\kappa - \kappa'} + \frac{\partial \Breve{Y}_{j,\kappa'}^{(1)}}{\partial \tilde{\zeta}} \Breve{f}_{w,\kappa - \kappa'} \right) \Biggr] \end{aligned} \\ \begin{align}\tag{23} \Breve{X}^{(1)} = & \sum_{\kappa \neq 0} \frac{i}{\kappa} f_{w,\kappa} e^{i \kappa \zeta}, \end{align}\\ \begin{align}\tag{24} \Breve{Z}_k^{(1)} = & \sum_{\kappa \neq 0} \frac{ i e^{ i \kappa \zeta}}{\kappa} \Biggl[ - \frac{\partial \tilde{s}_k^{(0)}}{\partial \tilde{P}_j} \Breve{Y}_{j, \kappa}^{(1)} + \\ & \sum_{\kappa' \neq 0} \left( \frac{\partial \Breve{Z}_{k, \kappa}^{(0)}}{\partial \tilde{P}_j} \Breve{F}_{j,\kappa - \kappa'} + \frac{\partial \Breve{Z}_{k,\kappa'}^{(0)}}{\partial \tilde{\zeta}} \Breve{f}_{w,\kappa - \kappa'} \right) \Biggr] \end{align} \end{align}\] We refer to the model using these 1PA equations of motion and together with the sub-leading transformations as the “1PA model".

There is some ambiguity whether these sub-leading transformations \(\Breve{Y}_j^{(2)}\), \(\Breve{X}^{(1)}\), and \(\Breve{Z}_k^{(1)}\) belong at 1PA order, or whether only the leading order transformations should be included. We find that without them, the 1PA model is no more accurate than the 0PA model. We therefore choose to include them in our 1PA model.

4.3 Second post-adiabatic terms↩︎

To reach acceptable accuracy for comparable-mass systems, we extend the calculation of the NIT equations of motion to second postadiabatic order (2PA). The resulting averaged terms are: \[\begin{align} \label{eq:2PA95Averaged95Terms95Ecc} \begin{aligned} \tilde{F}_j^{(2)} = & \left<\frac{\partial \Breve{Y}_j^{(2)}}{\partial \tilde{P}_k} \Breve{F}_{k} \right> + \left<\frac{\partial \Breve{Y}_j^{(2)}}{\partial \tilde{\zeta}} \Breve{f}_{w} \right> \\ & - \frac{1}{2} \frac{\partial^2 \tilde{F}_j^{(0)}}{\partial \tilde{P}_k \partial \tilde{P}_i} \left< \Breve{Y}_k^{(1)} \Breve{ Y}_i^{(1)} \right> , \end{aligned}\\ \begin{align} \tilde{f}_{w}^{(2)} = & \left<\frac{\partial \Breve{X}^{(1)} }{\partial \tilde{P}_j} F_j \right> + \left<\frac{\partial \Breve{X}^{(1)} }{\partial \tilde{\zeta}} f_{w} \right>, \end{align}\\ \begin{align} \tilde{s}_k^{(2)} = & \left<\frac{\partial \Breve{Z}_k^{(1)}}{\partial \tilde{P}_j} \Breve{F}_{j} \right> + \left<\frac{\partial \Breve{Z}_k^{(1)}}{\partial \tilde{\zeta}} \Breve{f}_{w} \right> \\ & - \frac{1}{2} \frac{\partial^2 \tilde{s}_k^{(0)}}{\partial \tilde{P}_i \partial \tilde{P}_j} \left< \Breve{Y}_i^{(1)} \Breve{ Y}_j^{(1)} \right> . \end{align} \end{align}\tag{25}\] We refer to the model that adds these 2PA contributions to the 1PA as the “2PA model". Any accuracy gain from 1PA to 2PA can be attributed to these terms.

At this order, the equations of motion become ill-behaved in the circular-orbit limit \(e \rightarrow 0\) where \(f_{w}\) diverges as \(1/e\). Moreover, the osculating-orbit eccentricity ceases to behave as a smooth orbit averaged evolution plus a rapidly oscillating contribution. This phenomenon will be explored in a follow-up paper [134].

4.4 Summary of results↩︎

Combining the above results, we obtain the averaged equations of motion for the NIT variables through second post-adiabatic order: \[\begin{align} \begin{aligned} \frac{d \tilde{p}}{d \chi} &= \epsilon \tilde{F}_p^{(0)}(\tilde{p},\tilde{e}) + \epsilon^2 \tilde{F}_p^{(1)}(\tilde{p},\tilde{e}) + \epsilon^3 \tilde{F}_p^{(2)}(\tilde{p},\tilde{e}) + \mathcal{O}( \epsilon^{ 4 } ) , \end{aligned} \\ \begin{align} \frac{d \tilde{e}}{d \chi} &= \epsilon \tilde{F}_e^{(0)}(\tilde{p},\tilde{e}) + \epsilon^2 \tilde{F}_e^{(1)}(\tilde{p},\tilde{e}) + \epsilon^3 \tilde{F}_e^{(2)}(\tilde{p},\tilde{e}) + \mathcal{O}( \epsilon^{ 4 } ), \end{align}\\ \begin{align} \frac{d \tilde{\zeta}}{d \chi} &=1 + \epsilon \tilde{f}_{w}^{(1)}(\tilde{p},\tilde{e}) + \epsilon^2 \tilde{f}_{w}^{(2)} (\tilde{p},\tilde{e}) + \mathcal{O}( \epsilon^{ 3 } ), \end{align} \\ \begin{aligned} \frac{d \tilde{t}}{d \chi} & = f_t^{(0)} (\tilde{p},\tilde{e}) + \epsilon \tilde{f}_t^{(1)} (\tilde{p},\tilde{e}) + \epsilon^2 \tilde{f}_t^{(2)} (\tilde{p},\tilde{e}) + \mathcal{O}( \epsilon^{ 3 } ) \end{aligned}\\ \begin{align} \frac{d \tilde{\phi}}{d \chi} & = f_\phi^{(0)} (\tilde{p},\tilde{e}) + \epsilon \tilde{f}_\phi^{(1)} (\tilde{p},\tilde{e}) + \epsilon^2 \tilde{f}_\phi^{(2)} (\tilde{p},\tilde{e}) + \mathcal{O}( \epsilon^{ 3 } ) \end{align} \end{align}\] One can then relate the NIT variables to the OOE variables to ensure consistent initial conditions via \[\begin{align} \label{eq:OOE95NIT95Initial95Conditions} \begin{aligned} \tilde{p} &= p + \epsilon Y_p^{(1)}(p,e,\zeta ) + \epsilon^2 Y_p^{(2)}(p,e,\zeta ) + \mathcal{O}( \epsilon^{ 3 } ) , \end{aligned} \\ \begin{align} \tilde{e} &= e + \epsilon Y_e^{(1)}(p,e,\zeta ) + \epsilon^2 Y_e^{(2)}(p,e,\zeta ) + \mathcal{O}( \epsilon^{ 3 } ) \end{align}\\ \begin{align} \tilde{\zeta} &= \zeta + \epsilon X^{(1)}(p,e,\zeta ) + \mathcal{O}( \epsilon^{ 2 } ), \end{align} \\ \begin{aligned} \tilde{t} & = t + Z_t^{(0)} (p,e,\zeta ) + \epsilon Z_t^{(1)} (p,e,\zeta ) + \mathcal{O}( \epsilon^{ 2 } ) \end{aligned}\\ \begin{align} \tilde{\phi} & = \phi + Z_\phi^{(0)} (p,e,\zeta ) + \epsilon Z_\phi^{(1)} (p,e,\zeta ) + \mathcal{O}( \epsilon^{ 2 } ) \end{align} \end{align}\tag{26}\] The OOE variables can then be recovered from the solutions to the NIT equations of motion using the inverse transformations of Eqs. 14 .

5 Waveform Generation↩︎

We now discuss how to generate waveforms from both the EOB and averaged dynamics. The waveform strain \(h\) can be split into polarizations \(h_+\) and \(h_\times\) which can be further decomposed into \(-2\) spin-weighted spherical-harmonic modes as follows: \[\label{eq:StrainDecomposition} h = h_+ - i h_\times = \sum_{\ell \geq 2} \sum_{m=-\ell}^{\ell} h_{\ell m} \,_{-2}Y_{\ell m}(\theta, \phi).\tag{27}\] Within the PN formalism, the waveform modes \(h_{\ell m}\) have been calculated up to 3PN order [135]–[138], and transformed into the EOB formalism in Ref. [62]. As a proof of principle, here we employ the 3PN instantaneous contributions to the modes from Ref. [62], which take the general form: \[h_{\ell m} = \mathcal{A}_{\ell m} \left(r,p_r,p_\phi \right) \exp \left( - i m \phi \right).\] We model only the dominant \((\ell,m) = (2,2)\) mode, so we compute only \(h_{22}\). We provide explicit expressions for the \(\mathcal{A}_{22}\) mode amplitudes up to 3PN order for nonspinning systems in Appendix 12.

5.1 Waveforms from EOB dynamics↩︎

To generate waveforms from the EOB dynamics, one first solves the EOB equations of motion for the variables \((r(t),p_{r_*}(t), \phi(t), p_\phi(t))\) using standard numerical ODE solvers such as adaptive step-size Runge-Kutta methods.

Once solved, the EOB dynamics yield the solutions \((r(t),p_{r_*}(t), \phi(t), p_\phi(t))\) on non-uniform time steps. At each time step, one substitutes these solutions into the waveform-mode expressions, here \(h_{2 2}\). Finally, the modes are interpolated onto a dense uniformly sampled time grid that can then be Fourier transformed for mismatch and likelihood calculations.

5.2 Waveforms from NIT dynamics↩︎

When using the NIT equations of motion, one obtains solutions for the NIT variables \((\tilde{p}(\chi), \tilde{e}(\chi), \tilde{\zeta}(\chi), \tilde{t}(\chi), \tilde{\phi}(\chi))\) on sparse, non-uniform steps in \(\chi\). Substituted directly into the waveform mode expressions, these are too sparse to resolve the orbital-timescale oscillations of the waveform. We therefore first sample the NIT solutions onto a uniform grid in \(\chi\), with spacing \(\Delta \chi\) small enough to resolve the orbital-timescale oscillations as set by the Nyquist theorem: \[\label{eq:Nyquest} \Delta \chi \leq \frac{2 \pi}{2 \kappa_\text{step} + 1}.\tag{28}\] Here \(\kappa_\text{step}\) is the highest Fourier mode present in the inverse transformation for the extrinsic NIT variable \(t\) given by Eq. 15 . We determine this by examining the ratio of the Fourier coefficients of the leading transformation until we find \[\left| \frac{Z_{t,\kappa_\text{step}}^{(0)}(p_0,e_0) }{ Z_{t,1}^{(0)} (p_0,e_0)} \right| \leq 10^{-16},\] as any smaller mode is of the same order of magnitude as machine precision. To ensure we never undersample our waveforms at lower eccentricities, we take the maximum of the above criterion and \(\kappa_\text{step} = 20\). Thus for lower eccentricities, we can sample at lower rates than for higher eccentricity systems. This condition can be relaxed to larger values of \(\Delta \chi\), further reducing waveform-generation cost at some loss of accuracy. However, since our priority is to test the accuracy of the averaged dynamics, we choose a conservative sampling criterion.

In principle, one should also transform from the intrinsic NIT variables \((\tilde{p}, \tilde{e}, \tilde{\zeta})\) to the osculating orbital element variables \((p, e, \zeta)\) on this uniform \(\chi\) grid using Eqs. 14 . We find that as in the extreme-mass-ratio case, this has negligible impact on waveform accuracy and can be safely neglected [123]–[125], [127], [128], [130]. Thus one can reconstruct the EOB variables \((r(\chi), p_{r_*}(\chi), \phi(\chi), p_\phi(\chi))\) using the relations of Sec. 3.1 with the approximation \((p, e, \zeta) \sim (\tilde{p}, \tilde{e}, \tilde{\zeta})\). Using the mode expressions in Appendix 12 one obtains the waveform on a non-uniform time grid. Finally, this is interpolated and upsampled onto a dense uniform time grid that can then be Fourier transformed for mismatch and likelihood calculations.

6 Implementation↩︎

We now outline the implementation of our Mathematica code for solving the dynamics and generating waveforms. This covers the interpolation grid, the initial conditions, the numerical integration, and the handover conditions for switching from NIT to EOB dynamics. We also describe how we compute mismatches between NIT/OOE and EOB waveforms.

6.1 Interpolation grid↩︎

To obtain smooth functions over the parameter space, we numerically calculate the Fourier coefficients via a fast Fourier transform (FFT) for each of the forcing functions \(F_j\), \(f_{w}\), and \(s_k\) on a grid in \(\nu\), \(p\) and \(e\). Motivated by the fact that the rate of change of \(p\) with respect to time scales as \(p^{-3}\) in the weak field limit [28], we define the quantity \(u = \log \left( p^{-3} \right)\). An equally spaced grid in \(u\) covers both the strong- and weak-field regimes with sufficient resolution.

We then construct a grid that has uniform spacing in \(\nu\), \(e\), and \(u\) given by

• \(\nu \in [0.01,0.25 ]\) with \(\Delta \nu = 0.01\),

• \(e \in [ 0.01, 0.9]\) with \(\Delta e = 0.01\),

• \(u \in [ -17, \log \left((6 + 2 e_\text{max})^{-3}\right)]\) (which corresponds to \(p \in [6 + 2 e_\text{max}, 289.069 ]\)) with \(\Delta u = 0.05\).

This results in a grid of \(25 \times 90 \times 217 = 488,250\) points. We then use cubic-spline interpolation to evaluate the Fourier coefficients at any point in the parameter space.

For the postadiabatic transformation functions and averaged terms, we combine the Fourier coefficients of the OOE forcing functions with the relationships derived in Sec. 4 on the same grid, then use cubic-spline interpolation to evaluate them rapidly at any point in the \((\nu,p,e)\) parameter space.

6.2 Initial conditions↩︎

To facilitate comparisons between the EOB, OOE, and NIT dynamics, we must specify consistent initial conditions for each system. The initial conditions for the OOE system are \((p_0, e_0, \zeta_0, \phi_0)\) at an initial time \(t_0\). This corresponds to the purely conservative orbit to which the inspiral is instantaneously tangential. Often rather than \(p_0\), we specify the initial orbit averaged waveform frequency \(f_0\), and root-find numerically for the value of \(p_0\) that satisfies \[f_0 = \frac{2}{2 \pi} \left\langle \Omega_\phi \right\rangle = \frac{1}{\pi} \left( \frac{\tilde{f}_\phi^{(0)}(\nu,p_0,e_0)}{\tilde{f}_t^{(0)}(\nu,p_0,e_0)} \right).\] Here \(\tilde{f}_\phi^{(0)}\) and \(\tilde{f}_t^{(0)}\) are the leading order terms in the NIT equations of motion for the extrinsic variables \(\phi\) and \(t\) respectively as given by Eq. 18 .

To find the corresponding initial conditions for the EOB system \((r_0, p_{r_* 0}, \phi_0, p_{\phi0})\), one can use the relations Eq. 4 , Eq. 5 , and the relation for the initial radial momentum, given in Appendix 9².

To set the initial conditions for the NIT system, we use Eqs. 26 to transform the OOE initial conditions into the corresponding NIT ones. We use only the leading-order terms for the 0PA model and the sub-leading terms for higher PA models. This ensures that the initial conditions of all systems are consistent and correspond to the same physical inspiral, enabling meaningful comparisons between the resulting waveforms.

6.3 Numerical integrator↩︎

To solve the equations of motion, we use the NDSolve function in Mathematica with the method set to ExplicitRungeKutta with AccuracyGoal set to \(9\) and PrecisionGoal set to \(\infty\) to ensure that numerical errors from the ODE solver do not dominate the phase error for the long inspiral trajectories computed here. For the EOB equations of motion, we set the DifferenceOrder to \(4\) (i.e., RK4) to ensure numerical stability for long, eccentric inspirals. For the OOE and NIT equations of motion, the system of equations is much less stiff allowing DifferenceOrder \(8\) (i.e., RK8) and hence larger time steps and faster runtimes. When transitioning back to the EOB dynamics for plunge, we again use RK4 for stability.

6.4 Handover conditions↩︎

The OOE and NIT equations of motion are valid only during the inspiral and break down as the system approaches plunge and merger. We therefore specify handover conditions for switching from NIT to EOB dynamics to model the late inspiral and plunge accurately.

First we use the adiabaticity parameter given by \(\mathcal{Q} = \left\langle \Omega_\phi \right\rangle^2 / \left\langle \dot{\Omega}_\phi \right\rangle \gg 1\). We find empirically that the NIT equations of motion become significantly less accurate when \(\mathcal{Q} \lesssim 100\) though the resulting waveform accuracy is not highly sensitive to this exact threshold. To evaluate this quickly at each time step, we use the following approximation \[\begin{align} \mathcal{Q} & \sim \frac{\left(\tilde{f}^{(0)}_\phi / \tilde{f}^{(0)}_t \right)^2}{ \frac{d}{dt} \left(\tilde{f}^{(0)}_\phi / \tilde{f}^{(0)}_t \right)} \\& \sim \left(\tilde{f}^{(0)}_\phi \right)^2 \tilde{f}^{(0)}_t \left[ \left( \frac{\partial \tilde{f}^{(0)}_\phi}{\partial \tilde{P}_j} \tilde{f}^{(0)}_t - \frac{\partial \tilde{f}^{(0)}_t}{\partial \tilde{P}_j} \tilde{f}^{(0)}_\phi \right) \tilde{F}_j^{(0)} \right]^{-1}. \end{align}\]

This condition is too permissive for low-eccentricity, large-mass-ratio systems, and so we add another handover condition when \(\tilde{p} = 8.5\).

Finally, we also impose a handover when \(\tilde{e} < 0.035\) since the NIT equations become ill-behaved as \(e \rightarrow 0\) where \(f_{w}\) diverges as \(1/e\). This condition will be replaced with a transformation to NIT variables that remain well-behaved in the small-eccentricity limit, derived in a follow-up paper [134].

In summary, when evolving the NIT equations of motion, we check at each time step if any of the following conditions are violated:

• \(\tilde{p} < 8.5\)

• \(\mathcal{Q} < 100\)

• \(\tilde{e} < 0.035\)

If any of these is violated, we switch from the NIT to the EOB equations of motion. We have chosen these conditions conservatively, to prioritize the accuracy of the NIT waveforms compared to pure EOB waveforms. More optimal handover conditions undoubtedly exist that retain this accuracy while allowing a longer NIT evolution, but we leave such optimization to future work.

6.5 Mismatch calculation↩︎

To quantify the accuracy of NIT and OOE waveforms relative to EOB waveforms, we require a measure of similarity between them. The standard measure used in GW data analysis is the noise-weighted inner product between two waveforms \(h_1\) and \(h_2\), \[\left<h_1 | h_2 \right> = 4 \Re \int_{f_\text{low}}^{f_\text{high}} \frac{\tilde{h}_1(f) \tilde{h}_2^*(f)}{S_n(f)} df,\] where \(\tilde{h}(f)\) is the Fourier transform of the waveform \(h(t)\), \(\tilde{h}^*(f)\) is the complex conjugate of \(\tilde{h}(f)\), \(S_n(f)\) is the one-sided power spectral density (PSD) of the detector noise, and \(f_\text{low}\) and \(f_\text{high}\) are the lower and upper frequency cutoffs [139].

To represent current ground based detectors, we use the A+ LIGO PSD with a low-frequency cutoff of 10 Hz, corresponding to the design sensitivity for the fifth observing run (O5) [140]. We also use the EinsteinTelescope and CosmicExplorer2 PSDs from Ref. [140] with a low-frequency cutoff of 5 Hz. For LISA, we use an analytic estimate of the design sensitivity curve with a galactic confusion noise for the nominal four-year mission lifetime and a low-frequency cutoff of \(10^{-4}\) Hz as given by Eq. (1) of Ref. [141]. We use this PSD only to gain intuition into the dependence of the mismatch with expected frequency sensitivity of LISA. A rigorous assessment of the accuracy required for LISA would necessitate computing the full instrument response and projection onto TDI channels which can be explored in future work.

The inner product is then used to define the faithfulness between two waveforms as \[\mathcal{F} = \max_{t_c, \phi_c} \frac{\left<h_1 | h_2 \right>}{\sqrt{\left<h_1 | h_1 \right> \left<h_2 | h_2 \right>}},\] where the maximization over time and phase shifts \(t_c\) and \(\phi_c\) accounts for the fact that the two waveforms may not be perfectly aligned in time or phase [142]. The mismatch is then defined as \[\mathcal{M} = 1 - \mathcal{F}.\] Since we compute only the \((\ell,m) = (2,2)\) mode of the waveform, we report \(\mathcal{M}_{22}\) for the EOB, OOE, and NIT waveforms. While the subdominant modes do not become significantly more important with eccentricity, they do for large-mass-ratio systems (see Fig. 11 of Ref. [50], Fig. 14 of Ref. [104], or Fig. 16 of [143]). However, since we are primarily interested in the accuracy of the NIT dynamics rather than the physical accuracy of the waveforms, this is an acceptable measure here.

For simulations of stellar-mass BBHs as observed by ground based detectors, we use a sample rate of 4096 Hz, which gives us a high frequency resolution up to 2048 Hz via the Nyquist theorem. Since we exclude the merger and ringdown where the GW frequency sharply rises, this is more than sufficient to avoid aliasing. For simulations of massive black-hole binaries (MBHBs) as observed by LISA, we use a sample rate of 0.2 Hz (\(\Delta t = 5s\)), allowing us to resolve frequencies up to 0.1 Hz.

For the mismatch calculations, we developed a Mathematica package GWMMatch [144] to compute time- and phase- optimized mismatches using the same methods as pyCBC [145]. Since we compare inspiral only waveforms (without the merger and ringdown), we apply a Planck taper to both ends of the waveform to mitigate spectral leakage [146]. We then zero-pad the waveforms to a common length, FFT them, and compute the noise-weighted inner product with the PSDs described above. Finally, we maximize over time and phase to find the faithfulness, and thus the mismatch. We use this mismatch to quantify the accuracy of NIT and OOE waveforms relative to EOB waveforms in the following sections.

7 Results↩︎

We now present the results from our implementation of NITs for eccentric EOB inspirals. We begin with a case study comparing EOB, OOE and NIT waveforms for a specific set of parameters. We then compare the accuracy of different postadiabatic orders across parameter space for total mass \(M = 60M_\odot\) and then more closely examine the accuracy of the 2PA dynamics for several different total masses. Finally, we examine the speed-up from the NIT dynamics and how it varies with both mass ratio and total mass.

7.1 Case study↩︎

We first demonstrate the accuracy of the NIT dynamics and waveform construction for a specific set of parameters. We choose a system with total mass \(M = 20 M_\odot\), mass ratio \(q = 1\), and initial eccentricity \(e_0 = 0.75\) at an orbit-averaged reference frequency of \(5\) Hz, corresponding to an initial semilatus rectum of \(p_0 = 33.92\). This system evolves to an eccentricity \(e \sim 0.625\) at 10 Hz. We evolve to \(r = 4 M\), where we terminate the dynamics and the waveform. Note that the innermost stable circular orbit for an equal-mass binary with our choice of effective Hamiltonian is at \(r \simeq 4.45M\), so the final segment always corresponds to a plunging trajectory.

By rescaling the waveform with the total mass, we also examine the case of an MBHB system with \(M = 10^6 M_\odot\) at same mass ratio and initial eccentricity with a reference frequency of \(10^{-4}\) Hz.

In each case we generate EOB, OOE and NIT (0PA,1PA,2PA) waveforms and compute mismatches between them for O5, CE, ET and LISA using the method described in Sec. 6.5.

In Fig. 1, we compare the EOB, OOE and NIT waveforms for the case study parameters. The left panels show the start and end of the inspiral where one can see where we swap back from the NIT or OOE dynamics (dashed orange) to the EOB dynamics (dashed green). The right panels show the accumulated phase difference of the \((\ell,m) = (2,2)\) mode throughout the inspiral. The 0PA and 1PA NIT waveforms accumulate a significant phase difference of \(\sim 1\) radian by \(r = 4M\). In contrast, both the 2PA and OOE waveforms are visually indistinguishable from the EOB waveform and accumulate a phase difference of less than \(10^{-2}\) radians by the end of the inspiral.

We further test the NIT waveforms by computing their mismatches against EOB waveforms for each postadiabatic order for the case study parameters for different mass ratios \(q = (1,3,10)\). We also report the walltimes for the dynamics alone and for the full waveform generation time. The results are shown in Table 1.

First, although the mismatches vary somewhat between detector PSDs, they agree to within an order of magnitude. Conclusions drawn from the O5 PSD should therefore be broadly applicable to other detectors as well.

Second, the NIT dynamics (at all postadiabatic orders) speed up the dynamics by \(\sim 10 - 50 \times\), scaling with mass ratio. The total waveform generation speed-up is \(\sim 6 - 8 \times\). It is smaller than the speed-up of the dynamics alone because the dynamics are no longer the bottleneck - waveform generation is. This is underscored by the OOE waveforms whose total walltime is comparable to that of the EOB waveforms despite the OOE dynamics being much slower. Our new method for sampling the waveform is much more efficient than the conventional method. We will explore how interpolating the waveform modes in terms of the NIT variables can further accelerate waveform generation and reduce this bottleneck in a follow-up paper [134].

Third, the 0PA and 1PA NIT waveforms have mismatches of \(\mathcal{M}_{22} \sim 10^{-3}\) to \(10^{-4}\), which decrease with mass ratio, as expected, since this increases the adiabaticity of the inspiral. The 2PA NIT waveforms reach \(\mathcal{M}_{22} \sim 10^{-6}\) to \(10^{-7}\) making 2PA the preferred order for comparable-mass systems. As a control, we include the OOE waveforms which quantify the accuracy lost by errors arising from interpolation, the ODE solver, and our method for sampling the waveform described in Sec.5.2. Their mismatches of \(\mathcal{M}_{22} \sim 10^{-8}\) indicate that the NIT mismatches are dominated by the order at which we truncate the NIT expansion. One could push to higher PA orders to achieve further accuracy. A conservative estimate for the SNR at which the EOB and 2PA NIT waveforms would be distinguishable is given by \(\sim 1 / \sqrt{2 \mathcal{M}}\) [147]. For \(\mathcal{M}_{22} \sim 10^{-6}\), this would be an SNR of \(\sim 700\). Thus going beyond 2PA order would be unnecessary for most comparable mass binaries, bar exceptionally loud MBHB signals that one might observe with LISA [10], [11], [148].

In the following sections, we will see that these results broadly hold across the eccentric BHB parameter space for current ground-based detectors.

7.2 Comparison of different postadiabatic orders↩︎

We now examine the accuracy of different postadiabatic orders across parameter space. We fix the total mass at \(60M_\odot\). We grid \(q \in [1,96]\) in steps of 5 and \(e \in [0.1, 0.85]\) in steps of 0.05 for 320 points in total. Each system is initialized at periastron (i.e., \(\zeta_0 = 0\)) with an orbit averaged frequency of \(f_0 = 10 \text{Hz}\). We compute \(\mathcal{M}_{22}\) between EOB and the NIT/OOE waveforms using the O5 design sensitivity curve as outlined in Sec. 6.5.

Note that a small subset of these initial conditions can yield very short inspirals with low adiabaticity \(\mathcal{Q} < 100\). Such systems start after the handover from NIT to EOB dynamics (Sec. 6.4). This is especially true for high eccentricity systems with \(e_0>0.75\), suggesting that this method could be ill-suited for such high eccentricities [149]. We therefore exclude systems that immediately transition to the EOB dynamics, since their waveforms would be identical by construction.

The results of this comparison are shown in Fig. 2. Adding more postadiabatic information significantly improves waveform accuracy. The 0PA waveforms have a median mismatch of \(\sim 2 \times 10^{-3}\) with some becoming as large as \(\sim 0.5\). The 1PA waveforms improve on this, with a median \(2 \times 10^{-5}\) and a tail extending to \(6 \times 10^{-3}\). The 2PA waveforms reach a median of \(\mathcal{M}_{22} \sim 5 \times 10^{-6}\) and a maximum of \(5 \times 10^{-5}\). Thus for comparable-mass systems, at least 2PA accuracy is required.

Figure 3 shows where these mismatches occur in parameter space. The 0PA and 1PA waveforms have the highest mismatches for comparable-mass ratio and high-eccentricity systems. These systems are the least adiabatic, where the NIT expansions are expected to perform worst. 1PA may suffice even for comparable-mass systems at low eccentricities, but more information is needed for higher eccentricities.

The 2PA waveforms have low mismatches throughout, with the highest mismatches appearing sporadically across the parameter space. These errors most likely stem from the transformations back to the EOB coordinates used for the plunge, which depend on \(\zeta\) as well as \((\nu,p,e)\). This could be remedied by a more optimal choice of handover conditions than those used in Sec. 6.4. However, these mismatches remain very low, so this does not significantly affect waveform accuracy.

Finally, the OOE waveforms have low mismatches across the entire parameter space, though they are slightly larger for smaller eccentricities, with seemingly no dependence on the mass ratio. Since these are generally smaller than the 2PA mismatches, the errors arising from the low-eccentricity limit of this choice of orbital element parametrization, interpolation, the ODE solver, and our method for sampling the waveform are all subdominant to the 2PA truncation error.

7.3 Second-postadiabatic order accuracy across parameter space↩︎

We now examine the 2PA NIT accuracy across parameter space more thoroughly. We explore slices of the parameter space avoiding waveforms that are too long to generate (e.g., very low-mass, high-mass-ratio systems). For each fixed total mass, we start at periastron with orbit averaged frequency of \(f_0 = 10\)Hz and sample equally in eccentricity from \(e \in [0.1,0.85]\) in steps of \(0.05\). We sample in the mass ratio \(q\) with steps \(\Delta q\) as follows:

• \(M = 100 M_\odot\), \(q \in [1,96]\), \(\Delta q = 5\),

• \(M = 50 M_\odot\), \(q \in [1,96]\), \(\Delta q = 5\),

• \(M = 20 M_\odot\), \(q \in [1,20]\), \(\Delta q = 1\),

• \(M = 10 M_\odot\), \(q \in [1,10]\), \(\Delta q = 0.5\).

We then compute \(\mathcal{M}_{22}\) between EOB and the 2PA NIT waveforms using the O5 design sensitivity curve as in Sec. 6.5. Again a small subset of these initial conditions (especially higher-mass, high-eccentricity, and low-mass-ratio systems) can yield very short inspirals that start after the NIT handover conditions and thus are excluded from the analysis.

The results are shown in Fig. 4. The 2PA NIT waveforms perform well across the entire parameter space, with median mismatch of \(\mathcal{M}_{22} = 1.25 \times 10^{-6}\) and a maximum \(\mathcal{M}_{22} = 8.05 \times 10^{-5}\). The lowest mismatches occur for low-mass systems with long inspirals, showing that the 2PA dynamics track the phase evolution for long-lived signals accurately. The worst mismatches occur for shorter-lived, high-mass signals, indicating that most of the error is accumulated in the late inspiral and transition to plunge, rather than the early inspiral.

A further point is that, at fixed physical frequency (here 10Hz), increasing the total mass raises the geometric frequency: a system with e.g., \(e = 0.5\) is more relativistic at \(100 M_\odot\) than at \(10 M_\odot\). The \(M = 100 M_\odot\) histogram in Fig. 4 thus samples much more relativistic systems than the \(M = 10 M_\odot\) one, which may also contribute to the loss of accuracy.

Figure 5 shows contour plots of the 2PA NIT mismatches for the four total masses. For low-mass systems, the highest mismatches occur at larger eccentricities and smaller mass-ratios, where the inspirals are less adiabatic. For higher mass systems, the highest mismatches occur sporadically across parameter space, likely due to our choice of handover conditions (Sec. 6.4); a better choice could remedy this. To mitigate the errors from the \(\zeta\) dependent transformation terms in the handover to EOB variables for the plunge, we produce the same plots for the inspiral section only in Appendix 13.

Finally, the mismatch depends weakly on the mass ratio, degrading as the mass ratio decreases. Since even the \(q = 1\) mismatches remain \(< 10^{-5}\), this confirms that the 2PA NIT waveforms are accurate enough for all mass ratios.

7.4 Speed-up↩︎

Now we explore the speed-up achieved by using NIT dynamics compared to EOB dynamics across parameter space. We report the timings for the waveforms of the 2PA accuracy comparison in Sec. 7.2 and take the ratio of the EOB to NIT walltime.

In Fig. 6, we show the speed-up of producing the dynamics alone (including the changeover to the EOB transition to plunge) relative to the EOB dynamics for the parameters of Sec. 7.2. Longer-lived, low-mass inspirals show the largest speed-up, since they are dominated by the inspiral. For systems with \(M = 10M_\odot\), the speed-up for the dynamics ranges from \(6.6 - 132\times\), with a median of \(\sim 50 \times\). The effect is less dramatic for heavier, shorter-lived systems where the dynamics cost is dominated by the transition to plunge. The dynamics speed-up ranges from \(1.25 - 16.7\times\) with a median of \(\sim 6\times\) at \(100M_\odot\).

In Fig. 7 we show the total walltime speed-up (dynamics plus waveform generation) of the NIT over the EOB approach, for the same parameters. The total speed-up depends mostly on total mass and mass ratio, ranging from \(1.8 - 7.8 \times\) with a median of \(\sim 6 \times\) over the parameter space explored. It is lower than the speed-up of the dynamics alone because waveform generation, not the dynamics, is now the computational bottleneck. This can be alleviated by interpolating the waveform modes as functions of the NIT variables, which we will implement in a follow-up paper [134].

We now study how the speed-up varies with total mass by examining the dynamics-only and total walltimes for a system starting at periastron with initial orbit-averaged frequency of \(10\)Hz and eccentricity of \(e_0 = 0.4\) as total mass is varied. In Fig. 8 the dynamics cost is almost independent of total mass and most of it comes from the final transition to plunge which uses the EOB dynamics.

In Fig. 9, the total walltime does increase with total mass, as expected since the waveform-generation cost scales with the length of the signal which grows inversely with total mass. Again, most of the NIT cost comes from the final EOB transition to plunge. Despite this, using the NIT dynamics for the inspiral still yields a significant speed-up of around \(1.5 - 7 \times\), with a median of \(\sim 4.5\) across the total masses examined.

Finally, we study how the speed-up varies with mass ratio for a system \(M = 60M_\odot\), started at periastron with initial orbit averaged frequency of \(10\)Hz and initial eccentricity of \(e_0 = 0.4\). In the top panel of Fig. 10, the EOB dynamics cost increases with mass ratio, while the cost of the NIT dynamics cost is independent of it. In the bottom panel, the total walltime of both the EOB and NIT waveforms increases with mass ratio. Again, this is because waveform generation, now the bottleneck, scales with signal length, which grows with mass ratio. Nonetheless, the NIT inspiral still gives a significant \(3 - 5.6\times\) speed-up with a median of \(\sim 5.3\times\) across the mass ratios examined. In both panels, we still see that the NIT timings are dominated by the final EOB transition to plunge.

8 Discussion↩︎

In this paper, we tackled the cost of eccentric EOB inspiral dynamics by recasting the equations of motion into quasi-Keplerian osculating orbital elements and then applying near-identity (averaging) transformations to obtain a new set of equations of motion that are independent of the orbital phases. These can be solved efficiently on the radiation-reaction timescale, decoupling the computational cost and the mass ratio. We combine this with an efficient waveform sampling method which further reduces computational cost. Finally, we transform back to the original EOB equations of motion to model the final transition to plunge. Comparing against waveforms generated in the standard way, we found that the NIT dynamics must be computed to at least second post-adiabatic (2PA) order to accurately model comparable mass systems. This yields mismatches \(\mathcal{M}_{22} \leq 8.05 \times 10^{-5}\) across the parameter space while reducing the dynamics cost by up to 2 orders of magnitude and achieving an overall waveform-generation speed-up of between \(\sim 1.5 - 8 \times\).

Several improvements remain possible. While the 2PA NIT waveforms are sufficiently accurate for current and future ground based detectors, the OOE waveforms suggest that higher postadiabatic orders would achieve even greater accuracy which might be necessary for the highest-SNR MBHB signals one may see with LISA. Additionally, we find that the highest mismatches occur sporadically across parameter space, most likely due to our transformation back to the EOB coordinates for the final transition to plunge. Moreover, most of the computational cost is due to this same region which uses the slower EOB dynamics. Both points warrant a more detailed study of optimal handover conditions. A more physically motivated eccentric transition to plunge could further improve the accuracy and speed of the resulting waveforms [150], [151].

Our current parametrization of the dynamics is ill-suited to the quasicircular limit. A follow-up paper will use a different set of osculating orbital elements that are better suited to this limit. For systems with moderate eccentricities that evolve to be less than our cutoff eccentricity of \(\tilde{e} < 0.035\), we will start the inspiral with the present parametrization and then transition to the new one for the remainder of the inspiral [134].

For long signals, we are now dominated by the waveform-generation cost. Combining the averaged dynamics with a transformation to monotonically evolving phases [152] would allow one to use a multi-voice decomposition of the waveform modes [153]. This can be trivially parallelized on both CPUs and GPUs to rapidly construct long inspiral signals as demonstrated by the FastEMRIWaveforms package [113]–[116]. Adapting this method to the variables presented here will also be explored in a follow-up paper [134].

Another approach would be to reduce the need for dense time-domain sampling to produce faithful frequency domain waveforms. Many ongoing works explore this, using nonuniform time sampling, frequency-domain transformations via the stationary phase approximation [154] or shifted uniform asymptotics [46], [155], or transforming the waveform to time-frequency-domain [156], [157]. All of these could be applied to our averaged dynamics to further reduce waveform-generation cost.

So far we have implemented only nonspinning dynamics and waveforms. Adding (anti-)aligned spin is a straightforward extension. The only difficulty is interpolating over the values of the primary and secondary spins, but efficient multidimensional interpolation methods are well established and should suffice [158], [159]. Adding spin precession is less trivial, but still possible, requiring transformations that average over both eccentric and spin precession oscillations [129], [152].

The averaged dynamics trivially relate to the eccentric dynamics for EMRIs, modeled with the gravitational self-force formalism. Many efforts have already used information from the test-mass limit to improve EOB waveforms at large mass ratios [143], [160]–[165]. Our framework can offer new calibration methods. For example, calibrating the test-mass limit of eccentric EOB radiation reaction force to eccentric self-force flux calculations [166] would greatly reduce the dominant error of EOB waveforms for extreme mass ratios stemming from missing adiabatic order information [159]. Conversely, one could supplement adiabatic self-force models with postadiabatic information from eccentric EOB models to produce eccentric hybrid models, similar to what has been done with post-Newtonian information [133], [167], [168].

Finally, the method presented here is agnostic to the specific EOB model, so long as the equations of motion take the form of Eqs. 3 . Having demonstrated it on a toy, inspiral-only EOB model, a natural next step would be to apply it to state-of-the-art eccentric EOB models, such as TEOBResumS-Dali [87]–[98], SEOBNRE [99]–[101], and SEOBNRv6EHM [104], [105]. This would enable efficient, accurate parameter estimation of eccentric BBHs for both current and future detectors. While more work remains, this method stands as an important first step toward efficient eccentric inspirals in the EOB framework.

9 Radial momentum expansion coefficients↩︎

To obtain an expression for \(p_{r_*}\) in terms of the orbital elements, we start from the EOB effective Hamiltonian given by Eq. 2 and rearrange it into a power series in \(p_{r_*}\): \[\begin{align} \hat{p}_{r*}^2 & = H_{\text{eff}}^2 - A \left(1 + \frac{p_\phi^2}{r^2} \right) \\& = p_{r_*}^2 + A Q_4 p_{r_*}^4 + A Q_6 p_{r_*}^6 + A Q_8 p_{r_*}^8 + \mathcal{O}(p_{r_*}^{10}), \end{align}\] where \(\hat{p}_{r*} = p_{r_*}\) when \(Q = 0\), and the potentials \(Q_4\), \(Q_6\), and \(Q_8\) are given explicitly in Appendix A of Ref. [83]. The quantity \(\hat{p}_{r*}\) is defined entirely in terms of quantities with explicit expressions in \((p,e,\zeta)\). We can thus perturbatively invert this expansion to obtain \(p_{r_*}\) as an approximate power series in \(\hat{p}_{r*}\): \[p_{r_*}(\nu,p,e,\zeta) \simeq \text{sgn}(\sin(\zeta)) \sum_{n=1}^{7} c_n(\nu,p,e,\zeta) \hat{p}_{r*}^n(\nu,p,e,\zeta),\] where the nonzero coefficients \(c_n\) are given by \[\begin{align} c_1 & = 1, \\ c_3 & = -\frac{1}{2} A Q_4, \\ c_5 & = \frac{1}{8} \left( 7 A^2 Q_4^2 - 4 A Q_6\right) \\ c_7 & = \frac{1}{16} \left( -33 A^3 Q_4^3 + 36 A^2 Q_4 Q_6 - 8 A Q_8 \right). \end{align}\]

10 Fourier coefficients of functions via their reciprocal↩︎

Here we summarize a method to compute the Fourier coefficients of a function \(A(\zeta)\) given those of its reciprocal \(A(\zeta)^{-1}\). This is particularly useful for the functions \(f_t(\zeta)\) and \(f_\phi(\zeta)\) which appear Eqs. 7 . In Fig. 11, we plot the absolute values of the Fourier coefficients of \(f_t(\zeta)\) and its reciprocal \(f_t(\zeta)^{-1}\) (normalized by the zeroth coefficient) with the parameters \((q,p_0,e_0) = (1,20,0.4)\). We see that even at moderate eccentricities, \(f_t(\zeta)\) would not be well approximated by a low-order Fourier expansion, while \(f_t(\zeta)^{-1}\) would. We now establish a method to exploit this observation.

Let the Fourier coefficients of a function \(A(\zeta)\) be denoted \(A_\kappa\) as in Eq. 9 . Let \(B(\zeta)\) be the reciprocal function \(A(\zeta)^{-1}\) with Fourier coefficients \(B_{\kappa'}\). Assuming \(A(\zeta) \neq 0\), so that \(A(\zeta) B(\zeta) = 1\), we can write: \[\sum_{\kappa} A_\kappa e^{i \kappa \zeta} \sum_{\kappa'} B_{\kappa'} e^{i \kappa' \zeta} = \sum_\kappa \sum_{\kappa'} A_{\kappa'} B_{\kappa - \kappa'} e^{i \kappa \zeta} = 1.\] Since the right-hand side is a constant that is independent of \(\zeta\), this must mean that the above sum is \(1\) for \(\kappa = 0\) and otherwise zero. Thus we have the following relationship between the Fourier coefficients of \(A\) and \(B\): \[\sum_{\kappa'} A_{\kappa'} B_{\kappa - \kappa'} = \delta_{\kappa,0},\] where \(\delta_{\kappa,0}\) is the Kronecker delta.

Note that we are restricted to a finite number of Fourier coefficients for both \(A\) and \(B\) such that \(-\kappa_\text{max} \leq \kappa \leq \kappa_\text{max}\) and \(-\kappa'_\text{max} \leq \kappa' \leq \kappa'_\text{max}\). We define the \(N = 2\kappa_\text{max} + 1\) dimensional column vector of Fourier coefficients for \(A\) as \(\mathbf{A} = (A_{-\kappa_\text{max}}, \ldots, A_{\kappa_\text{max}})^T\). From the coefficients of \(B_\kappa'\), we can construct the \(N \times N\) Toeplitz matrix \(\mathbf{G}\) as \[G_{i,j} = \begin{cases} B_{i - j} & |i-j| \leq \kappa_\text{max}' \\ 0 & |i-j| > \kappa_\text{max}'. \end{cases}\] Explicitly, the matrix \(\mathbf{G}\) has the following form: \[\mathbf{G} = \begin{bmatrix} B_0 & B_{-1} & B_{-2} & \cdots \\ B_1 & B_0 & B_{-1} & \cdots \\ B_2 & B_1 & B_0 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}.\]

This allows us to write the above relationship between the Fourier coefficients of \(A\) and \(B\) in matrix form as \(\mathbf{G} \mathbf{A} = \mathbf{e}_0\), where \(\mathbf{e}_0\) is a row vector of length \(N\) with a \(1\) at the index corresponding to the \(\kappa = 0\) entry and zeros everywhere else. Thus, so long as the matrix \(\mathbf{G}\) is invertible, we can solve for the Fourier coefficients of \(A\) given the Fourier coefficients of \(B\) via \(\mathbf{A} = \mathbf{e}_0 \mathbf{G}^{-1}\).

To demonstrate this method, we first compute the Fourier coefficients of \(f_t(\zeta)\) as outlined in Sec. 3.3 using \(\kappa_\text{max} = 30\). We then calculate the Fourier coefficients of \(f_t(\zeta)^{-1}\) using \(\kappa'_\text{max} = 10\) and then use the above method to compute the Fourier coefficients of \(f_t(\zeta)\) for \(-30 \le \kappa \le 30\). We show the results in Fig. 12. We plot the absolute value of the modes \(|f_{t,\kappa} |\) in the top panel. We see that the standard method reaches a numerical floor of \(\sim 10^{-9}\) at \(\kappa = \pm 18\). The reciprocal method, on the other hand, is able to compute Fourier coefficients up a numerical floor of \(\sim 10^{-15}\) which occurs at \(\kappa = \pm 28\), thus allowing us to resolve more of the harmonic content of \(f_t(\zeta)\). This is confirmed in the bottom panel where we plot the relative error between the original function and its Fourier series computed with each method. The \(L_2\) norm of the relative error is \(8.2 \times 10^{-10}\) for the standard method and \(7 \times 10^{-11}\) for the reciprocal method. Thus, by using only \(21\) Fourier coefficients of the reciprocal function, we are able to compute the \(61\) Fourier coefficients of the original function to higher accuracy than the standard method.

13 Second postadiabatic order mismatches for the inspiral only↩︎

In this section, we more closely examine the accuracy of the 2PA NIT model by plotting the mismatches for inspiral-only waveforms across parameter space, without handing over to the EOB dynamics to compute the transition to plunge.

In Fig. 13, we show the second post-adiabatic mismatches for a variety of eccentricities and mass ratios for four different total masses, with the same setup as in Fig. 5 but with the waveforms truncated when one of the three mismatch criteria is met.

This removes the errors that might occur when transforming from the NIT variables to the EOB variables when computing the transition to plunge. We see that the mismatches now vary smoothly across the parameter space, and are generally lower than in Fig. 5 with \(\mathcal{M}_{22} \lesssim 10^{-5}\). For all masses, we see slightly higher mismatches for the lowest eccentricities than at moderate eccentricities. We believe this is due to the poor convergence of the parametrization as one approaches the quasi-circular limit, especially for the most relativistic systems with \(M = 100M_\odot\).

Outside of this, we see that the mismatch generally follows the expected trend of getting worse for higher eccentricities and smaller mass ratios, corresponding to systems with lower adiabaticity.

References↩︎