# EFPE: Efficient fully precessing eccentric gravitational waveforms for
binaries with long inspirals

Antoine Klein  School of Physics and Astronomy & Institute for Gravitational
Wave Astronomy, University of Birmingham, Birmingham, B15 2TT, UK

###### Abstract

In this paper, we derive a set of equations of motions for binaries on
eccentric orbits undergoing spin-induced precession that can efficiently be
integrated on the radiation-reaction timescale. We find a family of solutions
with a computation cost improved by a factor $10$ \- $50$ down to $\sim 10$ ms
per waveform evaluation compared to waveforms obtained by directly integrating
the precession equations, that maintain a mismatch of the order $10^{-4}$ \-
$10^{-6}$ for waveforms lasting a million orbital cycles and a thousand spin-
induced precession cycles. We express it in terms of parameters that make the
solution regular in the equal-mass limit, thus bypassing a problem of previous
similar solutions. We point to ways in which the solution presented in this
paper can be perturbed to take into account effects such as general quadrupole
momenta and post-Newtonian corrections to the precession equations. This new
waveform, with its improved efficiency and its accuracy, makes possible
Bayesian parameter estimation using the full spin and eccentricity parameter
volume for long lasting inspiralling signals such as stellar-origin black hole
binaries observed by LISA.

###### pacs:

04.30.-w, 04.30.Tv

##  I Introduction

Gravitational wave (GW) observations have opened a new field of astronomy. The
ground-based network [1, 2, 3, 4, 5] has been able to constrain the population
of stellar-mass black holes [6, 7, 8] and detect neutron-star binaries [9],
allowing to place important theoretical constraints [10]. Pulsar timing array
observations [11] have recently been announced to have detected a very low
frequency stochastic GW background [12], potentially having implications on
its astrophysical source [13, 14]. The LISA detector [15] will, among other
observations, allow to vastly broaden the binary parameter space available for
detections, including massive black hole binary mergers [16, 17], extreme mass
ratio inspirals of stellar-mass compact objects around massive black holes
[18], galactic white dwarf and neutron star binaries [19, 20], and early-
inspiral stellar-mass black hole binaries [21].

The physics of inspiral binaries is rich, including effects such as spin-
induced precession [22, 23], tidal interactions [24], and orbital eccentricity
[25, 26]. The simultaneous measurement of those effects is an important
ingredient in discriminating between different black hole binary formation
scenarios [27, 28], and massive black hole binary evolution [29]. Efforts
towards building waveforms for preccessing and/or eccentric waveforms have
been carried out using various formalisms [23, 30, 31, 32, 33, 34, 35, 36].

Building waveforms for precessing binaries is complicated by the fact that
they experience variations on three different timescales: the longer
radiation-reaction timescale, the intermediate precession timescale, and the
faster orbital timescale [23]. Including eccentricity adds a fourth timescale,
the periastron precession timescale, of the same order as the precession
timescale [26]. It is possible to express the solutions in such a way as to
have to integrate the equations of motion only on the precession timescale,
inducing a significant speed-up in the computation of those waveforms [37,
38]. However, for systems with long inspirals such as neutron star or low-mass
black hole binaries in third-generation ground-based detectors, or stellar-
origin black hole binaries with LISA, using such a solution is still
computationally prohibitive. In this work, we use a previous expression of the
circular equations of motion that can be integrated analytically in the
absence of radiation reaction [39, 40], and slightly modify it in order to
regularize the equal-mass limit and include orbital eccentricity. We find a
complete analytical solution to these equations, that we can perturb in order
to take into account radiation-reaction effects. The additional equations can
be integrated on the radiation-reaction timescale, leading to a waveform with
a computational cost similar to ones describing nonprecessing binaries. Such a
waveform was recently used to demonstrate its applicability for Bayesian
parameter estimation in the full 17 dimensional parameter space of precessing
eccentric binaries [41].

While the precession equations had been solved in simplified situations such
as equal-mass binaries or binaries with a single nonvanishing spin [23, 42,
43, 30], and generic spin configuration waveforms had been mapped to simpler
ones [31], a general solution to these equations had eluded the modelling
community for a long time. Kesden  _et al._ [39, 40] were able to express the
problem in such a way as to retain only one dynamical parameter by expressing
the equations of precession in a non-inertial frame, thus finding a solution
in the absence of radiation reaction requiring the integration of a single
equation. This was later extended to account for radiation reaction effects,
including an analytical solution to this equation [44, 45]. One problem
remained in this formulation: the dynamical variable chosen to express the
solutions was the norm of the total spin $S$, which becomes constant in the
equal-mass limit, which introduces a singularity in the solution. While this
problem can safely be ignored in a large portion of the parameter space [46],
it can almost completely disappear by choosing a different parameter in which
to express the solutions, as we show in Sec. II.

Throughout this paper, we use geometric units with $G=c=1$, we write vectors
in boldface, and we write unit vectors with a hat. For convenience, we write
angular momenta as dimensionless quantities with e.g.
$\bm{S}=\bm{S}^{\text{phys}}/M^{2}$.

The paper is organized as follows. In Sec. II, we derive a general analytic
solution to the conservative equations of precession in the presence of
eccentricity. In Sec.III, we extend this solution to include radiation
reaction effects. In Sec. IV, we describe how to construct the two linear
gravitational wave polarizations in the Fourier domain using this solution. In
Sec. V, we describe the results of numerical simulations we performed to
assess the accuracy and the computational efficiency of these waveforms. In
Sec. VI, we discuss possible ways of extending these results to take into
account additional effects such as those involved in the description of
extended bodies. We conclude in Sec. VII.

##  II Solution of the conservative problem

We begin by writing the precession equations for black holes on eccentric
orbits, in the absence of radiation reaction, including leading post-
Newtionian (PN) order spin-orbit and spin-spin interactions [22, 38]:

$\displaystyle\mathcal{D}\bm{\hat{L}}$
$\displaystyle=-y^{6}\left(\bm{\Omega}_{1}+\bm{\Omega}_{2}\right),$ (1)
$\displaystyle\mathcal{D}\bm{s}_{1}$
$\displaystyle=\mu_{2}y^{5}\bm{\Omega}_{1},$ (2)
$\displaystyle\mathcal{D}\bm{s}_{2}$
$\displaystyle=\mu_{1}y^{5}\bm{\Omega}_{2},$ (3)

where we defined

$\displaystyle\mathcal{D}$
$\displaystyle=\frac{M}{\left(1-e^{2}\right)^{3/2}}\frac{d}{dt},$ (4)
$\displaystyle y$ $\displaystyle=\frac{(M\omega)^{1/3}}{\sqrt{1-e^{2}}},$ (5)
$\displaystyle\mu_{i}$ $\displaystyle=\frac{m_{i}}{M},$ (6)
$\displaystyle\bm{s}_{i}$ $\displaystyle=\frac{\bm{S}_{i}}{\mu_{i}},$ (7)
$\displaystyle L$ $\displaystyle=\frac{\nu}{y},$ (8)
$\displaystyle\bm{\Omega}_{i}$
$\displaystyle=\left[\frac{1}{2}\mu_{i}+\frac{3}{2}\left(1-y\bm{\hat{L}}\cdot\bm{s}\right)\right]\bm{\hat{L}}\times\bm{s}_{i}+\frac{1}{2}y\bm{s}_{j}\times\bm{s}_{i},$
(9)

$m_{i}$ are the individual masses, $\bm{S}_{i}$ are the individual spins,
$\bm{L}$ is the Newtonian angular momentum, $\mathcal{D}$ is a differential
operator related to the derivative with respect to the mean orbital phase,
$M=m_{1}+m_{2}$ is the total mass, $e$ is the orbital eccentricity, $y$ is a
PN parameter related to the norm of $\bm{L}$, $\omega$ is the mean orbital
frequency, $\mu_{i}$ are dimensionless mass parameters, $\nu=\mu_{1}\mu_{2}$
is the symmetric mass ratio, $\bm{s}_{i}$ are the reduced individual spins,
and $\bm{s}=\bm{s}_{1}+\bm{s}_{2}$ is the total reduced spin.

The precession equations for $\bm{L},\bm{s}_{1},\bm{s}_{2}$ contain a number
of conserved quantities: the norms of the vectors $L$, $s_{1}$, and $s_{2}$
are conserved, as well as the total angular momentum vector
$\bm{J}=\bm{L}+\mu_{1}\bm{s}_{1}+\mu_{2}\bm{s}_{2}$, leading to six conserved
quantities. The effective spin [47]

$\displaystyle\chi_{\text{eff}}$ $\displaystyle=\bm{\hat{L}}\cdot\bm{s}$ (10)

is a seventh conserved quantity, leaving only two dynamical variables needed
to completely characterize the three angular momenta.

###  II.1 Solution in a rotating frame

We can express the angular momenta in a frame where the $z$-axis is aligned
with $\bm{J}$, and where the orbital angular momentum $\bm{L}$ is
perpendicular to the $y$-axis, with $\bm{\hat{L}}\cdot\bm{\hat{x}}\geq 0$ (see
Fig 1). In this frame, the angular momenta can all be expressed in terms of
one variable. As we mentioned earlier, the choice of [39, 40] to use the norm
of the total spin $S$ is possible, but leads to a singularity in the equal-
mass limit. Instead we choose to use

$\displaystyle\delta\chi$
$\displaystyle=\bm{\hat{L}}\cdot\left(\bm{s}_{1}-\bm{s}_{2}\right).$ (11)

We can express the orbital angular momentum and the total spin in the
precessing frame as

$\displaystyle\bm{\hat{L}}$
$\displaystyle=\cos\theta_{L}\bm{\hat{z}}+\sin\theta_{L}\bm{\hat{x}},$ (12)
$\displaystyle\bm{S}$
$\displaystyle=S\cos\theta_{S}\bm{\hat{z}}-S\sin\theta_{S}\bm{\hat{x}},$ (13)
$\displaystyle S^{2}$
$\displaystyle=J^{2}-L^{2}-L(\chi_{\text{eff}}+\delta\mu\delta\chi),$ (14)
$\displaystyle\cos\theta_{L}$
$\displaystyle=\frac{1}{2J}\left(2L+\chi_{\text{eff}}+\delta\mu\delta\chi\right),$
(15) $\displaystyle\cos(\theta_{L}+\theta_{S})$
$\displaystyle=\frac{1}{2S}\left(\chi_{\text{eff}}+\delta\mu\delta\chi\right),$
(16)

where $J$ is the norm of the total angular momentum,
$\bm{S}=\bm{S}_{1}+\bm{S}_{2}$ is the total spin, and
$\delta\mu=\mu_{1}-\mu_{2}$ is the dimensionless mass difference. We can see
that $S^{2}$ and $\cos\theta_{L}$ are related to $\delta\chi$ by an affine
transformation, but they become conserved in the equal-mass limit where
$\delta\mu\to 0$, unlike $\delta\chi$. Therefore, a solution written in terms
of any one of those three variables is equivalent as long as $\delta\mu\neq
0$. The fact that $S$ and $\cos\theta_{L}$ are conserved in the equal-mass
limit is the cause of the singularity arising when using one of those
parameters to describe precession, as the individual spins $\bm{S}_{i}$
themselves are not conserved.

$\bm{\hat{z}}=\bm{\hat{J}}$
---
$\bm{\hat{x}}$
---
$\bm{\hat{y}}$
---
$\bm{S}_{1}$
---
$\bm{S}_{2}$
---
$\bm{L}$
---
$m_{1}$
---
$m_{2}$
---
$\theta_{L}$
---
Figure 1: Graphical representation of the rotating frame used to describe the
angular momenta. The $z$ axis is fixed onto the total angular momentum
direction $\bm{\hat{J}}$, and the $x$ axis is rotating together with the
orbital angular momentum $\bm{L}$ so that the latter stays in the $x$-$z$
plane, with $\bm{\hat{x}}\cdot\bm{\hat{L}}\geq 0$.

In order to express the individual spins in terms of these variables, we can
construct a frame described by $\bm{\hat{z}}^{\prime}=\bm{\hat{S}}$,
$\bm{\hat{y}}^{\prime}=\bm{\hat{y}}$, and
$\bm{\hat{x}}^{\prime}=\bm{\hat{y}}^{\prime}\times\bm{\hat{z}}^{\prime}$. If
we define the spherical angles of $\bm{s}_{1}$ in this frame to be
$(\theta^{\prime},\phi^{\prime})$, we can write

$\displaystyle\cos\theta^{\prime}$
$\displaystyle=\frac{1}{2SS_{1}}\left[S^{2}+S_{1}^{2}-S_{2}^{2}\right],$ (17)
$\displaystyle\cos\phi^{\prime}$
$\displaystyle=\frac{1}{4A_{L}A_{P}}\big{[}S^{2}(\delta\mu\chi_{\text{eff}}+\delta\chi)$
$\displaystyle-\left(S_{1}^{2}-S_{2}^{2}\right)(\chi_{\text{eff}}+\delta\mu\delta\chi)\big{]},$
(18) $\displaystyle A_{L}$
$\displaystyle=J\sin\theta_{L}=S\sin(\theta_{L}+\theta_{S}),$ (19)
$\displaystyle A_{P}$ $\displaystyle=SS_{1}\sin\theta^{\prime},$ (20)
$\displaystyle\text{sign}\left(\sin\phi^{\prime}\right)$
$\displaystyle=\text{sign}\left[\left(\bm{\hat{L}}\times\bm{s}_{1}\right)\cdot\bm{s}_{2}\right].$
(21)

We can then reconstruct the second spin with $\bm{S}_{2}=\bm{S}-\bm{S}_{1}$.

Using these relations together with the precession equations, we can find a
differential equation for $\delta\chi$

$\displaystyle\left(\mathcal{D}\delta\chi\right)^{2}$
$\displaystyle=\frac{9}{4}A^{2}y^{11}\left(\delta\mu\delta\chi^{3}+B\delta\chi^{2}+C\delta\chi+D\right),$
(22) $\displaystyle=-\frac{1}{y}\left(\frac{3}{2}Ay^{6}\right)^{2}$
$\displaystyle\times(\delta\chi-\delta\chi_{+})(\delta\chi-\delta\chi_{-})(\delta\chi_{3}-\delta\mu\delta\chi),$
(23) $\displaystyle A$ $\displaystyle=1-y\chi_{\text{eff}},$ (24)
$\displaystyle\delta\chi_{-}$
$\displaystyle\leq\delta\chi_{+}\leq\frac{\delta\chi_{3}}{\delta\mu},$ (25)

