# Efficient fully precessing gravitational waveforms for binaries with neutron
stars

Michael LaHaye  Department of Physics, University of Guelph, Guelph, Ontario,
N1G 2W1, Canada  Huan Yang  Perimeter Institute for Theoretical Physics,
Ontario, N2L 2Y5, Canada  Department of Physics, University of Guelph, Guelph,
Ontario, N1G 2W1, Canada  Béatrice Bonga  Institute for Mathematics,
Astrophysics and Particle Physics,Radboud University, 6525 AJ Nijmegen, The
Netherlands  Zhenwei Lyu  Kavli Institute for Astronomy and Astrophysics,
Peking University, Beijing 100871, China  Department of Physics, University of
Guelph, Guelph, Ontario, N1G 2W1, Canada

###### Abstract

We construct an efficient frequency domain waveform for generic circular
compact object binaries that include neutron stars. The orbital precession is
solved on the radiation reaction timescale (and then transformed to the
frequency domain), which is used to map the non-precessional waveform from the
source frame of the binary to the lab frame. The treatment of orbital
precession is different from that for precessional binary black holes, as
$\chi_{\rm eff}$ is no longer conserved due to the spin-induced quadrupole
moments of neutron stars. We show that the new waveform achieves $\leq
10^{-4}$ mismatch compared with waveforms generated by numerically evolved
precession for neutron star-black hole systems for $\geq 90\%$ configurations
with component mass/spin magnitude assumed in the analysis and randomized
initial spin directions. We expect this waveform to be useful to test the
nature of the mass-gap objects similar to the one discovered in GW 190814 by
measuring their spin-induced quadrupole moments, as it is possible that these
mass-gap objects are rapidly spinning. It is also applicable for the tests of
black hole mimickers in precessional binary black hole events, if the black
hole mimicker candidates have nontrivial spin-induced quadrupole moments.

##  I Introduction

Since the first gravitational wave (GW) observation in 2015, almost a hundred
binary coalescences have been detected LIGOScientific:2021djp ;
LIGOScientific:2020ibl ; LIGOScientific:2018mvr , including black hole
binaries (BHBH), neutron star binaries (NSNS) and neutron star-black hole
binaries (NSBH). In order to distinguish neutron stars from black holes in
these binaries, besides information from electromagnetic counterparts, a
natural way is to measure the tidal Love number of the compact object
Flanagan:2007ix , which is zero for black holes but could be
$\mathcal{O}(10^{2}-10^{3})$ for normal neutron stars. However, no definite
detection of nonzero Love number has been made yet. For example, in the first
binary neutron star event detected (GW 170817), the dimensionless tidal Love
number is constrained to be $\leq 800$ LIGOScientific:2017vwq ; Radice:2017lry
. In addition, the tidal Love number decreases sharply with increasing neutron
star mass. For heavy neutron stars near their maximum mass, the corresponding
dimensionless Love number is likely only $\mathcal{O}(1)$ Hinderer:2007mb ;
Hu:2020ujf  , which is at best detectable by the third-generation
gravitational wave detectors Castro:2022mpw . Without definite information
from tidal Love numbers and electromagnetic counterparts, a common practice is
to label compact objects with mass $\leq 3M_{\odot}$ as neutron stars, and
those with mass $\geq 5M_{\odot}$ as black holes. This is motivated by the
lack of black hole observation from X-ray binaries (the “mass-gap”) and
neutron star equation-of-state considerations. However, it remains an open
question whether there is a detectable population of low-mass black holes
(that is, black holes with comparable masses to neutron stars) Yang:2017gfb ,
possibly coming from delayed supernovae explosions, binary neutron star
mergers or primordial black holes. A recent study suggests that they may also
take place in extreme mass ratio inspirals suitable for space-borne
gravitational wave detections Pan:2021lyw , with an accelerated formation rate
through the interaction with accretion disks Pan:2021oob ; Pan:2021ksp . It is
both physically interesting and astrophysically important to unambiguously
identify the nature of (at least some of the) mass-gap objects.

The measurement of the spin-induced quadrupole moment may provide a promising
method to distinguish neutron stars from low-mass black holes in cases in
which this may not be possible otherwise. A key example being GW190814
LIGOScientific:2020zkf , which describes the coalescence of a $\sim
23M_{\odot}$ black hole with a compact object of $\sim 2.6M_{\odot}$; this
compact object may be a low-mass black hole or a heavy neutron star
Vattis:2020iuz ; Godzieba:2020tjn ; Most:2020bba  111If more exotic objects
such as boson stars exist, they may be candidates of mass-gap objects as well
with the corresponding tidal Love number potentially measurable Mendes:2016vdr
; Cardoso:2019rvt .. With the lighter object lying squarely inside the lower
mass gap, an extremely small tidal Love number (expected to be
$\leq\mathcal{O}(1)$) and no optical counterpart, it is difficult to determine
the identity of this object. A measurement of the spin-induced quadrupole
moment from the gravitational waveform could provide insight into the nature
of this object, if sufficiently different from that of a black hole. This
particular event is especially tantalizing: While neutron stars are not in
general expected to have large spins, if this object were a neutron star, it
is natural to expect that it would be rapidly spinning in order to support its
large mass. On the other hand, if it were a low-mass black hole, large spin is
also expected if it is formed in a binary neutron star merger or delayed
supernova with significant accretion. This potentially large spin would make
the effects of spin-precession more pronounced, making it a good candidate for
the measurement of the spin-induced quadrupole moment. In particular, we can
write the spin-induced quadrupole moment as

$\displaystyle Q=\kappa a^{2}m^{3}$ (1)

with $Q$ being the magnitude of the quadrupole moment, $m$ its mass and $a$
its spin. For black holes, the quadrupole constant $\kappa$ equals one (i.e.,
the Kerr metric) and for neutron stars $\kappa$ is a (larger-than-one) number
depending on the star’s mass and equation of state (see Fig. 1 of
Harry:2018hke ).

In order to determine the spin-induced quadrupole moment, accurate waveform
models are required. However, currently, there are no efficient methods of
generating waveforms that include the effects of the spin-induced quadrupole
moment on the precession of generic compact objects. Such methods only exist
for black hole binaries Chatziioannou:2017tdw ; Klein:2021jtd ;
Chatziioannou:2016ezg , but these have not been extended to include the
effects of the quadrupole moment of neutron stars. This is the goal of this
paper. We solve the spin dynamics component of the waveform, other parts of
the waveform generation are then kept the same as non precessional systems.

For precessing compact binaries, it is computationally expensive to track the
binary evolution on the precession timescale, which is generally longer than
the orbital timescale but shorter than the radiation reaction timescale. For
binary black hole systems, it has been shown that there are a sufficient
number of conserved quantities (in particular, $\chi_{\rm eff}$, as discussed
later) such that the spin evolution equation can be solved algebraically after
performing an average over the precession timescale Chatziioannou:2017tdw ;
Klein:2021jtd . However, the quantity $\chi_{\rm eff}$ is no longer conserved
for generic black hole-neutron star systems as neutron stars have $\kappa\neq
1$, so that the precessional binary black hole waveform no longer applies for
binaries with neutron stars, even without considering tidal effects.

The separation of the orbital, precession, and radiation reaction timescales
has been used in much of the literature starting with Apostolatos:1994mx  and
more recently Kesden:2014sla ; Gerosa:2015tea . We will work within the orbit
averaged equations, and take advantage of the separation between the
precession and radiation reaction (RR) timescales to separate our description
into a secular portion (determined by the RR) and a periodic portion
(determined by the conservative dynamics). This process requires first solving
the conservative problem, so we start by finding an approximate analytic
solution to the conservative dynamics. When solving the conservative problem
we have three spins each with three components, so that we have nine variables
describing the system. For black hole binaries there are seven conserved
quantities and two dynamical variables. As mentioned, one of these conserved
quantities, $\chi_{\rm eff}$, is no longer conserved when one of the compact
objects is not a black hole. This leaves us with a choice when solving the
conservative dynamics: either find a new conserved quantity to reduce the
problem back to two dynamical variables or work with three dynamical
variables.

To decide which is the preferable option we looked at the BHBH case for
insight. For black hole binaries, the conservative dynamics have been solved
for several different choices of dynamical variables Chatziioannou:2017tdw ;
Klein:2021jtd . We follow the definitions used in Klein (2021) Klein:2021jtd ,
namely using the sum ($\chi_{\text{eff}}$) and difference ($\delta\chi$) of
the spins projected onto the orbital angular momentum as two of our three
dynamical variables. For this choice, the conservative problem has an exact
solution, but this solution relies on the specific form of the evolution
equation for $\delta\chi$: because it is cubic, its roots can be found exactly
making the solution efficient. We find that the new conserved quantity in the
general case is quadratic in both $\chi_{\text{eff}}$ and $\delta\chi$. This
changes the cubic form of $\delta\chi$’s equation so that we can no longer use
the techniques that applied previously to solve the system. Thus, while using
the conserved quantity directly to reduce the number of dynamical variables
may be tempting at first, we opt to continue working with three dynamical
variables.

With this choice in mind we find an approximate analytic solution to the
conservative dynamics, which can be described in terms of an average,
amplitude, and precession phase. The average (secular) portion and amplitude
both evolve on the RR timescale. The precession phase obviously evolves on the
precession timescale, however, it can be evolved accurately on the RR
timescale because the frequency only evolves on the RR timescale. As a result,
one can evolve all quantities on the RR timescale, while still fully
accounting for the precession.

The organization of the paper is as follows. In Sec. II we lay out the
relevant definitions, and in the following section, Sec. III, we lay out the
relevant equations. In Sec. IV, we solve the conservative dynamics, and in
Sec. V we introduce radiation reaction. In Sec. VI we evaluate the accuracy of
this model, and compare the waveform with the binary black hole case to
illustrate the phase difference introduced by the spin-induced quadrupole
moment.

##  II Angular momenta and spin variables

A note about notation: we will use an arrow $\vec{V}$ to denote vectors, and a
hat $\hat{V}$ to denote unit vectors, and merely the symbol itself $V$ to
denote the magnitudes. From the total mass $M=m_{1}+m_{2}$, we can define the
reduced masses

$\begin{split}\mu_{i}&=\frac{m_{i}}{M}.\end{split}$ (2)

From this, it will be useful to define two combinations of the two reduced
masses, the first being their difference

$\begin{split}\delta\mu=\mu_{1}-\mu_{2},\end{split}$ (3)

which is small in the equal mass limit. The second combination is their
product (the symmetric mass ratio)

$\begin{split}\eta&=\mu_{1}\mu_{2}.\end{split}$ (4)

The magnitude of the orbital angular momentum can be related to the PN
parameter through this symmetric mass ratio

$\begin{split}L&=\frac{\eta}{y}.\end{split}$ (5)

This is one of the conserved quantities of the conservative dynamics. The next
quantities of interest are the dimensionless spin parameter

$\begin{split}\chi_{i}=\frac{S_{i}}{{m_{i}}^{2}}.\end{split}$ (6)

and the reduced spins

$\begin{split}\vec{s}_{i}=\frac{1}{\mu_{i}}\vec{S}_{i}\;.\end{split}$ (7)

The magnitudes of the spins, $S_{1}$ and $S_{2}$ are two more conserved
quantities on the precession timescale. From the reduced spins and the orbital
angular momentum, the total angular momentum is defined as

$\begin{split}\vec{J}=\vec{L}+\mu_{1}\vec{s}_{1}+\mu_{2}\vec{s}_{2}\;.\end{split}$
(8)

Its three components form three additional conserved quantities of the
conservative dynamics (and consequently its magnitude is another, but it is
not unique).

These spins evolve according to the precession equations

$\begin{split}\frac{d\hat{L}}{dt}&=-y^{6}(\Omega_{1}+\Omega_{2})\\\
\frac{d\vec{s}_{1}}{dt}&=\mu_{2}y^{5}\Omega_{1}\\\
\frac{d\vec{s}_{2}}{dt}&=\mu_{1}y^{5}\Omega_{2}\end{split}$ (9)

where

$\begin{split}\Omega_{i}=\left\\{\frac{1}{2}\mu_{i}+\frac{3}{2}\left[1-y\hat{L}\cdot(\kappa_{i}\vec{s}_{i}+\vec{s}_{j})\right]\right\\}\hat{L}\times\vec{s}_{i}+\frac{1}{2}y\vec{s}_{j}\times\vec{s}_{i}\;.\end{split}$
(10)

From these quantities, we denote the sum of the projections of the reduced
spins onto the orbital angular momentum as $\chi_{\rm eff}$

$\begin{split}\chi_{\text{eff}}=\hat{L}\cdot\vec{s}_{1}+\hat{L}\cdot\vec{s}_{2},\end{split}$
(11)

which is one of the conserved quantities of the BBH system first found by
Damour:2001tu  and later shown to be conserved by Racine:2008qv  (after orbit
averaging, it is not instantaneously conserved). The difference is
$\delta\chi$

$\begin{split}\delta\chi=\hat{L}\cdot\vec{s}_{1}-\hat{L}\cdot\vec{s}_{2}.\end{split}$
(12)

The angle between $\hat{L}$ and $\hat{J}$ is denoted

$\begin{split}\cos(\theta_{L})=\hat{L}\cdot\hat{J},\end{split}$ (13)

which can also be written in terms of the other projections using Eq. (8)-(12)

$\begin{split}\cos(\theta_{L})=\frac{1}{2J}(2L+\chi_{\text{eff}}+\delta\mu\delta\chi)\;.\end{split}$
(14)

##  III Equations of motion

The precession equations for compact objects on circular orbits including
leading post-Newtonian (PN) order spin-orbit and spin-spin interactions —
without radiation reaction terms — in principle require nine variables for the
description: three for the spin of each body and three for the Newtonian
angular momentum. At 2.5 post-Newtonian order one can show that there are six
straightforward conserved quantities (the magnitudes of each individual
angular momentum $S_{1}$, $S_{2}$, and $L$ as well as the three components of
the total angular momentum). As a result, we will need three additional
variables to complete the description of the angular momenta.

