# Fourier domain gravitational waveforms for precessing eccentric binaries

Antoine Klein, Yannick Boetzel, Achamveedu Gopakumar, Philippe Jetzer, Lorenzo de Vittori

► To cite this version:

Antoine Klein, Yannick Boetzel, Achamveedu Gopakumar, Philippe Jetzer, Lorenzo de Vittori. Fourier domain gravitational waveforms for precessing eccentric binaries. Physical Review D, 2018, 98 (10), pp.104043. 10.1103/PhysRevD.98.104043. hal-01707540

HAL Id: hal-01707540

https://hal.science/hal-01707540v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Fourier domain gravitational waveforms for precessing eccentric binaries

Antoine Klein,¹,²,³ Yannick Boetzel,⁴ Achamveedu Gopakumar,⁵ Philippe Jetzer,⁴ and Lorenzo de Vittori⁴

¹Department of Physics and Astronomy, The University of Mississippi, University, Mississippi 38677, USA

²CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa,
Avenida Rovisco Pais 1, 1049 Lisboa, Portugal

³CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Bd Arago, 75014 Paris, France

⁴Physik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland

⁵Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research,
Mumbai 400005, India

<figure><img src="image_2.png" /></figure>

(Received 25 January 2018; published 26 November 2018)

We build two families of inspiral waveforms for precessing binaries on eccentric orbits in the Fourier domain. To achieve this, we use a small eccentricity expansion of the waveform amplitudes in order to separate the periastron precession timescale from the orbital timescale, and use a shifted uniform asymptotics transformation to compute the Fourier transform in the presence of spin-induced precession. We show that the resulting waveforms can yield a median faithfulness above 0.993 when compared to an equivalent time domain waveform with an initial eccentricity of $e_0 \approx 0.3$. We also show that when the spins are large, using a circular waveform can potentially lead to significant biases in the recovery of the parameters, even when the system has fully circularized, particularly when the accumulated number of cycles is large. This is an effect of the residual eccentricity present when the objects forming the binary have nonvanishing spin components in the orbital plane.

DOI: 10.1103/PhysRevD.98.104043

# I. INTRODUCTION

The recent discoveries of gravitational wave (GW) signals by the LIGO and Virgo Collaborations have opened a new observation window on the Universe [1–9], through which the potential for new discoveries in astrophysics is truly tremendous. So far, those events have been analyzed with the assumption that the systems that produced them were evolving on circular orbits. Indeed, it has been a well-known fact that the emission of gravitational waves by binary systems has the tendency to circularize their orbits [10]. Nevertheless, it has been argued that certain astrophysical scenarios could lead to stellar-origin black hole binaries having high initial eccentricities [11–14], so they would still be measurable when the signal reaches the frequency window of the space-based GW detector LISA [15–17]. Furthermore, recent results have shown that eccentricity measurements by LIGO could be used to constrain stellar-mass black hole formation mechanisms [14,18–22]. It has been estimated that large biases in the recovery of the parameters of the first direct detection GW150914 could have occurred if the initial eccentricity in the detector was $e_0 \gtrsim 0.05$ [23]. Supermassive black hole binaries could also have important eccentricities in the late inspiral, if triple systems are a significant ingredient of supermassive black hole evolution [24–29]. Furthermore, in some spin configurations, it has been shown that the eccentricity of the system never truly vanishes, but reaches

a stationary value where it ceases to decrease through the emission of GWs [30].

This has motivated the development of waveforms that include the effects of a nonzero eccentricity in GW binary signals. The first steps towards this goal rely on the derivation of quasi-Keplerian equations describing the orbits [31], the derivation of the evolution equations for the orbital elements [32–36], and the derivation of GW polarization amplitudes [37]. The effects of individual spins were later added to this approach [30,38–42]. Using these solutions, several waveforms have been developed. Yunes et al. [43] proposed an analytic eccentric waveform in the post-Newtonian (PN) postcircular approximation, by solving for the Fourier phase of a binary signal analytically at Newtonian order using a small eccentricity expansion. Cornish and Key [44–46] and Gopakumar and Schäfer [47] independently developed a numerical waveform in the time domain by solving the 1.5PN equations of motion numerically, together with the spin-orbit precession equations, and using 1.5PN accurate amplitudes. Huerta et al. [48] expanded the analytical work of Yunes et al. by including the most important eccentricity-dependent terms up to 3.5PN order and at eighth order in the initial eccentricity for nonspinning systems. Tanay et al. [49] later computed the full 2PN Fourier phase for nonspinning systems at second order in the eccentricity. Moore et al. [50] then expanded this result to 3PN order. Huerta et al. [51,52] and Hinder et al. [53]
combined those results with numerical relativity to produce an eccentric inspiral-merger-ringdown waveform for nonspinning binaries. Recently, Hinderer and Babak [54] and Cao and Han [55] developed an eccentric waveform using a new approach in the effective one-body (EOB) formalism.

In this work, we further develop the formalism of post-Newtonian eccentric waveforms to include the effects of spin-induced precession in the Fourier domain. The advantage of Fourier domain waveforms over time domain ones is that they provide a much more computationally efficient way of computing a GW signal. Indeed, in order to produce a time domain waveform, one has to compute an equally spaced time series of the signal before computing its Fourier transform to use in detection or parameter estimation algorithms. The relevant timescale for this time series is the inverse of the maximum orbital frequency, which, being very short, makes this process computationally very expensive. Having a waveform available directly in the Fourier domain circumvents this problem and greatly reduces the computational cost of GW data analysis. In order to construct such a waveform, we solve the evolution equations for the orbital elements, together with the orbit-averaged spin precession equations numerically at 3PN order, including spin effects at 2PN order. Using a quasi-Keplerian description of the orbit, we employ instantaneous nonspinning amplitudes to construct the resulting GW polarizations. We then use a shifted uniform asymptotics (SUA) technique [56] to compute the waveforms in the Fourier domain. The resulting waveform has the advantage that the phasing is computed without any expansion for small eccentricities and thus can be very faithful compared to corresponding time domain waveforms for moderate to large eccentricities ($e \lesssim 0.4$). However, the amplitudes require a small-eccentricity expansion, and thus we do not expect the present waveforms to be faithful for arbitrarily large eccentricities.

In Sec. II, we derive two different families of eccentric waveforms. Due to the similarity between the orbital timescale and the periastron-to-periastron timescale, we derive the first family by expanding the Fourier domain waveform into combined harmonics of the mean orbital phase and of the mean anomaly. We then derive the second family by further expanding the resulting Fourier phase and time-frequency relations for small differences between the two similar phases. In Sec. III, we describe simulations that we performed to compute the faithfulness between our Fourier domain waveforms and a corresponding time domain waveform, including a detailed summary of how these different waveforms are constructed. We also compare a circular waveform to probe which domain of the parameter space allows for such circular waveforms to be effectively used for parameter estimation of binary signals. We give concluding remarks in Sec. IV. Throughout this paper, we use geometric units where $G = c = 1$.

## II. WAVEFORM

In the presence of spins, the orbit of a binary system is, in general, not restricted to an orbital plane [57]. Indeed, interactions between the spins and the orbit cause them to precess. However, in the post-Newtonian regime, the timescale on which this precession occurs is well separated from the other timescales present in the problem. We can therefore approximate the spin-orbital precession to be occurring much more slowly than the orbit, which allows us to describe it using a so-called quasi-Keplerian parametrization inside an orbital plane that stays perpendicular to the orbital angular momentum as the latter precesses. A quasi-Keplerian parametrization of the orbit of a spinning binary system is known at 3PN order for the nonspinning part [31,33] and at 2PN order for the spin-dependent part [30]. In this work, we restrict the quasi-Keplerian orbital description at 2PN for the computation of the polarization amplitudes. We can express the orbit at 2PN order as

$$r = a(1 - e_r \cos u) + f_r(v), \quad (1a)$$

$$\phi = (1 + k)v + f_{\phi}(v), \quad (1b)$$

$$\tan \frac{v}{2} = \sqrt{\frac{1 + e_{\phi}}{1 - e_{\phi}}} \tan \frac{u}{2}, \quad (1c)$$

$$l = u - e_t \sin u + f_t(u, v), \quad (1d)$$

$$\dot{l} = n, \quad (1e)$$

where $(r, \phi)$ is a polar representation of the separation vector in the orbital plane; $a$ is the semimajor axis; $u$ is the eccentric anomaly; $v$ is the true anomaly; $l$ is the mean anomaly; $n$ is the mean motion; $e_r$, $e_\phi$ and $e_t$ are eccentricity parameters; and the functions $f_i$ are general relativistic corrections given by [30–32]

$$f_r(v) = \sum_{i=0}^{2} b_{r,i} \cos(2v - 2\psi_i), \quad (2a)$$

$$f_{\phi}(v) = \sum_{k=2}^{3} a_{\phi,k} \sin(kv) + \sum_{k=1}^{2} \sum_{i=0}^{2} b_{\phi,k,i} \sin(kv - 2\psi_i), \quad (2b)$$

$$f_t(u, v) = g_t(v - u) + a_t \sin(v), \quad (2c)$$

where $\psi_i$ is the angle between the periastron line and the projection of spin $i$ onto the orbital plane (see Fig. 1), $\psi_0 = (\psi_1 + \psi_2)/2$, and the constants $a_A$, $b_A$ and $g_t$ are listed in Appendix B. We complemented the spinning solution of [30] by including quadrupole-monopole terms as described in Appendix A.
<figure><img src="image_3.png" /><figcaption>FIG. 1. Angles used in the definition of the relativistic corrections defined in Eq. (2). The orbital plane is perpendicular to the Newtonian orbital angular momentum $\mathbf{L}$, and the invariant plane is perpendicular to the conserved total angular momentum $\mathbf{J}$. The angle $\phi_e$ locates the periastron line with respect to the fixed invariant plane, and the angles $\psi_i$ are subtended by the periastron line and the projections of the spins onto the orbital plane.</figcaption></figure>

The orbital phase $\phi$ and the mean anomaly $l$ can be decomposed as the sum of a linearly growing part and a periodic part [32],

$$ \phi = \lambda + W_{\phi}, \qquad (3a) $$

$$ \dot{\lambda} = (1 + k)n, \qquad (3b) $$

$$ \dot{l} = n, \qquad (3c) $$

$$ W_{\phi} = (1 + k)(v - l) + f_{\phi}(v). \qquad (3d) $$

We choose to express our equations in terms of the post-Newtonian parameter $y$ and the eccentricity parameter $e$ defined by

$$ y = \frac{[M(1 + k)n]^{1/3}}{\sqrt{1 - e_t^2}}, \qquad (4a) $$

$$ e = e_t. \qquad (4b) $$

The constants in the quasi-Keplerian parametrization are given in terms of these parameters in Appendix B.

As the system emits gravitational waves, the orbital frequency and the eccentricity will evolve according to the following equations [30,32]:

$$ M \frac{dy}{dt} = (1 - e^2)^{3/2} \nu y^9 \left( a_0 + \sum_{n=2} a_n y^n \right), \qquad (5a) $$

$$ M \frac{de^2}{dt} = -(1 - e^2)^{3/2} \nu y^8 \left( b_0 + \sum_{n=2} b_n y^n \right), \qquad (5b) $$

where the constants $a_i$ and $b_i$ are given at 3PN order for nonspinning systems and at 2PN order for spinning systems in Appendix C. Here we also complemented the spinning solution of [30] by including quadrupole-monopole terms as described in Appendix A. We found that the minimum value for the eccentricity $e_{\min}^2$ found in [30] is unchanged by the addition of quadrupole-monopole effects, with

$$ \begin{aligned} e_{\min}^2 &= \frac{5y^4}{304} \sigma(-1, 1, 0, 2, -2, 0) \\ &= \frac{5y^4}{304} |\mathbf{s}_{\perp}^{(-)}|^2, \end{aligned} \qquad (6) $$

where the 2PN spin-spin coupling $\sigma$ can be found in Appendix A, and

$$ \mathbf{s}_{\perp}^{(-)} = [\mathbf{s}_1 - \hat{\mathbf{L}}(\hat{\mathbf{L}} \cdot \mathbf{s}_1)] - [\mathbf{s}_2 - \hat{\mathbf{L}}(\hat{\mathbf{L}} \cdot \mathbf{s}_2)], \qquad (7) $$

where $\hat{\mathbf{L}}$ is a normal to the orbital plane.

Note that we found a typo in [30], where the constant factor in $e_{\min}^2$ should read 5/304 instead of 5/340. This minimum eccentricity depends on the spin orientations: it is maximal when the spins lie inside the orbital plane and are opposite to one another, and it vanishes when the projections of $\mathbf{s}_1$ and $\mathbf{s}_2$ onto the orbital plane are equal to each other. The maximum value it can take is independent of the mass ratio; it is $e_{\min}^2 = 5y^4/304$, which evaluates to $e_{\min} \approx 0.021$ at the ISCO defined by $y = 6^{-1/2}$, and it is multiplied by a factor $(f/f_{\text{ISCO}})^{2/3}$ earlier in the inspiral. Note that this minimum eccentricity, being a spin effect, is unrelated to a similar effect observed in extreme mass-ratio inspirals around Schwarzschild black holes in [58], and also unrelated to another effect due to orbital effects derived in [59], which is of order $e_{\min}^2 \sim y^{10}$ and is independent of the spins.

The 2PN orbit-averaged equations of precession are given by [57,60]

$$ M \hat{\mathbf{L}} = -(1 - e^2)^{3/2} y^6 (\mathbf{\Omega}_1 + \mathbf{\Omega}_2), \qquad (8a) $$

$$ M \dot{\mathbf{s}}_1 = (1 - e^2)^{3/2} \mu_2 y^5 \mathbf{\Omega}_1, \qquad (8b) $$

$$ M \dot{\mathbf{s}}_2 = (1 - e^2)^{3/2} \mu_1 y^5 \mathbf{\Omega}_2, \qquad (8c) $$

where we defined the reduced spins $\mathbf{s}_i = \mathbf{S}_i/m_i$, the reduced individual masses $\mu_i = m_i/M$, and the precession vectors $\mathbf{\Omega}_i$ are given by