where the constant coefficients $B$, $C$, and $D$ are given in appendix B. A
similar equation was found in [44, 45] in terms of $S^{2}$, and we can find a
similar analytic solution in terms of the roots of the cubic polynomial:

$\displaystyle\delta\chi$
$\displaystyle=\delta\chi_{-}+(\delta\chi_{+}-\delta\chi_{-})\text{sn}^{2}(\psi_{p},m),$
(26) $\displaystyle m$
$\displaystyle=\frac{\delta\mu(\delta\chi_{+}-\delta\chi_{-})}{\delta\chi_{3}-\delta\mu\delta\chi_{-}},$
(27) $\displaystyle\mathcal{D}\psi_{\text{p}}$
$\displaystyle=\frac{3Ay^{6}}{4}\sqrt{\frac{1}{y}(\delta\chi_{3}-\delta\mu\delta\chi_{-})},$
(28)

where $\text{sn}(\psi_{p},m)$ is a Jacobi elliptic function with parameter
$m$.

In order to find the roots of the polynomial in terms of its coefficients, we
define

$\displaystyle p$
$\displaystyle=\frac{1}{y^{2}}\left(\frac{B^{2}}{3}-\delta\mu C\right),$ (29)
$\displaystyle q$
$\displaystyle=\frac{1}{y^{3}}\left(\frac{2B^{3}}{27}-\delta\mu\frac{BC}{3}+\delta\mu^{2}D\right).$
(30)

The solutions can then be expressed in terms of

$\displaystyle Y_{3}$
$\displaystyle=2\left|\frac{p}{3}\right|^{1/2}\cos\left[\frac{\arg\left(G\right)}{3}\right],$
(31) $\displaystyle Y_{\pm}$
$\displaystyle=2\left|\frac{p}{3}\right|^{1/2}\cos\left[\frac{\arg\left(G\right)\mp
2\pi}{3}\right],$ (32) $\displaystyle G$
$\displaystyle=-\frac{q}{2}+i\left[\left(\frac{p}{3}\right)^{3}-\left(\frac{q}{2}\right)^{2}\right]^{1/2},$
(33) $\displaystyle dY$ $\displaystyle=\frac{B}{3y},$ (34)

as

$\displaystyle\delta\chi_{3}$ $\displaystyle=y\left(Y_{3}-dY\right),$ (35)
$\displaystyle\delta\chi_{\pm}$
$\displaystyle=\frac{y}{\delta\mu}\left(Y_{\pm}-dY\right).$ (36)

Using those quantities, we can write

$\displaystyle\mathcal{D}\psi_{\text{p}}$
$\displaystyle=\frac{3Ay^{6}}{4}\sqrt{Y_{3}-Y_{-}},$ (37) $\displaystyle m$
$\displaystyle=\frac{Y_{+}-Y_{-}}{Y_{3}-Y_{-}}.$ (38)

Since $Y_{j}$ are regular in the equal-mass limit and $Y_{+}\leq Y_{3}$, the
solution described by the two equations above is regular in this limit as
well. A more detailed discussion of this solution is given in appendix B,
including the limits $\delta\mu\ll 1$ and $y\ll 1$.

###  II.2 Solutions in an inertial frame

In order to reconstruct the angular momenta in an inertial frame, we need to
compute the angular velocity of the rotation of the precessing frame with
respect to an inertial one. Since the precessing frame is defined so that
$\bm{\hat{L}}$ stays in the $x$-$z$ plane, and is related to an inertial frame
by a rotation around $\bm{\hat{J}}$ by some angle $\phi_{z}$, we can write

$\displaystyle\mathcal{D}\phi_{z}$
$\displaystyle=\frac{1}{\sin^{2}\theta_{L}}\left(\mathcal{D}\bm{\hat{L}}\right)\cdot\left(\bm{\hat{J}}\times\bm{\hat{L}}\right).$
(39)

Using the expression of the angular momenta and the precession equations, we
can write this in terms of the solution we found for $\delta\chi$ as

$\displaystyle\mathcal{D}\phi_{z}$
$\displaystyle=\frac{Jy^{6}}{2}+\frac{\mathcal{D}\psi_{p}}{\nu\sqrt{Y_{3}-Y_{-}}}\left[\frac{N_{+}}{D_{+}}+\frac{N_{-}}{D_{-}}\right],$
(40) $\displaystyle N_{+}$
$\displaystyle=(J+L)(J+L+2\nu\chi_{\text{eff}})-\delta\mu\left(S_{1}^{2}-S_{2}^{2}\right),$
(41) $\displaystyle N_{-}$
$\displaystyle=(J-L)(J-L-2\nu\chi_{\text{eff}})-\delta\mu\left(S_{1}^{2}-S_{2}^{2}\right),$
(42) $\displaystyle D_{+}$
$\displaystyle=2(J+L)+\chi_{\text{eff}}+\delta\mu\delta\chi,$ (43)
$\displaystyle=B_{+}-C_{+}\left[1-2\text{sn}^{2}(\psi_{p},m)\right],$ (44)
$\displaystyle D_{-}$
$\displaystyle=2(J-L)-\chi_{\text{eff}}-\delta\mu\delta\chi,$ (45)
$\displaystyle=B_{-}-C_{-}\left[1-2\text{sn}^{2}(\psi_{p},m)\right].$ (46)

We can integrate this equation analytically, and separate the secular part
from the periodic part:

$\displaystyle\phi_{z}$ $\displaystyle=\phi_{z,0}+\delta\phi_{z},$ (47)
$\displaystyle\mathcal{D}\phi_{z,0}$
$\displaystyle=\frac{Jy^{6}}{2}+\frac{3(1-y\chi_{\text{eff}})y^{6}}{4\nu
K(m)}\left[\frac{N_{+}}{B_{+}-C_{+}}\Pi\left(\frac{-2C_{+}}{B_{+}-C_{+}},m\right)+\frac{N_{-}}{B_{-}-C_{-}}\Pi\left(\frac{-2C_{-}}{B_{-}-C_{-}},m\right)\right],$
(48) $\displaystyle\delta\phi_{z}$
$\displaystyle=\frac{N_{+}}{\nu(B_{+}-C_{+})(Y_{3}-Y_{-})^{1/2}}\left\\{\Pi\left[\frac{-2C_{+}}{B_{+}-C_{+}};\text{am}(\hat{\psi}_{p},m);m\right]-\frac{\Pi\left(\frac{-2C_{+}}{B_{+}-C_{+}},m\right)}{K(m)}\hat{\psi}_{p}\right\\}$
$\displaystyle+\frac{N_{-}}{\nu(B_{-}-C_{-})(Y_{3}-Y_{-})^{1/2}}\left\\{\Pi\left[\frac{-2C_{-}}{B_{-}-C_{-}};\text{am}(\hat{\psi}_{p},m);m\right]-\frac{\Pi\left(\frac{-2C_{-}}{B_{-}-C_{-}},m\right)}{K(m)}\hat{\psi}_{p}\right\\},$
(49)

where $\Pi(n,m)$ is the complete elliptic integral of the third kind,
$\Pi(n;\phi;m)$ is the incomplete elliptic integral of the third kind, and
$K(m)$ is the complete elliptic integral of the first kind. The angle
$\hat{\psi}_{p}$ satisfies $-K(m)<\hat{\psi}_{p}\leq K(m)$, and
$\psi_{p}-\hat{\psi}_{p}=2nK(m)$, for some $n\in\mathbb{Z}$. The conventions
for the elliptic integrals used in this paper are given in appendix A.

In order to build a gravitational waveform, it is useful to compute the
solution for a related angle $\zeta$ satisfying
$\mathcal{D}\zeta=-\cos\theta_{L}\mathcal{D}\phi_{z}$ [48, 49]. With the
expressions computed earlier, we find the following solution, similarly
separated into a secular and a periodic part:

$\displaystyle\zeta$ $\displaystyle=\zeta_{0}+\delta\zeta,$ (50)
$\displaystyle\mathcal{D}\zeta_{0}$
$\displaystyle=-\frac{\left(2L+\chi_{\text{eff}}+\delta\mu\delta\chi_{\text{av}}\right)y^{6}}{4}+\frac{\delta\mu\delta\chi_{\text{diff}}\
y^{6}}{2m}\left[\frac{E(m)}{K(m)}-1+\frac{m}{2}\right]-\frac{3(L+\nu\chi_{\text{eff}})(1-y\chi_{\text{eff}})y^{6}}{2\nu}$
$\displaystyle+\frac{3(1-y\chi_{\text{eff}})y^{6}}{4\nu
K(m)}\left[\frac{N_{+}}{B_{+}-C_{+}}\Pi\left(\frac{-2C_{+}}{B_{+}-C_{+}},m\right)-\frac{N_{-}}{B_{-}-C_{-}}\Pi\left(\frac{-2C_{-}}{B_{-}-C_{-}},m\right)\right],$
(51) $\displaystyle\delta\zeta$
$\displaystyle=\frac{2\delta\mu\delta\chi_{\text{diff}}}{3m(1-y\chi_{\text{eff}})(Y_{3}-Y_{-})^{1/2}}\left\\{E\left[\text{am}(\hat{\psi}_{p},m);m\right]-\frac{E(m)}{K(m)}\hat{\psi}_{p}\right\\}$
$\displaystyle+\frac{N_{+}}{\nu(B_{+}-C_{+})(Y_{3}-Y_{-})^{1/2}}\left\\{\Pi\left[\frac{-2C_{+}}{B_{+}-C_{+}};\text{am}(\hat{\psi}_{p},m);m\right]-\frac{\Pi\left(\frac{-2C_{+}}{B_{+}-C_{+}},m\right)}{K(m)}\hat{\psi}_{p}\right\\}$
$\displaystyle-\frac{N_{-}}{\nu(B_{-}-C_{-})(Y_{3}-Y_{-})^{1/2}}\left\\{\Pi\left[\frac{-2C_{-}}{B_{-}-C_{-}};\text{am}(\hat{\psi}_{p},m);m\right]-\frac{\Pi\left(\frac{-2C_{-}}{B_{-}-C_{-}},m\right)}{K(m)}\hat{\psi}_{p}\right\\},$
(52)

where

$\displaystyle\lim_{m\to
0}\frac{1}{m}\left[\frac{E(m)}{K(m)}-1+\frac{m}{2}\right]=0,$ (53)

$E(m)$ is the complete elliptic integral of the second kind, $E(\phi;m)$ is
the incomplete elliptic integral of the second kind, $\text{am}(\psi,m)$ is
the Jacobi amplitude satisfying $\sin[\text{am}(\psi,m)]=\text{sn}(\psi,m)$,
$\delta\chi_{\text{av}}=(\delta\chi_{+}+\delta\chi_{-})/2$, and
$\delta\chi_{\text{diff}}=(\delta\chi_{+}-\delta\chi_{-})/2$.

###  II.3 Approximations of the solution

The full solution presented above can be sped up by making some
approximations. Since the solutions for $\phi_{z}$ and $\zeta$ are well
separated into a secular and a periodic part, one obvious simplification would
be to neglect the periodic part. This does save evaluation time as we will
discuss in Sec. V, but more time can be saved at a moderate cost in
faithfulness. As we show in appendix B, a PN expansion of the parameter $m$,
valid if $\delta\mu\neq 0$, yields $m=\mathcal{O}\left(y^{2}\right)$.
Furthermore, if $\delta\mu$ is small, we can write
$m=\mathcal{O}\left(\delta\mu\right)$. We therefore expect $m$ to be small in
a large number of cases, justifying the approximation $m\approx 0$. We show
distributions for $m$ taken from randomized angular momenta configurations at
different times in the inspiral in Fig. 2. In the limit $m\to 0$, we find the
simpler expressions

$\displaystyle\mathcal{D}\phi_{z,0}$
$\displaystyle=\frac{Jy^{6}}{2}+\frac{3(1-y\chi_{\text{eff}})y^{6}}{4\nu}\left[\frac{N_{+}}{\left(B_{+}^{2}-C_{+}^{2}\right)^{1/2}}+\frac{N_{-}}{\left(B_{-}^{2}-C_{-}^{2}\right)^{1/2}}\right],$
(54) $\displaystyle\mathcal{D}\zeta_{0}$
$\displaystyle=\frac{\left(2L+\chi_{\text{eff}}+\delta\mu\delta\chi_{\text{av}}\right)y^{6}}{4}+\frac{3(L+\nu\chi_{\text{eff}})(1-y\chi_{\text{eff}})y^{6}}{2\nu}-\frac{3(1-y\chi_{\text{eff}})y^{6}}{4\nu}\left[\frac{N_{+}}{\left(B_{+}^{2}-C_{+}^{2}\right)^{1/2}}-\frac{N_{-}}{\left(B_{-}^{2}-C_{-}^{2}\right)^{1/2}}\right],$
(55) $\displaystyle\delta\phi_{z}$
$\displaystyle=\frac{N_{+}}{\nu\left[(Y_{3}-Y_{-})\left(B_{+}^{2}-C_{+}^{2}\right)\right]^{1/2}}\left\\{\arctan\left[\left(\frac{B_{+}+C_{+}}{B_{+}-C_{+}}\right)^{1/2}\tan\hat{\psi}_{p}\right]-\hat{\psi}_{p}\right\\}$
$\displaystyle+\frac{N_{-}}{\nu\left[(Y_{3}-Y_{-})\left(B_{-}^{2}-C_{-}^{2}\right)\right]^{1/2}}\left\\{\arctan\left[\left(\frac{B_{-}+C_{-}}{B_{-}-C_{-}}\right)^{1/2}\tan\hat{\psi}_{p}\right]-\hat{\psi}_{p}\right\\},$
(56) $\displaystyle\delta\zeta$
$\displaystyle=-\frac{\delta\mu\delta\chi_{\text{diff}}\sin
2\hat{\psi}_{p}}{6(1-y\chi_{\text{eff}})(Y_{3}-Y_{-})^{1/2}}-\frac{N_{+}}{\nu\left[(Y_{3}-Y_{-})\left(B_{+}^{2}-C_{+}^{2}\right)\right]^{1/2}}\left\\{\arctan\left[\left(\frac{B_{+}+C_{+}}{B_{+}-C_{+}}\right)^{1/2}\tan\hat{\psi}_{p}\right]-\hat{\psi}_{p}\right\\}$
$\displaystyle+\frac{N_{-}}{\nu\left[(Y_{3}-Y_{-})\left(B_{-}^{2}-C_{-}^{2}\right)\right]^{1/2}}\left\\{\arctan\left[\left(\frac{B_{-}+C_{-}}{B_{-}-C_{-}}\right)^{1/2}\tan\hat{\psi}_{p}\right]-\hat{\psi}_{p}\right\\},$
(57)

