# Spin and quadrupolar effects in the secular evolution of precessing compact
binaries with black hole, neutron star, gravastar, or boson star components

Zoltán Keresztes†, Márton Tápai‡, László Árpád Gergely⋆ Institute of Physics,
University of Szeged, Dóm tér 9, Szeged 6720, Hungary
†E-mail<EMAIL_ADDRESS>‡E-mail:
<EMAIL_ADDRESS>⋆ E-mail<EMAIL_ADDRESS>
###### Abstract

We discuss precessing compact binaries on eccentric orbit with gravastar,
black hole, neutron star, or boson star components. We derive the secular
evolution equations to second post-Newtonian–order accuracy, with leading-
order spin-orbit, spin-spin, and mass quadrupole-monopole contributions
included. The emerging closed system of first-order differential equations
evolves the pairs of polar and azimuthal angles of the spin and orbital
angular momentum vectors together with the periastron angle. In contrast with
the instantaneous dynamics, the secular dynamics is autonomous. The validity
of the latter is confirmed numerically, showing that secular evolutions look
like smoothed-out instantaneous evolutions over timescales where radiation
reaction is negligible. The secular evolution of the spin polar angles and the
difference of their azimuthal angles generates a closed subsystem, which
despite the apparent singularity of spherical polar coordinates, remains well
defined through aligned configurations. We study analytically this system for
the particular cases of one spin dominating over the other and for black hole
- boson star binaries of equal masses. In the first case known large flip-
flops of the smaller spin are reproduced, when the larger spin is almost
coplanar with the orbit. We also find new, quadrupole-induced flip-flops
occurring when the neutron star with dominant spin has a quadrupolar parameter
$w_{1}\approx 3$. Finally we analyze the evolutions of the spin angles
numerically by comparing the cases when the black hole companion is either a
gravastar, another black hole or a boson star with identical mass. We find
that both the amplitude and period of the flip-flop are the largest, when the
companion is a black hole. In the case of a boson star companion the frequency
of the flip-flop increases significantly. Further, while in the case of
gravastars and black holes a swinging-type azimuthal evolution occurs, with
the spins of the components periodically overpassing each other, their
sequence is conserved when the companion is a boson star.

##  I Introduction

In the O1 and O2 runs of the Advanced LIGO Detector, the LIGO Scientific
Collaboration and the Virgo Collaboration have announced a total of 11
gravitational wave detections, ten of them produced by coalescing stellar mass
black holes LIGOcatO12 , while in one case the source was a coalescing neutron
star binary LIGO5NS . A much larger event rate has been seen in the more
sensitive O3 run, with 39 compact binary mergers identified during its first
half O3a , among them other plausible neutron star mergers and neutron star -
black hole mergers. The spin of the merging black holes was identified with
high significance in two cases during the O1–O2 runs LIGOcatO12 ; LIGOspin ,
while in the first half of the O3 run nine events were identified with the
effective spin parameter nonvanishing in the $90\%$ symmetric credible
interval O3a . The accurate description of both the orbital and spin dynamics
of compact binary systems is important for gravitational wave source modeling;
however, the imprint of spin effects also occurs in radio astronomy. From the
analysis of the Very Long Baseline Interferometry (VLBI) radio data of a
binary spanning over 18 years the spin precession of the dominant supermassive
black hole could be identified VLBIKun .

Whenever the black hole spins $\mathbf{S}_{\mathbf{1}}$,
$\mathbf{S}_{\mathbf{2}}$ and the Newtonian orbital angular momentum
$\mathbf{L}_{\mathbf{N}}$ of the binary are not aligned, they undergo
precession barkeroconnell2 ; barkeroconnell1 ; barkeroconnell . The total
angular momentum
$\mathbf{J}=\mathbf{L}+\mathbf{S}_{\mathbf{1}}+\mathbf{S}_{\mathbf{2}}$
composed of the total orbital angular momentum $\mathbf{L}$ and the individual
spins is conserved up to the second post-Newtonian (2PN) order KWW ; KIDDER52
; WILL ; WW , radiative dissipation appearing at 2.5PN orders Peters . Both
conservative and dissipative contributions to dynamics arising from leading
1.5 PN-order spin-orbit (SO) coupling have been thoroughly analyzed RS ; GPV3
, similarly the 2PN spin-spin (SS) GSS ; GSS2 ; KJ  and the 2PN mass
quadrupole - mass monopole (QM) contributions poisson ; GKQM ; Racine ; HLLR ;
VHF ; BFH . The 2PN self-interaction spin effect in the radiative losses,
representing the 2PN correction to the 1.5 PN order Lense-Thirring
approximation, was first identified in Refs. GSS ; GSS2  and explored later to
derive the respective contributions to the accumulated orbital phase MVG . The
dynamics of compact binary systems has been analyzed to 4PN-order accuracy
either in the nonspinning case 4PNb ; 4PNa , or in a conservative approach at
the Hamiltonian level LS4PN . Hamiltonian methods applied for compact binaries
generated several notable results; see Ref. HamPN  and references herein.

Compact binary dynamics in the inspiral regime exhibit three distinct
timescales. The shortest is the radial timescale, defined by the period of the
orbital motion of the reduced mass. Under the precessional timescale the
orbital angular momentum $\mathbf{L}_{\mathbf{N}}$ and the spins
$\mathbf{S}_{\mathbf{1}}$ and $\mathbf{S}_{\mathbf{2}}$ undergo a full
rotation about their precession axis. Over the gravitational radiation
reaction timescale the effects of gravitational dissipation become noticeable.
Averaging the dynamics over any of these timescales may be rewarding. When
precession-related effects are targeted, averaging over a radial period will
remove insignificant instantaneous effects, but keep the dominant precessional
evolution. Several interesting spin-related evolutions were identified by this
method. When the orbital angular momentum nearly cancels the total spin the
orbital plane changes significantly during a relatively short-lived
transitional precession period ACST . The direction of dominant spin relative
to the total angular momentum can change significantly over the radiation
reaction timescale in binaries where the components have significantly
different mass spinflip . This spin flip may explain the formation of X-shape
radio galaxies xshaped . By an additional averaging over the precession
timescale the evolution of the magnitude of total angular momentum over the
radiation reaction timescale was investigated in Refs. kesden1 ; kesden2 .

When the dominant spin vector is approximately perpendicular to
$\mathbf{L}_{\mathbf{N}}$ and the smaller spin is closely aligned with it, the
smaller spin slowly evolves to be antialigned with $\mathbf{L}_{\mathbf{N}}$,
then periodically changes back and forth on a timescale shorter than the
gravitational radiation reaction timescale flipflop . This effect was
investigated in a wider parameter range for binaries moving on circular orbits
flipflop1 ; flipflop2 . This spin flip-flop effect was first found
qualitatively in Ref. majar  as a harmonic "wooble" in the polar angle of the
spin, which evolves "from pole to pole". Recently, a parameter range has been
identified where the flip-flop happens on relatively short timescales, dubbed
as wide precession Gerosa2019 . Then over the period during which the
magnitude of the total spin changes from its minimum to its maximum and back
to the minimum value, one of the two spins evolves from complete alignment
with $\mathbf{L}_{\mathbf{N}}$ to complete antialignment.

Kidder has derived the orbit-averaged (secular) spin-precession equations
KIDDER52  for circular orbits, with SO and SS contributions included, but QM
contributions omitted. The QM couplings were included in the discussion of
angular evolutions in black hole binaries (where the quadrupole parameters are
$w_{i}=1$) by Racine Racine , who presented a new constant of motion $\lambda$
of the orbit-averaged dynamics (or equivalently $\xi=\lambda L_{N},$ the last
factor representing the magnitude of $\mathbf{L}_{\mathbf{N}}$).

Racine also solved analytically the averaged equations for equal masses and
derived approximate analytical solutions in the unequal mass case. His
analytical solution has been generalized for arbitrary masses and spins (but
still $w_{i}=1$ and circular orbits) in Ref. KGSBS , which identified three
distinct regimes in the orbit-averaged precessional dynamics: librations about
the configurations of the two spin projections to the orbital plane either
aligned or antialigned and a "circulating" configuration, when one of the
spins precesses much faster.

Because of the recent discovery of gravitational waves from mergers of neutron
stars LIGO5NS ; 2NS , the interest in their internal structure and equation of
state RPZB , implying a better constrained range of the parameter $w_{i}$ has
been revitalized. In this work we leave the quadrupole parameter unspecified,
including neutrons stars with $w_{i}\in$ $\left(2,14\right)$ NSw , NSUrbanecw
and also other exotic compact objects, like boson stars with
$w_{i}\in\left(10,150\right)$ BosonStarw  or gravastars with
$w_{i}\in\left(-0.8,1\right)$ GravaStarw  as compact binary components.
Allowing for $w_{i}\neq 1$ complicates the dynamics, as for example the
quantity conserved for black hole binaries, identified by Racine, becomes
dynamical, unless $w_{i}=1$.

Although gravitational radiation tends to circularize the orbit of the binary
Peters , significant eccentricity can be preserved at the end of the inspiral.
This happens for binaries in either dense galactic nuclei ecc01 ; ecc02  or
within accretion disks ecc03 ; ecc04 . Furthermore, because of the Kozai
mechanism, the relativistic orbital resonances in hierarchical triples can
also retain eccentricity ecc05 ; ecc06 ; ecc07 ; ecc08 . The interaction
between supermassive black hole binaries and their star populations results in
significant eccentricity toward the end of the inspiral too ecc09 ; ecc10 .
Hence allowing for eccentricity in the dynamics may be rewarding.

The instantaneous dynamics (including SO, SS and QM effects, also
eccentricity) in terms of dimensionless variables was discussed in Refs.
Inspiral1 ; Inspiral2 ; chameleon , based on earlier works on binary dynamics
of Refs. barkeroconnell2 ; barkeroconnell1 ; barkeroconnell ; KWW ;
damourderuelle ; blanchet ; OBKPSSZ ; KMG . The 2PN conservative dynamics of
compact binary systems was given by Eqs. (36)-(42) of Ref. chameleon  in terms
of dimensionless osculating orbital elements $\mathfrak{l}_{r}$, $e_{r}$,
$\psi_{p}$, $\alpha$, and $\phi_{n}$; spin polar and azimuthal angles
$\kappa_{i}$ and $\zeta_{i}$ ($i=1,2$); and true anomaly parameter $\chi_{p}$.
The time evolution of $\chi_{p}$ is governed by Eq. (43) of Ref. chameleon .
The polar angles $\kappa_{i}$ of the spin vectors are measured from
$\mathbf{L}_{\mathbf{N}}$. The azimuthal angles $\zeta_{i}$ are measured from
the Laplace-Runge-Lenz vector $\mathbf{A}_{\mathbf{N}}$ in the plane
perpendicular to $\mathbf{L}_{\mathbf{N}}$. The argument of the periastron,
$\psi_{p}$ is defined by
$\psi_{p}=\arccos\left(\mathbf{\hat{l}}\cdot\mathbf{\hat{A}}_{\mathbf{N}}\right)$,
with $\mathbf{\hat{l}}=\mathbf{\hat{J}}\times\mathbf{\hat{L}}_{\mathbf{N}}$,
where $\mathbf{\hat{J}}$ is the direction of the total angular momentum which
is conserved in the 2PN dynamics. The inclination $\alpha$ is the polar angle
of $\mathbf{\hat{L}}_{\mathbf{N}}$ measured from $\mathbf{\hat{J}}$. The last
angle is the longitude of the ascending node $-\phi_{n}$, spanned by the
inertial axis $\mathbf{\hat{x}}$ (arbitrarily chosen in the plane
perpendicular to $\mathbf{\hat{J}}$) and $\mathbf{\hat{l}}$. This angle is
related to the azimuthal angle $\left(\pi/2-\phi_{n}\right)$ of
$\mathbf{\hat{L}}_{\mathbf{N}}$, measured from $\mathbf{\hat{x}}$ in the plane
perpendicular to $\mathbf{\hat{J}}$. (The angles $\chi_{p}$, $\kappa_{i}$,
$\psi_{i}=\zeta_{i}-\psi_{p}$, $\psi_{p}$, $\alpha$, and $\phi_{n}$ are shown
in Fig. 2 in Ref. Inspiral1 .)

The dimensionless orbital angular momentum

$\mathfrak{l}_{r}=\frac{cL_{N}}{Gm\mu}\leavevmode\nobreak\ $ (1)

and the eccentricity

$e_{r}=\frac{A_{N}}{Gm\mu}\leavevmode\nobreak\ $ (2)

characterize the osculating ellipse of the orbit; hence, they are shape
variables. The total and reduced masses of the binary are denoted as
$m=m_{1}+m_{2}$ and $\mu=m_{1}m_{2}/m$, respectively. We also employ the mass
ratio $\nu=m_{2}/m_{1}\leq 1$ and the symmetric mass ratio $\eta=\mu/m$. The
gravitational constant, the speed of light, and the magnitude of
$\mathbf{A}_{\mathbf{N}}$ are denoted by $G$, $c$, and $A_{N}$ respectively.
The magnitudes of the spins are characterized by the dimensionless spin
parameters $\chi_{i}=cS_{i}/Gm_{i}^{2}$ $\left(i=1,2\right)$. The dot will
denote the derivative with respect to the dimensionless time variable
$\mathfrak{t}=tc^{3}/Gm$ (with time $t$) introduced in Ref. chameleon . The PN
parameter is defined as $\varepsilon=Gm/c^{2}r$.

In this paper, we will investigate precessing compact binary systems on
eccentric orbit subject to bound motion, first establishing the 2PN secular
dynamics in terms of the above-mentioned dimensionless variables, then
analyzing the spin evolutions with the methods of dynamical systems, with SO,
SS and QM contributions included. The PN equations of motion depend on the
choice of a spin supplementary condition (SSC)111A comparison of the three
widely used SSCs has been presented in Ref. mikoczissc , proving the SSC
dependence of the radiative multipole moments.. We employ the Newton-Wigner-
Price NWP1 ; NWP2  SSC, similarly as in Ref. chameleon . This system is valid
for eccentric orbits and for binaries composed of either black holes, neutron
stars or other exotic compact objects (boson stars or gravastars).

The instantaneous dynamics described in Ref. chameleon  is averaged over a
suitably defined radial period in Section II, obtaining the secular precessing
compact binary dynamics in terms of the dimensionless osculating orbital
elements $\mathfrak{l}_{r}$, $e_{r}$, $\psi_{p}$, $\alpha$, $\phi_{n}$ and
spin angles $\kappa_{i}$, $\zeta_{i}$ ($i=1,2$). These equations contain PN,
SO, SS, QM and 2PN contributions. These are generalized Lagrangian planetary
equations, which become singular for vanishing $e_{r}$, nevertheless the
singularity can be eliminated by a transformation of variables
LagrangePlanetarytrafo . For completeness, we also give in this section the
secular precession angular velocities and constraints relating the variables.

For the purpose of averaging, the PN expansion of the radial period is
required. For clarity of presentation we deferred the tedious but
straightforward bulk of calculations leading to it to the lengthy Appendix A.
There, first, we derive the radial period in terms of the variables evaluated
at the periastron (characterized by the value of the true anomaly parameter
$\chi_{p}=0$). The $\chi_{p}$ dependence of the shape variables is also
derived by integrating the corresponding first-order system given in Ref.
chameleon . Then, we compute their time-averaged values
$\bar{\mathfrak{l}}_{r}$, $\bar{e}_{r}$ over the radial period to 2PN
accuracy, with the inclusion of all spin and mass quadrupole effects to this
order. Next, we express the shape variables evaluated in the periastron in
terms of the corresponding averaged quantities. With this, we can write the
radial period as a PN expansion expressed in terms of averaged quantities. We
also give there a similar expansion in terms of the chosen variables of the PN
parameter.

In Sec. III, we analyze the role of eccentricity in the secular evolution by
comparing low-eccentricity and medium-eccentricity evolutions for the values
of the PN parameter $0.01$ and $0.0005$. We also prove that the secular
dynamics follows closely the instantaneous dynamics given in Ref. chameleon
over the conservative timescale.

The secular evolution of the spin angles generates a closed subsystem of three
variables, discussed in Sec. IV. Despite the apparent singularity of spherical
polar coordinates, we prove in Appendix B that the system remains well defined
through aligned configurations.

In Sec. V, we study analytically and numerically this closed subsystem for the
particular case of one spin dominating over the other, concentrating on the
flip-flop effect of the polar angle of the smaller spin. We identify a
diamond-shaped region in the parameter plane span by the dominant spin polar
angle and quadrupolar parameter. Along the horizontal axis of the diamond the
known flip-flop effect for black hole binaries is reproduced; however, the
vertical axis signifies mass quadrupole-induced flip-flops occurring for
neutron stars with a particular quadrupolar parameter.

In Sec. VI, we study the closed subsystem for the particular case of black
hole - boson star binaries. We investigate the spin angle dynamics both
analytically and numerically. We also compare it to typical evolutions in
black hole - black hole and black hole - gravastar binaries, pointing out the
differences.

In Sec. VII, we present the conclusions.

##  II Secular conservative dynamics of precessing compact binaries

In this section we present the orbital averaged evolution equations of the
dimensionless osculating orbital elements $\mathfrak{l}_{r}$, $e_{r}$,
$\psi_{p}$, $\alpha$, and $\phi_{n}$ and spin angles $\kappa_{i}$, $\zeta_{i}$
($i=1,2$) at 2PN accuracy, with spin and mass quadrupole effects included from
the instantaneous evolutions derived in Ref. chameleon .

###  II.1 Averaging method

For bounded motion the separation $r$ between the binary components can be
parametrized similarly to a Keplerian orbit chameleon :

$\frac{c^{2}r}{Gm}=\frac{\mathfrak{l}_{r}^{2}}{1+e_{r}\cos\chi_{p}}\leavevmode\nobreak\
,$ (3)

where $\chi_{p}$ is the true anomaly. However, unlike a Keplerian orbit, both
shape variables $\mathfrak{l}_{r}$ and $e_{r}$ are time dependent.

The dimensionless period $\mathfrak{T}$ is defined as the change in the
dimensionless time $\mathfrak{t}\equiv tc^{3}/Gm$ during the evolution of the
true anomaly over $\chi_{p}\in\left[0,2\pi\right]$ as

$\mathfrak{T}\equiv\int_{0}^{\mathfrak{T}}d\mathfrak{t}=\int_{0}^{2\pi}\frac{d\chi_{p}}{\dot{\chi}_{p}}\leavevmode\nobreak\
.$ (4)

The average $\bar{f}$ of any quantity $f\left(\mathfrak{t}\right)$ with
respect to $\mathfrak{t}$ is introduced by

$\mathfrak{T}\bar{f}=\int_{0}^{\mathfrak{T}}f\left(\mathfrak{t}\right)d\mathfrak{t}=\int_{0}^{2\pi}\frac{f\left(\chi_{p}\right)}{\dot{\chi}_{p}}d\chi_{p}\leavevmode\nobreak\
.$ (5)

Then, as described in Appendix A the period can be rewritten in terms of the
orbital averaged shape variables as

$\displaystyle\mathfrak{T}$ $\displaystyle=$
$\displaystyle\mathfrak{\tilde{T}}\left(1+\frac{1}{\mathfrak{\bar{l}}_{r}^{2}}\tilde{\tau}_{PN}+\frac{1}{\mathfrak{\bar{l}}_{r}^{3}}\tilde{\tau}_{SO}+\frac{1}{\mathfrak{\bar{l}}_{r}^{4}}\tilde{\tau}_{QM}\right.$
(6)
$\displaystyle\left.+\frac{1}{\mathfrak{\bar{l}}_{r}^{4}}\tilde{\tau}_{SS}+\frac{1}{\mathfrak{\bar{l}}_{r}^{4}}\tilde{\tau}_{2PN}\right)\leavevmode\nobreak\
,$

with the expressions $\mathfrak{\tilde{T}}$, $\tilde{\tau}_{PN}$,
$\tilde{\tau}_{SO}$, $\tilde{\tau}_{QM}$,$\leavevmode\nobreak\
\tilde{\tau}_{SS}$, and$\leavevmode\nobreak\ \tilde{\tau}_{2PN}$ also given
there. Here, $1/\mathfrak{l}_{r}^{2}$ stands for one relative PN order, as
indicated by Eq. (8) of Ref. chameleon .

###  II.2 Shape variables

Employing the averaging method described in the previous section, long but
straightforward calculations lead to the secular evolutions of the
dimensionless orbital angular momentum $\mathfrak{l}_{r}$ and orbital
eccentricity $e_{r}$ as

$\mathfrak{\bar{\dot{l}}}_{r}=\mathfrak{\bar{\dot{l}}}_{r}^{PN}=\mathfrak{\bar{\dot{l}}}_{r}^{SO}=\mathfrak{\bar{\dot{l}}}_{r}^{SS}=\mathfrak{\bar{\dot{l}}}_{r}^{QM}=\mathfrak{\bar{\dot{l}}}_{r}^{2PN}=0\leavevmode\nobreak\
,$ (7)
$\bar{\dot{e}}_{r}=\bar{\dot{e}}_{r}^{PN}=\bar{\dot{e}}_{r}^{SO}=\bar{\dot{e}}_{r}^{SS}=\bar{\dot{e}}_{r}^{QM}=\bar{\dot{e}}_{r}^{2PN}=0\leavevmode\nobreak\
.$ (8)

