# Disk induced binary precession:
Implications for dynamics and multi-messenger observations of black hole
binaries

Christopher Tiede Niels Bohr International Academy, Niels Bohr Institute,
Blegdamsvej 17, 2100 Copenhagen, Denmark  Daniel J. D’Orazio Niels Bohr
International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen,
Denmark  Lorenz Zwick Niels Bohr International Academy, Niels Bohr Institute,
Blegdamsvej 17, 2100 Copenhagen, Denmark  Paul C. Duffell Department of
Physics and Astronomy, Purdue University, 525 Northwestern Avenue, West
Lafayette, IN 47907, USA

(Received November 30, 2023)

###### Abstract

Many studies have recently documented the orbital response of eccentric
binaries accreting from thin circumbinary disks, characterizing the change in
binary semi-major axis and eccentricity. We extend these calculations to
include the precession of the binary’s longitude of periapse induced by the
circumbinary disk, and we characterize this precession continuously with
binary eccentricity $e_{b}$ for equal mass components. This disk-induced
apsidal precession is prograde with a weak dependence on binary eccentricity
when $e_{b}\lesssim 0.4$ and decreases approximately linearly for
$e_{b}\gtrsim 0.4$; yet at all $e_{b}$ binary precession is faster than the
rates of change to the semi-major axis and eccentricity by an order of
magnitude. We estimate that such precession effects are likely most important
for sub-parsec separated binaries with masses $\lesssim 10^{7}{M_{\odot}}$,
like LISA precursors. We find that accreting, equal-mass LISA binaries with
$M<10^{6}{M_{\odot}}$ (and the most massive $M\sim 10^{7}{M_{\odot}}$ binaries
out to $z\sim 3$) may acquire a detectable phase offset due to the disk-
induced precession. Moreover, disk-induced precession can compete with General
Relativistic precession in vacuum, making it important for observer-dependent
electromagnetic searches for accreting massive binaries—like Doppler boost and
binary self-lensing models—after potentially only a few orbital periods.

Accretion (14) — Orbits (1184) — Hydrodynamical simulations (767) —
Supermassive black holes (1663) — Gravitational waves (678)

††journal: ApJ

##  1 Introduction

Circumbinary accretion is astrophysically important for a variety of binaries
ranging from protoplanetary systems, to binary stars, to massive black holes
binaries. In recent years, much numerical work has been performed to determine
both the binary orbital and disk morphologic response to in-plane, prograde
accretion (MacFadyen & Milosavljević, 2008; Cuadra et al., 2009; Shi et al.,
2012; Shi & Krolik, 2015; D’Orazio et al., 2016; Muñoz & Lai, 2016; Miranda et
al., 2017; Tang et al., 2017; Moody et al., 2019; Tiede et al., 2020; Muñoz &
Lithwick, 2020; Duffell et al., 2020; Heath & Nixon, 2020; Dittmann & Ryan,
2022; Franchini et al., 2021). Many of these works have focused primarily on
equal-mass binaries with circular orbits, but the most recent studies have
begun characterizing accreting binary systems across binary orbital
eccentricities. In particular, both Zrake et al. (2021) and D’Orazio & Duffell
(2021) found that while near-circular binaries with eccentricity $e\lesssim
0.1$ have their eccentricity damped towards orbital circularity–where the
accretion flow causes them to expand their orbit–all other initial binary
orbital eccentricities $e\gtrsim 0.1$ are driven towards an equilibrium
eccentricity $e_{\rm eq}\sim 0.4-0.45$; and at $e_{\rm eq}$ the disk causes
the binary to shrink its semi-major axis. This is in agreement with three
values of the eccentricity studied in Muñoz et al. (2019). Siwek et al. (2023)
established that this general phenomenon holds true for all binary mass ratios
$q>0.1$ with an equilibrium eccentricity that can vary from $0.25\lesssim
e_{\rm eq}\lesssim 0.5$ (but, however, did not find a circularizing regime at
small $e$). Tiede & D’Orazio (2023) additionally investigated the orbital
response of eccentric binaries in retrograde disks and found, contrary to
prograde solutions, that the binary orbital eccentricity grows and the semi-
major axis shrinks at all eccentricities $e\leq 0.8$.

For all of these studies, measurements of the binary orbital response to
accretion from a circumbinary disk is focused on the disk mediated rate of
change to the binary semi-major axis, orbital eccentricity, and mass ratio
(Duffell et al., 2020; Dittmann & Ryan, 2023). However, in order to fully
characterize disk driven alterations to the binary orbit, one must also
consider binary apsidal precession induced by the gas forces. Such effects
have recently been noted in simulations of stellar-mass binaries embedded
inside AGN disks (Dittmann et al., 2023a; Calcino et al., 2023), but have not
been addressed in detail from full steady-state solutions of isolated binaries
accreting from thin disks. However, the presence of disk-induced binary
precession has the potential to alter existing solutions for the orbital
evolution of accreting binaries as well as to leave detectable effects in the
observations of these systems with electromagnetic and gravitational waves
(GWs).

Using existing data from D’Orazio & Duffell (2021), in this paper we calculate
the disk induced apsidal precession of the binary continuously for all
eccentricities $0<e<0.9$. In Section 2 we briefly describe the simulations
that underpin our analysis and lay out our framework for calculating and
contextualizing the induced binary precession. In Section 3 we present our
computations for the disk induced binary precession and an analysis on its
dynamical origin. Section 4 details how and when such precession may be
important for current and future modelling of accreting binary systems, and
Section 5 addresses how and when this precession might appear in observations
of massive black hole binaries with periodic light curve searches and the
space-based GW detector LISA.

##  2 Methods

The data used in this study is the same as that from D’Orazio & Duffell (2021,
hereafter DD21). We highlight the most germane aspects of the system setup and
numerical solution below and refer the reader to DD21 and Duffell et al.
(2020) for more specific details.

Data was generated using the grid-based, moving-mesh hydrodynamics code DISCO
to solve the 2D equations for viscous, locally isothermal hydrodynamics in the
presence of a time varying binary potential. The binary has equal-mass
components and is always fixed on a Keplerian orbit. The circumbinary disk is
treated in the thin-disk limit with aspect ratio $h/r\sim\mathcal{M}^{-1}=0.1$
(where $\mathcal{M}$ is the Mach number of the flow). The disk viscosity is
chosen as a constant kinematic viscosity $\nu=10^{-3}a^{2}\Omega$ where $a$ is
the binary semi-major axis and $\Omega$ is the binary orbital frequency. The
system is initially fixed on a circular orbit where it is run for 500 orbits
to reach a quasi-steady configuration; and then the eccentricity of the binary
$e_{b}$ is increased adiabatically up to $e_{b}=0.9$ over $2\times 10^{4}$
binary orbits.

As a diagnostic, DISCO outputs all forces $\mathbf{f}$ on the binary due to
the gas. The precession of the binary’s longitude of periapse can be
calculated from these forces (Murray, 1994) as

$\displaystyle\dot{\varpi}_{b}=\sqrt{\frac{a}{GM}\frac{(1-e_{b}^{2})}{e_{b}^{2}}}\left[-f_{r}\cos{\nu}+f_{\phi}\sin{\nu}\frac{2+e_{b}\cos{\nu}}{1+e_{b}\cos{\nu}}\right]\
.$ (1)

where $M$ is the total binary mass and $\nu$ is the orbital true anomaly.
There is no contribution to this precession from the accretion of mass
($\dot{M}$) and we demonstrate this in Appendix A. There is a small
contribution from the direct accretion of momentum, but these effects are
subdominant to the gravitational forces from the circumbinary material; so all
effects herein are calculated only from the gravitational forces on the
binary.

In order to translate this precession into meaningful units, we re-write
Equation 1 as

$\displaystyle\begin{split}\dot{\varpi}_{\rm
b}&=\frac{G\Sigma_{0}}{v_{b}}\sqrt{\frac{1-e_{b}^{2}}{e_{b}^{2}}}\left(\tilde{f}_{\phi}\sin{\nu}\frac{2+e_{b}\cos{\nu}}{1+e_{b}\cos{\nu}}-\tilde{f}_{r}\cos{\nu}\right)\\\
&=\frac{G\Sigma_{0}}{v_{b}}\dot{\varpi}_{\rm sim}\end{split}$ (2)

where $\Sigma_{0}$ is the density scale and $v_{b}=\sqrt{GM/a}$ is the average
binary orbital velocity. $\tilde{f}_{i}$ denotes the specific forces
(accelerations) measured in code-units where $GM=\Sigma_{0}=a=1$.
$\dot{\varpi}_{\rm sim}$ is the measured precession rate in code units.
Therefore, the physical binary precession rate from the circumbinary disk can
be written in terms of the binary orbital frequency
$\Omega_{b}=\sqrt{GM/a^{3}}$ as

