# Eccentricity evolution of spinning binaries

Sayak Datta<EMAIL_ADDRESS>Max-Planck-Institut für Gravitationsphysik
(Albert-Einstein-Institut), D-30167 Hannover, Germany  Leibniz Universität
Hannover, D-30167 Hannover, Germany

###### Abstract

We study the evolution of the eccentricity of an eccentric orbit with spinning
components. We develop a prescription to express the evolving eccentricity in
terms of initial eccentricity and frequency. For that purpose we considered
the spins to be perpendicular to the orbital plane. Using this we found an
analytical result for the contribution of spin in eccentricity evolution. As a
result, we expressed orbital eccentricity in a series of initial eccentricity
and gravitational wave frequency. The prescription developed here can easily
be used to find arbitrarily higher-order contributions of initial
eccentricity. With the eccentricity evolution at hand, we computed the
evolving energy and angular momentum fluxes for eccentric orbit with spinning
components. This result can be used to construct the waveforms of spinning
compact objects in an eccentric orbit.

##  I Introduction

In recent times, the detection of gravitational waves (GWs) from the
coalescence of compact binaries with ground-based detectors Aasi et al.
(2015); Acernese et al. (2015) has opened up a new era of astronomy Abbott et
al. (2019a, 2021a). Most of the sources are believed to be black hole (BH)
binaries. The merger of two neutron stars (NSs) was also observed in the event
GW170817 Abbott et al. (2017), and possibly also GW190425 Abbott et al.
(2020a). More recently, detections of GW200105 and GW200115 Abbott et al.
(2021b) were also made where it is believed that it is BH-NS binary.
Currently, the existing detectors are continuously being upgraded. Alongside,
there are proposals for several next-generation ground-based detectors such as
the Einstein telescope Maggiore et al. (2020) and cosmic explorer Reitze et
al. (2019). These detectors will be significantly more sensitive compared to
the current detectors. As a result, it will be possible to measure very small
features in the signals. Similarly, space-based detectors such as Laser
Interferometer Space Antenna (LISA) Amaro-Seoane et al. (2017) are also being
built. LISA will observe binaries comprising supermassive bodies. These
sources will be either very loud or will last very long in the detector for
the detector to measure very small features in the signal. Therefore modeling
the signals as accurately as possible has become a necessity.

In the context of GW astronomy, primarily the focus has been on the circular
orbits. This is reasonable as we expect the stellar mass binaries to have low
eccentricities Peters (1964); Peters and Mathews (1963). However, compact
binaries that formed via the dynamical interactions in dense stellar
environments or through the Kozai-Lidov processes Kozai (1962); Lidov (1962),
are expected to retain residual eccentricities $e_{0}\gtrsim.1$ during
observation with ground-based detectors Abbott et al. (2019b). Therefore, the
measurement of nonzero eccentricities of the binaries may shed light on our
understanding of the formation channels. On the other hand extreme mass ratio
inspirals are also expected to have large eccentricities in the observable
band Amaro-Seoane (2018). Hence, it is crucial to understand and model the
effects of eccentricity on a binary.

Interestingly, some of the recent analyses Abbott et al. (2020b); Kimball et
al. (2021); Romero-Shaw et al. (2021, 2022); O’Shea and Kumar (2021) support
the presence of eccentricity in the observed binary black hole (BBH) events.
Alongside, it was argued in Ref. Favata et al. (2022), that the presence of
even smaller eccentricities $e_{0}\sim.01-.05$ may induce systematic biases in
parameter estimation analyses. Since the next generation ground-based
detectors will have improved low-frequency sensitivity compared to the current
generation detectors, these detectors can make confident observations of
eccentric systems Lower et al. (2018).

There have been efforts in the past to model the inspiral waveforms from
compact binary mergers Mishra et al. (2015); Moore et al. (2016); Tanay et al.
(2016); Boetzel et al. (2019); Ebersold et al. (2019); Konigsdorffer and
Gopakumar (2006); Moore and Yunes (2019); Paul and Mishra (2022); Datta
(2023). Although the effect of spins has been modeled within the post-
Newtonian formalism Marsat et al. (2013); Hartung and Steinhoff (2011); Levi
and Steinhoff (2015); Bohe et al. (2013); Bohé et al. (2013); Marsat et al.
(2014); Bohé et al. (2015); Blanchet et al. (2011); Chatziioannou et al.
(2013); Henry et al. (2022); Arun et al. (2009a), a combined treatment
including spins and eccentricity is largely absent. There have been efforts
such as those of Ref. Kidder (1995); Buonanno et al. (2013); Majar and Vasuth
(2008); Klein and Jetzer (2010); Klein et al. (2018); Klein (2021); Blanchet
et al. (2008); Kidder (2008); Blanchet (2014); Phukon et al. (2019) that
attempt to address this concern to some extent. However, we still lack
waveforms as per the requirement.

Although there are some works on eccentric waveforms, mostly they do not
consider the effect of the interaction between spin and eccentricity. Keeping
this in mind, we will try to find the eccentricity evolution and fluxes due to
such interaction. In Ref. Klein and Jetzer (2010) the equations governing the
time evolution of the orbital parameters, and in particular the eccentricity,
including the spin-orbit and the spin-spin couplings were computed. These are
needed to accurately compute the post-Newtonian approximation up to 2PN
accuracy. However, these results are not sufficient to compute GW waveforms.
The computed results depend on the “instantaneous” eccentricity of the orbit.
However, in the inspiral time scale, the eccentricity of an orbit evolves.
Hence, it is important to express the instantaneous eccentricity in terms of
some “initial eccentricity” and GW frequency. In the leading order of initial
eccentricity for nonspinning orbit this was explored in a few works Arun et
al. (2009b); Moore et al. (2016); Moore and Yunes (2019). Some, higher order
terms of initial eccentricity were computed in ref. Yunes et al. (2009),
albeit for non-spinning binaries. However, in ref. Datta (2023) a more
comprehensive prescription was developed for non-spinning binaries, that can
be used to iteratively find the arbitrarily higher order terms in initial
eccentricities and GW frequency.

In the current work we extend the results of Ref. Datta (2023) to spinning
binaries. Therefore, it will be for the first time that the eccentricity
evolution of spinning binaries will be analytically expressed in terms of
initial eccentricity and GW frequency. This as a result provides eccentricity-
spin coupling terms in fluxes. We develop the prescription for spins
perpendicular to the orbital plane. The prescription discussed here can be
extended to arbitrarily high order in eccentricity, in principle. With the
newly found expressions, we also compute the leading order shifts in the
fluxes due to the eccentricity-spin coupling terms. Therefore these results
can be directly used to model the GW waveforms for spinning bodies in an
eccentric orbit.

##  II Equations of eccentric orbits

Our purpose in the current work is to formulate a method that can be used to
analytically compute the eccentricity of an orbit with spinning components
with respect to the emitted GW frequency. As a result, this can be used to
compute frequency-dependent fluxes and consequently the GW waveform. For this
purpose we will primarily follow the notations and discussions in Ref. Klein
and Jetzer (2010). The spin-orbit couplings appear at 1.5PN order and spin-
spin couplings appear at 2PN order. For this reason Ref. Klein and Jetzer
(2010) considered only the Newtonian and spin-coupling terms in the equations
of motion. For simplicity, we will use a system of units where $G=c=M=1$,
where $M$ is the total mass of the system. The boldface represents a
$3-$vector and a hat above represents a unit vector. The Lagrangian of the
system in the center of mass frame used in Ref. Kidder et al. (1993); Kidder
(1995):