where $\hat{\psi}_{p}$ is computed from $\psi_{p}$ using $K(m=0)=\pi/2$.

Figure 2: Distributions of the Jacobi elliptic parameter $m$ at different
points in the inspiral. The solid blue line is computed at $y\approx 0.027$,
corresponding to a $50M_{\odot}$-$50M_{\odot}$ system four years before
merger, the dashed orange line is computed at $y\approx 0.12$, corresponding
to the same system when it exits the LISA band at $2f_{\text{orb}}=1$ Hz, and
the dotted green line is computed at $y=6^{-1/2}$, when the binary enters the
merger stage. While $m$ does reach values close to its maximum
$m_{\text{max}}=1$, it is significantly smaller than unity in most cases in
the inspiral.

To produce Fig. 2, we computed $10^{6}$ different realizations of angular
momenta configurations, with all unit vectors uniformly distributed on the
sphere, both dimensionless spin magnitudes uniformly distributed between $0$
and $1$, and both masses uniformly distributed between $10M_{\odot}$ and
$100M_{\odot}$. We choose three different stages of the inspiral, to provide a
sense of the evolution of the Jacobi elliptic parameter $m$: $y_{0}\approx
0.027$, computed for a $50M_{\odot}$-$50M_{\odot}$ binary system on a circular
orbit four years before merger, at leading PN order; $y_{1}\approx 0.12$,
computed for the same binary system when its dominant harmonic exits the LISA
band, i.e. when the orbital frequency satisfies $2f_{\text{orb}}=1$ Hz; and
$y_{2}=6^{1/2}$, chosen as the boundary between the inspiral and the merger
stages. Looking at these distributions, we can see that $m$ tends to increase
as the binary gets closer to merger, as expected from the PN order estimation
provided in appendix B. Furthermore, these distributions justify the
approximation $m\approx 0$, especially in the early inspiral as in the case of
stellar-origin black hole binaries observed with LISA. This parameter however
does reach values close to its maximum of $m_{\text{max}}=1$, especially later
on in the inspiral.

##  III Addition of radiation reaction

In order to build a fast and accurate waveform for a precessing binary on an
elliptic orbit, one can take advantage of this analytical conservative
solution to construct a set of evolution equations that can be integrated on
the radiation reaction timescale. The solution presented in the previous
section is given in terms of several quantities that become dynamical when the
loss of energy of the system due to gravitational wave emission is taken into
account. The norms of the orbital angular momentum and of the total angular
momentum $L$ and $J$, the square eccentricity $e^{2}$, as well as the
direction of the total angular momentum $\bm{\hat{J}}$ become dynamical. The
norms of the individual spins $S_{1}$ and $S_{2}$ and the effective spin
$\chi_{\text{eff}}$ however are still conserved.

The evolution of $L$ and $e^{2}$ can be expressed as post-Newtonian series
[38], to which we can add evolution equations for the mean orbital phase
$\lambda$ and the argument of periastron $\delta\lambda$:

$\displaystyle\mathcal{D}y$ $\displaystyle=y^{9}\sum_{n\geq
0}a_{n}\left(y,e^{2},\bm{\hat{L}},\bm{s}_{1},\bm{s}_{2}\right)y^{n},$ (58)
$\displaystyle\mathcal{D}e^{2}$ $\displaystyle=y^{8}\sum_{n\geq
0}b_{n}\left(y,e^{2},\bm{\hat{L}},\bm{s}_{1},\bm{s}_{2}\right)y^{n},$ (59)
$\displaystyle\mathcal{D}\lambda$ $\displaystyle=y^{3},$ (60)
$\displaystyle\mathcal{D}\delta\lambda$
$\displaystyle=\frac{ky^{3}}{1+k},\quad k=y^{2}\sum_{n\geq
0}k_{n}\left(y,e^{2},\bm{\hat{L}},\bm{s}_{1},\bm{s}_{2}\right)y^{n},$ (61)

where $k$ is the periastron advance.

Those series are given in terms of $y$ and $e^{2}$ that vary on the radiation
reaction timescale, but include spin couplings that depend on the precession
timescale. In order to eliminate this dependency from the equations, a natural
choice is to substitute them with their time average, computed on the
precession timescale.

###  III.1 Spin couplings averaging

The 1.5PN spin-orbit couplings can be written in terms of

$\displaystyle\beta_{i}$ $\displaystyle=\bm{\hat{L}}\cdot\bm{s}_{i},$ (62)

with $i\in\\{1,2\\}$. We can compute their averages using our solution and
find

$\displaystyle\left\langle\beta_{1}\right\rangle$
$\displaystyle=\frac{1}{2}\left(\chi_{\text{eff}}+\left\langle\delta\chi\right\rangle\right),$
(63) $\displaystyle\left\langle\beta_{2}\right\rangle$
$\displaystyle=\frac{1}{2}\left(\chi_{\text{eff}}-\left\langle\delta\chi\right\rangle\right),$
(64) $\displaystyle\left\langle\delta\chi\right\rangle$
$\displaystyle=\delta\chi_{\text{av}}-\frac{2\delta\chi_{\text{diff}}}{m}\left[\frac{E(m)}{K(m)}-1+\frac{m}{2}\right].$
(65)

The 2PN spin-spin couplings can be expressed as linear combinations of three
different types of couplings:

$\displaystyle\sigma_{i}^{(1)}$ $\displaystyle=\bm{s}_{i}^{2},$ (66)
$\displaystyle\sigma_{i}^{(2)}$
$\displaystyle=\left(\bm{\hat{L}}\cdot\bm{s}_{i}\right)^{2},$ (67)
$\displaystyle\sigma_{i}^{(3)}$
$\displaystyle=\left|\bm{\hat{L}}\times\bm{s}_{i}\right|^{2}\cos 2\psi_{i},$
(68)

with $i\in\\{0,1,2\\}$, $\bm{s}_{0}=\bm{s}_{1}+\bm{s}_{2}$, and $\psi_{i}$
denotes the angle subtended by the periastron line and the projection of
$\bm{s}_{i}$ onto the orbital plane.

$\sigma_{1,2}^{(1)}=s_{1,2}^{2}$ and $\sigma_{0}^{(2)}=\chi_{\text{eff}}^{2}$
are conserved, and we can find the precession averages for the remaining
$\sigma_{i}^{(1,2)}$:

$\displaystyle\left\langle\sigma_{1}^{(2)}\right\rangle$
$\displaystyle=\frac{1}{4}\left(\chi_{\text{eff}}+\left\langle\delta\chi\right\rangle\right)^{2}+\frac{\delta\chi_{\text{diff}}^{2}}{8}+\frac{\delta\chi_{\text{diff}}^{2}}{m^{2}}\bigg{[}\frac{E(m)^{2}}{K(m)^{2}}$
$\displaystyle+\frac{2E(m)}{3K(m)}(m-2)+\frac{m^{2}}{8}-\frac{m}{3}+\frac{1}{3}\bigg{]},$
(69) $\displaystyle\left\langle\sigma_{2}^{(2)}\right\rangle$
$\displaystyle=\frac{1}{4}\left(\chi_{\text{eff}}-\left\langle\delta\chi\right\rangle\right)^{2}+\frac{\delta\chi_{\text{diff}}^{2}}{8}+\frac{\delta\chi_{\text{diff}}^{2}}{m^{2}}\bigg{[}\frac{E(m)^{2}}{K(m)^{2}}$
$\displaystyle+\frac{2E(m)}{3K(m)}(m-2)+\frac{m^{2}}{8}-\frac{m}{3}+\frac{1}{3}\bigg{]},$
(70) $\displaystyle\left\langle\sigma_{0}^{(1)}\right\rangle$
$\displaystyle=\frac{1}{\nu}\big{[}J^{2}-L^{2}$
$\displaystyle-\delta\mu\left(\mu_{1}s_{1}^{2}-\mu_{2}s_{2}^{2}\right)-L(\chi_{\text{eff}}+\delta\mu\left\langle\delta\chi\right\rangle)\big{]},$
(71)

where

$\displaystyle\lim_{m\to 0}\frac{1}{m^{2}}$
$\displaystyle\bigg{[}\frac{E(m)^{2}}{K(m)^{2}}+\frac{2E(m)}{3K(m)}(m-2)$
$\displaystyle+\frac{m^{2}}{8}-\frac{m}{3}+\frac{1}{3}\bigg{]}=0.$ (72)

In order to estimate the effect of $\sigma_{i}^{(3)}$, it is useful to express
the relevant quantities involved in a frame tied to the orbital plane. We can
find such a frame by rotating the precessing frame tied to $\bm{\hat{J}}$ by
an angle $\theta_{L}$ about $\bm{\hat{y}}$. In this frame,
$\bm{\hat{\bar{z}}}=\bm{\hat{L}}$,
$\bm{\hat{\bar{y}}}\propto\bm{\hat{J}}\times\bm{\hat{L}}$, and
$\bm{\hat{\bar{x}}}=\bm{\hat{\bar{y}}}\times\bm{\hat{\bar{z}}}$. We can write

$\displaystyle\bm{\hat{\bar{s}}}_{i}$
$\displaystyle=\left(\sin\bar{\theta}_{i}\cos\bar{\phi}_{i},\sin\bar{\theta}_{i}\sin\bar{\phi}_{i},\cos\bar{\theta}_{i}\right),$
(73) $\displaystyle\psi_{i}$ $\displaystyle=\bar{\phi}_{i}-\delta\lambda,$
(74)

where $\delta\lambda$ is the angle subtended by the periastron line and
$\bm{\hat{x}}^{\prime}$. Relativistic effects will cause this angle to acquire
a derivative given by the periastron precession rate

$\displaystyle\mathcal{D}\delta\lambda$
$\displaystyle=3y^{5}+\mathcal{O}\left(y^{6}\right).$ (75)

Spin precession effects will cause $\bar{\phi}_{i}$ to acquire a derivative.
We can compute it using

$\displaystyle\mathcal{D}\bar{\phi}_{i}$
$\displaystyle=\frac{1}{\sin^{2}\bar{\theta}_{i}}\left(\mathcal{D}\bm{\bar{s}}_{i}\right)\cdot\left(\bm{\hat{\bar{z}}}\times\bm{s}_{i}\right),$
(76)

where

$\displaystyle\mathcal{D}\bm{\bar{s}}_{i}$
$\displaystyle=\mathcal{D}\bm{s}_{i}-\left(\bm{\hat{\bar{x}}}\cdot\bm{s}_{i}\right)\mathcal{D}\bm{\hat{\bar{x}}}-\left(\bm{\hat{\bar{y}}}\cdot\bm{s}_{i}\right)\mathcal{D}\bm{\hat{\bar{y}}}-\left(\bm{\hat{\bar{z}}}\cdot\bm{s}_{i}\right)\mathcal{D}\bm{\hat{\bar{z}}}.$
(77)

We find

$\displaystyle\mathcal{D}\bar{\phi}_{1}$
$\displaystyle=y^{5}\left(\frac{3\mu_{2}}{2}+\frac{\nu}{2}\right)+\mathcal{O}\left(y^{6}\right),$
(78) $\displaystyle\mathcal{D}\bar{\phi}_{2}$
$\displaystyle=y^{5}\left(\frac{3\mu_{1}}{2}+\frac{\nu}{2}\right)+\mathcal{O}\left(y^{6}\right).$
(79)

We find that $\bar{\phi}_{1,2}$ and $\delta\lambda$ evolve on the same
timescale, but their derivatives cannot cancel each other, which eliminates
the possibility of a resonance. Therefore, we can make the approximation that
$\sigma_{1,2}^{(3)}$ average to zero, through the averaging of $\cos
2\psi_{1,2}$.

In order to look at the average of $\sigma_{0}^{(3)}$, we first note that if
we define $\bm{\hat{A}}$ to be a unit vector in the direction of the
periastron line, such as the Laplace-Runge-Lenz vector, and
$\bm{\hat{B}}=\bm{\hat{L}}\times\bm{\hat{A}}$, we can write

$\displaystyle\sigma_{i}^{(3)}$
$\displaystyle=\left(\bm{\hat{A}}\cdot\bm{s}_{i}\right)^{2}-\left(\bm{\hat{B}}\cdot\bm{s}_{i}\right)^{2}.$
(80)

We can see from that relation that we can rewrite

$\displaystyle\sigma_{0}^{(3)}$
$\displaystyle=\sigma_{1}^{(3)}+\sigma_{2}^{(3)}+2\left|\bm{\hat{L}}\times\bm{s}_{1}\right|\left|\bm{\hat{L}}\times\bm{s}_{2}\right|\cos(\psi_{1}+\psi_{2}).$
(81)

We can derive

$\displaystyle\mathcal{D}(\psi_{1}+\psi_{2})$
$\displaystyle=-y^{5}\left(\frac{9}{2}-\frac{\nu}{2}\right)+\mathcal{O}\left(y^{6}\right),$
(82)

which justifies neglecting the precession average of $\sigma_{0}^{(3)}$ as
well.

###  III.2 Total angular momentum

Through radiation reaction, the total angular momentum will acquire a
derivative that we can express as

$\displaystyle\mathcal{D}\bm{J}$ $\displaystyle=-|\mathcal{D}L|\bm{\hat{L}}$
(83) $\displaystyle=-\frac{\nu\mathcal{D}y}{y^{2}}\bm{\hat{L}}.$ (84)

Therefore, if we decompose $\bm{\hat{L}}$ into a parallel and a perpendicular
parts as