$\displaystyle\dot{\varpi}_{\rm b}=\frac{\Sigma_{0}a^{2}}{M}\dot{\varpi}_{\rm
sim}\Omega_{b}=q_{D}\dot{\varpi}_{\rm sim}\Omega_{b}$ (3)

where $q_{D}=\Sigma_{0}a^{2}/M$ is the disk-to-binary mass ratio. The
magnitude of the binary precession relative to other relevant timescales
(discussed below) depends on the magnitude of $q_{D}$. Because the underlying
simulations are scale free (i.e. are true for all mass and length scales $a$,
$M$, and $\Sigma_{0}$)111To the extent that the assumptions of a thin,
radiatively efficient, and gravitationally stable disk hold true; this is
discussed further in Section 4., derivatives of the binary orbital elements
are typically reported per unit mass accreted (in units of $\dot{M}$) as
opposed to per unit time. Thus, we can also express the precession rate per
accreted mass by noting that in steady-state $\Sigma_{0}=\dot{M}_{0}/3\pi\nu$
with $\dot{M}_{0}$ the accretion rate through the disk and
$\nu=\nu_{0}a^{2}\Omega_{b}$ such that

$\displaystyle\dot{\varpi}_{\rm b}=\frac{\dot{\varpi}_{\rm
sim}}{3\pi\nu_{0}}\frac{\dot{M}_{0}}{M}\quad\text{and}\quad
q_{D}=\frac{1}{3\pi\nu_{0}}\frac{\dot{M}_{0}}{M}\,\Omega_{b}^{-1}\ .$ (4)

##  3 Results

The primary finding of DD21 was the continuous characterization of how an
accreting binary changes its semi-major axis $\dot{a}/a$ and eccentricity
$\dot{e}_{b}$ as a function of binary eccentricity. In Figure 1 we plot these
same curves for $\dot{a}/a$ and $\dot{e}_{b}$ and add the new measurement of
the binary apsidal precession rate $\dot{\varpi}_{\rm b}/2\pi$ according to
Equation 4. Apsidal precession is shown by the _black_ curve, the change in
semi-major axis by the _purple_ , and the change to eccentricity in _orange_.
We observe that disk-induced binary precession is _prograde_ and approximately
constant for $e_{b}\lesssim 0.3$, peaks at $e_{b}\approx 0.4$, and decays
approximately linearly for $e_{b}\gtrsim 0.4$. In Appendix B we also show the
binary apsidal precession induced by a _retrograde_ CBD (Tiede & D’Orazio,
2023) noting that such scenarios only change the precession rate by at most a
factor of two222One might also consider inclined disks where the interplay of
binary and disk eccentricity with binary and disk inclination can follow more
complex, Kozai-Lidov-type oscillations that are tied with the disk-induced
precession (c.f. Martin et al., 2023).. We include approximate fitting
functions for these curves in Appendix C.

Of primary note, we see that the apsidal precession of the binary is a full
order of magnitude faster than corresponding changes to the binary
eccentricity and semi-major axis. We can cultivate an intuitive understanding
for this by more closely examining the forces on the binary. Specifically, the
disk always exerts a comparatively large outward radial force on the binary
(except for short times at pericenter when $e_{b}$ is large). This is evident
in Figure 2 which shows the average radial (_orange_ curve) and azimuthal
(_green_ curve) force at each binary phase (here given by the orbital true
anomaly $\nu$) for four representative eccentricities
$e_{b}=\\{0.1,0.3,0.5,0.7\\}$. Generally speaking, the radial force is largest
near binary apocenter and smallest near pericenter because these are when the
components are nearest and furthest, respectively, from the CBD cavity edge
and the bulk of the disk. To illustrate this, a generalized gravitational
force $f=K/(r_{c}-r_{b})^{2}$ between a point mass placed at the disk cavity
radius $r_{c}=3a$ and another particle on a Keplerian orbit with phase-
dependent radius $r_{b}=a(1-e^{2})/(1+e\cos{\nu})$; with the constant
arbitrarily chosen as $K=2$. However, at larger eccentricities, this
approximation breaks down because of increasing amounts of persistent intra-
orbit material.

Figure 1: Change in binary orbital elements per unit accreted mass: longitude
of periapse (_black_), semi-major axis (_purple_), and eccentricity
(_orange_). We see that binary precession occurs on an order of magnitude
faster timescale than changes to the binary semi-major axis or eccentricity.
We note that these values can be mapped to units that depend on the disk mass
$q_{D}$ and binary orbital frequency through Equation 4.  Figure 2:  Average
forces in the radial (_orange_) and azimuthal (_green_) directions at each
binary phase (true anomaly, $\nu$). At nearly all phases, there is a large
outward radial force from the bulk of the CBD which is the dominant
contribution to the total force (except for at pericenter when $e_{b}$ is
large; e.g. $e_{b}=0.7$). The _dashed-blue_ line shows the linear force in a
toy model between a point mass at the disk cavity radius $r_{c}=3\,a$ and a
particle on a Keplerian orbit with eccentricity $e_{b}$ and radius
$r_{b}=(1-e_{b}^{2})/(1+e_{b}\cos{\nu})$ as it sweeps through true anomalies
$\nu\in[0,2\pi]$.  Figure 3:  The cumulative sum of the radial force times the
sine (_top_) and cosine (_bottom_) of the binary phase. The evolution of the
binary semi-major axis and eccentricity are proportional to the former, and we
see that the sum over a full binary orbit is comparatively very small because
of the predominantly symmetric nature of the radial force. The precession of
the binary’s longitude of pericenter, however, is proportional to the latter,
which remains comparatively large when summed over a full orbit.  Figure 4:
The binary precession rate due to gravitational forces from the circumbinary
disk $\dot{\varpi}_{\rm b}$ compared with the quadrupolar precession of the
disk itself $\dot{\varpi}_{\rm Q}$ (_left_ panel; _dotted, colored_ curves) as
well as the empirically measured precession of the disk $\dot{\varpi}_{\rm D}$
(_right_ panel; _dashed, colored_ curves) at different disk-to-binary mass
ratios $q_{D}$. In both scenarios, $\dot{\varpi}_{\rm b}$ competes with the
precession of the disk when $q_{D}\sim\mathcal{O}(10^{-2})$; and disk-induced
precession becomes the dominant effect at all eccentricities when
$q_{D}\gtrsim 1/50$.

The large radial force does not contribute significantly to $\dot{a}$ and
$\dot{e}_{b}$ because these are only affected by torques ($T=rf_{\phi}$) and
components of the work $W\sim f_{r}v_{r}\propto f_{r}\sin{\nu}$ (see Appendix
A). The latter is almost always small because $v_{r}$ is anti-symmetric around
pericenter/apocenter ($\sin\nu$ dependence), whereas the large outward portion
of $f_{r}$ due to the bulk of the disk is symmetric around the same points
(with varying $r(\nu)$). Thus, only a small anti-symmetric part of $f_{r}$
contributes to $\dot{a}$ and $\dot{e}$. This can be seen in the top panel of
Figure 3 which shows the cumulative sum

$\displaystyle\int_{-\pi/2}^{\nu}\,\tilde{f}_{r}\,\sin{\nu^{\prime}}\,d\nu^{\prime}$

(in shorthand, $\int^{\nu}\tilde{f}_{r}\sin{\nu^{\prime}}$). One can see that
the orbit average of $f_{r}\sin\nu$ is nearly zero, such that the large radial
force from the bulk of the CBD plays a subdominant role in $\dot{a}$ and
$\dot{e}_{b}$. However, the precession rate is proportional to the binary
acceleration as opposed to the orbital velocity, $\dot{\varpi}_{\rm b}\propto
f_{r}\cos{\nu}$ (see Equation 3 and Appendix A). The bottom panel of Figure 3
illustrates the corresponding orbital average
$\int^{\nu}\tilde{f}_{r}\cos\nu^{\prime}$ which is notably larger than the
purely anti-symmetric component giving a substantial contribution to the
binary precession. Therefore, as the binary moves from the point of maximal
outward radial velocity, through apocenter, to the point of maximal negative
inward velocity, it is decelerating, and the outward radial force fights this;
hence delaying the apocenter turning point. This slows the radial oscillations
of the orbit compared to the azimuthal oscillations advancing the pericenter
angle. The opposite happens during the accelerating portion of the binary
orbit. The force near apocenter is stronger than that at pericenter because
the binary is closer to the circumbinary disk, and this mismatch leads to
pericenter advance prograde with the orbit. Interestingly the $f_{r}\cos{\nu}$
symmetric part of $f_{r}$ does no work on the binary, and rather—in analogy to
the magnetic portion of the electromagnetic field tensor—only changes its
orientation via $\dot{\varpi}_{\rm b}$.