$\begin{split}\mathcal{L}=&\frac{\nu}{2}v^{2}+\frac{\nu}{r}+\frac{\nu}{2}({\bf
v}\times{\bf a}).{\bf\xi}-\frac{2\nu}{r^{3}}{\bf(x\times
v).({\bf\zeta}+{\bf\xi})}\\\ &+\frac{1}{r^{3}}{\bf
S_{1}.S_{2}}-\frac{3}{r^{5}}{\bf(x.S_{1})(x.S_{2})},\end{split}$ (1)

where $\nu=m_{1}m_{2}$ and

$\displaystyle{\bf\zeta}=$ $\displaystyle{\bf S_{1}}+{\bf S_{2}}$ (2)
$\displaystyle{\bf\xi}=$ $\displaystyle\frac{m_{2}}{m_{1}}{\bf
S_{1}}+\frac{m_{1}}{m_{2}}{\bf S_{2}}.$ (3)

A quasi-Keplerian solution to the equation of motion derived from the
Lagrangian in Eq. (1) was found in Ref. Klein and Jetzer (2010)

$\displaystyle r=$ $\displaystyle a(1-e_{r}\cos{u})+f_{r}\cos{2(v-\psi)}$ (4)
$\displaystyle\phi=$
$\displaystyle(1+k)v+f_{\phi,1}\sin{(v-2\psi)}+f_{\phi,2}\sin{2(v-\psi)}$ (5)
$\displaystyle v=$ $\displaystyle
2\arctan\Big{(}\sqrt{\frac{1+e_{\phi}}{1-e_{\phi}}}\tan\frac{u}{2}\Big{)}$ (6)
$\displaystyle l=$ $\displaystyle n(t-t_{0})=u-e_{t}\sin{u},$ (7)

where $(r,\phi)$ is a polar coordinate system in the plane of motion, $n$ is
the mean motion, $u$, $v$, and $l$ are the eccentric, true, and mean
anomalies, a is the semi-major axis, $e_{t}$, $e_{r}$, and $e_{\phi}$ are
eccentricities, $k$ is the perihelion precession, and the $f_{i}$ are
constants Klein and Jetzer (2010).

In the current work, we will compute only leading order shifts in energy and
angular momentum fluxes due to spin-eccentricity couplings. Therefore, we only
require the following leading order results of reduced energy $(E)$ and
angular averaged reduced angular momentum $(L)$Klein and Jetzer (2010),

$\displaystyle E=$ $\displaystyle-\frac{x}{2}$ (8) $\displaystyle L=$
$\displaystyle\frac{\sqrt{1-e_{t}^{2}}}{x^{1/2}}\left[1-\frac{x^{3/2}\beta(4,3-e_{t}^{2})}{2(1-e_{t}^{2})^{3/2}}-\frac{x^{2}\gamma_{1}}{2(1-e_{t}^{2})^{2}}\right],$
(9)

where,

$\displaystyle\beta(a,b)=$ $\displaystyle a\beta_{a}+b\beta_{b}$ (10)
$\displaystyle\beta_{a}=$ $\displaystyle{\bf\hat{J}.{\bf\zeta}}$ (11)
$\displaystyle\beta_{b}=$ $\displaystyle{\bf\hat{J}.{\bf\xi}}$ (12)
$\displaystyle\gamma_{1}=$ $\displaystyle\frac{1}{\nu}\left[{\bf
S_{1}.S_{2}}-3({\bf\hat{J}.S_{1}})({\bf\hat{J}.S_{2}})\right]$ (13)
$\displaystyle\gamma_{2}=$ $\displaystyle\frac{1}{\nu}|\hat{\bf J}\times{\bf
S}_{1}||\hat{\bf J}\times{\bf S}_{2}|,$ (14)

where ${\bf S}_{i}$ are the individual spins and ${\bf J}$ is the reduced
orbital angular momentum discussed in Ref. Klein and Jetzer (2010). ${\bf J}$
helps encapsulate the effects of spin precession. It can be checked from Eq.
(10) in Ref. Klein and Jetzer (2010), when the spins of the bodies are either
aligned or antialigned with the direction perpendicular to the orbital plane,
say $\hat{\bf L}$, ${\bf J}$ is also aligned or antialigned with $\hat{\bf
L}$. In such a case therefore, $\hat{\bf J}||\hat{\bf L}$.

Spin and eccentricity dependence of several orbital quantities were found in
Ref. Klein and Jetzer (2010). These can also be expressed in a similar manner
as follows,

$\displaystyle k=$
$\displaystyle-\frac{x^{3/2}\beta(4,3)}{(1-e_{t}^{2})^{3/2}}-\frac{3x^{2}\gamma_{1}}{2(1-e_{t}^{2})^{2}}$
(15) $\displaystyle e_{r}=$ $\displaystyle
e_{t}\left[1-\frac{x^{3/2}\beta(2,1)}{\sqrt{1-e_{t}^{2}}}-\frac{x^{2}\gamma_{1}}{2(1-e_{t}^{2})}\right]$
(16) $\displaystyle e_{\phi}=$ $\displaystyle
e_{t}\left[1-\frac{x^{3/2}\beta(2,2)}{\sqrt{1-e_{t}^{2}}}-\frac{x^{2}\gamma_{1}}{(1-e_{t}^{2})}\right]$
(17) $\displaystyle f_{r}=$ $\displaystyle-\frac{x}{2(1-e_{t}^{2})}\gamma_{2}$
(18) $\displaystyle f_{\phi,1}=$
$\displaystyle-\frac{e_{t}x^{2}}{(1-e_{t}^{2})^{2}}\gamma_{2}$ (19)
$\displaystyle f_{\phi,2}=$
$\displaystyle-\frac{x^{2}}{4(1-e_{t}^{2})^{2}}\gamma_{2}.$ (20)

In the quasi-Keplerian formalism, PN expansions of elliptical orbit quantities
are performed most naturally in terms of the radial orbit angular frequency
$\omega_{r}\equiv n$. This frequency is the mean motion or periastron-to-
periastron angular frequency. However, the limitation of this frequency is
that correspondence with the circular orbit limit is not straightforward.
Unlike the eccentric case, the circular orbit quantities are more naturally
expanded in terms of the azimuthal or $\phi-$angular frequency
$\omega_{\phi}\equiv\xi_{\phi}/M$. For this purpose we use the variable $x$,
defined as $\xi_{\phi}=x^{3/2}$ Moore et al. (2016). The relationship between
these two frequencies is as follows,

$\begin{split}\frac{n}{\xi_{\phi}}=&\,1-\frac{3x}{1-e_{t}^{2}}-\left[18-28\nu+\Big{(}51-26\nu\Big{)}e_{t}^{2}\right]\frac{x^{2}}{4(1-e_{t}^{2})^{2}}\\\
\equiv&\,\mathcal{F}(x,e_{t}(x)).\end{split}$ (21)

##  III eccentricity evolution

In this section, our focus will be on computing the eccentricity of a binary
in terms of binary frequency. Currently, these results are available for non-
spinning binaries Datta (2023); Moore et al. (2016). In this section, we will
construct a prescription that can be used to find the eccentricity in the
power of initial eccentricity and frequency to an arbitrary power for a
spinning binary.

In Sec. II we discussed the parametric solution in the quasi-Keplerian
approach. This is done using the constants of the motions considering only the
conservative part. However, during orbital motion, the system emits GWs that
carry away energy and angular momentum. As a result, the conserved quantities
start to change in the inspiral time scale. This corresponds to the
dissipative part of the equation of motion. Conventionally it is tackled by
solving the evolution equation of the conserved quantities in the inspiral
time scale. For the current work, we only require the evolution equation of
the $e_{t}$ and $n$. Their evolution equation can be expressed as Arun et al.
(2009b); Klein and Jetzer (2010),