When the two compact objects are black holes, there is in fact an easily
identifiable seventh conserved quantity, $\chi_{\rm eff}$, which is a
projection of the sum of the reduced spins onto the direction of the Newtonian
orbital angular momentum $L$ (see Appendix. II for the mathematical
expression). This would leave two dynamical quantities to describe the system
completely. Previously for the binary black hole scenario, the square of the
magnitude of the sum of the spins, $S^{2}$, was chosen as the first dynamical
variable, while an angle describing how much the orbital angular momentum has
precessed around the total angular momentum, $\phi_{z}$, was chosen as the
second Chatziioannou:2017tdw ; Kesden:2014sla . Another equivalent choice is
to instead use the cosine of the angle between the orbital angular momentum
and the total angular momentum, $\cos(\theta_{L})$, as the first dynamical
variable (keeping $\phi_{z}$ as the second). Both choices have the
disadvantage that their evolution is singular in the equal mass limit. Klein
proposed a different choice for the first dynamical quantity: $\delta\chi$,
which is the difference of the projection of the reduced spins. With this
choice, the evolution is well defined in the equal mass limit. For this reason
this choice is preferred over using $S^{2}$ or $\cos(\theta_{L})$ and
$\phi_{z}$. The relation between these variables and the angular momenta of
the individual objects and orbital angular momentum is given in App. II (see
also Eqs.(10)-(11) in Klein:2021jtd , where a detailed description of their
evolution appears).

Here, we extend previous results describing the precession of black holes to
other compact objects by including the leading order finite-size effects
through the quadrupole moment constants $\kappa_{i}$. $\kappa_{i}$ is a
coefficient that appears in the quadrupole moment of compact objects,
$Q_{i}=-\kappa_{i}{\chi_{i}}^{2}{m_{i}}^{3}$, and is determined by their
properties, e.g. the equation of state. For rotating black holes $\kappa$ is
equal to one. For other compact objects, $\kappa\neq 1$. Specifically, for
neutron stars, $\kappa>1$ and it is larger for stiffer equations of state
Poisson:1997ha .

For a system where $\kappa_{i}\neq 1$, $\chi_{\rm eff}$ is no longer a
conserved quantity. This leaves us with a choice: either find a seventh
conserved quantity to reduce the number of dynamical variables back to two, or
work with three dynamical variables. While the former option may seem simpler,
as we will show in a later section, because of the more complicated form of
the conserved quantity it will turn out to be easier to work with three
dynamical variables. Therefore, we use $\chi_{\rm eff},\delta\chi$ and
$\phi_{z}$ as our set of dynamical variables. As noted in Klein:2021jtd , the
equations for $\chi_{\text{eff}}$ and $\delta\chi$ take the form:

$\begin{split}\left(\frac{d\delta\chi}{dt}\right)^{2}=\frac{9y^{11}}{4}A_{\delta\chi}^{2}\left(\delta\mu\delta\chi^{3}+B\delta\chi^{2}+C\delta\chi+D\right)\end{split}$
(15)

and

$\begin{split}\left(\frac{d\chi_{\text{eff}}}{dt}\right)^{2}=\frac{9y^{11}}{4}A_{\chi_{\text{eff}}}^{2}\left(\delta\mu\delta\chi^{3}+B\delta\chi^{2}+C\delta\chi+D\right),\end{split}$
(16)

where

$\begin{split}A_{\delta\chi}&=1+y\,A_{\delta\chi,\delta\chi}\delta\chi+y\,A_{\delta\chi,\chi_{\text{eff}}}\chi_{\text{eff}},\\\
A_{\chi_{\text{eff}}}&=y\,A_{\chi_{\text{eff}},\delta\chi}\delta\chi+y\,A_{\chi_{\text{eff}},\chi_{\text{eff}}}\chi_{\text{eff}},\\\
\end{split}$ (17)

with $y$ a post-Newtonian (PN) parameter related to the norm of the orbital
angular momentum (in particular, $y=(M\omega)^{\frac{1}{3}}$ with $M$ the
total mass, and $\omega$ the mean orbital frequency). The coefficients $B,C$,
and $D$ depend on the conserved quantities, the PN parameter, and $\chi_{\rm
eff}$, and are given in App. II (see also (Klein:2021jtd, , App. B)). The
coefficients in $A$ are

$\displaystyle A_{\delta\chi,\delta\chi}$
$\displaystyle=\frac{\kappa_{2}-\kappa_{1}}{4},$ (18a) $\displaystyle
A_{\delta\chi,\chi_{\rm eff}}$
$\displaystyle=-\frac{\kappa_{1}+\kappa_{2}+2}{4},$ (18b) $\displaystyle
A_{\chi_{\rm eff},\delta\chi}$
$\displaystyle=\frac{\kappa_{1}+\kappa_{2}-2}{4},$ (18c) $\displaystyle
A_{\chi_{\rm eff},\chi_{\rm eff}}$
$\displaystyle=\frac{\kappa_{1}-\kappa_{2}}{4}.$ (18d)

When $\kappa_{i}=1$, the equations reduce to those in Klein:2021jtd . To solve
the dynamical equations, we will separate the behaviour of $\chi_{\text{eff}}$
and $\delta\chi$ in the next section.

The amount the orbital angular momentum has precessed around the total angular
momentum $\phi_{z}$ evolves according to

$\begin{split}\frac{d\phi_{z}}{dt}=\frac{1}{\sin^{2}(\theta_{L})}\left[\frac{d\hat{L}}{dt}\cdot\left(\hat{J}\times\hat{L}\right)\right]\;,\end{split}$
(19)

where $\theta_{L}$ is the angle between the orbital angular momentum and the
total angular momentum, and the hats on the orbital angular momentum $L$ and
total angular momentum $J$ indicate that these are the direction vectors in
the directions of the corresponding angular momentum. While $\phi_{z}$ has a
nice physical interpretation, it will also be useful to define a related angle
$\zeta$, used in waveform generation Chatziioannou:2017tdw

$\begin{split}\frac{d\zeta}{dt}=-\cos(\theta_{L})\frac{d\phi_{z}}{dt}\;.\end{split}$
(20)

Once we include radiative effects, we need to take into account that the
conserved quantities are no longer constant. The radiation reaction behavior
is determined, largely, by two things, $dy/dt$ and $dJ/dt$, the latter of
which is given by Klein:2021jtd :

$\begin{split}\frac{dJ}{dt}=-\frac{L}{2Jy}\frac{dy}{dt}\left(2L+\chi_{\text{eff}}+\delta\mu\;\delta\chi\right),\end{split}$
(21)

where $\delta\mu$ is the difference of the reduced masses
$\delta\mu=\mu_{1}-\mu_{2}$.

##  IV Analytic Solution to the Precession Without Radiation Reaction

To solve the equations describing inspiralling precessing black hole or
neutron star binaries on circular orbits, we take advantage of the fact that
an (approximate) analytic solution exists for the conservative problem, as we
will show in this section. In Sec. V, we include the effects of radiation
reaction. To determine the conservative evolution, we separate the problem
into a secular and a periodic part. The secular part of the solution is
constant for $\delta\chi$, $\chi_{\text{eff}}$ and $J$ without considering
radiation reaction. In contrast, the remaining key quantities, that is,
$\phi_{z}$ and $\zeta$, have a secular part that evolves even in the absence
of radiation reaction. As such, for the former quantities, it is satisfactory
to focus only on their periodic parts. The evolution of the angles $\phi_{z}$
and $\zeta$ requires additional treatment. Therefore, we start by examining
the periodic part of $\delta\chi$ and $\chi_{\text{eff}}$, before we discuss
the more involved cases of $\phi_{z}$ and $\zeta$.

As mentioned, in the black hole binary case (for which
$\kappa_{1}=\kappa_{2}=1$), the effective spin $\chi_{\text{eff}}$ is a
constant of the conservative dynamics. This makes the process for solving the
equation for $\delta\chi$ straightforward. In particular, the coefficients in
the equation for $\delta\chi$ in Eq. (15) are all constant so that
$\delta\chi$ oscillates between its minimum and maximum values,
$\delta\chi_{-}$ and $\delta\chi_{+}$, respectively. The solution is then
easily obtained by treating the solution as an average and oscillatory part,
with $\delta\chi_{-}$ and $\delta\chi_{+}$ determining the average part of the
solution as well as the amplitude of the oscillatory part. One is then only
left to evolve the phase of the oscillatory part.

If one or both objects in the binary are neutron stars so that $\kappa_{i}\neq
1$, $\chi_{\text{eff}}$ becomes dynamical and oscillates. As a result, the
coefficients $B$, $C$, and $D$ appearing in Eq. (15) and Eq. (16) are now
dynamical and so are its minima and maxima, $\delta\chi_{-}$ and
$\delta\chi_{+}$. We overcome this complication by deriving an approximate
linear relation between $\delta\chi$ and $\chi_{\text{eff}}$, so that
$\chi_{\text{eff}}$ can be replaced by $\delta\chi$ in Eq. (15) and the
resulting equation for $\delta\chi$ becomes independent of
$\chi_{\text{eff}}$. Consequently, the roots on the right hand side of Eq.
(15) again correspond to the true maximum and minimum of $\delta\chi$. This
allows us to solve for the dynamics of $\delta\chi$. Then, using the linear
relation between $\delta\chi$ and $\chi_{\rm eff}$, the evolution of
$\chi_{\text{eff}}$ is trivially obtained.

In the remainder of this section, we first derive this important linear
relation in two different ways in Sec. IV.1 and IV.2. Next, we solve for the
dynamics of $\delta\chi$ and $\chi_{\rm eff}$ in Sec. IV.3-IV.5. In Sec. IV.6,
we discuss the dynamics for $\phi_{z}$ and $\zeta$.

###  IV.1 Approximate relation between $\delta\chi$ and $\chi_{\text{eff}}$

If the coefficients in Eq. (15) are constant, as is the case when both objects
in the binary are black holes, the _exact_ solution to this equation is given
by the Jacobi elliptic function

$\delta\chi(t)=\delta\chi_{-}+\left(\delta\chi_{+}-\delta\chi_{-}\right)\text{sn}^{2}(\psi^{\prime}(t),m)\;,$
(22)

where $\psi^{\prime}$ describes the phase evolution of the precession and the
parameter
$m=\delta\mu\left(\delta\chi_{+}-\delta\chi_{-}\right)/(\delta\chi_{3}-\delta\mu\;\delta\chi_{-})$
with $\delta\chi_{3}$ containing information of the largest root of the cubic
equation on the right hand side of Eq. (15) (see Eq. (43)). In the limit $m\to
0$, the Jacobi elliptic function reduces to the usual sine function. We will
use this insight to approximate the solutions for $\delta\chi$ and $\chi_{\rm
eff}$ in the generic case for which $\kappa_{i}\neq 1$ as

$\begin{split}\delta\chi&\approx\left<\delta\chi\right>+G_{\delta\chi}\sin(\psi),\\\
\chi_{\text{eff}}&\approx\left<\chi_{\text{eff}}\right>+G_{\chi_{\text{eff}}}\sin(\psi),\\\
\end{split}$ (23)

where

$\begin{split}\left<\delta\chi\right>&=\frac{1}{2}\left(\delta\chi_{+}+\delta\chi_{-}\right),\\\
G_{\delta\chi}&=\frac{1}{2}\left(\delta\chi_{+}-\delta\chi_{-}\right),\\\
\left<\chi_{\text{eff}}\right>&=\frac{1}{2}\left(\chi_{\text{eff},+}+\chi_{\text{eff},-}\right),\\\
G_{\chi_{\text{eff}}}&=\frac{1}{2}\left(\chi_{\text{eff},+}-\chi_{\text{eff},-}\right).\end{split}$
(24)

Here we have used the brackets to indicate that there is some time average
underlying these expressions, indeed the first expression corresponds to the
actual precession average of $\delta\chi$ in Eq. (22) in the $m\to 0$ limit.
This precession average is performed on the precession timescale, $T_{pr}\sim
O(y^{-5})$ which corresponds to the rate at which the precession phase,
$\psi$, evolves. This is in contrast to the radiation reaction timescale,
$T_{rr}\sim O(y^{-8})$, on which the entire evolution occurs. This separation
of scales allows one to disregard the changes in these quantities during the
precession averaging when adding radiation reaction.

In these approximate expressions, we have use $\psi$ instead of
$\psi^{\prime}$ that appears in the original expressions, because we have
simplified the $\sin^{2}$ term using the half angle formula and identified
$2\psi^{\prime}=\psi+\pi/2$ so that the final expression is expressible in
terms of a sine function again. This simple form of the solution matches well
with numerical evolutions of the coupled equations. Moreover, this form is
very powerful as it will allow us to relate $\delta\chi$ and $\chi_{\rm eff}$.
Note that due to the similarity in their derivatives in Eq. (15) and (16),
when one solution has reached its extrema the other must have as well. This is
the reason why both solutions oscillate with a single phase $\psi$. In analogy
to the solutions in Chatziioannou:2017tdw ; Klein:2021jtd , we will refer to
the ansatz in Eq. (23) as the “m=0” approximation.

Operating within this approximation, Eq. (15) and (16) can be related to
obtain

$\begin{split}A_{\chi_{\text{eff}}}^{2}\left(\frac{d\delta\chi}{dt}\right)^{2}=A_{\delta\chi}^{2}\left(\frac{d\chi_{\text{eff}}}{dt}\right)^{2}.\end{split}$
(25)

Substituting the approximate definitions (23) into the above relation then
gives

$(yA_{\chi_{\text{eff}},\delta\chi}\delta\chi+yA_{\chi_{\text{eff}},\chi_{\text{eff}}}\chi_{\text{eff}})^{2}G_{\delta\chi}^{2}=(1+yA_{\delta\chi,\delta\chi}\delta\chi+yA_{\delta\chi,\chi_{\text{eff}}}\chi_{\text{eff}})^{2}G_{\chi_{\text{eff}}}^{2}.$
(26)

Evaluating this at $\psi=0$ relates the amplitude of one quantity,
$G_{\chi_{\text{eff}}}$, to the amplitude of the other,
$G_{\delta\chi}$222Note that the inverse of this expression is not well-
defined in the black hole limit.

$\begin{split}\left|\frac{G_{\chi_{\text{eff}}}}{G_{\delta\chi}}\right|=\left|\frac{yA_{\chi_{\text{eff}},\delta\chi}\left<\delta\chi\right>+yA_{\chi_{\text{eff}},\chi_{\text{eff}}}\left<\chi_{\text{eff}}\right>}{1+yA_{\delta\chi,\delta\chi}\left<\delta\chi\right>+yA_{\delta\chi,\chi_{\text{eff}}}\left<\chi_{\text{eff}}\right>}\right|\;.\end{split}$
(27)

Next, solving for $\sin(\psi)$ in (23), we obtain the desired relation between
$\chi_{\rm eff}$ and $\delta\chi$:

$\begin{split}\chi_{\text{eff}}=\left<\chi_{\text{eff}}\right>-\frac{G_{\chi_{\text{eff}}}}{G_{\delta\chi}}\left<\delta\chi\right>+\frac{G_{\chi_{\text{eff}}}}{G_{\delta\chi}}\delta\chi.\end{split}$
(28)

Having already found the relation between the amplitudes
$G_{\chi_{\text{eff}}}/G_{\delta\chi}$ in terms of the averages, this equation
is entirely determined so long as the average values of $\chi_{\rm eff}$ and
$\delta\chi$ are known. Thus the desired result is achieved: an approximate
linear relation in the absence of radiation reaction between $\delta\chi$ and
$\chi_{\text{eff}}$ defined by their averages (which are known constants).
This relation will be used in Sec. IV.4 to obtain a new dynamical equation for
$\delta\chi$ that no longer depends on the dynamics of $\chi_{\rm eff}$ but
only on its average and amplitude.

###  IV.2 Another perspective: the amplitude relation from a conserved
quantity

The relation between $\chi_{\rm eff}$ and $\delta\chi$ in (28) holds
approximately because it was derived assuming that the solutions for
$\chi_{\rm eff}$ and $\delta\chi$ have the simple form given in Eq. (23).
Here, we show that the amplitude relation in Eq. (27) holds more generically
as long as the average/amplitude of $\chi_{\rm eff}$ and $\delta\chi$ can be
understood purely as a sum/difference between the maxima and minima.

Starting with (15) and (16), the chain rule gives:

$\begin{split}\left(\frac{d\chi_{\text{eff}}}{d\delta\chi}\right)^{2}=\left(\frac{yA_{\chi_{\text{eff}},\delta\chi}\delta\chi+yA_{\chi_{\text{eff}},\chi_{\text{eff}}}\chi_{\text{eff}}}{1+yA_{\delta\chi,\delta\chi}\delta\chi+yA_{\delta\chi,\chi_{\text{eff}}}\chi_{\text{eff}}}\right)^{2}.\end{split}$
(29)

Let us first discuss the special case with
$A_{\chi_{\text{eff}},\chi_{\text{eff}}}=0=A_{\delta\chi,\delta\chi}$ for
which one directly obtains a solution through its quadratures

$1+yA_{\delta\chi,\chi_{\text{eff}}}\chi_{\text{eff}}=\sqrt{1+A_{\chi_{\text{eff}},\delta\chi}A_{\delta\chi,\chi_{\text{eff}}}y^{2}(c_{1}\delta\chi^{2}+c_{2})},$
(30)

where $c_{1}=\pm 1$ arises from a choice of sign when taking the square root
of (29), and $c_{2}$ is a constant of integration. Squaring both sides gives

$2yA_{\delta\chi,\chi_{\text{eff}}}\chi_{\text{eff}}+y^{2}A_{\delta\chi,\chi_{\text{eff}}}^{2}\chi_{\text{eff}}^{2}=A_{\chi_{\text{eff}},\delta\chi}A_{\delta\chi,\chi_{\text{eff}}}y^{2}(c_{1}\delta\chi^{2}+c_{2}).$
(31)

While (23) is an approximate relation for generic values of $\psi$, by
definition of the average and amplitudes in (24), Eq. (23) holds exactly when
$\psi=\pm\pi/2$. As a result, we can substitute Eq. (23) evaluated at
$\psi=\pi/2$ into (31) without loss of generality. Subtracting off the same
equation evaluated at $\psi=-\pi/2$ gives:

$\begin{split}\left|\frac{G_{\chi_{\text{eff}}}}{G_{\delta\chi}}\right|=\left|\frac{yA_{\chi_{\text{eff}},\delta\chi}\left<\delta\chi\right>}{1+yA_{\delta\chi,\chi_{\text{eff}}}\left<\chi_{\text{eff}}\right>}\right|,\end{split}$
(32)

where the factor of $c_{1}$ has been replaced in favor of expressing this as a
relation in terms of the magnitudes of the amplitudes. This equation recovers
the result in Eq. (27) for the case at hand:
$A_{\chi_{\text{eff}},\chi_{\text{eff}}}=0=A_{\delta\chi,\delta\chi}$.

The general case in which $A_{\chi_{\text{eff}},\chi_{\text{eff}}}\neq 0$ and
$A_{\delta\chi,\delta\chi}\neq 0$ is obtained in the same manner.
Specifically, by making the substitution

$\begin{split}Z&=\chi_{\text{eff}}-A_{\chi_{\text{eff}},\delta\chi}/y(A_{\delta\chi,\delta\chi}A_{\chi_{\text{eff}},\chi_{\text{eff}}}-A_{\chi_{\text{eff}},\delta\chi}A_{\delta\chi,\chi_{\text{eff}}})\\\
X&=\delta\chi+A_{\chi_{\text{eff}},\chi_{\text{eff}}}/y(A_{\delta\chi,\delta\chi}A_{\chi_{\text{eff}},\chi_{\text{eff}}}-A_{\chi_{\text{eff}},\delta\chi}A_{\delta\chi,\chi_{\text{eff}}}),\end{split}$
(33)

and using the chain rule, we find

$\begin{split}\frac{dZ}{dX}=\frac{A_{\chi_{\text{eff}},\delta\chi}X+A_{\chi_{\text{eff}},\chi_{\text{eff}}}Z}{A_{\delta\chi,\delta\chi}X+A_{\delta\chi,\chi_{\text{eff}}}Z}.\end{split}$
(34)

The solution to this equation is implicitly given by

$4(A_{\delta\chi,\chi_{\text{eff}}}Z+A_{\delta\chi,\delta\chi}X)^{2}=-(\kappa_{1}\kappa_{2}-1)X^{2}+c_{2},$
(35)

where we have again denoted the constant of integration corresponding to the
conserved quantity as $c_{2}$, to highlight that it arises in a manner similar
as in Eq. (31). Again, evaluating this expression using (23) at $\psi=\pi/2$
and subtracting off the same equation evaluated at $\psi=-\pi/2$, one finds
the amplitude relation in Eq. (27).

As alluded to before, Eq. (35) points to a conserved quantity in the case with
$\kappa_{i}\neq 1$. In the black hole case, this conserved quantity was simply
$\chi_{\text{eff}}$, but here it is nonlinear in $\delta\chi$ and
$\chi_{\text{eff}}$. Of course, in the limit $\kappa_{i}=1$, this equation
simply states that $\chi_{\rm eff}$ is conserved. In Eq. (31) this limit is
slightly more obvious, $A_{\chi_{\text{eff}},\delta\chi}=0$ in this limit and
(after absorbing other conserved quantities into our definition) we find that
$\chi_{\rm eff}$ is conserved.

While the existence of this conserved quantity would reduce the number of
dynamic quantities from three to two, it is more useful to use the approximate
linear relation instead. This is because the linear relation maintains the
cubic nature of Eq. (15) in terms of $\delta\chi$, meaning its roots can be
found analytically. If instead the exact relation in Eq. (31) were used, it
would make the roots of the equation too difficult to find analytically.
Numerical root finding would make this method impractical computationally, for
little benefit in accuracy.

###  IV.3 Finding the initial averages of $\delta\chi$ and $\chi_{\text{eff}}$

In order to calculate the amplitudes of $\delta\chi$ and $\chi_{\rm eff}$,
from Eqs. (27) and (28), it is clear that we need to know their average
values. Consequently, the averages need to be calculated from known quantities
(such as $J$, $y$ and the initial values of $\chi_{\text{eff}}$ and
$\delta\chi$). To do this, first we compute the numerical value of the first
three derivatives for the initial values of $\chi_{\text{eff}}$ and
$\delta\chi$, evaluated using Eq. (15). We then relate these values to the
first three derivatives of the approximate relation Eq. (23), given as:

$\begin{split}\frac{d\delta\chi}{dt}&\approx\frac{d\psi}{dt}G_{\delta\chi}\cos(\psi),\\\
\frac{d^{2}\delta\chi}{dt^{2}}&\approx-\left(\frac{d\psi}{dt}\right)^{2}G_{\delta\chi}\sin(\psi),\\\
\frac{d^{3}\delta\chi}{dt^{3}}&\approx-\left(\frac{d\psi}{dt}\right)^{3}G_{\delta\chi}\cos(\psi),\end{split}$
(36)

where here we have not included the second and third derivatives of $\psi$
because they are zero in the absence of radiation reaction. In the presence of
radiation reaction they are not zero, but as a result of multi-scale analysis
they correspond to the evolution of $\psi$ on the radiation reaction timescale
and as a result their non-inclusion does not produce significant error. We
then solve the ensuing system of equations to obtain:

$\begin{split}\left<\delta\chi\right>&\approx\delta\chi-\delta\chi^{\prime\prime}\delta\chi/\delta\chi^{\prime\prime\prime},\\\
\tan(\psi)&\approx\left|\left(\delta\chi^{\prime\prime}/\delta\chi^{\prime\prime\prime}\right)\sqrt{-\delta\chi^{\prime\prime\prime}/\delta\chi^{\prime}}\right|,\\\
\left|\frac{d\psi}{dt}\right|&\approx\sqrt{-\delta\chi^{\prime\prime\prime}/\delta\chi^{\prime}},\\\
\left|G_{\delta\chi}\right|&\approx\sqrt{(\delta\chi^{\prime})^{2}/(d\psi/dt)^{2}+(\delta\chi^{\prime\prime})^{2}/(d\psi/dt)^{4}},\end{split}$
(37)

where the second equation is true up to a factor of $\pi$ and a plus or minus
sign depending on the signs of the initial derivatives, to get around this we
define the phase in such a way that the amplitude of $\delta\chi$ is always
positive. Since we are working in the conservative dynamics currently, the
average, initial amplitude, and precession frequency are constant, when we
include radiation reaction these will also evolve, the description of which is
given later.

###  IV.4 Finding the amplitudes of $\delta\chi$ and $\chi_{\text{eff}}$

Now, with the averages in hand, the relation (28) can be used to find the
amplitudes. To simplify the next set of calculations we define

$\begin{split}\chi_{\text{eff}}=N_{0}+N_{1}\delta\chi,\end{split}$ (38)

where (as a result of (28)) these coefficients are given by

$\begin{split}N_{0}&=\left<\chi_{\text{eff}}\right>-\frac{G_{\chi_{\text{eff}}}}{G_{\delta\chi}}\left<\delta\chi\right>,\\\
N_{1}&=\frac{G_{\chi_{\text{eff}}}}{G_{\delta\chi}}.\end{split}$ (39)

Substituting this into (15) gives a new equation of the form:

$\begin{split}\left(\frac{d\delta\chi}{dt}\right)^{2}=\frac{9y^{11}}{4}A_{\delta\chi}^{2}\left(X_{3}\delta\chi^{3}+X_{2}\delta\chi^{2}+X_{1}\delta\chi+X_{0}\right),\end{split}$
(40)

where

$\begin{split}X_{0}&=D_{3}N_{0}^{3}+D_{2}N_{0}^{2}+D_{1}N_{0}+D_{0},\\\
X_{1}&=3D_{3}N_{0}^{2}N_{1}+C_{2}N_{0}^{2}+2D_{2}N_{0}N_{1}+C_{1}N_{0}+D_{1}N_{1}+C_{0},\\\
X_{2}&=3D_{3}N_{0}N_{1}^{2}+2C_{2}N_{0}N_{1}+D_{2}N_{1}^{2}+B_{1}N_{0}+C_{1}N_{1}+B_{0},\\\
X_{3}&=D_{3}N_{1}^{3}+C_{2}N_{1}^{2}+B_{1}N_{1}+\delta\mu.\end{split}$ (41)

In the above expression, we have organized the coefficients of $B,C$ and $D$
that appear in Eq. (15) in the following way:

$\begin{split}B&=B_{0}+B_{1}\chi_{\text{eff}},\\\
C&=C_{0}+C_{1}\chi_{\text{eff}}+C_{2}\chi_{\text{eff}}^{2},\\\
D&=D_{0}+D_{1}\chi_{\text{eff}}+D_{2}\chi_{\text{eff}}^{2}+D_{3}\chi_{\text{eff}}^{3}.\end{split}$
(42)

The explicit form of these coefficients is given in App. II. As discussed,
because this cubic no longer implicates $\chi_{\text{eff}}$, the roots are now
the true maximum and minimum values of $\delta\chi$, $\delta\chi_{+}$ and
$\delta\chi_{-}$. The amplitude of the oscillatory part of $\delta\chi$ is
then obtained via its definition
$G_{\delta\chi}=\frac{1}{2}(\delta\chi_{+}-\delta\chi_{-})$. Finally, Eq. (28)
can then be used to find the amplitude for $\chi_{\text{eff}}$.

###  IV.5 Evolving the phase of $\delta\chi$ and $\chi_{\rm eff}$

With both the averages and amplitudes of $\delta\chi$ and $\chi_{\rm eff}$ in
hand, we are only left to determine the evolution of the phase $\psi$. For
reference we re-express the relevant derivative, Eq. (40), in the most
immediately useful form:

$\begin{split}\left(\frac{d\delta\chi}{dt}\right)^{2}=\frac{9y^{11}}{4}X_{3}A_{\delta\chi}^{2}\left(\delta\chi-\delta\chi_{+}\right)\left(\delta\chi-\delta\chi_{-}\right)\left(\delta\chi-\frac{\delta\chi_{3}}{\delta\mu}\right)\;,\end{split}$
(43)

where the roots of the cubic equation are ordered such that
$\delta\chi_{-}\leq\delta\chi_{+}\leq\delta\chi_{3}/\delta\mu$. To get an
expression for the phase evolution $\psi$, we substitute Eq. (23) into the
equation above and use the definition in Eq. (25) to simplify the expression.
The resulting equation is

$\begin{split}\left(\frac{d\psi}{dt}\right)^{2}=\frac{9y^{11}}{4}X_{3}A_{\delta\chi}^{2}\left(\frac{\delta\chi_{3}}{\delta\mu}-\left<\delta\chi\right>-G_{\delta\chi}\sin(\psi)\right).\end{split}$
(44)

This relation is not approximate, all of the nonlinearity in the original
equation is accounted for. Instead of this exact result, it is preferable to
use an averaged version of this equation so that the phase evolves at a fixed
rate. Intuitively, the times at which the phase would have evolved
faster/slower than its average rate correspond to neglected higher order modes
in its Fourier series. (One could include these higher order modes in the
definition in (23) to account for this; we will not do that here.) To simplify
the final expression of this average we first define