As expected, the average shape of the orbit does not change when dissipation
by gravitational waves is neglected.

###  II.3 Euler angles

These evolutions are nontrivial, as discussed below.

####  II.3.1 Inclination $\alpha$

The secular evolution of the inclination
$\alpha=\arccos\left(\mathbf{\hat{J}}\cdot\mathbf{\hat{L}}_{\mathbf{N}}\right)$
emerges as the sum of the contributions:

$\bar{\dot{\alpha}}^{PN}=0\leavevmode\nobreak\ ,$ (9)
$\displaystyle\bar{\dot{\alpha}}^{SO}$ $\displaystyle=$
$\displaystyle\frac{\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}\sum_{k=1}^{2}\left(4\nu^{2k-3}+3\right)$ (10)
$\displaystyle\times\chi_{k}\sin\kappa_{k}\cos\left(\psi_{p}+\zeta_{k}\right)\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\alpha}}^{SS}$ $\displaystyle=$
$\displaystyle-\frac{3\eta\pi}{\mathfrak{T\leavevmode\nobreak\
\bar{l}}_{r}^{4}\leavevmode\nobreak\ }\chi_{1}\chi_{2}$ (11)
$\displaystyle\times\left[\sin\kappa_{1}\cos\kappa_{2}\cos\left(\psi_{p}+\zeta_{1}\right)\right.$
$\displaystyle\left.+\cos\kappa_{1}\sin\kappa_{2}\cos\left(\psi_{p}+\zeta_{2}\right)\right]\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\alpha}}^{QM}$ $\displaystyle=$
$\displaystyle-\frac{3\eta\pi}{2\mathfrak{T\leavevmode\nobreak\
\bar{l}}_{r}^{4}\leavevmode\nobreak\
}\sum_{k=1}^{2}\nu^{2k-3}w_{k}\chi_{k}^{2}$ (12) $\displaystyle\times\sin
2\kappa_{k}\cos\left(\psi_{p}+\zeta_{k}\right)\leavevmode\nobreak\ .$
$\bar{\dot{\alpha}}^{2PN}=0\leavevmode\nobreak\ .$ (13)

As expected, the inclination only changes due to spin and quadrupolar effects.

####  II.3.2 Longitude of the ascending node $-\phi_{n}$

The longitude of the ascending node $-\phi_{n}$ is subtended by the inertial
axis $\mathbf{\hat{x}}$ and the ascending node
$\mathbf{\hat{l}}=\mathbf{\hat{L}}_{\mathbf{N}}\times\mathbf{\hat{J}}$. It has
the following contributions to its secular evolution:

$\bar{\dot{\phi}}_{n}^{PN}=0\leavevmode\nobreak\ ,$ (14)
$\displaystyle\bar{\dot{\phi}}_{n}^{SO}$ $\displaystyle=$
$\displaystyle-\frac{\eta\pi}{\mathfrak{T\leavevmode\nobreak\
\bar{l}}_{r}^{3}\sin\alpha}\sum_{k=1}^{2}\left(4\nu^{2k-3}+3\right)$ (15)
$\displaystyle\times\chi_{k}\sin\kappa_{k}\sin\left(\psi_{p}+\zeta_{k}\right)\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\phi}}_{n}^{SS}$ $\displaystyle=$
$\displaystyle\frac{3\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}\sin\alpha}\chi_{1}\chi_{2}$ (16)
$\displaystyle\times\left[\sin\kappa_{1}\cos\kappa_{2}\sin\left(\psi_{p}+\zeta_{1}\right)\right.$
$\displaystyle\left.+\cos\kappa_{1}\sin\kappa_{2}\sin\left(\psi_{p}+\zeta_{2}\right)\right]\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\phi}}_{n}^{QM}$ $\displaystyle=$
$\displaystyle\frac{3\eta\pi}{2\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}\sin\alpha}\sum_{k=1}^{2}\nu^{2k-3}w_{k}\chi_{k}^{2}$
(17) $\displaystyle\times\sin
2\kappa_{k}\sin\left(\psi_{p}+\zeta_{k}\right)\leavevmode\nobreak\ ,$
$\bar{\dot{\phi}}_{n}^{2PN}=0\leavevmode\nobreak\ .$ (18)

Again, only spin and quadrupolar effects contribute.

####  II.3.3 Argument of the periastron $\psi_{p}$

The secular evolution of $\psi_{p}$, the angle between the node line
($\mathbf{\hat{l}}$, perpendicular to both $\mathbf{L}_{\mathbf{N}}$ and
$\mathbf{J}$) and the periastron ($\mathbf{\hat{A}}_{\mathbf{N}}$), is the sum
of:

$\bar{\dot{\psi_{p}}}^{PN}=\frac{6\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{2}}\leavevmode\nobreak\ ,$ (19)
$\displaystyle\bar{\dot{\psi_{p}}}^{SO}$ $\displaystyle=$
$\displaystyle-\frac{\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}\sum_{k=1}^{2}\left(4\nu^{2k-3}+3\right)$ (20)
$\displaystyle\times\chi_{k}\left[2\cos\kappa_{k}\right.$
$\displaystyle\left.+\cot\alpha\sin\kappa_{k}\sin\left(\psi_{p}+\zeta_{k}\right)\right]\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\psi_{p}}}^{SS}$ $\displaystyle=$
$\displaystyle\frac{3\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}}\chi_{1}\chi_{2}$ (21)
$\displaystyle\times\left\\{\cot\alpha\left[\sin\kappa_{1}\cos\kappa_{2}\sin\left(\psi_{p}+\zeta_{1}\right)\right.\right.$
$\displaystyle\left.+\cos\kappa_{1}\sin\kappa_{2}\sin\left(\psi_{p}+\zeta_{2}\right)\right]+2\cos\kappa_{1}$
$\displaystyle\left.\times\cos\kappa_{2}-\sin\kappa_{1}\sin\kappa_{2}\cos\,\left(\zeta_{2}-\zeta_{1}\right)\right\\}\leavevmode\nobreak\
,$

$\displaystyle\bar{\dot{\psi_{p}}}^{QM}$ $\displaystyle=$
$\displaystyle\frac{3\eta\pi}{2\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}}\sum_{k=1}^{2}\nu^{2k-3}w_{k}\chi_{k}^{2}$ (22)
$\displaystyle\times\left[\cot\alpha\sin
2\kappa_{k}\sin\left(\psi_{p}+\zeta_{k}\right)\right.$
$\displaystyle\left.-3\sin^{2}\kappa_{k}+2\right]\leavevmode\nobreak\ ,$
$\bar{\dot{\psi_{p}}}^{2PN}=\frac{3\pi}{2\mathfrak{T\bar{l}}_{r}^{4}}\left[\allowbreak
33\bar{e}_{r}^{2}-4\eta-6\bar{e}_{r}^{2}\eta+2\right]\leavevmode\nobreak\ .$
(23)

All PN, spin and quadrupolar corrections lead to periastron precession.

###  II.4 Spin angles

The contributions to the secular evolutions of the spin polar angles
$\kappa_{i}$ and the azimuthal angles $\zeta_{i}$ are

$\bar{\dot{\kappa_{i}}}^{PN}=0\leavevmode\nobreak\ ,$ (24)
$\displaystyle\bar{\dot{\kappa_{i}}}^{SO}$ $\displaystyle=$
$\displaystyle\frac{\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}$ (25)
$\displaystyle\times\left(4\nu^{2j-3}+3\right)\chi_{j}\sin\kappa_{j}\sin\left(\zeta_{i}-\zeta_{j}\right)\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\kappa_{i}}}^{SS}$ $\displaystyle=$
$\displaystyle-\frac{\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}}\chi_{j}\sin\kappa_{j}$ (26)
$\displaystyle\times\sin\left(\zeta_{i}-\zeta_{j}\right)\left(\mathfrak{\bar{l}}_{r}\nu^{2j-3}+3\chi_{i}\cos\kappa_{i}\right)\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\kappa_{i}}}^{QM}$ $\displaystyle=$
$\displaystyle-\frac{3\eta\pi}{2\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}}\nu^{2j-3}w_{j}\chi_{j}^{2}$ (27)
$\displaystyle\times\sin
2\kappa_{j}\sin\left(\zeta_{i}-\zeta_{j}\right)\leavevmode\nobreak\ ,$
$\bar{\dot{\kappa_{i}}}^{2PN}=0\leavevmode\nobreak\ ,$ (28)
$\bar{\dot{\zeta_{i}}}^{PN}=-\bar{\dot{\psi_{p}}}^{PN}\leavevmode\nobreak\ ,$
(29) $\displaystyle\bar{\dot{\zeta_{i}}}^{SO}$ $\displaystyle=$
$\displaystyle\frac{\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}\left\\{\mathfrak{\bar{l}}_{r}\left(4+3\nu^{3-2i}\right)\right.$
(30)
$\displaystyle+3\left(4\nu^{2i-3}+3\right)\chi_{i}\cos\kappa_{i}+\left(4\nu^{2j-3}+3\right)\chi_{j}$
$\displaystyle\times\left.\left[2\cos\kappa_{j}+\cot\kappa_{i}\sin\kappa_{j}\cos\left(\zeta_{i}-\zeta_{j}\right)\right]\right\\}\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\zeta_{i}}}^{SS}$ $\displaystyle=$
$\displaystyle-\frac{\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}\nu^{2j-3}\chi_{j}\left[2\cos\kappa_{j}\right.$
(31)
$\displaystyle\left.+\cot\kappa_{i}\sin\kappa_{j}\cos\left(\zeta_{i}-\zeta_{j}\right)\right]$
$\displaystyle-\frac{3\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}}\chi_{i}\chi_{j}$
$\displaystyle\times\left\\{\cot\kappa_{i}\left[3\sin\kappa_{i}\cos\kappa_{j}\right.\right.$
$\displaystyle\left.+\cos\kappa_{i}\sin\kappa_{j}\cos\left(\zeta_{i}-\zeta_{j}\right)\right]$
$\displaystyle\left.-\sin\kappa_{i}\sin\kappa_{j}\cos\left(\zeta_{i}-\zeta_{j}\right)\right\\}\leavevmode\nobreak\
,$ $\displaystyle\bar{\dot{\zeta_{i}}}^{QM}$ $\displaystyle=$
$\displaystyle-\frac{3\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}w_{i}\chi_{i}\cos\kappa_{i}$ (32)
$\displaystyle-\frac{3\eta\pi}{2\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{4}}\sum_{k=1}^{2}w_{k}\nu^{2k-3}\chi_{k}^{2}$
$\displaystyle\times\left[2-3\sin^{2}\kappa_{k}\right.$
$\displaystyle\left.+\cot\kappa_{i}\sin\left(2\kappa_{k}\right)\cos\left(\zeta_{i}-\zeta_{k}\right)\right]\leavevmode\nobreak\
,$
$\bar{\dot{\zeta_{i}}}^{2PN}=-\bar{\dot{\psi_{p}}}^{2PN}\leavevmode\nobreak\
.$ (33)

Here, $i\neq jm,$ and $i=1$, $2$.

It is easy to check that in the binary black hole case the dynamics indeed
provides the additional constant of motion found by Racine Racine , which in
our notation becomes

$\xi=\frac{G\mu
m}{c}\left[\left(1+\nu^{-1}\right)\chi_{1}\cos\kappa_{1}+\left(1+\nu\right)\chi_{2}\cos\kappa_{2}\right]$
(34)

Its time derivative indeed vanishes as can be seen by employing Eqs. (24)–(28)
with $w_{i}=1$. However, there is no obvious way to generalize this constant
for arbitrary quadrupole parameter $w_{i}\neq 1$.

###  II.5 Secular precession angular velocities

The averaged precession angular velocities are calculated from Eqs. (31)-(33)
of Ref. chameleon 222In Eqs. (B34) of Ref. Inspiral2  the SS and QM terms have
typos: the $1/2$ factors should be removed. We thank Krisztina Kövér for
pointing this out. Because of this, the second term of Eq. (33) of Ref.
chameleon  contains unnecessary $1/2$ factors on the rhs (but otherwise all
conclusions remain unchanged). In the present paper, these have been
corrected, and both the instantaneous and secular dynamics are represented
correctly..

$\displaystyle\overline{\mathbf{\omega}_{\mathbf{i}}\cdot\mathbf{\hat{A}}_{\mathbf{N}}}$
$\displaystyle=$
$\displaystyle\frac{2\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}\left(\nu^{2j-3}\chi_{j}\sin\kappa_{j}\cos\zeta_{j}\right.$
(35) $\displaystyle\left.+3w_{i}\chi_{i}\sin\kappa_{i}\cos\zeta_{i}\right)$
$\displaystyle\overline{\mathbf{\omega}_{\mathbf{i}}\cdot\mathbf{\hat{Q}}_{\mathbf{N}}}$
$\displaystyle=$
$\displaystyle\frac{2\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}\left(\nu^{2j-3}\chi_{j}\sin\kappa_{j}\sin\zeta_{j}\right.$
(36) $\displaystyle\left.+3w_{i}\chi_{i}\sin\kappa_{i}\sin\zeta_{i}\right)$
$\displaystyle\overline{\mathbf{\omega}_{\mathbf{i}}\cdot\mathbf{\hat{L}}_{\mathbf{N}}}$
$\displaystyle=$ $\displaystyle\frac{\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{3}}\left[\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}\left(4+3\nu^{3-2i}\right)\right.$ (37)
$\displaystyle\left.+2\nu^{2j-3}\chi_{j}\cos\kappa_{j}\right]$

It is easy to see by checking the leading-order term of $\mathfrak{T}$ that as
$\bar{e}_{r}$ goes to $1$ the precession becomes increasingly small. It is
explained by the fact that on parabolic orbits, when $\bar{e}_{r}=1$, the
motion becomes unbound, and there is no well-defined period, thus no
precession.

###  II.6 Constraints

Figure 1: The panels represent the secular (green) and instantaneous (red)
evolutions of $\mathfrak{l}_{r}$, $e_{r}$, $\alpha$, $\phi_{n}$, $\psi_{p}$,
$\kappa_{i}$ and $\zeta_{i}$ ($i=1,2$) as functions of the number of orbital
periods $N$. The parameters and initial conditions for the instantaneous
evolutions are: total mass $m=20M_{\odot}$, mass ratio $\nu=0.5$,
dimensionless spin parameters $\chi_{1}=\chi_{2}=0.9982$, eccentricity
$e_{r}\left(0\right)=0.1$, spin polar angles $\kappa_{1}\left(0\right)=\pi/6$
and $\kappa_{2}\left(0\right)=\pi/4$, spin azimuthal angles
$\zeta_{1}\left(0\right)=0$ and $\zeta_{2}\left(0\right)=\pi/2$, longitude of
the ascending node $\phi_{n}\left(0\right)=0$ and post-Newtonian parameter
$\varepsilon\left(0\right)=0.01$. Matching initial conditions for the secular
dynamics have been chosen, as explained in the main text. On some of the plots
the two evolutions overlap at the chosen resolution. On each panel the first
two periods are also shown in order to illustrate that the match only occurs
on the longer run.  Figure 2: The effect of an increased eccentricity
$e_{r}\left(0\right)=0.5$, the rest of the parameters being identical to those
for Fig. 1. Figure 3: The evolutions at a much earlier stage of the inspiral,
characterized by post-Newtonian parameter $\varepsilon\left(0\right)=0.0005$,
the rest of the parameters being identical to those for Fig. 1.  Figure 4: The
effect of an increased eccentricity $e_{r}\left(0\right)=0.5$, the rest of the
parameters being identical to those for Fig. 3.

The 2PN dynamics of compact binaries (given by Eqs. (36)-(45) of Ref.
chameleon ) have four first integrals which are derived in Secs. V.A. and V.B.
of Ref. chameleon . The first integrals give the conservation of total energy
and total angular momentum. Two first integrals which express the direction of
total angular momentum and given by Eq. (53) and the ratio of Eqs. (54) and
(55) of Ref. chameleon  represent two constraints among the variables
occurring in the secular evolution equations. As a result of this, the secular
dynamics given above is subject to two constraints among the shape variables,
the Euler angles, and spin angles, which will be shown below.

The constraints determining the direction of the total angular momentum to
leading order (see Eqs. (58) and (62) of Ref. chameleon ) read as

$\mathfrak{\bar{l}}_{r}=\sum_{i=1}^{2}\nu^{2i-3}\chi_{i}\left[\sin\kappa_{i}\sin\left(\zeta_{i}+\psi_{p}\right)\cot\alpha-\cos\kappa_{i}\right]\
,$ (38)

and

$\sum_{i=1}^{2}\nu^{2i-3}\chi_{i}\sin\kappa_{i}\cos\left(\zeta_{i}+\psi_{p}\right)=0\
.$ (39)

Here $\mathfrak{l}_{r}\left(\chi_{p}\right)$ was changed to
$\mathfrak{\bar{l}}_{r}$, an approximation valid to leading order. These
expressions do not depend explicitly on the variable $\chi_{p}$ characterizing
the position of the body with reduced mass. Only the correction terms
$\mathcal{O}\left(\mathfrak{\bar{l}}_{r}^{-2}\right)$ to the above equations
exhibit explicit $\chi_{p}$ dependence. Therefore, we expect the constraints
(38) and (39) to be conserved under the secular evolution.

Since $\mathfrak{\bar{l}}_{r}$ is a constant, the time derivative of the left-
hand side of (38) is zero. A long but straightforward computation with the
application of the secular evolution equations shows that the time derivative
of the right-hand side also vanishes. Similarly we checked that the time
derivative of Eq. (39) vanishes to 2PN order. These prove that Eqs. (38) and
(39) are first integrals of the secular evolution equations given in Secs.
II.2–II.4. They can be used to reduce the order of the differential equation
system or for checking the accuracy of the numerical integration. In addition,
Eqs. (38) and (39) have to be fulfilled when setting up initial conditions.

##  III Secular and instantaneous dynamics of black hole binaries compared

###  III.1 Conservative timescale

Because of gravitational radiation, the PN parameter increases during the
inspiral. During the time

$\Delta\tau=\frac{5}{256\eta}\left(\varepsilon_{\left(in\right)}^{-4}-\varepsilon_{\left(out\right)}^{-4}\right)\frac{Gm}{c^{3}}\leavevmode\nobreak\
$ (40)

the PN parameter evolves from $\varepsilon_{\left(in\right)}$ to
$\varepsilon_{\left(out\right)}$ MVG . For numerical evolution purposes, we
define the dimensionless conservative timescale $\mathfrak{\tilde{T}}_{cons}$
as a 1% increase in the PN parameter:

$\mathfrak{\tilde{T}}_{cons}=\frac{5}{256\eta}\varepsilon_{\left(in\right)}^{-4}\left(1-\frac{1}{1.01^{4}}\right)\approx\frac{0.195}{256\eta}\varepsilon_{\left(in\right)}^{-4}\leavevmode\nobreak\
.$ (41)

Dividing this by (the leading-order dimensionless) orbital period
$\mathfrak{\tilde{T}}_{orb}=2\pi\varepsilon^{-3/2}$ gives the number

$N_{\max}=\frac{0.195}{512\pi\eta}\varepsilon_{\left(in\right)}^{-5/2}\leavevmode\nobreak\
$ (42)

of radial periods of the conservative timescale. Evolving numerically for
$N_{\max}$ periods keeps the error of disregarding the dissipation below
$1\%$.

###  III.2 Accuracy of the secular dynamics

We check the long-term accuracy of the secular dynamics by a numerical
comparison over $N_{\max}$ periods with the instantaneous evolutions given by
Eqs. (36)-(42) of Ref. chameleon . The results are shown on Figs. 1–4. All
figures are for $m=20M_{\odot}$, $\nu=1/2$, and high spin parameter values
$\chi_{1}=\chi_{2}=0.9982$. The initial value for the PN parameter is
$\varepsilon\left(0\right)=10^{-2}$ in Figs. 1 and 2, while
$\varepsilon\left(0\right)=5\times 10^{-4}$ in Figs. 3 and 4. The initial
value for the eccentricity is $e_{r}\left(0\right)=0.1$ in Figs. 1 and 3 and
$e_{r}\left(0\right)=0.5$ in Figs. 2 and 4.

For the instantaneous evolutions, the shape parameter
$\mathfrak{l}_{r}\left(0\right)$ emerged from the PN parameter
$\varepsilon\left(0\right)$ and eccentricity $e_{r}\left(0\right)$ cf. Eq. (8)
of Ref. chameleon . Then $\psi_{p}\left(0\right)$ and $\alpha\left(0\right)$
were computed from the constraints given by Eq. (53) and the ratio of Eqs.
(54) and (55) of Ref. chameleon .

The initial values $\bar{f}\left(0\right)$ for the secular dynamics were
extracted from the instantaneous evolution during the first orbit in the
following manner:

$\bar{f}\left(0\right)=\bar{f}-\frac{1}{2}\left[f\left(\mathfrak{T}\right)-f\left(0\right)\right]\leavevmode\nobreak\
,$ (43)