$\begin{split}\frac{dn}{dt}=&\frac{\nu
x^{11/2}}{(1-e_{t}^{2})^{7/2}}\Big{(}\mathcal{N}_{n}+x\mathcal{N}_{1}+x^{3/2}\mathcal{N}_{hered}-x^{3/2}\mathcal{N}_{\beta}-x^{2}\mathcal{N_{\sigma}}+x^{2}\mathcal{N}_{\tau}\Big{)}\\\
2e_{t}\frac{de_{t}}{dt}=&\frac{-\nu
x^{4}}{(1-e_{t}^{2})^{5/2}}\Big{(}2e_{t}^{2}\mathcal{E}_{n}+2e_{t}^{2}x\mathcal{E}_{1}+2e_{t}^{2}x^{3/2}\mathcal{E}_{hered}-x^{3/2}\mathcal{E}_{\beta}-x^{2}\mathcal{E_{\sigma}}+x^{2}\mathcal{E}_{\tau}\Big{)}.\end{split}$
(22)

To find an expression of $de_{t}/dx$ it is required to find $dx/dt$. With the
result at hand, we can compute using, $de_{t}/dx=(de_{t}/dt)/(dx/dt)$. To find
$dx/dt$ we used Eq. (22) along with Eq. (21). We find,

$\frac{dn}{dt}=\frac{dx}{dt}\Bigg{(}\frac{3}{2}x^{1/2}\mathcal{F}(x,e_{t})+x^{3/2}\Big{(}\frac{\partial\mathcal{F}(x,e_{t})}{\partial
x}+\frac{\partial\mathcal{F}(x,e_{t})}{\partial
e_{t}}\frac{de_{t}}{dx}\Big{)}\Bigg{)}.$ (23)

Once an expression of $de_{t}/dx$ is found, the right-hand side of the
equation can be expressed in a series expansion of $e_{t}$ where the
coefficients of each term depend on $x$. This equation then is solved to find
the expression of $e_{t}$ in terms of an initial eccentricity $e_{0}$ and $x$.
In this section, we will develop a prescription that can be used to find very
high-order powers in $e_{0}$.

To demonstrate the prescription, we will keep up to power $e_{t}^{3}$. Then
the eccentricity evolution can be expressed as follows,

$\frac{de_{t}}{dx}=e_{t}^{-1}f_{-1}(x)+e_{t}f_{1}(x)+e_{t}^{3}f_{3}(x)+\mathcal{O}(e_{t}^{4}).$
(24)

It will be demonstrated later that for aligned or antialigned spins
$f_{-1}(x)=0$. Hence, we will fix $f_{-1}(x)=0$ from here on. In such case,
the equation simplifies to,

$\frac{de_{t}}{dx}=e_{t}f_{1}(x)+e_{t}^{3}f_{3}(x)+\mathcal{O}(e_{t}^{4}).$
(25)

Note, we have not explicitly specified the functional expression of $f_{1}(x)$
and $f_{3}(x)$. Once the prescription is described, we will substitute them
with the required PN accuracy.

From Eq. (25) the resulting solution can be found by integrating the equation
on both sides. During the process, we identify that $e_{t}\rightarrow e_{0}$,
i.e. the initial eccentricity, when $x\rightarrow x_{0}$, i.e. the reference
frequency. we find,

$\begin{split}\int_{e_{0}}^{e_{t}}\frac{de_{t}}{e_{t}}=&\int_{x_{0}}^{x}dx\Big{(}f_{1}(x)+e_{t}^{2}f_{3}(x)+\mathcal{O}(e_{t}^{3})\Big{)}\\\
\ln(\frac{e_{t}}{e_{0}})=&\int_{x_{0}}^{x}dxf_{1}(x)+\int_{x_{0}}^{x}dxe_{t}^{2}f_{3}(x)+\mathcal{O}(e_{t}^{3})\\\
e_{t}=&e_{0}e^{\int_{x_{0}}^{x}dxf_{1}(x)}e^{\int_{x_{0}}^{x}dxe_{t}^{2}f_{3}(x)+\mathcal{O}(e_{t}^{3})}.\end{split}$
(26)

Note, that the first integral on the right-hand side is independent of the
eccentricity while the second one depends on it. Therefore although the second
integral can not be computed without the explicit knowledge of $e_{t}$, the
first integral can be computed. Hence the expression can be rearranged as,

$e_{t}=e_{0}\frac{e^{F_{1}{(x})}}{e^{F_{1}{(x_{0})}}}e^{\int_{x_{0}}^{x}dxe_{t}^{2}f_{3}(x)+\mathcal{O}(e_{t}^{3})}.$
(27)

The second integral has $e_{t}$ inside the integral. Therefore, without the
exact knowledge of $e_{t}$ in terms of $x$, this integral can not be computed
exactly. However, the eccentricity $e_{t}$ inside the second integral can be
replaced with the above equation and as a result, a leading order term of the
integral can be computed. Therefore, although an exact integral is not
computable, an approximate result can be found which is exact to a particular
order. This as a result can be used to find the next order term. This can be
continued for arbitrary powers of $e_{0}$. In this work, we will only compute
up to $e_{0}^{3}$ term as Eq. (24) keeps only up to $e_{t}^{3}$. But in
principle, this can be continued iteratively by considering higher order terms
in Eq. (24). After rearranging the expressions they can be expressed as,

$\begin{split}\frac{e_{t}e^{F_{1}{(x_{0})}}}{e_{0}e^{F_{1}{(x})}}=&e^{e_{0}^{2}\int_{x_{0}}^{x}dx\frac{e^{2F_{1}{(x})}}{e^{2F_{1}{(x_{0})}}}\Big{(}e^{2\int_{x_{0}}^{x}d\bar{x}e_{t}^{2}f_{3}(\bar{x})+\mathcal{O}(e_{t}^{3})}\Big{)}\,\,f_{3}(x)+\mathcal{O}(e_{t}^{3})}\\\
\frac{e_{t}e^{F_{1}{(x_{0})}}}{e_{0}e^{F_{1}{(x})}}=&\Big{(}1+e_{0}^{2}\int_{x_{0}}^{x}dx\frac{e^{2F_{1}{(x})}}{e^{2F_{1}{(x_{0})}}}\Big{(}e^{2\int_{x_{0}}^{x}d\bar{x}e_{t}^{2}f_{3}(\bar{x})+\mathcal{O}(e_{t}^{3})}\Big{)}\,\,f_{3}(x)+\mathcal{O}(e_{0}^{3})\Big{)}\\\
\frac{e_{t}e^{F_{1}{(x_{0})}}}{e_{0}e^{F_{1}{(x})}}=&\Big{(}1+e_{0}^{2}\int_{x_{0}}^{x}dx\frac{e^{2F_{1}{(x})}}{e^{2F_{1}{(x_{0})}}}\Big{(}1+\mathcal{O}(e_{0}^{2})\Big{)}\,\,f_{3}(x)+\mathcal{O}(e_{0}^{3})\Big{)}.\end{split}$
(28)

This as a result boils down to a series of $e_{0}$, where the coefficients of
the expansions are integrals in $x$. This approach can be used consecutively
after deriving individual coefficients of a particular order. Interestingly,
to find the schematic structure of the coefficients in integral form of an
arbitrary power of $e_{0}$ it is not required to know $f_{i}(x)$. The
expression can be derived in terms of the integrals of $f_{i}(x)$s.

The above results can be expressed in a further simplified and compact form
as,