$\begin{split}\left<A_{\delta\chi}\right>&=1+yA_{\delta\chi,\delta\chi}\left<\delta\chi\right>+yA_{\delta\chi,\chi_{\text{eff}}}\left<\chi_{\text{eff}}\right>,\\\
G_{A}&=yA_{\delta\chi,\delta\chi}G_{\delta\chi}+yA_{\delta\chi,\chi_{\text{eff}}}G_{\chi_{\text{eff}}}.\end{split}$
(45)

The average rate of change of $\psi$ is then given by

$\begin{split}\left(\frac{d\psi}{dt}\right)^{2}\approx\frac{9y^{11}}{8}X_{3}(\delta\chi_{3}/\delta\mu-\left<\delta\chi\right>)(2\left<A_{\delta\chi}\right>^{2}+G_{A}^{2})-G_{\delta\chi}G_{A}\left<A_{\delta\chi}\right>).\end{split}$
(46)

In the absence of radiation reaction, all quantities in the right hand side
are constant and thus this can be integrated exactly to give:

$\begin{split}\psi(t)\approx\frac{3y^{11/2}}{2^{3/2}}\sqrt{X_{3}(\delta\chi_{3}/\delta\mu-\left<\delta\chi\right>)(2\left<A_{\delta\chi}\right>^{2}+G_{A}^{2})-G_{\delta\chi}G_{A}\left<A_{\delta\chi}\right>)}\;\;t+\psi(0).\end{split}$
(47)

###  IV.6 Evolving $\phi_{z}$ and $\zeta$

Given $\chi$ and $\chi_{\rm eff}$, we now turn to the slightly more
complicated evolutions of $\phi_{z}$ and $\zeta$ as these variables have both
an evolving periodic and an evolving secular part, also in the absence of
radiation. The derivative of $\phi_{z}$ can be given approximately as (see
appendix B):

$\begin{split}\frac{d\phi_{z}}{dt}\approx&\frac{Jy^{6}}{2}\left\\{Q_{4}+Q_{5}\sin(\psi)+\frac{H_{0}+H_{1}\sin(\psi)+H_{2}\sin^{2}(\psi)+H_{3}\sin^{3}(\psi)}{(1+H_{-}\sin(\psi))(1+H_{+}\sin(\psi))}\right\\}.\end{split}$
(48)

While this can be integrated exactly, it is useful to split $\phi_{z}$ into
its secular and periodic part. To do this, we rewrite the derivative of
$\phi_{z}$ to separate terms with differing behavior:

$\begin{split}\frac{d\phi_{z}}{dt}\approx&\frac{Jy^{6}}{2}\left\\{Q_{4}+\frac{H_{2}H_{+}H_{-}-H_{3}H_{-}-H_{3}H_{+}}{{H_{+}}^{2}{H_{-}}^{2}}\right\\}+\frac{Jy^{6}}{2}\left\\{Q_{5}+\frac{H_{3}}{H_{+}H_{-}}\right\\}\sin(\psi)\\\
&+\frac{Jy^{6}}{2}\left\\{\frac{H_{0}{H_{+}}^{3}-H_{1}{H_{+}}^{2}+H_{2}{H_{+}}-H_{3}}{(H_{+}-H_{-}){H_{+}}^{2}(1+H_{+}\sin(\psi))}\right\\}-\frac{Jy^{6}}{2}\left\\{\frac{H_{0}{H_{-}}^{3}-H_{1}{H_{-}}^{2}+H_{2}{H_{-}}-H_{3}}{(H_{+}-H_{-}){H_{-}}^{2}(1+H_{-}\sin(\psi))}\right\\}.\end{split}$
(49)

The first term is purely secular, the second term is purely periodic and the
last two terms are mixed. To split the last two terms into secular and
periodic parts, we precession average these terms. The precession average then
contributes to the secular part, while the remaining part contributes to the
periodic behavior. The resulting secular part is

$\begin{split}\frac{d\left<\phi_{z}\right>}{dt}\approx\Phi_{0}+\Phi_{+}+\Phi_{-},\end{split}$
(50)

and the periodic part is

$\begin{split}\frac{d\phi_{z}}{dt}-\frac{d\left<\phi_{z}\right>}{dt}\approx\Phi_{s}\sin(\psi)+\Phi_{+}\left(\frac{\sqrt{1-{H_{+}}^{2}}}{1+H_{+}\sin(\psi)}-1\right)+\Phi_{-}\left(\frac{\sqrt{1-{H_{-}}^{2}}}{1+H_{-}\sin(\psi)}-1\right),\end{split}$
(51)

where $\Phi_{0}$, $\Phi_{s}$, $\Phi_{+}$, and $\Phi_{-}$ are independent of
$\phi_{z}$ and $\zeta$ and defined in App. B. The periodic part can be
integrated exactly to yield

$\begin{split}\phi_{z}-\left<\phi_{z}\right>\approx&-\frac{\Phi_{s}}{\dot{\psi}}\cos(\psi)+\frac{\Phi_{+}}{\dot{\psi}}\left(2\arctan\left(\frac{\tan(\psi/2)+H_{+}}{\sqrt{1-{H_{+}}^{2}}}\right)-\psi-\arcsin(H_{+})\right)\\\
&+\frac{\Phi_{-}}{\dot{\psi}}\left(2\arctan\left(\frac{\tan(\psi/2)+H_{-}}{\sqrt{1-{H_{-}}^{2}}}\right)-\psi-\arcsin(H_{-})\right).\end{split}$
(52)

The constant terms (involving arcsin) come from the constant of integration,
and are used to ensure that the average of the periodic part is zero.

The evolution for $\zeta$ is solved in a similar manner. Its derivative is
given by $d\zeta/dt=-\cos(\theta_{L})d\phi_{z}/dt$. Substituting Eq. (49) and
the expression for $\cos(\theta_{L})$ gives

$\begin{split}\frac{d\zeta}{dt}\approx-\left(\Theta_{0}+\Theta_{s}\sin(\psi)\right)\left(\Phi_{0}+\Phi_{s}\sin(\psi)+\frac{\Phi_{+}\sqrt{1-{H_{+}}^{2}}}{1+H_{+}\sin(\psi)}+\frac{\Phi_{-}\sqrt{1-{H_{-}}^{2}}}{1+H_{-}\sin(\psi)}\right)\end{split}$
(53)

Splitting this into a secular and periodic part, we find that the secular part
is

$\begin{split}\frac{d\left<\zeta\right>}{dt}\approx-\Theta_{0}(\Phi_{0}+\Phi_{+}+\Phi_{-})-\frac{\Theta_{s}\Phi_{s}}{2}-\frac{\Theta_{s}\Phi_{+}}{H_{+}}\left(\sqrt{1-{H_{+}}^{2}}-1\right)-\frac{\Theta_{s}\Phi_{-}}{H_{-}}\left(\sqrt{1-{H_{-}}^{2}}-1\right),\end{split}$
(54)

and the periodic part is

$\begin{split}\frac{d\zeta}{dt}-\frac{d\left<\zeta\right>}{dt}\approx&-\Theta_{0}\left(\frac{d\phi_{z}}{dt}-\frac{d\left<\phi_{z}\right>}{dt}\right)-\Theta_{s}\Phi_{0}\sin(\psi)+\frac{\Theta_{s}\Phi_{s}}{2}\cos(2\psi)\\\
&+\frac{\Theta_{s}\Phi_{+}}{H_{+}}\left(\frac{\sqrt{1-{H_{+}}^{2}}}{1+H_{+}\sin(\psi)}-1\right)+\frac{\Theta_{s}\Phi_{-}}{H_{+}}\left(\frac{\sqrt{1-{H_{-}}^{2}}}{1+H_{-}\sin(\psi)}-1\right).\end{split}$
(55)

The periodic part can be integrated exactly to yield

$\begin{split}\zeta-\left<\zeta\right>\approx&-\Theta_{0}\left(\phi_{z}-\left<\phi_{z}\right>\right)-\frac{\Theta_{s}\Phi_{0}}{\dot{\psi}}\cos(\psi)+\frac{\Theta_{s}\Phi_{s}}{4\dot{\psi}}\sin(2\psi)\\\
&+\frac{\Theta_{s}\Phi_{+}}{H_{+}}\left(2\arctan\left(\frac{\tan(\psi/2)+H_{+}}{\sqrt{1-{H_{+}}^{2}}}\right)-\psi-\arcsin(H_{+})\right)\\\
&+\frac{\Theta_{s}\Phi_{-}}{H_{+}}\left(2\arctan\left(\frac{\tan(\psi/2)+H_{-}}{\sqrt{1-{H_{-}}^{2}}}\right)-\psi-\arcsin(H_{-})\right).\end{split}$
(56)

##  V Adding Radiation Reaction

We have discussed the complete conservative dynamics of spins in Sec. IV. With
gravitational radiation reaction included, the orbital frequency increases as
a function of time, which can be obtained as an expansion in Post-Newtonian
orders. The total angular momentum $J$ and the average part of $\chi_{\rm
eff},\delta\chi$ become time-dependent, with corresponding evolution equations
discussed below. We do not discuss the evolution of the secular part of
$\phi_{z}$ and $\zeta$ here as the previous expressions for their derivatives
can simply be evaluated and evolved on the radiation reaction timescale, but
with updated values for the averages of $J$, $\delta\chi$ and
$\chi_{\text{eff}}$ at each step.

###  V.1 Evolving $J$

Restating it here for simplicity, the radiation reaction equation for $J$ is

$\begin{split}\frac{dJ}{dt}=\frac{-L}{2Jy}\frac{dy}{dt}\left(2L+\delta\mu\delta\chi+\chi_{\text{eff}}\right).\end{split}$
(57)

Similarly to before, this can be expressed in the “m=0” approximation,
expanding the solution into a secular part and periodic part. Here, the
solution is slightly out of phase with $\delta\chi$/$\chi_{\text{eff}}$, so an
additional term will be required to account for this. Thus, the solution for
$J$ should take the following form

$\begin{split}J\approx\left<J\right>+G_{J,s}\sin(\psi)+G_{J,c}\cos(\psi)\end{split}$
(58)

and we are left to determine $\left<J\right>,G_{J,s}$ and $G_{J,c}$. Because
the secular part is much smaller than the periodic part,
$G_{J,s}/\left<J\right>$ and $G_{J,c}/\left<J\right>$ are assumed to be small.
Substituting Eqs. (23) and (58) into Eq. (72) and applying this approximation
gives

$\begin{split}\frac{d\left<J\right>}{dt}+G_{J,s}\frac{d\psi}{dt}\cos(\psi)-G_{J,c}\frac{d\psi}{dt}\sin(\psi)&\approx\frac{-L}{2\left<J\right>y}\frac{dy}{dt}\left(2L+\delta\mu\delta\chi+\chi_{\text{eff}}\right)\\\
&+\frac{L}{2\left<J\right>^{2}y}\frac{dy}{dt}\left(2L+\delta\mu\delta\chi+\chi_{\text{eff}}\right)\left(G_{J,s}\sin(\psi)+G_{J,c}\cos(\psi)\right).\end{split}$
(59)

$d\left<J\right>/dt$ is found by averaging both sides of this equation with
respect to $\psi$ over one cycle:

$\begin{split}\frac{d\left<J\right>}{dt}=\frac{-L}{2\left<J\right>y}\frac{dy}{dt}\left(2L+\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>\right)+\frac{L}{4\left<J\right>^{2}y}\frac{dy}{dt}G_{J,s}\left(\delta\mu
G_{\delta\chi}+G_{\chi_{\text{eff}}}\right),\end{split}$ (60)

where we have treated $d\psi/dt$ as a constant for the purpose of the
precession average, because it only evolves on the radiation reaction
timescale. Another, alternative way of thinking about this is that the
correction to this equation resulting from accounting for the change in
$d\psi/dt$ will always be multiplied by one of the two amplitudes of $J$, this
correction being small and the amplitudes of $J$ also being small means the
resulting term is sufficiently small so that it can safely be neglected. A
system of equations for $G_{J,s}$ and $G_{J,c}$ can be obtained by multiplying
by either $\sin(\psi)$ or $\cos(\psi)$ and averaging over $\psi$:

$\displaystyle
G_{J,s}\frac{d\psi}{dt}=\frac{L}{2\left<J\right>^{2}y}\frac{dy}{dt}G_{J,c}\left(2L+\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>\right)$
(61) $\displaystyle-
G_{J,c}\frac{d\psi}{dt}=\frac{-L}{2\left<J\right>y}\frac{dy}{dt}\left(\delta\mu
G_{\delta\chi}+G_{\chi_{\text{eff}}}\right)+\frac{L}{2\left<J\right>^{2}y}\frac{dy}{dt}G_{J,s}\left(2L+\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>\right).$
(62)

Eqs. (60)-(62) form a system of equations for the desired unknowns:
$d\left<J\right>/dt$, $G_{J,s}$, and $G_{J,c}$. In previous works
Chatziioannou:2017tdw ; Klein:2021jtd  the value of $J$ is then used to
calculate the roots of $\delta\chi$, instead here we use this derivative of
the average of $J$ to find the derivatives of the averages of $\delta\chi$ and
$\chi_{\text{eff}}$.

###  V.2 Evolving $\left<\delta\chi\right>$ and
$\left<\chi_{\text{eff}}\right>$