$$ \mathbf{\Omega}_i = \left\{ \left[ 2\mu_i + \frac{3}{2}\mu_j - \frac{3}{2}y\hat{\mathbf{L}} \cdot (\mathbf{s} + (q_i - 1)\mathbf{s}_i) \right] \hat{\mathbf{L}} + \frac{1}{2}y\mathbf{s}_j \right\} \times \mathbf{s}_i, \qquad (9) $$

where $i, j \in \{1, 2\}$, $i \neq j$, and the $q_i$ are quadrupole parameters defined in such a way that $q_i = 1$ for black holes.
The gravitational waveform emitted by a binary system on such an orbit has been computed at 3PN order for nonspinning binaries, omitting tail effects [37]. The result has the following structure:

$$h(t) = F_{+}(t)h_{+}(t) + F_{\times}(t)h_{\times}(t), \quad (10a)$$

$$h_{+, \times}(t) = \sum_{n \in \mathbb{Z}} H_{+, \times}^{(n)}(y, e, e \cos u, e \sin u) e^{in(\phi + \phi_T)}, \quad (10b)$$

where $F_+$ and $F_\times$ are antenna pattern functions [61], and the Thomas precession phase $\phi_T$ is given by [62]

$$\dot{\phi}_T = \frac{\hat{\mathbf{L}} \cdot \hat{\mathbf{N}}}{1 - (\hat{\mathbf{L}} \cdot \hat{\mathbf{N}})^2} (\hat{\mathbf{L}} \times \hat{\mathbf{N}}) \cdot \dot{\hat{\mathbf{L}}}, \quad (11)$$

with respect to a given sky location vector $\hat{\mathbf{N}}$.

In order to compute the Fourier transform of this signal, we need to separate the orbital timescale from the precessional one and express the orbital timescale dependence in terms of linearly growing phases. To do so, we follow [63] and compute an inversion of the PN-accurate Kepler equation (1d) as

$$u = l + \sum_{s=1}^{\infty} A_s \sin(sl), \quad (12)$$

with the Fourier coefficients $A_s$ given by

$$A_s = \frac{2}{s} J_s(se) + \sum_{j=1}^{\infty} \alpha_j [J_{s+j}(se) - J_{s-j}(se)]. \quad (13)$$

The PN-accurate constants $\alpha_j$ can be computed from [31] and are given in Eq. (18) of [63]. Similarly, we can find a Fourier expansion of the true anomaly $v$ and the orbital phase $\phi$ in terms of the mean anomaly $l$:

$$v = l + \sum_{s=1}^{\infty} B_s \sin(sl), \quad (14a)$$

$$\phi = \lambda + \sum_{s=1}^{\infty} C_s \sin(sl). \quad (14b)$$

The Fourier coefficients $A_s$, $B_s$ and $C_s$ can be found up to $O(y^4, e^5)$ in Appendix D. Using this solution, we can then express

$$e^{iku} = \sum_{p \in \mathbb{Z}} e^{ku_p} e^{-ipl}, \quad (15a)$$

$$e^{ikv} = \sum_{p \in \mathbb{Z}} e^{kv_p} e^{-ipl}, \quad (15b)$$

$$e^{ik\phi} = e^{ik\lambda} \sum_{s \in \mathbb{Z}} P_p^{k\phi} e^{-ipl}, \quad (15c)$$

where the coefficients $e^{ku_p}$, $e^{kv_p}$ and $P_p^{k\phi}$ are given as a Taylor expansion in both $e$ and $y$. We refer to Eqs. (30), (34), (E11) of [63] for how to calculate these Fourier coefficients.

This small eccentricity expansion allows us to express the waveform as

$$h_{+, \times}(t) = \sum_{n \in \mathbb{Z}} \sum_{p \in \mathbb{Z}} H_{+, \times}^{(p,n)} e^{-i(n\lambda + pl)}, \quad (16)$$

where we included the Thomas phase $\phi_T$ into the amplitudes $H_{+, \times}^{(p,n)}$, which vary on the spin-precession timescale. To separate the periastron precession timescale from the orbital timescale, we define $\delta\lambda = \lambda - l$, such that

$$M\dot{\lambda} = (1 - e^2)^{3/2}y^3, \quad (17a)$$

$$\delta\dot{\lambda} = \frac{k}{1+k}\dot{\lambda}. \quad (17b)$$

This new angle defines the periastron precession timescale, which is similar to the spin precession timescale since $\delta\dot{\lambda}/\dot{\lambda} = O(y^2)$.

Using this, we can then further simplify the waveform with

$$h_{+, \times}(t) = \sum_{n \in \mathbb{Z}} H_{+, \times}^{(n)} e^{-in\lambda}, \quad (18a)$$

$$H_{+, \times}^{(n)} = \frac{M\nu y^2}{d_L} \sum_{m \in \mathbb{Z}} G_{+, \times}^{(m,n)} e^{-im\delta\lambda} e^{-i(n+m)\phi_T}. \quad (18b)$$

The amplitudes $G_{+, \times}^{(m,n)}$ are given in Appendix E at order $O(y^2, e)$ [64].

## A. Fourier transform approximations

Before we compute an approximation of the Fourier transform of our signal, let us introduce two useful techniques.

Let us first assume that we have a signal of the form

$$h(t) = A(t)e^{-i\phi(t)}, \quad (19)$$

with $\dot{\phi}(t)$ a positive and monotonically increasing function of time, and that we want to compute its Fourier transform

$$\tilde{h}(f) = \int_{-\infty}^{\infty} h(t) e^{2\pi ift} dt. \quad (20)$$

The stationary phase approximation (SPA) of this Fourier transform consists in Taylor expanding the amplitude $A(t)$ and phase $\phi(t)$ around the stationary point $t_0$ defined by the relation

$$2\pi f = \dot{\phi}(t_0), \quad (21)$$
keeping only the constant term in the expansion of the amplitude and up to the quadratic order in the expansion of the phase. We get

$$h(t) \approx A(t_0) \exp \left\{ -i \left[ \phi(t_0) + \dot{\phi}(t_0)(t-t_0) + \frac{1}{2} \ddot{\phi}(t_0)(t-t_0)^2 \right] \right\}. \quad (22)$$

We can compute the Fourier transform of this approximate signal analytically, and we get

$$\tilde{h}(f) \approx \sqrt{\frac{2\pi}{|\dot{\phi}(t_0)|^2}} A(t_0) e^{i(2\pi ft_0 - \phi(t_0) - \pi/4)}. \quad (23)$$

This approximation will be accurate if $|\dot{A}/A| \ll |\ddot{\phi}|^{1/2}$ around the stationary point, and if the quadratic approximation is accurate around the stationary point. For a formal derivation, see e.g., [65].

Let us now suppose that our signal is of the form

$$h(t) = A e^{-i[\phi_C + B \sin \beta]}, \quad (24)$$

with $\dot{A}/A = \mathcal{O}(y^8)$, $\dot{\phi}_C = \mathcal{O}(y^3)$, $B = \mathcal{O}(y)$, $\dot{\beta} = \mathcal{O}(y^5)$, and that each additional time derivative adds a factor $\mathcal{O}(y^8)$ to the various quantities present in the signal, with y a small expansion parameter. This is the simplified form of a GW signal that we expect from a binary system undergoing spin-induced orbital precession, with y being a PN expansion parameter. The SPA cannot be directly used in this case, because the two terms in the second time derivative of the signal phase

$$\ddot{\phi} = \ddot{\phi}_C - B\dot{\beta}^2 \sin\beta + \mathcal{O}(y^{13}) \quad (25)$$

are of the same PN order and can cancel each other. The shifted uniform asymptotics (SUA) method [56] offers an approximation of the Fourier transform of such a signal by first expanding the signal using Bessel functions as

$$h(t) = A \sum_k J_k(B) e^{-i(\phi_C + k\beta)}, \quad (26)$$

so its Fourier transform can be approximated by a series of SPA, since $\ddot{\beta} \ll \ddot{\phi}_C$. The Fourier transform then becomes

$$\tilde{h}(f) \approx \sum_k A \sqrt{\frac{2\pi}{\ddot{\phi}_C + k\ddot{\beta}}} \times \exp \left[ i \left( 2\pi ft_k - \phi_C - k\beta - \frac{\pi}{4} \right) \right], \quad (27)$$

where the various functions are evaluated for each $k \in \mathbb{Z}$ at the stationary time $t_k$ defined by

$$2\pi f = \dot{\phi}_C(t_k) + k\dot{\beta}(t_k). \quad (28)$$

The different stationary times can be related to each other by Taylor expanding their defining equations around $t_0$ and solving for the difference order by order:

$$t_k - t_0 = - \frac{k\dot{\beta}(t_0)}{\ddot{\phi}_C(t_0)} + \mathcal{O}(y^{-4}). \quad (29)$$

By Taylor expanding Eq. (27) around $t_0$, and keeping only the leading PN order amplitude and the phase accurate to order $\mathcal{O}(y^0)$, we obtain

$$\tilde{h}(f) \approx \tilde{h}_0(f) \tilde{h}_{\text{corr}}(f), \quad (30a)$$

$$\tilde{h}_0(f) = \sqrt{\frac{2\pi}{|\dot{\phi}(t_0)|^2}} A(t_0) e^{i(2\pi ft_0 - \phi_C(t_0) - \pi/4)}, \quad (30b)$$

$$\tilde{h}_{\text{corr}}(f) = \sum_k J_k[B(t_0)] \exp \left[ -k\beta(t_0) + \frac{1}{2} T^2 k^2 \dot{\beta}^2(t_0) \right], \quad (30c)$$

$$T = \sqrt{\frac{1}{\ddot{\phi}_C(t_0)}}. \quad (30d)$$

After some manipulation, we can resum the Bessel functions in $\tilde{h}_{\text{corr}}(f)$ as

$$\tilde{h}_{\text{corr}}(f) = \sum_{p \ge 0} \frac{(-i T^2)^p}{2^p p!} \partial_t^{2p} e^{-iB \sin \beta}, \quad (31)$$

where the functions are evaluated at $t = t_0$. Truncating this series at some order $p = k_{\text{max}}$ and using a stencil around $t_0$ to approximate the different order time derivatives, we obtain the SUA approximation

$$\tilde{h}_{\text{corr}}(f) \approx \sum_{k=-k_{\text{max}}}^{k_{\text{max}}} a_{k,k_{\text{max}}} e^{-iB \sin \beta(t_0+kT)}, \quad (32)$$

where the constants $a_{k,k_{\text{max}}}$ satisfy the following linear system of equations:

$$\sum_{k=1}^{k_{\text{max}}} a_{k,k_{\text{max}}} + \frac{1}{2} a_{0,k_{\text{max}}} = 1, \quad (33a)$$

$$\sum_{k=1}^{k_{\text{max}}} a_{k,k_{\text{max}}} \frac{k^{2p}}{(2p)!} = \frac{(-i)^p}{2^p p!}, \quad p \in \{1, \dots, k_{\text{max}}\}, \quad (33b)$$

$$a_{-k,k_{\text{max}}} = a_{k,k_{\text{max}}}. \quad (33c)$$
To summarize, if we are able to separate the spin-precessional timescale effects from a carrier phase $\phi_C$ that satisfies $\dot{\phi}_C > 0$ and $\ddot{\phi}_C > 0$ as

$$h(t) = A(t)e^{-i\phi_C(t)}, \quad (34)$$

where all spin-precessional timescale effects are included in $A(t)$, then we can write the SUA approximation of its Fourier transform:

$$\tilde{h}(f) = \sqrt{\frac{2\pi}{\ddot{\phi}_C}} e^{i(2\pi ft_0 - \phi_C(t_0) - \pi/4)} \times \sum_{k=-k_{\max}}^{k_{\max}} a_{k,k_{\max}} A(t_0 + kT), \quad (35)$$

with the constants $a_{k,k_{\max}}$ satisfying the linear system of Eqs. (33), and $T = [\ddot{\phi}_C(t_0)]^{-1/2}$.

## B. Periastron precession effects

Let us first derive a waveform in the Fourier domain that is valid for nonprecessing spins, and add the effects of spin-orbit precession later. Putting aside spin-orbit precession, the signal in the time domain can be expressed as in Eqs. (18):

$$h_{+,x}(t) = \frac{M\nu y^2}{d_L} \sum_{n \in \mathbb{Z}} \sum_{m \in \mathbb{Z}} G_{+,x}^{(m,n)} e^{-i(n\lambda + m\delta\dot{\lambda})}. \quad (36)$$

Using the SPA, we can approximate its Fourier transform by

$$\begin{aligned} \tilde{h}_{+,x}(f) &= \int h_{+,x}(t)e^{2\pi i f t} dt \\ &= \frac{M\nu y^2}{d_L} \sum_{n \ge 1} \sum_{m \in \mathbb{Z}} \sqrt{\frac{2\pi}{n\ddot{\lambda} + m\delta\dot{\lambda}}} G_{+,x}^{(m,n)} \\ & \times e^{i[2\pi ft_{n,m} - n\lambda(t_{n,m}) - m\delta\dot{\lambda}(t_{n,m}) - \pi/4]}, \quad (37a) \\ 2\pi f &= n\dot{\lambda}(t_{n,m}) + m\delta\dot{\lambda}(t_{n,m}), \quad (37b) \end{aligned}$$

where each of the harmonics $(n, m)$ has to be evaluated at a different time. It is worth noting here that we assumed that $n\dot{\lambda} + m\delta\dot{\lambda} > 0$, which is not necessarily true for every $(n, m)$ pair during the whole inspiral. However, for this assumption to break down, one needs negative and sufficiently large $m$, since $\delta\dot{\lambda}/\dot{\lambda} = O(y^2)$, and the corresponding amplitude will be suppressed by a factor $e^m$. We verified that ignoring this fact does not lead to high inaccuracies, at least for initial eccentricities $e_0 \lesssim 0.4$.

In order to simplify the expression of the Fourier domain waveform and to improve its computational efficiency, we look for an expression of the following form:

$$\tilde{h}_{+,x}(f) = \frac{M\nu y^2}{d_L} \sum_{n \ge 1} \tilde{h}_{n,0}(f) \sum_{m \in \mathbb{Z}} \tilde{h}_{n,m}^{\text{PP}}(f), \quad (38a)$$