$\begin{split}e_{t}=&\frac{e^{F_{1}{(x})}}{e^{F_{1}{(x_{0})}}}(e_{0}+e_{0}^{3}\int_{x_{0}}^{x}dx\frac{e^{2F_{1}{(x})}}{e^{2F_{1}{(x_{0})}}}f_{3}(x)+\mathcal{O}(e_{0}^{4}))\\\
e_{t}=&e_{0}\frac{\mathcal{A}(x)}{\mathcal{A}(x_{0})}+e_{0}^{3}\Big{(}\frac{\mathcal{A}(x)\mathcal{B}(x)}{\mathcal{A}(x_{0})\mathcal{A}_{2}(x_{0})}-\frac{\mathcal{A}(x)\mathcal{B}(x_{0})}{\mathcal{A}(x_{0})\mathcal{A}_{2}(x_{0})}\Big{)}\\\
\equiv&\epsilon_{1}e_{0}+\epsilon_{3}e_{0}^{3},\end{split}$ (29)

where,

$\displaystyle\mathcal{A}(x)=$ $\displaystyle e^{F_{1}{(x})}$
$\displaystyle\mathcal{A}_{2}(x)=$ $\displaystyle e^{2F_{1}{(x})}$ (30)
$\displaystyle\mathcal{B}(x)-\mathcal{B}(x_{0})=$
$\displaystyle\int_{x_{0}}^{x}dxe^{2F_{1}{(x})}f_{3}(x).$

These expressions provide us with the eccentricity evolution with respect to
the frequency. These can be used for further computations for fluxes and
waveforms. In the later sections, we will derive equations for the $f_{i}$s.
Using the derived expression finally, we will find the explicit expressions
for the evolving eccentricity. Once the eccentricity in terms of $x$ is known
we will use it to compute the shifts in the fluxes.

##  IV aligned and antialigned orbits

In the last section, we constructed a prescription that can be used to compute
eccentricity evolution in terms of initial eccentricity and $x$. Using this we
have found the general expression for $e_{t}$ up to $e_{0}^{3}$. Although, we
limited ourselves to order $e_{0}^{3}$, in principle this can be easily
extended further.

In this work, we want to compute the leading order shifts due to spin-
eccentricity coupling. Hence we will use 2PN order expressions of $f_{1}$ and
$f_{3}$. However, in 2PN order we will keep only the spin-dependent terms. The
required expressions are as follows Arun et al. (2009b); Klein and Jetzer
(2010):