$\displaystyle\bm{\hat{L}}$
$\displaystyle=\bm{\hat{L}}_{\parallel}+\bm{\hat{L}}_{\perp},$ (85)
$\displaystyle\bm{\hat{L}}_{\parallel}$
$\displaystyle=\cos\theta_{L}\bm{\hat{J}},$ (86)

we can separate the angular momentum loss contribution to $J$ from the one to
$\bm{\hat{J}}$. We find

$\displaystyle\mathcal{D}J$
$\displaystyle=-\frac{L\mathcal{D}y}{2Jy}\left(2L+\chi_{\text{eff}}+\delta\mu\delta\chi\right),$
(87) $\displaystyle\mathcal{D}\bm{\hat{J}}$
$\displaystyle=-\frac{L\mathcal{D}y}{y}\sin\theta_{L}\left(\cos\phi_{z}\bm{\hat{x}}+\sin\phi_{z}\bm{\hat{y}}\right),$
(88)

where $\bm{\hat{x}}$ and $\bm{\hat{y}}$ are part of an inertial triad together
with $\bm{\hat{J}}^{(0)}$, a vector aligned with $\bm{\hat{J}}$ at some
reference time. Since $\phi_{z}$ is increasing on the spin precession
timescale, much shorter than the radiation reaction timescale, we can neglect
the evolution of the total angular momentum direction [23]. We can use the
conservative solution for $\delta\chi$ in order to find a leading-order
multiple-scale analysis solution for $J$. Separating it into a secular and a
periodic parts, we find

$\displaystyle J$ $\displaystyle=J_{0}+\delta J,$ (89)
$\displaystyle\frac{\mathcal{D}J_{0}}{\mathcal{D}y}$
$\displaystyle=-\frac{L}{2J_{0}y}\left(2L+\chi_{\text{eff}}+\delta\mu\left\langle\delta\chi\right\rangle\right)+\mathcal{O}\left(y\right),$
(90) $\displaystyle\delta J$
$\displaystyle=\frac{\nu\delta\mu\delta\chi_{\text{diff}}}{m\sqrt{Y_{3}-Y_{-}}}\left(\frac{64}{15}+\frac{56}{15}e^{2}\right)y^{2}$
$\displaystyle\times\left\\{E[\text{am}(\hat{\psi}_{p},m);m]-\frac{E(m)}{K(m)}\hat{\psi}_{p}\right\\}+\mathcal{O}\left(y^{3}\right).$
(91)

Since $J_{0}=\mathcal{O}\left(y^{-1}\right)$ and $\delta
J=\mathcal{O}\left(y^{2}\right)$, we can neglect the contribution of $\delta
J$ in our solution.

Note that, provided a PN expansion for $\left\langle\delta\chi\right\rangle$,
this leaves us with an equation that we can integrate

$\displaystyle J_{0}^{2}$
$\displaystyle=\frac{\nu^{2}}{y^{2}}+\frac{\nu}{y}\left[\chi_{\text{eff}}+\delta\mu\left\langle\delta\chi\right\rangle^{(0)}\right]$
$\displaystyle+C_{0}^{2}-\nu\delta\mu\left\langle\delta\chi\right\rangle^{(1)}\log
y+\mathcal{O}\left(y\right),$ (92)
$\displaystyle\left\langle\delta\chi\right\rangle$ $\displaystyle=\sum_{n\geq
0}\left\langle\delta\chi\right\rangle^{(n)}y^{n},$ (93)

with $C_{0}$ an arbitrary constant. However, we choose not to use this
solution in the following since finding the PN expansion for
$\left\langle\delta\chi\right\rangle$ is not trivial in practice, and
integrate Eq. (90) which varies on the radiation reaction timescale.

###  III.3 Effects of the evolution of the elliptic parameter

Through radiation reaction, the coefficients of the cubic equation for
$\delta\chi$ given in Eq. (22) evolve. Therefore, the parameter $m$ of the
Jacobi elliptic function describing the solution does as well. Since this
parameter influences the period $2K(m)$ of this function, the phase $\psi_{p}$
accumulated early on in the evolution will correspond to a different number of
cycles later on when the period has changed, causing a dephasing in the
solution. We can remedy this by not tracking the phase of the Jacobi elliptic
function directly, but a quantity proportional to the corresponding
accumulated number of cycles. For example, we can define $\bar{\psi}_{p}$ such
that

$\displaystyle\mathcal{D}\bar{\psi}_{p}$
$\displaystyle=\frac{\pi}{2K(m)}\mathcal{D}\psi_{p},$ (94)

and recover the correct phase of the Jacobi elliptic function when we need to
evaluate it, with

$\displaystyle\psi_{p}(t)$
$\displaystyle=\frac{2K[m(t)]}{\pi}\bar{\psi}_{p}(t).$ (95)

###  III.4 Averaged spin parameters for the characterization of precessing
binaries

It is useful for the characterization of binary observations through GWs to
describe relevant parameters, the observation of which is indicative of the
presence of some physical effect influencing the radiation. In the case of the
binary components’ spins, one linear combination of the spins’ projection onto
the orbital angular momentum enters the frequency evolution at leading order
in the eccentricity and at lowest PN order. Thus, it forms a natural parameter
to describe the presence of nonzero spins:

$\displaystyle\beta$
$\displaystyle=\frac{94}{113}\chi_{\text{eff}}+\frac{19}{113}\delta\mu\delta\chi.$
(96)

As a sign of the presence of precession, the parameter $\chi_{p}$ was
introduced in [50]. It was recently argued that using a precession averaged
version of this parameter could be even more useful [51]. Following those
works, we define

$\displaystyle\chi_{p}^{2}$
$\displaystyle=\frac{1}{\Omega_{1}^{2}}\left|\mathcal{D}\bm{\hat{L}}\right|^{2},$
(97) $\displaystyle\Omega_{1}$
$\displaystyle=\frac{\mu_{1}y^{6}}{2}[\mu_{1}+3(1-y\chi_{\text{eff}})].$ (98)

Using our solution, we can find an analytic expression for the precession
average of those parameters. We find

$\displaystyle\left\langle\beta\right\rangle$
$\displaystyle=\frac{94}{113}\chi_{\text{eff}}+\frac{19}{113}\delta\mu\left\langle\delta\chi\right\rangle,$
(99)

where $\left\langle\delta\chi\right\rangle$ is given in Eq. (65), and

$\displaystyle\left\langle\left|\mathcal{D}\bm{\hat{L}}\right|^{2}\right\rangle$
$\displaystyle=\frac{y^{10}}{16}\left[c-b\,y\,\delta\mu\left\langle\delta\chi\right\rangle-y^{2}\delta\mu^{2}\left\langle\delta\chi\right\rangle^{2}\right]$
$\displaystyle+\frac{y^{12}\delta\mu^{2}\delta\chi_{\text{diff}}^{2}}{4m^{2}}\bigg{[}\frac{E(m)^{2}}{K(m)^{2}}$
$\displaystyle+\frac{2E(m)}{3K(m)}(m-2)-\frac{m-1}{3}\bigg{]},$ (100)

where the constants $b$ and $c$ are given in appendix C. Note that we were
able to find an analytic expression for $\langle|\chi_{p}|^{2}\rangle^{1/2}$,
where $\chi_{p}$ is defined as the “generalized” parameter introduced in [51],
which slightly differs from $\langle\chi_{p}\rangle$ as was proposed there.

###  III.5 Summary of the full solution

We are now ready to present a full set of evolution equations that vary
exclusively on the radiation reaction timescale, that we can integrate to
construct a full gravitational waveform valid over many precessional cycles of
evolution:

$\displaystyle\mathcal{D}y$ $\displaystyle=y^{9}\sum_{n\geq
0}a_{n}\left(y,e^{2}\right)y^{n},$ (101) $\displaystyle\mathcal{D}e^{2}$
$\displaystyle=y^{8}\sum_{n\geq 0}b_{n}\left(y,e^{2}\right)y^{n},$ (102)
$\displaystyle\mathcal{D}\lambda$ $\displaystyle=y^{3},$ (103)
$\displaystyle\mathcal{D}\delta\lambda$
$\displaystyle=\frac{ky^{3}}{1+k},\quad k=y^{2}\sum_{n\geq
0}k_{n}\left(y,e^{2}\right)y^{n},$ (104) $\displaystyle\mathcal{D}J_{0}$
$\displaystyle=-\frac{L\mathcal{D}y}{2J_{0}y}\left(2L+\chi_{\text{eff}}+\delta\mu\left\langle\delta\chi\right\rangle\right),$
(105) $\displaystyle\mathcal{D}\bar{\psi}_{p}$
$\displaystyle=\frac{3(1-y\chi_{\text{eff}})y^{6}}{4}\frac{\pi}{2K(m)}\sqrt{Y_{3}-Y_{-}},$
(106) $\displaystyle\mathcal{D}\phi_{z,0}$
$\displaystyle=\frac{Jy^{6}}{2}+\frac{3(1-y\chi_{\text{eff}})y^{6}}{4\nu
K(m)}\left(P_{+}+P_{-}\right),$ (107) $\displaystyle\mathcal{D}\zeta_{0}$
$\displaystyle=-\frac{\left(2L+\chi_{\text{eff}}+\delta\mu\left\langle\delta\chi\right\rangle\right)y^{6}}{4}$
$\displaystyle-\frac{3(L+\nu\chi_{\text{eff}})(1-y\chi_{\text{eff}})y^{6}}{2\nu}$
$\displaystyle+\frac{3(1-y\chi_{\text{eff}})y^{6}}{4\nu
K(m)}\left(P_{+}-P_{-}\right),$ (108) $\displaystyle P_{\pm}$
$\displaystyle=\frac{N_{\pm}}{B_{\pm}-C_{\pm}}\Pi\left(\frac{-2C_{\pm}}{B_{\pm}-C_{\pm}},m\right),$
(109)

where we have made use of the averaged spin couplings derived earlier. In all
the expressions for the coefficients of the conservative solution, we further
substitute $J\to J_{0}$. The absence of the spin-precession timescale in these
equations allows one to integrate them efficiently, significantly reducing the
computational cost of an implementation of a solution to the precession
equations, as we will discuss in Sec. V.3.

##  IV Waveform construction

Using our solution to the evolution equations, we can describe the
gravitational radiation emitted in the Fourier domain by expanding the wave
amplitudes for small eccentricities [52, 53, 54, 55]. We first express the
waveform in a fixed polarization basis as

$\displaystyle h_{ab}(t)$
$\displaystyle=h_{+}(t)e_{ab}^{+}+h_{\times}(t)e_{ab}^{\times},$ (110)
$\displaystyle e_{ab}^{+}$ $\displaystyle=\epsilon_{ab}^{+}\cos
2\psi-\epsilon_{ab}^{\times}\sin 2\psi,$ (111) $\displaystyle e_{ab}^{\times}$
$\displaystyle=\epsilon_{ab}^{+}\sin 2\psi+\epsilon_{ab}^{\times}\cos 2\psi,$
(112) $\displaystyle\epsilon_{ab}^{+}$
$\displaystyle=\hat{p}^{a}\hat{p}^{b}-\hat{q}^{a}\hat{q}^{b},$ (113)
$\displaystyle\epsilon_{ab}^{\times}$
$\displaystyle=\hat{p}^{a}\hat{q}^{b}+\hat{q}^{a}\hat{p}^{b},$ (114)

where $\bm{\hat{p}}$ and $\bm{\hat{q}}$ form an orthonormal triad together
with the wave propagation vector $\bm{\hat{k}}$. We can express the wave
polarizations $h_{+}$ and $h_{\times}$ through a rotated spin-weighted
spherical harmonic decomposition. Using an inertial frame aligned with the
constant vector $\bm{\hat{J}}^{(0)}$, we can write [48, 49]

$\displaystyle h_{+}-ih_{\times}$ $\displaystyle=\sum_{l\geq
2}\sum_{m=-l}^{l}H_{lm}Y_{lm}$ (115) $\displaystyle Y_{lm}$
$\displaystyle=\sum_{m^{\prime}=-l}^{l}D_{m^{\prime},m}^{l}(\phi_{z},\theta_{L},\zeta)\left.{}_{-2}Y_{lm^{\prime}}(\theta_{s},\phi_{s})\right.,$
(116)

where $H_{l,-m}=(-1)^{l}H_{lm}^{*}$, $(\theta_{L},\phi_{z})$ are the spherical
angles of $\bm{\hat{L}}$ in the inertial frame, $(\theta_{s},\phi_{s})$ are
the spherical angles of $\bm{\hat{k}}$ in the same frame, and $\zeta$
satisfies $\dot{\zeta}=-\cos\theta_{L}\dot{\phi}_{z}$.
$\left.{}_{-2}Y_{lm^{\prime}}(\theta,\phi)\right.$ are the spin-weighted
spherical harmonics, and $D_{m^{\prime},m}^{l}(\alpha,\beta,\gamma)$ are the
Wigner D-matrices. The conventions used in this work are given in appendix A.

Similarly as in [38], we can expand the wave amplitudes $h_{lm}$ with a
Fourier decomposition in the mean orbital phase $\lambda$ and the argument of
periastron $\delta\lambda$ as

$\displaystyle H_{lm}$
$\displaystyle=\sum_{n\in\mathbb{Z}}G_{lm}^{(m-n)}e^{i(m-n)\delta\lambda}e^{in\lambda},$
(117) $\displaystyle G_{lm}^{(p)}$
$\displaystyle=\mathcal{O}\left(e^{|p|}\right).$ (118)

We can construct an eccentric precessing waveform based on the solution
described in Sec. III.5, along the same principles as laid out in [38]:

$\displaystyle\tilde{h}_{+,\times}(f)=$ $\displaystyle\;\sum_{n\geq
1}\tilde{h}_{n}^{(0)}(f)\tilde{h}_{n}^{\text{PP,$+,\times$}}(f)\,,$ (119)
$\displaystyle\tilde{h}_{n}^{(0)}(f)=$
$\displaystyle\;\sqrt{2\pi}\;T_{n}\exp\Big{[}i\Big{(}2\pi
ft_{n}-n\lambda(t_{n})-\frac{\pi}{4}\Big{)}\Big{]}\,,$ (120) $\displaystyle
2\pi f=$ $\displaystyle\;n\dot{\lambda}(t_{n})\,,$ (121) $\displaystyle
T_{n}=$ $\displaystyle\;\left[n\ddot{\lambda}(t_{n})\right]^{-1/2}\,,$ (122)
$\displaystyle\tilde{h}_{n}^{\text{PP,$+,\times$}}(f)=$
$\displaystyle\;\sum_{p=-P}^{P}e^{i\Delta\Psi_{n,p}}$
$\displaystyle\times\sum_{k=-k_{\text{max}}}^{k_{\text{max}}}a_{k,k_{\text{max}}}\mathcal{A}_{n,p}^{+,\times}(t_{n}+\Delta
t_{n,p}+kT_{n})\,,$ (123)

where

$\displaystyle\Delta t_{n,p}=\;$
$\displaystyle\sum_{q=1}^{P}\frac{1}{q!}\left(-\frac{p}{n}\right)^{q}D^{q-1}\left(\frac{\delta\dot{\lambda}^{q}}{\ddot{\lambda}}\right)\,,$
(124) $\displaystyle\Delta\Psi_{n,p}=\;$
$\displaystyle-p\delta\lambda+n\sum_{q=2}^{P+1}\frac{1}{q!}\left(-\frac{p}{n}\right)^{q}D^{q-2}\left(\frac{\delta\dot{\lambda}^{q}}{\ddot{\lambda}}\right)\,,$
(125) $\displaystyle D$
$\displaystyle=\frac{1}{\ddot{\lambda}}\frac{d}{dt}\,,$ (126)

$P\geq 0$ and $k_{\text{max}}\geq 0$ are arbitrary integers,
$a_{k,k_{\text{max}}}$ satisfy the linear system of equations

$\displaystyle\sum_{k=1}^{k_{\text{max}}}a_{k,k_{\text{max}}}+\frac{1}{2}a_{0,k_{\text{max}}}=1\,,$
(127)
$\displaystyle\sum_{k=1}^{k_{\text{max}}}a_{k,k_{\text{max}}}\frac{k^{2q}}{(2q)!}=\frac{(-i)^{q}}{2^{q}q!}\,,\quad
q\in\\{1,\ldots,k_{\text{max}}\\}\,,$ (128) $\displaystyle
a_{-k,k_{\text{max}}}=a_{k,k_{\text{max}}}\,,$ (129)

and the amplitudes are given by

$\displaystyle\mathcal{A}_{n,p}^{+}$
$\displaystyle=\frac{1}{2}\sum_{l\geq|n-p|}\left.G_{l,p-n}^{(p)}\right.^{*}\left[(-1)^{l}Y_{l,n-p}+Y_{l,p-n}^{*}\right],$
(130) $\displaystyle\mathcal{A}_{n,p}^{\times}$
$\displaystyle=\frac{i}{2}\sum_{l\geq|n-p|}\left.G_{l,p-n}^{(p)}\right.^{*}\left[(-1)^{l}Y_{l,n-p}-Y_{l,p-n}^{*}\right].$
(131)

In particular, we have implemented this waveform to describe inspiralling
binaries observed with LISA. In order to do so, we follow [56] to produce the
rigid adiabatic approximation [57] of the $X$, $Y$, and $Z$ TDI variables [58]
in terms of the Fourier-domain waveform polarizations
$\tilde{h}_{+,\times}(f)$ at the Solar System barycenter given in Eq. (119).

##  V Waveform comparisons

We performed match and computational efficiency comparisons between different
waveforms, including the one presented in this paper and the one derived in
[38] which will serve as a reference. The difference between those two
waveforms is the way in which the equations of motion are computed: in [38]
the equations of precession are integrated directly, resulting in the
precession timescale dictating the computational efficiency of the waveform.
In addition, we include a version of these waveforms with aligned spins for
comparison, and an aligned-spins TaylorF2 [59, 60] circular waveform at 3.5PN
order for computational cost comparisons. We have implemented the LISA
response of these waveforms, and use them as an example for efficiency and
accuracy evaluation. We expect the results to be similar for implementations
applicable to ground-based detectors.

###  V.1 Description of the waveform models

In these two comparison studies, we compare four different waveforms:

  * •

A direct precession solution waveform, hereafter referred to as DPS, for which
we directly solve the equations of precession (1-3), together with Eqs.
(101-104) in order to build the Solar System barycenter Fourier-domain
waveform polarizations $\tilde{h}_{+,\times}(f)$ in Eq.(119).

  * •

A family of efficient fully precessing eccentric waveforms, hereafter referred
to as EFPE-($i,j$) and the main result of the present work that we want to
assess, for which we solve Eqs. (101-105) together with different versions of
Eqs. (106-108) in order to build the Solar System barycenter Fourier-domain
waveform polarizations $\tilde{h}_{+,\times}(f)$ in Eq.(119).

The two integers $i\in\\{0,1\\}$ and $j\in\\{0,1,2\\}$ indicate the version of
Eqs. (106-108) used, thus controlling the accuracy of the $\phi_{z}$ and
$\zeta$ solutions.

When $i=0$, we use an $m=0$ approximation for $\phi_{z,0}$ and $\zeta_{0}$
described in Eqs. (54) and (55)).

When $i=1$, we use the general solutions for $\phi_{z,0}$ and $\zeta_{0}$ in
Eqs. (48) and (51)).

When $j=0$, we ignore the precession corrections $\delta\phi_{z}$ and
$\delta\zeta$ and do not solve Eq. (106).

When $j=1$, we use an $m=0$ approximation for $\delta\phi_{z}$ and
$\delta\zeta$ in Eqs. (56) and (57).

When $j=2$, we use the general solutions for $\delta\phi_{z}$ and
$\delta\zeta$ in Eqs. (49) and (52).

We thus defined here six different EFPE waveforms. However, we feel not
sensible to use the approximate versions of $\phi_{z,0}$ and $\zeta_{0}$ while
using the general solution for $\delta\phi_{z}$ and $\delta\zeta$, and
therefore ignore the waveform EFPE-($0,2$) in the following discussion.

  * •

A spin-aligned eccentric waveform, hereafter referred to as SAE, for which we
solve Eqs. (101-104) without spin-induced precession, in order to build the
Solar System barycenter Fourier-domain waveform polarizations
$\tilde{h}_{+,\times}(f)$ in Eq.(119).

  * •

A circular, spin-aligned, TaylorF2 waveform, hereafter referred to as TF2,
with a standard stationary phase approximation used to build the Solar System
barycenter Fourier-domain waveform polarizations $\tilde{h}_{+,\times}(f)$.

We use the following parameters for the eccentric precessing waveforms:

  * •

$\log(m_{1})$ and $\log(m_{2})$, the logarithm of the redshifted masses of the
binary components.

  * •

$\sin b$ and $l$, respectively the sine ecliptic latitude and longitude of the
source.

  * •

$\sin b_{L}$ and $l_{L}$, respectively the initial sine ecliptic latitude and
initial ecliptic longitude of the orbital angular momentum of the source.

  * •

$\chi_{1}$, $\sin b_{1}$ and $l_{1}$, respectively the dimensionless
magnitude, initial sine ecliptic latitude and initial ecliptic longitude of
the spin of the primary binary component.

  * •

$\chi_{2}$, $\sin b_{2}$ and $l_{2}$, respectively the dimensionless
magnitude, initial sine ecliptic latitude and initial ecliptic longitude of
the spin of the secondary binary component.

  * •

$\lambda_{0}$ and $\delta\lambda_{0}$, respectively the initial mean orbital
phase and initial argument of periastron.

  * •

$\log(e^{2}_{0})$, the logarithm of the initial square orbital eccentricity.

  * •

$t_{M}$, a merger time parameter, defined as the time where the PN parameter
$y$ reaches $y_{M}=6^{-1/2}$ using a leading-order PN solution for $y(t)$ with
constant $e^{2}=e^{2}_{0}$ in order to compute the initial $y(t=0)$. While
this is an approximation of the true merger time, we felt that it was
sufficient for the purposes of this study.

Note that since we only compute normalized matches, the luminosity distance of
the source $d_{L}$ does not affect the results, and we therefore ignore this
parameter. All other parameters follow a flat distribution in our simulations.

$\chi_{1,\ell}$ and $\chi_{2,\ell}$, the spin parameters of the binary
components of the spin-aligned waveforms, are defined as the initial
projections of the individual dimensionless spins onto the orbital angular
momentum of the corresponding precessing waveform. The orbital angular
momentum of the spin-aligned waveform is defined as the initial orbital
angular momentum of the corresponding precessing waveform.

We ran three different sets of simulations, each investigating the efficiency
and accuracy of our waveforms for systems at different stages of their
inspiral. For each system, we start our simulation at $t=0$, corresponding to
the start of the LISA mission, and stop either at the end of the inspiral,
i.e. when the PN parameter $y$ has reached $y=6^{-1/2}$, or when the system
exits the LISA band, i.e. when the orbital frequency $f$ has reached $f=1$ Hz,
whichever comes first. This excludes the merger from our simulations and lets
us concentrate on the inspiral.

In common for each set, we randomize the sine ecliptic latitudes between $-1$
and $1$, the ecliptic longitudes between $0$ and $2\pi$, the spin magnitudes
between $0$ and $1$, $\lambda_{0}$ and $\delta\lambda_{0}$ between $0$ and
$2\pi$, $\log(e_{0}^{2})$ between $\log(10^{-8})$ and $\log(0.1)$, and $t_{M}$
between $2$ and $4$ years.

The three sets are distinguished by the bounds for $\log(m_{1})$ and
$\log(m_{2})$ as

  * •

SOBHB: $10M_{\odot}<m_{1,2}<100M_{\odot}$, a mass distribution corresponding
to stellar-origin black hole binaries.

  * •

IMBHB: $10^{3}M_{\odot}<m_{1,2}<10^{4}M_{\odot}$, a mass distribution
corresponding to intermediate mass black hole binaries.

  * •

MBHB: $10^{5}M_{\odot}<m_{1,2}<10^{6}M_{\odot}$, a mass distribution
corresponding to massive black hole binaries.

Note that despite the fact that the existence of IMBHs lacks supporting
evidence, we feel that inculding this mass distribution is important as the
results are likely to be similar for SOBHBs observed with intermediate band
detectors like DECIGO [61].

Masses | $\mathcal{N}_{\text{orb}}$ | $\mathcal{N}_{\text{spin}}$ | $\mathcal{N}_{\delta\chi}$ | $\mathcal{N}_{\text{ecc}}$
---|---|---|---|---
SOBHB | $1.6\times 10^{6}$ | $1200$ | $530$ | $5800$
IMBHB | $9.6\times 10^{4}$ | $230$ | $100$ | $1100$
MBHB | $5600$ | $42$ | $18$ | $200$
Table 1: For each mass distribution, average number of orbital cycles, spin-
induced precession cycles, $\delta\chi$ cycles, and periastron precession
cycles.

We present in Table 1 the average number of orbital cycles
$\mathcal{N}_{\text{orb}}=\langle\Delta\lambda\rangle/2\pi$, the average
number of spin-induced precession cycles
$\mathcal{N}_{\text{spin}}=\langle\Delta\phi_{z,0}\rangle/2\pi$, the average
number of $\delta\chi$ cycles
$\mathcal{N}_{\delta\chi}=\langle\Delta\bar{\psi}_{p}\rangle/2\pi$, and the
average number of periastron precession cycles
$\mathcal{N}_{\text{ecc}}=\langle\Delta\delta\lambda\rangle/2\pi$ for each
mass distribution. We can see from this table that a faithful modelling of
these sources can be expected to be challenging, especially for SOBHBs which
go on average through over a million orbital cycles, and over a thousand spin
precession cycles.

###  V.2 Faithfulness comparisons

In order to determine the accuracy of the different waveforms, we first define
the match as

$\displaystyle M(h_{1},h_{2})$
$\displaystyle=\frac{(h_{1}|h_{2})}{(h_{1}|h_{1})^{1/2}(h_{2}|h_{2})^{1/2}},$
(132) $\displaystyle(a|b)$ $\displaystyle=4\Re\sum_{k}\int
df\frac{\tilde{a}_{k}(f)\tilde{b}_{k}(f)^{*}}{S_{n,k}(f)},$ (133)

where $k$ denotes one of the three noise-independent TDI channels $A$, $E$ and
$T$, and $S_{n,k}(f)$ is the one-sided power spectral density of the channel
$k$ [62, 63, 64].

The two waveforms $h_{1}$ and $h_{2}$ describe systems with all but two
parameters equal to each other; the waveform $h_{2}$ is allowed to take an
arbitrary initial orbital phase $\lambda_{0}$ and the initial conditions are
defined at an arbitrary initial time $t_{0}$, not necessarily equal to those
defining $h_{1}$. We then define the unfaithfulness as

$\displaystyle U(h_{1},h_{2})$
$\displaystyle=1-\max_{\lambda_{0},t_{0}}M(h_{1},h_{2}).$ (134)

The unfaithfulness can provide a comparison between the modelling bias and the
statistical errors in a parameter estimation simulation through the equation
(see [45] for a derivation)

$\displaystyle E(M)$ $\displaystyle\approx 1-\frac{(D-1)}{2\rho^{2}},$ (135)

where $E(M)$ is the expectation value of the match between a sample from the
posterior distribution and the maximum likelihood sample, $D$ is the
dimensionality of the parameter space, and $\rho=(h|h)^{1/2}$ is the signal-
to-noise ratio (SNR). We can expect the modelling bias to be smaller than the
statistical errors when $U<1-E(M)$. In our case, with $D=17$ we find the
requirement $U<3.6\times 10^{-2}$ for $\rho=15$, $U<1.3\times 10^{-2}$ for
$\rho=25$, and $U<3.2\times 10^{-3}$ for $\rho=50$.

For each mass distribution SOBHB, IMBHB, and MBHB, we simulated 10000 random
systems and computed the unfaithfulness between the reference DPS waveform
$h_{1}$ and different waveforms $h_{2}$, chosen as EFPE-($i,j$) or SAE. To
maximize over phase and time shifts, we used a Newton-Raphson method. As a
convergence criterion, we stopped the iteration when
$|U^{(i+1)}-U^{(i)}|<10^{-8}$, or when $|t_{0}^{(i+1)}-t_{0}^{(i)}|<0.01$ s.