$$\tilde{h}_{n,0}(f) = \sqrt{\frac{2\pi}{n\ddot{\lambda}}} e^{i[2\pi ft_n - n\dot{\lambda}(t_n) - \pi/4]}, \quad (38b)$$

$$2\pi f = n\dot{\lambda}(t_n), \quad (38c)$$

where $\tilde{h}_{n,0}(f)$ is a waveform harmonic without any periastron precession effects and $\tilde{h}_{n,m}^{\text{PP}}(f)$ are corrections to it. In order to evaluate $\tilde{h}_{n,m}^{\text{PP}}(f)$, we define

$$\Delta t_{n,m} = t_{n,m} - t_n, \quad (39)$$

and Taylor expand Eq. (37b) around $t_n$:

$$2\pi f = \sum_{p \ge 0} \frac{\Delta t_{n,m}^p}{p!} \left. \frac{d^p}{dt^p} (n\dot{\lambda} + m\delta\dot{\lambda}) \right|_{t=t_n}. \quad (40)$$

We can use this together with Eq. (38c) to solve for the PN expansion of $\Delta t_{n,m}$ order by order, and we obtain

$$\Delta t_{n,m} = \sum_{p=1}^{P} \frac{1}{p!} \left( -\frac{m}{n} \right)^p D^{p-1} \left( \frac{\delta \dot{\lambda}^p}{\ddot{\lambda}} \right), \quad (41)$$

where the differential operator $D$ is given by

$$D = \frac{1}{\ddot{\lambda}} \frac{d}{dt}, \quad (42)$$

and every function of time is evaluated at $t = t_n$. We have checked that this expression remains valid at least up to $P = 6$.

Using this, we can then Taylor expand the phase in Eq. (37a) around $t_n$ to compute

$$\tilde{h}_{n,m}^{\text{PP}}(f) = G_{+,x}^{(m,n)} e^{i\Delta\Psi_{n,m}}, \quad (43a)$$

$$\Delta\Psi_{n,m} = -m\delta\dot{\lambda} + n \sum_{p=2}^{P+1} \frac{1}{p!} \left(-\frac{m}{n}\right)^p D^{p-2} \left(\frac{\delta\dot{\lambda}^p}{\ddot{\lambda}}\right), \quad (43b)$$

where all functions are once again evaluated at the stationary time $t_n$ defined by Eq. (38c), and we checked that the latter equation is valid at least up to $P = 6$. Equations (41) and (43b) are PN expansions in the sense that each increasing order in $m$ is multiplied by a factor of PN order $(\delta\dot{\lambda}/\ddot{\lambda})(d/d_t) = O(y^2)$, as both $\delta\dot{\lambda}$ and $\ddot{\lambda}$ evolve on the radiation reaction timescale. This implies that the formal expansion in $m$ in these two equations coincides with a PN expansion at order $2P$ beyond leading order.

## C. Complete waveform

We can now add spin precession by using a SUA transformation [56] instead of a SPA. We start by
noting that we can express the waveform in the time domain by

$$h(t) = \sum_{n,m} A_{n,m}(t) e^{-i(n\lambda+m\delta\lambda)}, \quad (44)$$

where all spin-precession timescale effects are included in the amplitudes

$$A_{n,m}(t) = \frac{M\nu y^2}{d_L} [F_+(t)G_+^{(n,m)}(t) + F_\times(t)G_\times^{(n,m)}(t)]e^{-i(n+m)\phi_r}. \quad (45)$$

This allows us to directly use a SUA transformation. If we restrict the amplitudes to $O(y^N, e^M)$, we then obtain

$$\tilde{h}(f) = \sum_{n=\max(1,2-N)}^{2+N} \sum_{m=-M}^{M} \tilde{h}_{n,m}(f), \quad (46a)$$

$$\tilde{h}_{n,m}(f) = \tilde{h}_{n,m}^{(0)}(f) \tilde{h}_{n,m}^{\text{SP}}(f), \quad (46b)$$

$$\tilde{h}_{n,m}^{(0)}(f) = \sqrt{2\pi T_{n,m}} \exp[i(2\pi f t_{n,m} - n\lambda(t_{n,m}) - m\delta\lambda(t_{n,m}) - \pi/4)], \quad (46c)$$

$$2\pi f = n\dot{\lambda}(t_{n,m}) + m\delta\dot{\lambda}(t_{n,m}), \quad (46d)$$

$$T_{n,m} = [n\ddot{\lambda}(t_{n,m}) + m\delta\ddot{\lambda}(t_{n,m})]^{-1/2}, \quad (46e)$$

$$h_{n,m}^{\text{SP}}(f) = \sum_{k=-k_{\max}}^{k_{\max}} a_{k,k_{\max}} A_{n,m}(t_{n,m} + kT_{n,m}), \quad (46f)$$

where the constants $a_{k,k_{\max}}$ satisfy the linear system of equations defined in Eq. (33). This waveform is in the Fourier domain and consistently includes the effects of spin-induced precession and periastron precession. As we will see in the next section, it allows for large matches with time domain waveforms with eccentricities $e \lesssim 0.3$ that we can consider as moderate in the modeling sense, because only the amplitudes, not the phasing, have been expanded for small eccentricities.

The waveform defined by Eq. (46a) suffers from the fact that it includes a double sum, and therefore, its computational cost rises quickly as the precision of the amplitudes increases. However, in order to increase its computational efficiency, we can use a similar strategy as described in the previous subsection and expand the waveform in powers of $\delta\dot{\lambda}/\dot{\lambda}$.

First, we can approximate the SUA timescale in Eq. (46e) by $T_{n,m} \approx T_n = T_{n,0}$. Next, we can use Eqs. (41) and (43b) to define $\Delta t_{n,m}$ and $\Delta \Psi_{n,m}$ at order P:

$$\Delta t_{n,m} = \sum_{p=1}^{P} \frac{1}{p!} \left(-\frac{m}{n}\right)^p D^{p-1} \left(\frac{\delta \dot{\lambda}^p}{\dot{\lambda}}\right), \quad (47a)$$

$$\Delta\Psi_{n,m} = -m\delta\lambda + n \sum_{p=2}^{P+1} \frac{1}{p!} \left(-\frac{m}{n}\right)^p D^{p-2} \left(\frac{\delta\dot{\lambda}^p}{\dot{\lambda}}\right). \quad (47b)$$

We can use Eqs. (5a) and (5b) together with the chain rule

$$\frac{d}{dt} f(y, e^2) = \frac{\partial f}{\partial y} \frac{dy}{dt} + \frac{\partial f}{\partial e^2} \frac{de^2}{dt}, \quad (48)$$

to get the necessary derivatives of $\lambda$ and $\delta\lambda$ as PN expanded functions. Thus, we can simplify the waveform as

$$\tilde{h}(f) = \sum_{n=\max(1,2-N)}^{2+N} \tilde{h}_n^{(0)}(f) \tilde{h}_n^{\text{PP}}(f), \quad (49a)$$

$$\tilde{h}_n^{(0)}(f) = \sqrt{2\pi T_n} \exp \left[ i \left( 2\pi f t_n - n\lambda(t_n) - \frac{\pi}{4} \right) \right], \quad (49b)$$

$$2\pi f = n\dot{\lambda}(t_n), \quad (49c)$$

$$T_n = [n\ddot{\lambda}(t_n)]^{-1/2}, \quad (49d)$$

$$\tilde{h}_n^{\text{PP}}(f) = \sum_{m=-M}^{M} e^{i\Delta\Psi_{n,m}} \\ \times \sum_{k=-k_{\max}}^{k_{\max}} a_{k,k_{\max}} A_{n,m}(t_n + \Delta t_{n,m} + kT_n). \quad (49e)$$

Equation (49a) presents a further expanded waveform, and can possibly be made more efficient than the one defined by Eq. (46a), especially for amplitudes of high ($N, M$) order. Thus, we get a family of Fourier domain gravitational waveforms for spin-precessing binaries on eccentric orbits characterized by the expansion orders ($P, N, M$).

### III. COMPARISONS

We have run different sets of simulations in order to probe under what circumstances our waveforms defined in Eqs. (46) and (49) are sufficiently faithful to equivalent waveforms obtained in the time domain. For all the waveforms used in our comparisons, we use nonspinning amplitudes at 2PN order omitting tail terms [37], and we use evolution equations for $y$ and $e^2$ at 3PN nonspinning order and 2PN spinning order, including tail terms, as described in Appendix C. For all Fourier domain waveforms, we use a SUA transformation as in [56] with $k_{\max} = 3$.

We use as a reference time domain waveform $h_R$ obtained in the following way:

(i) Equations (5a)–(8c) are solved numerically together with Eqs. (17a), (17b), and (11) in order to yield solutions for $y(t)$, $e^2(t)$, $\hat{\mathbf{L}}(t)$, $\mathbf{s}_1(t)$, $\mathbf{s}_2(t)$, $\lambda(t)$, $\phi_T(t)$, and $\delta\dot{\lambda}(t)$, from an initial time $t_0$ to a maximum time $t_{\max}$ defined by the ISCO-like condition $M\dot{\lambda}(t_{\max}) = 6^{-3/2}$.
(ii) Time is equally sampled between $t_0$ and $t_{max}$ using a sampling time $\Delta t = 2\pi/[24\dot{\lambda}(t_{max})]$, in order to ensure that the first 12 waveform harmonics fall below the Nyquist frequency. Equations (1b)–(1d) are solved at each step to get the orbital phase $\phi$ and the eccentric anomaly $u$. Equation (1d) is inverted numerically to yield $u(l = \dot{\lambda} - \delta\dot{\lambda}, e = e_i)$.

(iii) A waveform signal is constructed using Eqs. (10a) and (10b), and the solutions for $y(t)$, $e(t)$, $\phi(t)$, $\phi_T(t)$, and $u(t)$. The antenna pattern functions are chosen in the low-frequency limit, for a static detector [61]. The waveform amplitudes are included at 2PN order, with the omission of spin effects and tail terms.

(iv) A Tukey window is introduced in order to reduce spectral leakage, and a discrete Fourier transform of the signal is taken to yield the waveform in the Fourier domain.

We compare different waveforms to the reference one:

(i) A nonexpanded eccentric one (NE$_M$) defined by Eq. (46a) and $N = 4$, $M \in \{0, ..., 6\}$, i.e., with amplitudes at $N/2 = 2$PN order and amplitudes expanded at $M$th order in $e$.

(ii) An expanded eccentric one (EE$_{M,P}$) defined by Eq. (49a) and $N = 4$, $M \in \{0, ..., 6\}$, $P \in \{0, ..., 3\}$, i.e., with amplitudes at $N/2 = 2$ PN order, amplitudes expanded at $M$th order in $e$, and the waveform expanded at $P$th order in $\delta\dot{\lambda}/\dot{\lambda}$ as in Eqs. (47).

(iii) A circular one (C) with amplitudes at 2PN order, taken from [56].

Note that Eqs. (47) imply that the waveforms NE$_0$ and EE$_{0,P}$ are identical for any P.

In order to make our comparisons, we compute the faithfulness $F = \max \mathcal{M}$, defined by the match $\mathcal{M}$ maximized over some of the waveform parameters, with

$$ \mathcal{M} = \frac{(h, h_R)}{\sqrt{(h, h)(h_R, h_R)}}, \qquad (50a) $$

$$ (a, b) = \int_{f_{\min}}^{f_{\max}} \tilde{a}(f) \tilde{b}^*(f) df, \qquad (50b) $$

where we chose a white detector noise in order to make as few assumptions about the detector as possible. For the eccentric waveforms, since they use the same phasing as the reference one, we do not maximize over any parameters, while for the circular waveform, we maximize the match over time and orbital phase shifts to obtain the faithfulness. We compare the faithfulness obtained this way to a fiducial value of $F = 0.993$, corresponding to a faithfulness level at which we can estimate that the errors in the recovered parameters due to mismodeling are smaller than the statistical errors coming from the detector noise, for $D = 10$ intrinsic parameters and a signal-to-noise ratio (SNR) of $\rho = 25$ [66]. The relation between the faithfulness and the SNR at which

the mismodeling error becomes likely to exceed the statistical error in a GW detection is [66]

$$ F \approx 1 - \frac{D-1}{2\rho^2}. \qquad (51) $$

We ran two different sets of simulations: one to study systems in the late inspiral as observed by the LIGO/Virgo network and by LISA in the case of massive black hole binaries [denoted by (Xa)], and the other to study systems in the early inspiral as observed by LISA for stellar-origin black hole binaries [denoted by (Xb)]. We made six different runs in order to probe the faithfulness of our waveforms as a function of the starting eccentricity in different situations:

(I) We randomize the initial eccentricity with a log-flat distribution $10^{-5} < e_0 < 0.5$ and the spin magnitudes with a flat distribution $0 < \chi_i < 1$.

(II) We randomize the initial eccentricity with a log-flat distribution $10^{-5} < e_0 < 0.5$ and spin magnitudes with a flat distribution $0 < \chi_i < 0.1$.

(III) We randomize the initial eccentricity with a log-flat distribution $10^{-5} < e_0 < 0.5$ and spin magnitudes set to the maximum value $\chi_i = 1$.

(IV) We start with zero initial eccentricity and spin magnitudes with a flat distribution $0 < \chi_i < 1$.

(V) We start with zero initial eccentricity and spin magnitudes with a flat distribution $0 < \chi_i < 0.1$.

(VI) We start with zero initial eccentricity and spin magnitudes set to the maximum value $\chi_i = 1$.

We thus have twelve runs: (Ia)–(VIa) in the late inspiral case and (Ib)–(VIb) in the early inspiral case.