After using the results of Sec. IV.3 we can obtain the averages
$\left<\delta\chi\right>$ and $\left<\chi_{\text{eff}}\right>$ at the initial
time. In the presence of radiation reaction we must then evolve these
averages. First we evolve the roots ($\delta\chi_{+}$, $\delta\chi_{-}$,
$\chi_{\text{eff},+}$, $\chi_{\text{eff},-}$) and then use the relation
between these roots and the averages (Eq. (24) to evolve the averages. To find
these derivatives, we start with their definitions as the roots of the
following equations:

$\displaystyle\delta\mu{\delta\chi_{+}}^{3}+(B_{0}+B_{1}\chi_{\text{eff},+}){\delta\chi_{+}}^{2}+(C_{0}+C_{1}\chi_{\text{eff},+}+C_{2}{\chi_{\text{eff},+}}^{2}){\delta\chi_{+}}$
$\displaystyle+(D_{0}+D_{1}\chi_{\text{eff},+}+D_{2}{\chi_{\text{eff},+}}^{2}+D_{3}{\chi_{\text{eff},+}}^{3})=0,$
(63)
$\displaystyle\delta\mu{\delta\chi_{-}}^{3}+(B_{0}+B_{1}\chi_{\text{eff},-}){\delta\chi_{-}}^{2}+(C_{0}+C_{1}\chi_{\text{eff},-}+C_{2}{\chi_{\text{eff},-}}^{2}){\delta\chi_{-}}$
$\displaystyle+(D_{0}+D_{1}\chi_{\text{eff},-}+D_{2}{\chi_{\text{eff},-}}^{2}+D_{3}{\chi_{\text{eff},-}}^{3})=0,.$
(64)

Using the implicit function theorem we can get two equations, one satisfied by
the positive roots $\delta\chi_{+}$ and $\chi_{\text{eff},+}$ (the derivative
of Eq. (V.2)), and one satisfied by the negative roots $\delta\chi_{-}$ and
$\chi_{\text{eff},-}$ (the derivative of Eq. (V.2)):

$\displaystyle
B_{0}^{\prime}\delta\chi_{+}^{2}+\left(C_{0}^{\prime}+C_{1}^{\prime}\chi_{\text{eff},+}\right)\delta\chi_{+}+\left(D_{0}^{\prime}+D_{1}^{\prime}\chi_{\text{eff},+}+D_{2}^{\prime}\chi_{\text{eff},+}\right)+\left(3\delta\mu\delta\chi_{+}^{2}+2B\delta\chi_{+}+C\delta\chi_{+}\right)\delta\chi_{+}^{\prime}$
$\displaystyle+\left(B_{1}\chi_{\text{eff},+}\delta\chi_{+}^{2}+(C_{1}+2C_{2}\chi_{\text{eff},+})\delta\chi_{+}+(D_{1}+2D_{2}\chi_{\text{eff},+}+3D_{3}\chi_{\text{eff},+}^{2})\right)\chi_{\text{eff},+}^{\prime}=0,$
(65) $\displaystyle
B_{0}^{\prime}\delta\chi_{-}^{2}+\left(C_{0}^{\prime}+C_{1}^{\prime}\chi_{\text{eff},-}\right)\delta\chi_{-}+\left(D_{0}^{\prime}+D_{1}^{\prime}\chi_{\text{eff},-}+D_{2}^{\prime}\chi_{\text{eff},-}\right)+\left(3\delta\mu\delta\chi_{-}^{2}+2B\delta\chi_{-}+C\delta\chi_{-}\right)\delta\chi_{-}^{\prime}$
$\displaystyle+\left(B_{1}\chi_{\text{eff},-}\delta\chi_{-}^{2}+(C_{1}+2C_{2}\chi_{\text{eff},-})\delta\chi_{-}+(D_{1}+2D_{2}\chi_{\text{eff},-}+3D_{3}\chi_{\text{eff},-}^{2})\right)\chi_{\text{eff},-}^{\prime}=0,$
(66)

where the primes denote time derivatives. These can be obtained using the
definitions of the coefficients, $B$, $C$, and $D$ in conjunction with the
time derivatives of $J$ and $y$, as an example $B_{0}^{\prime}$ is given as

$\begin{split}B_{0}^{\prime}=\frac{\partial B_{0}}{\partial
J}\frac{d\left<J\right>}{dt}+\frac{\partial B_{0}}{\partial
y}\frac{dy}{dt}\end{split}$ (67)

A second set of relations follow from the derivatives directly, giving a
complete set of equations:

$\displaystyle\chi_{\text{eff},+}^{\prime}$
$\displaystyle=\frac{yA_{\chi_{\text{eff}},\delta\chi}\delta\chi_{+}+yA_{\chi_{\text{eff}},\chi_{\text{eff}}}\chi_{\text{eff},+}}{1+yA_{\delta\chi,\delta\chi}\delta\chi_{+}+yA_{\delta\chi,\chi_{\text{eff}}}\chi_{\text{eff},+}}\delta\chi_{+}^{\prime},$
(68) $\displaystyle\chi_{\text{eff},-}^{\prime}$
$\displaystyle=\frac{yA_{\chi_{\text{eff}},\delta\chi}\delta\chi_{-}+yA_{\chi_{\text{eff}},\chi_{\text{eff}}}\chi_{\text{eff},-}}{1+yA_{\delta\chi,\delta\chi}\delta\chi_{-}+yA_{\delta\chi,\chi_{\text{eff}}}\chi_{\text{eff},-}}\delta\chi_{-}^{\prime}.$
(69)

The derivative of the averages are then straightforwardly obtained from the
solutions to these via:

$\displaystyle\frac{d\left<\delta\chi\right>}{dt}$
$\displaystyle=\frac{1}{2}\left(\delta\chi_{+}^{\prime}+\delta\chi_{-}^{\prime}\right)$
(70a) $\displaystyle\frac{d\left<\chi_{\text{eff}}\right>}{dt}$
$\displaystyle=\frac{1}{2}\left(\chi_{\text{eff},+}^{\prime}+\chi_{\text{eff},-}^{\prime}\right)$
(70b)

In the absence of radiation reaction the derivatives of $y$ and $J$ are zero
and consequently the derivatives of the coefficients ($B_{0}^{\prime}$,
$B_{1}^{\prime}$, $C_{0}^{\prime}$, etc) are zero. The solution to the
resulting system of equations is
$\delta\chi_{+}^{\prime}=\delta\chi_{-}^{\prime}=\chi_{\text{eff},+}^{\prime}=\chi_{\text{eff},-}^{\prime}=0$.
Thus, we nicely recover that in the absence of radiation reaction the averages
do not evolve.

###  V.3 Summarizing

At this point, we have discussed all the equation necessary for describing the
spin dynamics under the influence of radiation reaction. Our method is
summarized in Fig. 1 in the form of a flow chart. The general idea can be
broken up into several blocks, which can be described generally as:
initializing the averages (red), updating the averages (yellow), and
calculating the periodic part (green). The first two deal with the secular
parts, while the latter is critical for accounting for the effects of
precession. Once the evolution of these key quantities are known, the waveform
is generated using standard techniques (see e.g. (Klein:2021jtd, , Sec. IV)).
Following Chatziioannou:2017tdw , we refer to the resulting wavefor as a
frequency domain waveform because, while the equations are evolved in the time
domain, these time domain solutions can be transformed directly to frequency
domain waveforms via the method of shifted uniform asymptotics (as opposed to
constructing a time domain waveform then transforming this to the frequency
domain). We do not present the equations to produce this waveform here because
they are identical to those found in (Klein:2021jtd, , Sec. IV). In the
analysis here, we only use the leading term in $dy/dt=(32\eta/5)y^{9}$ unless
otherwise specified. When constructing realistic waveforms for actual analysis
of gravitational waves, this should be updated with the most current PN
expression for $dy/dt$. Thus, to include the effects of higher PN results for
the nonspinning part one should adjust $dy/dt$ as well as modify the
expressions for the nonspinning part appearing in the waveform generation.

Start Initialize $\left<\delta\chi\right>(0)$,
$\left<\chi_{\text{eff}}\right>(0)$, and $\psi(0)$ (Section IV.3) Initialize
$\left<J\right>$(0) (Section V.1) Initialize $\left<\phi_{z}\right>(0)$ and
$\left<\zeta\right>$(0) (Section IV.6) Evolve $\psi$ (Section IV.5) Evolve
$\left<J\right>$ (Section V.1) Evolve $\left<\delta\chi\right>$ and
$\left<\delta\chi\right>$ (Section V.2) Evolve $\left<\phi_{z}\right>$ and
$\left<\zeta\right>$ (Section IV.6) Get periodic part of $\delta\chi$ and
$\chi_{\text{eff}}$ (Section IV.4) Get periodic part of $J$ (Section V.1) Get
periodic part of $\phi_{z}$ and $\zeta$ (Section IV.6) Are we finished? noend
yes Figure 1: This flow chart summarizes the key steps in solving for the
precession equations of compact objects (with possibly $\kappa\neq 1$) on
circular orbits. The parts in red initialize the averages, the parts in yellow
update the averages, and the parts in green evolve the periodic parts of the
solutions.

##  VI Comparing models

In this section, we first establish the validity of our method in Sec. VI.1
before we compare waveforms of various binaries in Sec. VI.2. In Sec. VI.1 we
compute the mismatch between the waveform constructed according to the
discussion in Sec. IV,V and the waveform following numerical evolution of
spins. The latter section describe the key results of this paper: its shows
that the quadrupole $\kappa_{i}$ has an observable imprint on the waveform, at
least for some of the spin/orbital configurations.

###  VI.1 Model validation

To evaluate the effectiveness of the proposed method, we compare the method to
the fully numerically evolved precession equations Eq. (9). Our method is
better than simply numerically evolving the precession equations directly
because of the timescales on which they can be accurately evolved (while still
accounting for the full effects of precession). To numerically evolve the
precession equations one must evolve them on the precession timescale,
otherwise they run the risk of diverging. These divergences usually occur as a
result of one point lying outside the upper or lower bound on $\delta\chi$ or
$\chi_{\text{eff}}$, the solution then diverges secularly. To avoid this
behavior, one typically has to use an overly expensive number of points in the
evolution to keep the error low. Our solution has no such issues. Since we are
evolving the average values, there is no chance of such an issue happening in
the first place and it can be accurately evolved on the radiation reaction
timescale. Because of the different timescales on which the system must be
evolved, our method takes on average 5.1ms evolving from 10-100Hz with a
single CPU for the system configurations in this section, while evolving the
precession equations numerically (while keeping enough points to avoid the
aforementioned divergence) takes on average 54ms evolving from 10-100Hz with a
single CPU for the same configurations. For high SNR events with better
detector sensitivity, the waveform will be longer so that this speed-up factor
will be larger. Our method is only constrained by how many frequency points
one wants to evaluate, while the full numerical evolution is constrained by
the precession timescale.

As concrete examples for comparison, we examine two systems: a binary
comprised of two neutron stars and a binary consisting of one neutron star and
one black hole. For these comparisons we will use $\chi_{1}=0.4$ and
$\chi_{2}=0.7$ for the spins, we set the angles of $\vec{s}_{1}$ and
$\vec{s}_{2}$ relative to $\hat{L}$ as $\theta_{s1}=\pi/20$, $\phi_{s1}=0$,
$\theta_{s2}=\pi/4$, $\phi_{s2}=\pi/10$. The only difference between the NSNS
and NSBH systems we analyze are the masses and spin-induced quadrupole moment
constants:

  * •

For NSNS binaries, we will use
$M_{1}=2.6M_{\odot},M_{2}=1.5M_{\odot},\kappa_{1}=2.5$, $\kappa_{2}=3.5$.

  * •

For NSBH binaries, we will use
$M_{1}=23.0M_{\odot},M_{2}=2.6M_{\odot},\kappa_{1}=1.0$ and $\kappa_{2}=2.5$.

Here we have chosen values of $\kappa$ that are smaller than the typical range
for neutron stars ($\kappa\sim 4-8$ Steinhoff:2021dsn ) to be conservative in
examining whether this is potentially measurable. To facilitate the
comparison, the tidal effects of NSs are not included here, although their
implementation is straightforward to add.

To validate the effectiveness of our method, we start by looking at several
spin variables as a function of frequency to gain insight into the error
resulting from the approximations we made to see how we can improve the
accuracy in the future. To facilitate this, we separate the error into secular
and periodic errors. We start with $\delta\chi$ and $\chi_{\text{eff}}$, since
they form two of our three dynamical variables, then we examine the variables
more relevant to waveform generation: $\phi_{z}$ and $\theta_{L}$ (we exclude
$\zeta$ because its evolution is similar to that of $\phi_{z}$). Finally we
randomize spin configurations and compute the fidelity/mismatch for various
NSBH/NSNS systems, using the mismatch as a measure of the accuracy of our spin
evolution scheme.

In the logarithmic error plots below, one can easily determine whether the
secular evolution or the periodic evolution dominates the error budget. The
presence of a series of sharp dips means that the error approaches zero
repeatedly, implying that the periodic part is the main source of error. When
these sharp dips are absent, the secular error dominates. When there is a
transition from the periodic part to the secular part being the dominant
source of error, every other peak will gradually decrease in magnitude until
there are no more sharp dips. This is because two consecutive peaks correspond
to periods where the quantity is below/above what it should be. As an example,
if the secular error diverges to positive values, the quantity will gradually
spend less time being “too negative” and more time being “too positive”, until
it is never “too negative”, as a result the peaks that correspond to being
“too negative” will diminish in size until they vanish altogether.

We begin by comparing the numerical evolution of $\delta\chi$ and
$\chi_{\text{eff}}$ to our semi-analytic approach in Fig. 2. First, note that
there is good agreement between the precession frequency of the numerically
evolved spins and the spins evolved using the methods in this paper. The error
is largest for $\chi_{\rm eff}$ in the NSNS binary when the frequencies are
high, but in most cases, the difference is not visible. Obtaining an accurate
precession frequency is critical in reducing the error of the final solution,
because small phase errors produce secular errors of the same order in the
oscillation amplitude, which can introduce large periodic errors in the final
result. Second, while the agreement between the numerical and our method is
good at the extrema of the oscillations, there is some disagreement in
between. This error is a result of approximating the solutions by their first
mode only (i.e., the “m=0” approximation). Finally, there is a slight error in
the amplitude and average of the oscillations. These latter error sources are
the main contributions to the error in the final results. To see this, we
start by examining $\phi_{z}$.

(a) NSNS binary

(b) NSBH binary

(c) NSNS binary

(d) NSBH binary

Figure 2: Plot of $\chi_{\text{eff}}/M^{2}$ vs frequency for the NSNS (NSBH)
binary in panel a (b), and $\delta\chi/M^{2}$ for the NSNS (NSBH) binary in
panel c (d). In each plot the black curve is the numerically evolved system,
while the orange curves are the solution found using the method outlined in
this paper. The second plot in each panel is the absolute error.

The absolute error for $\phi_{z}$ at low frequencies is $O(10^{-2})$ radians
and at higher frequencies $O(10^{-1})$ radians for the NSNS system, while for
the NSBH system is maximally $O(10^{-2})$ radians, as shown in Fig. 3. This
error is dominated by its secular part, because of the complete lack of sharp
dips (except at the lowest frequency, but this corresponds to the frequency at
which we initialized the system). The main source of this secular error is not
obvious; there are two potential candidates. The first is the “m=0”
approximation, by its very nature it produces errors of order $O(m)$ (which is
equivalent to $O(y^{2})$). The second is the observed secular error in
$\delta\chi$ and $\chi_{\text{eff}}$, which are also suppressed by $O(y^{2})$.
Since both error sources contribute at the same order, it is difficult to
ascertain which one dominates this error budget.

(a) NSNS Binary $\phi_{z}$

(b) NSBH Binary $\phi_{z}$

Figure 3: Plots of the $\phi_{z}$ for the NSNS binary (a) and the NSBH binary
(b), in black is the numerically evolved system, while orange is the system
evolved with our method. The top panel is the actual values of $\phi_{z}$
while the bottom panels are the absolute error between the value produced by
the numerical evolution and that produced by our method.

Finally, the error for $\theta_{L}$ is plotted in Fig. 4. This error is
periodic and consequently one might be tempted to assume that its main source
is the non-inclusion of higher terms in the Fourier series directly. However,
this is not the case, as the error has the same frequency as the precession,
if they were from the non inclusion of these terms directly they would
correspond to oscillations at double or triple the frequency. Hence, the major
contributing source of error is small errors in the amplitude of $\delta\chi$.

(a) NSNS Binary $\theta_{L}$

(b) NSBH Binary $\theta_{L}$

Figure 4: Plots of $\theta_{L}$ as a function of frequency for the NSNS binary
(a) and the NSBH binary (b). The black line shows the numerically evolved
system, while orange shows the system evolved with our method. The top panel
displays the actual values of $\theta_{L}$ while the bottom panels are the
absolute error between the value produced by the numerical evolution and that
produced by our method.

To evaluate the accuracy of this waveform we evaluate the mismatch for two
sets of systems:

The first set to evaluate the effectiveness for general NSNS systems, we
evenly space the first mass between $[1.8,2.6]M_{\odot}$, and for the second
we distribute its mass evenly between $1.0M_{\odot}$ and $0.9m_{1}$ (to ensure
that the mass ratio is not equal to one). We set the first spin to be
$\chi_{1}=0.7$, and the second to be $\chi_{2}=0.2$, we then evenly distribute
the angles $\theta_{s1}$, $\theta_{s2}$, and $\phi_{s2}$ every $20^{\circ}$
over the sphere, evolving every system from 10-100Hz.

The second set of systems is more restricted in order to explore how much
error is produced for models similar to GW190814, since this is the system of
interest to us. In particular, we restrict the masses to $m_{1}=23M_{\odot}$
and $m_{2}=2.6M_{\odot}$, and the spins to $\chi_{1}=\chi_{2}=0.6$. Similarly
we distribute the angles $\theta_{s1}$, $\theta_{s2}$, and $\phi_{s2}$ every
$20^{\circ}$ over the sphere, evolving every system from 10-100Hz.

In the following we use the overlap, defined here as

$\begin{split}{\rm
overlap}(h_{1},h_{2})&=\frac{(h_{1},h_{2})}{\sqrt{(h_{1},h_{1})(h_{2},h_{2})}},\\\
(h_{1},h_{2})&=4\rm
Re\int_{f_{min}}^{f_{max}}h_{1}(f)h_{2}^{*}(f)df,\end{split}$ (71)

where the waveform phase is aligned at the initial frequency. This is done,
since it will be a more conservative estimate of the accuracy than performing
matched filterring. The cumulative distribution functions (CDF) of the
mismatch (one minus the overlap) between the waveform produced via our method
and the waveform produced via the full numerical evolution are shown in Fig.
5.

(a) NSNS Binary

(b) NSBH Binary

Figure 5: Plots of the CDF for mismatch between the numerically evolved
NSNS/NSBH systems and the same systems evolved with our method.

For the set of binary neutron star systems we examine we find that only
1.35%(7.44%) have a fidelity worse than 0.965(0.994). This is the same order
of magnitude as in Chatziioannou:2017tdw . We have investigated several of the
higher mismatch cases, and find that the main cause of these high mismatch
cases was the non-inclusion of higher terms in the Fourier series of our “m=0”
approximation. This is because in these higher mismatch cases, $\theta_{L}$
drifts close to zero at its minimum. A small error in the value of
$\theta_{L}$ at this point then introduces larger effects in the secular
evolution of $\phi_{z}$. When the system precesses away from $\theta_{L}$
being close to zero, this then poses an issue and in turn contributes to the
mismatch. Adding more terms in the approximation would make this error
smaller, in principal allowing us to reduce this mismatch. In the NSBH case no
systems we examined have overlap worse than 0.999%.

###  VI.2 Differences between BNS/NSBH and BBH Waveforms

Having established the excellent accuracy of our method with mismatch less
than 3.5% for 98.65% in the more extreme case of NSNS systems (and an order of
magnitude better accuracy in BHNS systems), here we study the qualitative
effects of the spin-induced quadrupole moment on the evolution. In order to do
so, we compare the two systems in Sec. VI.1 to a binary with two black holes
(BHBH) with the same parameters, except with $\kappa_{1}=1=\kappa_{2}$.

Here we follow a similar procedure as in the previous section, using the same
two sample configurations to evolve the NSNS/NSBH systems and then evolve the
same system as though they were a BHBH system (with the other parameters all
the same).

First, the effect of the quadrupole moment on $\chi_{\text{eff}}$ and
$\delta\chi$ is shown in Fig. 6. Qualitatively the largest difference is the
introduction of oscillations to $\chi_{\text{eff}}$, which are absent in the
BHBH case. So one would expect the corresponding correction to be the same PN
order as the amplitude of $\chi_{\text{eff}}$’s oscillations, which are
$O(y^{2})$. If the amplitude alone were modified this would introduce a
relatively small periodic error. As we will see, the main contributor to the
error is the modification to the precession frequency, observed via the
difference in oscillation frequency between the blue and black curves. Even
small changes in frequency can introduce large differences in the final
result. Also note the absence of large secular differences in either
$\chi_{\text{eff}}$ or $\delta\chi$, this will be important for later.

(a) NSNS binary

(b) NSBH binary

(c) NSNS binary

(d) NSBH binary

Figure 6: Plot of $\chi_{\text{eff}}/M^{2}$ vs frequency for the NSNS (NSBH)
binary in panel a (b), and $\delta\chi/M^{2}$ for the NSNS (NSBH) binary in
panel c (d). In each plot the black curve is the numerically evolved neutron
star system, while the blue curves are the corresponding BBH systems.

Next, we turn to the angle $\phi_{z}$ in Fig. 7. The difference between the
NSNS and BHBH case is large enough that it is directly visible, and the
difference accumulates to several radians. The source of this secular effect
is the corrections to the terms in the derivative of the average, Eq. (50),
involving $Q_{11}$-$Q_{16}$. Thus, the effects of the spin-induced quadrupole
moment are important and introduce potentially huge errors to the waveform if
left unaccounted for, even for smaller quadrupole moment constants between
1.1-1.5

The difference between the NSBH and BHBH systems is not secular but periodic
in nature, characterized by the sharp, repeating dips. This is because the
secular part for a NSBH system is suppressed by the mass ratio. Intuitively,
this makes sense: the larger body’s spin-induced quadrupole moment affects the
waveform more. In the NSNS case, both masses are similar, whereas the masses
can differ significantly in a NSBH binary. For the NSBH case we considered,
the mass ratio is large so that the secular part of the evolution of
$\phi_{z}$ is little affected by the quadrupole moment of the light neutron
star.

(a) NSNS Binary $\phi_{z}$

(b) NSBH Binary $\phi_{z}$

Figure 7: Plots of the $\phi_{z}$ for the NSNS binary (left) and the NSBH
binary (right), in black is the numerically evolved system, while blue is the
equivalent BBH system. The top panels show the actual values of $\phi_{z}$,
while the bottom panels show the absolute difference between the NSNS/NSBH and
BHBH binaries.

Finally, we investigate the influence of the quadrupole moment on $\theta_{L}$
in Fig. 8. The quadrupole moment leaves an imprint that is periodic. This is
as expected, since we observe no significant secular imprint of $\kappa\neq 1$
on either $\chi_{\text{eff}}$ or $\delta\chi$. The only other possible
quantity that could introduce a secular effect to $\theta_{L}$ is $J$. However
the secular effect in $J$ is largely determined by that of $\chi_{\text{eff}}$
and $\delta\chi$, so if the secular effect in $\chi_{\text{eff}}$ or
$\delta\chi$ is already small, the secular effect in $J$ will also be small.
As a result, since the secular effect of the quadrupole moment is small in any
of the quantities that determine $\theta_{L}$, there is little secular effect
on $\theta_{L}$.

Since $\phi_{z}$ and $\theta_{L}$ play an important role in waveform
generation, these results imply that precession is very important for
parameter estimation.

(a) NSNS Binary $\theta_{L}$

(b) NSBH Binary $\theta_{L}$

Figure 8: Plots of the $\theta_{L}$ for the NSNS binary (a) and the NSBH
binary (b), in black is the numerically evolved system, while blue is the
equivalent BBH system. The top panels are the actual values of $\theta_{L}$
while the bottom panels are the absolute difference between the value produced
by the neutron star system and that of a BBH.

To accurately assess the measurability of these effects we will use a more
accurate measure of the mismatch that reflects what would be used in practice.
First, we define the fidelity:

$\begin{split}{\rm
F}(h_{1},h_{2})&=\max_{\phi}\frac{(h_{1},h_{2})}{\sqrt{(h_{1},h_{1})(h_{2},h_{2})}},\\\
(h_{1},h_{2})&=4\rm
Re\int_{f_{min}}^{f_{max}}\frac{h_{1}(f)h_{2}^{*}(f)}{S_{n}(f)}df,\end{split}$
(72)

where now we perform matched filtering, maximizing the fidelity over the
initial phase (or minimizing the mismatch), and we also include the advanced
LIGO zero-detuned, high power noise sensitivity LIGOScientific:2014pky . So
that we can accurately determine the magnitude of these effects with
everything accounted for, we include the higher order terms in $dy/dt$ and the
contributions of the quadrupole moment constant to the non-spinning part of
the waveform Harry:2018hke . The cumulative distribution functions (CDF) of
the mismatch (one minus the fidelity) between the BHNS/BNS waveform and the
corresponding BBH waveform are shown in Fig. fig:effect-kappa-cdf. We fix the
masses and quadrupole moment constants to those indicated in each panel in the
figure, the dimensionless spins are fixed to $\chi_{1}=0.7$ and $\chi_{2}=0.6$
for all samples. We take 1000 samples with these inputs with randomly
distributed isotropic spin angles and inclination angle.

(a) NSNS Binary

(b) NSBH Binary

Figure 9: Plots of the CDF for mismatch between the numerically evolved
NSNS/NSBH systems and the corresponding numerically evolved BHBH system.

For the neutron star binary, this effect produces a large mismatch, but for
the NSBH system this effect is more difficult to measure. The explanation is
that the secular effect is suppressed by the mass ratio, which is relatively
small in the higher mass ratio BHNS cases we examine. We can see this directly
in the second panel of Fig. fig:effect-kappa-cdf, where for gradually larger
mass ratios the average mismatch between the BNS/NSBH system and the
corresponding BBH system decreases.

By examining which cases have higher/lower mismatch, we find the quadrupole
moment constant is more measurable for somewhat aligned/anti-aligned spins.
This can be understood intuitively from the precession equations (Eq. (9)),
while the precession is larger for spins perpendicular to the orbital angular
momentum, this is not the hypothetical best place to find the _difference_ in
the precession. This is because the correction terms scale with the projection
of the spin onto the orbital angular momentum. This means that the best
systems to measure the difference between the black hole precession and
neutron star precession are systems for which their spins are misaligned
enough for the precession to be measurable, but somewhat aligned otherwise to
maximize the difference between the neutron star and black hole precession.

These results are consistent with those reported in studies of the effects on
the waveform Krishnendu:2017shb ; Krishnendu:2019tjp ; Chia:2022rwc  in
several qualitative ways. Firstly, the spin precession effect is larger for
more aligned spins, but as we have shown, if there is a slight misalignment in
the spins the effects on the precession cannot be ignored if one wants
accurate waveforms. Secondly, the effects are more noticeable when the first
spin is larger, because the secondary’s spin affects the waveform less.
Finally, for the dimensionless spin magnitudes used here, this effect is in
principle measurable.

The tidal heating effect, which comes from the energy and angular momentum on
the black hole horizon, contributes at $2.5$ PN order to the waveform if the
black hole is spinning Tagoshi:1997jy , and $4$PN if the black hole is non-
spinning Poisson:1994yf . Since neutron stars do not have tidal heating
effects, in principle the tidal heating effect can also be used to probe the
nature of mass-gap objects. However, the calculation in Alvi:2001mx  shows
that the tidal heating effect from an equal mass black hole (with high spin
$a=0.998$) binary contributes to $\sim 0.05$ gravitational wave cycle, which
is approximately $0.3$ rad in phase modulation. For unequal mass ratio
binaries, the flux absorbed by the less massive black hole $M_{2}$ is roughly
$M_{2}/M_{1}$ times smaller than the flux absorbed by $M_{1}$, if both black
holes are spinning, and $(M_{2}/M_{1})^{2}$, if both are non-spinning. As a
result, replacing $M_{2}$ by a neutron star in a GW 190814-like binary has
negligible impact on the gravitational wave phase if only the tidal heating
effect is considered.

##  VII Conclusion

We have presented a new method to construct frequency domain waveforms for
circular compact object binaries that include neutron stars. Because of the
analytical treatment of the spin evolution equations, the spin variables are
evolved on the radiation reaction timescale, which is convenient to transform
to the frequency domain. The new waveform is able to achieve more than an
order of magnitude speed-up compared to the one with fully numerical evolution
of spins on the precession timescale. For generic mass ratio and spin
configurations, the mismatch between the new waveform and the one with fully
numerical evolution of spins is $\leq 3.5\%$ for 98.65% of configurations
examined. We have also investigated the difference between BHBH, NSNS and NSBH
waveforms due to different spin-induced quadrupole moment constant values,
assuming the same component mass and spins and neglecting possible tidal
effects. We find that the maximum difference occurs for system with somewhat
aligned or anti-aligned spins. In such cases the mismatch between BHBH and BNS
waveforms can be approximately $\sim 10\%$, which is promising for detection
with Advanced LIGO/Virgo LIGOScientific:2014pky ; VIRGO:2014yos  and LIGO A+
KAGRA:2013rdx  . For other spin configurations, the difference between the
waveforms are rather small, which may require third-generation gravitational
wave detectors to measure the spin-induced quadrupole moment constant. For BNS
systems, if the neutron stars have significant spins, the mismatch may reach
$5\%$ for a good fraction of spin-configurations. This is particularly
interesting if in the future we discover BNS systems with two mass-gap
objects.

There are several avenues for further developments. First, improving this
waveform’s speed and accuracy via a semi-analytic evolution. This can be done
by expanding the averages and amplitudes of $\delta\chi$ and
$\chi_{\text{eff}}$ in a series in terms of the PN parameter $y$:

$\begin{split}\left<\delta\chi\right>&=\left<\delta\chi\right>_{0}+\left<\delta\chi\right>_{1}y+\left<\delta\chi\right>_{2}y^{2},\\\
G_{\delta\chi}&=G_{\delta\chi,0}+G_{\delta\chi,1}y+G_{\delta\chi,2}y^{2},\\\
\left<\chi_{\text{eff}}\right>&=\left<\chi_{\text{eff}}\right>_{0}+\left<\chi_{\text{eff}}\right>_{1}y+\left<\chi_{\text{eff}}\right>_{2}y^{2},\\\
G_{\chi_{\text{eff}}}&=G_{\chi_{\text{eff}},0}+G_{\chi_{\text{eff}},1}y+G_{\chi_{\text{eff}},2}y^{2}.\end{split}$
(73)

This approach has the advantage that while there is a higher up front cost to
initialize this system (requiring solving a system of linear equations), it
makes the cost of calculating the amplitudes and averages at each step
significantly lower and removes the need for root-finding altogether. The
system of equations for these coefficients would be formed from the initial
values of the amplitudes and averages, which can be found using a similar set
of derivatives of $\delta\chi$ and $\chi_{\text{eff}}$ used to calculate the
initial averages here. The rest of the equations would be formed from
derivatives of the roots of $\delta\chi$ and $\chi_{\text{eff}}$; we have
discussed in an earlier section how to find such derivatives via the implicit
function theorem. To keep high accuracy in the waveform, $\psi$ and
$\left<\phi_{z}\right>$ could potentially still be numerically evolved if
necessary.

Normally speed versus accuracy is a balancing act: improving one typically
harms the other. The approach taken in this paper provides the potential for
improving both simultaneously. The speed improvement we mentioned already and
the accuracy can be further improved by including the next terms in the
Fourier series in the “m=0” approximation. The higher order terms improve the
accuracy in two ways. First, including them make the estimation of the average
and amplitude of $\delta\chi$ and $\chi_{\text{eff}}$ more accurate, thereby
lowering the periodic effects in the error plots of $\theta_{L}$ which
correspond to $1\times$ the precession frequency. On top of this, the
inclusion of the higher order terms should remove some of the error with
double this frequency that is buried beneath this error. This improvement in
$\theta_{L}$ will also improve the accuracy of $\phi_{z}$’s evolution, since
the main source of error in $\phi_{z}$ is due to the error in the minimum of
$\theta_{L}$.

Currently, there is no waveform model that can simultaneously handle the
evolution of eccentric orbits and precession. For neutron stars binaries this
task is particularly complicated as tidal excitations may include multiple
harmonics with eccentric orbits Yang:2018bzx ; Yang:2019kmf . However, this is
(astro)physically important as some of the dynamically formed binaries may
carry nonzero eccentricity and non-negligible precession at the same time. It
is a promising direction to consider whether eccentric orbit evolution can be
incorporated into the scheme discussed here.

Finally, it will be interesting to perform Bayesian parameter estimation with
this waveform, either with the real data (such as GW 190814) or with
artificial data (detector noise plus injected signal). Given that this is
probably the best method to probe the nature of the mass-gap objects similar
to the one found in GW 190814 —-if no electromagnetic counterparts are present
— it is important to assess the ability of using these precession waveforms to
measure the spin-induced quadrupole moment constant in compact binaries with
respect to various detector sensitivities. We are currently performing such
analysis, the results will be reported in a future publication.

Acknowledgement. We thank Mohammed Khalili for reading over the manuscript and
providing many helpful comments. We thank Reed Essick and the Perimeter LSC
discussion group for useful discussions. M. L., Z. L. and H. Y. are supported
by the Natural Sciences and Engineering Research Council of Canada and in part
by Perimeter Institute for Theoretical Physics. Research at Perimeter
Institute is supported in part by the Government of Canada through the
Department of Innovation, Science and Economic Development Canada and by the
Province of Ontario through the Ministry of Colleges and Universities.

##  Appendix A Coefficients of $d\delta\chi/dt$ with behaviour separated

The coefficients in Eq. (15)-(16) can be found in Klein Klein:2021jtd , but
are repeated here:

$\displaystyle B$
$\displaystyle=\frac{y}{2\eta^{2}}\left[-2\eta(J^{2}-L^{2}-L\chi_{\text{eff}})+\delta\mu(S_{1}^{2}-S_{2}^{2})-\delta\mu^{2}(2L^{2}+S_{1}^{2}+S_{2}^{2})\right],$
(74a) $\displaystyle C$
$\displaystyle=\frac{y}{2\eta^{2}}\left\\{(1+\delta\mu^{2})\chi_{\text{eff}}(S_{1}^{2}-S_{2}^{2})+2\delta\mu\left[2L(J^{2}-L^{2}-L\chi_{\text{eff}})-(2L+\chi_{\text{eff}})(S_{1}^{2}+S_{2}^{2})-\eta
L\chi_{\text{eff}}^{2}\right]\right\\}$ (74b) $\displaystyle D$
$\displaystyle=\frac{y}{2\eta^{2}}\left\\{-2(J^{2}-L^{2}-L\chi_{\text{eff}})\left[J^{2}-L^{2}-L\chi_{\text{eff}}-2(S_{1}^{2}+S_{2}^{2})-\eta\chi_{\text{eff}}^{2}\right]\right.$
(74c)
$\displaystyle+\left.(S_{1}^{2}-S_{2}^{2})\left[\delta\mu\chi_{\text{eff}}^{2}-2(S_{1}^{2}-S_{2}^{2})\right]-\chi_{\text{eff}}^{2}(S_{1}^{2}+S_{2}^{2})\right\\}$

From these it follows that the coefficients with the behavior of
$\chi_{\text{eff}}$ separated (see Eqs. (42)) are

$\displaystyle B_{0}$
$\displaystyle=\frac{y}{2\eta^{2}}\left[-2\eta(J^{2}-L^{2})+\delta\mu(S_{1}^{2}-S_{2}^{2})-\delta\mu^{2}(2L^{2}+S_{1}^{2}+S_{2}^{2})\right]$
(75a) $\displaystyle B_{1}$ $\displaystyle=1$ (75b) $\displaystyle C_{0}$
$\displaystyle=\frac{2\delta\mu}{\eta}\left[J^{2}-L^{2}-S_{1}^{2}-S_{2}^{2}\right]$
(75a) $\displaystyle C_{1}$
$\displaystyle=\frac{y}{2\eta^{2}}\left[(1+\delta\mu^{2})(S_{1}^{2}-S_{2}^{2})-2\delta\mu(2L^{2}+S_{1}^{2}+S_{2}^{2})\right]$
(75b) $\displaystyle C_{2}$ $\displaystyle=-\delta\mu$ (75c) $\displaystyle
D_{0}$
$\displaystyle=\frac{-y}{\eta^{2}}\left(J^{2}-L^{2}-S_{1}^{2}-2S_{1}S_{2}-S_{2}^{2}\right)\left(J^{2}-L^{2}-S_{1}^{2}+2S_{1}S_{2}-S_{2}^{2}\right)$
(75a) $\displaystyle D_{1}$
$\displaystyle=\frac{2}{\eta}\left(J^{2}-L^{2}-S_{1}^{2}-S_{2}^{2}\right)$
(75b) $\displaystyle D_{2}$
$\displaystyle=\frac{y}{2\eta^{2}}\left[2\eta(J^{2}-L^{2})-(2L^{2}+S_{1}^{2}+S_{2}^{2})+\delta\mu(S_{1}^{2}-S_{2}^{2})\right]$
(75c) $\displaystyle D_{3}$ $\displaystyle=-1$ (75d)

##  Appendix B Derivative of $\phi_{z}$ definitions

We restate the derivative of $\phi_{z}$ here:

$\begin{split}\frac{d\phi_{z}}{dt}=\frac{1}{\sin^{2}(\theta_{L})}\left[\frac{d\hat{L}}{dt}\cdot\left(\hat{J}\times\hat{L}\right)\right].\end{split}$
(76)

By substituting the derivative $d\hat{L}/dt$, we can separate the behaviour
into the terms that come from the quadrupole moment constant being $\neq 1$
and those that correspond to the black hole portion, doing so gives:

$\begin{split}\frac{d\phi_{z}}{dt}=&-\frac{y^{6}}{2\sin^{2}(\theta_{L})}\left[\hat{L}\times(\mu_{1}\vec{s}_{1}+\mu_{2}\vec{s}_{2})\right]\cdot(\hat{J}\times\hat{L})\\\
&-\frac{3y^{6}}{2\sin^{2}(\theta_{L})}\left[1-y\chi_{\text{eff}}\right]\left\\{\hat{L}\times\vec{s}_{1}+\hat{L}\times\vec{s}_{2}\right\\}\cdot(\hat{J}\times\hat{L})\\\
&+\frac{3y^{7}}{2\sin^{2}(\theta_{L})}\left\\{(\kappa_{1}-1)(\hat{L}\cdot\vec{s}_{1})(\hat{L}\times\vec{s}_{1})+(\kappa_{2}-1)(\hat{L}\cdot\vec{s}_{2})(\hat{L}\times\vec{s}_{2})\right\\}\cdot(\hat{J}\times\hat{L}).\end{split}$
(77)

By relying on the definition of $\vec{J}$ given in Eq. (8), in addition to the
definition of $\theta_{L}$ as the angle between $\vec{L}$ and $\vec{J}$ the
first term can be simplified substantially

$\begin{split}\frac{d\phi_{z}}{dt}=&\frac{Jy^{6}}{2}-\frac{3y^{6}}{2\sin^{2}(\theta_{L})}\left[1-y\chi_{\text{eff}}\right]\left\\{\hat{L}\times\vec{s}_{1}+\hat{L}\times\vec{s}_{2}\right\\}\cdot(\hat{J}\times\hat{L})\\\
&+\frac{3y^{7}}{2\sin^{2}(\theta_{L})}\left\\{(\kappa_{1}-1)(\hat{L}\cdot\vec{s}_{1})(\hat{L}\times\vec{s}_{1})+(\kappa_{2}-1)(\hat{L}\cdot\vec{s}_{2})(\hat{L}\times\vec{s}_{2})\right\\}\cdot(\hat{J}\times\hat{L}).\end{split}$
(78)

The latter terms correspond to neutron stars, while the former terms (alone)
compose the entirety of the black hole case. To simplify the derivative of
$\phi_{z}$ we start by looking at its expression in the black hole case.
Substituting Eq. (23) gives

$\begin{split}\left[\frac{d\phi_{z}}{dt}\right]^{\text{BH}}\approx&\frac{Jy^{6}}{2}\left\\{Q_{1}+Q_{2}\sin(\psi)+\frac{(Q_{3}+Q_{4}\sin(\psi))(Q_{5}+Q_{6}\sin(\psi)+Q_{7}\sin^{2}(\psi))}{Q_{8}Q_{9}(1-\frac{Q_{10}}{Q_{8}}\sin(\psi))(1+\frac{Q_{10}}{Q_{9}}\sin(\psi))}\right\\},\end{split}$
(79)

where

$\displaystyle Q_{1}$
$\displaystyle=1+\frac{3}{2\eta}\left(1-y\left<\chi_{\text{eff}}\right>\right)$
(80a) $\displaystyle Q_{2}$
$\displaystyle=-\frac{3}{2\eta}yG_{\chi_{\text{eff}}}$ (80b) $\displaystyle
Q_{3}$
$\displaystyle=-\frac{3(1+q)}{2q}\left(1-y\left<\chi_{\text{eff}}\right>\right)=\frac{1-Q_{1}}{1+q}$
(80c) $\displaystyle Q_{4}$
$\displaystyle=\frac{3(1+q)}{2q}yG_{\chi_{\text{eff}}}=\frac{-Q_{2}}{1+q}$
(80d) $\displaystyle Q_{5}$
$\displaystyle=4(1-q)(S_{1}^{2}-S_{2}^{2})-(1+q)(\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>)(\delta\mu\left<\delta\chi\right>+(1-4\eta)\left<\chi_{\text{eff}}\right>)$
(80e) $\displaystyle Q_{6}$
$\displaystyle=-(1+q)\left[(\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>)(\delta\mu
G_{\delta\chi}+(1-4\eta)G_{\chi_{\text{eff}}})+(\delta\mu
G_{\delta\chi}+G_{\chi_{\text{eff}}})(\delta\mu\left<\delta\chi\right>+(1-4\eta)\left<\chi_{\text{eff}}\right>)\right]$
(80f) $\displaystyle Q_{7}$ $\displaystyle=-(1+q)(\delta\mu
G_{\delta\chi}+G_{\chi_{\text{eff}}})(\delta\mu
G_{\delta\chi}+(1-4\eta)G_{\chi_{\text{eff}}})$ (80g) $\displaystyle Q_{8}$
$\displaystyle=2\left<J\right>-\delta\mu\left<\delta\chi\right>-\left<\chi_{\text{eff}}\right>-2L$
(80h) $\displaystyle Q_{9}$
$\displaystyle=2\left<J\right>+\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>+2L=4\left<J\right>-Q_{8}$
(80i) $\displaystyle Q_{10}$ $\displaystyle=\delta\mu
G_{\delta\chi}+G_{\chi_{\text{eff}}}.$ (80j)

The rest of the terms, corresponding to the spin induced quadrupole moment of
the neutron star are given by

$\begin{split}\frac{d\phi_{z}}{dt}=\left[\frac{d\phi_{z}}{dt}\right]^{\text{BH}}+\frac{Jy^{6}}{2}\left\\{\frac{Q_{11}+Q_{12}\sin(\psi)+(Q_{13}+Q_{14}\sin(\psi))(Q_{15}+Q_{16}\sin(\psi)+Q_{17}\sin^{2}(\psi))}{Q_{8}Q_{9}(1-\frac{Q_{10}}{Q_{8}}\sin(\psi))(1+\frac{Q_{10}}{Q_{9}}\sin(\psi))}\right\\},\end{split}$
(81)

where

$\displaystyle Q_{11}$
$\displaystyle=-3y\frac{J^{2}-L^{2}}{L\eta}\left[(\kappa_{2}-1)\mu_{1}(\left<\chi_{\text{eff}}\right>-\left<\delta\chi\right>)+(\kappa_{1}-1)\mu_{2}(\left<\chi_{\text{eff}}\right>+\left<\delta\chi\right>)\right]$
$\displaystyle+3y\frac{S_{1}^{2}-S_{2}^{2}}{L\eta}\left[(\kappa_{2}-1)\mu_{1}(\left<\chi_{\text{eff}}\right>-\left<\delta\chi\right>)-(\kappa_{1}-1)\mu_{2}(\left<\chi_{\text{eff}}\right>+\left<\delta\chi\right>)\right]$
(82a) $\displaystyle Q_{12}$
$\displaystyle=-3y\frac{J^{2}-L^{2}}{L\eta}\left[(\kappa_{2}-1)\mu_{1}(G_{\chi_{\text{eff}}}-G_{\delta\chi})+(\kappa_{1}-1)\mu_{2}(G_{\chi_{\text{eff}}}+G_{\delta\chi})\right]$
$\displaystyle+3y\frac{S_{1}^{2}-S_{2}^{2}}{L\eta}\left[(\kappa_{2}-1)\mu_{1}(G_{\chi_{\text{eff}}}-G_{\delta\chi})-(\kappa_{1}-1)\mu_{2}(G_{\chi_{\text{eff}}}+G_{\delta\chi})\right]$
(82b) $\displaystyle Q_{13}$
$\displaystyle=\frac{3}{2}(\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>)$
(82c) $\displaystyle Q_{14}$ $\displaystyle=\frac{3}{2}(\delta\mu
G_{\delta\chi}+G_{\chi_{\text{eff}}})$ (82d) $\displaystyle Q_{15}$
$\displaystyle=4yA_{\chi_{\text{eff}},\delta\chi}\left<\delta\chi\right>^{2}+2(4yA_{\chi_{\text{eff}},\chi_{\text{eff}}}\left<\chi_{\text{eff}}\right>+\mu_{2}\kappa_{1}-\mu_{1}\kappa_{2}+\delta\mu)\left<\delta\chi\right>$
$\displaystyle+4yA_{\chi_{\text{eff}},\delta\chi}\left<\chi_{\text{eff}}\right>^{2}+2(\mu_{2}\kappa_{1}+\mu_{1}\kappa_{2}-1)\left<\chi_{\text{eff}}\right>$
(82e) $\displaystyle Q_{16}$
$\displaystyle=2G_{\delta\chi}(4y(A_{\chi_{\text{eff}},\chi_{\text{eff}}}\left<\chi_{\text{eff}}\right>+A_{\chi_{\text{eff}},\delta\chi}\left<\delta\chi\right>)-\mu_{1}\kappa_{2}+\mu_{2}\kappa_{1}+\delta\mu)$
$\displaystyle+2G_{\chi_{\text{eff}}}(4y(A_{\chi_{\text{eff}},\chi_{\text{eff}}}\left<\chi_{\text{eff}}\right>+A_{\chi_{\text{eff}},\delta\chi}\left<\delta\chi\right>)+\mu_{1}\kappa_{2}+\mu_{2}\kappa_{1}-1)$
(82f) $\displaystyle Q_{17}$
$\displaystyle=4y(A_{\chi_{\text{eff}},\delta\chi}(G_{\delta\chi}^{2}+G_{\chi_{\text{eff}}}^{2})+2A_{\chi_{\text{eff}},\chi_{\text{eff}}}G_{\delta\chi}G_{\chi_{\text{eff}}}).$
(82g)

Using these definitions the entire derivative can be rewritten as

$\begin{split}\frac{d\phi_{z}}{dt}\approx&\frac{Jy^{6}}{2}\left\\{Q_{1}+Q_{2}\sin(\psi)+\frac{H_{0}+H_{1}\sin(\psi)+H_{2}\sin^{2}(\psi)+H_{3}\sin^{3}(\psi)}{(1+H_{-}\sin(\psi))(1+H_{+}\sin(\psi))}\right\\},\end{split}$
(83)

where

$\displaystyle H_{0}$
$\displaystyle=\frac{Q_{3}Q_{5}+Q_{11}+Q_{13}Q_{15}}{Q_{8}Q_{9}}$ (84a)
$\displaystyle H_{1}$
$\displaystyle=\frac{Q_{3}Q_{6}+Q_{5}Q_{4}+Q_{12}+Q_{13}Q_{16}+Q_{14}Q_{15}}{Q_{8}Q_{9}}$
(84b) $\displaystyle H_{2}$
$\displaystyle=\frac{Q_{3}Q_{7}+Q_{4}Q_{6}+Q_{13}Q_{17}+Q_{14}Q_{16}}{Q_{8}Q_{9}}$
(84c) $\displaystyle H_{3}$
$\displaystyle=\frac{Q_{4}Q_{7}+Q_{14}Q_{17}}{Q_{8}Q_{9}}$ (84d)
$\displaystyle H_{-}$ $\displaystyle=-\frac{Q_{10}}{Q_{8}}$ (84e)
$\displaystyle H_{+}$ $\displaystyle=\frac{Q_{10}}{Q_{9}}.$ (84f)

The angles appearing in the solutions for $\phi_{z}$ and $\zeta$ are related
to the above expressions in the following way

$\displaystyle\Phi_{0}$
$\displaystyle=\frac{Jy^{6}}{2}\left\\{Q_{1}+\frac{H_{2}H_{+}H_{-}-H_{3}H_{-}-H_{3}H_{+}}{{H_{+}}^{2}{H_{-}}^{2}}\right\\}$
(85a) $\displaystyle\Phi_{s}$
$\displaystyle=\frac{Jy^{6}}{2}\left\\{Q_{2}+\frac{H_{3}}{H_{+}H_{-}}\right\\}$
(85b) $\displaystyle\Phi_{+}$
$\displaystyle=\frac{Jy^{6}}{2}\left\\{\frac{H_{0}{H_{+}}^{3}-H_{1}{H_{+}}^{2}+H_{2}{H_{+}}-H_{3}}{(H_{+}-H_{-}){H_{+}}^{2}\sqrt{1-{H_{+}}^{2}}}\right\\}$
(85c) $\displaystyle\Phi_{-}$
$\displaystyle=-\frac{Jy^{6}}{2}\left\\{\frac{H_{0}{H_{-}}^{3}-H_{1}{H_{-}}^{2}+H_{2}{H_{-}}-H_{3}}{(H_{+}-H_{-}){H_{-}}^{2}\sqrt{1-{H_{-}}^{2}}}\right\\}$
(85d) $\displaystyle\Theta_{0}$
$\displaystyle=\frac{2L+\delta\mu\left<\delta\chi\right>+\left<\chi_{\text{eff}}\right>}{2J}$
(85e) $\displaystyle\Theta_{s}$ $\displaystyle=\frac{\delta\mu
G_{\delta\chi}+G_{\chi_{\text{eff}}}}{2J}.$ (85f)

## References

  * [1] R. Abbott et al.  GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run.  11 2021.
  * [2] R. Abbott et al.  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.
  * [3] 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(3):031040, 2019.
  * [4] Eanna E. Flanagan and Tanja Hinderer.  Constraining neutron star tidal Love numbers with gravitational wave detectors.  Phys. Rev. D, 77:021502, 2008.
  * [5] B. P. Abbott et al.  GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral.  Phys. Rev. Lett., 119(16):161101, 2017.
  * [6] David Radice, Albino Perego, Francesco Zappa, and Sebastiano Bernuzzi.  GW170817: Joint Constraint on the Neutron Star Equation of State from Multimessenger Observations.  Astrophys. J. Lett., 852(2):L29, 2018.
  * [7] Tanja Hinderer.  Tidal Love numbers of neutron stars.  Astrophys. J., 677:1216–1220, 2008.
  * [8] Jinniu Hu, Shishao Bao, Ying Zhang, Ken’ichiro Nakazato, Kohsuke Sumiyoshi, and Hong Shen.  Effects of symmetry energy on the radius and tidal deformability of neutron stars in the relativistic mean-field model.  PTEP, 2020(4):043D01, 2020.
  * [9] Gonçalo Castro, Leonardo Gualtieri, Andrea Maselli, and Paolo Pani.  Impact and detectability of spin-tidal couplings in neutron star inspirals.  Phys. Rev. D, 106(2):024011, 2022.
  * [10] Huan Yang, William E. East, and Luis Lehner.  Can we distinguish low mass black holes in neutron star binaries?  Astrophys. J., 856(2):110, 2018.  [Erratum: Astrophys.J. 870, 139 (2019)].
  * [11] Zhen Pan, Zhenwei Lyu, and Huan Yang.  Mass-gap extreme mass ratio inspirals.  Phys. Rev. D, 105(8):083005, 2022.
  * [12] Zhen Pan, Zhenwei Lyu, and Huan Yang.  Wet extreme mass ratio inspirals may be more common for spaceborne gravitational wave detection.  Phys. Rev. D, 104(6):063007, 2021.
  * [13] Zhen Pan and Huan Yang.  Formation Rate of Extreme Mass Ratio Inspirals in Active Galactic Nuclei.  Phys. Rev. D, 103(10):103018, 2021.
  * [14] R. Abbott et al.  GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object.  Astrophys. J. Lett., 896(2):L44, 2020.
  * [15] Kyriakos Vattis, Isabelle S. Goldstein, and Savvas M. Koushiappas.  Could the 2.6 $M_{\odot}$ object in GW190814 be a primordial black hole?  Phys. Rev. D, 102(6):061301, 2020.
  * [16] Daniel A. Godzieba, David Radice, and Sebastiano Bernuzzi.  On the maximum mass of neutron stars and GW190814.  Astrophys. J., 908(2):122, 2021.
  * [17] Elias R. Most, L. Jens Papenfort, Lukas R. Weih, and Luciano Rezzolla.  A lower bound on the maximum mass if the secondary in GW190814 was once a rapidly spinning neutron star.  Mon. Not. Roy. Astron. Soc., 499(1):L82–L86, 2020.
  * [18] Raissa F. P. Mendes and Huan Yang.  Tidal deformability of boson stars and dark matter clumps.  Class. Quant. Grav., 34(18):185001, 2017.
  * [19] Vitor Cardoso and Paolo Pani.  Testing the nature of dark compact objects: a status report.  Living Rev. Rel., 22(1):4, 2019.
  * [20] Ian Harry and Tanja Hinderer.  Observing and measuring the neutron-star equation-of-state in spinning binary neutron star systems.  Class. Quant. Grav., 35(14):145010, 2018.
  * [21] Katerina Chatziioannou, Antoine Klein, Nicolás Yunes, and Neil Cornish.  Constructing Gravitational Waves from Generic Spin-Precessing Compact Binary Inspirals.  Phys. Rev. D, 95(10):104004, 2017.
  * [22] Antoine Klein.  EFPE: Efficient fully precessing eccentric gravitational waveforms for binaries with long inspirals.  6 2021.
  * [23] Katerina Chatziioannou, Antoine Klein, Neil Cornish, and Nicolas Yunes.  Analytic Gravitational Waveforms for Generic Precessing Binary Inspirals.  Phys. Rev. Lett., 118(5):051101, 2017.
  * [24] Theocharis A. Apostolatos, Curt Cutler, Gerald J. Sussman, and Kip S. Thorne.  Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries.  Phys. Rev. D, 49:6274–6297, 1994.
  * [25] Michael Kesden, Davide Gerosa, Richard O’Shaughnessy, Emanuele Berti, and Ulrich Sperhake.  Effective potentials and morphological transitions for binary black-hole spin precession.  Phys. Rev. Lett., 114(8):081103, 2015.
  * [26] Davide Gerosa, Michael Kesden, Ulrich Sperhake, Emanuele Berti, and Richard O’Shaughnessy.  Multi-timescale analysis of phase transitions in precessing black-hole binaries.  Phys. Rev. D, 92:064016, 2015.
  * [27] Thibault Damour.  Coalescence of two spinning black holes: an effective one-body approach.  Phys. Rev. D, 64:124013, 2001.
  * [28] Etienne Racine.  Analysis of spin precession in binary black hole systems including quadrupole-monopole interaction.  Phys. Rev. D, 78:044021, 2008.
  * [29] Eric Poisson.  Gravitational waves from inspiraling compact binaries: The Quadrupole moment term.  Phys. Rev. D, 57:5287–5290, 1998.
  * [30] Jan Steinhoff, Tanja Hinderer, Tim Dietrich, and Francois Foucart.  Spin effects on neutron star fundamental-mode dynamical tides: Phenomenology and comparison to numerical simulations.  Phys. Rev. Res., 3(3):033129, 2021.
  * [31] J. Aasi et al.  Advanced LIGO.  Class. Quant. Grav., 32:074001, 2015.
  * [32] N. V. Krishnendu, K. G. Arun, and Chandra Kant Mishra.  Testing the binary black hole nature of a compact binary coalescence.  Phys. Rev. Lett., 119(9):091101, 2017.
  * [33] N. V. Krishnendu, M. Saleem, A. Samajdar, K. G. Arun, W. Del Pozzo, and Chandra Kant Mishra.  Constraints on the binary black hole nature of GW151226 and GW170608 from the measurement of spin-induced quadrupole moments.  Phys. Rev. D, 100(10):104019, 2019.
  * [34] Horng Sheng Chia, Thomas D. P. Edwards, Richard N. George, Aaron Zimmerman, Adam Coogan, Katherine Freese, Cody Messick, and Christian N. Setzer.  Dimensionally Reduced Waveforms for Spin-Induced Quadrupole Searches.  10 2022.
  * [35] Hideyuki Tagoshi, Shuhei Mano, and Eiichi Takasugi.  PostNewtonian expansion of gravitational waves from a particle in circular orbits around a rotating black hole: Effects of black hole absorption.  Prog. Theor. Phys., 98:829–850, 1997.
  * [36] Eric Poisson and Misao Sasaki.  Gravitational radiation from a particle in circular orbit around a black hole. 5: Black hole absorption and tail corrections.  Phys. Rev. D, 51:5753–5767, 1995.
  * [37] Kashif Alvi.  Energy and angular momentum flow into a black hole in a binary.  Phys. Rev. D, 64:104020, 2001.
  * [38] F. Acernese et al.  Advanced Virgo: a second-generation interferometric gravitational wave detector.  Class. Quant. Grav., 32(2):024001, 2015.
  * [39] B. P. Abbott et al.  Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA.  Living Rev. Rel., 21(1):3, 2018.
  * [40] Huan Yang, William E. East, Vasileios Paschalidis, Frans Pretorius, and Raissa F. P. Mendes.  Evolution of Highly Eccentric Binary Neutron Stars Including Tidal Effects.  Phys. Rev. D, 98(4):044007, 2018.
  * [41] Huan Yang.  Inspiralling eccentric binary neutron stars: Orbital motion and tidal resonance.  Phys. Rev. D, 100(6):064023, 2019.