The trend towards constant time averaged $\dot{\varpi}_{\rm b}$ at small $e$
can be understood by expanding the dominant $f_{r}$ term in Equation 1 in
terms of the mean anomaly $\tilde{M}(t)$,

$-\frac{f_{r}}{e}\cos{\nu}=-\frac{f_{r}\cos{\tilde{M}}}{e}-f_{r}\left(\cos{2\tilde{M}}-1\right)-\mathcal{O}(e).$

Averaging the above over an orbit will yield the time averaged value of the
precession rate $\left<\dot{\varpi}_{\rm b}\right>_{P}$. As seen in the top
panel of Figure 2, at small $e_{b}$ we can approximate
$f_{r}=A_{r}(1-e\cos{\tilde{M}})$ with $A_{r}$ constant. The orbital average,
then, is

$\left<\dot{\varpi}_{\rm
b}\right>_{P}\approx\sqrt{\frac{a}{GM}}\left(\frac{3}{2}A_{r}\right)\approx
0.75q_{D}\Omega_{b},$ (5)

which for $A_{r}\sim 0.5$, from the top panel of Figure 2, corresponds to the
value plotted in Figure 1 at $e\lesssim 0.2$ (by Equation 4,
$13/2\pi\times\dot{M}_{0}/M\approx 0.77\,q_{D}\Omega_{b}$; see Figure 4).

##  4 Implications for CBD modelling

###  4.1 Binary vs. disk precession

Figure 5: The critical value of the binary-to-disk mass ratio $1/q_{D}$ where
the binary precession rate $\dot{\varpi}_{\rm b}$ is equal to the quadrupolar
precession rate (_green, dashed)_ at the disk cavity edge $r_{c}=3.25\,a$ and
the measured disk precession (_gold, solid)_. Below each curve
$\dot{\varpi}_{\rm b}$ dominates. Between $0.18\lesssim e_{b}\lesssim 0.4$,
the critical binary-to-disk mass ratio $1/q_{D}$ goes to infinity because the
disk no longer precesses. In this regime, the bounding limit is a viscous time
at the disk’s cavity edge which yields a critical value $1/q_{D}\approx
4.8\times 10^{3}$.

Most prior solutions for eccentric binaries accreting from steady-state disks
assume that the binary is on a fixed Keplerian orbit and does not
precess.333Although solutions from Smoothed Particle Hydrodynamics studies
that directly integrate the binary orbit ought to observe this effect for
large enough disk masses. A primary finding for circular (and low-$e_{b}$)
systems is that the CBD becomes eccentric and precesses around the binary with
a frequency comparable to that induced by the binary potential’s quadrupole
moment (e.g. MacFadyen & Milosavljević, 2008; Shi et al., 2012; Muñoz &
Lithwick, 2020)

$\dot{\varpi}_{\rm
Q}\approx\frac{3}{4}\frac{q}{(1+q)^{2}}\left(1+\frac{3}{2}e_{b}^{2}\right)\left(\frac{r_{c}}{a}\right)^{-7/2}\Omega_{b}\
$ (6)

where $q$ is the binary mass-ratio and $r_{c}$ is the approximate radius of
the disk cavity. DD21 found that this holds for eccentric binaries
$e_{b}\lesssim 0.2$, but that for $0.2\lesssim e_{b}\lesssim 0.4$ the disk
undergoes a state-transition to a predominantly symmetric, non-precessing
configuration; and for $e_{b}\gtrsim 0.5$ the disk becomes highly eccentric
and precesses around the binary slower than the associated frequency from the
binary quadrupole moment. However, it is conceivable that such behavior could
be altered if the binary is also precessing at a comparable rate.

In order to determine the relative importance of binary precession compared to
the disk precession itself, we compare the binary precession rate to that of a
test particle in the quadrupole potential of the binary and to the measured
disk precession rate $\dot{\varpi}_{\rm D}$ from DD21. This is illustrated in
Figure 4 where the black curve is equivalent to that in Figure 1, except given
in units according to Equation 3. The colored curves in the left panel show
the precession rate from the binary quadrupole moment with $q=1$ and
$r_{c}=3.25a$ scaled into equivalent units through $q_{D}$. The curves in the
right panel show the same comparison to the measured disk precession rate,
again scaled to equivalent units through $q_{D}$. The displayed values of
$q_{D}$ were chosen to highlight values of the disk-to-binary mass ratio where
the disk-induced binary precession might compete with the precession of the
disk itself. In particular we see that for disk masses $M_{D}\lesssim M/400$
the induced binary precession is always subdominant to the precession of the
disk itself (and quadrupolar precession at $r_{c}=3.25a$). However, for disk
masses $M_{D}\gtrsim M/400$ the induced binary precession can compete with the
precession of the disk at most eccentricities, and for masses $M_{D}\gtrsim
M/100$ the binary precession is always faster than that of the disk. In these
regimes, the effects of binary precession may be significant for the full
hydrodynamics solution including the disk precession itself. The exceptions to
this are for $0.2\lesssim e_{b}\lesssim 0.4$ when the disk ceases to precess,
causing the binary precession to always dominate in this regime.

For intermediate values of the disk mass $1/400\lesssim q_{D}\lesssim 1/50$ we
can solve for the critical value where the binary precession rate equals that
of the quadrupolar rate and disk rate respectively; these values are shown as
a function of eccentricity in Figure 6. The _green-dashed_ curve shows the
critical values of $1/q_{D}$ for quadrupolar precession at $r_{c}=3.25a$, and
the _gold-solid_ curve illustrates the critical values for the empirical disk
precession. The shaded regions below each curve show disk masses where binary
precession dominates. We see, again, that because the disk ceases precession
for $0.2\lesssim e_{b}\lesssim 0.4$ the critical value of $1/q_{D}$ diverges
and binary precession dominates always. In this case, the bounding timescale
would be the viscous time in the disk; i.e. if binary precession is fast and
the orientation of the gravitational potential rotates faster than the disk
can relax in response to their changing relative orientation.

###  4.2 The viscous limit

Similar to Section 4.1, we can compute an estimate of the disk mass where
binary precession (and the other evolution effects) may influence the full
solution by equating their characteristic timescales to the viscous time in
the disk. That is to say, we compute the disk mass where the binary changes
its orbit elements faster than the disk can relax and communicate these
changes viscously. We compute the characteristic number of orbits
$\tau_{\chi}$ required for an order-unity change to each of the binary orbital
elements $\chi=\\{\dot{a}/a,\dot{e}_{b},\dot{\varpi}_{\rm b}/2\pi\\}$ through
Equation 3 (and as the inverse of the curves shown in Figure 1). We note that
we can generalize Equation 3 to $\dot{\chi}=q_{d}\dot{\chi}_{\rm
sim}\Omega_{b}$, such that the timescale associated with each $\dot{\chi}$ is
set by $q_{D}$. Thus, in Figure 6 we compute the value of $q_{D}$ where the
timescale for changing each orbital element equals the viscous time at the
cavity wall $\tau_{\chi}=t_{\nu}^{\rm cav}$ as a function of the binary
eccentricity. We see that the limit for the disk’s ability to relax viscously
to changes in $\dot{a}/a$ and $\dot{e}_{b}$ occurs at disk masses
$M_{D}\approx 10^{-2}M$; but that the binary precession begins out-running the
viscous time for order-of-magnitude less massive disks, $M_{D}\approx
10^{-3}M$. Therefore, for disk masses $M_{D}\gtrsim 10^{-3}M$ it may be
important to include the disk-induced binary precession (i.e. by integrating
the binary orbit in response to forces from the CBD) in order to fully
characterize the solution. For completeness, we also overlay the limit where
the binary precession equals the empirical disk precession from Section 4.1 as
the _blue-dotted_ curve.