To get the binary parameters used in our runs, we randomize all vector directions with a flat distribution on the sphere. Since the distance does not affect the match $\mathcal{M}$ in Eq. (50), we fix it at some fiducial value. We randomize the initial orbital phase and the initial periastron-ascending node angle $\phi_e$ (see Fig. 1) with a flat distribution in $[0, 2\pi]$. Whenever the randomized initial eccentricity is lower than the minimal value given in Eq. (6), we set $e_0 = e_{\min}$. Note that cases (IV) to (VI) correspond to fully circularized binaries, but Eq. (6) prevents them from having truly zero eccentricity unless the reduced spins have exactly equal support in the orbital plane. In each late inspiral run, we start our simulations with an initial eccentricity $e_0$ and at a frequency $M\dot{\lambda}_{\text{start}} = 6^{-3/2}/10$, and stop at $M\dot{\lambda}_{\text{stop}} = 6^{-3/2}$. We also randomize the mass ratio $q = m_2/m_1$ with a log-flat distribution between 1 and 1/30, and use a fixed total mass $M = 100 M_{\odot}$, taking advantage of the white detector noise. In each early inspiral run, we randomize the two masses with a log-flat distribution $10 M_{\odot} < m_i < 100 M_{\odot}$. We then use the Newtonian time-frequency relation and the initial eccentricity to determine the starting frequency such that the system will evolve to have an orbital frequency of $f_{\text{end}} = 1$ Hz after $T = 4$ yr:
$$y_{\text{start}} = \left[ \frac{5M y_{\text{end}}^8}{5M + 32\nu T(1 - e_0^2)^{3/2}(8 + 7e_0^2)y_{\text{end}}^8} \right]^{1/8}, \quad (52a)$$

$$y_{\text{end}} = \left[ \frac{2\pi M f_{\text{end}}}{(1 - e_0^2)^{3/2}} \right]^{1/3}. \quad (52b)$$

We then let the system evolve and stop after four years, and set the maximum frequency $f_{\text{max}} = 1$ Hz in Eq. (50b).

## A. Late inspiral systems

We present in Fig. 2 the results from late inspiral run (Ia), with starting eccentricity $10^{-5} < e_0 < 0.5$ and

<figure><img src="image_4.png" /><figcaption>FIG. 2. Results from late inspiral run (Ia), with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $0 < \chi_i < 1$. On the top, we show median faithfulness as a function of the starting eccentricity for the circular waveform C, the expanded eccentric waveforms $EE_{6,p}$, with amplitudes at 6th order in the eccentricity and at $P$th order in $\delta\dot{\lambda}/\dot{\lambda}$, $0 \le P \le 3$ with increasing $P$, and the nonexpanded waveform with amplitudes at 6th order in the eccentricity $NE_6$. At the bottom, we show median faithfulness as a function of the starting eccentricity for the expanded eccentric waveforms $EE_{M,2}$, with amplitudes at $M$th order in the eccentricity, $0 \le M \le 6$ and increasing $M$, and at second order in $\delta\dot{\lambda}/\dot{\lambda}$. The left axis shows the unfaithfulness $1-F$, and the right axis shows the corresponding threshold SNR $\rho$ above which we can expect mismodeling errors to exceed the accuracy in a measurement. In both panels, the thin horizontal black line shows a fiducial faithfulness of 0.993 or a corresponding SNR of 25.</figcaption></figure>

spin magnitudes $0 < \chi_i < 1$. In it, we compare the mean faithfulness as a function of the initial eccentricity for different waveforms. The top panel shows a comparison between the results for the circular waveform C, the nonexpanded eccentric waveform $NE_6$, and the expanded eccentric waveforms $EE_{6,p}$, $0 \le P \le 3$, and the bottom panel shows a comparison between the expanded eccentric waveforms $EE_{M,2}$, $0 \le M \le 6$. We can see in the top panel that the circular waveform stays above the fiducial faithfulness only for initial eccentricities of $e_0 \lesssim 0.008$. Furthermore, the results for the expanded eccentric waveform become very close to the nonexpanded version starting at second order in $\delta\dot{\lambda}/\dot{\lambda}$, and leads to a faithfulness above the fiducial threshold for eccentricities below $e_0 \lesssim 0.3$. On the bottom panel, we can see the effects of the expansion of the waveform amplitudes for small eccentricities. We can see that the largest starting eccentricity for which the median faithfulness stays above the threshold increases with increasing order in the expansion. Furthermore, we can see that below a certain starting eccentricity depending on the specific order, increasing the expansion order has no effect on the faithfulness, as the errors due to this approximation become subdominant.

We present in Fig. 3 the results from late inspiral run (IIa), with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $0 < \chi_i < 0.1$. These results are very similar to

<figure><img src="image_5.png" /><figcaption>FIG. 3. Same as Fig. 2, for late inspiral run (IIa) with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $0 < \chi_i < 0.1$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 4. Same as Fig. 2, for late inspiral run (IIIa) with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $\chi_i = 1$.</figcaption></figure>

the results of run (Ia), but due to the reduced spin magnitudes, the starting eccentricities reach smaller values. On the top panel, we can see that below a starting eccentricity of $e_0 \lesssim 10^{-3}$, the loss of faithfulness using circular waveforms with respect to our eccentric models becomes negligible.

We present in Fig. 4 the results from late inspiral run (IIIa), with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $\chi_i = 1$. The results are similar to the ones shown in Figs. 2 and 3, but the increased magnitudes of the spins slightly reduce the performance of the circular waveform. Comparing this figure to Figs. 2 and 3, we can conclude that the value of the spin magnitudes has little effect on the faithfulness, other than on the limiting residual eccentricity.

We present in Fig. 5 the results from late inspiral run (IVa), with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 1$. We can see here an effect due to the residual eccentricity. Indeed, the circular waveform performs poorly in some cases, even when the binaries are fully circularized. In our simulations, 7% of the faithfulness for the circular waveform was below the threshold line, while virtually no faithfulness was found below it for waveforms that used eccentric phasing, even with the lowest order amplitudes. While this does not represent a

<figure><img src="image_7.png" /><figcaption>FIG. 5. Results from late inspiral run (IVa) with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 1$. The systems simulated here correspond to highly spinning fully circularized binaries. The blue line corresponds to the circular waveform C, the red line to the lowest-order expanded eccentric waveform EE$_{0,0}$, and the green line to the highly accurate expanded eccentric waveform EE$_{6,2}$. The bottom axis shows the unfaithfulness $1-F$, and the top axis shows the corresponding threshold SNR $\rho$ above which we can expect mismodeling errors to exceed the accuracy in a measurement. The thin vertical line corresponds to a fiducial faithfulness of 0.993 or a corresponding SNR of 25. Note that due to the eccentricity being taken into account in the phasing, even the lowest-order eccentric waveform EE$_{0,0}$ performs better than the circular waveform C.</figcaption></figure>

large proportion of binaries, this number will only increase when considering binaries with higher SNRs.

We present in Fig. 6 the results from late inspiral run (Va), with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 0.1$. Comparing with the results shown in Fig. 5, we can see that assuming lower spins prevents the circular waveforms from having faithfulness below the threshold line. Thus, eccentricity effects in the inspiral can

<figure><img src="image_8.png" /><figcaption>FIG. 6. Same as Fig. 5, for late inspiral run (Va) with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 0.1$. The systems simulated here correspond to slowly spinning fully circularized binaries.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Same as Fig. 5, for late inspiral run (VIa) with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $\chi_i = 1$. The systems simulated here correspond to maximally spinning fully circularized binaries.</figcaption></figure>

be safely ignored when only the last part of it is visible. This further shows that the starting eccentricity is the most important factor to influencing the accuracy of our waveforms in the late inspiral.

We present in Fig. 7 the results from late inspiral run (VIa), with starting eccentricity $e_0 = e_{min}$ and spin

magnitudes $\chi_i = 1$. The results here are similar to the ones shown in Fig. 5, but more pronounced. The proportion of binaries for which the circular waveform has a faithfulness lying below the threshold line increases to 25%, indicating that the inclusion of eccentricity effects might be important even for fully circularized binaries in the last stages of their inspiral when their spins are large.

## B. Early inspiral systems

We present in Fig. 8 the results from early inspiral run (Ib), with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $0 < \chi_i < 1$. We can see that, in this case, using circular waveforms will likely result in large biases even when the starting eccentricity is below $10^{-3}$. The large number of orbital cycles accumulated is such that the small difference in the frequency evolution induces very low faithfulness even for very low eccentricities. On the other hand, the eccentric waveforms perform better than in the late inspiral case. In the top panel, we can see that the low-order $EE_{6,0}$ waveform stays above the faithfulness threshold for $e_0 \lesssim 0.05$, while the high-order one $EE_{6,2}$ is above the threshold for the whole parameter space that we investigated. In the bottom panel, we can see that the waveform with circular amplitudes $EE_{0,2}$ stays above the

<figure><img src="image_10.png" /><figcaption>FIG. 8. Same as Fig. 2, for early inspiral run (Ib) with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $0 < \chi_i < 1$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. Same as Fig. 2, for early inspiral run (Ib) with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $0 < \chi_i < 0.1$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Same as Fig. 2, for early inspiral run (IIb) with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $\chi_i = 1$.</figcaption></figure>

threshold for $e_0 \lesssim 0.1$, while the waveforms $EE_{M,2}$, $M \ge 2$ do so for $e_0 \lesssim 0.3$.

We present in Fig. 9 the results from early inspiral run (IIb), with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $0 < \chi_i < 0.1$. We can see that, for initial eccentricities $e_0 \gtrsim 10^{-4}$, circular waveforms yield a faithfulness below $F = 0.9$. Thus, even if they are slowly spinning, the use of circular waveforms for parameter estimation for such binaries is likely to yield important biases. Using eccentric waveforms for early inspiral systems is therefore crucial in order to ensure accurate parameter recovery, even with initial eccentricities as low as $e_0 \sim 10^{-4}$. In the bottom panel, similarly to run (Ib), we can see that the waveform with circular amplitudes $EE_{0,2}$ stays above the threshold for $e_0 \lesssim 0.1$, while the waveforms $EE_{M,2}$, $M \ge 2$ do so for $e_0 \lesssim 0.3$.

We present in Fig. 10 the results from early inspiral run (IIb), with starting eccentricity $10^{-5} < e_0 < 0.5$ and spin magnitudes $\chi_i = 1$. While the results for the eccentric waveforms are similar to the ones shown in Figs. 8 and 9, the circular waveform never reached a median faithfulness above the threshold line above an initial eccentricity of $e_0 = 3 \times 10^{-5}$. This indicates that highly spinning systems in the early inspiral will require the use of an eccentric model irrespective of their initial eccentricity.

<figure><img src="image_13.png" /><figcaption>FIG. 11. Same as Fig. 5, for early inspiral run (IVb) with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 1$. The systems simulated here correspond to highly spinning fully circularized binaries.</figcaption></figure>

We present in Fig. 11 the results from early inspiral run (IVb), with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 1$. We can see that for these systems, including the eccentricity in the phasing is important, but the order used in other effects matters very little. Indeed, the faithfulness distributions for the two eccentric waveforms $EE_{0,0}$ and $EE_{6,2}$ are indistinguishable and have support almost exclusively above the faithfulness threshold, whereas the faithfulness distribution for the circular waveforms has 46% support below the threshold.

We present in Fig. 12 the results from early inspiral run (Vb), with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 0.1$. We can see that for these systems, circular waveforms have a faithfulness distribution almost identical to those of eccentric waveforms, indicating that when the spins are small and the binaries have fully

<figure><img src="image_14.png" /><figcaption>FIG. 12. Same as Fig. 5, for early inspiral run (Vb) with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $0 < \chi_i < 0.1$. The systems simulated here correspond to slowly spinning fully circularized binaries.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIG. 13. Same as Fig. 5, for early inspiral run (VIb) with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $\chi_i = 1$. The systems simulated here correspond to maximally spinning fully circularized binaries.</figcaption></figure>

circularized, the use of circular waveforms may be suffi-
cient for unbiased parameter estimation.

We present in Fig. 13 the results from early inspiral run (VIb), with starting eccentricity $e_0 = e_{min}$ and spin magnitudes $\chi_i = 1$. We observe in this figure that the faithfulness distribution for the circular waveforms has 94% support below the threshold line, indicating that for highly spinning binaries, the use of eccentric waveforms will be crucial for unbiased parameter estimation. However, the distributions for the two eccentric waveforms $EE_{0,0}$ and $EE_{6,2}$ are indistinguishable also in this case, indicating that the precision of the waveform amplitude is of little importance. Thus, eccentricity and spins will be important to include in the analysis of stellar-origin black hole binaries with LISA to account for the possibility of high spins, even if the binaries have fully circularized. However, using accurate amplitudes might be unnecessary for those sources.

Comparing the different results, we find that the median
faithfulness for each waveform is mainly influenced by the
initial eccentricity and the stage in the inspiral that they
find themselves in. We summarize in Table I the initial

<table><caption>TABLE I. For a few select waveforms, maximum initial eccentricity $e_0^{\text{max},1}$ for which the median faithfulness in the late inspiral runs stays above the faithfulness threshold, and the same for the early inspiral runs $e_0^{\text{max},e}$.</caption>
  <thead>
    <tr>
      <th>Waveform</th>
      <th>e<sub>0</sub><sup>max,1</sup></th>
      <th>e<sub>0</sub><sup>max,e</sup></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>C</td>
      <td>0.0078</td>
      <td>&lt;7 × 10<sup>-5</sup></td>
    </tr>
    <tr>
      <td>EE<sub>0,0</sub></td>
      <td>0.056</td>
      <td>0.086</td>
    </tr>
    <tr>
      <td>EE<sub>6,0</sub></td>
      <td>0.036</td>
      <td>0.04</td>
    </tr>
    <tr>
      <td>EE<sub>0,2</sub></td>
      <td>0.056</td>
      <td>0.086</td>
    </tr>
    <tr>
      <td>EE<sub>2,2</sub></td>
      <td>0.23</td>
      <td>0.32</td>
    </tr>
    <tr>
      <td>EE<sub>4,2</sub></td>
      <td>0.29</td>
      <td>&gt;0.4</td>
    </tr>
    <tr>
      <td>EE<sub>6,2</sub></td>
      <td>0.32</td>
      <td>&gt;0.4</td>
    </tr>
    <tr>
      <td>NE<sub>6</sub></td>
      <td>0.33</td>
      <td>&gt;0.4</td>
    </tr>
  </tbody>
</table>

