Disk-induced Binary Precession: Implications for Dynamics and Multimessenger Observations of Black Hole Binaries

Many studies have recently documented the orbital response of eccentric binaries accreting from thin circumbinary disks, characterizing the change in the binary semimajor 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 the binary eccentricity when e b ≲ 0.4 and decreases approximately linearly for e b ≳ 0.4; yet at all e b binary precession is faster than the rates of change to the semimajor axis and eccentricity by an order of magnitude. We estimate that such precession effects are likely most important for subparsec separated binaries with masses ≲107 M ⊙, like LISA precursors. We find that accreting, equal-mass LISA binaries with M < 106 M ⊙ (and the most massive M ∼ 107 M ⊙ binaries out to z ∼ 3) may acquire a detectable phase offset due to the disk-induced precession. Moreover, disk-induced precession can compete with general relativistic precession in a 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.


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 equalmass 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) andD'Orazio &Duffell (2021) found that while near-circular binaries with eccentricity e ≲ 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 ≳ 0.1 are driven towards an equilibrium eccentricity e eq ∼ 0.4−0.45;and at e 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 ≲ e eq ≲ 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 ≤ 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.

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, movingmesh 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 ∼ M −1 = 0.1 (where M is the Mach number of the flow).The disk viscosity is chosen as a constant kinematic viscosity ν = 10 −3 a 2 Ω where a is the binary semi-major axis and Ω 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 × 10 4 binary orbits.
As a diagnostic, DISCO outputs all forces 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 where M is the total binary mass and ν is the orbital true anomaly.There is no contribution to this precession from the accretion of mass ( Ṁ ) 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 1as where Σ 0 is the density scale and v b = GM/a is the average binary orbital velocity.fi denotes the specific forces (accelerations) measured in code-units where GM = Σ 0 = a = 1.π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 where q D = Σ 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 Σ 0 ) 1 , derivatives of the binary orbital elements are typically reported per unit mass accreted (in units of Ṁ ) as opposed to per unit time.Thus, we can also express the precession rate per accreted mass by noting that in steady-state Σ 0 = Ṁ0 /3πν with Ṁ0 the accretion rate through the disk and ν = ν 0 a 2 Ω b such that The primary finding of DD21 was the continuous characterization of how an accreting binary changes its semimajor axis ȧ/a and eccentricity ėb as a function of binary eccentricity.In Figure 1 we plot these same curves for ȧ/a and ėb and add the new measurement of the binary apsidal precession rate πb /2π 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 diskinduced binary precession is prograde and approximately constant for e b ≲ 0.3, peaks at e b ≈ 0.4, and decays approximately linearly for e b ≳ 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 two 2 .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 semimajor 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 ν) 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 ν) with the constant arbitrarily chosen as K = 2.However, at larger eccentricities, 2 One 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).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 qD and binary orbital frequency through Equation 4.
this approximation breaks down because of increasing amounts of persistent intra-orbit material.
The large radial force does not contribute significantly to ȧ and ėb because these are only affected by torques (T = rf ϕ ) and components of the work W ∼ f r v r ∝ f r sin ν (see Appendix A).The latter is almost always small because v r is anti-symmetric around pericenter/apocenter (sin ν 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(ν)).Thus, only a small anti-symmetric part of f r contributes to ȧ and ė.This can be seen in the top panel of Figure 3 which shows the cumulative sum ν fr sin ν ′ ).One can see that the orbit average of f r sin ν is nearly zero, such that the large radial force from the bulk of the CBD plays a subdominant role in ȧ and ėb .However, the precession rate is proportional to the binary acceleration as opposed to the orbital velocity, πb ∝ f r cos ν (see Equation 3and Appendix A).The bottom panel of Figure 3 illustrates the corresponding orbital average ν fr cos ν ′ 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  Average forces in the radial (orange) and azimuthal (green) directions at each binary phase (true anomaly, ν).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  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.
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 ν 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 πb .
The trend towards constant time averaged πb at small e can be understood by expanding the dominant f r term in Equation 1 in terms of the mean anomaly M (t), Binary v. quadrupole precession Averaging the above over an orbit will yield the time averaged value of the precession rate ⟨ πb ⟩ P .As seen in the top panel of Figure 2, at small e b we can approximate f r = A r (1 − e cos M ) with A r constant.The orbital average, then, is which for A r ∼ 0.5, from the top panel of Figure 2, corresponds to the value plotted in Figure 1 at e ≲ 0.2 (by Equation 4, 13/2π × Ṁ0 /M ≈ 0.77 q D Ω b ; see Figure 4).Lastly, we note that because the precession rate is dominated by the outward radial pull of the CBD, we expect our results to be comparatively insensitive to system uncertainties like the disk thermodynamics (unlike the other orbital elements; see e.

disk precession
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. 3A 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 3 Although solutions from Smoothed Particle Hydrodynamics studies that directly integrate the binary orbit ought to observe this effect for large enough disk masses.
. The critical value of the binary-to-disk mass ratio 1/qD where the binary precession rate πb is equal to the quadrupolar precession rate (green, dashed) at the disk cavity edge rc = 3.25 a and the measured disk precession (gold, solid).Below each curve πb dominates.Between 0.18 ≲ e b ≲ 0.4, the critical binary-to-disk mass ratio 1/qD 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/qD ≈ 4.8 × 10 3 .
the binary potential's quadrupole moment (e.g.Mac-Fadyen & Milosavljević 2008;Shi et al. 2012 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 ≲ 0.2, but that for 0.2 ≲ e b ≲ 0.4 the disk undergoes a state-transition to a predominantly symmetric, non-precessing configuration; and for e b ≳ 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 π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 ≲ 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 ≳ M/400 the induced binary precession can compete with the precession of the disk at most eccentricities, and for masses M D ≳ 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 ≲ e b ≲ 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 ≲ q D ≲ 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 5.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 ≲ e b ≲ 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.

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 τ χ required for an order-unity change to each of the binary orbital elements χ = { ȧ/a, ėb , πb /2π} through Equation 3 (and as the inverse of the curves shown in Figure 1).We note that we can generalize Equation 3 to χ = q d χsim Ω b , such that the timescale associated with each χ 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 τ χ = t cav ν as a function of the binary eccentricity.We see that the limit for the disk's ability to relax viscously to changes in ȧ/a and ėb occurs at disk masses M D ≈ 10 −2 M ; but that the binary precession begins out-running the viscous time for order-of-magnitude less massive disks, M D ≈ 10 −3 M .Therefore, for disk masses M D ≳ 10 −3 M it may be important to include the diskinduced 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.

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 = Σ k a 2 /M for two steady-state solutions: α-disk where the disk viscosity is proportional to the total pressure (gas + radiation), and β-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 electronscattering, a middle-region that is gas-pressure dominated with electron-scattering opacity, and an outerregion that is gas-pressure dominated with free-free opacity.The corresponding disk density is then taken as Σ k = Σ (α, β) (r = a).Figure 7 shows values of the diskto-binary mass ratio q D for binary's accreting at the Eddington rate with mass M/M ⊙ ∈ [10 3 − 10 10 ] and separations a ∈ [10 −5 , 1] parsec (where we've chosen the disk viscosity according to the α-prescription with α = 0.3 Shakura & Sunyaev 1973).The left panel shows solutions for α-disks and the right for β-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 ⊙ (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.

OBSERVATIONAL IMPLICATIONS
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 πb = πGR , with 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 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 πb / π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 (ν) r equals the orbital decay rate due to gravitational wave radiation ȧ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 semimajor 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 ∼ 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 where we note that π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).  2009a, eqs. 4-14) accreting at their Eddington limit with Σ k = Σ (α, β) (a).The dashed-cyan lines shows qD = 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 ≲ 10 7 M⊙, which are likely progenitors of mergers in the LISA band.

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. 2020Chen et al. , 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 topto-bottom, P b ∈ [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 selflensing 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 ≲ 10 −3 when πb is slow compared to the other relevant timescales, but it is possible that the charac-teristic frequencies of accretion variability change in the limit q D ≳ 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 ∼ 10 −3 may induce a comparable effect for massive binaries 104 ≲ M/M ⊙ ≲ 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 flare 4 .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 9we equate the accumulated precession angle δθ ≈ 2πq D N 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, 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 .

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 ≲ M/M ⊙ ≲ 10 7 ) can stay coupled to their circumbinary disks anywhere from ∼ 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 diskinduced 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 eccentricity5 .Apsidal precession is then best understood as a change in orbital frequency of magnitude ∼ πb .However, for a monochromatic source with an observed GW frequency f GW , there would be no simple way to distinguish between a vacuum source emitting at the "Keplerian" frequency f GW /2 = f k , or at a slightly perturbed frequency where ∆f = (2πP k ) −1 ⟨δθ⟩ P ≈ q D f k , assuming the binary eccentricity is small and π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 ĖGW ∼ f 6 GW a 4 ∼ ḟ such that a binary observed at f GW that is perturbed by surrounding gas will exist at a slightly modified separation a p (f GW ) ∼ (2πf p ) −2/3 (compared to that on a Keplerian orbit in vacuum, a k ).Thus, the modified chirp evolution can be expressed The total accumulated phase of a GW signal is ϕ(f GW ) = 2π f / ḟp df such that the difference in accumulated phase between the perturbed source and one in vacuum (for small ∆f ) is where ϕ k ∼ 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 δϕ ≳ 10 × SNR −1 with SNR the signalto-noise ratio of the event (Kocsis et al. 2011).Thus, Equation 13gives a detectability condition for the disk mass as a function of the event SNR and the total number of orbits in-band Equal mass mergers in the LISA band will have typical SNRs of order O(10 1 −10 3 ) and will spend N ∼ 10−10 3 orbits in band (with the exception of the most massive 10 7 M ⊙ binaries which will have SNR ∼ 10 − 100 and N ∼ 1 − 10).
Figure 9 illustrates the critical disk mass q D,crit necessary to satisfy Equation 14 as a function of the total source frame binary mass M and redshift z by calculating SNR(M, z) and ϕ 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 ≳ 10 −5 ; and the most nearby sources for q D ≳ 10 −6 .However, this estimate serves only as a lower bound on the disk mass required to induce a detectable dephasing in LISA 6 because we have ignored degeneracies in the induced chirp behavior (e.g.Garg et al. 2022), have assumed each source spends the maximum amount of time in-band possible for a 4 year mission lifetime, and have neglected the disk mass evolution with shrinking semi-major axis; and for the lowest values of q D this signal 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.

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 π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 π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 ≳ 1/300 with the exception of binaries with 0.2 ≲ e b ≲ 0.4 where the circumbinary disk is empirically found not to precess.We compared the timescale for πb to the viscous time in the disk (Figure 6) and found that for disk masses q D ≳ 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 α and β steady-state disk solutions (Figure 7) and determined that circumbinary accretion solutions and orbital evolution modelling for binaries with masses up-to ∼ 10 7 M ⊙ and at subparsec separations may need to consider the effects of disk-induced binary precession.
For observational purposes, we compared the diskinduced 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 ≳ 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.
Figure 2.Average forces in the radial (orange) and azimuthal (green) directions at each binary phase (true anomaly, ν).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 rc = 3 a and a particle on a Keplerian orbit with eccentricity e b and radius r b = (1 − e 2 b )/(1 + e b cos ν) as it sweeps through true anomalies ν ∈ [0, 2π].
Figure 2.Average forces in the radial (orange) and azimuthal (green) directions at each binary phase (true anomaly, ν).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 rc = 3 a and a particle on a Keplerian orbit with eccentricity e b and radius r b = (1 − e 2 b )/(1 + e b cos ν) as it sweeps through true anomalies ν ∈ [0, 2π].