Figure 6: Disk-to-binary mass ratios where $\tau_{\chi}$ equals the viscous
time at the cavity edge $t_{\nu}^{\rm cav}$ as a function of binary
eccentricity. The _dashed-gold_ line shows the critical value of $q_{D}$ where
disk-induced binary precession rate equals the precession rate of the
eccentric disk $\dot{\varpi}_{\rm b}=\dot{\varpi}_{\rm D}$.  Figure 7:
Contours of $q_{D}=\Sigma_{k}a^{2}/M$ from steady-state solutions for
$\alpha$\- and $\beta$-disks (_left_ and _right_ respectively) (Haiman et al.,
2009a, eqs. 4-14) accreting at their Eddington limit with
$\Sigma_{k}=\Sigma_{(\alpha,\,\beta)}(a)$. The _dashed-cyan_ lines shows
$q_{D}=10^{-3}$ near the "viscous limit" determined in Section 4.2. Above this
line, precession effects may alter the full solution. The grey shaded region
in the upper portion of the plot illustrates the limit where disks are no-
longer stable against their own gravity. The white region in the lower right
is where the semi-major axis is within the ISCO of the binary components. We
see that disk-induced binary precession from a gravitationally stable disk
could be important for the dynamics of binaries with $M\lesssim
10^{7}{M_{\odot}}$, which are likely progenitors of mergers in the LISA band.

###  4.3 Disk mass from steady-state models

Lastly, steady-state solutions for thin accretion disks around single black
holes allow one to directly calculate the mass of the surrounding disk. Thus,
using Equations 4-14 from Haiman et al. (2009a) we calculate the disk-to-
binary mass ratio $q_{D}=\Sigma_{k}a^{2}/M$ for two steady-state solutions:
_$\alpha$ -disk_ where the disk viscosity is proportional to the total
pressure (gas + radiation), and _$\beta$ -disk_ where the viscosity is only
proportional to the gas pressure. Such solutions consider a disk with three
distinct regions: an inner-region that is radiation pressure dominated with
opacity given by electron-scattering, a middle-region that is gas-pressure
dominated with electron-scattering opacity, and an outer-region that is gas-
pressure dominated with free-free opacity. The corresponding disk density is
then taken as $\Sigma_{k}=\Sigma_{(\alpha,\,\beta)}(r=a)$. Figure 7 shows
values of the disk-to-binary mass ratio $q_{D}$ for binary’s accreting at the
Eddington rate with mass $M/{M_{\odot}}\in[10^{3}-10^{10}]$ and separations
$a\in[10^{-5},1]$ parsec (where we’ve chosen the disk viscosity according to
the $\alpha$-prescription with $\alpha=0.3$ Shakura & Sunyaev 1973). The
_left_ panel shows solutions for $\alpha$-disks and the _right_ for
$\beta$-disks. The _hatched, grey_ region illustrates where such a steady-
state solution is no longer gravitationally stable and the _white_ region in
the bottom right corner shows where the ISCO of the binary components is
larger than the semi-major axis. The _dashed, cyan_ line shows the “viscous
limit” $q_{D}=10^{-3}$ where binary precession may outpace the ability of the
circumbinary disk to relax to the changing binary apsides. Therefore, for
binaries of mass up-to $10^{7}M_{\odot}$ (likely LISA progenitors) that
experience an active, gas accretion phase during their evolution through the
sub-parsec regime, the full description of the binary evolution and associated
observational signatures may require a solution that includes such precession
effects.

##  5 Observational implications

Figure 8: Contours of binary mass and separation where $\dot{\varpi}_{\rm
b}/\dot{\varpi}_{\rm GR}=1$ for 5 values of the binary eccentricity with a
disk mass $q_{D}=10^{-3}$. The lightly-shaded regions above the curves
indicate where disk-mediated binary precession would dominate over GR-
precession for each eccentricity. The black dashed line shows binary
separations where orbital decay from the emission of GWs causes the binary to
decouple from the CBD (e.g. is faster than the viscous inflow of the disk at
its inner-edge, $r_{c}=3.25a$). The _blue, dash-dotted_ lines show binary
periods $P_{b}$ of 10 years, 1 year, and 1 month from top-to-bottom.

Binary precession can leave detectable imprints in both electromagnetic and
gravitational wave emission from accreting or inspiraling binaries, especially
in so far as it competes with precession from general relativity. We explore
these effects by comparing disk-induced precession rates to those from general
relativity in vacuum. The lines in Figure 8 show the combinations of binary
mass $M$ and semi-major axis $a$ where the rate of disk-induced precession is
equal to the orbital precession from GR in vacuum $\dot{\varpi}_{\rm
b}=\dot{\varpi}_{\rm GR}$, with

$\dot{\varpi}_{\rm GR}=\frac{3GM}{a(1-e_{b}^{2})c^{2}}\Omega_{b}\ ,$ (7)

for varying binary eccentricity and a disk mass $q_{D}=10^{-3}$. The shaded
regions show where binary precession dominates, and the white regions where GR
precession is fastest. The ratio of binary-to-GR precession can be written as

$\frac{\dot{\varpi}_{\rm b}}{\dot{\varpi}_{\rm
GR}}=\frac{a(1-e_{b}^{2})c^{2}q_{D}}{3GM}\dot{\varpi}_{\rm sim},\ $ (8)

such that the colored curves in Figure 8 simply depict the linear relationship
between mass and semi-major axis at fixed $e_{b}$ when this ratio
$\dot{\varpi}_{\rm b}/\dot{\varpi}_{\rm GR}=1$; and the eccentricity
dependence is set by the relationship in Figure 1. The dashed black line shows
the binary-disk decoupling radius where the viscous inflow rate at the disk’s
inner edge $v^{(\nu)}_{r}$ equals the orbital decay rate due to gravitational
wave radiation $\dot{a}_{\rm GW}$ (c.f. Armitage & Natarajan, 2002). For
$q_{D}=10^{-3}$, at the limit where the disk has sufficient time to viscously
relax to the changing binary orbit, GR precession becomes the dominant source
of binary precession markedly before the binary decouples from the
circumbinary disk. For smaller values of $q_{D}$ this transition recedes to
larger binary semi-major axes, while for larger $q_{D}$—although in this
regime our solutions may no longer be applicable (see Section 4.2 )—disk-
induced precession encroaches on the decoupling radius; and it dominates over
GR precession all the way down to the decoupling radius when $q_{D}\sim 0.1$.

The primary observable effect of apsidal precession, whether or not it
dominates over other forms of precession, is the accumulation of extra orbital
phase relative to an unperturbed orbit. We show below that this can arise in
observations of periodic lightcurves from eccentric accreting binaries or in
the dephasing of GW emission from supermassive black hole binaries. To lowest
order in orbital eccentricity, the magnitude of this change over a period is

$\left<\delta\theta\right>_{P}\approx\dot{\varpi}_{\rm b}P_{b}\sim 2\pi q_{D}$
(9)

where we note that $\dot{\varpi}_{\rm b}$ and this dephasing do not
necessarily vanish for very small or zero eccentricities because external
forces can still alter the effective gravitational potential and the orbital
velocity of the binary (see Appendix D).

###  5.1 Periodic Light Curve Searches for Black Hole Binaries

The relative contributions of disk-induced precession vs. general relativistic
precession can be significant for electromagnetic searches for massive
binaries in galactic nuclei. A subset of search methods aim to identify
periodic features in target system lightcurves that can be connected to the
underlying orbital period of an accreting binary (e.g., Graham et al., 2015;
Charisi et al., 2016; Liu et al., 2019; Chen et al., 2020, 2022). Orbital
periods targeted in these searches range from weeks to years (e.g., Haiman et
al., 2009b; Xin & Haiman, 2021; Haiman et al., 2023). For comparison, Figure 8
shows _blue, dashed-dotted_ lines for three binaries periods, from top-to-
bottom, $P_{b}\in$ [10 years, 1 year, 1 month]; like those that might be
identifiable in electromagnetic searches.

Binary models for the origin of periodic variability employ either
hydrodynamic variability caused by variations in the binary accretion rate
(e.g. D’Orazio et al., 2013; Farris et al., 2014; Tang et al., 2018; Dittmann
& Ryan, 2022; Westernacher-Schneider et al., 2022; Gutiérrez et al., 2022) or
observer-dependent relativistic Doppler Boosts (D’Orazio et al., 2015;
D’Orazio & Haiman, 2017; Charisi et al., 2018) and binary self-lensing events
(D’Orazio & Di Stefano, 2018; Hu et al., 2020; Kelley et al., 2021; Davelaar &
Haiman, 2022; Major Krauth et al., 2023b)—both modulated by the binary orbit.
Hydrodynamic accretion-rate variability would likely not be sensitive to disk-
induced binary precession for $q_{D}\lesssim 10^{-3}$ when $\dot{\varpi}_{\rm
b}$ is slow compared to the other relevant timescales, but it is possible that
the characteristic frequencies of accretion variability change in the limit
$q_{D}\gtrsim 10^{-3}$ (i.e., Figure 6); we leave exploration of this effect
to future work.