where $\bar{f}$ represents the first orbital average. This method corrects for
the periodic component of the instantaneous motion leading to the proper
initial condition for the secular evolution representing the orbital average.

In the secular dynamics the shape variables $\bar{e}_{r}$ and
$\mathfrak{\bar{l}}_{r}$ are conserved. Finally, the initial values of
$\psi_{p}$ and $\alpha$ were computed from the constraints (38) and (39).

The smaller pictures zoom into the first two radial periods of the evolution.
On the conservative timescale the secular evolution follows closely the
instantaneous evolutions. Certain evolutions even overlap making them hard to
distinguish.

###  III.3 Effects of the eccentricity and PN parameter

The comparison of Figs. (1) and (2), and also (3) and (4), yields to the
remark that $N_{\max}$ decreases with increasing eccentricity, whereas
comparing the evolutions of distant binaries in Figs. (2) and (4) to the close
binary evolutions (1) and (3) yields that $N_{\max}$ decreases with increasing
PN parameter.

While Figs. (3) and (4) represent the early stages of the inspiral, Figs. (1)
and (2) refer to the very end of it. While the secular evolution of the spin
azimuthal angles $\zeta_{1,2}$ indicate more than $60$ precessional cycles at
low eccentricity and more than $20$ for high eccentricity in the early
inspiral regime, at the end of it, $N_{\max}$ will not reach even one single
precessional cycle. Similarly, the number of nutations (represented by the
number of cycles of the polar angles $\kappa_{i}$ between their maximal and a
minimal values) is $18$ at low eccentricity and $5$ at high eccentricity in
the early inspiral, while at the end of it, $N_{\max}$ covers only a fraction
of the nutational period. The polar angle $\alpha$ and azimuthal angle
$\phi_{n}$ of the orbital angular momentum exhibit similar behavior to that of
the spins.

##  IV Closed system for the secular spin angle evolutions

Remarkably, the secular evolution of the spin angles $\kappa_{1}$ and
$\kappa_{2}$ and $\Delta\zeta\equiv\zeta_{1}-\zeta_{2}$ discussed in Sec. II
form a closed subset:

$\displaystyle\frac{1}{R}\overline{\frac{d\kappa_{1}}{d\mathfrak{t}}}$
$\displaystyle=$ $\displaystyle\left(1+\nu-x_{1}\cos\kappa_{1}-\nu
w_{2}x_{2}\cos\kappa_{2}\right)$ (44) $\displaystyle\times
x_{2}\sin\kappa_{2}\sin\Delta\zeta\leavevmode\nobreak\ ,$
$\displaystyle\frac{1}{R}\overline{\frac{d\kappa_{2}}{d\mathfrak{t}}}$
$\displaystyle=$
$\displaystyle-\left(1+\nu^{-1}-x_{2}\cos\kappa_{2}-\nu^{-1}w_{1}x_{1}\cos\kappa_{1}\right)$
(45) $\displaystyle\times
x_{1}\sin\kappa_{1}\sin\Delta\zeta\leavevmode\nobreak\ ,$
$\displaystyle\frac{1}{R}\overline{\frac{d\Delta\zeta}{d\mathfrak{t}}}$
$\displaystyle=$ $\displaystyle\nu-\nu^{-1}+\left(1+2\nu^{-1}-w_{1}\right.$
(46)
$\displaystyle\left.-w_{1}\nu^{-1}x_{1}\cos\kappa_{1}\right)x_{1}\cos\kappa_{1}$
$\displaystyle-\left(1+2\nu-w_{2}-w_{2}\nu
x_{2}\cos\kappa_{2}\right)x_{2}\cos\kappa_{2}$
$\displaystyle-\left(1+\nu^{-1}-w_{1}\nu^{-1}x_{1}\cos\kappa_{1}\right)$
$\displaystyle\times x_{1}\cot\kappa_{2}\sin\kappa_{1}\cos\Delta\zeta$
$\displaystyle+\left(1+\nu-w_{2}\nu x_{2}\cos\kappa_{2}\right)$
$\displaystyle\times x_{2}\cot\kappa_{1}\sin\kappa_{2}\cos\Delta\zeta$
$\displaystyle-
x_{1}x_{2}\left(\frac{\sin\kappa_{2}}{\sin\kappa_{1}}-\frac{\sin\kappa_{1}}{\sin\kappa_{2}}\right)\cos\Delta\zeta\leavevmode\nobreak\
,$

where we have introduced the notations

$R=\frac{3\eta\pi}{\mathfrak{T}\leavevmode\nobreak\
\mathfrak{\bar{l}}_{r}^{2}}\leavevmode\nobreak\ ,\leavevmode\nobreak\
x_{i}=\frac{\chi_{i}}{\mathfrak{\bar{l}}_{r}}\leavevmode\nobreak\ .$ (47)

In Appendix B, we prove that the evolutions remain regular across manifest
coordinate singularities.

In the next two sections, we will discuss applications for this closed system.

##  V Spin flip-flops when one spin dominates over the other

In this section, we consider the case of one spin dominating over the other,
thus $\chi_{2}\ll\chi_{1}$. Under these conditions we will recover
configurations with large spin flip-flop already discussed in the literature
and an additional flip-flop induced by a particular value of the quadrupole
parameter, relevant for neutron star binaries.

Since

$\mathcal{O}\left(\overline{\frac{d\kappa_{1}}{d\mathfrak{t}}}\right)/\mathcal{O}\left(\overline{\frac{d\kappa_{2}}{d\mathfrak{t}}}\right)\approx\frac{\nu\chi_{2}}{\chi_{1}}\ll
1\leavevmode\nobreak\ ,$ (48)

in this case, $\kappa_{1}$ behaves as a quasiconstant. The system governing
$\Delta\zeta$ and $\kappa_{2}$ simplifies to

$\overline{\frac{d\Delta\zeta}{d\mathfrak{t}}}=A_{S_{1}}+B_{S_{1}}\cot\kappa_{2}\cos\Delta\zeta\leavevmode\nobreak\
,$ (49)
$\overline{\frac{d\kappa_{2}}{d\mathfrak{t}}}=B_{S_{1}}\sin\Delta\zeta\leavevmode\nobreak\
,$ (50)

with coefficients

$\displaystyle\frac{A_{S_{1}}}{R}$ $\displaystyle=$
$\displaystyle\nu-\frac{1}{\nu}+\left(1+\frac{2}{\nu}-w_{1}\right.$ (51)
$\displaystyle\left.-\frac{w_{1}x_{1}}{\nu}\cos\kappa_{1}\right)x_{1}\cos\kappa_{1}\leavevmode\nobreak\
,$

and

$\frac{B_{S_{1}}}{R}=-\left(1+\frac{1}{\nu}-\frac{w_{1}x_{1}}{\nu}\cos\kappa_{1}\right)x_{1}\sin\kappa_{1}\leavevmode\nobreak\
.$ (52)

For the values $\kappa_{1}=\left\\{0,\pi\right\\}$, the angle $\kappa_{2}$
becomes a constant, and $\Delta\zeta=A_{S_{1}}\mathfrak{t}+$constant.

For generic $\kappa_{1}$, the system has fix points given by
$\overline{d\Delta\zeta}/d\mathfrak{t}=0$ and
$\overline{d\kappa_{2}}/d\mathfrak{t}=0$, resulting in either

$\Delta\zeta=0\leavevmode\nobreak\ ,\leavevmode\nobreak\
\kappa_{2}=\arctan\left(-\frac{B_{S_{1}}}{A_{S_{1}}}\right)\leavevmode\nobreak\
$ (53)

or

$\Delta\zeta=\pi\leavevmode\nobreak\ ,\leavevmode\nobreak\
\kappa_{2}=\arctan\left(\frac{B_{S_{1}}}{A_{S_{1}}}\right)\leavevmode\nobreak\
.$ (54)

In the rest of the cases [other than $\kappa_{1}=\left\\{0,\pi\right\\}$ or
the fix points (53) or (54)], we derive the following second-order
differential equation from Eqs. (49)–(50),

$\overline{\frac{d^{2}}{d\mathfrak{t}^{2}}}\left(\sin\kappa_{2}\sin\Delta\zeta\right)+\Omega_{S_{1}}^{2}\sin\kappa_{2}\sin\Delta\zeta=0\leavevmode\nobreak\
,$ (55)

where

$\Omega_{S_{1}}=\sqrt{A_{S_{1}}^{2}+B_{S_{1}}^{2}}\leavevmode\nobreak\ .$ (56)

Equation (55) represents a harmonic oscillator, and its solution reads

$\sin\kappa_{2}\sin\Delta\zeta=K_{1}\cos\left(\Omega_{S_{1}}\mathfrak{\bar{t}}+D_{S_{1}}\right)\leavevmode\nobreak\
,$ (57)

with integration constants $\left|K_{1}\right|\leq 1$ and $D_{S_{1}}$. Then,
from (50), we obtain

$\cos\kappa_{2}=-\frac{K_{1}B_{S_{1}}}{\Omega_{S_{1}}}\sin\left(\Omega_{S_{1}}\mathfrak{\bar{t}}+D_{S_{1}}\right)+K_{2}\leavevmode\nobreak\
,$ (58)

where $\left|K_{2}\right|\leq 1$ is an additional integration constant. Since
the system (49)–(50) of two first-order differential equations admits only two
integration constants, $K_{1}$, $K_{2}$, and $D_{S_{1}}$ are not independent.
Indeed, by taking the derivative of Eq. (57) and using Eqs. (49)–(50) and
(58), we find that

$A_{S_{1}}\sin\kappa_{2}\cos\Delta\zeta=\frac{A_{S_{1}}^{2}}{B_{S_{1}}}\cos\kappa_{2}-\frac{\Omega_{S_{1}}^{2}}{B_{S_{1}}}K_{2}\leavevmode\nobreak\
.$ (59)

Hence $A_{S_{1}}=0$ also implies $K_{2}=0$. We introduce a new constant $C$ as

$K_{2}=A_{S_{1}}C\leavevmode\nobreak\ $ (60)

with arbitrary value for $A_{S_{1}}=0$ and from (59) and (57)–(58)

$K_{1}^{2}=1-\Omega_{S_{1}}^{2}C^{2}$ (61)

otherwise. Hence the solutions (57) and (58) can be rewritten as

$\cos\kappa_{2}=A_{S_{1}}C+\frac{\epsilon
B_{S_{1}}}{\Omega_{S_{1}}}\sqrt{1-\Omega_{S_{1}}^{2}C^{2}}\sin\left(\Omega_{S_{1}}\mathfrak{\bar{t}}+D_{S_{1}}\right)\leavevmode\nobreak\
,$ (62)
$\sin\kappa_{2}\sin\Delta\zeta=-\epsilon\sqrt{1-\Omega_{S_{1}}^{2}C^{2}}\cos\left(\Omega_{S_{1}}\mathfrak{\bar{t}}+D_{S_{1}}\right)\leavevmode\nobreak\
,$ (63)

where $\epsilon=\pm 1$ and $-\epsilon$ gives the sign of $K_{1}$.333Note that
the reverse case $\chi_{1}\ll\chi_{2}$ can be obtained by the following
notational changes: $\chi_{1}\rightarrow\chi_{2}$, $\nu\rightarrow\nu^{-1}$,
$w_{1}\rightarrow w_{2}$, $\kappa_{1}\rightarrow\kappa_{2}$, and
$\Delta\zeta\rightarrow-\Delta\zeta$.

From Eq. (62), the maximal variation of $\cos\kappa_{2}$ during the evolution
is

$\displaystyle\left|\Delta\cos\kappa_{2}\right|$ $\displaystyle=$
$\displaystyle\left|\frac{2B_{S_{1}}}{\Omega_{S_{1}}}\sqrt{1-\Omega_{S_{1}}^{2}C^{2}}\right|$
(64) $\displaystyle=$
$\displaystyle\left|2\sqrt{\frac{B_{S_{1}}^{2}}{A_{S_{1}}^{2}+B_{S_{1}}^{2}}-B_{S_{1}}^{2}C^{2}}\right|\leavevmode\nobreak\
.$

This is the largest when

$A_{S_{1}}^{2}<<B_{S_{1}}^{2}\leavevmode\nobreak\ .$ (65)

As $x_{1}=\chi_{1}/\mathfrak{\bar{l}}_{r}$ and $1/\mathfrak{\bar{l}}_{r}$
represents 1/2PN order (Ref. chameleon ), a necessary condition for Eq. (65)
is $\nu\approx 1$, since in this case the leading-order term of
$A_{S_{1}}^{2}$ is negligible. Further the condition (65), implying a large
change in the polar angle of the smaller spin cf. Eq. (64), holds in two
cases,

$i)\leavevmode\nobreak\
\cos\kappa_{1}=\mathcal{O}\left(\frac{1}{\mathfrak{\bar{l}}_{r}}\right)\leavevmode\nobreak\
,$ (66)

and

$ii)\leavevmode\nobreak\
1+\frac{2}{\nu}-w_{1}-\frac{w_{1}x_{1}}{\nu}\cos\kappa_{1}=\mathcal{O}\left(\frac{1}{\mathfrak{\bar{l}}_{r}}\right)\leavevmode\nobreak\
.$ (67)

The condition (66) means that the larger spin vector is almost perpendicular
to $\mathbf{L}_{N}$ (with a deviation of 1/4PN order allowed). This
configuration, which results in large flip-flops of the smaller spin was
analyzed first in Ref. flipflop  for quasicircular orbits. By contrast, as
$\kappa_{1}$ is contained in a
$\mathcal{O}\left(1/\mathfrak{\bar{l}}_{r}\right)$ term of Eq. (67), the
second condition could hold for a wide range of angles $\kappa_{1}$ (i.e.
almost independently of the direction of the dominant spin), provided the
binary component with dominant spin has a quadrupole moment $w_{1}\approx 3$.
This situation can be relevant for neutron star - neutron star binary systems
where one of the binary components is spinning much faster than the other.
Both flip-flop situations arising under i) and ii) are represented on Fig. 5
in a combined fashion as the red diamond-shaped regions. Case i) occurs along
the horizontal axis, while case ii) along the vertical one.

Figure 5: The color map represents the flip-flop angle $\Delta\kappa_{2}$ as a
function of the initial polar angle $\kappa_{1}\left(0\right)$ and the
quadrupole parameter $w_{1}$ of the dominant spin ($\chi_{1}\gg\chi_{2}$).
From the top left to lower right panel, $\Delta\zeta\left(0\right)$ is ${0}$,
${\pi/4}$, ${\pi/2}$, and ${3\pi/4}$, respectively. The additional parameters
of the binary are $m_{1}=m_{2}=1.4M_{\odot}$, $\bar{e}_{r}=0.1$,
$\bar{\varepsilon}=0.0001$, $\chi_{1}=0.95$, $\chi_{2}=0.05$, and
$\kappa_{2}\left(0\right)=\pi/10$. In the small white rectangle regions,
$\Delta\kappa_{2}$ does not reach its maximum within the conservative
timescale. The flip-flop is large in the red diamond-shaped regions.

##  VI Binaries with a black hole and a gravastar, another black hole or boson
star companion

###  VI.1 Black hole - boson star binaries

In this subsection we discuss analytically the case of black hole - boson star
binaries with $w_{1}=1$ and $w_{2}\approx 100$, also equal masses and at large
separations, such that $x_{1}\approx x_{2}\ll 1$ (as they are of the order of
1/2PN order; also we assume $w_{2}x_{2}\ll 1$). The closed system (44)-(46)
can be approximated as

$\frac{1}{R}\overline{\frac{d\kappa_{1}}{d\mathfrak{t}}}=2x_{2}\sin\kappa_{2}\sin\Delta\zeta\leavevmode\nobreak\
,$ (68)
$\frac{1}{R}\overline{\frac{d\kappa_{2}}{d\mathfrak{t}}}=-2x_{1}\sin\kappa_{1}\sin\Delta\zeta\leavevmode\nobreak\
,$ (69)
$\displaystyle\frac{1}{R}\overline{\frac{d\Delta\zeta}{d\mathfrak{t}}}$
$\displaystyle=$ $\displaystyle 2x_{1}\cos\kappa_{1}+w_{2}x_{2}\cos\kappa_{2}$
(70) $\displaystyle-2x_{1}\cot\kappa_{2}\sin\kappa_{1}\cos\Delta\zeta$
$\displaystyle+2x_{2}\cot\kappa_{1}\sin\kappa_{2}\cos\Delta\zeta\leavevmode\nobreak\
.$

The first two equations give

$\overline{\frac{d\cos\kappa_{1}}{d\cos\kappa_{2}}}=-\frac{x_{2}}{x_{1}}\leavevmode\nobreak\
,$ (71)

or

$\cos\kappa_{2}=\frac{x_{1}}{x_{2}}\left(A-\cos\kappa_{1}\right)\leavevmode\nobreak\
,$ (72)

with $A$ a suitable constant to render $\cos\kappa_{2}$ in the allowed range.
If $\cos\kappa_{2}\neq 0$, the first term of Eq. (70) can also be dropped.

In the particular case $A=0\leavevmode\nobreak\ $and $x_{2}=x_{1}$ Eq. (72)
yields $\kappa_{2}=\pi-\kappa_{1}$.

This is a highly symmetrical configuration, with the evolutions of the two
spin polar angles occurring symmetrically to the orbital plane. In this setup,
Eqs. (68)–(70) reduce to

$\frac{1}{R}\overline{\frac{d\kappa_{1}}{d\mathfrak{t}}}=2x_{1}\sin\kappa_{1}\sin\Delta\zeta\leavevmode\nobreak\
,$ (73)
$\frac{1}{R}\overline{\frac{d\Delta\zeta}{d\mathfrak{t}}}=-w_{2}x_{1}\cos\kappa_{1}\leavevmode\nobreak\
.$ (74)

These suitably combined integrate into

$\cos\Delta\zeta=\ln\left|C\sin^{\frac{w_{2}}{2}}\kappa_{1}\right|\leavevmode\nobreak\
,$ (75)

with $C$ an integration constant. Inserting this into Eq. (73), we obtain an
ordinary differential equation

$\displaystyle\frac{1}{R}\overline{\frac{d\kappa_{1}}{d\mathfrak{t}}}$
$\displaystyle=$ $\displaystyle\pm 2x_{1}\sin\kappa_{1}$ (76)
$\displaystyle\times\left[1-\left(\ln\left|C\sin^{\frac{w_{2}}{2}}\kappa_{1}\right|\right)^{2}\right]^{1/2}\leavevmode\nobreak\
,$

with formal solution

$\displaystyle 2x_{1}R\left(\mathfrak{t-t}_{0}\right)$ (77) $\displaystyle=$
$\displaystyle\pm\int\frac{d\kappa_{1}}{\left[1-\left(\ln\left|C\sin^{\frac{w_{2}}{2}}\kappa_{1}\right|\right)^{2}\right]^{1/2}\sin\kappa_{1}}\leavevmode\nobreak\
,$

where $\mathfrak{t}_{0}$ is another constant. Hence, the time evolution of the
spin polar angles is given by Eq. (77), while the evolution of the difference
of their azimuthal angles is given by Eq. (75).

###  VI.2 Comparing flip-flops in black hole - black hole, black hole -
gravastar and black hole - boson star binaries of equal mass

We wish to study here the effect of the quadrupole parameter of the companion
compact object to a black hole in the spin flip-flop. The masses of all
compact objects were chosen equal; hence, $\nu=1$. We monitored the evolution
at PN parameter $\bar{\varepsilon}=0.0001$ and eccentricity $\bar{e}_{r}=0.1$,
leading to $\mathfrak{\bar{l}}_{r}=99.5$. We chose the spin magnitudes
$\chi_{1}=\chi_{2}=0.95$, generating the parameters $x\equiv
x_{1}=x_{2}=0.01$. For initial values of the spin angles we picked
$\kappa_{1}\left(0\right)=\kappa_{2}\left(0\right)=\pi/2$ (hence, the spins
lying in the orbital plane), separated by an azimuthal angle difference of
$\Delta\zeta\left(0\right)=3\pi/4$. With these values we monitored the spin
angle evolutions for the three distinct binary systems. The results are
represented in Fig. 6.

Figure 6: The evolution of the spin polar angles (left column) and of the
difference of their azimuthal angle (right column) as function of the
dimensionless time $\mathfrak{t}$ in units of $xR$, for binaries consisting of
a black hole and either a gravastar with $w=-0.8$ (first row), another black
hole (second row) or a boson star with $w=100$ (third row). The mass ratio is
$\nu=1$ in all cases. The plots are for the PN parameter
$\bar{\varepsilon}=0.0001$, eccentricity $\bar{e}_{r}=0.1$, the parameters
$x\equiv x_{1}=x_{2}=0.01$, and initial spin orientations
$\kappa_{1}\left(0\right)=\kappa_{2}\left(0\right)=\pi/2$,
$\Delta\zeta\left(0\right)=3\pi/4$. The polar angles exhibit a flip-flopping
behavior, with both the frequency and amplitude depending on the quadrupole
parameter of the compact companion. The difference of the azimuthal angles is
also sensitive to $w$.