Figure 3 .
Figure3.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 .
Figure 4.The binary precession rate due to gravitational forces from the circumbinary disk πb compared with the quadrupolar precession of the disk itself πQ (left panel; dotted, colored curves) as well as the empirically measured precession of the disk πD (right panel; dashed, colored curves) at different disk-to-binary mass ratios qD.In both scenarios, πb competes with the precession of the disk when qD ∼ O(10 −2 ); and disk-induced precession becomes the dominant effect at all eccentricities when qD ≳ 1/50.

Figure 6 .
Figure 6.Disk-to-binary mass ratios where τχ equals the viscous time at the cavity edge t cav ν as a function of binary eccentricity.The dashed-gold line shows the critical value of qD where disk-induced binary precession rate equals the precession rate of the eccentric disk πb = πD.

Figure 7 .
Figure7.Contours of qD = Σ k a 2 /M from steady-state solutions for αand β-disks (left and right respectively,Haiman et al.  2009a, eqs.4-14)  accreting at their Eddington limit with Σ k = Σ (α, β) (a).The dashed-cyan lines shows qD = 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 ≲ 10 7 M⊙, which are likely progenitors of mergers in the LISA band.

Figure 8 .
Figure8.Contours of binary mass and separation where πb / πGR = 1 for 5 values of the binary eccentricity with a disk mass qD = 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, rc = 3.25a).The blue, dash-dotted lines show binary periods P b of 10 years, 1 year, and 1 month from topto-bottom.

Figure 9 .
Figure9.The critical disk mass qD, 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 qD ∼ 10 −3 determined in Section 4.2.

Figure A1 .
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.