Masses | $h_{2}$ | $q_{0.1}$ | $q_{0.5}$ | $q_{0.9}$
---|---|---|---|---
SOBHB | SAE | $0.77$ | $0.19$ | $8.2\times 10^{-3}$
SOBHB | EFPE-($0,0$) | $7.6\times 10^{-4}$ | $5.0\times 10^{-5}$ | $2.3\times 10^{-6}$
SOBHB | EFPE-($0,1$) | $5.0\times 10^{-4}$ | $5.6\times 10^{-6}$ | $3.8\times 10^{-8}$
SOBHB | EFPE-($1,0$) | $4.7\times 10^{-4}$ | $4.8\times 10^{-5}$ | $2.3\times 10^{-6}$
SOBHB | EFPE-($1,1$) | $1.4\times 10^{-4}$ | $2.8\times 10^{-6}$ | $4.0\times 10^{-8}$
SOBHB | EFPE-($1,2$) | $1.4\times 10^{-4}$ | $2.8\times 10^{-6}$ | $4.0\times 10^{-8}$
IMBHB | SAE | $0.89$ | $0.16$ | $3.0\times 10^{-3}$
IMBHB | EFPE-($0,0$) | $3.4\times 10^{-3}$ | $1.4\times 10^{-4}$ | $3.4\times 10^{-6}$
IMBHB | EFPE-($0,1$) | $2.5\times 10^{-3}$ | $9.9\times 10^{-6}$ | $8.0\times 10^{-8}$
IMBHB | EFPE-($1,0$) | $1.2\times 10^{-3}$ | $1.1\times 10^{-4}$ | $3.4\times 10^{-6}$
IMBHB | EFPE-($1,1$) | $2.3\times 10^{-5}$ | $7.0\times 10^{-7}$ | $5.4\times 10^{-8}$
IMBHB | EFPE-($1,2$) | $1.8\times 10^{-5}$ | $6.7\times 10^{-7}$ | $5.3\times 10^{-8}$
MBHB | SAE | $0.91$ | $0.13$ | $4.4\times 10^{-3}$
MBHB | EFPE-($0,0$) | $1.3\times 10^{-2}$ | $7.5\times 10^{-4}$ | $3.6\times 10^{-5}$
MBHB | EFPE-($0,1$) | $7.4\times 10^{-3}$ | $1.0\times 10^{-4}$ | $3.7\times 10^{-6}$
MBHB | EFPE-($1,0$) | $6.7\times 10^{-3}$ | $6.8\times 10^{-4}$ | $3.6\times 10^{-5}$
MBHB | EFPE-($1,1$) | $9.6\times 10^{-4}$ | $9.3\times 10^{-5}$ | $4.0\times 10^{-6}$
MBHB | EFPE-($1,2$) | $8.5\times 10^{-4}$ | $9.0\times 10^{-5}$ | $4.0\times 10^{-6}$
Table 2: For each mass distribution and each test waveform $h_{2}$, 10th
percentile $q_{0.1}$, median $q_{0.5}$, and 90th percentile $q_{0.9}$ of the
distributions of the unfaithfulness $U(h_{1},h_{2})$ with the reference direct
precession solution (DPS) waveform $h_{1}$, for the spin-aligned eccentric
(SAE) and the different efficient full precessing eccentric (EFPE) waveforms.
Figure 3: Distributions of the unfaithfulness $U(h_{1},h_{2})$ between the
reference direct precession solution (DPS) waveform $h_{1}$ and the test
waveforms $h_{2}$ SAE (spin-aligned eccentric, solid blue line), EFPE-($0,0$)
(fastest efficient full precessing eccentric, dashed orange line), and
EFPE-($1,2$) (most accurate efficient full precessing eccentric, dotted green
line), for SOBHB systems. Figure 4: Distributions of the unfaithfulness
$U(h_{1},h_{2})$ between the reference direct precession solution (DPS)
waveform $h_{1}$ and the test waveforms $h_{2}$ SAE (spin-aligned eccentric,
solid blue line), EFPE-($0,0$) (fastest efficient full precessing eccentric,
dashed orange line), and EFPE-($1,2$) (most accurate efficient full precessing
eccentric, dotted green line), for IMBHB systems. Figure 5: Distributions of
the unfaithfulness $U(h_{1},h_{2})$ between the reference direct precession
solution (DPS) waveform $h_{1}$ and the test waveforms $h_{2}$ SAE (spin-
aligned eccentric, solid blue line), EFPE-($0,0$) (fastest efficient full
precessing eccentric, dashed orange line), and EFPE-($1,2$) (most accurate
efficient full precessing eccentric, dotted green line), for MBHB systems.

We show in Table 2 the 10th percentile $q_{0.1}$, median $q_{0.5}$, and 90th
percentile $q_{0.9}$ of the resulting unfaithfulness distributions, and we
show the unfaithfulness distributions for the SAE, the EFPE-($0,0$), and the
EFPE-($1,2$) waveforms in Fig. 3 for the SOBHB systems, in Fig. 4 for the
IMBHB systems, and in Fig. 5 for the MBHB systems. From these quantiles, we
can note that the typical unfaithfulness for the SAE waveforms is above 0.1
for all mass distributions, which is well above the requirement even for an
SNR of 15. The distributions have a long tail with 90th percentile below
$10^{-2}$, which we postulate is due to systems with small spin projections
onto the orbital plane, relatively inaffected by spin-induced precession
effects. The EFPE-($i,j$) waveforms, on the other hand, all pass the
unfaithfulness requirement for an SNR of at least 25, even when considering
the 10th percentile. We also see that the inclusion of the precession
corrections $\delta\phi_{z}$ and $\delta\zeta$ decreases the median
unfaithfulness by a factor of $\sim 10$ for most waveforms, and by a factor of
$\sim 100$ for the EFPE-($1,j$) waveforms modelling IMBHB systems. Comparing
the different EFPE-($i,j$) waveforms, we see little difference between
EFPE-($0,0$) and EFPE-($1,0$), as well as between EFPE-($1,1$) and
EFPE-($1,2$). For IMBHB systems, we do see an improvement in the median
unfaithfulness by a factor of $\sim 10$ from EFPE-($0,1$) to EFPE-($1,1$).
Using the median results, the best performing EFPE-($1,2$) waveform passes the
unfaithfulness requirement for an SNR of $\rho\lesssim 1400$ for SOBHB
systems, $\rho\lesssim 3500$ for IMBHB systems, and $\rho\lesssim 300$ for
MBHB systems, while the simplest EFPE-($0,0$) for an SNR of $\rho\lesssim 400$
for SOBHB systems, and $\rho\lesssim 100$ for IMBHB and MBHB systems. We
verified that the median unfaithfulnesses were independent on the initial
eccentricity.

We stress that the waveforms used in this study are identical, except for the
treatment of spin-induced precession. The goal of these comparisons is to
evaluate the amount of mismatch appearing in our waveform due to this
difference, while other effects, such as the truncation of PN series, the
approximation of the Fourier transform, finite armlength effects, etc., will
be additional sources of mismatch.

###  V.3 Computational efficiency comparisons

For each mass distribution, we ran a set of simulations aimed at evaluating
the computational efficiency of each waveform model. In order to do so, we
pre-randomized $10^{5}$ sets of binary parameters and started evaluating each
waveform model using these. We timed the waveform evaluations and stopped
after 20 minutes, or when all waveforms had been evaluated, whichever came
first. A waveform evaluation consisted in solving the equations of motion if
applicable, and evaluate the waveform on a fixed set of 301 frequencies,
equally spaced in log between $10$ mHz and $1$ Hz for SOBHBs, between $1$ mHz
and $1$ Hz for IMBHBs, and between $0.05$ and $20$ mHz for MBHBs. We chose
fixed grids in order to reduce to a maximum the computation time not used to
construct waveforms during the timed evaluation, and set their bounds so that
all simulated waves emitted in a significant portion of the window.

| SOBHB | IMBHB | MBHB
---|---|---|---
DPS | 449 | 71.4 | 13.3
TF2 | 1.06 | 0.64 | 0.44
SAE | 5.22 | 2.00 | 1.23
EFPE-($0,0$) | 10.5 | 4.23 | 2.62
EFPE-($0,1$) | 11.0 | 4.63 | 2.84
EFPE-($1,0$) | 23.9 | 8.15 | 4.86
EFPE-($1,1$) | 24.5 | 8.53 | 5.03
EFPE-($1,2$) | 27.5 | 10.3 | 6.48
Table 3: Average evaluation time in milliseconds for each waveform and for
each mass distribution. In the SOBHB and IMBHB cases, the reference direct
precession solution waveform (DPS) is too computationally expensive for
Bayesian parameter estimation, unlike the EFPE-($i,j$) models which still
capture precession and eccentricity effects. In the MBHB case, while the DPS
waveform is fast enough, the EFPE-($i,j$) models still provide a speed-up of a
factor 2-5.

We show in Table 3 the results of the timed evaluations. We can notice that
the computational costs of the EFPE-($i,j$) models depend very weakly on the
treatment of $\delta\phi_{z}$ and $\delta\zeta$. We find that the EFPE-($1,j$)
models are about twice as costly as the EFPE-($0,j$) models, and the latter
are only about twice as costly as the SAE model. We expect the computational
cost of the TF2 model to be largely independent of the mass distribution,
which can provide us with an estimate of the accuracy of the timing estimates.
Comparing to the DPS model, while the EFPE-($0,0$) is about a factor of five
more efficient for MBHBs, the difference increases as the length of the
waveform increases, being about a factor of $20$ for IMBHBs, and a factor of
$50$ for SOBHBs. While an average evaluation time of $449$ ms renders the DPS
model impractical for Bayessian parameter estimation of SOBHB systems, the
EFPE-($i,j$) models with an average of $10$-$25$ ms per evaluation makes it
feasible. Finally, we note that the results of the TF2 waveform leads us to
believe that developing a purely analytic version of the waveform we presented
here might still lead to substantial computational efficiency improvements.

##  VI Further extensions of the solution

The solutions presented here are based on a few assumptions that it might be
preferable to relax in some situations. For example, when extended objects
such as neutron stars are considered, the equations of precession are modified
in such a way that the effective spin $\chi_{\text{eff}}$ is no longer
conserved [22, 65]. Furthermore, the solution presented in this work is based
on neglecting certain terms in the spin precession equations, such as next-to-
leading order spin-orbit and spin-spin terms [66, 67], or cubic-in-spin terms
[68]. The extent to which this approximation is valid remains to be assessed,
and to do so we need to construct a waveform that includes those terms. The
present work can be used as a basis to accomplish such tasks, since such
effects appear at higher PN order than the one we included. In the following,
we briefly sketch a possible avenue to construct such an extension, by using
quadrupole-monopole terms [22, 65] as an example.

The equations of precession acquire a modification when applied to extended
bodies:

$\bm{\Omega}_{i}=\left\\{\frac{1}{2}\mu_{i}+\frac{3}{2}\left[1-y\bm{\hat{L}}\cdot\left(\kappa_{i}\bm{s}_{i}+\bm{s}_{j}\right)\right]\right\\}\bm{\hat{L}}\times\bm{s}_{i}\\\
+\frac{1}{2}y\bm{s}_{j}\times\bm{s}_{i},$ (136)

where $\kappa_{i}$ is the quadrupole parameter of body $i$. In the case of
black holes, we have $\kappa_{i}=1$, so that we can write

$\displaystyle\bm{\Omega}_{i}$
$\displaystyle=\bm{\Omega}_{i}^{\text{BH}}-\frac{3}{2}y\delta\kappa_{i}\left(\bm{\hat{L}}\cdot\bm{s}_{i}\right)\bm{\hat{L}}\times\bm{s}_{i},$
(137)

with $\delta\kappa_{i}=\kappa_{i}-1$.

Eqs. (14-21) being purely geometrical, they are not affected by a modification
of the precession equations. On the other hand, we find as the derivatives of
$\delta\chi$ and $\chi_{\text{eff}}$

$\displaystyle\mathcal{D}\delta\chi$
$\displaystyle=6y^{6}\left[1-y\chi_{\text{eff}}-\frac{y}{2}\left(\delta\kappa_{1}\chi_{1}+\delta\kappa_{2}\chi_{2}\right)\right]V,$
(138) $\displaystyle\mathcal{D}\chi_{\text{eff}}$
$\displaystyle=3y^{7}\left[\delta\kappa_{1}\chi_{1}-\delta\kappa_{2}\chi_{2}\right]V,$
(139)

where we defined $\chi_{i}=\bm{\hat{L}}\cdot\bm{s}_{i}$ and
$V=\left(\bm{\hat{L}}\times\bm{s}_{1}\right)\cdot\bm{s}_{2}$.

We can observe that

$\displaystyle\mathcal{D}\chi_{\text{eff}}$
$\displaystyle=\mathcal{O}\left(y\right)\times\mathcal{D}\delta\chi.$ (140)

Furthermore, since the form of Eq. (22) follows from the unmodified Eqs.
(14-21), we can still write

$\displaystyle\left(\mathcal{D}\delta\chi\right)^{2}$
$\displaystyle=\frac{9}{4}A^{2}y^{11}\left(\delta\mu\delta\chi^{3}+B\delta\chi^{2}+C\delta\chi+D\right),$
(141)

with

$\displaystyle A$
$\displaystyle=1-y\chi_{\text{eff}}-\frac{y}{2}(\delta\kappa_{1}\chi_{1}+\delta\kappa_{2}\chi_{2}),$
(142)

and $B$, $C$, and $D$ unmodified with respect to the black hole solutions.

This presents a few complications from the fact that $A$, $B$, $C$, and $D$
are no longer constants. However, Eqs. (140), (142), and (155-157) lead us to
postulate that if we write

$\displaystyle\chi_{\text{eff}}$ $\displaystyle=\sum_{k\geq
0}\chi_{\text{eff}}^{(k)}y^{k},$ (143) $\displaystyle\delta\chi$
$\displaystyle=\sum_{k\geq 0}\delta\chi^{(k)}y^{k},$ (144)