In all cases, the spin polar angles exhibit the flip-flopping behavior, its
amplitude and frequency being affected by the nature of the companion (through
the value of the quandrupolar parameter $w$). We found that the largest
amplitude and period occur when the companion is an identical black hole. Both
for gravastar and boson star companions, the amplitude and period of the flip-
flop decrease. The flip-flop frequency is largely increased for boson star
companions. As concerning the difference of the azimuthal angles
$\Delta\zeta\,$, the evolutions are similar for low values of $w$, thus for
gravastar and black hole companions. The spins of the binary components
exhibit precessions during which they overpass each other periodically
($\Delta\zeta$ evolves through a sequence of positive and negative values).
Nevertheless when the companion is a boson star, this swinging behavior
disappears ($\Delta\zeta$ remains positive). The amplitude of the periodical
evolution in $\Delta\zeta$ is also decreased.

##  VII Conclusions

In this paper we derived the conservative secular evolution of precessing
compact binaries on eccentric orbits, to second post-Newtonian–order accuracy,
with leading-order spin-orbit, spin-spin and mass quadrupole-monopole
contributions included. Our approach relies on employing a properly chosen set
of dimensionless variables, advanced in Ref. chameleon  and a method of
averaging over the radial period. The secular dynamics emerged by applying
this to the instantaneous dynamics discussed in Ref. chameleon . The inclusion
of the mass quadrupole parameter allows to apply the formalism for binaries
with arbitrary compact components, like black holes, neutron stars, boson
stars or gravastars.

The derived secular dynamics generalizes previous results from the literature.
In Ref. GPV3  the dynamics was only expressed to 1.5 PN order by employing
different shape variables. In Ref. ACST , the precession was examined with
leading-order SO and SS effects for orbits with negligible eccentricity.
Reference damourderuelle  investigated the PN dynamics with eccentricity, but
without spins. Reference WOEdyn  discussed the dynamics of small mass ratio
binaries with only the smaller body having spin. The seminal review on
gravitational radiation from compact binary sources by Blanchet blanchet
discusses the SO effect in its last section, but omits the SS and QM
contributions to the dynamics. The secular dynamics where the leading-order
SO, SS, and QM coupling are included is investigated analitically in Ref.
Racine  only for black holes.

The secular evolution equations emerged as a closed system of first-order
differential equations, which in contrast with the instantaneous evolutions
presented in Ref. chameleon , is autonomous. The dependent variables are the
polar ($\kappa_{1}$ and $\kappa_{2}$) and the azimuthal angles ($\zeta_{1}$
and $\zeta_{2}$) of the spin vectors, the angles $\alpha$ and $\phi_{n}$
giving the orientation of the orbital angular momentum vector, together with
the periastron angle $\psi_{p}$, the dimensionless magnitude of the orbital
angular momentum $\mathfrak{l}_{r}$ and the eccentricity $e_{r}$. Over the
conservative timescale the secular dynamics can be regarded as some sort of
smoothed-out intantaneous evolution, as illustrated on Figs. 1–4.

Moreover we have shown that the spin polar angles and the difference of their
azimuthal angles in the system defined by the orbital plane evolve according
to a closed subsystem of the secular dynamics. In spite of the apparent
singularity of spherical polar coordinates, the evolutions remains well
defined through aligned configurations. We studied in detail this closed
subsystem in two significant cases. First we assumed that the masses are
comparable, but one of the spins dominates over the other. In this case we i)
derived analytically that large flip-flops of the smaller spin emerge when the
larger spin is almost coplanar with the orbit (a known result) and ii) found
new flip-flop configurations arising for the quadrupole parameter
$w_{1}\approx 3$ of the neutron star with dominant spin.

We also studied black hole - boson star binaries. In this case, the huge
quadrupolar parameter of the boson star allows for significant simplification
of the closed subsystem, allowing us to derive a formal analytical solution.

Finally, we analyzed the evolutions of the spin angles numerically by
comparing the cases when the black hole companion is either a gravastar,
another black hole or a boson star with identical mass. We found that the
amplitude and period of the flip-flop is maximal, when the companion is a
black hole. In the case of a boson star companion the frequency of the flip-
flop increased significantly. The precession of the spins is also sensitive to
the quadrupolar parameter of the companion. While in the case of gravastars
and black holes a swinging-type evolution occurs, when the spins of the
components regularly overpass each other, their sequence is conserved when the
companion is a boson star.

In a related paper KGstab  we will discuss the equilibrium of the spin
configurations and their linear stability in precessing compact binaries with
black hole, neutron star, gravastar, or boson star components.

## Acknowledgements

This work was supported by the Hungarian National Research Development and
Innovation Office (NKFIH) in the form of the Grant No. 123996 and has been
carried out in the framework of COST actions CA16104 (GWverse) and CA18108
(QG-MM) supported by COST (European Cooperation in Science and Technology). Z.
K. was further supported by the János Bolyai Research Scholarship of the
Hungarian Academy of Sciences and by the ÚNKP-20-5–New National Excellence
Program of the Ministry for Innovation and Technology through its National
Research, Development and Innovation Fund. In the early stages of this work L.
Á. G. was supported by the European Union and the State of Hungary, cofinanced
by the European Social Fund in the framework of TÁMOP 4.2.4.
A/2-11-/1-2012-0001 “National Excellence Program.”

##  Appendix A The PN expansion of the radial period

When averaging instantaneous variables over one radial orbit in order to
obtain their secular counterpart, the expression of $\dot{\chi}_{p}$ given by
Eq. (43) of Ref. chameleon  is needed, as explained in Sec. II.1, with both
its Newtonian and PN contributions expressible by the true anomaly $\chi_{p}$
and the shape variables $\mathfrak{l}_{r}\left(\chi_{p}\right)$ and
$e_{r}\left(\chi_{p}\right)$ alone. Hence
$\mathfrak{l}_{r}\left(\chi_{p}\right)$ and $e_{r}\left(\chi_{p}\right)$ will
be required to 2PN-order accuracy, while the rest of the orbital elements
$\psi_{p}$, $\alpha$, and $\phi_{n}$, and spin angles $\kappa_{i}$ and
$\zeta_{i}$ ($i=1,2$) only to leading order, where they are constants.

In this Appendix first we derive the $\chi_{p}$ dependence of the
dimensionless orbital angular momentum $\mathfrak{l}_{r}$ and the
dimensionless orbital eccentricity $e_{r}$ to 2PN order, in terms of their
values at the periastron (characterized by the true anomaly $\chi_{p}=0$).
Next, employ these expressions to compute the radial period to 2PN accuracy.
The derivation of the time-averaged values $\bar{\mathfrak{l}}_{r}$ and
$\bar{e}_{r}$ over the radial period follows, again to 2PN accuracy, with the
inclusion of all spin and mass quadrupole effects to this order. This enables
us to express the shape variables evaluated at the periastron in terms of the
corresponding averaged quantities. This is employed for rewriting the radial
period as a PN expansion in terms of averaged quantities. At the end of the
Appendix we also give a similar expansion of the averaged PN parameter.

###  A.1 $\chi_{p}$ dependence of $\mathfrak{l}_{r}$

A lower index $0$ indicates values taken at $\chi_{p}=0$:

$\displaystyle\mathfrak{l}_{r}\left(\chi_{p}=0\right)$ $\displaystyle=$
$\displaystyle\mathfrak{l}_{r0}\leavevmode\nobreak\ ,$ (78) $\displaystyle
e_{r}\left(\chi_{p}=0\right)$ $\displaystyle=$ $\displaystyle
e_{r0}\leavevmode\nobreak\ .$ (79)

The expressions $\mathfrak{\dot{l}}_{r}$ and $\dot{\chi}_{p}$ given by Eqs.
(36) and (43) of Ref. chameleon  allow for deriving

$\displaystyle\mathfrak{l}_{r}\left(\chi_{p}\right)$ $\displaystyle=$
$\displaystyle\mathfrak{l}_{r0}+\int_{0}^{\chi_{p}}\frac{\mathfrak{\dot{l}}_{r}}{\dot{\chi}_{p}}d\chi_{p}$
(80) $\displaystyle=$
$\displaystyle\mathfrak{l}_{r0}+\frac{\mathfrak{l}_{rPN}\left(\chi_{p}\right)}{\mathfrak{l}_{r0}}+\frac{\mathfrak{l}_{rSO}\left(\chi_{p}\right)}{\mathfrak{l}_{r0}^{2}}$
$\displaystyle+\frac{\mathfrak{l}_{rQM}\left(\chi_{p}\right)}{\mathfrak{l}_{r0}^{3}}+\frac{\mathfrak{l}_{rSS}\left(\chi_{p}\right)}{\mathfrak{l}_{r0}^{3}}$
$\displaystyle+\frac{\mathfrak{l}_{r2PN}\left(\chi_{p}\right)}{\mathfrak{l}_{r0}^{3}}\leavevmode\nobreak\
,$

with

$\mathfrak{l}_{rPN}\left(\chi_{p}\right)=2\left(2-\eta\right)e_{r0}\left(1-\cos\chi_{p}\right)\leavevmode\nobreak\
,$ (81)
$\mathfrak{l}_{r2PN}\left(\chi_{p}\right)=\sum_{k=0}^{4}\sum_{l=0}^{2}L_{kl}^{2PN}\sin^{2l}\chi_{p}\cos^{k}\chi_{p}\leavevmode\nobreak\
,$ (82) $\displaystyle\mathfrak{l}_{rSO}\left(\chi_{p}\right)$
$\displaystyle=$ $\displaystyle-\frac{\eta
e_{r0}}{2}\left(\cos\chi_{p}-1\right)$ (83)
$\displaystyle\times\sum_{k=1}^{2}\left(4^{2k-3}+3\right)\chi_{k}\cos\kappa_{k}\leavevmode\nobreak\
,$ $\displaystyle\mathfrak{l}_{rSS}\left(\chi_{p}\right)$ $\displaystyle=$
$\displaystyle\eta\chi_{1}\chi_{2}\sin\kappa_{1}\sin\kappa_{2}$ (84)
$\displaystyle\times\left[\cos\zeta_{+}\sum_{k=0}^{3}L_{k}^{SS}\cos^{k}\chi_{p}\right.$
$\displaystyle\left.+\sin\zeta_{+}\sin\chi_{p}\sum_{k=0}^{2}K_{k}^{SS}\cos^{k}\chi_{p}\right]\leavevmode\nobreak\
,$ $\displaystyle\mathfrak{l}_{rQM}\left(\chi_{p}\right)$ $\displaystyle=$
$\displaystyle\frac{\eta}{2}\sum_{i=1}^{2}\chi_{i}^{2}w_{i}\nu^{2i-3}\sin^{2}\kappa_{i}$
(85) $\displaystyle\times\left[\cos
2\zeta_{i}\sum_{k=0}^{3}L_{k}^{QM}\cos^{k}\chi_{p}\right.$
$\displaystyle\left.+\sin
2\zeta_{i}\sin\chi_{p}\sum_{k=0}^{2}K_{k}^{QM}\cos^{k}\chi_{p}\right]\leavevmode\nobreak\
,$

where the coefficients $L_{kl}^{2PN}$, $L_{k}^{SS}$, $K_{k}^{SS}$,
$L_{k}^{QM}$, and $K_{k}^{QM}$ are enlisted in Table 1.

Table 1: The coefficients of $\mathfrak{l}_{r}\left(\chi_{p}\right)$ in Eqs. (82), (84), and (85). $Coefficient$ | $Expression$
---|---
$L_{00}^{2PN}$ | $\frac{1}{96}\left[432e_{r0}^{3}\eta+e_{r0}^{2}\left(-117\eta^{2}+54\eta+48\right)\right.$
| $+32e_{r0}\left(2\eta^{2}-83\eta+50\right)$
| $\left.-48\left(\eta^{2}-5\eta+6\right)\right]$
$L_{10}^{2PN}$ | $\frac{e_{r0}}{8}\left[\left(2\eta-33\right)\eta e_{r0}^{2}\right.$
| $\left.-116\eta^{2}+256\eta-160\right]$
$L_{20}^{2PN}$ | $\frac{1}{8}\left[e_{r0}^{2}\left(9\eta^{2}-3\eta-4\right)\right.$
| $\left.+4\left(\eta^{2}-5\eta+6\right)\right]$
$L_{30}^{2PN}$ | $\frac{e_{r0}}{24}\left[-3\left(2\eta+3\right)\eta e_{r0}^{2}\right.$
| $\left.+32\eta^{2}-104\eta+80\right]$
$L_{40}^{2PN}$ | $\frac{3e_{r0}^{2}(\eta-2)\eta}{32}$
$L_{01}^{2PN}$ | $\frac{1}{8}\left[e_{r0}^{2}\left(-9\eta^{2}+3\eta+4\right)\right.$
| $\left.-4\left(\eta^{2}-5\eta+6\right)\right]$
$L_{11}^{2PN}$ | $\frac{e_{r0}}{8}\left[3\left(2\eta+3\right)\eta e_{r0}^{2}\right.$
| $\left.-32\eta^{2}+104\eta-80\right]$
$L_{21}^{2PN}$ | $-\frac{9e_{r0}^{2}(\eta-2)\eta}{16}$
$L_{31}^{2PN}$ | $0$
$L_{41}^{2PN}$ | $0$
$L_{02}^{2PN}$ | $\frac{3e_{r0}^{2}(\eta-2)\eta}{32}$
$L_{12}^{2PN}$ | $0$
$L_{22}^{2PN}$ | $0$
$L_{32}^{2PN}$ | $0$
$L_{42}^{2PN}$ | $0$
$L_{0}^{SS}$ | $2e_{r0}+3$
$L_{1}^{SS}$ | $0$
$L_{2}^{SS}$ | $-3$
$L_{3}^{SS}$ | $-2e_{r0}$
$K_{0}^{SS}$ | $-e_{r0}$
$K_{1}^{SS}$ | $-3$
$K_{2}^{SS}$ | $-2e_{r0}$
$L_{0}^{QM}$ | $-\left(2e_{r0}+3\right)$
$L_{1}^{QM}$ | $0$
$L_{2}^{QM}$ | $-3$
$L_{3}^{QM}$ | $-2e_{r0}$
$K_{0}^{QM}$ | $-e_{r0}$
$K_{1}^{QM}$ | $-3$
$K_{2}^{QM}$ | $-2e_{r0}$

###  A.2 $\chi_{p}$ dependence of $e_{r}$

From $\dot{e}_{r}$ and $\dot{\chi}_{p}$ given by Eqs. (37) and (43) of Ref.
chameleon , respectively, we find

$\displaystyle e_{r}\left(\chi_{p}\right)$ $\displaystyle=$ $\displaystyle
e_{r0}+\int_{0}^{\chi_{p}}\frac{\dot{e}_{r}}{\dot{\chi}_{p}}d\chi_{p}\leavevmode\nobreak\
=$ (86) $\displaystyle=$ $\displaystyle
e_{r0}+\frac{1}{\mathfrak{l}_{r0}^{2}}e_{rPN}\left(\chi_{p}\right)+\frac{1}{\mathfrak{l}_{r0}^{3}}e_{rSO}\left(\chi_{p}\right)$
$\displaystyle+\frac{1}{\mathfrak{l}_{r0}^{4}}e_{rQM}\left(\chi_{p}\right)+\frac{1}{\mathfrak{l}_{r0}^{4}}e_{rSS}\left(\chi_{p}\right)$
$\displaystyle+\frac{1}{\mathfrak{l}_{r0}^{4}}e_{r2PN}\left(\chi_{p}\right)\leavevmode\nobreak\
,$

with

$e_{rPN}\left(\chi_{p}\right)=\sum_{k=0}^{3}E_{k}^{PN}\cos^{k}\chi_{p}\leavevmode\nobreak\
,$ (87) $\displaystyle e_{rSO}\left(\chi_{p}\right)$ $\displaystyle=$
$\displaystyle\frac{\eta}{2}\left(1-e_{r0}^{2}\right)\left(1-\cos\chi_{p}\right)\times$
(88)
$\displaystyle\times\sum_{i=1}^{2}\left(4^{2k-3}+3\right)\chi_{i}\cos\kappa_{i}\leavevmode\nobreak\
,$
$e_{r2PN}\left(\chi_{p}\right)=\sum_{l=0}^{3}\sum_{k=0}^{6}E_{kl}^{2PN}\cos^{k}\chi_{p}\sin^{2l}\chi_{p}\leavevmode\nobreak\
,$ (89) $\displaystyle
e_{rSS}\left(\chi_{p}\right)=\eta\chi_{1}\chi_{2}\left[\sum_{k=0}^{5}E_{k}^{SS}\cos^{k}\chi_{p}\right.$
$\displaystyle\left.+\sin\kappa_{1}\sin\kappa_{2}\sin\zeta_{+}\sin\chi_{p}\sum_{k=0}^{4}F_{k}^{SS}\cos^{k}\chi_{p}\right]\leavevmode\nobreak\
,$ (90) $\displaystyle e_{rQM}\left(\chi_{p}\right)$ $\displaystyle=$
$\displaystyle\frac{\eta}{2}\sum_{i=1}^{2}\chi_{i}^{2}\nu^{2i-3}w_{i}$ (91)
$\displaystyle\times\left[\sum_{k=0}^{6}E_{k}^{QM}\cos^{k}\chi_{p}\right.$
$\displaystyle+\sin 2\zeta_{i}\sin^{2}\kappa_{i}\sin\chi_{p}$
$\displaystyle\left.\times\sum_{k=0}^{4}F_{k}^{QM}\cos^{k}\chi_{p}\right]\leavevmode\nobreak\
.$

The coefficients $E_{kl}^{2PN}$, $E_{k}^{SS}$, $F_{k}^{SS}$, $E_{k}^{QM}$, and
$F_{k}^{QM}\leavevmode\nobreak\ $of $e_{r}\left(\chi_{p}\right)$ for the PN,
SS, and QM contributions are collected in Table 2. The coefficients
$E_{k}^{2PN}$ and $F_{k}^{2PN}$ are enlisted in Table 3.