eccentricity below which the median faithfulness falls
above the threshold line for a few of the waveforms
compared in our simulations. Interestingly, we find that
waveform EE<sub>0,0</sub> performs slightly better than waveform
EE<sub>6,0</sub>. We find the same to be true by comparing EE<sub>0,0</sub> to
any waveform EE<sub>M,0</sub> or EE<sub>M,1</sub> with M > 0. We thus remark
that in order for the inclusion of beyond-circular effects in
the amplitudes to increase the accuracy of the waveform,
one also needs to include periastron precession effects at
least at second order.

IV. CONCLUSION

We have constructed two families of Fourier domain
waveforms for spin-precessing binaries on eccentric orbits.
These include phasing at the third nonspinning post-
Newtonian order, including leading-order spin-orbit and
spin-spin interactions, such as instantaneous amplitudes at
second post-Newtonian order as small eccentricity expan-
sions. In this work, we have used amplitudes up to O(e⁶),
but the extension to higher orders in the eccentricity
would be trivial though lengthy. Through comparisons
with a complete time domain waveform at consistent
post-Newtonian order, we find that our new waveforms
faithfully reproduce their Fourier transform for initial eccen-
tricities up to e₀ ~ 0.3 for systems in the late inspiral, and at
least up to e₀ ~ 0.4 for systems in the early inspiral such as
stellar-origin black hole binaries as observed by LISA.

Comparing results, we find that using circular wave-
forms would likely lead to significant biases in parameter
recovery, even for fully circularized binaries with a signal-
to-noise ratio around 25, provided they are highly spinning.
Indeed, a 2PN spin effect prevents the eccentricity of a
binary system from vanishing completely unless the pro-
jections of the reduced spins in the orbital plane are exactly
equal to each other. We find that the use of circular
waveforms can cause biases if fully circularized systems
with large spin magnitudes and random orientations are
observed in the late inspiral, but not if the spin magnitudes
are small. This situation is made worse if binary systems are
observed in the early inspiral, and we expect large biases
with circular waveforms irrespective of the initial eccen-
tricity for highly spinning systems, even if they are fully
circularized. However, if the spins are sufficiently small and
the binaries have circularized below an eccentricity of 10⁻⁴
when the observations start, we expect the use of circular
waveforms to be appropriate for parameter estimation.
Overall, we expect circular waveforms to be safe to use
for parameter estimation in the late inspiral if the initial
eccentricity falls below 10⁻² and in the early inspiral when
it falls below 10⁻⁴, but we would recommend the use of
eccentric phasing in the waveform to describe highly
spinning systems, even if they have fully circularized.

Those waveforms provide a step towards the inclusion of
the eccentricity in gravitational wave data analysis such as
that performed by the LIGO/Virgo community. We argue
from the simulations described in this paper that the inclusion of spins and eccentricity might be of importance for reducing potential biases in the parameter recovery of binaries, even when they are fully circularized. While circular templates might be appropriate to describe slowly spinning systems, it can be important to include in the modeling of highly spinning systems. It is worth noting that the faithfulness measurements described in this work are not suitable to estimate the loss of events due to mismodelling, or the measurability of binary parameters, including the initial eccentricity. We leave those questions open for future work.

Some assumptions made in this work, particularly the neglect of orbital timescale effects in the spin-orbit precession dynamics, have to be more closely investigated. Furthermore, the inclusion of the merger and ringdown signals in our waveforms is also very important work for the future, and will have to be taken into account in the construction of waveform templates to use in current and future detectors. The waveform that we have presented in this work, while useful to describe inspiral-dominated signals such as stellar-origin black hole binaries in LISA or neutron star binaries in the LIGO/Virgo network, is inspiral-only and therefore cannot be used alone in the characterization of merger-dominated signals such as black hole binaries as observed by the LIGO/Virgo network.

## ACKNOWLEDGMENTS

We thank Katerina Chatziioannou and Eliu Huerta for useful comments. We thank our referees for useful comments. A. K. is supported by NSF CAREER Grant No. PHY-1055103, by FCT Contract No. IF/00797/2014/CP1214/CT0012 under the IF2014 Programme, and by H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. This work was supported by the Centre National d'Études Spatiales. Y. B. and L. d. V. are supported by the Swiss National Science Foundation. This work has made use of the Horizon Cluster, hosted by the Institut d'Astrophysique de Paris. We thank Stéphane Rouberol for smoothly running this cluster for us.

## APPENDIX A: QUADRUPOLE-MONOPOLE EFFECTS

The 2PN part of the quasi-Keplerian parametrization found in [30] is based upon the reduced Lagrangian

$$ \frac{\mathcal{L}}{\nu} = \frac{\dot{r}^2}{2} + \frac{1}{r} + \frac{1}{r^3} s_1 \cdot s_2 - \frac{3}{r^3} (\hat{\mathbf{r}} \cdot s_1)(\hat{\mathbf{r}} \cdot s_2), \quad (A1) $$

where the reduced spins $s_i = S_i/m_i$. The quadrupole-monopole part of the reduced Lagrangian is [42,67]

$$ \frac{\mathcal{L}_{\text{QM}}}{\nu} = \frac{1}{2r^3} \sum_{i=1}^{2} q_i [s_i^2 - 3(\hat{\mathbf{r}} \cdot s_i)^2], \quad (A2) $$

where the quadrupole parameter $q_i$ is defined in such a way that $q_i = 1$ for black holes. The total Lagrangian can then be written as

$$ \frac{\mathcal{L}}{\nu} = \frac{\dot{\mathbf{r}}^2}{2} + \frac{1}{r} + \frac{1}{2r^3} s^2 - \frac{3}{2r^3} (\hat{\mathbf{r}} \cdot s)^2 \\ + \frac{1}{2r^3} \sum_{i=1}^{2} (q_i - 1) [s_i^2 - 3(\hat{\mathbf{r}} \cdot s_i)^2], \quad (A3) $$

where $s = s_1 + s_2$.

Thus, a quasi-Keplerian description of the orbit including quadrupole-monopole terms can be found by adding the 2PN terms of [30], using the substitutions ($s_1 \to s/\sqrt{2}$, $s_2 \to s/\sqrt{2}$), ($s_1 \to s_1\sqrt{(q_1-1)/2}$, $s_2 \to s_1\sqrt{(q_1-1)/2}$), and ($s_1 \to s_2\sqrt{(q_2-1)/2}$, $s_2 \to s_2\sqrt{(q_2-1)/2}$). It reads

$$ r = a(1 - e_r \cos u) + f_r(v), \quad (A4a) $$

$$ \phi = (1 + k)v + f_{\phi,1}(v) + f_{\phi,2}(v), \quad (A4b) $$

$$ \tan \frac{v}{2} = \sqrt{\frac{1 + e_{\phi}}{1 - e_{\phi}}} \tan \frac{u}{2}, \quad (A4c) $$

$$ l = u - e \sin u, \quad (A4d) $$

$$ \dot{l} = n, \quad (A4e) $$

with

$$ a = \frac{1}{(1 - e^2)y^2} \left[ 1 + \frac{1}{2}(1 + e^2)\gamma_1 y^4 \right], \quad (A5a) $$

$$ e_r^2 = e^2[1 + (1 - e^2)\gamma_1 y^4], \quad (A5b) $$

$$ k = \frac{3}{2}\gamma_1 y^4, \quad (A5c) $$

$$ e_{\phi}^2 = e^2[1 + 2(1 - e^2)\gamma_1 y^4], \quad (A5d) $$

$$ n = (1 - e^2)^{3/2} y^3 \left(1 - \frac{3}{2}\gamma_1 y^4\right), \quad (A5e) $$

$$ \gamma_1 = \frac{1}{2} \left\{ 3(\hat{\mathbf{L}} \cdot \mathbf{s})^2 - \mathbf{s}^2 + \sum_{i=1}^{2} (q_i - 1) [3(\hat{\mathbf{L}} \cdot \mathbf{s}_i)^2 - s_i^2] \right\}, \quad (A5f) $$

$$ f_r(v) = -\frac{y^2}{4} \left[ |\hat{\mathbf{L}} \times \mathbf{s}|^2 \cos(2v - 2\psi) + \sum_{i=1}^{2} (q_i - 1) |\hat{\mathbf{L}} \times \mathbf{s}_i|^2 \cos(2v - 2\psi_i) \right], \quad (A5g) $$

$$ f_{\phi,1}(v) = -\frac{y^4 e}{2} \left[ |\hat{\mathbf{L}} \times \mathbf{s}|^2 \sin(v - 2\psi) + \sum_{i=1}^{2} (q_i - 1) |\hat{\mathbf{L}} \times \mathbf{s}_i|^2 \sin(v - 2\psi_i) \right], \quad (A5h) $$
$$f_{\phi,2}(v) = -\frac{y^4}{8} \left[ |\hat{\mathbf{L}} \times \mathbf{s}|^2 \sin(2v - 2\psi) + \sum_{i=1}^{2} (q_i - 1) |\hat{\mathbf{L}} \times \mathbf{s}_i|^2 \sin(2v - 2\psi_i) \right], \quad (\text{A5i})$$

where $\psi$ is the angle subtended by the total reduced spin $\mathbf{s}$ and the periastron line, $\psi_i$ is the angle

subtended by the individual reduced spin $\mathbf{s}_i$ and the periastron line (see Fig. 1), and the periastron line is defined by the equation $v = u = l = 2p\pi$, $p \in \mathbb{Z}$. We can then use this representation of the orbit together with the orbit averaged evolution equations for the energy and orbital angular momentum computed in [40] to find

$$M \frac{dy}{dt} = \nu(1 - e^2)^{3/2} \left\{ \left( \frac{32}{5} + \frac{28}{5} e^2 \right) y^9 + \sigma \left[ -\frac{84}{5} - \frac{228}{5} e^2 - \frac{33}{5} e^4, \frac{242}{5} + \frac{654}{5} e^2 + \frac{381}{20} e^4, -\frac{447}{10} e^2 - \frac{93}{10} e^4, \right. \right. \\ \left. \left. \frac{88}{5} - 16q + \left( 48 - \frac{216}{5}q \right) e^2 + \left( \frac{69}{10} - \frac{63}{10}q \right) e^4, -\frac{244}{5} + 48q + \left( -132 + \frac{648}{5}q \right) e^2 + \left( -\frac{96}{5} + \frac{189}{10}q \right) e^4, \right. \right. \\ \left. \left. (1-q) \left( \frac{447}{10} e^2 + \frac{93}{10} e^4 \right) \right] y^{13} \right\}, \quad (\text{A6a})$$

$$M \frac{de^2}{dt} = -\nu(1 - e^2)^{3/2} \left\{ \left( \frac{608}{15} e^2 + \frac{242}{15} e^4 \right) y^8 + \sigma \left[ \frac{2}{3} - \frac{1961}{15} e^2 - \frac{2527}{12} e^4 - \frac{157}{8} e^6, -\frac{2}{3} + \frac{5623}{15} e^2 + \frac{2393}{4} e^4 \right. \right. \\ + \frac{447}{8} e^6, -\frac{5527}{30} e^2 - \frac{10117}{30} e^4 - \frac{5507}{160} e^6, -\frac{4}{3} + \left( \frac{682}{5} - \frac{1876}{15} q \right) e^2 + \left( \frac{1337}{6} - \frac{595}{3} q \right) e^4 + \left( \frac{83}{4} - \frac{37}{2} q \right) e^6, \\ \left. \left. \frac{4}{3} + \left( -\frac{5618}{15} + \frac{1876}{5} q \right) e^2 + \left( -\frac{1203}{2} + 595q \right) e^4 + \left( -\frac{225}{4} + \frac{111}{4} q \right) e^6, \right. \right. \\ \left. \left. \left( \frac{2764}{15} - \frac{921}{5} q \right) e^2 + \left( \frac{1687}{5} - \frac{5056}{15} q \right) e^4 + \left( \frac{551}{16} - \frac{172}{5} q \right) e^6 \right] y^{12} \right\}, \quad (\text{A6b})$$

where

$$\sigma(a, b, c, a_1 + a_2 q, b_1 + b_2 q, c_1 + c_2 q) = as^2 + b(\hat{\mathbf{L}} \cdot \mathbf{s})^2 + c|\hat{\mathbf{L}} \times \mathbf{s}|^2 \cos 2\psi + \sum_{i=1}^{2} \left[ (a_1 + a_2 q_i)\mathbf{s}_i^2 + (b_1 + b_2 q_i)(\hat{\mathbf{L}} \cdot \mathbf{s}_i)^2 + (c_1 + c_2 q_i)|\hat{\mathbf{L}} \times \mathbf{s}_i|^2 \cos 2\psi_i \right]. \quad (\text{A7})$$

We thus find the residual eccentricity found in [30] unchanged by quadrupole-monopole effects.

## APPENDIX B: QUASI-KEPLERIAN PARAMETRIZATION

A full quasi-Keplerian parametrization of the orbit at 2PN order in harmonic coordinates is [30,31]

$$r = a(1 - e_r \cos u) + f_r(v), \quad (\text{B1a})$$

$$\phi = (1 + k)v + f_{\phi}(v), \quad (\text{B1b})$$

$$\tan \frac{v}{2} = \sqrt{\frac{1 + e_{\phi}}{1 - e_{\phi}}} \tan \frac{u}{2}, \quad (\text{B1c})$$

$$l = u - e \sin u + f_l(u, v), \quad (\text{B1d})$$

$$\dot{i} = n, \quad (\text{B1e})$$
with