$\displaystyle\mathcal{N}_{n}=$
$\displaystyle\frac{1}{5}(96+292e_{t}^{2}+37e_{t}^{4})$ (31)
$\displaystyle\mathcal{N}_{1}=$
$\displaystyle\frac{1}{(1-e_{t}^{2})}(-\frac{4846}{35}-\frac{264\nu}{5}+e_{t}^{2}(\frac{5001}{35}-570\nu)+e_{t}^{4}(\frac{2489}{4}-\frac{5061\nu}{10})+e_{t}^{6}(\frac{11717}{280}-\frac{148\nu}{5})$
(32) $\displaystyle\mathcal{N}_{hered}=$ $\displaystyle
4\pi\frac{96(1-e_{t}^{2})^{7/2}}{5}(1+\frac{2335}{192}e_{t}^{2}+\frac{42955}{768}e_{t}^{4})$
(33) $\displaystyle\mathcal{N}_{\beta}=$
$\displaystyle\frac{1}{10(1-e_{t}^{2})^{3/2}}\beta(3088+15528e_{t}^{2}+7026e_{t}^{4}+195e_{t}^{6},2160+11720e_{t}^{2}+5982e_{t}^{4}+207e_{t}^{6})$
(34) $\displaystyle\mathcal{N}_{\sigma}=$
$\displaystyle\frac{1}{160(1-e_{t}^{2})^{2}}\sigma(21952+128544e_{t}^{2}+73752e_{t}^{4}+3084e_{t}^{6},64576+373472e_{t}^{2}+210216e_{t}^{4}$
(35)
$\displaystyle+8532e_{t}^{6},131344e_{t}^{2}+127888e_{t}^{4}+7593e_{t}^{6})$
$\displaystyle\mathcal{N}_{\tau}=$
$\displaystyle\frac{1}{320(1-e_{t}^{2})^{2}}\tau(448+4256e_{t}^{2}+3864e_{t}^{4}+252e_{t}^{6},64+608e_{t}^{2}+552e_{t}^{4}+36e_{t}^{6},16e_{t}^{2}+80e_{t}^{4}+9e_{t}^{6}).$
(36)

$\displaystyle\mathcal{E}_{n}=$
$\displaystyle(\frac{304}{15}+\frac{121}{15}e_{t}^{2})$ (37)
$\displaystyle\mathcal{E}_{1}=$
$\displaystyle\frac{1}{(1-e_{t}^{2})}(-\frac{939}{35}-\frac{4084\nu}{45}+e_{t}^{2}(\frac{29917}{105}-\frac{7753}{30}\nu)+e_{t}^{4}(\frac{13929}{280}-\frac{1664\nu}{45})$
(38) $\displaystyle\mathcal{E}_{hered}=$
$\displaystyle\frac{32}{5}\frac{985(1-e_{t}^{2})^{5/2}}{48}\pi(1+\frac{21729}{3940}e_{t}^{2}-\frac{7007}{788}e_{t}^{4})$
(39) $\displaystyle\mathcal{E}_{\beta}=$
$\displaystyle\frac{e_{t}^{2}}{15(1-e_{t}^{2})^{3/2}}\beta(13048+12000e_{t}^{2}+789e_{t}^{4},9208+10026e_{t}^{2}+835e_{t}^{4})$
(40) $\displaystyle\mathcal{E}_{\sigma}=$
$\displaystyle\frac{1}{240(1-e_{t}^{2})^{2}}\sigma(-320+101664e_{t}^{2}+116568e_{t}^{4}+9420e_{t}^{6},-320+296672e_{t}^{2}+333624e_{t}^{4}$
(41)
$\displaystyle+26820e_{t}^{6},88432e_{t}^{2}+161872e_{t}^{4}+16521e_{t}^{6})$
$\displaystyle\mathcal{E}_{\tau}=$
$\displaystyle\frac{1}{480(1-e_{t}^{2})^{2}}\tau(-320+2720e_{t}^{2}+5880e_{t}^{4}+540e_{t}^{6},-320-160e_{t}^{2}+1560e_{t}^{4}+180e_{t}^{6},16e_{t}^{2}+80e_{t}^{4}+9e_{t}^{6}).$
(42)

$\sigma(a,b,c)$ and $\tau(a,b,c)$ are spin dependent functions defined as
follows,

$\displaystyle\sigma(a,b,c)=$ $\displaystyle
a\sigma_{a}+b\sigma_{b}+c\sigma_{c},$ (43) $\displaystyle\tau(a,b,c)=$
$\displaystyle a\tau_{a}+b\tau_{b}+c\tau_{c}.$ (44)

The functions defined above $\sigma_{i}$ and $\tau_{i}$, where $i=a,b,c$,
depend on the spin magnitude of both of the bodies in the binary and their
directions. The expressions of these functions in terms of individual spins
and their components can be expressed as follows,

$\displaystyle\sigma_{a}=$ $\displaystyle\frac{1}{\nu}{\bf S_{1}.S_{2}}$ (45)
$\displaystyle\sigma_{b}=$
$\displaystyle-\frac{1}{\nu}({\bf\hat{J}.S_{1}})({\bf\hat{J}.S_{2}})$ (46)
$\displaystyle\sigma_{c}=$ $\displaystyle\frac{1}{\nu}|{\bf\hat{J}\times
S_{1}}||{\bf\hat{J}\times S_{2}}|\cos 2\psi$ (47) $\displaystyle\tau_{a}=$
$\displaystyle\sum_{i=1}^{2}\frac{1}{m_{i}^{2}}{\bf S_{i}^{2}}$ (48)
$\displaystyle\tau_{b}=$
$\displaystyle-\sum_{i=1}^{2}\frac{1}{m_{i}^{2}}({\bf\hat{J}.S_{i}})^{2}$ (49)
$\displaystyle\tau_{c}=$
$\displaystyle\sum_{i=1}^{2}\frac{1}{m_{i}^{2}}|{\bf\hat{J}\times
S_{i}}|^{2}\cos 2\psi_{i},$ (50)

where, $\psi_{i}$ is the angle subtended of ${\bf S}_{i}$ in the plane of
motion and the periastron line. When the spins of the bodies are aligned or
antialigned with $\hat{\bf L}$, $\hat{\bf J}||\hat{\bf L}$. As a result,
$\sigma_{c}=\tau_{c}=0$, $\sigma_{a}=-\sigma_{b}=S_{1}S_{2}/\nu$,
$\tau_{a}=-\tau_{b}=\sum_{i=1}^{2}\frac{S_{i}^{2}}{m_{i}^{2}}$, and
$\gamma_{2}=0$, where $S_{i}={\bf S}_{i}.\hat{{\bf L}}$. Using the above
expressions it is straightforward to compute $f_{1}(x)$ and $f_{3}(x)$. Then
using Eq. (III) further computation can be done. We find that $f_{-1}$
contributes at 2-PN order for the spinning case. For non-spinning cases, it
vanishes. The expression is as follows,

$\begin{split}f_{-1}(x)=&-\frac{5}{192}x(2\sigma_{a}+2\sigma_{b}-\tau_{a}-\tau_{b}).\end{split}$
(51)

This behaviour of $f_{-1}$ arises because $e_{t}^{-1}$ comes primarily due to
$\mathcal{E}_{\sigma}$ and $\mathcal{E}_{\tau}$. Note in the above expression
$f_{-1}(x)$ vanishes if $\sigma_{a}=-\sigma_{b}$ and $\tau_{a}=-\tau_{b}$.
Therefore, this term does not contribute when the spins of the bodies are
aligned or antialigned with $\hat{\bf L}$. In such a case we can apply the
prescription constructed in the last section, and eccentricity evolution can
be found. We find,

$\begin{split}f_{1}(x)=&-\frac{19}{12x}\left(1+\frac{(2833-5516\nu)x}{3192}+\frac{1}{456}x^{3/2}(-2452\beta_{a}-1776\beta_{b}+1131\pi)+\frac{x^{2}\left(6618528\sigma_{a}+91728\tau_{a}\right)}{1072512}\right)\\\
f_{3}(x)=&\frac{3323}{576x}\left(1+\frac{(472943-653228\nu)x}{159504}+\frac{x^{3/2}(-353345\beta_{a}-270423\beta_{b}+159321\pi)}{39876}+\frac{x^{2}\left(3432250080\sigma_{a}+39238416\tau_{a}\right)}{375153408}\right).\end{split}$
(52)

With this, we can compute the eccentricity evolution as follows,

$\begin{split}x^{19/12}\mathcal{A}(x)=&1+\frac{(5516\nu-2833)x}{2016}+\frac{x^{3/2}}{432}(2452\beta_{a}+1776\beta_{b}-1131\pi)+\frac{x^{2}(-39711168\sigma_{a}-550368\tau_{a})}{8128512}.\end{split}$
(53)

$\begin{split}x^{19/6}\mathcal{A}_{2}(x)=&1+\frac{(5516\nu-2833)x}{1008}+x^{3/2}\Big{(}\frac{613\beta_{a}}{54}+\frac{74\beta_{b}}{9}-\frac{377\pi}{72}\Big{)}+\frac{x^{2}(-19855584\sigma_{a}-275184\tau_{a})}{2032128}.\end{split}$
(54)

$\begin{split}\mathcal{B}(x)=&-\frac{3323}{1824x^{19/6}}\Big{(}1+\frac{19(576485\nu+64718)x}{5443074}-\frac{19x^{3/2}(890535\pi-2(893893\beta_{a}+517017\beta_{b}))}{7177680}\\\
&+\frac{19x^{2}(-1049901048\sigma_{a}-52036236\tau_{a})}{11817332352}\Big{)}.\end{split}$
(55)

$\begin{split}\mathcal{A}(x)\mathcal{B}(x)=&-\frac{3323}{1824x^{19/4}}\Big{(}1+\frac{(137845708\nu-34236165)x}{29029728}+\frac{x^{3/2}(12451319\beta_{a}+8192481\beta_{b}-5951955\pi)}{1196280}\\\
&+\frac{x^{2}(-3(598429191360\sigma_{a}+13780542048\tau_{a}))}{273111681024}\Big{)}.\end{split}$
(56)

##  V Orbital quantities

In the previous section, we formulated a prescription to compute eccentricity
evolution. We expressed eccentricity in a series expansion of $e_{0}$. This
result considers the spins of the components to be (anti)aligned with
$\hat{\bf L}$. As discussed in Sec. II, several conservative quantities of an
eccentric orbit can be expressed in terms of the PN parameter $x$ and
eccentricity $e_{t}$. These quantities also evolve as the orbit evolves in the
inspiral time scale. With the knowledge of $e_{t}$ in a PN expansion, it is
possible to find expressions for these quantities also. Since,
$f_{r},\,f_{\phi,1},\,f_{\phi,2}$ are proportional to $\gamma_{2}$, they all
vanish. By using Eq. (29) in Eq. (15) and simplifying we find the eccentricity
and spin-dependent shifts as follows,

$\displaystyle\delta_{\chi,e_{0}}k=$
$\displaystyle-\frac{3}{2}e_{0}^{2}\Big{\\{}\epsilon_{1}^{2}x^{3/2}\big{(}4\beta_{a}+3\beta_{b})+2\epsilon_{1}^{2}x^{2}\gamma_{1}\Big{\\}}+e_{0}^{4}\Big{\\{}-\frac{3}{8}x^{3/2}(4\beta_{a}+3\beta_{b})\big{(}5\epsilon_{1}^{4}+8\epsilon_{1}\epsilon_{3}\big{)}-\frac{3}{2}\gamma_{1}x^{2}\big{(}3\epsilon_{1}^{4}+4\epsilon_{1}\epsilon_{3}\big{)}\Big{\\}}$
(57) $\displaystyle\delta_{\chi,e_{0}}e_{r}=$
$\displaystyle-\frac{e_{0}\epsilon_{1}}{2}\big{(}-2+4\beta_{a}x^{3/2}+2\beta_{b}x^{3/2}+\gamma_{1}x^{2}\big{)}$
(58)
$\displaystyle-\frac{e_{0}^{3}}{2}\big{(}-2\epsilon_{3}+2\beta_{a}\epsilon_{1}^{3}x^{3/2}+\beta_{b}\epsilon_{1}^{3}x^{3/2}+4\beta_{a}\epsilon_{3}x^{3/2}+2\beta_{b}\epsilon_{3}x^{3/2}+\epsilon_{1}^{3}\gamma_{1}x^{2}+\epsilon_{3}\gamma_{1}x^{2}\big{)}$
$\displaystyle\delta_{\chi,e_{0}}e_{\phi}=$ $\displaystyle-
e_{0}\epsilon_{1}\big{(}-1+2\beta_{a}x^{3/2}+2\beta_{b}x^{3/2}+\gamma_{1}x^{2}\big{)}$
(59) $\displaystyle-
e_{0}^{3}\big{(}-\epsilon_{3}+\beta_{a}\epsilon_{1}^{3}x^{3/2}+\beta_{b}\epsilon_{1}^{3}x^{3/2}+2\beta_{a}\epsilon_{3}x^{3/2}+2\beta_{b}\epsilon_{3}x^{3/2}+\epsilon_{1}^{3}\gamma_{1}x^{2}+\epsilon_{3}\gamma_{1}x^{2}\big{)}.$

##  VI Energy and angular momentum Fluxes

In the context of GW astronomy, the most crucial thing is to model the GW
waveform emitted by a system. This boils down to modeling both the amplitude
and the phase of the waveform as these quantities evolve. These results are
usually computed from the knowledge of the energy flux $(\dot{E})$ and the
angular momentum flux $(\dot{L})$ from the GW sources. These fluxes can be
separated into explicit spin-independent and spin-dependent parts.

$\displaystyle\dot{E}=$
$\displaystyle\nu(\dot{E}_{N}+\dot{E}_{SO}+\dot{E}_{SS})$ (60)
$\displaystyle\dot{L}=$
$\displaystyle\nu(\dot{L}_{N}+\dot{L}_{SO}+\dot{L}_{SS}),$ (61)

where SO and SS represent the spin-orbit and spin-spin coupling respectively.
These expressions depend on orbital energy and angular momentum as follows,

$\displaystyle\dot{E}_{N}=$
$\displaystyle-\frac{(-2E)^{3/2}}{15L^{7}}(96+292A^{2}+37A^{4})$ (62)
$\displaystyle\dot{E}_{SO}=$
$\displaystyle\frac{(-2E)^{3/2}}{10L^{10}}\beta(2704+7320A^{2}+2490A^{4}+65A^{6},1976+5096A^{2}+1569A^{4}+32A^{6})$
(63) $\displaystyle\dot{E}_{SS}=$
$\displaystyle\frac{(-2E)^{3/2}}{960L^{11}}\Big{[}2\sigma(42048+154272A^{2}+75528A^{4}+3084A^{6},124864+450656A^{2}+215544A^{4}+8532A^{6},$
(64) $\displaystyle
131344A^{2}+127888A^{4}+7593A^{6})-\tau(448+4256A^{2}+3864A^{4}+252A^{6},64+608A^{2}+552A^{4}+36A^{6},$
$\displaystyle 16A^{2}+80A^{4}+9A^{6})\Big{]},$

$\displaystyle\dot{L}_{N}=$
$\displaystyle-\frac{4(-2E)^{3/2}}{5L^{4}}(8+7A^{2})$ (65)
$\displaystyle\dot{L}_{SO}=$
$\displaystyle\frac{(-2E)^{3/2}}{15L^{7}}\beta(2264+2784A^{2}+297A^{4},1620+1852A^{2}+193A^{4})$
(66) $\displaystyle\dot{L}_{SS}=$
$\displaystyle\frac{(-2E)^{3/2}}{20L^{8}}\Big{[}2\sigma(552+996A^{2}+132A^{4},1616+2868A^{2}+381A^{4},894A^{2}+186A^{4})$
(67) $\displaystyle-(8+24A^{2}+3A^{4})\tau(2,1,0)\Big{]},$

where $A=\sqrt{1+2EL^{2}}$.

Using Eq. (8) and the expression of $A$ it is possible to find a PN expression
for the fluxes. We use $\delta_{\chi,e_{0}}$ to represent the sift of
$\dot{E}$ and $\dot{L}$ due to both the nonzero spin and eccentricity. We
compute the leading order eccentricity-spin coupling terms arising from spin-
orbit and spin-spin coupling.

$\displaystyle\delta_{\chi,e_{0}}\dot{E}_{N}=$
$\displaystyle\mathcal{E}_{N,\chi,2}e_{0}^{2}+\mathcal{E}_{N,\chi,4}e_{0}^{4}$
(68) $\displaystyle\delta_{\chi,e_{0}}\dot{L}_{N}=$
$\displaystyle\mathcal{L}_{N,\chi,2}e_{0}^{2}+\mathcal{L}_{N,\chi,4}e_{0}^{4}$
(69) $\displaystyle\delta_{\chi,e_{0}}\dot{E}_{SO}=$
$\displaystyle\mathcal{E}_{SO,2}e_{0}^{2}+\mathcal{E}_{SO,4}e_{0}^{4}$ (70)
$\displaystyle\delta_{\chi,e_{0}}\dot{E}_{SS}=$
$\displaystyle\mathcal{E}_{SS,2}e_{0}^{2}+\mathcal{E}_{SS,4}e_{0}^{4}$ (71)
$\displaystyle\delta_{\chi,e_{0}}\dot{L}_{SO}=$
$\displaystyle\mathcal{L}_{SO,2}e_{0}^{2}+\mathcal{L}_{SO,4}e_{0}^{4}$ (72)
$\displaystyle\delta_{\chi,e_{0}}\dot{L}_{SS}=$
$\displaystyle\mathcal{L}_{SS,2}e_{0}^{2}+\mathcal{L}_{SS,4}e_{0}^{4}.$ (73)

In the above expression, we have decomposed the energy and angular momentum
fluxes in a series expansion with respect to initial eccentricity. We have
kept only up to $e_{0}^{4}$. With more computation of eccentricity evolution,
this can be extended further. The coefficients of these expansions are
themselves series expansion in PN parameter $x$. The corresponding expressions
are computed and demonstrated in the following,

$\displaystyle\mathcal{E}_{SO,2}=$
$\displaystyle\frac{4x^{10/3}(2605\beta_{a}+1872\beta_{b})}{5\mathcal{A}_{0}^{2}}$
(74) $\displaystyle\mathcal{L}_{SO,2}=$
$\displaystyle\frac{2x^{11/6}(5354\beta_{a}+3761\beta_{b})}{15\mathcal{A}_{0}^{2}}$
(75) $\displaystyle\mathcal{E}_{SS,2}=$
$\displaystyle-\frac{2x^{23/6}(3916\sigma_{a}+15\tau_{a})}{5\mathcal{A}_{0}^{2}}$
(76) $\displaystyle\mathcal{L}_{SS,2}=$
$\displaystyle-\frac{2x^{7/3}(1532\sigma_{a}+7\tau_{a})}{5\mathcal{A}_{0}^{2}}.$
(77)

$\begin{split}\mathcal{E}_{N,\chi,2}=&-\frac{x^{11/6}}{1905120\mathcal{A}_{0}^{2}}\Big{(}-2352x^{3/2}(-1236868\beta_{a}-883848\beta_{b})+x^{2}(540800064\gamma_{1}+157(-4963896\sigma_{a}-68796\tau_{a}))\Big{)}.\end{split}$
(78)

$\begin{split}\mathcal{L}_{N,\chi,2}=&-\frac{x^{1/3}}{635040\mathcal{A}_{0}^{2}}\Big{(}-2352x^{3/2}(-131996\beta_{a}-92580\beta_{b})+x^{2}(50295168\gamma_{1}-23\left(4963896\sigma_{a}+68796\tau_{a}\right))\Big{)}.\end{split}$
(79)

$\displaystyle\mathcal{E}_{SO,4}=$
$\displaystyle\frac{x^{1/6}}{1140\mathcal{A}_{0}^{4}\mathcal{A}_{2}(x_{0})}\Big{(}-8656415\mathcal{A}_{0}^{2}\beta_{a}-6220656\mathcal{A}_{0}^{2}\beta_{b}+9080100\mathcal{A}_{2}(x_{0})\beta_{a}+6462546\mathcal{A}_{2}(x_{0})\beta_{b}\Big{)}$
(80) $\displaystyle\mathcal{L}_{SO,4}=$
$\displaystyle\frac{1}{6840\mathcal{A}_{0}^{4}\mathcal{A}_{2}(x_{0})x^{4/3}}\Big{(}-17791342\mathcal{A}_{0}^{2}\beta_{a}-12497803\mathcal{A}_{0}^{2}\beta_{b}+12708720\mathcal{A}_{2}(x_{0})\beta_{a}+8861220\mathcal{A}_{2}(x_{0})\beta_{b}\Big{)}$
(81) $\displaystyle\mathcal{E}_{SS,4}=$
$\displaystyle\frac{x^{2/3}}{2280\mathcal{A}_{0}^{4}\mathcal{A}_{2}(x_{0})}\Big{(}13012868\mathcal{A}_{0}^{2}\sigma_{a}+49845\mathcal{A}_{0}^{2}\tau_{a}-15439704\mathcal{A}_{2}(x_{0})\sigma_{a}-71820\mathcal{A}_{2}(x_{0})\tau_{a}\Big{)}$
(82) $\displaystyle\mathcal{L}_{SS,4}=$
$\displaystyle\frac{1}{2280\mathcal{A}_{0}^{4}\mathcal{A}_{2}(x_{0})x^{5/6}}\Big{(}5090836\mathcal{A}_{0}^{2}\sigma_{a}+23261\mathcal{A}_{0}^{2}\tau_{a}-4189956\mathcal{A}_{2}(x_{0})\sigma_{a}-20406\mathcal{A}_{2}(x_{0})\tau_{a}\Big{)}.$
(83)

$\begin{split}x^{4/3}\mathcal{E}_{N,\chi,4}=&\frac{1}{3162194380800\mathcal{A}_{0}^{4}\mathcal{A}_{2}(x_{0})}\Bigg{[}\mathcal{A}_{0}^{2}\Big{[}428064x^{3/2}(46434089278\beta_{a}+32454792462\beta_{b})+40x^{2}\Big{\\{}81767076876576\gamma_{1}\\\
&-157\Big{(}880186480992\sigma_{a}+16828189560\tau_{a}\Big{)}\Big{\\}}\Big{]}+4357080\mathcal{A}_{2}(x_{0})\Big{[}1344x^{3/2}(-3426812\beta_{a}-2390010\beta_{b})\\\
&-x^{2}\Big{\\{}759773952\gamma_{1}-605\Big{(}2836512\sigma_{a}+39312\tau_{a}\Big{)}\Big{\\}}\Big{]}\Bigg{]}.\end{split}$
(84)

$\begin{split}x^{17/6}\mathcal{L}_{N,\chi,4}=&\frac{336x^{3/2}(-364012\beta_{a}-249321\beta_{b})+x^{2}\Big{\\{}19\Big{(}2836512\sigma_{a}+39312\tau_{a}\Big{)}-18724608\gamma_{1}\Big{\\}}}{90720\mathcal{A}_{0}^{4}}\\\
&+\frac{1}{131758099200\mathcal{A}_{0}^{2}\mathcal{A}_{2}(x_{0})}\Bigg{[}53508x^{3/2}(5167489562\beta_{a}+3528306258\beta_{b})+5x^{2}\Big{\\{}7604453368512\gamma_{1}\\\
&-23\Big{(}880186480992\sigma_{a}+16828189560\tau_{a}\Big{)}\Big{\\}}\Bigg{]},\end{split}$
(85)

where $\mathcal{A}_{0}=\mathcal{A}(x=x_{0})$. Here we have computed the shift
in fluxes only due to the leading order spin-orbit and spin-spin interaction.
These fluxes can further be used to compute the shift in the emitted GWs.
Therefore these expressions will be necessary to construct GW waveforms for
spinning bodies in an eccentric orbit.

##  VII discussion

The main result of the paper is the formulation of the prescription that can
be used to compute the eccentricity evolution of a spinning binary to very
high order. To our knowledge, this is for the first time where spin and
$e_{0}$ couplings have been computed in the expression of eccentricity
evolution. It provides us with a frequency-dependent evolution of eccentricity
in terms of initial eccentricity. We also discussed how it can be extended to
higher orders iteratively. We considered the spins to be (anti)aligned to
simplify the calculations. As a result in the current case $f_{-1}(x)=0$. With
the computed PN expansion of eccentricity, we added several orbital quantities
and the fluxes of energy and angular momentum. These fluxes are crucial for
the computation of the GW waveform. Hence, these results can be used to model
GW waveforms considering eccentricity-spin coupling.

In general $f_{-1}(x)\neq 0$. Considering such cases, need more studies which
we will pursue in the future. In such a case different approaches combining
analytical and numerical fitting may be required. In the presence of
nonvanishing $f_{-1}(x)$ integrating the equation analytically becomes
challenging. To address this a semi-analytical approach can be constructed. By
multiplying both side with $2e_{t}$ in Eq. (24) it can be found,

$\begin{split}\frac{de_{t}^{2}}{dx}=&2f_{-1}(x)+2e_{t}^{2}f_{1}(x)+2e_{t}^{4}f_{3}(x)+\mathcal{O}(e_{t}^{5})\\\
\equiv&f_{0}(x)+e_{t}^{2}f_{2}(x)+e_{t}^{4}f_{4}(x)+\mathcal{O}(e_{t}^{5})\end{split}$
(86)

The semi-analytic approach can be constructed by defining
$\tilde{e}_{t}^{2}=\bar{e}_{t}^{2}+\Delta e_{t}^{2}$. Then rather than solving
the exact differential equation satisfied by $e_{t}$, we can solve for two
approximated differential equations such as

$\frac{d\Delta e_{t}^{2}}{dx}=f_{0}(x)$ (87)

$\frac{d\bar{e}_{t}^{2}}{dx}=\sum_{n=2}^{\infty}\bar{e}_{t}^{n}f_{n}(x).$ (88)

Eq. (87) can be solved analytically. Similarly, Eq. (88) can be solved
analytically in an iterative manner as discussed in the current paper. This
will provide us with an analytical expression for $\tilde{e}_{t}$. Note,
$\tilde{e}_{t}\neq e_{t}$ in general. However, as long as
$\frac{\tilde{e}_{t}-e_{t}}{e_{t}}$ is lesser than the observable statistical
error, this systematic error can be ignored. In such a case by comparing a
numerically found solution for $e_{t}$ from the differential equation and
comparing it with $\tilde{e}_{t}$ we can estimate the systematic error and its
significance compared to the statistical error. This needs to be explored in
detail.

## Acknowledgements

We want to thank K. G. Arun, Antoine Klein, Chandra Kant Mishra, and Khun Sang
Phukon for the very useful discussions. We also like to thank Chandra Kant
Mishra, and Khun Sang Phukon for reading an earlier draft and providing us
with valuable comments.

## References

  * Aasi et al. (2015) J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), eprint 1411.4547.
  * Acernese et al. (2015) F. Acernese et al. (VIRGO), Class. Quant. Grav. 32, 024001 (2015), eprint 1408.3978.
  * Abbott et al. (2019a) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X9, 031040 (2019a), eprint 1811.12907.
  * Abbott et al. (2021a) R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), eprint 2010.14527.
  * Abbott et al. (2017) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017), eprint 1710.05832.
  * Abbott et al. (2020a) B. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020a), eprint 2001.01761.
  * Abbott et al. (2021b) R. Abbott et al. (LIGO Scientific, KAGRA, VIRGO), Astrophys. J. Lett. 915, L5 (2021b), eprint 2106.15163.
  * Maggiore et al. (2020) M. Maggiore, C. Van Den Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, M. Branchesi, S. Clesse, S. Foffa, J. García-Bellido, et al., Journal of Cosmology and Astroparticle Physics 2020, 050 (2020).
  * Reitze et al. (2019) D. Reitze et al., Bull. Am. Astron. Soc. 51, 035 (2019), eprint 1907.04833.
  * Amaro-Seoane et al. (2017) P. Amaro-Seoane et al. (LISA) (2017), eprint 1702.00786.
  * Peters (1964) P. C. Peters, Phys. Rev. 136, B1224 (1964).
  * Peters and Mathews (1963) P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963).
  * Kozai (1962) Y. Kozai, Astron. J. 67, 591 (1962).
  * Lidov (1962) M. L. Lidov, Planetary and Space Science 9, 719 (1962).
  * Abbott et al. (2019b) B. P. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. 883, 149 (2019b), eprint 1907.09384.
  * Amaro-Seoane (2018) P. Amaro-Seoane, Living Rev. Rel. 21, 4 (2018), eprint 1205.5240.
  * Abbott et al. (2020b) R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 900, L13 (2020b), eprint 2009.01190.
  * Kimball et al. (2021) C. Kimball et al., Astrophys. J. Lett. 915, L35 (2021), eprint 2011.05332.
  * Romero-Shaw et al. (2021) I. M. Romero-Shaw, P. D. Lasky, and E. Thrane, Astrophys. J. Lett. 921, L31 (2021), eprint 2108.01284.
  * Romero-Shaw et al. (2022) I. M. Romero-Shaw, P. D. Lasky, and E. Thrane, Astrophys. J. 940, 171 (2022), eprint 2206.14695.
  * O’Shea and Kumar (2021) E. O’Shea and P. Kumar (2021), eprint 2107.07981.
  * Favata et al. (2022) M. Favata, C. Kim, K. G. Arun, J. Kim, and H. W. Lee, Phys. Rev. D 105, 023003 (2022), eprint 2108.05861.
  * Lower et al. (2018) M. E. Lower, E. Thrane, P. D. Lasky, and R. Smith, Phys. Rev. D 98, 083028 (2018), eprint 1806.05350.
  * Mishra et al. (2015) C. K. Mishra, K. G. Arun, and B. R. Iyer, Phys. Rev. D 91, 084040 (2015), eprint 1501.07096.
  * Moore et al. (2016) B. Moore, M. Favata, K. G. Arun, and C. K. Mishra, Phys. Rev. D 93, 124061 (2016), eprint 1605.00304.
  * Tanay et al. (2016) S. Tanay, M. Haney, and A. Gopakumar, Phys. Rev. D 93, 064031 (2016), eprint 1602.03081.
  * Boetzel et al. (2019) Y. Boetzel, C. K. Mishra, G. Faye, A. Gopakumar, and B. R. Iyer, Phys. Rev. D 100, 044018 (2019), eprint 1904.11814.
  * Ebersold et al. (2019) M. Ebersold, Y. Boetzel, G. Faye, C. K. Mishra, B. R. Iyer, and P. Jetzer, Phys. Rev. D 100, 084043 (2019), eprint 1906.06263.
  * Konigsdorffer and Gopakumar (2006) C. Konigsdorffer and A. Gopakumar, Phys. Rev. D 73, 124012 (2006), eprint gr-qc/0603056.
  * Moore and Yunes (2019) B. Moore and N. Yunes, Class. Quant. Grav. 36, 185003 (2019), eprint 1903.05203.
  * Paul and Mishra (2022) K. Paul and C. K. Mishra (2022), eprint 2211.04155.
  * Datta (2023) S. Datta (2023), eprint 2305.03771.
  * Marsat et al. (2013) S. Marsat, A. Bohe, G. Faye, and L. Blanchet, Class. Quant. Grav. 30, 055007 (2013), eprint 1210.4143.
  * Hartung and Steinhoff (2011) J. Hartung and J. Steinhoff, Annalen Phys. 523, 783 (2011), eprint 1104.3079.
  * Levi and Steinhoff (2015) M. Levi and J. Steinhoff, JHEP 09, 219 (2015), eprint 1501.04956.
  * Bohe et al. (2013) A. Bohe, S. Marsat, G. Faye, and L. Blanchet, Class. Quant. Grav. 30, 075017 (2013), eprint 1212.5520.
  * Bohé et al. (2013) A. Bohé, S. Marsat, and L. Blanchet, Class. Quant. Grav. 30, 135009 (2013), eprint 1303.7412.
  * Marsat et al. (2014) S. Marsat, A. Bohé, L. Blanchet, and A. Buonanno, Class. Quant. Grav. 31, 025023 (2014), eprint 1307.6793.
  * Bohé et al. (2015) A. Bohé, G. Faye, S. Marsat, and E. K. Porter, Class. Quant. Grav. 32, 195010 (2015), eprint 1501.01529.
  * Blanchet et al. (2011) L. Blanchet, A. Buonanno, and G. Faye, Phys. Rev. D 84, 064041 (2011), eprint 1104.5659.
  * Chatziioannou et al. (2013) K. Chatziioannou, A. Klein, N. Yunes, and N. Cornish, Phys. Rev. D 88, 063011 (2013), eprint 1307.4418.
  * Henry et al. (2022) Q. Henry, S. Marsat, and M. Khalil, Phys. Rev. D 106, 124018 (2022), eprint 2209.00374.
  * Arun et al. (2009a) K. G. Arun, A. Buonanno, G. Faye, and E. Ochsner, Phys. Rev. D79, 104023 (2009a), [Erratum: Phys. Rev.D84,049901(2011)], eprint 0810.5336.
  * Kidder (1995) L. E. Kidder, Phys. Rev. D 52, 821 (1995), eprint gr-qc/9506022.
  * Buonanno et al. (2013) A. Buonanno, G. Faye, and T. Hinderer, Phys. Rev. D 87, 044009 (2013), eprint 1209.6349.
  * Majar and Vasuth (2008) J. Majar and M. Vasuth, Phys. Rev. D 77, 104005 (2008), eprint 0806.2273.
  * Klein and Jetzer (2010) A. Klein and P. Jetzer, Phys. Rev. D 81, 124001 (2010), eprint 1005.2046.
  * Klein et al. (2018) A. Klein, Y. Boetzel, A. Gopakumar, P. Jetzer, and L. de Vittori, Phys. Rev. D 98, 104043 (2018), eprint 1801.08542.
  * Klein (2021) A. Klein (2021), eprint 2106.10291.
  * Blanchet et al. (2008) L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, Class. Quant. Grav. 25, 165003 (2008), [Erratum: Class.Quant.Grav. 29, 239501 (2012)], eprint 0802.1249.
  * Kidder (2008) L. E. Kidder, Phys. Rev. D 77, 044016 (2008), eprint 0710.0614.
  * Blanchet (2014) L. Blanchet, Living Rev. Rel. 17, 2 (2014), eprint 1310.1528.
  * Phukon et al. (2019) K. S. Phukon, A. Gupta, S. Bose, and P. Jain, Phys. Rev. D 100, 124008 (2019), eprint 1904.03985.
  * Arun et al. (2009b) K. G. Arun, L. Blanchet, B. R. Iyer, and S. Sinha, Phys. Rev. D 80, 124018 (2009b), eprint 0908.3854.
  * Yunes et al. (2009) N. Yunes, K. G. Arun, E. Berti, and C. M. Will, Phys. Rev. D 80, 084001 (2009), [Erratum: Phys.Rev.D 89, 109901 (2014)], eprint 0906.0313.
  * Kidder et al. (1993) L. E. Kidder, C. M. Will, and A. G. Wiseman, Phys. Rev. D 47, R4183 (1993), eprint gr-qc/9211025.