Table 2: The coefficients of $e_{r}\left(\chi_{p}\right)$ in Eqs. (87), (90), and (91).  $Coefficient$ | $Expression$
---|---
$E_{0}^{PN}$ | $3-\eta+\left(5-4\eta\right)e_{r0}+e_{r0}^{2}\left(7-6\eta\right)$
$E_{1}^{PN}$ | $-\left[3-\eta+e_{r0}^{2}\left(7-\frac{11}{2}\eta\right)\right]$
$E_{2}^{PN}$ | $-\left(5-4\eta\right)e_{r0}$
$E_{3}^{PN}$ | $\frac{\eta}{2}e_{r0}^{2}$
$E_{0}^{SS}$ | $-\cos\kappa_{1}\cos\kappa_{2}\left(e_{r0}^{2}+3e_{r0}+3\right)$
|
$+\frac{1}{2}\sin\kappa_{1}\sin\kappa_{2}\left[\left(e_{r0}^{2}+3e_{r0}+3\right)\cos\zeta_{-}\right.$
| $+\left.\left(7e_{r0}^{2}+15e_{r0}+5\right)\cos\zeta_{+}\right]$
$E_{1}^{SS}$ | $\frac{3}{2}\sin\kappa_{1}\sin\kappa_{2}\left(3\cos\zeta_{+}-\cos\zeta_{-}\right)$
| $+3\cos\kappa_{1}\cos\kappa_{2}$
$E_{2}^{SS}$ | $\frac{3}{2}e_{r0}\sin\kappa_{1}\sin\kappa_{2}\left(\cos\zeta_{+}-\cos\zeta_{-}\right)$
| $+3e_{r0}\cos\kappa_{1}\cos\kappa_{2}$
$E_{3}^{SS}$ | $-\frac{1}{2}e_{r0}^{2}\sin\kappa_{1}\sin\kappa_{2}\left(\cos\zeta_{+}+\cos\zeta_{-}\right)$
| $-7\sin\kappa_{1}\sin\kappa_{2}\cos\zeta_{+}$
| $+e_{r0}^{2}\cos\kappa_{1}\cos\kappa_{2}$
$E_{4}^{SS}$ | $-9e_{r0}\sin\kappa_{1}\sin\kappa_{2}\cos\zeta_{+}$
$E_{5}^{SS}$ | $-3e_{r0}^{2}\sin\kappa_{1}\sin\kappa_{2}\cos\zeta_{+}$
$F_{0}^{SS}$ | $1-e_{r0}^{2}$
$F_{1}^{SS}$ | $-3e_{r0}$
$F_{2}^{SS}$ | $-\left(2e_{r0}^{2}+7\right)$
$F_{3}^{SS}$ | $-9e_{r0}$
$F_{4}^{SS}$ | $-3e_{r0}^{2}$
$E_{0}^{QM}$ | $\left(7e_{r0}^{2}+15e_{r0}+5\right)\cos^{2}\zeta_{i}\sin^{2}\kappa_{i}$
| $+\left(2e_{r0}^{2}+3e_{r0}-2\right)\cos^{2}\kappa_{i}$
| $-3e_{r0}^{2}-6e_{r0}-1$
$E_{1}^{QM}$ | $9\cos^{2}\kappa_{i}\sin^{2}\zeta_{i}+3\left(3\cos^{2}\zeta_{i}-2\right)$
$E_{2}^{QM}$ | $3e_{r0}\left(2-\cos^{2}\zeta_{i}\right)\cos^{2}\kappa_{i}$
| $-3e_{r0}\sin^{2}\zeta_{i}$
$E_{3}^{QM}$ | $e_{r0}^{2}\cos^{2}\kappa_{i}$
|
$+\left[-\left(e_{r0}^{2}+14\right)\cos^{2}\zeta_{i}+7\right]\sin^{2}\kappa_{i}$
$E_{4}^{QM}$ | $-9e_{r0}\sin^{2}\kappa_{i}\cos 2\zeta_{i}$
$E_{5}^{QM}$ | $-3e_{r0}^{2}\sin^{2}\kappa_{i}\cos 2\zeta_{i}$
$F_{0}^{QM}$ | $\left(1-e_{r0}^{2}\right)$
$F_{1}^{QM}$ | $-3e_{r0}$
$F_{2}^{QM}$ | $-\left(2e_{r0}^{2}+7\right)$
$F_{3}^{QM}$ | $-9e_{r0}$
$F_{4}^{QM}$ | $-3e_{r0}^{2}$
Table 3: The coefficients of $e_{r}\left(\chi_{p}\right)$ in Eq. (89). $Coefficient$ | $Expression$
---|---
$E_{00}^{2PN}$ | $\frac{1}{1920e_{r0}}\left[1920e_{r0}^{5}\eta(3\eta+8)\right.$
| $+e_{r0}^{4}\left(-1845\eta^{2}+8880\eta+1800\right)$
| $+32e_{r0}^{3}\left(232\eta^{2}-2825\eta+1740\right)$
| $-180e_{r0}^{2}\left(29\eta^{2}+89\eta+60\right)$
| $+160e_{r0}\left(8\eta^{2}-187\eta-60\right)$
| $\left.-480(\eta-3)^{2}\right]$
$E_{10}^{2PN}$ | $-\frac{1}{64}\left[e_{r0}^{4}\eta\left(161\eta+477\right)\right.$
| $+4e_{r0}^{2}\left(136\eta^{2}-849\eta+564\right)$
| $\left.+16\eta\left(8\eta-85\right)\right]$
$E_{20}^{2PN}$ | $\frac{1}{256e_{r0}}\left[e_{r0}^{4}\left(269\eta^{2}-1312\eta-256\right)\right.$
| $+32e_{r0}^{2}\left(5\eta^{2}+109\eta+20\right)$
| $\left.+64\left(\eta-3\right)^{2}\right]$
$E_{30}^{2PN}$ | $\frac{1}{384}\left[-3e_{r0}^{4}\eta\left(53\eta+73\right)\right.$
| $+8e_{r0}^{2}\left(208\eta^{2}-269\eta+300\right)$
| $\left.+128\left(4\eta^{2}-17\eta+15\right)\right]$
$E_{40}^{2PN}$ | $\frac{e_{r0}}{128}\left[e_{r0}^{2}\left(-13\eta^{2}+64\eta+8\right)\right.$
| $\left.+268\eta^{2}-676\eta+400\right]$
$E_{50}^{2PN}$ | $-\frac{3e_{r0}^{2}\eta}{640}\left[5e_{r0}^{2}(3\eta-1)-64\eta+80\right]$
$E_{60}^{2PN}$ | $\frac{3e_{r0}^{3}\eta^{2}}{256}$
$E_{01}^{2PN}$ | $-E_{20}^{2PN}$
$E_{11}^{2PN}$ | $-3E_{30}^{2PN}$
$E_{21}^{2PN}$ | $-\frac{3}{2}E_{40}^{2PN}$
$E_{31}^{2PN}$ | $-\frac{1}{10}E_{50}^{2PN}$
$E_{41}^{2PN}$ | $-9E_{60}^{2PN}$
$E_{51}^{2PN}$ | $0$
$E_{61}^{2PN}$ | $0$
$E_{02}^{2PN}$ | $E_{40}^{2PN}$
$E_{12}^{2PN}$ | $\frac{1}{5}E_{50}^{2PN}$
$E_{22}^{2PN}$ | $9E_{60}^{2PN}$
$E_{32}^{2PN}$ | $0$
$E_{42}^{2PN}$ | $0$
$E_{52}^{2PN}$ | $0$
$E_{62}^{2PN}$ | $0$
$E_{03}^{2PN}$ | $-E_{60}^{2PN}$
$E_{13}^{2PN}$ | $0$
$E_{23}^{2PN}$ | $0$
$E_{33}^{2PN}$ | $0$
$E_{43}^{2PN}$ | $0$
$E_{53}^{2PN}$ | $0$
$E_{63}^{2PN}$ | $0$

###  A.3 Dimensionless 2PN orbital period

We insert the expressions of $\mathfrak{l}_{r}\left(\chi_{p}\right)$ and
$e_{r}\left(\chi_{p}\right)$ into the integral (4) and expand it in Taylor
series to 2PN order.

The integral (4) leads to the PN expansion:

$\displaystyle\mathfrak{T}$ $\displaystyle=$
$\displaystyle\mathfrak{T}_{0}\left(1+\frac{\tau_{0PN}}{\mathfrak{l}_{r}^{2}}+\frac{\tau_{0SO}}{\mathfrak{l}_{r}^{3}}+\frac{\tau_{0SS}}{\mathfrak{l}_{r}^{4}}\right.$
(92)
$\displaystyle\left.+\frac{\tau_{0QM}}{\mathfrak{l}_{r}^{4}}+\frac{\tau_{02PN}}{\mathfrak{l}_{r}^{4}}\right)\leavevmode\nobreak\
,$

where the lower index $0$ stands for $\chi_{p}=0$. As explained in the main
text $\tau_{0PN}/\mathfrak{l}_{r}^{2}$ and $\tau_{02PN}/\mathfrak{l}_{r}^{4}$
give PN- and 2PN-order contributions, respectively.

The terms $\mathfrak{T}_{0}$, $\tau_{0PN}$, $\tau_{0SO}$, $\tau_{0SS}$,
$\tau_{0QM}$ and $\tau_{02PN}$ are found by exploring the expressions
$\mathfrak{l}_{r}\left(\chi_{p}\right)$ and $e_{r}\left(\chi_{p}\right)$
derived above. They read as

$\mathfrak{T}_{0}=\frac{2\pi\mathfrak{l}_{r0}^{3}}{\left(1-e_{r0}^{2}\right)^{3/2}}\leavevmode\nobreak\
,$ (93)
$\tau_{0PN}=-\frac{\left(1-e_{r0}^{2}\right)\left(e_{r0}^{2}(7\eta-6)+2e_{r0}(5\eta-3)+4\eta-18\right)}{2(e_{r0}-1)^{2}}\leavevmode\nobreak\
,$ $\tau_{0SO}=0\leavevmode\nobreak\ ,$ (94) $\displaystyle\tau_{02PN}$
$\displaystyle=$
$\displaystyle\frac{1}{40(1-e_{r0})^{2}e_{r0}^{4}}\left[\sum_{k=0}^{10}U_{k}e_{r0}^{k}\right.$
(95)
$\displaystyle\left.-\frac{\left(1-e_{r0}^{2}\right)^{3/2}}{2(1-e_{r0})}\sum_{k=0}^{7}V_{k}e_{r0}^{k}\right]\leavevmode\nobreak\
,$ $\displaystyle\tau_{0SS}$ $\displaystyle=$
$\displaystyle-\frac{3\chi_{1}\chi_{2}\left(1+e_{r0}\right)^{2}\eta}{\left(1-e_{r0}\right)\mathfrak{l}_{r0}^{4}}\left[\cos\kappa_{1}\cos\kappa_{2}+\sin\kappa_{1}\sin\kappa_{2}\right.$
(96)
$\displaystyle\left.\times\left(\sin\zeta_{1}\sin\zeta_{2}-2\cos\zeta_{1}\cos\zeta_{2}\right)\right]\leavevmode\nobreak\
,$ $\displaystyle\tau_{0QM}$ $\displaystyle=$
$\displaystyle-\frac{3\eta\left(1+e_{r0}\right)^{2}}{2\left(1-e_{r0}\right)\mathfrak{l}_{r0}^{4}}\sum_{i=1}^{2}\chi_{i}^{2}\nu^{2i-3}w_{i}$
(97)
$\displaystyle\times\left(1-3\sin^{2}\kappa_{i}\cos^{2}\zeta_{i}\right)\leavevmode\nobreak\
.$

The coefficients $U_{k}$ and $V_{k}$ of $\tau_{02PN}$ are enlisted in Table 4.

Table 4: The coefficients $U_{k}$ and $V_{k}$ of $\tau_{02PN}$ in Eq. (95). $Coefficient$ | $Expression$
---|---
$U_{10}$ | $-105\left(\eta+1\right)\eta$
$U_{9}$ | $10\left(-559\eta+297\eta^{2}+228\right)$
$U_{8}$ | $5\left(-2674\eta+1289\eta^{2}+1336\right)$
$U_{7}$ | $-4\left(-235\eta+186\eta^{2}-280\right)$
$U_{6}$ | $-482\eta^{2}+25\eta+1120$
$U_{5}$ | $2\left(-17\,245\eta+5496\eta^{2}+12\,200\right)$
$U_{4}$ | $2\left(-2240\eta+79\eta^{2}+3350\right)$
$U_{3}$ | $-4\left(-7075\eta+2137\eta^{2}+5300\right)$
$U_{2}$ | $2\left(-3915\eta+1317\eta^{2}+2050\right)$
$U_{1}$ | $40\left(-254\eta+67\eta^{2}+210\right)$
$U_{0}$ | $-20\left(-238\eta+65\eta^{2}+180\right)$
$V_{7}$ | $5\left(-932\eta+509\eta^{2}+504\right)$
$V_{6}$ | $-5\left(-2700\eta+1399\eta^{2}+1352\right)$
$V_{5}$ | $1427\eta^{2}+920\eta-4240$
$V_{4}$ | $12079\eta^{2}-35880\eta+25680$
$V_{3}$ | $-4\left(-6875\eta+2606\eta^{2}+3650\right)$
$V_{2}$ | $-4\left(-4735\eta+998\eta^{2}+4850\right)$
$V_{1}$ | $40\left(-746\eta+199\eta^{2}+600\right)$
$V_{0}$ | $-40\left(-238\eta+65\eta^{2}+180\right)$

###  A.4 Secular shape variables $\mathfrak{\bar{l}}_{r}$ and $\bar{e}_{r}$

Next we calculate the time averages of the shape variables $\mathfrak{l}_{r}$
and $e_{r}$ during one radial period

$\displaystyle\mathfrak{\bar{l}}_{r}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathfrak{T}}\int_{0}^{2\pi}\frac{\mathfrak{l}_{r}\left(\chi_{p}\right)}{\dot{\chi}_{p}}d\chi_{p}\leavevmode\nobreak\
,$ (98) $\displaystyle\bar{e}_{r}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathfrak{T}}\int_{0}^{2\pi}\frac{e_{r}\left(\chi_{p}\right)}{\dot{\chi}_{p}}d\chi_{p}\leavevmode\nobreak\
,$ (99)

in terms of their initial values at $\chi_{p}=0$, with the decomposition of
$\mathfrak{T}$ given in Eq. (92). Their PN expansion can be formally written
as

$\displaystyle\mathfrak{\bar{l}}_{r}$ $\displaystyle\\!=\\!$
$\displaystyle\frac{\mathfrak{l}_{r0N}\\!+\\!\mathfrak{\bar{l}}_{rPN}\\!+\\!\mathfrak{\bar{l}}_{rSO}\\!+\\!\mathfrak{\bar{l}}_{rSS}\\!+\\!\mathfrak{\bar{l}}_{rQM}\\!+\\!\mathfrak{\bar{l}}_{r2PN}}{\mathfrak{T}}\leavevmode\nobreak\
,$ (100) $\displaystyle\bar{e}_{r}$ $\displaystyle\\!=\\!$
$\displaystyle\frac{e_{r0N}\\!+\\!\bar{e}_{rPN}\\!+\\!\bar{e}_{rSO}\\!+\\!\bar{e}_{rSS}\\!+\\!\bar{e}_{rQM}\\!+\\!\bar{e}_{r2PN}}{\mathfrak{T}}\leavevmode\nobreak\
.$ (101)

The contributions to the integral in Eq. (98) become

$\mathfrak{\bar{l}}_{rN}=\frac{2\pi\mathfrak{l}_{r0}^{4}}{\left(1-e_{r0}^{2}\right)^{3/2}}\leavevmode\nobreak\
,$ (102) $\displaystyle\mathfrak{\bar{l}}_{rPN}$ $\displaystyle=$
$\displaystyle-\frac{\pi\mathfrak{l}_{r0}^{2}}{(1-e_{r0})^{2}\sqrt{1-e_{r0}^{2}}}\left[e_{r0}^{2}(3\eta+2)\right.$
(103)
$\displaystyle\left.+14e_{r0}(\eta-1)+4\eta-18\right]\leavevmode\nobreak\ ,$
$\displaystyle\mathfrak{\bar{l}}_{rSO}$ $\displaystyle=$
$\displaystyle-\frac{\pi
e_{r0}\eta\mathfrak{l}_{r0}}{(1-e_{r0})\sqrt{1-e_{r0}^{2}}}$ (104)
$\displaystyle\times\sum_{k=1}^{2}\left(4^{2k-3}+3\right)\chi_{k}\cos\kappa_{k}\leavevmode\nobreak\
,$ $\displaystyle\mathfrak{\bar{l}}_{r2PN}$ $\displaystyle=$
$\displaystyle\frac{\pi}{120e_{r0}^{4}}\sum_{k=0}^{4}\bar{L}_{k}^{2PN}e_{r0}^{k}$
(105)
$\displaystyle-\frac{\pi\sqrt{1-e_{r0}^{2}}}{60(e_{r0}-1)^{4}e_{r0}^{4}}\sum_{k=0}^{8}\bar{K}_{k}^{2PN}e_{r0}^{k}\leavevmode\nobreak\
,$ $\displaystyle\mathfrak{\bar{l}}_{rSS}$ $\displaystyle=$
$\displaystyle\frac{\pi\chi_{1}\chi_{2}\eta}{16\left(1-e_{r0}\right)e_{r0}^{2}\left(1-e_{r0}^{2}\right)^{3/2}}$
(106)
$\displaystyle\times\left(\cos\kappa_{1}\cos\kappa_{2}\bar{L}^{SS}\right.$
$\displaystyle\left.+\sin\kappa_{1}\sin\kappa_{2}\bar{K}^{SS}\right)\leavevmode\nobreak\
,$ $\displaystyle\mathfrak{\bar{l}}_{rQM}$ $\displaystyle=$
$\displaystyle\frac{\pi\left(e_{r0}^{2}+2\right)\eta}{256e_{r0}\left(1-e_{r0}^{2}\right)^{5/2}}\sum_{k=1}^{2}\chi_{k}^{2}\nu^{2k-3}w_{k}$
(107)
$\displaystyle\times\left[-4\left(47e_{r0}^{3}+1050e_{r0}^{2}+488e_{r0}\right.\right.$
$\displaystyle\left.+480\right)\sin^{2}\kappa_{k}\cos 2\zeta_{k}$
$\displaystyle+16\left(5e_{r0}^{3}+21e_{r0}^{2}+15e_{r0}+6\right)$
$\displaystyle\left.\times\left(3\cos
2\kappa_{k}+1\right)\right]\leavevmode\nobreak\ .$

The coefficients $\bar{L}_{k}^{2PN}$, $\bar{K}_{k}^{2PN}$, $\bar{L}^{SS}$, and
$\bar{K}^{SS}$ are collected in Table 5.

Table 5: The coefficients of $\mathfrak{\bar{l}}_{r}$ in Eqs. (105) and (106). $Coefficient$ | $Expression$
---|---
$\bar{L}_{4}^{2PN}$ | $15\left(467\eta^{2}-580\eta+296\right)$
$\bar{L}_{3}^{2PN}$ | $480\left(4\eta^{2}-3\eta+5\right)$
$\bar{L}_{2}^{2PN}$ | $-4\left(3001\eta^{2}-9445\eta+6610\right)$
$\bar{L}_{1}^{2PN}$ | $-480\left(\eta^{2}-8\eta+15\right)$
$\bar{L}_{0}^{2PN}$ | $120\left(65\eta^{2}-238\eta+180\right)$
$\bar{K}_{8}^{2PN}$ | $15\eta(29\eta-3)$
$\bar{K}_{7}^{2PN}$ | $-60\left(129\eta^{2}-188\eta+74\right)$
$\bar{K}_{6}^{2PN}$ | $-15\left(116\eta^{2}-711\eta+304\right)$
$\bar{K}_{5}^{2PN}$ | $2\left(5516\eta^{2}-15155\eta+5420\right)$
$\bar{K}_{4}^{2PN}$ | $-4\left(5347\eta^{2}-13720\eta+6400\right)$
$\bar{K}_{3}^{2PN}$ | $8\left(\eta^{2}+2735\eta-4265\right)$
$\bar{K}_{2}^{2PN}$ | $20308\eta^{2}-81610\eta+71380$
$\bar{K}_{1}^{2PN}$ | $-720\left(22\eta^{2}-82\eta+65\right)$
$\bar{K}_{0}^{2PN}$ | $60\left(65\eta^{2}-238\eta+180\right)$
$\bar{L}^{SS}$ | $4e_{r0}\left(4e_{r0}^{4}+29e_{r0}^{3}\right.$
| $\left.+30e_{r0}^{2}+48e_{r0}+24\right)$
$\bar{K}^{SS}$ | $\cos\left(\zeta_{1}+\zeta_{2}\right)\left[96e_{r0}^{3}-236e_{r0}^{2}\right.$
| $-171e_{r0}^{4}+95e_{r0}^{5}+56e_{r0}-32$
|
$\left.-32\sqrt{1-e_{r0}^{2}}\left(e_{r0}^{2}-e_{r0}^{3}+e_{r0}-1\right)\right]$
|
$-2e_{r0}\cos\left(\zeta_{1}-\zeta_{2}\right)\left(4e_{r0}^{4}+29e_{r0}^{3}\right.$
| $\left.+30e_{r0}^{2}+48e_{r0}+24\right)$

The integral of Eq. (99) results in

$\bar{e}_{rN}=\frac{2\pi
e_{r0}\mathfrak{l}_{r0}^{3}}{\left(1-e_{r0}^{2}\right)^{3/2}}\leavevmode\nobreak\
,$ (108) $\displaystyle\bar{e}_{rPN}$ $\displaystyle=$
$\displaystyle\frac{\pi\mathfrak{l}_{r0}}{e_{r0}}\left\\{2(3\eta-5)\right.$
(109)
$\displaystyle+\frac{\sqrt{1-e_{r0}^{2}}}{(1-e_{r0})^{3}(e_{r0}+1)}\left[4e_{r0}^{4}\left(\eta-2\right)\right.$
$\displaystyle-20e_{r0}^{3}\left(\eta-1\right)+e_{r0}^{2}\left(22-9\eta\right)$
$\displaystyle\left.\left.+2e_{r0}\left(5\eta-7\right)-6\eta+10\right]\right\\}\leavevmode\nobreak\
,$ $\displaystyle\bar{e}_{rSO}$ $\displaystyle=$
$\displaystyle\frac{\pi(e_{r0}+1)\eta}{\sqrt{1-e_{r0}^{2}}}$ (110)
$\displaystyle\times\sum_{k=1}^{2}\left(4^{2k-3}+3\right)\chi_{k}\cos\kappa_{k}\leavevmode\nobreak\
,$ $\displaystyle\bar{e}_{r2PN}$ $\displaystyle=$
$\displaystyle\frac{\pi}{480e_{r0}^{3}\mathfrak{l}_{r0}}\sum_{k=0}^{4}\bar{E}_{k}^{2PN}e_{r0}^{k}$
(111)
$\displaystyle-\frac{\pi\sqrt{1-e_{r0}^{2}}}{60(e_{r0}-1)^{4}e_{r0}^{3}\mathfrak{l}_{r0}}\sum_{k=0}^{8}\bar{F}_{k}^{2PN}e_{r0}^{k}\leavevmode\nobreak\
,$ $\displaystyle\bar{e}_{rSS}$ $\displaystyle=$
$\displaystyle\frac{3\pi\chi_{1}\chi_{2}\eta}{16\left(1-e_{r0}\right)e_{r0}^{3}\left(1-e_{r0}^{2}\right)^{2}\mathfrak{l}_{r0}}$
(112)
$\displaystyle\times\left(\cos\kappa_{1}\cos\kappa_{2}\bar{E}^{SS}\right.$
$\displaystyle\left.+\sin\kappa_{1}\sin\kappa_{2}\bar{F}^{SS}\right)\leavevmode\nobreak\
,$ $\displaystyle\bar{e}_{rQM}$ $\displaystyle=$
$\displaystyle\frac{\pi\left(e_{r0}^{2}+2\right)\eta}{128\left(1-e_{r0}^{2}\right)^{5/2}\mathfrak{l}_{r0}}\sum_{k=1}^{2}\chi_{k}^{2}\nu^{2k-3}w_{k}$
(113)
$\displaystyle\times\left[-4\left(10e_{r0}^{3}+609e_{r0}^{2}+260e_{r0}\right.\right.$
$\displaystyle\left.+336\right)\sin^{2}\kappa_{k}\cos
2\zeta_{k}+8\left(4e_{r0}^{3}\right.$
$\displaystyle\left.+27e_{r0}^{2}+15e_{r0}+12\right)$
$\displaystyle\left.\times\left(3\cos
2\kappa_{k}+1\right)\right]\leavevmode\nobreak\ .$