$$a = \frac{1}{(1 - e^2)y^2} \left\{ \begin{aligned} & 1 + \left[-1 + \frac{\nu}{3} + \left(3 - \frac{\nu}{3}\right)e^2\right]y^2 + \beta\left(\frac{2}{3} + 2e^2, 1 + e^2\right)y^3 + \left[5 + \frac{11}{4}\nu + \frac{\nu^2}{9}\right. \\ & \left. + \left(\frac{21}{2} - \frac{73}{6}\nu - \frac{2}{9}\nu^2\right)e^2 + \left(1 + \frac{5}{12}\nu + \frac{\nu^2}{9}\right)e^4 + (1 - e^2)^{3/2}(-5 + 2\nu) + \frac{\gamma_1}{2}(1 + e^2)\right]y^4, \end{aligned} \right. \quad (B2a)$$

$$e_r^2 = e^2 \left\{ 1 + (1 - e^2) \left\{ (8 - 3\nu)y^2 + \beta(4, 2)y^3 + \left[ 32 - \frac{467}{12}\nu + 4\nu^2 + \left(-40 + \frac{371}{12}\nu - 4\nu^2\right)e^2 + \sqrt{1 - e^2}(15 - 6\nu) + \gamma_1 \right] y^4 \right\} \right\}, \quad (B2b)$$

$$k = 3y^2 + \beta(4, 3)y^3 + \left[\frac{27}{2} - 7\nu + \left(\frac{51}{4} - \frac{13}{2}\nu\right)e^2 + \frac{3}{2}\gamma_1\right]y^4, \quad (B2c)$$

$$e_{\phi}^2 = e^2 \left\{ 1 + (1 - e^2) \left\{ (8 - 2\nu)y^2 + \beta(4, 4)y^3 + \left[ 42 - \frac{113}{12}\nu + \frac{11}{12}\nu^2 + \left(-40 + \frac{1043}{48}\nu - \frac{89}{48}\nu^2\right)e^2 + \sqrt{1 - e^2}(15 - 6\nu) + 2\gamma_1 \right] y^4 \right\} \right\}, \quad (B2d)$$

$$n = (1 - e^2)^{3/2} y^3 \left\{ 1 - 3y^2 - \beta(4, 3)y^3 + \left[ -\frac{9}{2} + 7\nu + \left(-\frac{51}{4} + \frac{13}{2}\nu\right)e^2 - \frac{3}{2}\gamma_1 \right] y^4 \right\}, \quad (B2e)$$

where

$$\beta(a, b) = -[(a\mu_1 + b\mu_2)s_1 + (b\mu_1 + a\mu_2)s_2] \cdot \hat{\mathbf{L}}, \quad (B3a)$$

$$\gamma_1 = \frac{1}{2} \left\{ 3(\hat{\mathbf{L}} \cdot \mathbf{s})^2 - \mathbf{s}^2 + \sum_{i=1}^{2} (q_i - 1)[3(\hat{\mathbf{L}} \cdot \mathbf{s}_i)^2 - \mathbf{s}_i^2] \right\}. \quad (B3b)$$

The functions $f_r$, $f_\phi$, $f_t$, and $f_n$ are given by

$$f_r(v) = \sum_{i=0}^{2} b_{r,i} \cos(2v - 2\psi_i), \quad (B4a)$$

$$f_{\phi}(v) = \sum_{k=2}^{3} a_{\phi,k} \sin(kv) + \sum_{k=1}^{2} \sum_{i=0}^{2} b_{\phi,k,i} \sin(kv - 2\psi_i), \quad (B4b)$$

$$f_t(u, v) = g_t(u - v) + a_t \sin(v), \quad (B4c)$$

with

$$b_{r,i} = -\frac{y^2}{4} F_i |\hat{\mathbf{L}} \times \mathbf{s}_i|^2, \quad (B5a)$$

$$a_{\phi,2} = e^2 \left(\frac{\nu}{8} - \frac{3}{8}\nu^2\right) y^4, \quad (B5b)$$

$$a_{\phi,3} = -e^3 \frac{3}{32} \nu^2 y^4, \quad (B5c)$$
$$b_{\phi,1,i} = -\frac{e}{2} F_i |\hat{L} \times \mathbf{s}_i|^2 y^4, \qquad (B5d)$$

$$b_{\phi,2,i} = -\frac{1}{8} F_i |\hat{L} \times \mathbf{s}_i|^2 y^4, \qquad (B5e)$$

$$g_t = (1 - e^2)^{3/2} \left( \frac{15}{2} - 3\nu \right) y^4, \qquad (B5f)$$

$$a_t = e(1 - e^2)^{3/2} \left( -\frac{\nu}{2} - \frac{\nu^2}{8} \right) y^4, \qquad (B5g)$$

where we defined, for convenience, $s_0 = s$, $\psi_0 = \psi$, $F_0 = 1$, $F_1 = q_1 - 1$, and $F_2 = q_2 - 1$.

## APPENDIX C: EVOLUTION EQUATIONS

The evolution equations of y and e are given at 3PN order by [30,34–36]

$$M \frac{dy}{dt} = (1 - e^2)^{3/2} \nu y^9 \left( a_0 + \sum_{n=2}^{6} a_n y^n \right), \qquad (C1a)$$

$$M \frac{de^2}{dt} = -(1 - e^2)^{3/2} \nu y^8 \left( b_0 + \sum_{n=2}^{6} b_n y^n \right), \qquad (C1b)$$

where

$$a_0 = \frac{32}{5} + \frac{28}{5} e^2, \qquad (C2a)$$

$$a_2 = -\frac{1486}{105} - \frac{88}{5}\nu + \left(\frac{12296}{105} - \frac{5258}{45}\nu\right)e^2 + \left(\frac{3007}{84} - \frac{244}{9}\nu\right)e^4, \qquad (C2b)$$

$$a_3 = \frac{128\pi}{5}\phi_y + \beta\left(\frac{904}{15} + \frac{2224}{15}e^2 + \frac{99}{5}e^4, 40 + \frac{1916}{15}e^2 + \frac{314}{15}e^4\right), \qquad (C2c)$$

$$\begin{aligned} a_4 = & \frac{34103}{2835} + \frac{13661}{315}\nu + \frac{944}{45}\nu^2 + \left(-\frac{256723}{945} - \frac{173587}{315}\nu + \frac{147443}{270}\nu^2\right)e^2 + \left(\frac{2095517}{7560} - \frac{589507}{504}\nu + \frac{34679}{45}\nu^2\right)e^4 \\ & + \left(\frac{53881}{2520} - \frac{7357}{90}\nu + \frac{9392}{135}\nu^2\right)e^6 + \frac{e^2}{1 - e^2}\left(\frac{85}{6} + \frac{1445}{6}\nu\right) + \frac{1 - \sqrt{1 - e^2}}{\sqrt{1 - e^2}}\left[16 - \frac{32}{5}\nu + \left(266 - \frac{532}{5}\nu\right)e^2\right. \\ & \left. + \left(-\frac{859}{2} + \frac{859}{5}\nu\right)e^4 + (-65 + 26\nu)e^6\right] + \sigma\left[-\frac{84}{5} - \frac{228}{5}e^2 - \frac{33}{5}e^4, \frac{242}{5} + \frac{654}{5}e^2 + \frac{381}{20}e^4, -\frac{447}{10}e^2 - \frac{93}{10}e^4, \right. \\ & \left. \frac{88}{5} - 16q + \left(48 - \frac{216}{5}q\right)e^2 + \left(\frac{69}{10} - \frac{63}{10}q\right)e^4, -\frac{244}{5} + 48q + \left(-132 + \frac{648}{5}q\right)e^2 + \left(-\frac{96}{5} + \frac{189}{10}q\right)e^4, \right. \\ & \left. (1 - q)\left(\frac{447}{10}e^2 + \frac{93}{10}e^4\right)\right], \qquad (C2d) \end{aligned}$$

$$a_5 = \pi \left( -\frac{4159}{105}\psi_y - \frac{756}{5}\nu\zeta_y \right), \qquad (C2e)$$
$$a_6 = \frac{16447322263}{21829500} - \frac{54784}{525} \gamma_E + \frac{512}{15} \pi^2 + \left(-\frac{56198689}{34020} + \frac{902}{15} \pi^2\right) \nu + \frac{541}{140} \nu^2 - \frac{1121}{81} \nu^3 + \left[\frac{247611308999}{87318000} - \frac{392048}{525} \gamma_E + \frac{3664}{15} \pi^2 + \left(-\frac{2828420479}{680400} + \frac{477}{4} \pi^2\right) \nu + \frac{1070903}{315} \nu^2 - \frac{392945}{324} \nu^3\right] e^2 + \left[-\frac{236637777001}{58212000} - \frac{93304}{175} \gamma_E + \frac{872}{5} \pi^2 + \left(\frac{2963572847}{453600} - \frac{53131}{960} \pi^2\right) \nu + \frac{44123941}{6048} \nu^2 - \frac{2198212}{405} \nu^3\right] e^4 + \left[-\frac{28913792717}{6468000} - \frac{4922}{175} \gamma_E + \frac{46}{5} \pi^2 + \left(\frac{107275139}{30240} - \frac{369}{80} \pi^2\right) \nu + \frac{5155951}{1512} \nu^2 - \frac{44338}{15} \nu^3\right] e^6 + \left(-\frac{243511057}{887040} + \frac{4179523}{15120} \nu + \frac{83701}{3780} \nu^2 - \frac{1876}{15} \nu^3\right) e^8 + \frac{e^2}{1-e^2} \left[\frac{91284763}{378000} + \left(\frac{19505077}{5040} - \frac{595}{8} \pi^2\right) \nu - \frac{48569}{12} \nu^2 - \frac{730168}{23625(1+\sqrt{1-e^2})}\right] + \frac{1-\sqrt{1-e^2}}{\sqrt{1-e^2}} \left\{-\frac{1425319}{3375} + \left(\frac{9874}{315} - \frac{41}{30} \pi^2\right) \nu + \frac{632}{15} \nu^2 + \left[\frac{2385427}{1050} + \left(-\frac{274234}{45} + \frac{4223}{240} \pi^2\right) \nu + \frac{70946}{45} \nu^2\right] e^2 + \left[\frac{8364697}{4200} + \left(\frac{1900517}{630} - \frac{32267}{960} \pi^2\right) \nu - \frac{47443}{90} \nu^2\right] e^4 + \left[-\frac{167385119}{25200} + \left(\frac{4272491}{504} - \frac{123}{160} \pi^2\right) \nu - \frac{43607}{18} \nu^2\right] e^6 + \left(-\frac{65279}{168} + \frac{510361}{1260} \nu - \frac{5623}{45} \nu^2\right) e^8\right\} + \frac{1284}{175} \kappa_y + \left(\frac{54784}{525} + \frac{392048}{525} e^2 + \frac{93304}{175} e^4 + \frac{4922}{175} e^6\right) \log\left[\frac{1+\sqrt{1-e^2}}{8y(1-e^2)^{3/2}}\right], \quad (C2f)$$

$$b_0 = \frac{608}{15} e^2 + \frac{242}{15} e^4, \quad (C2g)$$

$$b_2 = \left(-\frac{1878}{35} - \frac{8168}{45} \nu\right) e^2 + \left(\frac{59834}{105} - \frac{7753}{15} \nu\right) e^4 + \left(\frac{13929}{140} - \frac{3328}{45} \nu\right) e^6, \quad (C2h)$$

$$b_3 = \frac{788\pi e^2}{3} \phi_e + \beta \left( \frac{19688}{45} e^2 + \frac{28256}{45} e^4 + \frac{263}{5} e^6, \frac{1448}{5} e^2 + \frac{1618}{3} e^4 + \frac{167}{3} e^6 \right), \quad (C2i)$$

$$b_4 = \left(-\frac{952397}{945} + \frac{5937}{7}\nu + \frac{1504}{5}\nu^2\right) e^2 + \left(-\frac{3113989}{1260} - \frac{388419}{140}\nu + \frac{64433}{20}\nu^2\right) e^4 + \left(\frac{4656611}{1512} - \frac{13057267}{2520}\nu + \frac{127411}{45}\nu^2\right) e^6 + \left(\frac{420727}{1680} - \frac{362071}{1260}\nu + \frac{1642}{9}\nu^2\right) e^8 + \sqrt{1-e^2} \left[\left(\frac{2672}{3} - \frac{5344}{15}\nu\right) e^2 + \left(2321 - \frac{4642}{5}\nu\right) e^4 + \left(\frac{565}{3} - \frac{226}{3}\nu\right) e^6\right] + \sigma \left[\frac{2}{3} - \frac{1961}{15}e^2 - \frac{2527}{12}e^4 - \frac{157}{8}e^6, -\frac{2}{3} + \frac{5623}{15}e^2 + \frac{2393}{4}e^4 + \frac{447}{8}e^6, -\frac{5527}{30}e^2 - \frac{10117}{30}e^4 - \frac{5507}{160}e^6, -\frac{4}{3} + \left(\frac{682}{5} - \frac{1876}{15}q\right) e^2 + \left(\frac{1337}{6} - \frac{595}{3}q\right) e^4 + \left(\frac{83}{4} - \frac{37}{2}q\right) e^6, \frac{4}{3} + \left(-\frac{5618}{15} + \frac{1876}{5}q\right) e^2 + \left(-\frac{1203}{2} + 595q\right) e^4 + \left(-\frac{225}{4} + \frac{111}{4}q\right) e^6, \left(\frac{2764}{15} - \frac{921}{5}q\right) e^2 + \left(\frac{1687}{5} - \frac{5056}{15}q\right) e^4 + \left(\frac{551}{16} - \frac{172}{5}q\right) e^6\right], \quad (C2j)$$

$$b_5 = \pi \left( -\frac{55691}{105} \psi_e - \frac{610144}{315} \nu \zeta_e \right) e^2, \quad (C2k)$$
$$
\begin{aligned}
b_6 = & \left[ \frac{61655211971}{4365900} - \frac{2633056}{1575} \gamma_E + \frac{24608}{45} \pi^2 + \left( \frac{43386337}{56700} + \frac{1017}{5} \pi^2 \right) \nu - \frac{4148897}{1260} \nu^2 - \frac{61001}{243} \nu^3 \right] e^2 \\
& + \left[ \frac{64020009407}{21829500} - \frac{9525568}{1575} \gamma_E + \frac{89024}{45} \pi^2 + \left( \frac{770214901}{12600} - \frac{15727}{96} \pi^2 \right) \nu - \frac{80915371}{7560} \nu^2 - \frac{86910509}{9720} \nu^3 \right] e^4 \\
& + \left[ -\frac{1167012417073}{58212000} - \frac{4588588}{1575} \gamma_E + \frac{42884}{45} \pi^2 + \left( \frac{8799500893}{453600} - \frac{295559}{960} \pi^2 \right) \nu + \frac{351962207}{10080} \nu^2 - \frac{2223241}{90} \nu^3 \right] e^6 \\
& + \left[ \frac{120660628321}{12936000} - \frac{20437}{175} \gamma_E + \frac{191}{5} \pi^2 + \left( -\frac{91818931}{5040} - \frac{6519}{320} \pi^2 \right) \nu + \frac{2495471}{126} \nu^2 - \frac{11792069}{1215} \nu^3 \right] e^8 \\
& + \left( \frac{302322169}{887040} - \frac{1921387}{5040} \nu + \frac{41179}{108} \nu^2 - \frac{386792}{1215} \nu^3 \right) e^{10} + \sqrt{1-e^2} \left\{ \left[ -\frac{22713049}{7875} + \left( -\frac{11053982}{945} + \frac{8323}{90} \pi^2 \right) \nu \right. \right. \\
& \left. \left. + \frac{108664}{45} \nu^2 \right] e^2 + \left[ \frac{178791374}{7875} + \left( -\frac{38295557}{630} + \frac{94177}{480} \pi^2 \right) \nu + \frac{681989}{45} \nu^2 \right] e^4 + \left[ \frac{5321445613}{189000} \right. \right. \\
& \left. \left. + \left( -\frac{26478311}{756} + \frac{2501}{1440} \pi^2 \right) \nu + \frac{450212}{45} \nu^2 \right] e^6 + \left[ \frac{186961}{168} - \frac{289691}{252} \nu + \frac{3197}{9} \nu^2 \right] e^8 \right\} + \frac{1460336}{23625} (1-\sqrt{1-e^2}) \\
& + \frac{428}{1575} e^2 \kappa_e + \left( \frac{2633056}{1575} e^2 + \frac{9525568}{1575} e^4 + \frac{4588588}{1575} e^6 + \frac{20437}{175} e^8 \right) \log \left[ \frac{1+\sqrt{1-e^2}}{8y(1-e^2)^{3/2}} \right],
\end{aligned}
\quad (C21)
$$

with the tail terms given, in terms of the functions found in [34,36], by

$$ \phi_y = (1 - e^2)^{7/2} \tilde{\phi}, \quad (C3a) $$

$$ \phi_e = \frac{192(1 - e^2)^{9/2}}{985e^2} (\sqrt{1 - e^2}\phi - \tilde{\phi}), \quad (C3b) $$

$$ \psi_y = (1 - e^2)^{9/2} \left( -\frac{8064}{4159} \sqrt{1 - e^2} \phi + \frac{4032}{4159} \tilde{\phi} + \frac{8191}{4159} \tilde{\psi} \right), \quad (C3c) $$

$$ \zeta_y = (1 - e^2)^{7/2} \left[ \frac{160(1 - e^2)^{3/2}}{567} \phi + \left( -\frac{176}{567} + \frac{80}{567} e^2 \right) \tilde{\phi} + \frac{583(1 - e^2)}{567} \tilde{\zeta} \right], \quad (C3d) $$

$$ \psi_e = \frac{16382(1 - e^2)^{9/2}}{55691e^2} \left[ \left( \frac{9408}{8191} - \frac{14784}{8191} e^2 \right) \sqrt{1 - e^2} \phi + \left( -\frac{9408}{8191} + \frac{4032}{8191} e^2 \right) \tilde{\phi} + (1 - e^2)(\sqrt{1 - e^2}\psi - \tilde{\psi}) \right], \quad (C3e) $$

$$ \zeta_e = \frac{12243(1 - e^2)^{9/2}}{76268e^2} \left[ -\frac{16(1 - e^2)^{3/2}}{53} \phi + \left( \frac{16}{53} - \frac{80}{583} e^2 \right) \tilde{\phi} + (1 - e^2)(\sqrt{1 - e^2}\zeta - \tilde{\zeta}) \right], \quad (C3f) $$

$$ \kappa_y = -\frac{934088(1 - e^2)^5}{33705}(\tilde{\kappa} - \tilde{F}), \quad (C3g) $$

$$ \kappa_e = -\frac{5604528(1 - e^2)^6}{3745e^2} \left[ \sqrt{1 - e^2}(\kappa - F) - (\tilde{\kappa} - \tilde{F}) \right]. \quad (C3h) $$

We chose to only include in the 3PN enhancement functions $\kappa_i$ the terms proportional to $\log n$, as there is only a finite number of other terms, and they can be combined with nontail terms. Using the formalism developed in [34,36], we give them here at tenth order in the eccentricity:

$$ \phi_y = 1 + \frac{97}{32} e^2 + \frac{49}{128} e^4 - \frac{49}{18432} e^6 - \frac{109}{147456} e^8 - \frac{2567}{58982400} e^{10} + \mathcal{O}(e^{12}), \quad (C4a) $$

$$ \phi_e = 1 + \frac{5969}{3940} e^2 + \frac{24217}{189120} e^4 + \frac{623}{4538880} e^6 - \frac{96811}{363110400} e^8 - \frac{5971}{4357324800} e^{10} + \mathcal{O}(e^{12}), \quad (C4b) $$
$$ \psi_y = 1 - \frac{207671}{8318} e^2 - \frac{8382869}{266176} e^4 - \frac{8437609}{4791168} e^6 + \frac{10075915}{306634752} e^8 - \frac{38077159}{15331737600} e^{10} + \mathcal{O}(e^{12}), \qquad (C4c) $$

$$ \zeta_y = 1 + \frac{113002}{11907} e^2 + \frac{6035543}{762048} e^4 + \frac{253177}{571536} e^6 - \frac{850489}{877879296} e^8 - \frac{1888651}{10973491200} e^{10} + \mathcal{O}(e^{12}), \qquad (C4d) $$

$$ \psi_e = 1 - \frac{9904271}{891056} e^2 - \frac{101704075}{10692672} e^4 - \frac{217413779}{513248256} e^6 + \frac{35703577}{6843310080} e^8 - \frac{3311197679}{9854366515200} e^{10} + \mathcal{O}(e^{12}), \qquad (C4e) $$

$$ \zeta_e = 1 + \frac{11228233}{2440576} e^2 + \frac{37095275}{14643456} e^4 + \frac{151238443}{1405771776} e^6 - \frac{118111}{611205120} e^8 - \frac{407523451}{26990818099300} e^{10} + \mathcal{O}(e^{12}), \qquad (C4f) $$

$$ \kappa_y = 244 \log 2 \left( e^2 - \frac{18881}{1098} e^4 + \frac{6159821}{39528} e^6 - \frac{16811095}{19764} e^8 + \frac{446132351}{123525} e^{10} \right) - 243 \log 3 \left( e^2 - \frac{39}{4} e^4 + \frac{2735}{64} e^6 \right. \\ \left. + \frac{25959}{512} e^8 - \frac{638032239}{409600} e^{10} \right) - \frac{48828125 \log 5}{5184} \left( e^6 - \frac{83}{8} e^8 + \frac{12637}{256} e^{10} \right) - \frac{4747561509943 \log 7}{33177600} e^{10} + \mathcal{O}(e^{12}), \qquad (C4g) $$

$$ \kappa_e = 6536 \log 2 \left( 1 - \frac{22314}{817} e^2 + \frac{7170067}{19608} e^4 - \frac{10943033}{4128} e^6 + \frac{230370959}{15480} e^8 - \frac{866124466133}{8823600} e^{10} \right) \\ - 6561 \log 3 \left( 1 - \frac{49}{4} e^2 + \frac{4369}{64} e^4 + \frac{214449}{512} e^6 - \frac{623830739}{81920} e^8 + \frac{76513915569}{1638400} e^{10} \right) \\ - \frac{48828125 \log 5}{64} \left( e^4 - \frac{293}{24} e^6 + \frac{159007}{2304} e^8 - \frac{6631171}{27648} e^{10} \right) - \frac{4747561509943 \log 7}{245760} \left( e^8 - \frac{259}{20} e^{10} \right) + \mathcal{O}(e^{12}). \qquad (C4h) $$

It can be noted that those enhancement functions converge much more quickly than the ones presented in [34,36]. Indeed, because of the inclusion of factors of $\sqrt{1 - e^2}$ in them, the enhancement functions seem to converge in the parabolic limit $e \to 1$. We believe it to be related to the fact that the PN parameter $y$ we used here is related to the Newtonian orbital angular momentum and thus is finite and nonzero in this limit. In contrast, the PN parameter $(M\omega)^{1/3}$ is related to the energy and thus vanishes in this limit. In that case, in order for the tail effects to stay nonzero, the enhancement functions are forced to diverge.

## APPENDIX D: TRUE AND ECCENTRIC ANOMALY EXPANSION

The Fourier coefficients of the eccentric anomaly, true anomaly and orbital phase are given to order $\mathcal{O}(y^4, e^5)$ by

$$ A_1 = e - \frac{e^3}{8} + \frac{e^5}{192} + y^4 \left( e^3 \left( \frac{105}{8} - \frac{51\eta}{64} - \frac{19\eta^2}{64} \right) + e \left( -\frac{15}{2} + \frac{9\eta}{8} + \frac{\eta^2}{8} \right) + e^5 \left( -\frac{735}{128} - \frac{489\eta}{512} + \frac{111\eta^2}{512} \right) \right), $$

$$ A_2 = \frac{e^2}{2} - \frac{e^4}{6} + y^4 \left( e^4 \left( \frac{75}{4} - \frac{5\eta}{32} - \frac{47\eta^2}{96} \right) + e^2 \left( -\frac{75}{8} + \frac{15\eta}{16} + \frac{3\eta^2}{16} \right) \right), $$

$$ A_3 = \frac{3e^3}{8} - \frac{27e^5}{128} + y^4 \left( e^5 \left( \frac{6825}{256} + \frac{705\eta}{1024} - \frac{775\eta^2}{1024} \right) + e^3 \left( -\frac{95}{8} + \frac{49\eta}{64} + \frac{17\eta^2}{64} \right) \right), $$

$$ A_4 = \frac{e^4}{3} + e^4 y^4 \left( -\frac{975}{64} + \frac{35\eta}{64} + \frac{71\eta^2}{192} \right), $$

$$ A_5 = \frac{125e^5}{384} + e^5 y^4 \left( -\frac{5049}{256} + \frac{1167\eta}{5120} + \frac{523\eta^2}{1024} \right), \qquad (D1a) $$
$$B_1 = 2e^{-\frac{e^3}{4}} + \frac{5e^5}{96} + y^2 \left( e(4-\eta) + e^5 \left( \frac{17}{48} - \frac{17\eta}{192} \right) + e^3 \left( -\frac{7}{2} + \frac{7\eta}{8} \right) \right) \\ + y^4 \left( e^3 \left( -\frac{39}{4} + \frac{197\eta}{32} - \frac{7\eta^2}{32} \right) + e \left( 13 - \frac{31\eta}{12} + \frac{\eta^2}{12} \right) + e^5 \left( -\frac{401}{192} - \frac{569\eta}{144} + \frac{107\eta^2}{576} \right) \right),$$

$$B_2 = \frac{5e^2}{4} - \frac{11e^4}{24} + y^2 \left( e^2(4-\eta) + e^4 \left( -\frac{14}{3} + \frac{7\eta}{6} \right) \right) \\ + y^4 \left( e^4 \left( -\frac{217}{24} + \frac{3187\eta}{288} - \frac{313\eta^2}{288} \right) + e^2 \left( \frac{91}{8} - \frac{101\eta}{24} + \frac{11\eta^2}{24} \right) \right),$$

$$B_3 = \frac{13e^3}{12} - \frac{43e^5}{64} + y^2 \left( e^3 \left( \frac{9}{2} - \frac{9\eta}{8} \right) + e^5 \left( -\frac{207}{32} + \frac{207\eta}{128} \right) \right) \\ + y^4 \left( e^5 \left( -\frac{1011}{128} + \frac{9471\eta}{512} - \frac{1161\eta^2}{512} \right) + e^3 \left( \frac{43}{4} - \frac{51\eta}{8} + \frac{7\eta^2}{8} \right) \right),$$

$$B_4 = \frac{103e^4}{96} + e^4 y^2 \left( \frac{16}{3} - \frac{4\eta}{3} \right) + e^4 y^4 \left( \frac{1969}{192} - \frac{2711\eta}{288} + \frac{205\eta^2}{144} \right),$$

$$B_5 = \frac{1097e^5}{960} + e^5 y^2 \left( \frac{625}{96} - \frac{625\eta}{384} \right) + e^5 y^4 \left( \frac{3641}{384} - \frac{315821\eta}{23040} + \frac{10039\eta^2}{4608} \right), \qquad (\text{D1b})$$

$$C_1 = 2e^{-\frac{e^3}{4}} + \frac{5e^5}{96} + y^2 \left( e(10-\eta) + e^5 \left( \frac{49}{96} - \frac{17\eta}{192} \right) + e^3 \left( -\frac{17}{4} + \frac{7\eta}{8} \right) \right) \\ + y^4 \left( e^5 \left( -\frac{317}{96} + \frac{2497\eta}{2304} - \frac{1579\eta^2}{2304} \right) + e \left( 52 - \frac{235\eta}{12} + \frac{\eta^2}{12} \right) + e^3 \left( \frac{13}{8} - \frac{231\eta}{32} + \frac{17\eta^2}{32} \right) \right)$$

$$C_2 = \frac{5e^2}{4} - \frac{11e^4}{24} + y^2 \left( e^2 \left( \frac{31}{4} - \eta \right) + e^4 \left( -\frac{145}{24} + \frac{7\eta}{6} \right) \right) \\ + y^4 \left( e^2 \left( \frac{323}{8} - \frac{163\eta}{12} + \frac{\eta^2}{12} \right) + e^4 \left( -\frac{331}{24} + \frac{\eta}{18} + \frac{25\eta^2}{36} \right) \right)$$

$$C_3 = \frac{13e^3}{12} - \frac{43e^5}{64} + y^2 \left( e^3 \left( \frac{31}{4} - \frac{9\eta}{8} \right) + e^5 \left( -\frac{543}{64} + \frac{207\eta}{128} \right) \right) \\ + y^4 \left( e^3 \left( \frac{313}{8} - \frac{1205\eta}{96} + \frac{\eta^2}{32} \right) + e^5 \left( -\frac{749}{32} + \frac{7217\eta}{1536} + \frac{567\eta^2}{512} \right) \right)$$

$$C_4 = \frac{103e^4}{96} + e^4 y^2 \left( \frac{821}{96} - \frac{4\eta}{3} \right) + e^4 y^4 \left( \frac{1975}{48} - \frac{118\eta}{9} - \frac{11\eta^2}{144} \right)$$

$$C_5 = \frac{1097e^5}{960} + e^5 y^2 \left( \frac{9541}{960} - \frac{625\eta}{384} \right) + e^5 y^4 \left( \frac{10813}{240} - \frac{338987\eta}{23040} - \frac{1211\eta^2}{4608} \right). \qquad (\text{D1c})$$

## APPENDIX E: WAVEFORM AMPLITUDES EXPANSION

The amplitudes $G_{+,x}^{(m,n)}$ in Eq. (18b) are given to order $O(y^2, e)$ for $n < 0$, with $C_i = \cos i = \hat{L} \cdot \hat{N}$ and $S_i = \sin i$, by

$$G_{+}^{(1,-1)} = e \left( \frac{1}{2} - \frac{C_i^2}{2} \right) + ey^2 \left( -\frac{35}{16} + \frac{9C_i^2}{4} - \frac{C_i^4}{16} + \eta \left( -\frac{5}{48} - \frac{C_i^2}{12} + \frac{3C_i^4}{16} \right) \right),$$

$$G_{+}^{(1,-2)} = ey\delta \left( -\frac{3S_i}{2} + \frac{1}{2}C_i^2S_i \right),$$

$$G_{+}^{(1,-3)} = e \left( -\frac{9}{4} - \frac{9C_i^2}{4} \right) + ey^2 \left( -\frac{9}{8} - \frac{63C_i^2}{8} + \eta \left( -\frac{75}{8} + \frac{87C_i^2}{8} \right) \right),$$

$$G_{+}^{(1,-4)} = ey\delta(4S_i + 4C_i^2S_i),$$
$$G_{+}^{(1,-5)} = ey^2 \left( -\frac{625S_i^2}{96} - \frac{625}{96}C_i^2S_i^2 + \eta \left( \frac{625S_i^2}{32} + \frac{625}{32}C_i^2S_i^2 \right) \right),$$

$$G_{+}^{(0,-1)} = y\delta\left(-\frac{5S_i}{8} - \frac{1}{8}C_i^2S_i\right),$$

$$G_{+}^{(0,-2)} = -1 - C_i^2 + y^2 \left( \frac{19}{6} + \frac{3C_i^2}{2} - \frac{C_i^4}{3} + \eta \left( -\frac{19}{6} + \frac{11C_i^2}{6} + C_i^4 \right) \right),$$

$$G_{+}^{(0,-3)} = y\delta\left(\frac{9S_i}{8} + \frac{9}{8}C_i^2S_i\right),$$

$$G_{+}^{(0,-4)} = y^2 \left( -\frac{4S_i^2}{3} - \frac{4}{3}C_i^2S_i^2 + \eta(4S_i^2 + 4C_i^2S_i^2) \right),$$

$$G_{+}^{(-1,-1)} = e\left(\frac{3}{4} + \frac{3C_i^2}{4}\right) + ey^2\left(\frac{143}{24} + \frac{53C_i^2}{8} + \frac{C_i^4}{12} + \eta\left(\frac{7}{8} - \frac{9C_i^2}{8} - \frac{C_i^4}{4}\right)\right),$$

$$G_{+}^{(-1,-2)} = ey\delta\left(-\frac{3S_i}{2} - \frac{3}{2}C_i^2S_i\right),$$

$$G_{+}^{(-1,-3)} = ey^2 \left( \frac{81S_i^2}{32} + \frac{81}{32}C_i^2S_i^2 + \eta \left( -\frac{243S_i^2}{32} - \frac{243}{32}C_i^2S_i^2 \right) \right), \qquad (E1a)$$

$$G_{\times}^{(1,-1)} = iey^2 \left( -\frac{1}{4}C_iS_i^2 + \frac{3}{4}\eta C_iS_i^2 \right),$$

$$G_{\times}^{(1,-2)} = iey\delta C_iS_i,$$

$$G_{\times}^{(1,-3)} = i \left( \frac{9eC_i}{2} + ey^2 \left( 9C_i - \frac{27}{8}C_iS_i^2 + \eta \left( -\frac{3C_i}{2} + \frac{81}{8}C_iS_i^2 \right) \right) \right),$$

$$G_{\times}^{(1,-4)} = -8i\epsilon\gamma\delta C_iS_i,$$

$$G_{\times}^{(1,-5)} = iey^2 \left( \frac{625}{48} C_i S_i^2 - \frac{625}{16} \eta C_i S_i^2 \right),$$

$$G_{\times}^{(0,-1)} = \frac{3}{4}i\gamma\delta C_iS_i,$$

$$G_{\times}^{(0,-2)} = i \left( 2C_i + y^2 \left( -\frac{13C_i}{3} - \frac{4}{3}C_iS_i^2 + \eta \left( \frac{C_i}{3} + 4C_iS_i^2 \right) \right) \right),$$

$$G_{\times}^{(0,-3)} = -\frac{9}{4}i\gamma\delta C_iS_i,$$

$$G_{\times}^{(0,-4)} = i\gamma^2 \left( \frac{8}{3}C_iS_i^2 - 8\eta C_iS_i^2 \right),$$

$$G_{\times}^{(-1,-1)} = i \left( -\frac{3eC_i}{2} + ey^2 \left( -\frac{38C_i}{3} + \frac{11}{24}C_iS_i^2 + \eta \left( \frac{C_i}{2} - \frac{11}{8}C_iS_i^2 \right) \right) \right),$$

$$G_{\times}^{(-1,-2)} = 3i\epsilon\gamma\delta C_iS_i,$$

$$G_{\times}^{(-1,-3)} = i\epsilon\gamma^2 \left( -\frac{81}{16}C_iS_i^2 + \frac{243}{16}\eta C_iS_i^2 \right). \qquad (E1b)$$

Note that $G_{+,x}^{(m,-n)} = \bar{G}_{+,x}^{(-m,n)}$.
[1] J. Aasi et al., Classical Quantum Gravity **32**, 115012 (2015).

[2] T. Accadia et al., J. Instrum. **7**, P03012 (2012).

[3] H. Grote, Classical Quantum Gravity **27**, 084003 (2010).

[4] B. P. Abbott et al., Phys. Rev. Lett. **116**, 061102 (2016).

[5] B. P. Abbott et al., Phys. Rev. Lett. **116**, 241103 (2016).

[6] B. P. Abbott et al., Phys. Rev. Lett. **118**, 221101 (2017).

[7] B. P. Abbott et al., Phys. Rev. Lett. **119**, 141101 (2017).

[8] B. P. Abbott et al., Phys. Rev. Lett. **119**, 161101 (2017).

[9] B. P. Abbott et al., Astrophys. J. Lett. **851**, L35 (2017).

[10] P. C. Peters, Phys. Rev. **136**, B1224 (1964).

[11] K. A. Postnov and L. R. Yungelson, Living Rev. Relativity **17**, 3 (2014).

[12] B. J. Shappee and T. A. Thompson, Astrophys. J. **766**, 64 (2013).

[13] F. Antonini, S. Chatterjee, C. L. Rodriguez, M. Morscher, B. Pattabiraman, V. Kalogera, and F. A. Rasio, Astrophys. J. **816**, 65 (2016).

[14] F. Antonini, S. Toonen, and A. S. Hamers, Astrophys. J. **841**, 77 (2017).

[15] A. Nishizawa, E. Berti, A. Klein, and A. Sesana, Phys. Rev. D **94**, 064020 (2016).

[16] A. Nishizawa, A. Sesana, E. Berti, and A. Klein, Mon. Not. R. Astron. Soc. **465**, 4375 (2017).

[17] K. Breivik, C. L. Rodriguez, S. L. Larson, V. Kalogera, and F. A. Rasio, Astrophys. J. **830**, L18 (2016).

[18] C. Petrovich and F. Antonini, Astrophys. J. **846**, 146 (2017).

[19] J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, Astrophys. J. **853**, 140 (2018).

[20] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, and F. A. Rasio, Phys. Rev. Lett. **120**, 151101 (2018).

[21] B.-M. Hoang, S. Naoz, B. Kocsis, F. A. Rasio, and F. Dosopoulou, Astrophys. J. **856**, 140 (2018).

[22] J. Samsing, D. J. D'Orazio, A. Askar, and M. Giersz, arXiv:1802.08654.

[23] B. P. Abbott et al., Classical Quantum Gravity **34**, 104002 (2017).

[24] O. Blaes, M. H. Lee, and A. Socrates, Astrophys. J. **578**, 775 (2002).

[25] L. Hoffman and A. Loeb, Mon. Not. R. Astron. Soc. **377**, 957 (2007).

[26] P. Amaro-Seoane, A. Sesana, L. Hoffman, M. Benacquista, C. Eichhorn, J. Makino, and R. Spurzem, Mon. Not. R. Astron. Soc. **402**, 2308 (2010).

[27] M. Bonetti, E. Barausse, G. Faye, F. Haardt, and A. Sesana, Classical Quantum Gravity **34**, 215004 (2017).

[28] M. Bonetti, F. Haardt, A. Sesana, and E. Barausse, Mon. Not. R. Astron. Soc. **477**, 3910 (2018).

[29] M. Bonetti, A. Sesana, E. Barausse, and F. Haardt, Mon. Not. R. Astron. Soc. **477**, 2599 (2018).

[30] A. Klein and P. Jetzer, Phys. Rev. D **81**, 124001 (2010).

[31] R.-M. Memmesheimer, A. Gopakumar, and G. Schäfer, Phys. Rev. D **70**, 104011 (2004).

[32] T. Damour, A. Gopakumar, and B. R. Iyer, Phys. Rev. D **70**, 064028 (2004).

[33] C. Königsdörffer and A. Gopakumar, Phys. Rev. D **73**, 124012 (2006).

[34] K. G. Arun, L. Blanchet, B. R. Iyer, and M. S. S. Qusailah, Phys. Rev. D **77**, 064034 (2008).

[35] K. G. Arun, L. Blanchet, B. R. Iyer, and M. S. S. Qusailah, Phys. Rev. D **77**, 064035 (2008).

[36] K. G. Arun, L. Blanchet, B. R. Iyer, and S. Sinha, Phys. Rev. D **80**, 124018 (2009).

[37] C. K. Mishra, K. G. Arun, and B. R. Iyer, Phys. Rev. D **91**, 084040 (2015).

[38] L. Á. Gergely, Z. I. Perjés, and M. Vasúth, Phys. Rev. D **58**, 124001 (1998).

[39] L. Á. Gergely, Phys. Rev. D **61**, 024035 (1999).

[40] L. Á. Gergely and Z. Keresztes, Phys. Rev. D **67**, 024020 (2003).

[41] B. Mikóczi, M. Vasúth, and L. Á. Gergely, Phys. Rev. D **71**, 124043 (2005).

[42] Z. Keresztes, B. Mikóczi, and L. Á. Gergely, Phys. Rev. D **72**, 104022 (2005).

[43] N. Yunes, K. G. Arun, E. Berti, and C. M. Will, Phys. Rev. D **80**, 084001 (2009).

[44] N. J. Cornish and J. S. Key, Phys. Rev. D **82**, 044028 (2010).

[45] N. J. Cornish and J. S. Key, Phys. Rev. D **84**, 029901(E) (2011).

[46] J. S. Key and N. J. Cornish, Phys. Rev. D **83**, 083001 (2011).

[47] A. Gopakumar and G. Schäfer, Phys. Rev. D **84**, 124007 (2011).

[48] E. A. Huerta, P. Kumar, S. T. McWilliams, R. O'Shaughnessy, and N. Yunes, Phys. Rev. D **90**, 084016 (2014).

[49] S. Tanay, M. Haney, and A. Gopakumar, Phys. Rev. D **93**, 064031 (2016).

[50] B. Moore, M. Favata, K. G. Arun, and C. K. Mishra, Phys. Rev. D **93**, 124061 (2016).

[51] E. A. Huerta, P. Kumar, B. Agarwal, D. George, H.-Y. Schive, H. P. Pfeiffer, R. Haas, W. Ren, T. Chu, M. Boyle, D. A. Hemberger, L. E. Kidder, M. A. Scheel, and B. Szilagyi, Phys. Rev. D **95**, 024038 (2017).

[52] E. A. Huerta, C. J. Moore, P. Kumar, D. George, A. J. K. Chua, R. Haas, E. Wessel, D. Johnson, D. Glennon, A. Rebei, A. M. Holgado, J. R. Gair, and H. P. Pfeiffer, Phys. Rev. D **97**, 024031 (2018).

[53] I. Hinder, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. D **98**, 044015 (2018).

[54] T. Hinderer and S. Babak, Phys. Rev. D **96**, 104048 (2017).

[55] Z. Cao and W.-B. Han, Phys. Rev. D **96**, 044028 (2017).

[56] A. Klein, N. Cornish, and N. Yunes, Phys. Rev. D **90**, 124029 (2014).

[57] B. M. Barker and R. F. O'Connell, Phys. Rev. D **12**, 329 (1975).

[58] C. Cutler, D. Kennefick, and E. Poisson, Phys. Rev. D **50**, 3816 (1994).

[59] N. Loutrel, S. Liebersbach, N. Yunes, and N. Cornish, arXiv:1801.09009.

[60] E. Racine, Phys. Rev. D **78**, 044021 (2008).

[61] C. Cutler, Phys. Rev. D **57**, 7089 (1998).

[62] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D **49**, 6274 (1994).
[63] Y. Boetzel, A. Susobhanan, A. Gopakumar, A. Klein, and P. Jetzer, Phys. Rev. D **96**, 044011 (2017).

[64] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevD.98.104043 for a MATHEMATICA version of all amplitudes to order $O(y^4, e^{10})$ or upon request from boetzel@physik.uzh.ch.

[65] C. M. Bender and S. A. Orszag, *Advanced Mathematical Methods for Scientists and Engineers 1, Asymptotic Methods and Perturbation Theory* (Springer, New York, 1999).

[66] K. Chatziioannou, A. Klein, N. Yunes, and N. J. Cornish, Phys. Rev. D **95**, 104004 (2017).

[67] E. Poisson, Phys. Rev. D **57**, 5287 (1998).