Relativistic boosting and lensing, however, are sensitive to—and straight-
forwardly depend on—the orbital reconstruction. Specifically, the shape of
periodic modulations caused by the orbital Doppler boost along with the shape,
magnification, and crucially, timing of lensing flares depend on the orbital
eccentricity and argument of pericenter with respect to the observer’s line of
sight (see Fig. 7 in D’Orazio & Charisi, 2023); and thus, also on changes to
these quantities over time, such as the advance of pericenter. These models
sometimes include general relativistic precession in their modelling, but
Figure 8 demonstrates that a circumbinary disk sourcing the binary accretion
with $q_{D}\sim 10^{-3}$ may induce a comparable effect for massive binaries
$10^{4}\lesssim M/M_{\odot}\lesssim 10^{7}$ with orbital periods of order
months-to-years.

To estimate when precession from a circumbinary disk would significantly
affect the shape and timing of Doppler+lensing signatures, we calculate the
number of binary orbits needed for the accumulated precession angle to equal
the width of a binary self-lensing flare444This assumes at least one lensing
flare has been temporally resolved and so deviation from the non-precessing
model can be verified by measuring a later flare that is offset by greater
than one flare width.. This limit guarantees at least an observable change in
flare timing, if not also a change in lightcurve shape due to an altered
azimuthal viewing angle caused by apsidal precession. From Equation 9 we
equate the accumulated precession angle $\delta\theta\approx 2\pi q_{D}N_{\rm
orb}$ with the duration of the lensing event divided by the orbital period
(e.g., Equation (3) of D’Orazio & Di Stefano (2018)). Then the required number
of orbits needed to accumulate significant precession is,

$\displaystyle N_{\rm orb}$ $\displaystyle\gtrsim$
$\displaystyle\frac{(1+q_{b})^{-1/2}}{\pi^{5/3}q_{D}}\left(\frac{2GM}{c^{3}P_{\rm
b}}\right)^{1/3}$ $\displaystyle\approx$ $\displaystyle
1.2\left(\frac{q_{D}}{10^{-3}}\right)^{-1}\left(\frac{M}{10^{7}{M_{\odot}}}\right)^{1/3}\left(\frac{P_{\rm
b}}{2\mathrm{yr}}\right)^{-1/3}.$

In the second line we evaluate for $q_{b}=1$ and a binary at the edge of the
disk-precession-dominated regime in Figure 8. Hence, before relativistic
apsidal precession becomes dominant, the disk-induced precession can cause
significant changes to the Doppler-lensing periodic modulations over only a
few orbital periods. This motivates including disk orbital precession when
modelling periodic variability due to binary self-lensing. Moreover, such
lightcurve models with the inclusion of disk precession allow measurement of
the local disk mass and can provide otherwise lacking constraints.

Forces from the disk which cause pericenter precession also alter the orbital
velocity of the binary away from the Keplerian value. In our case, the
dominant contribution is from outward radial disk forces which do not strongly
affect orbital evolution (see Section 3), but do act to slow down the binary
compared to the Keplerian value for a given separation and binary mass. For
Doppler+lensing signatures this effect would bias the binary parameters that
one would recover assuming an orbit in vacuum by, an albeit, very small amount
of order $q_{D}$.

Figure 9: The critical disk mass $q_{D,\,\rm crit}$ required to induce a
detectable dephasing in a LISA source with source frame total mass $M$ at
redshift $z$. The _blue-dashed_ line indicates the viscous limit $q_{D}\sim
10^{-3}$ determined in Section 4.2.

###  5.2 Dephasing of Gravitational Waves in LISA

Disk-induced precession in the late stages of massive binary evolution could
also potentially have important consequences for near-equal-mass binary
inspiral events that will be detected by LISA. Binaries in the LISA band
($10^{3}\lesssim M/M_{\odot}\lesssim 10^{7}$) can stay coupled to their
circumbinary disks anywhere from $\sim 1$ year down to minutes prior to merger
(Dittmann et al., 2023b; Major Krauth et al., 2023a); or they may not decouple
at all (Avara et al., 2023). Therefore, even though the disk-induced
precession in-band is likely small compared to that from GR, the responsible
disk-forces could still feasibly alter gravitational wave emission and could
leave a detectable imprint on LISA signals.

In accordance with LISA expectations, we focus on sources with very low
eccentricity555This is consistent with estimates for the in-band eccentricity
$e_{b}\sim 10^{-3}-10^{-4}$ for equal-mass binaries that have had their
eccentricity driven near the equilibrium eccentricity by a prograde
circumbinary disk (Zrake et al., 2021; Garg et al., 2023). Apsidal precession
is then best understood as a change in orbital frequency of magnitude
$\sim\dot{\varpi}_{\rm b}$. However, for a monochromatic source with an
observed GW frequency $f_{\rm GW}$, there would be no simple way to
distinguish between a vacuum source emitting at the “Keplerian” frequency
$f_{\mathrm{GW}}/2=f_{k}$, or at a slightly perturbed frequency

$\displaystyle\frac{f_{\mathrm{GW}}}{2}=f_{p}=f_{k}\left(1+\frac{\Delta
f}{f_{k}}\right),$ (11)

where $\Delta f=(2\pi P_{k})^{-1}\left<\delta\theta\right>_{P}\approx
q_{D}f_{k}$, assuming the binary eccentricity is small and $\dot{\varpi}_{\rm
sim}\sim 1$. The two sources would generally only be distinguishable by their
frequency evolution (or chirping) because they radiate gravitational energy at
slightly different rates. Here we provide a simple estimate of this effect.

The radiated energy flux is $\dot{E}_{\mathrm{GW}}\sim
f_{\mathrm{GW}}^{6}a^{4}\sim\dot{f}$ such that a binary observed at
$f_{\mathrm{GW}}$ that is perturbed by surrounding gas will exist at a
slightly modified separation $a_{p}(f_{\mathrm{GW}})\sim(2\pi f_{p})^{-2/3}$
(compared to that on a Keplerian orbit in vacuum, $a_{k}$). Thus, the modified
chirp evolution can be expressed

$\displaystyle\frac{\dot{f}_{p}}{\dot{f}_{k}}\sim\frac{a_{p}(f_{\rm
GW})^{4}}{a_{k}(f_{\rm GW})^{4}}\sim\left(1+q_{D}\right)^{-8/3}\ .$ (12)

The total accumulated phase of a GW signal is $\phi(f_{\mathrm{GW}})=2\pi\int
f/\dot{f}_{p}df$ such that the difference in accumulated phase between the
perturbed source and one in vacuum (for small $\Delta f$) is

$\displaystyle\delta\phi=\phi_{\mathrm{GW}}-\phi_{k}\approx\frac{16\pi}{3}q_{D}\phi_{k}$
(13)

where $\phi_{k}\sim f^{-5/3}$ is the total phase accumulated by a source in
vacuum. Such a dephasing is considered detectable in the LISA band (without
accounting for degeneracies) if $\delta\phi\gtrsim 10\times\textrm{SNR}^{-1}$
with SNR the signal-to-noise ratio of the event (Kocsis et al., 2011). Thus,
Equation 13 gives a detectability condition for the disk mass as a function of
the event SNR and the total number of orbits in-band $N\sim\phi_{k}/2\pi$,

$\displaystyle
q_{D}\gtrsim\frac{3}{32\pi^{2}}\frac{10}{\textrm{SNR}}N^{-1}\approx\frac{0.1}{N\times
SNR}\ .$ (14)

Equal mass mergers in the LISA band will have typical SNRs of order
$\mathcal{O}(10^{1}-10^{3})$ and will spend $N\sim 10-10^{3}$ orbits in band
(with the exception of the most massive $10^{7}{M_{\odot}}$ binaries which
will have $\textrm{SNR}\sim 10-100$ and $N\sim 1-10$).

Figure 9 illustrates the critical disk mass $q_{D,{\rm crit}}$ necessary to
satisfy Equation 14 as a function of the total source frame binary mass $M$
and redshift $z$ by calculating $\textrm{SNR}(M,z)$ and $\phi_{k}(M,z)$ from
PhenomA model waveforms (Ajith et al., 2007) and an approximate LISA
sensitivity curve (Robson et al., 2019). The _blue-dashed_ line indicates the
$q_{D}=10^{-3}$ viscous limit determined from Section 4.2 such that the
majority of LISA sources are detectable for $q_{D}\gtrsim 10^{-5}$; and the
most nearby sources for $q_{D}\gtrsim 10^{-6}$. However, this estimate serves
only as a lower bound on the disk mass required to induce a detectable
dephasing in LISA because we have ignored any degeneracies in the induced
chirp behavior, have assumed each source spends the maximum amount of time in-
band possible for a 4 year mission lifetime, and for the lowest values of
$q_{D}$ it may start to compete with the frequency resolution of the detector.
Furthermore, a more rigorous treatment ought to self-consistently connect the
circumbinary disk forces to the time evolution of the binary quadrupole moment
and the resultant chirp behavior. We reserve this calculation for future work,
but emphasize that because disk-induced precession dominates over the other
orbital perturbations, it is the most likely candidate for detectable disk-
induced signals in LISA detections of near equal-mass binaries.