The coefficients $\bar{E}_{k}^{2PN}$, $\bar{F}_{k}^{2PN}$, $\bar{E}^{SS}$, and
$\bar{F}^{SS}$ are collected in Table 6.

Table 6: The coefficients of $\bar{e}_{r}$ in Eqs. (111) and (112). $Coefficient$ | $Expression$
---|---
$\bar{E}_{4}^{2PN}$ | $15\left(1111\eta^{2}-1624\eta+528\right)$
$\bar{E}_{3}^{2PN}$ | $4800\left(4\eta^{2}-7\eta+4\right)$
$\bar{E}_{2}^{2PN}$ | $-4\left(151\eta^{2}-7550\eta+5640\right)$
$\bar{E}_{1}^{2PN}$ | $2880(\eta-3)\eta$
$\bar{E}_{0}^{2PN}$ | $8\left(1501\eta^{2}-8090\eta+7260\right)$
$\bar{F}_{8}^{2PN}$ | $120(\eta-3)\eta$
$\bar{F}_{7}^{2PN}$ | $-60\left(44\eta^{2}-107\eta+44\right)$
$\bar{F}_{6}^{2PN}$ | $-15\left(317\eta^{2}-537\eta+96\right)$
$\bar{F}_{5}^{2PN}$ | $4\left(457\eta^{2}-1430\eta-960\right)$
$\bar{F}_{4}^{2PN}$ | $-5362\eta^{2}+18335\eta-4380$
$\bar{F}_{3}^{2PN}$ | $-6\left(714\eta^{2}-4315\eta+5420\right)$
$\bar{F}_{2}^{2PN}$ | $21\left(391\eta^{2}-2110\eta+2100\right)$
$\bar{F}_{1}^{2PN}$ | $-4\left(1411\eta^{2}-7820\eta+7260\right)$
$\bar{F}_{0}^{2PN}$ | $1501\eta^{2}-8090\eta+7260$
$\bar{E}^{SS}$ | $16e_{r0}^{3}(1+e_{r0})^{2}\sqrt{1-e_{r0}^{2}}\left(e_{r0}^{2}+2\right)$
$\bar{F}^{SS}$ | $-\left[8(e_{r0}+1)^{2}\sqrt{1-e_{r0}^{2}}\right.$
| $\times\left(e_{r0}^{2}+2\right)e_{r0}^{3}\cos(\zeta_{1}-\zeta_{2})$
|
$+\cos\left(\zeta_{1}+\zeta_{2}\right)\left[3176e_{r0}^{3}-3176e_{r0}^{2}\right.$
| $+1552e_{r0}^{4}-1552e_{r0}^{5}-35e_{r0}^{6}$
| $+\sqrt{1-e_{r0}^{2}}\left(2376e_{r0}^{2}-2328e_{r0}^{3}\right.$
| $-468e_{r0}^{4}+606e_{r0}^{5}-35e_{r0}^{6}$
| $\left.\left.+92e_{r0}^{7}+1600e_{r0}-1600\right)\right]$

###  A.5 Initial shape variables in terms of the secular shape variables

The contributions to the averaged shape variables in terms of
$\mathfrak{l}_{r0}$ and $e_{r0}$, Eqs. (100) and (101), were presented in the
previous subsection. Here, we invert these relations to generate
$\mathfrak{l}_{r0}$ and $e_{r0}$ in terms of $\mathfrak{\bar{l}}_{r}$ and
$\bar{e}_{r}$.

We do this in two steps. First, we take the perturbations to linear order.
Using Eqs. (100), (102), and (93) for $\mathfrak{l}_{r0}$, we find

$\displaystyle\mathfrak{l}_{r0}$ $\displaystyle=$
$\displaystyle\mathfrak{\bar{l}}_{r}\frac{\mathfrak{T}}{\mathfrak{T}_{0}}-\frac{\mathfrak{1}}{\mathfrak{T}_{0}}\left(\mathfrak{\bar{l}}_{rPN}+\mathfrak{\bar{l}}_{rSO}\right.$
(114)
$\displaystyle\left.+\mathfrak{\bar{l}}_{rSS}+\mathfrak{\bar{l}}_{rQM}\right)\leavevmode\nobreak\
.$

Using Eqs (101), (108), and (93) for $e_{r0}$, we get

$\displaystyle e_{r0}$ $\displaystyle=$
$\displaystyle\bar{e}_{r}\frac{\mathfrak{T}}{\mathfrak{T}_{0}}-\frac{\mathfrak{1}}{\mathfrak{T}_{0}}\left(\bar{e}_{rPN}+\bar{e}_{rSO}\right.$
(115)
$\displaystyle\left.+\bar{e}_{rSS}+\bar{e}_{rQM}\right)\leavevmode\nobreak\ .$

In the perturbation terms, we can insert the leading-order terms of
$\mathfrak{l}_{r0}$ and $e_{r0}$, which are

$\displaystyle\mathfrak{l}_{r0}$ $\displaystyle=$
$\displaystyle\mathfrak{\bar{l}}_{r}\leavevmode\nobreak\ ,$ (116)
$\displaystyle e_{r0}$ $\displaystyle=$
$\displaystyle\bar{e}_{r}\leavevmode\nobreak\ .$ (117)

The results are

$\mathfrak{l}_{r0PN}=\frac{2\bar{e}_{r}(\bar{e}_{r}+1)(\eta-2)}{\mathfrak{\bar{l}}_{r}}\leavevmode\nobreak\
,$ (118) $\displaystyle\mathfrak{l}_{r0SO}$ $\displaystyle=$
$\displaystyle-\frac{\bar{e}_{r}\left(1-\bar{e}_{r}^{2}\right)\eta}{2(\bar{e}_{r}-1)\mathfrak{\bar{l}}_{r}^{2}}$
(119)
$\displaystyle\times\sum_{k=1}^{2}\left(4^{2k-3}+3\right)\chi_{k}\cos\kappa_{k}\leavevmode\nobreak\
,$ $\displaystyle\mathfrak{l}_{r0SS}$ $\displaystyle=$
$\displaystyle-\frac{\chi_{1}\chi_{2}\eta}{32(1-\bar{e}_{r})\bar{e}_{r}^{2}\mathfrak{\bar{l}}_{r}^{3}}\left(L_{0}^{SS}\cos\kappa_{1}\cos\kappa_{2}\right.$
(120)
$\displaystyle\left.-K_{0}^{SS}\sin\kappa_{1}\sin\kappa_{2}\right)\leavevmode\nobreak\
,$ $\displaystyle\mathfrak{l}_{r0QM}$ $\displaystyle=$
$\displaystyle\frac{\eta}{512\bar{e}_{r}\left(\bar{e}_{r}^{2}-1\right)\mathfrak{\bar{l}}_{r}^{3}}\sum_{k=1}^{2}\chi_{k}^{2}\nu^{2k-3}w_{k}$
(121)
$\displaystyle\times\left[-4\left(47\bar{e}_{r}^{5}+1338\bar{e}_{r}^{4}+1446\bar{e}_{r}^{3}\right.\right.$
$\displaystyle\left.+3444\bar{e}_{r}^{2}+1264\bar{e}_{r}+960\right)$
$\displaystyle\times\sin^{2}\kappa_{k}\cos
2\zeta_{k}+16\left(5\bar{e}_{r}^{5}+33\bar{e}_{r}^{4}\right.$
$\displaystyle\left.+61\bar{e}_{r}^{3}+84\bar{e}_{r}^{2}+42\bar{e}_{r}+12\right)$
$\displaystyle\left.\times\left(3\cos 2\kappa_{k}+\allowbreak
1\right)\right]\leavevmode\nobreak\ ,$

where the coefficients $L_{0}^{SS}$ and $K_{0}^{SS}$ are listed in Table 7.
For $e_{r0}$, we have

$\displaystyle e_{r0PN}$ $\displaystyle=$
$\displaystyle\frac{1}{\bar{e}_{r}\mathfrak{\bar{l}}_{r}^{2}}\left\\{-\left(1-\bar{e}_{r}^{2}\right)^{3/2}\left(3\eta-5\right)\right.$
(122)
$\displaystyle+\frac{\left(1+\bar{e}_{r}\right)^{2}}{2}\left[\bar{e}_{r}^{2}\left(11\eta-14\right)\right.$
$\displaystyle\left.\left.-2\bar{e}_{r}\left(5\eta-7\right)+6\eta-10\right]\right\\}\leavevmode\nobreak\
,$ $\displaystyle e_{r0SO}$ $\displaystyle=$
$\displaystyle-\frac{(\bar{e}_{r}+1)\left(1-\bar{e}_{r}^{2}\right)\eta}{2\mathfrak{\bar{l}}_{r}^{3}}$
(123)
$\displaystyle\times\sum_{k=1}^{2}\left(4^{2k-3}+3\right)\chi_{k}\cos\kappa_{k}\leavevmode\nobreak\
,$ $\displaystyle e_{r0SS}$ $\displaystyle=$
$\displaystyle-\frac{3\chi_{1}\chi_{2}\eta}{32(1-\bar{e}_{r})\bar{e}_{r}^{3}\mathfrak{\bar{l}}_{r}^{4}}\left(E_{0}^{SS}\cos\kappa_{1}\cos\kappa_{2}\right.$
(124)
$\displaystyle\left.-F_{0}^{SS}\sin\kappa_{1}\sin\kappa_{2}\right)\leavevmode\nobreak\
,$ $\displaystyle e_{r0QM}$ $\displaystyle=$
$\displaystyle\frac{\eta}{256\left(\bar{e}_{r}^{2}-1\right)\mathfrak{\bar{l}}_{r}^{4}}\sum_{k=1}^{2}\chi_{k}^{2}\nu^{2k-3}w_{k}$
(125)
$\displaystyle\times\left[-4\left(10\bar{e}_{r}^{5}+753\bar{e}_{r}^{4}+712\bar{e}_{r}^{3}\right.\right.$
$\displaystyle\left.+1986\bar{e}_{r}^{2}+664\bar{e}_{r}+672\right)$
$\displaystyle\times\sin^{2}\kappa_{k}\cos
2\zeta_{k}+8\left(4\bar{e}_{r}^{5}+39\bar{e}_{r}^{4}\right.$
$\displaystyle+59\bar{e}_{r}^{3}+102\bar{e}_{r}^{2}+42\bar{e}_{r}$
$\displaystyle\left.\left.+24\right)\left(3\cos 2\kappa_{k}+\allowbreak
1\right)\right]\leavevmode\nobreak\ .$

where the coefficients $E_{0}^{SS}$ and $F_{0}^{SS}$ can be found in Table 8.

Table 7: The coefficients in Eqs. (120) and (128). $Coefficient$ | $Expression$
---|---
$L_{0,2}^{2PN}$ | $3\left(75\eta^{2}-176\eta+216\right)$
$L_{0,1}^{2PN}$ | $48\left(3\eta^{2}-11\eta+10\right)$
$L_{0,0}^{2PN}$ | $86\eta^{2}-260\eta+176$
$K_{0,4}^{2PN}$ | $12\left(12\eta^{2}-39\eta+28\right)$
$K_{0,3}^{2PN}$ | $-6\left(18\eta^{2}-63\eta+64\right)$
$K_{0,2}^{2PN}$ | $3\left(17\eta^{2}-7\eta+28\right)$
$K_{0,1}^{2PN}$ | $-2\left(7\eta^{2}+2\eta-32\right)$
$K_{0,0}^{2PN}$ | $43\eta^{2}-130\eta+88$
$L_{0}^{SS}$ | $4\bar{e}_{r}\left(4\bar{e}_{r}^{4}+53\bar{e}_{r}^{3}\right.$
| $\left.+78\bar{e}_{r}^{2}+72\bar{e}_{r}+24\right)$
$K_{0}^{SS}$ | $\left[-95\bar{e}_{r}^{5}+171\bar{e}_{r}^{4}-96\bar{e}_{r}^{3}\right.$
| $+236\bar{e}_{r}^{2}-56\bar{e}_{r}+32$
|
$\left.-32\left(1-\bar{e}_{r}\right)\left(1-\bar{e}_{r}^{2}\right)^{3/2}\right]$
| $\times\cos(\zeta_{1}+\zeta_{2})+96\bar{e}_{r}^{2}$
|
$\times\left(1+\bar{e}_{r}\right)^{2}\left(2\cos\zeta_{1}\cos\zeta_{2}\right.$
|
$\left.-\sin\zeta_{1}\sin\zeta_{2}\right)+2\bar{e}_{r}\left(4\bar{e}_{r}^{4}\right.$
| $\left.+29\bar{e}_{r}^{3}+30\bar{e}_{r}^{2}+48\bar{e}_{r}+24\right)$
| $\times\cos(\zeta_{1}-\zeta_{2})$
Table 8: The coefficients in Eqs. (124) and (129). $Coefficient$ | $Expression$
---|---
$E_{0,4}^{2PN}$ | $15\left(2915\eta^{2}-8904\eta+6192\right)$
$E_{0,3}^{2PN}$ | $960\left(36\eta^{2}-102\eta+65\right)$
$E_{0,2}^{2PN}$ | $4\left(7673\eta^{2}-20110\eta+15360\right)$
$E_{0,1}^{2PN}$ | $960\left(2\eta^{2}-11\eta+15\right)$
$E_{0,0}^{2PN}$ | $-8\left(4559\eta^{2}-13390\eta+9540\right)$
$F_{0,6}^{2PN}$ | $15\left(383\eta^{2}-989\eta+560\right)$
$F_{0,5}^{2PN}$ | $-30\left(51\eta^{2}-261\eta+224\right)$
$F_{0,4}^{2PN}$ | $-15\left(68\eta^{2}-301\eta+232\right)$
$F_{0,3}^{2PN}$ | $2\left(893\eta^{2}-3505\eta+3360\right)$
$F_{0,2}^{2PN}$ | $-2\left(1756\eta^{2}-6425\eta+5250\right)$
$F_{0,1}^{2PN}$ | $9358\eta^{2}-28100\eta+20880$
$F_{0,0}^{2PN}$ | $-4559\eta^{2}+13390\eta-9540$
$E_{0}^{SS}$ | $16\bar{e}_{r}^{3}\left(1+\bar{e}_{r}\right)\left(\bar{e}_{r}^{3}+3\bar{e}_{r}^{2}+4\bar{e}_{r}+2\right)$
$F_{0}^{SS}$ | $\cos\left(\zeta_{1}+\zeta_{2}\right)\left[92\bar{e}_{r}^{7}-35\bar{e}_{r}^{6}\right.$
| $+606\bar{e}_{r}^{5}-468\bar{e}_{r}^{4}-2328\bar{e}_{r}^{3}$
| $+2376\bar{e}_{r}^{2}+1600\bar{e}_{r}-1600$
|
$\left.+8\left(\bar{e}_{r}+1\right)\left(3\bar{e}_{r}^{2}+200\right)\left(1-\bar{e}_{r}\right)^{5/2}\right]$
| $+32(\bar{e}_{r}+1)^{2}\bar{e}_{r}^{4}$
| $\times\left(2\cos\zeta_{1}\cos\zeta_{2}-\sin\zeta_{1}\sin\zeta_{2}\right)$
| $+8(\bar{e}_{r}+1)^{2}\left(\bar{e}_{r}^{2}+2\right)\bar{e}_{r}^{3}$
| $\times\cos(\zeta_{1}-\zeta_{2})$

The second step is to derive the 2PN terms. For this, we use Eqs. (100),
(102), and (93) for $\mathfrak{l}_{r0}$ and Eqs. (101), (108), and (93) for
$e_{r0}$:

$\displaystyle\mathfrak{l}_{r0}$ $\displaystyle=$
$\displaystyle\mathfrak{\bar{l}}_{r}\frac{\mathfrak{T}}{\mathfrak{T}_{0}}-\frac{\mathfrak{1}}{\mathfrak{T}_{0}}\left(\mathfrak{\bar{l}}_{rPN}+\mathfrak{\bar{l}}_{r2PN}\right)\leavevmode\nobreak\
,$ (126) $\displaystyle e_{r0}$ $\displaystyle=$
$\displaystyle\bar{e}_{r}\frac{\mathfrak{T}}{\mathfrak{T}_{0}}-\frac{\mathfrak{1}}{\mathfrak{T}_{0}}\left(\bar{e}_{rPN}+\bar{e}_{r2PN}\right)\leavevmode\nobreak\
.$ (127)

This time, in order to get the 2PN terms, we need the previously calculated
1PN expressions of $\mathfrak{l}_{r0}$ and $e_{r0}$. After the Taylor
expansion to 2PN order, we find

$\displaystyle\mathfrak{l}_{r02PN}$ $\displaystyle=$
$\displaystyle-\frac{(1-\bar{e}_{r}^{2})^{3/2}}{24\bar{e}_{r}^{2}\mathfrak{\bar{l}}_{r}^{3}}\sum_{k=0}^{2}L_{0,k}^{2PN}\bar{e}_{r}^{k}$
(128)
$\displaystyle+\frac{(\bar{e}_{r}+1)^{2}}{12\bar{e}_{r}^{2}\mathfrak{\bar{l}}_{r}^{3}}\sum_{k=0}^{4}K_{0,k}^{2PN}\bar{e}_{r}^{k}\leavevmode\nobreak\
,$ $\displaystyle e_{r02PN}$ $\displaystyle=$
$\displaystyle-\frac{\left(1-\bar{e}_{r}^{2}\right)^{3/2}}{960\bar{e}_{r}^{3}\mathfrak{\bar{l}}_{r}^{4}}\sum_{k=0}^{4}E_{0,k}^{2PN}\bar{e}_{r}^{k}$
(129)
$\displaystyle+\frac{(1+\bar{e}_{r})^{2}}{120\bar{e}_{r}^{3}\mathfrak{\bar{l}}_{r}^{4}}\sum_{k=0}^{6}F_{0,k}^{2PN}\bar{e}_{r}^{k}\leavevmode\nobreak\
.$

The coefficients $L_{0,k}^{2PN}$ and $K_{0,k}^{2PN}$ are given in Table 7,
while the terms $E_{0,k}^{2PN}$ and $F_{0,k}^{2PN}$ can be found in Table 8.

The full expressions of $\mathfrak{l}_{r0}$ and $e_{r0}$ are the sum of the
corresponding above contributions Eqs. (116), (118)–(121), (128) and Eqs.
(117), (122)–(125), (129), respectively.

###  A.6 Radial period in terms of time averages of shape variables

Here, we finally are able to express the radial period in terms of averaged
quantities by replacing the initial values at the periastron with time
averages over $\chi_{p}\in\left[0,2\pi\right]$. The various order
contributions to Eq. (6) become

$\displaystyle\mathfrak{\tilde{T}}$ $\displaystyle=$
$\displaystyle\frac{2\mathfrak{\bar{l}}_{r}^{3}\pi}{\left(1-\bar{e}_{r}^{2}\right)^{3/2}}\leavevmode\nobreak\
,$ (130) $\displaystyle\tilde{\tau}_{PN}$ $\displaystyle=$
$\displaystyle\sqrt{1-\bar{e}_{r}^{2}}(15-9\eta)$ (131)
$\displaystyle+\left(1-\bar{e}_{r}^{2}\right)(7\eta-6)\leavevmode\nobreak\ ,$
$\displaystyle\tilde{\tau}_{SO}$ $\displaystyle=$ $\displaystyle
0\leavevmode\nobreak\ ,$ (132) $\displaystyle\tilde{\tau}_{2PN}$
$\displaystyle=$
$\displaystyle\frac{\sqrt{1-\bar{e}_{r}^{2}}}{64\bar{e}_{r}^{4}}\sum_{k=0}^{6}\bar{U}_{k}\bar{e}_{r}^{k}$
(133)
$\displaystyle-\frac{(\bar{e}_{r}+1)}{8\bar{e}_{r}^{4}}\sum_{k=0}^{7}\bar{V}_{k}\bar{e}_{r}^{k}\leavevmode\nobreak\
,$ $\displaystyle\tilde{\tau}_{QM}$ $\displaystyle=$
$\displaystyle\frac{3\eta}{512\bar{e}_{r}\left(1-\bar{e}_{r}^{2}\right)^{2}}\sum_{k=1}^{2}\chi_{k}^{2}\nu^{2k-3}w_{k}$
(134) $\displaystyle\times\left[\bar{U}^{QM}\sin^{2}\kappa_{k}\cos
2\zeta_{k}\right.$ $\displaystyle\left.+\bar{V}^{QM}\left(3\cos
2\kappa_{k}\allowbreak+1\right)\right]\leavevmode\nobreak\ ,$
$\displaystyle\bar{\tau}_{SS}$ $\displaystyle=$
$\displaystyle\frac{3\chi_{1}\chi_{2}\eta}{8(1-\bar{e}_{r})^{2}\bar{e}_{r}}\left(\cos\kappa_{1}\cos\kappa_{2}\bar{U}^{SS}\right.$
(135)
$\displaystyle\left.+\sin\kappa_{1}\sin\kappa_{2}\bar{V}^{SS}\right)\leavevmode\nobreak\
.$

The coefficients in the above expressions are listed in Table 9.

Table 9: The coefficients in Eqs. (133)–(135). $Coefficient$ | $Expression$
---|---
$\bar{U}_{6}$ | $-437\eta^{2}+3336\eta-1008$
$\bar{U}_{5}$ | $-64\left(8\eta^{2}-6\eta-5\right)$
$\bar{U}_{4}$ | $-8\left(211\eta^{2}-159\eta+336\right)$
$\bar{U}_{3}$ | $64\left(4\eta^{2}+11\eta-5\right)$
$\bar{U}_{2}$ | $-8\left(79\eta^{2}-600\eta+528\right)$
$\bar{U}_{1}$ | $-128\left(\eta^{2}-8\eta+15\right)$
$\bar{U}_{0}$ | $32\left(65\eta^{2}-238\eta+180\right)$
$\bar{V}_{7}$ | $224\eta^{2}-690\eta+360$
$\bar{V}_{6}$ | $2\left(64\eta^{2}-11\eta-12\right)$
$\bar{V}_{5}$ | $139\eta^{2}-410\eta+452$
$\bar{V}_{4}$ | $-179\eta^{2}+266\eta-308$
$\bar{V}_{3}$ | $-27\eta^{2}+28\eta+8$
$\bar{V}_{2}$ | $67\eta^{2}-4\eta+72$
$\bar{V}_{1}$ | $-12\left(23\eta^{2}-90\eta+80\right)$
$\bar{V}_{0}$ | $260\eta^{2}-952\eta+720$
$\bar{U}^{QM}$ | $-4\left(27\bar{e}_{r}^{7}-72\bar{e}_{r}^{6}+263\bar{e}_{r}^{5}\right.$
| $-1674\bar{e}_{r}^{4}-1702\bar{e}_{r}^{3}-4116\bar{e}_{r}^{2}$
| $\left.-1360\bar{e}_{r}-960\right)$
$\bar{V}^{QM}$ | $16\left(\bar{e}_{r}^{7}-2\bar{e}_{r}^{6}+9\bar{e}_{r}^{5}-43\bar{e}_{r}^{4}\right.$
| $\left.-69\bar{e}_{r}^{3}-108\bar{e}_{r}^{2}-46\bar{e}_{r}-12\right)$
$\bar{U}^{SS}$ | $-8\bar{e}_{r}^{5}+21\bar{e}_{r}^{4}-15\bar{e}_{r}^{3}$
| $-38\bar{e}_{r}^{2}-56\bar{e}_{r}-24$
$\bar{V}^{SS}$ | $\frac{\cos(\zeta_{1}-\zeta_{2})}{2}\allowbreak\left(8\bar{e}_{r}^{5}-21\bar{e}_{r}^{4}\right.$
| $\left.+15\bar{e}_{r}^{3}+38\bar{e}_{r}^{2}+56\bar{e}_{r}+24\right)$
|
$+\frac{\cos(\zeta_{1}+\zeta_{2})}{4\bar{e}_{r}^{2}(\bar{e}_{r}+1)}\left[371\bar{e}_{r}^{7}-276\bar{e}_{r}^{6}\right.$
|
$+\bar{e}_{r}^{5}\left(-48\bar{e}_{r}+104\sqrt{1-\bar{e}_{r}^{2}}+1771\right)$
|
$-\bar{e}_{r}^{4}\left(48\bar{e}_{r}+104\sqrt{1-\bar{e}_{r}^{2}}+1517\right)$
|
$+8\bar{e}_{r}^{3}\left(24\bar{e}_{r}+583\sqrt{1-\bar{e}_{r}^{2}}-854\right)$
|
$+4\bar{e}_{r}^{2}\left(48\bar{e}_{r}-1166\sqrt{-\bar{e}_{r}^{2}+1}+1881\right)$
| $-8\left(596\sqrt{1-\bar{e}_{r}^{2}}-593\right)\bar{e}_{r}$
| $\left.+4768\left(\sqrt{1-\bar{e}_{r}^{2}}-1\right)\right]$

###  A.7 Expansion of the averaged PN parameter

The PN parameter associated to the averaged dynamics is given by

$\bar{\varepsilon}=\bar{\varepsilon}_{0PN}+\bar{\varepsilon}_{1PN}+\bar{\varepsilon}_{2PN}+\bar{\varepsilon}_{SO}+\bar{\varepsilon}_{SS}+\bar{\varepsilon}_{QM}\leavevmode\nobreak\
,$ (136)

with

$\bar{\varepsilon}_{0PN}=\frac{2\mathfrak{\bar{l}}_{r}\pi}{\mathfrak{T}\left(1-\bar{e}_{r}^{2}\right)^{1/2}}\leavevmode\nobreak\
,$ (137)
$\bar{\varepsilon}_{1PN}=-\frac{\pi}{\mathfrak{T\bar{l}}_{r}}\left[\left(6-7\eta\right)\left(1-\bar{e}_{r}^{2}\right)^{1/2}-9\left(2-\eta\right)\right]\leavevmode\nobreak\
,$ (138) $\displaystyle\bar{\varepsilon}_{SO}$ $\displaystyle=$
$\displaystyle\frac{3\pi}{2\mathfrak{T\bar{l}}_{r}^{2}}\left(1-\frac{2}{3}\sqrt{\frac{1+\bar{e}_{r}}{1-\bar{e}_{r}}}\right)$
(139)
$\displaystyle\times\sum_{i=1}^{2}\chi_{i}\cos\kappa_{i}\left(4\nu^{2i-3}+3\right)\leavevmode\nobreak\
,$ $\displaystyle\bar{\varepsilon}_{SS}$ $\displaystyle=$
$\displaystyle\frac{\pi\chi_{1}\chi_{2}\eta}{16\bar{e}_{r}^{2}(\bar{e}_{r}-1)^{3}(\bar{e}_{r}+1)^{2}\mathfrak{\bar{l}}_{r}^{3}\mathfrak{T}}$
(140)
$\displaystyle\times\left\\{16\bar{e}_{r}\left(\bar{e}_{r}+1\right)\cos\kappa_{1}\cos\kappa_{2}M^{SS}\right.$
$\displaystyle\left.+\sin(\kappa_{1})\sin(\kappa_{2})N^{SS}\right\\}$
$\displaystyle\bar{\varepsilon}_{QM}$ $\displaystyle=$
$\displaystyle\frac{\pi\eta}{256\bar{e}_{r}\left(1-\bar{e}_{r}^{2}\right)^{5/2}\mathfrak{\bar{l}}_{r}^{3}\mathfrak{T}}\sum_{i=1}^{2}\nu^{2i-3}\chi_{i}^{2}w_{i}$
(141)
$\displaystyle\times\left\\{M^{QM}\left[\cos\left(2\kappa_{i}-2\zeta_{i}\right)+\cos\left(2\kappa_{i}+2\zeta_{i}\right)\right.\right.$
$\displaystyle\left.\left.-2\cos 2\zeta_{i}\right]+N^{QM}\right\\}$

The coefficients of $\bar{\varepsilon}_{SS}$ and $\bar{\varepsilon}_{QM}$ are
given in Table 10.

Table 10: The coefficients $M^{SS}$,$\leavevmode\nobreak\ M^{QM}$,$\leavevmode\nobreak\ N^{SS}$, and $N^{QM}$ of $\mathfrak{\bar{\varepsilon}}$ in Eqs. (140) and (141). $Coefficient$ | $Expression$
---|---
$M^{SS}$ | $-2\bar{e}_{r}+4\bar{e}_{r}^{2}-4\bar{e}_{r}^{4}+2\bar{e}_{r}^{5}$
|
$-\left(16\bar{e}_{r}^{2}+190+96-19\bar{e}_{r}^{5}-203\bar{e}_{r}^{4}\right.$
| $\left.-110\bar{e}_{r}^{3}\right)\sqrt{1-\bar{e}_{r}^{2}}$
$N^{SS}$ | $\cos\left(\zeta_{1}+\zeta_{2}\right)\left\\{\allowbreak 512+\allowbreak 2048\bar{e}_{r}^{2}\right.$
| $-\allowbreak 3072\bar{e}_{r}^{4}+2048\bar{e}_{r}^{6}-512\bar{e}_{r}^{8}$
| $+\left[5312-\allowbreak 5696\bar{e}_{r}-\allowbreak
1608\bar{e}_{r}^{2}\right.$
| $+\allowbreak 6248\bar{e}_{r}^{3}+\allowbreak
220\bar{e}_{r}^{4}\allowbreak-1994\bar{e}_{r}^{5}$
| $-5031\bar{e}_{r}^{6}+1244\bar{e}_{r}^{7}$
| $+\left(-4800+\allowbreak 4728\bar{e}_{r}^{2}+72\bar{e}_{r}^{4}\right)$
|
$\left.\left.\times\left(1-\bar{e}_{r}\right)^{3/2}\right]\sqrt{1-\bar{e}_{r}^{2}}\right\\}$
| $+8\bar{e}_{r}(\bar{e}_{r}+1)\cos\left(\zeta_{1}-\zeta_{2}\right)$
|
$\left[-2\bar{e}_{r}\allowbreak-4\bar{e}_{r}^{2}-4\bar{e}_{r}^{4}+2\bar{e}_{r}^{5}\right.$
| $\left(-96-190\bar{e}_{r}-176\bar{e}_{r}^{2}-82\bar{e}_{r}^{3}\right.$
|
$\left.\left.+203\bar{e}_{r}^{4}+19\bar{e}_{r}^{5}\right)\sqrt{1-\bar{e}_{r}^{2}}\right]$
$M^{QM}$ | $27\bar{e}_{r}^{7}-40\bar{e}_{r}^{6}+167\bar{e}_{r}^{5}-2122\bar{e}_{r}^{4}$
| $-1990\bar{e}_{r}^{3}-3700\bar{e}_{r}^{2}-976\bar{e}_{r}-960$
$N^{QM}$ | $16\bar{e}_{r}^{7}+416\bar{e}_{r}^{6}+1104\bar{e}_{r}^{5}-816\bar{e}_{r}^{4}$
| $-2064\bar{e}_{r}^{3}-2048\bar{e}_{r}^{2}-736\bar{e}_{r}-192$
|
$+\left(64\bar{e}_{r}+64\bar{e}_{r}^{5}-128\bar{e}_{r}^{3}\right)\sqrt{1-\bar{e}_{r}^{2}}$
| $-256\left(7\bar{e}_{r}^{4}+15\bar{e}_{r}^{3}-2\bar{e}_{r}^{2}\right.$
|
$\left.-15\bar{e}_{r}-5\right)\bar{e}_{r}^{2}\sin^{2}\kappa_{i}\cos^{2}\zeta_{i}$
|
$+16\left\\{3\bar{e}_{r}^{7}-34\bar{e}_{r}^{6}-121\bar{e}_{r}^{4}-304\bar{e}_{r}^{2}\right.$
| $-138\bar{e}_{r}-33\bar{e}_{r}^{5}-\allowbreak 147\bar{e}_{r}^{3}-36$
|
$\left.+\left(-24+12\bar{e}_{r}+12\bar{e}_{r}^{5}\right)\sqrt{1-\bar{e}_{r}^{2}}\right]\cos
2\kappa_{i}$

$\displaystyle\bar{\varepsilon}_{2PN}$ $\displaystyle=$
$\displaystyle\frac{-\pi}{4915200\left(1-\bar{e}_{r}^{2}\right)^{5/2}\bar{e}_{r}^{4}\mathfrak{\bar{l}}^{3}\mathfrak{T}}$
(142)
$\displaystyle\times\left\\{15\bar{e}_{r}^{5}M_{1}^{2PN}-2457600\left(1+\bar{e}_{r}\right)^{3}\bar{e}_{r}^{3}M_{2}^{2PN}\right.$
$\displaystyle-10240\bar{e}_{r}^{2}\left(1+\bar{e}_{r}\right)^{2}M_{3}^{2PN}$
$\displaystyle+2\left(1+\bar{e}_{r}\right)\left[8\left(1-\bar{e}_{r}-\bar{e}_{r}^{2}+\bar{e}_{r}^{3}\right)\sqrt{1-\bar{e}_{r}^{2}}\right.$
$\displaystyle\left.-8+8\bar{e}_{r}+12\bar{e}_{r}^{2}-14\bar{e}_{r}^{3}-\bar{e}_{r}^{4}\right]M_{4}^{2PN}$
$\displaystyle+5\left[15\bar{e}_{r}^{4}-20\bar{e}_{r}^{2}-8\left(1-\bar{e}_{r}^{2}\right)^{5/2}\right.$
$\displaystyle\left.\left.+8\right]M_{5}^{2PN}\right\\}$

The coefficents of $\bar{\varepsilon}_{2PN}$ are enlisted in Table 11.

Table 11: The coefficients $M_{k}^{2PN}$ of $\mathfrak{\bar{\varepsilon}}$ in Eq. (142). $Coefficient$ | $Expression$
---|---
$M_{1}^{2PN}$ | $-15\eta(98967\eta+118819)\bar{e}_{r}^{6}$
| $+1920(\eta(8399\eta-17134)+12016)\bar{e}_{r}^{5}$
| $+4(\eta(5837763\eta-15119020)+8076560)\bar{e}_{r}^{4}$
| $-5120(\eta(2431\eta-13254)+662)\bar{e}_{r}^{3}$
| $+160(41\eta(3071\eta-11629)+502740)\bar{e}_{r}^{2}$
| $+40960(\eta(68\eta+1067)-408)\bar{e}_{r}$
| $+1482240(\eta-3)^{2}$
$M_{2}^{2PN}$ | $8\left(4\eta^{2}-11\eta+6\right)\bar{e}_{r}^{5}$
| $+\left(161\eta^{2}-394\eta+256\right)\bar{e}_{r}^{4}$
| $+\left(49\eta^{2}-166\eta+148\right)\bar{e}_{r}^{3}$
| $+2\left(-25\eta^{2}+59\eta-48\right)\bar{e}_{r}^{2}$
| $-2\left(-33\eta^{2}+83\eta-52\right)\bar{e}_{r}$
| $+4\left(9\eta^{2}-36\eta+35\right)$
|
$+\sqrt{1-\bar{e}_{r}^{2}}\left[\left(-12\eta^{2}+56\eta-64\right)\bar{e}_{r}^{4}\right.$
| $+\left(24\eta^{2}-112\eta+128\right)\bar{e}_{r}^{3}$
| $+2\left(39\eta^{2}-95\eta+48\right)\bar{e}_{r}^{2}$
| $-2\left(27\eta^{2}-51\eta+10\right)\bar{e}_{r}$
| $\left.-4\left(9\eta^{2}-36\eta+35\right)\right]$
$M_{3}^{2PN}$ | $120\left(133\eta^{2}-349\eta+168\right)\bar{e}_{r}^{8}$
| $-240\left(33\eta^{2}-155\eta+96\right)\bar{e}_{r}^{7}$
| $+15\left(-2304\eta^{2}+5856\eta-\allowbreak 3680\right)\bar{e}_{r}^{6}$
| $-2\left(17\,256\eta^{2}+37\,520\eta+27\,520\right)\bar{e}_{r}^{5}$
| $+\left(-14\,856\eta^{2}+73\,720\eta-\allowbreak
47\,360\right)\bar{e}_{r}^{4}$
| $-8\left(-8312\eta^{2}+24050\eta-18\,080\right)\bar{e}_{r}^{3}$
| $+4\left(-6554\eta^{2}+12260\eta-\allowbreak 8720\right)\bar{e}_{r}^{2}$
| $-16\left(3709\eta^{2}-11310\eta+8660\right)\bar{e}_{r}$
| $+8\left(3049\eta^{2}-8490\eta+5660\right)$
|
$+\sqrt{1-\bar{e}_{r}^{2}}\left[15\left(587\eta^{2}-2696\eta+1264\right)\right.$
| $-2\left(8805\eta^{2}-40\,440\eta+21\,360\right)\bar{e}_{r}^{5}$
| $+\left(6177\eta^{2}-17\,840+27\,280\right)\bar{e}_{r}^{4}$
| $-8\left(5583\eta^{2}-12110\eta+12\,280\right)\bar{e}_{r}^{3}$
| $+4\left(3085\eta^{2}-1610\eta+360\right)\bar{e}_{r}^{2}$
| $+16\left(3709\eta^{2}-11310\eta+8660\right)\bar{e}_{r}$
| $\left.-8\left(3049\eta^{2}-8490\eta+5660\right)\right]$
$M_{4}^{2PN}$ | $1200\eta(851\eta-3617)\bar{e}_{r}^{6}$
| $+15\left(1731801\eta^{2}-3339132\eta+1884800\right)\bar{e}_{r}^{5}$
| $+800\left(37875\eta^{2}-119180\eta+33904\right)\bar{e}_{r}^{4}$
| $-96\left(460104\eta^{2}-1355255\eta+1016100\right)\bar{e}_{r}^{3}$
| $-25600\left(1177\eta^{2}+5202\eta-2962\right)\bar{e}_{r}^{2}$
| $-9600\left(399\eta^{2}-1814\eta+1851\right)\bar{e}_{r}$
| $-153600\left(43\eta^{2}-174\eta+135\right)$
$M_{5}^{2PN}$ | $15\eta(98967\eta+118819)\bar{e}_{r}^{6}$
| $-1920\left(8399\eta^{2}-17134\eta+12016\right)\bar{e}_{r}^{5}$
| $-4\left(5837763\eta^{2}-15119020\eta+8076560\right)\bar{e}_{r}^{4}$
| $+5120\left(2431\eta^{2}-13254\eta+662\right)\bar{e}_{r}^{3}$
| $-160\left(125911\eta^{2}-476789\eta+502740\right)\bar{e}_{r}^{2}$
| $-40960\left(68\eta^{2}+1067\eta-408\right)\bar{e}_{r}$
| $-1482240(\eta-3)^{2}$

##  Appendix B Regular evolution despite a jump in $\Delta\zeta$ when one of
the spins crosses the orbital angular momentum

We note that when either of the spins becomes perpendicular to the orbit Eq.
(46) blows up due to the $\sin^{-1}\kappa_{i}$ and $\cot\kappa_{i}$ factors,
the angle $\Delta\zeta$ becoming ill defined. We show in this subsection that
this is but a coordinate singularity.

To illustrate this, we assume that $\sin\kappa_{2}\ll 1$, i.e.

$\kappa_{2}\left(\mathfrak{t}\right)=\kappa_{\left(0\right)2}+\delta\kappa_{2}\left(\mathfrak{t}\right)\leavevmode\nobreak\
,$ (143)

with $\kappa_{\left(0\right)2}\in\left\\{0,\pi\right\\}$ and
$\left|\delta\kappa_{2}\right|\ll 1$. As said before, we are interested in the
evolution of $\mathbf{S}_{\mathbf{2}}$ across the orbital angular momentum
direction. During this $\kappa_{1}$ behaves as a quasiconstant (since
$\dot{\kappa}_{1}\propto\sin\kappa_{2}$):

$\kappa_{1}\left(\mathfrak{t}\right)=\kappa_{\left(0\right)1}+\delta\kappa_{1}\left(\mathfrak{t}\right)\leavevmode\nobreak\
,$ (144)

with $\kappa_{\left(0\right)1}$ a constant value and $\delta\kappa_{1}\ll 1$.
In addition we assume that $\mathbf{S}_{\mathbf{1}}$ points away from
$\mathbf{L}_{\mathbf{N}}$; thus,
$\sin\kappa_{\left(0\right)1}>>\cos\kappa_{\left(0\right)1}\delta\kappa_{1}$.
Then, the evolution equations for $\Delta\zeta$, $\delta\kappa_{1}$, and
$\delta\kappa_{2}$ to leading order become

$\overline{\frac{d\Delta\zeta}{d\mathfrak{t}}}=A_{\Delta\zeta}-B\frac{\cos\Delta\zeta}{\delta\kappa_{2}}\leavevmode\nobreak\
,$ (145) $\overline{\frac{d\delta\kappa_{1}}{d\mathfrak{t}}}=\epsilon
A_{\delta\kappa_{1}}\delta\kappa_{2}\sin\Delta\zeta\leavevmode\nobreak\ ,$
(146)
$\overline{\frac{d\delta\kappa_{2}}{d\mathfrak{t}}}=-B\sin\Delta\zeta\leavevmode\nobreak\
,$ (147)

with $\epsilon=\cos\kappa_{\left(0\right)2}=\pm 1$ and the coefficients

$\displaystyle\frac{A_{\Delta\zeta}}{R}$ $\displaystyle=$
$\displaystyle\nu-\frac{1}{\nu}-\epsilon\left(2\nu+1\right)x_{2}$ (148)
$\displaystyle+\left(\epsilon+\nu
x_{2}\right)w_{2}x_{2}+\left(\frac{2}{\nu}+1-w_{1}\right.$
$\displaystyle\left.-\frac{w_{1}x_{1}}{\nu}\cos\kappa_{\left(0\right)1}\right)x_{1}\cos\kappa_{\left(0\right)1}\leavevmode\nobreak\
,$ $\displaystyle\frac{B}{R}$ $\displaystyle=$
$\displaystyle\left(1+\frac{1}{\nu}-\epsilon x_{2}\right.$ (149)
$\displaystyle\left.-\frac{w_{1}x_{1}}{\nu}\cos\kappa_{\left(0\right)1}\right)x_{1}\sin\kappa_{\left(0\right)1}\leavevmode\nobreak\
,$ $\displaystyle\frac{A_{\delta\kappa_{1}}}{R}$ $\displaystyle=$
$\displaystyle\left(1+\nu-\epsilon\nu w_{2}x_{2}\right.$ (150)
$\displaystyle\left.-x_{1}\cos\kappa_{\left(0\right)1}\right)x_{2}\leavevmode\nobreak\
.$

Note that for notational simplicity we omitted the overbar form the secular
time derivatives.

From (145) and (147), we find

$\frac{d^{2}}{d\mathfrak{t}^{2}}\left[\delta\kappa_{2}\sin\Delta\zeta\right]+A_{\Delta\zeta}^{2}\delta\kappa_{2}\sin\Delta\zeta=0\leavevmode\nobreak\
,$ (151)

which gives

$\delta\kappa_{2}\sin\Delta\zeta=Q_{1}\cos\left(A_{\Delta\zeta}\mathfrak{t}+G\right)\leavevmode\nobreak\
,$ (152)

with constants $Q_{1}>0$ and $G$. Then, Eq. (147) results in

$\left(\delta\kappa_{2}\right)^{2}=-\frac{2BQ_{1}}{A_{\Delta\zeta}}\sin\left(A_{\Delta\zeta}\mathfrak{t}+G\right)+Q_{2}^{2}\leavevmode\nobreak\
,$ (153)

with an integration constant $Q_{2}^{2}$. By substituting the solutions (152)
and (153) into Eqs. (145) and (147), we find the following relation:

$Q_{2}^{2}=Q_{1}^{2}+\frac{B^{2}}{A_{\Delta\zeta}^{2}}\leavevmode\nobreak\ .$
(154)

With the notation

$C_{1}=\frac{A_{\Delta\zeta}}{B}Q_{1}\leavevmode\nobreak\ ,$ (155)

Eqs. (152) and (153) become

$\left(\delta\kappa_{2}\right)^{2}=\frac{B^{2}}{A_{\Delta\zeta}^{2}}\left[1+C_{1}^{2}-2C_{1}\sin\left(A_{\Delta\zeta}\mathfrak{t}+G\right)\right]\leavevmode\nobreak\
,$ (156)
$\sin\Delta\zeta=\frac{BC_{1}}{A_{\Delta\zeta}}\frac{\cos\left(A_{\Delta\zeta}\mathfrak{\bar{t}}+G\right)}{\delta\kappa_{2}}\leavevmode\nobreak\
,$ (157)

The minimum value of $\left(\delta\kappa_{2}\right)^{2}$ is given by
$B^{2}\left(1-C_{1}\right)^{2}/A_{\Delta\zeta}^{2}$ and when this does not
vanish, we find from Eq. (157) the following restriction for the integration
constant $C_{1}$:

$\left(\frac{C_{1}}{C_{1}-1}\right)^{2}<1\leavevmode\nobreak\ .$ (158)

In the other case when $\delta\kappa_{2}$ can vanish, the expression (157) is
regular for $\delta\kappa_{2}\rightarrow 0$ only if
$\cos\left(A_{\Delta\zeta}\mathfrak{\bar{t}}+G\right)\rightarrow 0$ at the
same time. From Eq. (156), we find that these conditions can be satisfied only
if $\sin\left(A_{\Delta\zeta}\mathfrak{t}+G\right)\rightarrow 1$ for
$\delta\kappa_{2}\rightarrow 0$ and the integration constant $C_{1}$ is $1$.
Then, the solutions read as

$\left(\delta\kappa_{2}\right)^{2}=\frac{2B^{2}}{A_{\Delta\zeta}^{2}}\left[1-\sin\left(A_{\Delta\zeta}\mathfrak{t}+G\right)\right]\leavevmode\nobreak\
,$ (159)
$\sin\Delta\zeta=\frac{B}{A_{\Delta\zeta}}\frac{\cos\left(A_{\Delta\zeta}\mathfrak{t}+G\right)}{\delta\kappa_{2}}\leavevmode\nobreak\
.$ (160)

With this expression, from (146), we have

$\delta\kappa_{1}=\epsilon\frac{A_{\delta\kappa_{1}}B}{A_{\Delta\zeta}^{2}}\sin\left(A_{\Delta\zeta}\mathfrak{t}+G\right)+D\leavevmode\nobreak\
,$ (161)

with $D$ an integration constant.

For $A_{\Delta\zeta}\mathfrak{t}+G\rightarrow M\pi/2$, where the integer $M$
is chosen such that $\delta\kappa_{2}\rightarrow 0$ at the same time, we find

$\left(\delta\kappa_{2}\right)^{2}=\frac{B^{2}}{A_{\Delta\zeta}^{2}}y^{2}\left(1-\frac{y^{2}}{12}+\mathcal{O}\left(y^{4}\right)\right)\leavevmode\nobreak\
,$ (162)
$\sin^{2}\Delta\zeta=1-\frac{y^{2}}{4}+\mathcal{O}\left(y^{4}\right)\leavevmode\nobreak\
,$ (163)

with

$y=A_{\Delta\zeta}\mathfrak{t}+G-M\frac{\pi}{2}\leavevmode\nobreak\ .$ (164)

Equations (162) and (163) show that, with the exception of the case when
$\kappa_{1}$ vanishes but otherwise for general configurations,
$\Delta\zeta\rightarrow\pm\pi/2$ as$\leavevmode\nobreak\
\kappa_{2}\rightarrow\left\\{0,\pi\right\\}$. According to the definition of
polar spin angle, $\delta\kappa_{2}$ does not change sign when the spin
crosses the axis. Thus, from Eq. (160), we find that $\sin\Delta\zeta$ must
change sign as $\delta\kappa_{2}$ vanishes, implying a jump of $\Delta\zeta$
by $\pi$ whenever $\mathbf{S}_{2}$ goes through the axis defined by
$\mathbf{L}_{\mathbf{N}}$.

Thus, we have proven that both $\cos\Delta\zeta\propto y$ and
$\delta\kappa_{2}\propto y$; thus, those terms in Eq. (46) which contain a
factor of $\cos\Delta\zeta/\sin\kappa_{2}$ remain finite as $\sin\kappa_{2}$
vanishes.

The reverse case, when $\kappa_{1}$ is close to $0$ or $\pi$ but
$\sin\kappa_{2}\ncong 0$, can be obtained by interchanging the indices
$1\leftrightarrow 2$ and $\nu\leftrightarrow\nu^{-1}$.

This results shows that the dynamics of the spin angles is well described even
if one of $\sin\kappa_{i}$ evolves through zero.

## References

  * (1) LIGO Scientific and Virgo Collaborations, GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019). [arXiv:1811.12907]
  * (2) LIGO Scientific and Virgo Collaborations, GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017). [arXiv:1710.05832]
  * (3) LIGO Scientific and Virgo Collaborations, GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, (2020) [arXiv:2010.14527]
  * (4) C. Kimball, C. P. L. Berry, and V. Kalogera, What GW170729’s exceptional mass and spin tells us about its family tree, Res. Notes AAS 4, 2 (2020). [arXiv:1903.07813]
  * (5) E. Kun, K. É. Gabányi, M. Karouzos, S. Britzen, and L. Á. Gergely, A spinning supermassive black hole binary model consistent with VLBI observations of the S5 1928+738 jet, Mon. Not. Royal Astron. Soc. 445 (2), 1370 (2014). [arXiv:1402.2644]
  * (6) B. M. Barker and R. F. O’Connell, Derivation of the Equations of Motion of a Gyroscope from the Quantum Theory of Gravitation, Phys. Rev. D 2, 1428 (1970).
  * (7) B. M. Barker and R. F. O’Connell, Gravitational two-body problem with arbitrary masses, spins, and quadrupole moments, Phys. Rev. D 12, 329 (1975)
  * (8) B. M. Barker and R. F. O’Connell, The gravitational interaction: Spin, rotation, and quantum effects-a review, Gen. Relativ. Gravit. 11, 149 (1979).
  * (9) L. E. Kidder, C. M. Will, and A. G. Wiseman, Spin effects in the inspiral of coalescing compact binaries, Phys. Rev. D 47, R4183 (1993). [gr-qc/9211025]
  * (10) L. E. Kidder, Coalescing binary systems of compact objects to $($post$)^{5/2}$-Newtonian order. V. Spin effects, Phys. Rev. D 52, 821 (1995). [arXiv:gr-qc/9506022]
  * (11) C. M. Will, Post-Newtonian gravitational radiation and equations of motion via direct integration of the relaxed Einstein equations. III. Radiation reaction for binary systems with spinning bodies, Phys. Rev. D 71, 084027 (2005). [arXiv:gr-qc/0502039]
  * (12) H. Wang and C. M. Will, Post-Newtonian gravitational radiation and equations of motion via direct integration of the relaxed Einstein equations. IV. Radiation reaction for binary systems with spin-spin coupling, Phys. Rev. D 75, 064017 (2007). [arXiv:gr-qc/0701047]
  * (13) P. C. Peters, Gravitational Radiation and the Motion of Two Point Masses, Phys. Rev. 136, B1224 (1964).
  * (14) R. Rieth and G. Schäfer, Spin and tail effects in the gravitational-wave emission of compact binaries, Class. Quantum Grav. 14, 2357 (1997).
  * (15) L. Á. Gergely, Z. I. Perjés, and M. Vasúth, Spin effects in gravitational radiation backreaction III. Compact binaries with two spinning components, Phys. Rev. D 58, 124001 (1998). [arXiv:gr-qc/9808063]
  * (16) L. Á. Gergely, Spin-spin effects in radiating compact binaries, Phys. Rev. D 61, 024035 (1999). [arXiv:gr-qc/9911082]
  * (17) L. Á. Gergely, Second post-Newtonian radiative evolution of the relative orientations of angular momenta in spinning compact binaries, Phys. Rev. D 62, 024007 (2000). [arXiv:gr-qc/0003037]
  * (18) A. Klein and P. Jetzer, Spin effects in the phasing of gravitational waves from binaries on eccentric orbits, Phys. Rev. D 81, 124001 (2010). [arXiv:1005.2046]
  * (19) E. Poisson, Gravitational waves from inspiraling compact binaries: The quadrupole-moment term, Phys. Rev. D 57, 5287 (1998). [arXiv:gr-qc/9709032]
  * (20) L. Á. Gergely and Z. Keresztes, Gravitational radiation reaction in compact binary systems: Contribution of the quadrupole-monopole interaction, Phys. Rev. D 67, 024020 (2003). [arXiv:gr-qc/0211027]
  * (21) É. Racine, Analysis of spin precession in binary black hole systems including quadrupole-monopole interaction, Phys. Rev. D 78, 044021 (2008). [arXiv:0803.1820]
  * (22) T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral, Phys. Rev. D 81, 123016 (2010). [arXiv:0911.3535]
  * (23) J. Vines, É. É. Flanagan, and T. Hinderer, Post-1-Newtonian tidal effects in the gravitational waveform from binary inspirals, Phys. Rev. D 83, 084051 (2011). [arXiv:1101.1673]
  * (24) A. Buonanno, G. Faye, and T. Hinderer, Spin effects on gravitational waves from inspiraling compact binaries at second post-Newtonian order, Phys. Rev. D 87, 044009 (2013). [arXiv:1209.6349]
  * (25) B. Mikóczi, M. Vasúth, and L. Á. Gergely, Self-interaction spin effects in inspiralling compact binaries, Phys. Rev. D 71, 124043 (2005). [arXiv:astro-ph/0504538]
  * (26) T. Marchand, L. Bernard, L. Blanchet, and G. Faye, Ambiguity-Free Completion of the Equations of Motion of Compact Binary Systems at the Fourth Post-Newtonian Order, Phys. Rev. D 97, 044023 (2018). [arXiv:1707.09289]
  * (27) L. Bernard, L. Blanchet, G. Faye, and T. Marchand, Center-of-Mass Equations of Motion and Conserved Integrals of Compact Binary Systems at the Fourth Post-Newtonian Order, Phys. Rev. D 97, 044037 (2018). [arXiv:1711.00283]
  * (28) M. Levi and J. Steinhoff, Complete conservative dynamics for inspiralling compact binaries with spins at fourth post-Newtonian order, (2016). [arXiv:1607.04252]
  * (29) G. Schäfer and P. Jaranowski, Hamiltonian formulation of general relativity and post-Newtonian dynamics of compact binaries, Liv. Rev. Rel. 21, 7 (2018). [arXiv:1805.07240]
  * (30) T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin-induced orbital precession and its modulation of the gravitational waveforms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * (31) L. Á. Gergely and P. L. Biermann, The spin-flip phenomenon in supermassive black hole binary mergers, Astrophys. J. 697, 1621 (2009). [arXiv:0704.1968]
  * (32) Gopal-Krishna, P.L. Biermann, L. Á. Gergely, and P. J. Wiita, On the origin of X-shaped radio galaxies, Res. Astron. Astrophys. 12, 127 (2012). [arXiv:1008.0789]
  * (33) M. Kesden, D. Gerosa, R. O’Shaughnessy, E. Berti, and U. Sperhake, Effective potentials and morphological transitions for binary black-hole spin precession, Phys. Rev. Lett. 114, 081103 (2015). [arXiv:1411.0674]
  * (34) D. Gerosa, M. Kesden, U. Sperhake, E. Berti, and R. O’Shaughnessy, Multi-timescale analysis of phase transitions in precessing black-hole binaries, Phys. Rev. D 92, 064016 (2015). [arXiv:1506.03492]
  * (35) C. O. Lousto and J. Healy, Flip-flopping binary black holes, Phys. Rev. Lett. 114, 141101 (2015). [arXiv:1410.3830]
  * (36) C. O. Lousto and J. Healy, Unstable flip-flopping spinning binary black holes, Phys.Rev. D 93, 124074 (2016). [arXiv:1601.05086]
  * (37) C. O. Lousto, J. Healy, and H. Nakano, Spin flips in generic black hole binaries, Phys. Rev. D 93, 044031 (2016). [arXiv:1506.04768]
  * (38) J. Majár, Spin-spin interaction in the spin-precession equations, Phys. Rev. D 80, 104028 (2009).
  * (39) D. Gerosa, A. Lima, E. Berti, U. Sperhake, M. Kesden, and R. O’Shaughnessy, Wide nutation: binary black-hole spins repeatedly oscillating from full alignment to full anti-alignment, Class. Quantum Grav. 36, 105003 (2019). [arXiv:1811.05979]
  * (40) M. Kesden, D. Gerosa, R. O’Shaughnessy, E. Berti, and U. Sperhake, Effective Potentials and Morphological Transitions for Binary Black Hole Spin Precession, Phys. Rev. Lett. 114, 081103 (2015). [arXiv:1411.0674]
  * (41) LIGO Scientific and Virgo Collaborations, GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$, Astrophys. J. Lett. 892, L3 (2020). [arXiv:2001.01761]
  * (42) D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, GW170817: Joint Constraint on the Neutron Star Equation of State from Multimessenger Observations, Astrophys. J. Lett. 852, L29 (2018). [arXiv:1711.03647]
  * (43) W. G. Laarakkers and E. Poisson, Quadrupole Moments of Rotating Neutron Stars, Astrophys. J. 512, 282 (1999). [arXiv:gr-qc/9709033]
  * (44) M. Urbanec, J. C. Miller, and Z. Stuchlik, Quadrupole moments of rotating neutron stars and strange stars, Mon. Not. Roy. Astron. Soc. 433, 1903 (2013). [arXiv:1301.5925]
  * (45) F. D. Ryan, Spinning boson stars with large self-interaction, Phys. Rev. D 55, 6081 (1997).
  * (46) N. Uchikata, S. Yoshida, and P. Pani, Tidal deformability and I-Love-Q relations for gravastars with polytropic thin shells, Phys. Rev. D 94, 064015 (2016). [arXiv:1607.03593]
  * (47) R. M. O’Leary, B. Kocsis, and A. Loeb, Gravitational waves from scattering of stellar-mass black holes in galactic nuclei, Mon. Not. R. Astron. Soc. 395, 2127 (2009). [arXiv:0807.2638]
  * (48) B. Kocsis and J. Levin, Repeated bursts from relativistic scattering of compact objects in galactic nuclei, Phys. Rev. D 85, 123005 (2012). [arXiv:1109.4170]
  * (49) J. Cuadra, P. J. Armitage, R. D. Alexander, and M. C. Begelman, Massive black hole binary mergers within subparsec scale gas discs, Mon. Not. R. Astron. Soc. 393, 1423 (2009). [arXiv:0809.0311]
  * (50) A. Sesana, Self consistent model for the evolution of eccentric massive black hole binaries in stellar environments: implications for gravitational wave observations, Astrophys. J. 719, 851 (2010). [arXiv:1006.0730]
  * (51) L. Wen, On the Eccentricity Distribution of Coalescing Black Hole Binaries Driven by the Kozai Mechanism in Globular Clusters, Astrophys. J. 598, 419 (2003). [arXiv:astro-ph/0211492]
  * (52) L. Hoffman and L. Loeb, Dynamics of triple black hole systems in hierarchically merging massive galaxies, Mon. Not. R. Astron. Soc. 377, 957 (2007). [arXiv:astro-ph/0612517]
  * (53) N. Seto, Relativistic resonant relations between massive black hole binary and extreme mass ratio inspiral, Phys. Rev. D 85, 064037 (2012). [arXiv:1202.4761]
  * (54) S. Naoz, B. Kocsis, A. Loeb, and N. Yunes, Resonant Post-Newtonian Eccentricity Excitation in Hierarchical Three-body Systems, Astrophys. J. 773, 187 (2013). [arXiv:1206.4316]
  * (55) M. Preto, I. Berentzen, P. Berczik, D. Merritt, and R. Spurzem, Merger of massive black holes using N-body simulations with post-Newtonian corrections, J. Phys. Conf. Ser. 154, 012049 (2009). [arXiv:0811.3501]
  * (56) U. Löckmann and H. Baumgardt, Tracing intermediate-mass black holes in the Galactic Centre, Mon. Not. R. Astron. Soc. 384, 323 (2008). [arXiv:0711.1326]
  * (57) L. Á. Gergely, Spinning compact binary inspiral: Independent variables and dynamically preserved spin configurations, Phys. Rev. D 81, 084025 (2010). [arXiv:0912.0459]
  * (58) L. Á. Gergely, Spinning compact binary inspiral. II. Conservative angular dynamics, Phys. Rev. D 82, 104031 (2010). [arXiv:1005.5330]
  * (59) L. Á. Gergely and Z. Keresztes, Spinning compact binary dynamics and chameleon orbits, Phys. Rev. D 91, 024012 (2015). [arXiv:1411.4057]
  * (60) T. Damour and N. Deruelle, General relativistic celestial mechanics of binary systems. I. The post-Newtonian motion, Ann. Inst. Henri Poincaré (Phys. Théor.) 43, 107 (1985).
  * (61) L. Blanchet, Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries, Living Rev. Relativity 17, 2 (2014). [arXiv:1310.1528]
  * (62) S. Ossokine, M. Boyle, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Comparing Post-Newtonian and Numerical-Relativity Precession Dynamics, Phys. Rev. D 92, 104028 (2015). [arXiv:1502.01747]
  * (63) Z. Keresztes, B. Mikóczi, and L. Á. Gergely, Kepler equation for inspiralling compact binaries, Phys. Rev. D 72, 104022 (2005). [arXiv:astro-ph/0510602]
  * (64) B. Mikóczi, Spin supplementary conditions for spinning compact binaries, Phys. Rev. D 95, 064023 (2017). [arXiv:1609.01536]
  * (65) T. D. Newton and E. P. Wigner, Localized States for Elementary Systems, Rev. Mod. Phys. 21, 400 (1949).
  * (66) M. H. L. Pryce, The Mass center in the restricted theory of relativity and its connection with the quantum theory of elementary particles, Proc. R. Soc. A 195, 62 (1948).
  * (67) J. Vinti, D. Ger, and A. Bonavito, Orbital and Celestial Mechanics, American Institute of Aeronautics and Astronautics, vol. 177, Reston (1998).
  * (68) N. Warburton, T. Osburn, and C. R. Evans, Evolution of small-mass-ratio binaries with a spinning secondary, Phys. Rev. D 96, 084057 (2017). [arXiv:1708.03720]
  * (69) Z. Keresztes and L. Á. Gergely, Stability analysis of the spin evolution fix points in inspiraling compact binaries with black hole, neutron star, gravastar or boson star components, Phys. Rev. D 103, to appear (2021).