we can solve order by order by using for $\chi_{\text{eff}}^{(0)}$ and
$\delta\chi^{(0)}$ the solutions of the black hole case derived in Sec. II.
The new solutions that we will then derive for $\chi_{\text{eff}}^{(1)}$ and
$\delta\chi^{(1)}$ will be linear in $\delta\kappa_{i}$, and so on. Therefore,
we can in principle derive new solutions as perturbations of the solution
presented in this paper, still varying on the radiation reaction timescale
only. We argue that this will continue to be true as long as the modifications
of the precession equations can be written as perturbations of Eqs. (1-3).

We stress that we intend the discussion presented in this section as a mere
sketch of a possible avenue for including a wide range of extra effects into
our solution. If more work is needed to verify its applicability, we believe
that its importance makes it worth mentioning.

##  VII Conclusion

In this work, we have constructed a new class of gravitational waveform models
for inspiralling precessing binary systems on eccentric orbits, by taking
advantage of an analytic solution of the conservative problem. This solution
allows for the modification of the complete set of equations of motion, such
that it only varies on the radiation reaction timescale. This permits the
construction of eccentric precessing waveforms with computational costs
comparable to ones valid for spin-aligned systems. While we observe a modest
efficiency gain of a factor $2$-$5$ for binaries in the late inspiral, we
observe a substantial gain of a factor $20$-$50$ for binaries in the early
inspiral such as stellar-origin binaries observed with LISA. This improvement
comes at a modest cost in accuracy, not likely to lead to unwanted biases in
parameter recovery for systems with signal-to-noise ratios lower than several
hundreds. This solution opens the way for Bayesian parameter estimation
studies of eccentric, precessing stellar-origin binaries with LISA, an
important class of sources for the study of black hole binary formation
mechanisms. Such waveforms have been implemented to model the LISA time-delay
interferometry response, and have been shown to require a computational cost
of about $10$-$30$ ms per waveform evaluation, well inside the range of
possibilities for Bayesian parameter estimation studies.

We point to possible extensions of this model based on perturbations of the
spin-induced precession equations. Important such extensions include
quadrupole-monopole interactions relevant to the modelling of binaries
involving neutron stars, higher post-Newtonian order spin interactions, and
higher order in spin interactions. While these interactions are likely to have
relatively small effects on the precession physics, they can still be relevant
for stellar-origin binaries in the LISA band as they acquire thousands of
spin-induced precession cycles. Studies determining their importance therefore
need to be conducted, and we point to future work for their implementation.

###### Acknowledgements.

We would like to thank Karthik Srinivasan, Chandra Mishra and Guillaume Faye
for interesting discussions that inspired this investigation.

##  Appendix A Conventions

In this paper, we use the following conventions for the elliptic integrals:

$\displaystyle F(\phi;m)$
$\displaystyle=\int_{0}^{\phi}\frac{dt}{\sqrt{1-m\sin^{2}t}},$ (145)
$\displaystyle E(\phi;m)$
$\displaystyle=\int_{0}^{\phi}dt\sqrt{1-m\sin^{2}t},$ (146)
$\displaystyle\Pi(n;\phi;m)$
$\displaystyle=\int_{0}^{\phi}\frac{dt}{\left(1-n\sin^{2}t\right)\sqrt{1-m\sin^{2}t}},$
(147) $\displaystyle K(m)$ $\displaystyle=F\left(\frac{\pi}{2};m\right),$
(148) $\displaystyle E(m)$ $\displaystyle=E\left(\frac{\pi}{2};m\right),$
(149) $\displaystyle\Pi(n,m)$
$\displaystyle=\Pi\left(n;\frac{\pi}{2};m\right),$ (150)
$\displaystyle\text{am}(\psi,m)$
$\displaystyle=\arcsin[\text{sn}(\psi,m)],\quad-K(m)<\psi<K(m).$ (151)

We use the following conventions for the spin-weighted spherical harmonics and
the Wigner D-matrices:

$\displaystyle\left.{}_{s}Y_{lm}(\theta,\phi)\right.$
$\displaystyle=(-1)^{s}e^{im\phi}\sqrt{\frac{2l+1}{4\pi}}\sqrt{(l+m)!(l-m)!(l+s)!(l-s)!}$
$\displaystyle\times\sum_{k=\max(0,m+s)}^{\min(l+m,l+s)}\frac{(-1)^{k}\sin^{2k-s-m}\left(\frac{\theta}{2}\right)\cos^{2(l-k)+m+s}\left(\frac{\theta}{2}\right)}{k!(k-m-s)!(l+m-k)!(l+s-k)!},$
(152) $\displaystyle D^{l}_{m^{\prime},m}(\alpha,\beta,\gamma)$
$\displaystyle=(-1)^{m+m^{\prime}}e^{-i(m^{\prime}\alpha+m\gamma)}\sqrt{(l+m^{\prime})!(l-m^{\prime})!(l+m)!(l-m)!}$
$\displaystyle\times\sum_{k=\max(0,m-m^{\prime})}^{\min(l+m,l-m^{\prime})}\frac{(-1)^{k}\sin^{2k+m^{\prime}-m}\left(\frac{\beta}{2}\right)\cos^{2(l-k)+m-m^{\prime}}\left(\frac{\beta}{2}\right)}{k!(k-m+m^{\prime})!(l+m-k)!(l-m^{\prime}-k)!}.$
(153)

We have the following relation:

$\displaystyle\left.{}_{s}Y_{lm}(\theta,\phi)\right.$
$\displaystyle=(-1)^{m}\sqrt{\frac{2l+1}{4\pi}}e^{-is\psi}D^{l}_{-m,s}(\phi,\theta,-\psi)$
(154)

##  Appendix B Discussion of the $\delta\chi$ solution

The coefficients of Eq. (22) are given by

$\displaystyle B$
$\displaystyle=\frac{y}{2\nu^{2}}\left[-2\nu\left(J^{2}-L^{2}-L\chi_{\text{eff}}\right)+\delta\mu\left(S_{1}^{2}-S_{2}^{2}\right)-\delta\mu^{2}\left(2L^{2}+S_{1}^{2}+S_{2}^{2}\right)\right],$
(155) $\displaystyle C$
$\displaystyle=\frac{y}{2\nu^{2}}\left\\{\left(1+\delta\mu^{2}\right)\chi_{\text{eff}}\left(S_{1}^{2}-S_{2}^{2}\right)+2\delta\mu\left[2L\left(J^{2}-L^{2}-L\chi_{\text{eff}}\right)-(2L+\chi_{\text{eff}})\left(S_{1}^{2}+S_{2}^{2}\right)-\nu
L\chi_{\text{eff}}^{2}\right]\right\\},$ (156) $\displaystyle D$
$\displaystyle=\frac{y}{2\nu^{2}}\big{\\{}-2\left(J^{2}-L^{2}-L\chi_{\text{eff}}\right)\left[J^{2}-L^{2}-L\chi_{\text{eff}}-2\left(S_{1}^{2}+S_{2}^{2}\right)-\nu\chi_{\text{eff}}^{2}\right]$
$\displaystyle+\left(S_{1}^{2}-S_{2}^{2}\right)\left[\delta\mu\chi_{\text{eff}}^{2}-2\left(S_{1}^{2}-S_{2}^{2}\right)\right]-\chi_{\text{eff}}^{2}\left(S_{1}^{2}+S_{2}^{2}\right)\big{\\}}.$
(157)

We can verify that $\delta\chi_{\pm}$ are indeed regular in the equal-mass
limit. Expanding the solutions of the cubic polynomial in terms of the
coefficients $B$, $C$, and $D$ for small $\delta\mu$, we get

$\displaystyle\delta\chi_{3}$
$\displaystyle=-B+\mathcal{O}\left(\delta\mu\right),$ (158)
$\displaystyle\delta\chi_{\pm}$
$\displaystyle=-\frac{C\pm\sqrt{C^{2}-4BD}}{2B}+\mathcal{O}\left(\delta\mu\right).$
(159)

We can thus see that $\delta\chi_{\pm}$ are regular in the equal-mass limit,
and we recognize the quadratic formula in the solutions at the limit. We can
further find

$\displaystyle m$
$\displaystyle=\frac{\sqrt{C^{2}-4BD}}{B^{2}}\delta\mu+\mathcal{O}\left(\delta\mu^{2}\right),$
(160) $\displaystyle\sqrt{Y_{3}-Y_{-}}$
$\displaystyle=\sqrt{-\frac{B}{y}}+\mathcal{O}\left(\delta\mu\right)$ (161)
$\displaystyle=\sqrt{\frac{1}{\nu}\left[J^{2}-L(L+\chi_{\text{eff}})\right]}+\mathcal{O}\left(\delta\mu\right).$
(162)

We can also PN expand the coefficients of the cubic polynomial. In order to do
that, we should realize that $J=L+\mathcal{O}\left(y^{0}\right)$. We can thus
write

$\displaystyle J$ $\displaystyle=L+\sum_{n\geq 0}j_{n}y^{n},$ (163)

and find

$\displaystyle\delta\chi_{\pm}$
$\displaystyle=\frac{1}{\delta\mu}\left(2j_{0}-\chi_{\text{eff}}\right)+\mathcal{O}\left(y\right),$
(164) $\displaystyle m$ $\displaystyle=\mathcal{O}\left(y^{2}\right),$ (165)
$\displaystyle\mathcal{D}\psi_{\text{p}}$ $\displaystyle=\frac{3\delta\mu
y^{5}}{2}+\mathcal{O}\left(y^{6}\right),$ (166)
$\displaystyle\delta\chi_{\text{diff}}$
$\displaystyle=\mathcal{O}\left(y\right),$ (167)
$\displaystyle\left\langle\mathcal{D}\phi_{z}\right\rangle$
$\displaystyle=\frac{y^{5}(7-\delta\mu^{2})}{8}+\mathcal{O}\left(y^{6}\right).$
(168)

Therefore, we expect the parameter $m$ to be small in the large separation
limit, and to increase as the binary components get closer and closer
together. Note that, while the solution that we presented is perfectly regular
in the equal-mass limit, the leading order PN term for $\mathcal{D}\psi_{p}$
vanishes in this limit, rendering it difficult to analytically integrate, as
was mentioned in [45].

##  Appendix C Coefficients entering the average $\chi_{p}$

The coefficients of Eq. (100) are given by

$\displaystyle b$
$\displaystyle=49-70y\chi_{\text{eff}}+24y^{2}\chi_{\text{eff}}^{2}-\delta\mu^{2},$
(169) $\displaystyle c$
$\displaystyle=\nu\frac{J^{2}}{L^{2}}\left[\left(7-6y\chi_{\text{eff}}\right)^{2}-\delta\mu^{2}\right]-\frac{1}{4}\left(7+8y\chi_{\text{eff}}-12y^{2}\chi_{\text{eff}}^{2}-\delta\mu^{2}\right)^{2}+3\delta\mu^{2}\left(3-8y\chi_{\text{eff}}+5y^{2}\chi_{\text{eff}}^{2}\right)$
$\displaystyle-6\delta\mu
y^{2}(1-y\chi_{\text{eff}})\left[\left(7-6y\chi_{\text{eff}}+\delta\mu\right)\frac{s_{1}^{2}}{\mu_{2}}-\left(7-6y\chi_{\text{eff}}-\delta\mu\right)\frac{s_{2}^{2}}{\mu_{1}}\right].$
(170)