##  6 Conclusions

We have calculated from viscous hydrodynamics simulations the apsidal
precession of eccentric binaries accreting from circumbinary disks. Such
simulations have previously been used to compute the disk-driven rate of the
change to the binary semi-major axis and eccentricity, and we have found that
the induced precession effect is an order of magnitude faster than these
changes (Figure 1). Moreover, we identify the primary source of this
precession and the fundamental reason for its comparative significance as the
relatively strong and symmetric outward radial gravitational force from the
bulk of the circumbinary disk (Figure 2 & 3).

The degree to which $\dot{\varpi}_{\rm b}$ competes with other timescales in
the problem scales with the disk-to-binary mass ratio $q_{D}$, and for
sufficiently large disk masses it is possible this precession could alter the
full solutions on which previous orbital-evolution results and our own
findings are based. We estimated these disk masses by first comparing
$\dot{\varpi}_{\rm b}$ to the precession of a test particle in the quadrupole
moment of the binary’s potential and of the circumbinary disk itself (Figure
4). We found that the disk-induced binary precession generally dominates over
these timescales when $q_{D}\gtrsim 1/300$ with the exception of binaries with
$0.2\lesssim e_{b}\lesssim 0.4$ where the circumbinary disk is empirically
found not to precess. We compared the timescale for $\dot{\varpi}_{\rm b}$ to
the viscous time in the disk (Figure 6) and found that for disk masses
$q_{D}\gtrsim 10^{-3}$ the induced binary precession may occur faster than the
disk can viscously relax to the changing binary apsides, and that such
situations may warrant more dedicated study to determine if the precession of
the binary alters the full solution. We additionally computed physical values
of $q_{D}$ from $\alpha$ and $\beta$ steady-state disk solutions (Figure 7)
and determined that circumbinary accretion solutions and orbital evolution
modelling for binaries with masses up-to $\sim 10^{7}M_{\odot}$ and at sub-
parsec separations may need to consider the effects of disk-induced binary
precession.

For observational purposes, we compared the disk-induced binary precession to
precession rates from general relativity (Figure 8) and discussed implications
for both electromagnetic and gravitational wave searches for accreting and
coalescing massive binaries. Specifically, we determined that disk-induced
precession may be significant over only a few orbital periods for modelling
sources with Doppler-lensing periodic modulations, and that the precession may
source an extra phase accumulation in accreting LISA systems that is generally
detectable (in the best case scenario and ignoring degeneracies) for equal
mass binaries with disk-binary mass ratios $q_{D}\gtrsim 10^{-5}$. Thus, we
have concluded that the effects of precession on existing solutions for binary
orbital evolution and circumbinary accretion signatures warrants future, more
detailed investigation, and that it should be considered in both current and
future observations of accreting massive binaries in electromagnetic surveys
and gravitational wave experiments.

##  7 Acknowledgments

The authors sincerely thank Bin Liu, Alessandro Trani, Andrew MacFadyen, and
members of the NBIA theoretical astrophysics group for insightful discussions.
D.J.D. and C.T. received support from the Danish Independent Research Fund
through Sapere Aude Starting Grant No. 121587. P.D. acknowledges support from
the National Science Foundation under grant AAG-2206299.

## References

  * Ajith et al. (2007) Ajith P., et al., 2007, Classical and Quantum Gravity, 24, S689
  * Armitage & Natarajan (2002) Armitage P. J., Natarajan P., 2002, ApJ, 567, L9
  * Avara et al. (2023) Avara M. J., Krolik J. H., Campanelli M., Noble S. C., Bowen D., Ryu T., 2023, arXiv e-prints, p. arXiv:2305.18538
  * Calcino et al. (2023) Calcino J., Dempsey A. M., Dittmann A. J., Li H., 2023, arXiv e-prints, p. arXiv:2311.13727
  * Charisi et al. (2016) Charisi M., Bartos I., Haiman Z., Price-Whelan A. M., Graham M. J., Bellm E. C., Laher R. R., Marka S., 2016, preprint,  (arXiv:1604.01020)
  * Charisi et al. (2018) Charisi M., Haiman Z., Schiminovich D., D’Orazio D. J., 2018, MNRAS, 476, 4617
  * Chen et al. (2020) Chen Y.-C., et al., 2020, MNRAS, 499, 2245
  * Chen et al. (2022) Chen Y.-J., et al., 2022, arXiv e-prints, p. arXiv:2206.11497
  * Cuadra et al. (2009) Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, MNRAS, 393, 1423
  * D’Orazio & Charisi (2023) D’Orazio D. J., Charisi M., 2023, arXiv e-prints, p. arXiv:2310.16896
  * D’Orazio & Di Stefano (2018) D’Orazio D. J., Di Stefano R., 2018, MNRAS, 474, 2975
  * D’Orazio & Duffell (2021) D’Orazio D. J., Duffell P. C., 2021, ApJ, 914, L21
  * D’Orazio & Haiman (2017) D’Orazio D. J., Haiman Z., 2017, MNRAS, 470, 1198
  * D’Orazio et al. (2013) D’Orazio D. J., Haiman Z., MacFadyen A., 2013, MNRAS, 436, 2997
  * D’Orazio et al. (2015) D’Orazio D. J., Haiman Z., Schiminovich D., 2015, Nature, 525, 351
  * D’Orazio et al. (2016) D’Orazio D. J., Haiman Z., Duffell P., MacFadyen A., Farris B., 2016, Monthly Notices of the Royal Astronomical Society, 459, 2379
  * Davelaar & Haiman (2022) Davelaar J., Haiman Z., 2022, Phys. Rev. D, 105, 103010
  * Dittmann & Ryan (2022) Dittmann A. J., Ryan G., 2022, MNRAS, 513, 6158
  * Dittmann & Ryan (2023) Dittmann A. J., Ryan G., 2023, arXiv e-prints, p. arXiv:2310.07758
  * Dittmann et al. (2023a) Dittmann A. J., Dempsey A. M., Li H., 2023a, arXiv e-prints, p. arXiv:2310.03832
  * Dittmann et al. (2023b) Dittmann A. J., Ryan G., Miller M. C., 2023b, ApJ, 949, L30
  * Duffell et al. (2020) Duffell P. C., D’Orazio D., Derdzinski A., Haiman Z., MacFadyen A., Rosen A. L., Zrake J., 2020, ApJ, 901, 25
  * Farris et al. (2014) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ, 783, 134
  * Franchini et al. (2021) Franchini A., Sesana A., Dotti M., 2021, MNRAS, 507, 1458
  * Garg et al. (2023) Garg M., Tiwari S., Derdzinski A., Baker J., Marsat S., Mayer L., 2023, arXiv e-prints, p. arXiv:2307.13367
  * Graham et al. (2015) Graham M. J., et al., 2015, Nature, 518, 74
  * Gutiérrez et al. (2022) Gutiérrez E. M., Combi L., Noble S. C., Campanelli M., Krolik J. H., López Armengol F., García F., 2022, ApJ, 928, 137
  * Haiman et al. (2009a) Haiman Z., Kocsis B., Menou K., Lippai Z., Frei Z., 2009a, Classical and Quantum Gravity, 26, 094032
  * Haiman et al. (2009b) Haiman Z., Kocsis B., Menou K., 2009b, ApJ, 700, 1952
  * Haiman et al. (2023) Haiman Z., et al., 2023, arXiv e-prints, p. arXiv:2306.14990
  * Heath & Nixon (2020) Heath R. M., Nixon C. J., 2020, A&A, 641, A64
  * Hu et al. (2020) Hu B. X., D’Orazio D. J., Haiman Z., Smith K. L., Snios B., Charisi M., Di Stefano R., 2020, MNRAS, 495, 4061
  * Kelley et al. (2021) Kelley L. Z., D’Orazio D. J., Di Stefano R., 2021, MNRAS, 508, 2524
  * Kocsis et al. (2011) Kocsis B., Yunes N., Loeb A., 2011, Phys. Rev. D, 84, 024032
  * Liu et al. (2019) Liu T., et al., 2019, ApJ, 884, 36
  * MacFadyen & Milosavljević (2008) MacFadyen A. I., Milosavljević M., 2008, ApJ, 672, 83
  * Major Krauth et al. (2023a) Major Krauth L., Davelaar J., Haiman Z., Westernacher-Schneider J. R., Zrake J., MacFadyen A., 2023a, arXiv e-prints, p. arXiv:2304.02575
  * Major Krauth et al. (2023b) Major Krauth L., Davelaar J., Haiman Z., Westernacher-Schneider J. R., Zrake J., MacFadyen A., 2023b, arXiv e-prints, p. arXiv:2310.19766
  * Martin et al. (2023) Martin R. G., Lepp S., Zhang B., Nixon C. J., Childs A. C., 2023, MNRAS,
  * Miranda et al. (2017) Miranda R., Muñoz D. J., Lai D., 2017, MNRAS, 466, 1170
  * Moody et al. (2019) Moody M. S. L., Shi J.-M., Stone J. M., 2019, The Astrophysical Journal, 875, 66
  * Muñoz & Lai (2016) Muñoz D. J., Lai D., 2016, ApJ, 827, 43
  * Muñoz & Lithwick (2020) Muñoz D. J., Lithwick Y., 2020, ApJ, 905, 106
  * Muñoz et al. (2019) Muñoz D. J., Miranda R., Lai D., 2019, ApJ, 871, 84
  * Murray (1994) Murray C. D., 1994, Icarus, 112, 465
  * Murray & Dermott (2000) Murray C. D., Dermott S. F., 2000, Solar System Dynamics, doi:10.1017/CBO9781139174817.
  * Robson et al. (2019) Robson T., Cornish N. J., Liu C., 2019, Classical and Quantum Gravity, 36, 105011
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shi & Krolik (2015) Shi J.-M., Krolik J. H., 2015, ApJ, 807, 131
  * Shi et al. (2012) Shi J.-M., Krolik J. H., Lubow S. H., Hawley J. F., 2012, The Astrophysical Journal, 749, 118
  * Siwek et al. (2023) Siwek M., Weinberger R., Hernquist L., 2023, MNRAS, 522, 2707
  * Tang et al. (2017) Tang Y., MacFadyen A., Haiman Z., 2017, MNRAS, 469, 4258
  * Tang et al. (2018) Tang Y., Haiman Z., MacFadyen A., 2018, MNRAS, 476, 2249
  * Tiede & D’Orazio (2023) Tiede C., D’Orazio D. J., 2023, arXiv e-prints, p. arXiv:2307.03775
  * Tiede et al. (2020) Tiede C., Zrake J., MacFadyen A., Haiman Z., 2020, ApJ, 900, 43
  * Westernacher-Schneider et al. (2022) Westernacher-Schneider J. R., Zrake J., MacFadyen A., Haiman Z., 2022, Phys. Rev. D, 106, 103010
  * Xin & Haiman (2021) Xin C., Haiman Z., 2021, MNRAS, 506, 2408
  * Zrake et al. (2021) Zrake J., Tiede C., MacFadyen A., Haiman Z., 2021, ApJ, 909, L13

##  Appendix A $\dot{\varpi}_{\rm b}$ derivation with accretion

Starting with the elliptic equation

$\displaystyle r=\frac{a(1-e_{b}^{2})}{1+e_{b}\cos{\nu}}$ (A1)

with $\nu=\theta-\varpi_{b}$ the true anomaly, $\theta$ the position angle
measured from the line of nodes, $\varpi_{b}$ the longitude of pericenter, $a$
the semi-major axis, and $e_{b}=\sqrt{1+2\ell^{2}E\mu^{-2}}$ the orbital
eccentricity. $\ell=\sqrt{\mu a(1-e_{b}^{2})}$ is the specific angular
momentum, $E=-\mu/2a$ is the specific orbital energy, and $\mu=GM$ with $M$
the total mass. Therefore, we write

$\displaystyle\begin{split}\ell^{2}&=\mu r(1+e_{b}\cos{\nu})\\\ &=\mu
r\left[1+\sqrt{1+2\ell^{2}E\mu^{-2}}\cos{(\theta-\varpi_{b})}\right]\
.\end{split}$ (A2)

Differentiating with respect to time while keeping the instantaneous radius
fixed (e.g. because of some external force) yields

$\displaystyle\begin{split}2\ell\dot{\ell}=\dot{\mu}r(1&+e_{b}\cos{\nu})+\mu
re_{b}^{-1}\left[\frac{E}{\mu^{2}}2\ell\dot{\ell}+\frac{\ell^{2}}{\mu^{2}}\dot{E}-\frac{2\ell^{2}E}{\mu^{2}}\frac{\dot{\mu}}{\mu}\right]\cos{\nu}+\mu
re_{b}\left[-\dot{\theta}+\dot{\varpi}_{b}\right]\sin{\nu}\ .\end{split}$ (A3)

The energy and angular momentum derivatives are given as

$\displaystyle\begin{split}&\dot{E}=\frac{\mu}{2a^{2}}\dot{a}-\frac{\dot{\mu}}{2a}\\\
2\ell&\dot{\ell}=\dot{\mu}a(1-e_{b}^{2})+\mu\dot{a}(1-e_{b}^{2})-2\mu
a\dot{e}_{b}\end{split}$ (A4)

such that all terms proportional to $\dot{\mu}$ in Equation A3 cancel.
Therefore, as expected, there is no contribution to the binary precession from
the accretion of mass. Solving for $\dot{\varpi}$, then, gives

$\displaystyle\dot{\varpi}$
$\displaystyle=\dot{\theta}+\frac{2\ell\dot{\tilde{\ell}}}{\mu
re_{b}\sin{\nu}}-\frac{\cot{\nu}}{e_{b}^{2}\mu^{2}}\left[2E\ell\dot{\tilde{\ell}}+\ell^{2}\dot{\tilde{E}}\right]$
$\displaystyle=\dot{\theta}+\left(\frac{r^{-1}-E(e\mu)^{-1}\cos{\nu}}{\mu
e_{b}\sin{\nu}}\right)2\ell\dot{\tilde{\ell}}-\frac{\ell^{2}\cot{\nu}}{\mu^{2}e_{b}^{2}}\dot{\tilde{E}}\
.$ (A5)

where $\dot{\tilde{\ell}}$ and $\dot{\tilde{E}}$ are the change in energy and
specific angular momentum due only to an external specific force $\mathbf{f}$

$\displaystyle\begin{split}\dot{\tilde{E}}&=f_{r}v_{r}+f_{\phi}v_{\phi}=\sqrt{\frac{\mu}{a(1-e_{b}^{2})}}\bigg{[}f_{r}e_{b}\sin{\nu}+f_{\phi}(1+e_{b}\cos{\nu})\bigg{]}\\\
\dot{\tilde{\ell}}&=rf_{\phi}\end{split}$ (A6)

where we’ve used the fact that $v_{r}=\sqrt{\mu
a^{-1}/(1-e_{b}^{2})}e_{b}\sin\nu$ and $v_{\phi}=\sqrt{\mu
a^{-1}/(1-e_{b}^{2})}(1+e_{b}\cos\nu)$. Lastly, we take $\dot{\theta}=0$
because we are only considering in-plane forces, fixing the binary longitude
of ascending node (Murray & Dermott, 2000). Plugging in to Equation A5 yields

$\displaystyle\dot{\varpi}_{b}=e_{b}^{-1}\sqrt{a\mu^{-1}(1-e_{b}^{2})}\left[-f_{r}\cos{\nu}+f_{\phi}\sin{\nu}\frac{2+e\cos{\nu}}{1+e\cos{\nu}}\right]\
.$ (A7)

##  Appendix B Retrograde CBDs

Tiede & D’Orazio (2023) also explored the orbital response of equal-mass
binaries accreting from retrograde CBDs. For completeness, we measure the
disk-induced precession for retrograde solutions; this is shown in Figure A1
as the _blue_ curve alongside $\dot{\varpi}_{\rm b}$ from the prograde
configuration. Despite the comparative symmetry and lack of resonances in the
retrograde scenario (see Tiede & D’Orazio 2023 for a detailed discussion), a
retrograde CBD still drives _prograde_ apsidal precession in the binary. This
disk-induced precession in retrograde configurations is $\sim 2\times$ faster
than in prograde disks at low binary eccentricity $e_{b}\lesssim 0.15$; but
decreases monotonically with growing $e_{b}$ such that the effect is
comparable by $e_{b}\gtrsim 0.4$. The comparatively large precession rate in
retrograde solutions is consistent with the relatively large magnitudes
measured for $\dot{a}/a$ and $\dot{e}_{b}$, and we posit that these effects
are due to stronger gravitational forcing associated with a less truncated CBD
(especially at low $e_{b}$) that is more tightly coupled to the binary.
Nevertheless, all analysis presented in this paper would change by at most a
factor of two if considering retrograde instead of prograde accretion
scenarios.