## References

  * Abbott  _et al._ [2018] B. P. Abbott _et al._ , Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Rev. Rel. 21, 3 (2018), arXiv:1304.0670 [gr-qc] .
  * LIGO Scientific Collaboration [2015] LIGO Scientific Collaboration, Advanced LIGO, Classical and Quantum Gravity 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ , Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Somiya [2012] K. Somiya, Detector configuration of KAGRA-the Japanese cryogenic gravitational-wave detector, Class. Quant. Grav. 29, 124007 (2012), arXiv:1111.7185 [gr-qc] .
  * Aso  _et al._ [2013] Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto, Interferometer design of the KAGRA gravitational wave detector, Phys. Rev. D 88, 043007 (2013), arXiv:1306.6747 [gr-qc] .
  * Abbott  _et al._ [2019a] B. P. Abbott _et al._ , Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo, Astrophys. J. Lett. 882, L24 (2019a), arXiv:1811.12940 [astro-ph.HE] .
  * Abbott  _et al._ [2019b] B. P. Abbott _et al._ , GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019b), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Gwtc-2: Compact binary coalescences observed by ligo and virgo during the first half of the third observing run, Phys. Rev. X 11, 021053 (2021), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ , GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2019c] B. P. Abbott _et al._ , Tests of general relativity with the binary black hole signals from the LIGO-Virgo catalog GWTC-1, Phys. Rev. D 100, 104036 (2019c), arXiv:1903.04467 [gr-qc] .
  * Perera  _et al._ [2019] B. B. P. Perera _et al._ , The International Pulsar Timing Array: second data release, Mon. Not. R. Astro. Soc. 490, 4666 (2019), arXiv:1909.04534 [astro-ph.HE] .
  * Arzoumanian  _et al._ [2020] Z. Arzoumanian _et al._ , The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background, Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Bian  _et al._ [2021] L. Bian, R.-G. Cai, J. Liu, X.-Y. Yang, and R. Zhou, Evidence for different gravitational-wave sources in the NANOGrav dataset, Phys. Rev. D 103, L081301 (2021), arXiv:2009.13893 [astro-ph.CO] .
  * Middleton  _et al._ [2021] H. Middleton, A. Sesana, S. Chen, A. Vecchio, W. Del Pozzo, and P. A. Rosado, Massive black hole binary systems and the NANOGrav 12.5 yr results, Mon. Not. R. Astro. Soc. 502, L99 (2021), arXiv:2011.01246 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ , Laser Interferometer Space Antenna, arXiv e-prints , arXiv:1702.00786 (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Klein  _et al._ [2016] A. Klein, E. Barausse, A. Sesana, A. Petiteau, E. Berti, S. Babak, J. Gair, S. Aoudia, I. Hinder, F. Ohme, and B. Wardell, Science with the space-based interferometer eLISA: Supermassive black hole binaries, Phys. Rev. D 93, 024003 (2016), arXiv:1511.05581 [gr-qc] .
  * Mangiagli  _et al._ [2020] A. Mangiagli, A. Klein, M. Bonetti, M. L. Katz, A. Sesana, M. Volonteri, M. Colpi, S. Marsat, and S. Babak, Observing the inspiral of coalescing massive black hole binaries with LISA in the era of multimessenger astrophysics, Phys. Rev. D 102, 084056 (2020), arXiv:2006.12513 [astro-ph.HE] .
  * Babak  _et al._ [2017] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Science with the space-based interferometer LISA. V. Extreme mass-ratio inspirals, Phys. Rev. D 95, 103012 (2017), arXiv:1703.09722 [gr-qc] .
  * Nissanke  _et al._ [2012] S. Nissanke, M. Vallisneri, G. Nelemans, and T. A. Prince, Gravitational-wave Emission from Compact Galactic Binaries, Astrophys. J. 758, 131 (2012), arXiv:1201.4613 [astro-ph.GA] .
  * Korol  _et al._ [2017] V. Korol, E. M. Rossi, P. J. Groot, G. Nelemans, S. Toonen, and A. G. A. Brown, Prospects for detection of detached double white dwarf binaries with Gaia, LSST and LISA, Mon. Not. R. Astro. Soc. 470, 1894 (2017), arXiv:1703.02555 [astro-ph.HE] .
  * Sesana [2016] A. Sesana, Prospects for Multiband Gravitational-Wave Astronomy after GW150914, Phys. Rev. Lett. 116, 231102 (2016), arXiv:1602.06951 [gr-qc] .
  * Barker and Oconnell [1979] B. M. Barker and R. F. Oconnell, The gravitational interaction: spin, rotation, and quantum effects - a review., Gen. Rel. Grav. 11, 149 (1979).
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin-induced orbital precession and its modulation of the gravitational waveforms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * Hinderer  _et al._ [2010] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral, Phys. Rev. D 81, 123016 (2010), arXiv:0911.3535 [astro-ph.HE] .
  * Peters and Mathews [1963] P. C. Peters and J. Mathews, Gravitational Radiation from Point Masses in a Keplerian Orbit, Phys. Rev. 131, 435 (1963).
  * Damour  _et al._ [2004] T. Damour, A. Gopakumar, and B. R. Iyer, Phasing of gravitational waves from inspiralling eccentric binaries, Phys. Rev. D 70, 064028 (2004), arXiv:gr-qc/0404128 [astro-ph] .
  * Nishizawa  _et al._ [2017] A. Nishizawa, A. Sesana, E. Berti, and A. Klein, Constraining stellar binary black hole formation scenarios with eLISA eccentricity measurements, Mon. Not. R. Astro. Soc. 465, 4375 (2017), arXiv:1606.09295 [astro-ph.HE] .
  * Zevin  _et al._ [2021] M. Zevin, I. M. Romero-Shaw, K. Kremer, E. Thrane, and P. D. Lasky, Implications of Eccentric Observations on Binary Black Hole Formation Channels, arXiv e-prints , arXiv:2106.09042 (2021), arXiv:2106.09042 [astro-ph.HE] .
  * Bonetti  _et al._ [2019] M. Bonetti, A. Sesana, F. Haardt, E. Barausse, and M. Colpi, Post-Newtonian evolution of massive black hole triplets in galactic nuclei - IV. Implications for LISA, Mon. Not. R. Astro. Soc. 486, 4044 (2019), arXiv:1812.01011 [astro-ph.GA] .
  * Lundgren and O’Shaughnessy [2014] A. Lundgren and R. O’Shaughnessy, Single-spin precessing gravitational waveform in closed form, Phys. Rev. D 89, 044021 (2014), arXiv:1304.3332 [gr-qc] .
  * Hannam  _et al._ [2014] M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Khan  _et al._ [2019] S. Khan, K. Chatziioannou, M. Hannam, and F. Ohme, Phenomenological model for the gravitational-wave signal from precessing binary black holes with two-spin effects, Phys. Rev. D 100, 024059 (2019), arXiv:1809.10113 [gr-qc] .
  * Taracchini  _et al._ [2014] A. Taracchini, A. Buonanno, Y. Pan, T. Hinderer, M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, B. Szilágyi, N. W. Taylor, and A. Zenginoglu, Effective-one-body model for black-hole binaries with generic mass ratios and spins, Phys. Rev. D 89, 061502 (2014), arXiv:1311.2544 [gr-qc] .
  * Hinderer and Babak [2017] T. Hinderer and S. Babak, Foundations of an effective-one-body model for coalescing binaries on eccentric orbits, Phys. Rev. D 96, 104048 (2017), arXiv:1707.08426 [gr-qc] .
  * Huerta  _et al._ [2018] 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, Eccentric, nonspinning, inspiral, Gaussian-process merger approximant for the detection and characterization of eccentric binary black hole mergers, Phys. Rev. D 97, 024031 (2018), arXiv:1711.06276 [gr-qc] .
  * Yun  _et al._ [2021] Q. Yun, W.-B. Han, X. Zhong, and C. A. Benavides-Gallego, A surrogate model for gravitational waveforms of spin-aligned binary black holes with eccentricities, arXiv e-prints , arXiv:2104.03789 (2021), arXiv:2104.03789 [gr-qc] .
  * Cornish and Key [2010] N. J. Cornish and J. S. Key, Computing waveforms for spinning compact binaries in quasi-eccentric orbits, Phys. Rev. D 82, 044028 (2010), arXiv:1004.5322 [gr-qc] .
  * Klein  _et al._ [2018] A. Klein, Y. Boetzel, A. Gopakumar, P. Jetzer, and L. de Vittori, Fourier domain gravitational waveforms for precessing eccentric binaries, Phys. Rev. D 98, 104043 (2018), arXiv:1801.08542 [gr-qc] .
  * Kesden  _et al._ [2015] M. Kesden, D. Gerosa, R. O’Shaughnessy, E. Berti, and U. Sperhake, Effective Potentials and Morphological Transitions for Binary Black Hole Spin Precession, Phys. Rev. Lett. 114, 081103 (2015), arXiv:1411.0674 [gr-qc] .
  * Gerosa  _et al._ [2015] D. Gerosa, M. Kesden, U. Sperhake, E. Berti, and R. O’Shaughnessy, Multi-timescale analysis of phase transitions in precessing black-hole binaries, Phys. Rev. D 92, 064016 (2015), arXiv:1506.03492 [gr-qc] .
  * Buscicchio  _et al._ [2021] R. Buscicchio, A. Klein, E. Roebber, C. J. Moore, D. Gerosa, E. Finch, and A. Vecchio, Bayesian parameter estimation of stellar-mass black-hole binaries with LISA, arXiv e-prints , arXiv:2106.05259 (2021), arXiv:2106.05259 [astro-ph.HE] .
  * Klein  _et al._ [2013] A. Klein, N. Cornish, and N. Yunes, Gravitational waveforms for precessing, quasicircular binaries via multiple scale analysis and uniform asymptotics: The near spin alignment case, Phys. Rev. D 88, 124015 (2013), arXiv:1305.1932 [gr-qc] .
  * Chatziioannou  _et al._ [2013] K. Chatziioannou, A. Klein, N. Yunes, and N. Cornish, Gravitational waveforms for precessing, quasicircular compact binaries with multiple scale analysis: Small spin expansion, Phys. Rev. D 88, 063011 (2013), arXiv:1307.4418 [gr-qc] .
  * Chatziioannou  _et al._ [2017a] K. Chatziioannou, A. Klein, N. Cornish, and N. Yunes, Analytic Gravitational Waveforms for Generic Precessing Binary Inspirals, Phys. Rev. Lett. 118, 051101 (2017a), arXiv:1606.03117 [gr-qc] .
  * Chatziioannou  _et al._ [2017b] K. Chatziioannou, A. Klein, N. Yunes, and N. Cornish, Constructing gravitational waves from generic spin-precessing compact binary inspirals, Phys. Rev. D 95, 104004 (2017b), arXiv:1703.03967 [gr-qc] .
  * Gerosa  _et al._ [2017] D. Gerosa, U. Sperhake, and J. Vošmera, On the equal-mass limit of precessing black-hole binaries, Class. Quant. Grav. 34, 064004 (2017), arXiv:1612.05263 [gr-qc] .
  * Racine [2008] É. Racine, Analysis of spin precession in binary black hole systems including quadrupole-monopole interaction, Phys. Rev. D 78, 044021 (2008), arXiv:0803.1820 [gr-qc] .
  * Schmidt  _et al._ [2011] P. Schmidt, M. Hannam, S. Husa, and P. Ajith, Tracking the precession of compact binaries from their gravitational-wave signal, Phys. Rev. D 84, 024046 (2011), arXiv:1012.2879 [gr-qc] .
  * Boyle  _et al._ [2011] M. Boyle, R. Owen, and H. P. Pfeiffer, Geometric approach to the precession of compact binaries, Phys. Rev. D 84, 124011 (2011), arXiv:1110.2965 [gr-qc] .
  * Schmidt  _et al._ [2015] P. Schmidt, F. Ohme, and M. Hannam, Towards models of gravitational waveforms from generic binaries: II. Modelling precession effects with a single effective precession parameter, Phys. Rev. D 91, 024043 (2015), arXiv:1408.1810 [gr-qc] .
  * Gerosa  _et al._ [2021] D. Gerosa, M. Mould, D. Gangardt, P. Schmidt, G. Pratten, and L. M. Thomas, A generalized precession parameter $\chi$p to interpret gravitational-wave data, Phys. Rev. D 103, 064067 (2021), arXiv:2011.11948 [gr-qc] .
  * Yunes  _et al._ [2009] N. Yunes, K. G. Arun, E. Berti, and C. M. Will, Post-circular expansion of eccentric binary inspirals: Fourier-domain waveforms in the stationary phase approximation, Phys. Rev. D 80, 084001 (2009), arXiv:0906.0313 [gr-qc] .
  * Mishra  _et al._ [2015] C. K. Mishra, K. G. Arun, and B. R. Iyer, Third post-Newtonian gravitational waveforms for compact binary systems in general orbits: Instantaneous terms, Phys. Rev. D 91, 084040 (2015), arXiv:1501.07096 [gr-qc] .
  * Boetzel  _et al._ [2019] Y. Boetzel, C. K. Mishra, G. Faye, A. Gopakumar, and B. R. Iyer, Gravitational-wave amplitudes for compact binaries in eccentric orbits at the third post-Newtonian order: Tail contributions and postadiabatic corrections, Phys. Rev. D 100, 044018 (2019), arXiv:1904.11814 [gr-qc] .
  * Ebersold  _et al._ [2019] M. Ebersold, Y. Boetzel, G. Faye, C. K. Mishra, B. R. Iyer, and P. Jetzer, Gravitational-wave amplitudes for compact binaries in eccentric orbits at the third post-Newtonian order: Memory contributions, Phys. Rev. D 100, 084043 (2019), arXiv:1906.06263 [gr-qc] .
  * Marsat and Baker [2018] S. Marsat and J. G. Baker, Fourier-domain modulations and delays of gravitational-wave signals,  (2018), arXiv:1806.10734 [gr-qc] .
  * Rubbo  _et al._ [2004] L. J. Rubbo, N. J. Cornish, and O. Poujade, Forward modeling of space-borne gravitational wave detectors, Phys. Rev. D 69, 082003 (2004), arXiv:gr-qc/0311069 [gr-qc] .
  * Tinto and Dhurandhar [2021] M. Tinto and S. V. Dhurandhar, Time-delay interferometry, Living Rev. Rel. 24, 1 (2021), arXiv:gr-qc/0409034 [gr-qc] .
  * Damour  _et al._ [2001] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Comparison of search templates for gravitational waves from binary inspiral, Phys. Rev. D 63, 044023 (2001), arXiv:gr-qc/0010009 [gr-qc] .
  * Damour  _et al._ [2005] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Erratum: Comparison of search templates for gravitational waves from binary inspiral: 3.5PN update [Phys. Rev. D 66, 027502 (2002)], Phys. Rev. D 72, 029901 (2005).
  * Kawamura  _et al._ [2020] S. Kawamura _et al._ , Current status of space gravitational wave antenna DECIGO and B-DECIGO, arXiv e-prints , arXiv:2006.13545 (2020), arXiv:2006.13545 [gr-qc] .
  * Prince  _et al._ [2002] T. A. Prince, M. Tinto, S. L. Larson, and J. W. Armstrong, LISA optimal sensitivity, Phys. Rev. D 66, 122002 (2002), arXiv:gr-qc/0209039 [gr-qc] .
  * LISA Science Study Team [2018] LISA Science Study Team,  _LISA Science Requirements Document_ , Tech. Rep. ESA-L3-EST-SCI-RS-001 (ESA, 2018) www.cosmos.esa.int/web/lisa/lisa-documents/.
  * Karnesis  _et al._ [2021] N. Karnesis, S. Babak, M. Pieroni, N. Cornish, and T. Littenberg, Characterization of the stochastic signal originating from compact binaries populations as measured by LISA, arXiv e-prints , arXiv:2103.14598 (2021), arXiv:2103.14598 [astro-ph.IM] .
  * Poisson [1998] E. Poisson, Gravitational waves from inspiraling compact binaries: The quadrupole-moment term, Phys. Rev. D 57, 5287 (1998), arXiv:gr-qc/9709032 [gr-qc] .
  * Bohé  _et al._ [2013] A. Bohé, S. Marsat, G. Faye, and L. Blanchet, Next-to-next-to-leading order spin-orbit effects in the near-zone metric and precession equations of compact binaries, Class. Quant. Grav. 30, 075017 (2013), arXiv:1212.5520 [gr-qc] .
  * Bohé  _et al._ [2015] A. Bohé, G. Faye, S. Marsat, and E. K. Porter, Quadratic-in-spin effects in the orbital dynamics and gravitational-wave energy flux of compact binaries at the 3PN order, Class. Quant. Grav. 32, 195010 (2015), arXiv:1501.01529 [gr-qc] .
  * Marsat [2015] S. Marsat, Cubic-order spin effects in the dynamics and gravitational wave energy flux of compact object binaries, Classical and Quantum Gravity 32, 085008 (2015), arXiv:1411.4118 [gr-qc] .