##  Appendix C $\dot{\varpi}_{\rm b}$ fitting functions

Figure A1 also includes as _dashed_ lines approximate fitting functions for
$\dot{\varpi}_{\rm b}$ in both prograde (_black_) and retrograde (_blue_)
configurations. The piecewise fitting function for prograde solutions is taken
as

$\displaystyle\dot{\varpi}_{\rm
fit}=\frac{\dot{M}_{0}}{M}\begin{cases}B_{0}&e<e_{0}\\\
B_{0}+m_{1}(e-e_{0})&e_{0}\leq e<e_{1}\\\ B_{1}+m_{2}(e-e_{1})&e_{1}\leq
e<e_{2}\\\ B_{2}+m_{3}(e-e_{2})&e>e_{2}\end{cases}$ (A1)

where $B_{0}=0.770$, $B_{1}=0.878,B_{2}=0.680$, $m_{1}=0.600$,
$m_{2}=-5.351,m_{3}=-0.860$, and $e_{0}=0.198$, $e_{1}=0.378$, $e_{2}=0.415$.
For retrograde solutions we approximate the disk-induced precession as

$\displaystyle\dot{\varpi}_{\rm
fit\\_r}=\frac{\dot{M}_{0}}{M}\begin{cases}C_{0}+n_{0}e&e<e_{0}^{\prime}\\\
C_{1}+n_{1}e&e\geq e_{0}^{\prime}\end{cases}$ (A2)

where $C_{0}=1.810$, $C_{1}=1.503$, $n_{0}=-2.738$, $n_{1}=-1.242$, and
$e_{0}^{\prime}=0.212$.

Figure A1: Binary precession from a retrograde CBD compared to the prograde
solution. Fitting functions provided in Appendix C are shown by the associated
_dashed_ lines.

##  Appendix D Phase Accumulation from Osculating Orbits

The angle tracking the true position of the orbit in the sky is
$\theta=\nu+\varpi_{b}$, where the true anomaly $\nu$ denotes the angle with
respect to the argument of pericenter $\dot{\varpi}_{\rm b}$. The true anomaly
under perturbations can be expressed,

$\dot{\nu}=\dot{\nu}_{k}-\dot{\varpi}_{\rm b},$ (A1)

where $\dot{\nu}_{k}$ is the Keplerian value of the true anomaly advance.
However, by definition

$\dot{\theta}=\dot{\nu}+\dot{\varpi}_{\rm b},$ (A2)

so that $\dot{\theta}=\dot{\nu}_{k}$. This implies that the only variation in
$\dot{\theta}$ is caused by the the perturbed Kepler equation. Thus, under
perturbed orbital elements $\chi_{p}=\chi_{0}+\chi_{1}$, the rate at which the
binary sweeps out this true longitude can be expressed

$\dot{\theta}=\dot{\nu}_{k}(\chi_{p})=\sqrt{\frac{\mu}{p^{3}_{p}}}\left(1+e_{p}\cos{\nu_{p}}\right)^{2},$
(A3)

where $p=a(1-e^{2})$, and $\chi_{1}=\int^{t}_{t_{0}}\dot{\chi}dt$ is
proportional to the magnitude of the perturbing force and so assumed small
over an orbit. We expand the true anomaly by writing
$\dot{\nu}=\dot{\nu_{k0}}-\dot{\varpi}_{\rm b}$ and integrate—without yet
specifying the time dependence of $\dot{\varpi}_{\rm b}$—so that
$\nu_{p}(t)=\nu_{k0}(t)-\varpi_{b}(t)$. Here $\nu_{k0}$ is shorthand for
$\nu_{k}(\chi_{0})$. We expand in the limit of slow precession and small
perturbations, keeping only quantities to first order in these small
quantities (e.g. $a_{p}=a_{0}+a_{1}$ and so forth) such that Equation A3
becomes

$\displaystyle\dot{\theta}=\dot{\nu}_{k0}+\dot{\nu}_{k0}\left[\left(-\frac{3}{2}\frac{a_{1}}{a_{0}}+\frac{3e_{0}e_{1}}{1-e^{2}_{0}}\right)+\frac{2}{1+e_{0}\cos{\nu_{k0}}}\left(e_{1}\cos{\nu_{k0}}+e_{0}\varpi_{b}(t)\sin{\nu_{k0}}\right)\right].$
(A4)

Dephasing is from the second term, and $\varpi_{b}(t)$ enters only in the last
term. Expanding the phase-dependent $\cos$ and $\sin$ terms in time via the
mean anomaly $\tilde{M}=\Omega_{b}t$ to 1st order in $e_{0}$, we write

$\delta\dot{\theta}(t)=2\sqrt{\frac{\mu}{p^{3}_{0}}}\left[e_{0}\varpi_{b}(t)\sin{\tilde{M}}-\frac{3}{4}\frac{a_{1}}{a_{0}}\left(1+2e_{0}\cos{\tilde{M}}\right)+e_{0}e_{1}\left(1+\frac{3}{2}\cos{2\tilde{M}}\right)\right].$
(A5)

Now consider the contribution from each of the perturbed quantities after one
orbit, $\varpi_{b}(P_{b})$, $a_{1}(P_{b})$, and $e_{1}(P_{b})$, in the small
$e_{0}$ limit where orbital dephasing is most relevant for LISA sources. To
see the contribution from $\varpi_{b}(P_{b})$ we expand
$e_{0}\dot{\varpi}_{\rm b}$ as,

$\displaystyle e_{0}\dot{\varpi}_{\rm b}$
$\displaystyle\approx\sqrt{\frac{p_{0}}{\mu}}\left[-f_{r}\cos{\nu_{k0}}+2f_{\phi}\sin{\nu_{k0}}\left(2-e_{0}\cos{\nu_{k0}}\right)\right]+\mathcal{O}(e^{2}_{0})$
$\displaystyle\approx\sqrt{\frac{p_{0}}{\mu}}\left[-f_{r}\cos{\tilde{M}}+4f_{\phi}\sin{\tilde{M}}\right]+\mathcal{O}(e_{0})$
(A6)

where we’ve once again expanded the oscillating functions of true anomaly in
mean anomaly (time) to linear order in $e_{0}$. Then the leading order
contribution to $\varpi_{b}$ and the dephasing will depend on the functional
dependence of the forces.

For prograde, coplanar accretion the radial force dominates $f_{r}\gg
f_{\phi}$ and for small $e_{0}$ is almost constant in time. For larger
$e_{0}$, $f_{r}$ has significant $\cos{\tilde{M}}$ components (see Figure 2).
Therefore, we consider two specific scenarios: that of constant radial force
$f_{r}=A_{r}$ and that where $f_{r}=-A_{r}\cos{\tilde{M}}$. For each of these
scenarios $a_{1}/a_{0},\,e_{0}e_{1}\sim f_{\phi}+\mathcal{O}(e_{0})$ (as
discussed in §3), and the dominant term in $\delta\dot{\theta}$ is the
precession term.

For the constant force, the largest contributions to $\varpi_{b}$ will be
those that do not vanish in the integral over an orbit. To leading order,

$\displaystyle
e_{0}\varpi_{b}(t)=-\frac{a^{2}_{0}}{\mu}A_{r}\sin{\tilde{M}}+\mathcal{O}(e_{0},f_{\phi})$
(A7)

and the dephasing over an orbit is given by

$\displaystyle\delta{\theta}(P_{b})$ $\displaystyle\approx$ $\displaystyle
2\sqrt{\frac{\mu}{p^{3}_{0}}}\int^{P_{b}}_{0}e_{0}\varpi_{b}(t)\sin{\tilde{M}}dt$
(A8) $\displaystyle\approx$
$\displaystyle-2\pi\frac{A_{r}a^{2}_{0}}{\mu}\propto-2\pi q_{D},$

where we’ve used $A_{r}/(\mu/a_{0}^{2})\sim q_{D}$ because the force is
proportional to the disk mass.

For, instead, $f_{r}=-A_{r}\cos{\tilde{M}}$,

$\displaystyle
e_{0}\varpi_{b}(t)=\frac{1}{4}\frac{a^{2}_{0}}{\mu}\left[A_{r}\left(2\tilde{M}+\sin{2\tilde{M}}\right)\right]+\mathcal{O}(e_{0},f_{\phi})$
(A9)

and

$\displaystyle\delta{\theta}(P_{b})\approx-2\pi\frac{A_{r}a^{2}_{0}}{\mu}\propto-2\pi
q_{D}.$ (A10)

Hence, for both scenarios, the dephasing from apsidal precession is to leading
order given by $2\pi q_{D}\sim\dot{\varpi}_{\rm b}P_{b}$.