# Eccentric Binaries in Retrograde Disks

Christopher Tiede1 & Daniel J. D’Orazio1
1Niels Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen, Denmark E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted November 13, 2023. Received November 13, 2023; in original form July
10, 2023)

###### Abstract

Modern numerical hydrodynamics tools have recently enabled detailed
examinations of binaries accreting from prograde circumbinary disks. These
have re-framed the current understanding of binary-disk interactions and disk
driven orbital evolution. We present the first full-domain grid-based
hydrodynamics simulations of equal-mass, eccentric binaries accreting from
retrograde circumbinary disks. We study binary eccentricities that span
$e=0.0$ to $e=0.8$ continuously, and explore the influence of retrograde
accretion on the binary orbital response, disk morphology, and observational
properties. We find that, at all eccentricities, retrograde accretion shrinks
the binary semi-major axis and pumps its eccentricity leading to the
previously identified possibility of highly eccentric mergers. Contrary to
past studies and models, we observe gravitational forces to dominate the
binary’s orbital evolution as opposed to the physical accretion of mass and
momentum. Retrograde accretion variability also differs strongly from prograde
solutions. Preeminently, binaries with $e>0.55$ reveal a unique two-period,
double-peaked accretion signature that has not previously been identified. We
additionally find evidence for the emergence of retrograde Lindblad resonances
at large eccentricities in accordance with predictions from linear theory. Our
results suggest that some astrophysical binaries for which retrograde
accretion is possible will experience factors-of-a-few times faster orbital
decay than in prograde disks and will have their eccentricities pumped beyond
the limits found from prograde solutions. Such effects could lead to rapid
inward migration for some young stellar binaries, the detection of highly-
eccentric LISA mergers, and the tentatively observed turnover at the low-
frequency end of the gravitational wave background.

###### keywords:

hydrodynamics – software:simulations – black hole mergers – gravitational
waves – quasars:general

††pubyear: 2024††pagerange: Eccentric Binaries in Retrograde
Disks–Gravitational Softening

##  1 Introduction

Accretion onto a binary from a surrounding circumbinary disk (CBD) is
important for the evolution and observation of many types of astrophysical
binaries ranging from protoplanetary systems, to binary stars, to massive
black hole binaries (Kley & Nelson, 2012; Orosz et al., 2012; Barnes &
Hernquist, 1996; Mayer, 2013). A plethora of work has been completed in recent
years to determine the effect of prograde CBD’s on the orbital evolution of
the inner binary and on the associated observational characteristics
(MacFadyen & Milosavljević, 2008; Cuadra et al., 2009; Shi et al., 2012;
D’Orazio et al., 2013; Shi & Krolik, 2015; Muñoz & Lai, 2016; Miranda et al.,
2017; Tang et al., 2017; Moody et al., 2019; Tiede et al., 2020; Duffell et
al., 2020; Zrake et al., 2021; D’Orazio & Duffell, 2021; Dittmann & Ryan,
2022; Franchini et al., 2022; Siwek et al., 2023a). A primary result from the
majority of these studies is that for the fiducial setup of an equal mass,
circular binary embedded in an $\alpha$-disk (Shakura & Sunyaev, 1973) with
characteristic scale height $h/r=0.1$ and turbulent viscosity parameter
$\alpha=0.1$, the CBD delivers angular momentum to the binary and drives
binary outspiral (Muñoz et al., 2019). The most recent studies have largely
focused on subsequently filling out the parameter space beyond this reference
model by varying the binary mass ratio (Duffell et al., 2020; Siwek et al.,
2023a), the binary eccentricity (Zrake et al., 2021; D’Orazio & Duffell, 2021;
Siwek et al., 2023b), the disk scale-height (Tiede et al., 2020; Heath &
Nixon, 2020; Dittmann & Ryan, 2022; Penzlin et al., 2022; Franchini et al.,
2022), and the disk cooling (Sudarshan et al., 2022; Wang et al., 2022, 2023).
Nearly all of these studies have found that both the direction and magnitude
of binary inspiral depends non-trivially on each of these parameters; see Lai
& Muñoz (2022) for a recent review.

Many of these works have also detailed observational characteristics of these
disks through measured accretion rates onto the binary as a proxy for the
emitted accretion flux (e.g. Farris et al., 2015; Tang et al., 2018). A
primary result from prograde disks is that binaries accrete at a rate that is
on average equal to that of a single central object. The underlying accretion
modulations occur near the binary orbital frequency (and its harmonics) as
well as on longer timescales at $\sim 5\times$ the orbital period due to
periodic over-feeding from a lopsided cavity formed in circular and near-
circular binaries (MacFadyen & Milosavljević, 2008; Shi et al., 2012; D’Orazio
et al., 2013; Farris et al., 2014; Tang et al., 2017; Duffell et al., 2020;
Dittmann & Ryan, 2022). Eccentric binaries in prograde disks, however, lose
the longer $\sim 5$ orbit period and instead have their accretion variability
dominated by the binary orbital frequency (Hayasaki et al., 2007; Dunhill et
al., 2015; Miranda et al., 2017; Zrake et al., 2021; Westernacher-Schneider et
al., 2022).

Of primary relevance to this study, both D’Orazio & Duffell (2021) and Zrake
et al. (2021) found that equal-mass binaries that start with small
eccentricity $e\lesssim 0.1$ have their eccentricity damped towards circular
orbits—at which they are driven apart by the gas—but all other initial binary
eccentricities $e\gtrsim 0.1$ are driven towards an equilibrium eccentricity
$e\sim 0.4-0.45$. Further, Siwek et al. (2023b) demonstrated that this general
behavior holds for all binary mass ratios $q>0.1$, and the equilibrium
eccentricity varies between $0.25\lesssim e_{\rm eq}\lesssim 0.5$. Such an
equilibrium eccentricity could manifest in populations of stellar binaries
that have undergone gas accretion phases as well as in residual eccentricity
measurements of merging super-massive black hole binaries (SMBHB’s) in
gravitational waves with the space-based interferometer LISA.

Nearly all modern numerical studies have focused on prograde CBD’s; but in
some astrophysical systems it is not clear a priori what the angular momentum
of the CBD ought to be relative to the binary. In particular, SMBHB’s that
undergo active accretion phases in the late-stages of their evolution
following a galaxy merger may receive gas injections isoptropically, and it is
possible that retrograde gas feeding onto a SMBHB is comparably likely to
prograde configurations (King & Pringle, 2006; Nixon et al., 2011b; Hobbs et
al., 2011). It is also possible that a non-negligible fraction of young
stellar binaries born in dense, chaotic star forming regions acquire
retrograde CBD’s (e.g., $\sim 10\%$ in Elsender et al. 2023; Bate et al.
2010). This has motivated past numerical studies of retrograde-CBD’s around
circular binaries with viscous, smoothed particle hydrodynamics (Nixon et al.,
2011a; Nixon & Lubow, 2015) and inviscid 3D magnetohydrodynamics (Bankert et
al., 2015). Roedig & Sesana (2014) have also simulated retrograde, self-
gravitating disks around eccentric binaries for a few eccentricities up to
$e=0.8$, and Amaro-Seoane et al. (2016) examined the effect of different sink
prescriptions for binaries with mass ratio $q=1/10$. These results have also
motivated a series of approximate analytic models for retrograde binary
accretion (Nixon et al., 2011a; Roedig & Sesana, 2014; Schnittman & Krolik,
2015; Amaro-Seoane et al., 2016). We review these works and their findings in
more detail in §2, but generally speaking, they all find that a retrograde-CBD
causes accreting binaries at all eccentricities and mass ratios to shrink
their separation (inspiral) and to pump their eccentricity (above some small
threshold).

That said, the simulations performed for these studies were either focused on
a few specific parameter setups or did not simulate the full multi-scale
domain including both the entire CBD and the binary with its intra-orbit
material111This has been demonstrated to be integral to converging on a
solution in prograde scenarios Farris et al. 2014; Tang et al. 2017; Muñoz et
al. 2019; with possible exceptions in confined regions of parameter space
e.g., Tiede et al. 2022.. Moreover, the developed analytic models of
retrograde binary orbital evolution relied on these simulations and typically
assumed the dominant contribution to be inelastic collisions with parcels of
gas at binary apocenter.

In addition to being equally astrophysically viable as prograde solutions,
retrograde circumbinary systems are analytically interesting because at low
binary eccentricity they lack the Lindblad resonances that strongly affect—and
possibly dominate (see, however, Mahesh et al., 2023)—the prograde scenario
(Goldreich & Tremaine, 1979, 1980). Thus, retrograde CBD’s offer an
informative insight into the relevant dynamics of accretion onto eccentric
binaries by comparison to their prograde counterparts.

In this paper we present the first full-domain, grid-based simulations of
equal-mass binaries accreting from retrograde CBD’s and expand upon recent
works exploring the orbital evolution of eccentric binaries in prograde
configurations – this work can be directly compared with the recent prograde
results of D’Orazio & Duffell (2021); Zrake et al. (2021). §2 reviews the
previous literature on retrograde binary accretion. §3 lays out the numerical
techniques used in this study. §4 includes our numerical results for binary
evolution (§4.1), details of the disk morphology and how it changes with
eccentricity (§4.2), an analysis of the disk-mediated eccentricity evolution
and its primary drivers (§4.3, §4.4), the appearance of retrograde resonances
(§4.5), and variability signatures from retrograde accretion (§4.6). We
discuss how our results compare with previous works and explore observational
implications in §5, and summarize in §6. The Appendix includes an additional
study on how our results depend on the specifics of the sink prescription and
gravitational softening.

##  2 Retrograde Results in the Literature

The initial motivation for the study of retrograde CBD’s was as an
astrophysically plausible mechanism to shepherd SMBHB’s from the canonical
dynamical friction stalling radius (Milosavljević & Merritt, 2003) down to the
sub-parsec separations where gravitational waves could merge the binary within
a Hubble time. At the time, analytic and numerical works had suggested that
prograde circumbinary disks would absorb angular momentum from the binary,
becoming either decretion disks (Webbink, 1976; Pringle, 1991) (that transfer
mass outwards) or unstable (Lodato et al., 2009) halting binary migration.

The angular momentum transport in the prograde interaction was thought to be
dominated by Lindblad resonances where the disk angular velocity $\Omega(r)$
was equal to the combination of binary orbital frequency $\omega_{b}$ and
positive integer mode number $m=1,2,...$ as $\Omega(r)=m\,\omega_{b}/(m\pm
1)$. However, this is only valid for $|\Omega|>|\omega_{b}|$ such that these
resonances are not present when $\Omega$ and $\omega_{b}$ have opposite sign.
Nixon et al. (2011a) explored the possibility that such retrograde CBD’s could
be a promising mode for shrinking an accreting SMBHB toward a GW dominated
regime and coalescence. They employ a toy model where the binary-disk
interaction is assumed to only occur at apocenter and pericenter where the
binary semi-major axis $a$ and eccentricity $e$ evolve assuming inelastic
collisions between the smaller, secondary binary component of mass $M_{2}$ and
Keplerian gas parcels of some mass $\Delta M$. Their model finds that all
binaries shrink their semi-major axis due to the loss of energy and angular
momentum and note that eccentricity is pumped at binary apocenter and damped
at pericenter. As a result, binaries that start with small eccentricity tend
to remain circular because of the relative balance between effects at
pericenter and apocenter; but binaries that start with eccentricity greater
than the characteristic scale-height of the disk $e\gtrsim h/r$ will have
their evolution dominated by interactions at apocenter and have their
eccentricity pumped until it is near unity and the binary coalesces. They
compare their toy model with a limited set of 3D viscous smoothed particle
hydrodynamics (SPH) simulations and confirm that secondary-gas interactions at
apocenter drive eccentricity and at pericenter damp eccentricity for a binary
mass ratio of $q=1/10$. They also include one simulation with $q=0.5$ and very
small initial eccentricity and demonstrate that it stays near circular; but
they do not include an example simulation of the case with initial
eccentricity $e_{0}>h/r$ where they would expect rapid eccentricity growth
towards binary coalescence. Nixon et al. (2011a) also point out that such an
interaction ought to drive eccentricity in the interacting gas parcels and
could potentially drive an eccentric CBD.

Roedig & Sesana (2014) simulated a more extensive set of eccentric binaries in
retrograde CBD’s using 3D SPH with $\beta$-cooling and self-gravity (instead
of a typical viscosity). All of their simulations were for near-equal mass
binaries with $q=1/3$, and similar to Nixon et al. (2011a), they find that
near-circular binaries in retrograde CBD’s remain nearly circular while
shrinking their semi-major axis (at a rate comparable to a prograde, $e=0$
control). For eccentricities above a critical value $e>0.04$, they find that
both $\\{\dot{a},\,\dot{e}\\}(e)\propto e$ such that eccentricity grows
exponentially until the binary merges at eccentricities near unity.
Additionally, they amend the toy model of Nixon et al. (2011a) based on an
impact-interaction at binary apocenter to include terms $\propto 1+q$ for when
$q\ll 1$ is not satisfied and find an overall qualitative agreement; although
they note that their model slightly underestimates both the binary semi-major
axis shrinking and the eccentricity growth for $e\gtrsim 0.3$.

It is worth mentioning that one major discrepancy between the simulations of
Nixon et al. (2011a) and Roedig & Sesana (2014) was that the former observed
the presence of material around each binary component in the form of circum-
single “minidisks”, while the latter did not. Roedig & Sesana (2014)
attributed this to their cooling prescription.

The first grid based simulations of retrograde circumbinary accretion were
performed by Bankert et al. (2015) in 3D Newtonian MHD. These simulations were
for equal mass binaries $q=1$ fixed on a circular orbit where the grid had an
inner boundary at $r=0.8a$ (with $a$ the binary separation) such that
hydrodynamics in direct proximity of the binary components themselves were not
resolved. This study used time- and azimuthally-averaged radial torque density
profiles from $r>0.8a$ to conclude that very little angular momentum is
transferred gravitationally between the binary and the disk and that the
dominant effect in the orbital evolution of the binary is from the direct
accretion of mass and retrograde angular momentum; in general agreement with
the modelling of Nixon et al. (2011a); Roedig & Sesana (2014). Bankert et al.
(2015) also found the the retrograde CBD does not become eccentric and
maintains its general axisymmetry unlike its prograde counterpart, but that
the retrograde MHD disk around a circular binary does still exhibit spiral
density perturbations despite the lack of standard Linblad and co-rotation
resonances.

Schnittman & Krolik (2015) used the conclusion from Bankert et al. (2015) that
the binary evolution is dominated by the direct accretion of mass and momentum
to develop an analytic model for binary separation and eccentricity evolution
as a function of both binary eccentricity and mass-ratio. This model was again
based around the “impact approximation” that all energy and angular momentum
exchange occurs at apocenter; however, their model is built around a value of
the specific angular momentum exchange measured from Bankert et al. (2015) and
includes the possibility of differential accretion between the two binary
components. Consistent with others, their model predicted semi-major axis
shrinking for all eccentricities and mass ratios, but in contrast with
previous modelling posited that eccentricity driving is maximal at small $e$
and decreases with growing eccentricity.

Using 2D viscous hydrodynamics, Amaro-Seoane et al. (2016) studied the
interaction between a binary with mass ratio $q=1/10$ and a retrograde CBD.
They considered both circular setups and configurations with an initial
eccentricity $e=0.6$. They point out that because of the large relative
velocities between the secondary and retrograde gas, it is important to use a
sink radius—for the removal of gas at the location of the black hole (see Sec.
3 for more precise definition and details)—that is smaller than the black
hole’s sphere of influence; and determine a sink radius of 2% the secondaries
Roche Lobe radius to be sufficient for a $q=1/10$ secondary. Notably, they
point out that this is because gas near to the secondary black hole can exert
strong gravitational torques that alter the orbit, and too large a sink radius
removes this gas from the domain resulting in erroneous estimates of the
evolutionary effects. Accordingly, they develop an analytic model for the
evolution of binary semi-major axis and eccentricity for small mass ratios
$M_{2}\ll M_{1}$ based on both accretion and gas-dynamical friction from
nearby, but unbound material. The outcome of this model depends on the
steepness of the CBD density profile, but they find that most eccentric
binaries will continue to have their eccentricity pumped $e\rightarrow 1$ by
the retrograde CBD. However, for near-circular binaries they find that
relatively flat density profiles give eccentricity growth, but slightly
steeper profiles with $\rho\propto r^{-n}$ and $n>3/2$ near-circular binaries
stay near-circular while shrinking their semi-major axis. They find this model
to be consistent with their circular and $e=0.6$ simulations, but do not
perform an extensive comparison.

Each of the studies above has noted that because of the lack of binary-disk
resonances, the retrograde CBD is not truncated away from the binary, but
rather extends all the way down until it interacts directly with the binary
orbit at apocenter. Nixon & Lubow (2015) pointed out, however, that while this
lack of Lindblad resonances is true for circular binaries, there exist
components of the potential expansion for eccentric binaries that rotate
retrograde to the binary allowing for resonant torques in eccentric,
retrograde binaries. They find that these torques are generally weak, but for
sufficiently high eccentricity, can become strong enough to drive spiral waves
in the CBD and possibly carve a cavity (similar to prograde solutions). They
corroborate their analytic calculations with a suite of SPH simulations and
find a critical eccentricity $e\sim 0.6$ where retrograde Lindblad resonances
become comparable to viscous torques in the CBD for $\alpha$-viscosity value
$\alpha=0.05$. They observe these resonances to drive spiral density waves
into the CBD and to possibly carve a circumbinary-cavity (although these can
be challenging to resolve in SPH simulations because of the extreme low-
densities associated with circumbinary cavities).

##  3 Numerical Methods

For comparison and increased confidence in results, this study utilizes two
separate grid-based Newtonian hydrodynamics codes: the moving-mesh code Disco
in cylindrical coordinates (Duffell, 2016), and the GPU-accelerated code
Sailfish in Cartesian coordinates (see Westernacher-Schneider et al., 2022,
2023). Both codes solve the the Navier-Stokes equations for viscous, locally
isothermal hydrodynamics in 2-dimensions. We refer the reader to an upcoming
code-comparison paper Duffel et al. (2023) for implementation-specific details
as well as a direct comparison of prograde CBD solutions.

###  3.1 Problem setup

We initialize the disk in two ways: In Sailfish, the vertically-averaged
surface density is initially uniform $\Sigma/\Sigma_{0}=1$ with a Keplerian
rotation profile $v_{\phi}=\sqrt{GM/\mathcalligra{r}\,}$. In Disco, a depleted
central cavity of minimum density $\delta_{0}=10^{-5}$ is imposed on the
otherwise uniform initial surface density profile
$\Sigma/\Sigma_{0}=(1-\delta_{0})e^{-(2.5/r)^{12}}+\delta_{0}$ and the initial
angular velocity of the gas is set by a Keplerian profile plus binary
quadrupole and pressure corrections far from the binary, denoted $\Omega_{0}$
(see, e.g., Miranda et al., 2017), down to $r=a$ at which it saturates to
$\Omega=\Omega_{b}\equiv\sqrt{GM/a^{3}}$ for $r<a$:
$\Omega(r)=\left[\Omega_{0}(r)^{-4}+\Omega_{b}^{-4}\right]^{-1/4}$. As both
sets of initial conditions capture steady state disk solutions far from the
binary, and because we allow the disk to relax for $\gtrsim 500$ binary orbits
before measurement (see § 3.3), these small difference in set up do not have a
significant impact on the results and diagnostics discussed below.

The disks are subject to the time-varying gravitational potential of a binary
with separation $a$, orbital frequency $\Omega_{b}$, total mass $M$, and mass
ratio $q$. $\Sigma_{0}$ is an arbitrary density scaling (=1 in code units),
and $\mathcalligra{r}\,=\sqrt{r^{2}+r_{s}^{2}}$ is the softened radial
coordinate to prevent divergences at the origin (the softening radius is
fiducially set to 5% the binary separation, $r_{s}=0.05\,a$). The disk has
constant kinematic viscosity $\nu=10^{-3}\,a^{2}\Omega_{b}$ which induces
radial inflow in the disk and implies a steady-state accretion rate at
infinity $\dot{M}_{0}=3\pi\nu\Sigma_{0}$.

We consider the disk to be in the thin-disk limit with Mach number
$\mathcal{M}=10$ implying a disk aspect ratio $h/r\sim\mathcal{M}^{-1}=0.1$.
The locally isothermal condition is applied via the sound speed definition
$c_{s}^{2}=-\Phi_{g}/\mathcal{M}^{2}$ with $\Phi_{g}$ the binary potential

$\Phi_{g}=\Phi_{1}+\Phi_{2}\,,\qquad\Phi_{j}=-\frac{GM_{j}}{\mathcalligra{r}\,_{j}},$
(1)

and $\mathcalligra{r}\,_{j}$ the smoothed distance to binary component
$j\in\\{1,2\\}$. For this study we only consider equal mass binaries
$M_{1}=M_{2}=M/2$ (or equivalently, binary mass ratio $q=1$). Moreover, we
take the disk mass $M_{d}\sim\Sigma_{0}a^{2}$ to be much smaller than the
total binary mass $M_{d}\ll M$ such that the timescale for altering the binary
orbit is long compared to the orbital time, and that we may ignore the disks
self-gravity (e.g., the disk Toomre parameter
$Q\sim\mathcal{M}^{-1}(M/M_{d})\gg 1$).

###  3.2 Source terms and boundary conditions

The mass and momentum conservation equations solved in both codes can be
written as

$\displaystyle\partial_{t}\Sigma$
$\displaystyle+\mathbf{\nabla}\cdot(\Sigma\mathbf{v})=S_{\Sigma}$ (2)
$\displaystyle\partial_{t}(\Sigma\mathbf{v})$
$\displaystyle+\mathbf{\nabla}\cdot(\Sigma\mathbf{v}\mathbf{v}+P\mathbf{I}-\mathbf{\tau})=\mathbf{S_{J}}+\mathbf{F_{g}}$
(3)

where $v_{i}$ is the mid-plane fluid velocity, $P=c_{s}^{2}\Sigma$ is the
vertically averaged gas pressure, $\tau$ is the viscous stress tensor, $F_{g}$
is the gravitational force density from the potential in Eq. 1, and
$S_{\\{\Sigma,J\\}}$ are mass- and momentum-sinks respectively. The viscous
stress tensor is calculated as

$\tau_{ij}=\nu\Sigma(\nabla_{i}v_{j}+\nabla_{j}v_{i}-\delta_{ij}\nabla_{k}v_{k})$
(4)

with the covariant derivatives taken in the relevant coordinate system. The
sink terms are included to mimic the accretion of material onto each binary
component because we do not resolve our solutions down to the physical
accretion boundaries. The mass and momentum sinks are given respectively as

$\displaystyle S_{\Sigma}$
$\displaystyle=-\gamma\Omega_{b}\Sigma\sum_{i}w_{i}$ (5)
$\displaystyle\mathbf{S_{J}}$
$\displaystyle=-\gamma\Omega_{b}\Sigma\sum_{i}\mathbf{v}_{i}^{\star}\,w_{i}$
(6)

where $\gamma$ is a dimensionless sink-rate which—unless specified
otherwise—is set to 1, and $\Omega_{b}=\sqrt{GM/a^{3}}$ is the Keplerian
orbital frequency of the binary. $w_{i}$ is a window function defining the
strength of the sink for a gas parcel distance $r_{i}$ away from black hole
$i$ in terms of some characteristic sink radius that is fiducially set equal
to the softening radius $r_{s}$:

$\displaystyle w_{i}=\exp{\left[-(r_{i}/r_{s})^{4}\right]}\ .$ (7)

Different sink behaviors are achieved through calculation of the adjusted gas-
velocity $\mathbf{v_{i}^{\star}}$ associated with each fluid element removed.
For nearly all runs presented in this study we have adopted “torque-free”
(Dempsey et al., 2020)—or _spinless_ —sinks such that no spin angular momentum
is accreted by each binary component—only orbital angular momentum. To
accomplish this, the adjusted velocity is calculated as

$\mathbf{v_{i}^{\star}}=(\mathbf{v}-\mathbf{v_{i}})\cdot\mathbf{\hat{r}_{i}}+\mathbf{v_{i}},$
(8)

where $\mathbf{v}$ is the gas velocity in the inertial frame, $\mathbf{v_{i}}$
is the binary component velocity, and $\mathbf{\hat{r}_{i}}$ is the radial
unit-vector in a coordinate system centered on binary component $i$. In the
appendix we briefly discuss a separate “acceleration free” sink prescription
(also sometimes referred to as a “standard sink” Dittmann & Ryan 2021) in
which $\mathbf{v_{i}^{\star}}=\mathbf{v}$ and the momentum sink reduces to
$\mathbf{S_{j}}=S_{\Sigma}\mathbf{v}$.

In both codes the outer boundary is set so to enforce the steady-state
accretion rate $\dot{M}_{0}$ so that the solution mimics that of an infinite
accretion disk. In Disco this is accomplished with outer ghost-cells fixed to
the setup initial conditions. In Sailfish this is attained via an additional
“buffer” source term that drives the solution back towards the initial
condition. This “buffer” prevents artifacts from the square domain propagating
into the solution. The details of the buffer prescription can be found in §3
of Westernacher-Schneider et al. (2022).

###  3.3 Adiabatic eccentricity variation

The primary calculations presented in this study fix the binary mass ratio
$q=1$ and slowly vary the binary eccentricity from $e_{0}=0$ up to
$e_{f}=0.8$. In the process we measure the rate of change of both the binary
semi-major axis $\dot{a}(e)$ and eccentricity $\dot{e}(e)$ as in D’Orazio &
Duffell (2021) (see also Duffell et al. 2020 for the same procedure performed
on $q$ at fixed-$e$, and Wang et al. 2022 on the disk cooling timescale). The
runs that vary the binary eccentricity are initialized from the output of a
simulation with eccentricity fixed at $e=0$ for 500 binary orbits. The
eccentricity growth is performed linearly over $n$ binary orbits as

$e(t)=e_{0}+\frac{e_{f}}{2\pi\,n}\,t\,.$ (9)

We take $n=5\times 10^{3}\,\text{and}\,10^{4}$ and have verified that further
increasing $n$ does not meaningfully change our results. The time rate of
change of the binary orbital elements $\dot{a}$, $\dot{e}$ are calculated in
the same way as is detailed in D’Orazio & Duffell (2021), except with
additional full accounting for all momentum accreted by the sinks (the mass
deposition effect is included as previously).

The fiducial grid resolution in Sailfish is taken to be $\delta=0.0067a$. In
Disco, the grid is hybrid-logarithmic in radius such that the resolution is
$\delta=0.0127a$ at $r=0.5\,a$, slightly higher inside $r<0.5\,a$, and
decreases as one moves toward the outer boundary at $r=50a$. In Sailfish, the
resolution is constant everywhere and the outer boundary is taken to be
$r=10a$ because initial simulations (and the previous numerical work discussed
in §2) showed that the retrograde CBD becomes near completely axisymmetric at
$r>\text{few}\times a$, with possible exception at large eccentricity to be
discussed later.

##  4 Numerical Results

###  4.1 Orbital evolution: $\dot{a}(e)$ and $\dot{e}(e)$

Figure 1: Time rates of change of binary semi-major axis (_top_) and
eccentricity (_bottom_) from a retrograde circumbinary disk as a function of
binary eccentricity. Results from Sailfish are shown in pink and from Disco in
green. The grey crosses show single, fixed-eccentricity runs computed with
Sailfish for comparison to the results from the eccentricity-varying runs.

Figure 1 shows calculations of $\dot{a}(e)$ (_top panel_) and $\dot{e}(e)$
(_bottom panel_) from the simulations that vary eccentricity adiabatically
from $e=0$ to $e=0.8$. Results from Disco are displayed in green and results
from Sailfish are given in pink. Both codes find that the retrograde CBD
shrinks the binary semi-major axis $\dot{a}<0$ and pumps the binary
eccentricity $\dot{e}>0$ for all eccentricities. The shrinking of the binary
orbit and pumping of the binary eccentricity is in general agreement with
previous studies of binary evolution in retrograde CBD’s, but unlike Nixon et
al. (2011a) and Roedig & Sesana (2014) these simulations show that near-
circular binaries are not in fact driven back towards circularity, but rather
have their eccentricity driven up by the retrograde CBD. The eccentricity
driving rate grows linearly towards a maximum of $de/d\log{M}\approx 3$ at
$e\approx 0.15$ and settle to a near constant growth rate $de/d\log{M}\approx
2.25$ by $e\gtrsim 0.4$. The shrinking of the binary $\dot{a}$ shows almost no
dependence on binary eccentricity and always retains a value
$d\log{a}/d\log{M}\approx-10$. 222Curiously, for a circular binary
$d\log{a}/d\log{M}\approx-10$ corresponds to an accretion eigenvalue (see Lai
& Muñoz 2022) $\ell_{0}=\dot{L}/\dot{M}=-1$. This value, however, is likely a
coincidence because it remains constant with binary eccentricity when changes
to the semi-major axis are no longer defined only by the flow of angular
momentum in the disk.

To sanity check these results we have also run a series of high-resolution
$\delta=0.005\,a$ fixed-eccentricity runs in Sailfish shown as grey crosses.
These runs were done for 2000 binary orbits and we report the rate of change
to the orbital elements over the final 500 orbits. These show near exact
agreement with the Sailfish variable eccentricity run and very close agreement
with the Disco run. At high eccentricities $e>0.7$ we observe a notable growth
in the fluctuations of $\dot{a}$ and $\dot{e}$ in the Sailfish runs alone. We
attribute this growth in variation in one code and not the other to the fact
that the Cartesian grid of Sailfish evolves the fluid linear momenta, and
thus, angular momentum is not perfectly conserved. Disco by contrast
explicitly evolves and conserves the fluid angular momentum by construction.
We have confirmed this effect by lowering the resolution of the Sailfish run
and verifying that the variations become more extreme, that they set in at
lower eccentricity, and that the solutions begin to diverge at large
eccentricities.

Apart from these variations, both simulations show a notable increase in
signal variability at $e\sim 0.55$. We attribute this growth in variability to
the emergence of a two-orbit periodic switching in the circumbinary disk
structure (discussed further in §4.2). This growth in variability also appears
to be coincident with the emergence of two-armed spiral density waves that
extend from the CBD cavity edge suggesting the possible realization of
retrograde circumbinary resonances in accordance with the predictions of Nixon
& Lubow (2015).

###  4.2 Disk morphology

Figure 2: Snapshot of the steady-state flow pattern for a circular binary in a
retrograde CBD. The arrows show the direction of the fluid velocity. The
binary is orbiting counter-clockwise. The minidisks are retrograde in
accordance with the bulk flow. The minidisks have persistent wakes that also
feed low-angular momentum material into the standing-bridge between the binary
components. Figure 3: Density maps at $e=\\{0.0,0.1,0.3,0.5,0.7\\}$ from the
eccentricity varying runs with both Disco (_top row_) and Sailfish (_middle
row_) as well as the fixed-eccentricity runs (_bottom row_), all taken at
pericenter.

A single snapshot for a circular binary ($e=0$) after 2000 orbits is shown in
Figure 2 that highlights many of the general features of retrograde solutions.
The colormap shows the log-density of the accretion flow, and the overlayed
arrows show the direction of the gas velocity. Similar to previous studies,
the circumbinary disk extends all the way to the binary orbital radius, and
remains almost completely axisymmetric everywhere except for the inner-most
$r<a$ (in contrast with the large, lopsided cavities observed in circular,
prograde scenarios). The binary carves a low-density cavity inside of its
orbital radius, but does not entirely expel it of material, and rather is
always orbitting in some ambient medium. The binary components capture
retrograde material into circum-single “minidisks” that retain opposite
angular momentum to the orbital angular momentum of the binary. Rather than
becoming perfectly circular, the minidisks have bow-shock-like structures on
their leading edge from the ram-pressure of incident counter-rotating
material, as well as trailing “tadpole wakes” of low angular momentum material
that falls almost radially toward the binary barycenter in both directions.
This creates a persistent standing “retrograde-bridge” of material between the
binary components. This retrograde-bridge is quasi-steady—it does wobble
slightly—in circular retrograde solutions, but becomes more transitory for
eccentric binaries as will be demonstrated later. We also note that all future
surface density plots will use the same colorbar indicated in Figure 2 and the
same $[-2a,2a]$ Cartesian stretch (unless otherwise indicated).

Figure 4: Phase dependent disk morphology for an $e=0.3$ orbit at four values
of the true anomaly $f=\\{0,\pi/2,\pi,3\pi/2\\}$. Overlayed arrows show the
direction of the fluid velocity. Figure 5: Phase dependent disk morphology for
two contiguous $e=0.7$ orbits at four values of the true anomaly showing the
two-orbit periodicity in the flow behavior. The first row illustrates the
first orbit in this two-orbit behavior, and the second row shows the second.
Overlayed arrows show the direction of the fluid velocity.

Figure 3 displays snapshots of the disk $\log$-surface-density for
eccentricities $e=\\{0.0,0.1,0.3,0.5,0.7\\}$ with Cartesian extent $[-3a,3a]$.
The _top two rows_ sample from the eccentricity varying runs while the _bottom
row_ uses results of the fixed-eccentricity runs. All snapshots are taken at
binary pericenter. There is generally very good agreement in disk morphology
between the eccentricity varying runs (Disco sweep and Sailfish sweep; _rows 1
and 2_) and with the fixed-eccentricity runs (_row 3_). The only discernible
differences are in slight variations to the minidisk density and the exact
shape of the retrograde-bridge between the binary components. We investigate
the sensitivity of these wakes and bridges to numerical choices in § Appendix
A: Dependence on Numerical Parameters.

The time-dependence of the disk evolution has three separate regimes: (_i_)
steady-state in the co-rotating frame near circularity, (_ii_) the driving of
axisymmetric density waves for eccentric binaries $0.025\lesssim e\lesssim
0.55$, and (_iii_) the forcing of non-axisymmetric spiral density waves at
large eccentricities $e\gtrsim 0.55$. Circular binaries (_i_) with $e\lesssim
0.025$ resemble Figure 2 and show no phase dependence in their evolution. The
systems in regime (_ii_) acquire a “breathing-mode” that drives axisymmetric
density waves into the CBD. This process is shown in Figure 4. At pericenter,
the binary has carved a fully-depleted, axisymmetric cavity. As it approaches
apocenter (true anomaly $0<f<\pi$; here-on we refer to this as “binary
waxing”), the minidisks circularize in the relative-vacuum of the cavity, but
CBD material encroaches on the cavity as the orbit expands and tidally
redirects material through the domain center. This redirected material from
each component collides forming the retrograde-bridge, and by binary apocenter
the cavity has been replenished with low-density gas. The ram-pressure of this
ambient material disrupts the minidisks compressing the leading edge and
stripping off a spiral wake. As the binary accelerates towards pericenter
($\pi<f<2\pi$; “binary waning”), it expels the material inside of its orbit
driving an axisymmetric density ring into the CBD. This ring can be seen
forming in the final panel of Figure 4 ($f=3\pi/2$), as circularized in the
first panel (at pericenter), and propagating through the CBD as a sound wave
in the second and third panels. We illustrate $e=0.3$ as an example, but all
eccentricities $0.025\lesssim e\lesssim 0.55$ were observed to have the same
phase-dependent behavior. We note that the driving of an axisymmetric density
wave sets in at very small eccentricity $e\approx 0.025$, but the carving of a
depleted cavity around pericenter, however, occurs more gradually. A fully-
depleted cavity is effectively formed once per orbit starting at $e\approx
0.1$.

Binaries in regime (_iii_) become qualitatively different as they acquire a
two-orbit periodicity in the phase-dependence of their flow. These two orbits
are illustrated in Figure 5. Similar to regime (_ii_), at the first pericenter
(selected arbitrarily, but shown at the top-left panel of Figure 5) the binary
has carved a fully-depleted, mostly axisymmetric cavity with the exception of
two weak spiral arms extending from the cavity wall. In the first orbit (_top
row_), during binary waxing, the minidisks similarly circularize and material
is tidally re-oriented into the depleted cavity. At apocenter, the cavity is
replenished with gas, the retrograde-bridge has formed, and the minidisks are
strongly perturbed from the collisions with encroaching CBD material. In the
approach to pericenter, the accelerating binary again expels gas from within
its orbit, but instead of driving an axisymmetric density wave (as in the
lower $e$ case), it propels two spiral density waves that propagate into the
disk. Because of this, at second-pericenter (bottom left panel), the binary
has carved a depleted cavity, but the cavity wall is no longer circular as it
is dominated by the $m=2$ spirals. In the second orbit (_bottom row_), the
same processes occur, but the non-axisymmetric cavity at pericenter is less
efficient at refilling the cavity during binary waxing. The minidisks are
less-perturbed at second-apocenter (in the bottom row) as a result of these
lower cavity densities; and while the binary still creates a two-armed spiral
structure upon expelling its intra-orbit material during binary waning, the
lower densities weaken the response, and the resulting pericenter cavity is
left essentially circular. This two orbit periodicity also appears as power in
the accretion rate time series and periodogram (see Figures 10 and 11).

We emphasize that at all eccentricities, these 2D Newtonian, isothermal
hydrodynamics simulations show the formation of binary minidisks and a
retrograde-bridge. The presence of these minidisks—and their visible
asymmetries throughout the binary orbit—preempt the importance of including
the binary and the central most $r<a$ regions of the accretion flow in order
to accurately ascertain the gravitational forces on the binary and its
resultant orbital evolution.

###  4.3 Gravity vs. accretion

A major component of most previous studies of orbital evolution from
retrograde CBD’s was that the evolution is dominated by the direct accretion
of both mass and angular momentum due to collisions between retrograde fluid
elements and the binary at apocenter; with the exception of Amaro-Seoane et
al. (2016) who pointed out that gas near the binary orbit can exert very
strong gravitational forces before physically accreting. In order to quantify
the relative effects of gravitational and accretion forces, we separated both
$\dot{a}$ and $\dot{e}$ into their components due solely to gravitational
forces and those due only to the accretion of mass and momentum. Figure 6
shows this decomposition by plotting the total $\dot{a}$ (_purple_) and
$\dot{e}$ (_orange_), and the components due to gravitational forces alone
$\dot{a}_{\rm grav}$, $\dot{e}_{\rm grav}$ (illustrated as _light-purple_ and
_light-orange_). The top panel shows this decomposition for the Disco run and
the bottom panel shows it for the Sailfish simulation. The primary conclusion
from these breakdowns is that—contrary to previous studies and
models—gravitational forces are the dominant component of the binary’s orbital
evolution for both $\dot{a}$ and $\dot{e}$ as their components due to gravity
alone near perfectly describe the full orbital evolution curves.

The effect from the physical accretion of mass and momentum is almost entirely
negligible for the eccentricity evolution in both codes and only represents a
$\sim 10\%$ contribution to the change in semi-major axis. We posit that
this—just as was pointed out in Amaro-Seoane et al. (2016) for $q=1/10$
binaries—is because material captured from the CBD at binary apocenter is not
directly added to the binary via accretion, but rather, is transferred onto
orbits around the individual binary components. Moreover, these circum-single,
“minidisk” structures are not symmetric around each binary component, but
instead have small wakes that trail each component exerting gravitational
torques on and removing energy from the binary orbit. Even larger wakes and
non-symmetric features are formed when the minidisks are partially disrupted
from the impact with the CBD at each apocenter passage as discussed in §4.2.

To compare the extent to which gravitational forces from this intra-orbit
material (material within $r\leq a$) dominate over those from the outer CBD
($r>a$), Figure 7 shows the average magnitudes of the unit-less torque
$|(\mathbf{r_{b}}\times\mathbf{f_{g}})/\sqrt{GMa(1-e^{2})}|$ and power
$|(\mathbf{f_{g}}\cdot\mathbf{v_{b}})/(GM/2a)|$, with $r_{b}$ the binary
separation, exerted on the binary from each of the fixed-eccentricity
comparison runs. The average torque is shown by crosses and the average power
by circles. The components from the intra-orbit material with $r\leq a$ are
given in red and those from the outer CBD in blue. At all eccentricities both
the gravitational torque and power exerted on the binary are dominated by the
intra-orbit material. Further, the components of the average torque and power
from the outer-CBD $r>a$ are almost entirely negligible at eccentricities
below the threshold for retrograde resonances to drive spiral density waves
into the CBD, $e\lesssim 0.55$. At $e=0.7$ we once again see the effect of
persistent spiral density waves in the CBD as the average torque from the
outer-disk is no longer negligible—and instead accounts for approximately 25%
of the total torque on the binary.

Figure 6: Orbital evolution from Disco (_top_) and Sailfish (_bottom_)
separated into the total effect and the component from gravitational forces
alone. The semi-major axis and eccentricity evolution are both dominated by
the gravitational pull of the gas. Figure 7: Magnitude of gravitational torque
and power for fixed-eccentricity runs with Sailfish. We see that gravitational
forces from the inner-most region of the flow dominate the binary evolution at
all eccentricities. At high-eccentricities $e\gtrsim 0.5$, the binary drives
spiral wakes into the outer-CBD, leading to growth in the torque component
from $r>a$; although this component remains subdominant by a factor of a few.

###  4.4 Eccentricity driving

Figure 8: Averaged change in eccentricity $\Delta e$ per binary phase (mean
anomaly) at four eccentricities $e=\\{0.02,0.1,0.3,0.7\\}$. The top
panel—representative of small eccentricities $e\sim\mathcal{O}(10^{-2})$ in a
_wake-driven_ mode—shows different phase-dependent eccentricity driving
behavior than all other eccentricities $e\gtrsim 0.1$ in a _cavity-mediated_
mode. Figure 9: Time averaged strength of azimuthal density modes with mode
numbers $m=\\{1,2,3,4\\}$. We see the emergence of only even-$m$ azimuthal
modes in accordance with linear theory of binary-disk resonances for equal
mass binaries. The strong growth of the $m=2$ mode at $e=0.7$ reflects the
observed $m=2$ spiral density waves in the associated surface density plots
and is evidence for the emergence of retrograde resonances. The non-zero
even-$m$ mode strengths at $e<0.55$ are likely due to the temporary non-
axisymmetry of the CBD around binary apocenter once every orbit.

A primary component of past modelling for binaries accreting from retrograde
disks was the assumption that eccentricity pumping is focused around binary
apocenter (and damping—if included—occurs around pericenter). We examine this
hypothesis by plotting the average change in binary eccentricity $\Delta
e(\theta)=\langle\dot{e}\delta t\rangle_{\theta}$ at each binary phase
(orbital mean anomaly, $\theta$) in Figure 8 with $\dot{e}$ the instantaneous
rate of change, $\delta t$ the timestep, and $\langle\rangle_{\theta}$
denoting the average at mean anomaly $\theta$. $\Delta e$ includes the fact
that the binary spends different amounts of time at each orbital phase, and
integrating the curve over all phases would yield the average per-orbit change
in eccentricity. This tells us the relative contribution—on average—of each
binary phase to the net change in binary eccentricity per orbit. We show the
total effect from both gravitational forces and the accretion of mass and
momentum as solid black curves as well as the decomposed contributions from
gravitational torque (_orange dashed_ curves) and gravitational power (_green
dash-dotted_ curves) given respectively as

$\displaystyle\dot{e}_{\rm
power}=-\frac{\dot{\epsilon}}{\epsilon}\left(\frac{1-e^{2}}{2e}\right)$
$\displaystyle=\frac{\mathbf{f_{g}}\cdot\mathbf{v_{b}}}{GM/a}\left(\frac{1-e^{2}}{e}\right)$
(10) $\displaystyle\dot{e}_{\rm
torque}=-2\frac{\dot{\ell}}{\ell}\left(\frac{1-e^{2}}{2e}\right)$
$\displaystyle=-\frac{\mathbf{r_{b}}\times\mathbf{f_{g}}}{\sqrt{GMa(1-e^{2})}}\left(\frac{1-e^{2}}{e}\right),$
(11)

such that $\dot{e}_{\rm grav}=\dot{e}_{\rm power}+\dot{e}_{\rm torque}$ (see
D’Orazio & Duffell 2021 Equations 6-8 for more detailed discussion of these
terms). These phased eccentricity driving curves are shown for four values of
eccentricity $e=\\{0.02,0.1,0.3,0.7\\}$ to encapsulate each of the three
morphologic disk regimes: $e=0.02$ in regime (_i_), $e=0.1,\,0.3$ in regime
(_ii_), and $e=0.7$ in regime (_iii_) (see §4.2). Of primary note,
eccentricity pumping is found to occur during binary waxing (pericenter
$\rightarrow$ apocenter) for all eccentricities, peaking near a mean anomaly
of $\pi/2$ (or equivalently at $P_{b}/4$) at small eccentricities and shifting
nearer to pericenter with growing $e$. Eccentricity is correspondingly damped
during binary waning (apocenter $\rightarrow$ pericenter), and the effect at
pericenter and apocenter is minimal at all eccentricities considered. The net
integrated effect over one full orbit yields the $\dot{e}_{\rm grav}$ curves
shown in Figure 6.

One major distinguishing characteristic between regime (_i_) and regimes (_ii_
, _iii_) is that in regime (_i_) the binary is always orbiting through ambient
intra-orbit material and as a result forms persistent tadpole wakes that trail
each binary component. In regimes (_ii_ & _iii_), the binary expels all intra-
orbit material once per orbit during its waning phase such that it evolves in
a fully-depleted cavity for at least half its orbit until material is once
again redirected within the orbital radius around apocenter (see Figures 4 &
5). We identify two separate eccentricity driving modes associated with these
disk morphologic regimes: (1) a _wake-driven_ mode when the binary evolves
through persistent intra-orbit material and never carves a fully depleted
cavity (regime _i_ ; see Figure 2). The associated eccentricity driving is
shown in the top-panel of Figure 8 for $e=0.02$ where torques from the tadpole
wakes pump eccentricity at all binary phases. The associated power from the
wakes always acts with opposite effect, damping the binary’s eccentricity; but
the torque-mediated pumping wins out on aggregate.

The second eccentricity driving mode, (2) a _cavity-mediated_ mode, occurs
once the binary has begun carving a fully depleted cavity and orbits in
relative vacuum for a significant portion of its orbit (regimes _ii_ & _iii_).
This mode is observed in the $e=0.1,\,0.3,\,0.7$ panels of Figure 8 where the
total effect on $\Delta e$ is dominated by, and near perfectly tracks, the
power contribution. Once the binary begins carving a depleted cavity,
gravitational power continues to damp eccentricity during binary waning as the
binary expels its intra-orbit material. However, the power effect switches to
eccentricity pumping during waxing when the binary and its minidisks evolve
through relative vacuum. This transition occurs gradually, but has mostly set
in by $e\approx 0.1$. The effect of torque on $\Delta e$ during this mode is
negligible at all phases, except near apocenter when the binary has refilled
its cavity and temporarily evolves through ambient material with the
associated tadpole wakes. This effect, though, becomes less significant with
growing eccentricity.

We attribute the initially small eccentricity driving at small eccentricities
$e\sim\mathcal{O}(10^{-2})$ in Figure 1 to the small net difference between
the eccentricity changing effects of gravitational torque and power from the
wakes before the binary begins to carve a depleted cavity once per orbit. The
slow linear growth in $\dot{e}$ with $e$ reflects the gradual transition to a
cavity-carving binary; and once the binary is effectively depleting its orbit
of material around pericenter—despite the amplitude of all eccentricity
modifying effects decreasing with $e$—the net balance remains relatively
constant with $e$.

###  4.5 Retrograde resonances

As discussed in §2, Nixon & Lubow (2015) pointed out that eccentric binaries
have components of the potential expansion that rotate retrograde to the
binary admitting retrograde Lindblad resonances; and their simulations find
that these resonances become strong enough to drive persistent spiral density
waves into the disk at eccentricities $e\geq 0.6$. The simulations presented
in this study also suggest evidence for the emergence of retrograde resonances
at eccentricities $e\gtrsim 0.55$. Figure 3 shows the first signs of non-
axisymmetric density waves at $e=0.5$, and at $e=0.7$ reveals the presence of
two-armed, $m=2$ spiral density waves originating at the CBD inner edge and
propagating into the CBD.

In traditional linear analyses of resonant torques on a CBD, the forces
associated with a harmonic mode $m$ of the gravitational potential will drive
spiral density waves into the CBD with the same azimuthal mode number.
Moreover, for equal-mass binaries, odd components of the potential expansion
disappear, so only even-$m$ resonant torques should act on the CBD333and the
presence of odd azimuthal modes—the $m=1$ mode in particular—has been
suggested as evidence for non-resonant effects in the cavity-carving process
for equal mass binaries (Mahesh et al., 2023).. To check, albeit indirectly,
for the presence of such resonances, we calculate timeseries of the azimuthal
density mode $m=\\{1,2,3,4\\}$ as

$\Psi_{m}(t)=\int_{a}^{r_{\rm
out}}dr\,\int_{0}^{2\pi}rd\phi\,\Sigma(t)\,e^{im\phi}.$ (12)

Figure 9 shows the time averaged strength of these azimuthal density modes
$\bar{\Psi}_{m}$ from the fixed-eccentricity runs. We observe the presence of
an $m=2$ azimuthal mode consistent with the observation of $m=2$ spiral
density waves at $e>0.55$. However, this mode appears to emerge smoothly with
growing eccentricity as opposed to an abrupt appearance at large $e$. The next
even mode, $m=4$, does appear only at large eccentricity $e\geq 0.5$. The
strength of the odd-modes is consistent with zero for all eccentricities. This
non-existence of odd azimuthal modes, the presence and growing strength of the
even modes, and the observed spiral density waves at large $e$, are consistent
with the emergence of retrograde resonances in our simulations. However, these
features could also emerge from non-resonant interactions with the time
varying potential of the equal mass binary and conclusively disentangling the
two would require a more detailed analysis.

###  4.6 Accretion rate

To determine periodicity structure in the accretion rate, we compute a 2D
periodogram of the accretion rate time series as measured onto both binary
components. To do so we utilized the entire $10^{4}$ binary orbit accretion
rate time series from Disco, spanning from $e=0-0.8$ (similar results are
found with Sailfish, see below). As in Duffell et al. (2020), we convolve a
Gaussian window in time (and so also binary eccentricity) with the inner
product of the accretion rate and Fourier vector with angular frequency
$\omega$ in frequency space,

$\mathcal{P}(e,\omega)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\int^{t(e_{f})}_{t(e_{0})}{\mathrm{e}^{-\frac{1}{2}\frac{\left(t(e)-\tau\right)^{2}}{\sigma^{2}}}\dot{M}(\tau)\mathrm{e}^{i\omega\tau}d\tau},$
(13)

where eccentricity dependence comes through $t(e)$, the inverse of Eq. (9). We
choose $\sigma=30P_{b}$ and compute the norm of Eq. (13) over a $300\times
300$ grid of values of $e$ ranging from $0.0$ to $0.8$ and $\omega$ ranging
from $0.1\Omega_{b}$ to $2.5\Omega_{b}$, corresponding to variability
timescales between $(0.1-2.5)P_{\rm{orb}}$.

Figure 10 plots contours of the log-base-10 power computed over this grid. The
most prominent power is concentrated in narrow bands at the orbital time and
its higher frequency harmonics. In addition, a wider band is concentrated at
twice the binary orbital period (half the orbital frequency), appearing for
eccentricities above $e\approx 0.55$. The latter appearance of a twice-orbital
periodicity is notably coincident with the appearance of retrograde resonances
as argued in § 4.5. The twice-orbital periodicity derives from the alternating
process of cavity-clearing described in § 4.2.

In Figure 11 we further explore the time series accretion rates for
representative binary eccentricities. In this figure we plot the Disco results
in thick black lines and also the Sailfish results in grey lines for
comparison. We first describe the Disco accretion-rate time series. For $e=0$,
the accretion rate is steady. As $e$ increases from zero, the accretion rate
becomes strongly modulated at the orbital period with peaks just following
apocenter, and with amplitude growing with $e$ until $e\sim 0.2$. From $e\sim
0.2-0.5$, the amplitude of accretion rate modulations saturates to $\sim 50\%$
of the mean, and the time of peak accretion rate drifts from just after
apocenter towards just before pericenter with increasing $e$. At $e\gtrsim
0.3$, a second peak in the accretion rate develops, occurring just after
apocenter and before the other peak, which, for $0.3\lesssim e\lesssim 0.5$,
occurs just before pericenter. For $e\gtrsim 0.55$, the twice-per-orbit
periodicity appears , and the amplitude of modulations begins to again grow
with $e$. For $e=0.6$ and $e=0.7$, the twice-orbital-period periodicity
manifests as the double peaked accretion rate modulation becoming a factor of
$\sim 2$ higher every other orbit, but otherwise similar in shape, due again
to the alternating cavity structure seen in Figure 5 and described in § 4.2.
For $e=0.8$, the higher accretion-rate modulation occurring every other orbit
is punctuated by a large, factor of $3.5$ accretion rate spike just following
apocenter.

The Sailfish accretion-rate time series are similar to the Disco results in
periodicity structure and in the magnitude of accretion rate modulations, but
exhibit a number of differences. Primarily, the double peaked structure does
not manifest until $e\gtrsim 0.6$ for the Sailfish runs. Rather, at
intermediate eccentricities, what appears as a double peaked structure in the
Disco runs, appears instead as a low accretion rate kink in the grey Sailfish
time series. At high eccentricities Sailfish also exhibits large spikes
following apocenter, but starting at $\sim 0.7$ compared to $e\approx 0.8$ for
Disco. We note that the Sailfish accretion rate time series in the final
$e=0.8$ panel of Figure 11, is likely exhibiting spurious behavior due to the
non-explicit conservation of angular momentum. Hence, while periodicity
timescales and the magnitude of accretion-rate variations are robust between
the codes, small differences in the shape of the accretion rate time series
such as the observed double peaked structure are not, and caution should be
taken in applying these to observable features of accreting binaries (in
addition to further complications in conversion of accretion rates to
luminosities).

$\begin{array}[]{c}\includegraphics[width=433.62pt]{Figures/Pdgm_Disco_q1_retro.png}\end{array}$

Figure 10:  Power in the total binary accretion-rate periodogram computed via
Equation 13. The dominant power at all eccentricities is at the orbital period
and its higher frequency harmonics. For eccentricities $e\gtrsim 0.55$ power
also emerges in a wide band around two-times the orbital period.

$\begin{array}[]{ccc}\includegraphics[scale={0.33}]{Figures/e0p0_SDcompare_MdotTS_pa}&\includegraphics[scale={0.33}]{Figures/e0p1_SDcompare_MdotTS_pa}&\includegraphics[scale={0.33}]{Figures/e0p2_SDcompare_MdotTS_pa}\\\
\includegraphics[scale={0.33}]{Figures/e0p3_SDcompare_MdotTS_pa}&\includegraphics[scale={0.33}]{Figures/e0p4_SDcompare_MdotTS_pa}&\includegraphics[scale={0.33}]{Figures/e0p5_SDcompare_MdotTS_pa}\\\
\includegraphics[scale={0.33}]{Figures/e0p6_SDcompare_MdotTS_pa}&\includegraphics[scale={0.33}]{Figures/e0p7_SDcompare_MdotTS_pa}&\includegraphics[scale={0.33}]{Figures/e0p8_SDcompare_MdotTS_pa}\end{array}$

Figure 11:  The total accretion rate measured onto the binary components from
the Disco sweep (black) and for comparison, the Sailfish sweep (grey). In each
panel, accretion rate time series are shown for 10 orbits starting from where
the eccentricity sweep (Equation 9) reaches the binary orbital eccentricity
$e$, indicated in each panel. Vertical dotted blue (red) lines denote
pericenter (apocenter). While there are small differences between codes,
especially at high eccentricities, a few robust qualitative features persist:
For lower eccentricities ($e\lesssim 0.55$) the accretion rate times series is
modulated at the orbital period, with a double peaked (for Disco) or kinked
(for Sailfish) structure at intermediate eccentricities. For $e\gtrsim 0.55$,
a twice-orbital period modulation arises, as indicated by Figure 10. For the
highest eccentricities probed in this study, a large accretion rate spike is
induced following apocenter.

##  5 Discussion

###  5.1 Comparison with previous studies

The results presented in this study generally agree qualitatively with
previous numerical and analytic works on retrograde binary accretion, but
there exist pertinent differences. Analytic modelling of retrograde accretion
scenarios from Nixon et al. (2011b); Roedig & Sesana (2014); Schnittman &
Krolik (2015) were all generally founded upon the assumption that binary
evolution is dominated by the physical accretion of counter-rotating gas at
binary apocenter (and sometimes at pericenter). Generally, these models
accurately predict that the binary will shrink its semi-major axis and grow
its eccentricity at nearly all eccentricities. However, Nixon et al. (2011a)
forecasted that at small eccentricities $e<h/r$, the binary would actually
shrink its eccentricity and remain near circular due to opposing accretion
effects at binary apocenter and pericenter. This study does not corroborate
this prediction and finds that the binary grows its eccentricity always. On
the contrary, the modelling of Schnittman & Krolik (2015) anticipates binary
eccentricity growth at all eccentricities, but quantitatively predicts the
effect to be largest at small eccentricities and to decline as eccentricity
grows. The results in Figure 1 show the opposite. Eccentricity growth is
smallest at small $e$ and grows to a peak value at $e\sim 0.15$ after which it
stabilizes to $de/d\log{M}\sim 2$.

Figure 12: Upper- and lower-bounds (_left_ and _right_ respectively) to the
eccentricity achieved by MBHB’s accreting from a retrograde disk determined
according to Equation 15. The light shaded region in the upper-right of each
panel shows the region where a binary of that mass and initial semi-major axis
does not fit in a gravitationally stable disk. The dark grey shaded region in
the bottom right corner shows where the binary does not fit outside the
component ISCO’s for a given mass. The dashed black lines illustrate initial
semi-major axes $a_{0}$ (in parsec) associated with the initial binary periods
$P_{0}$ shown on the y-axes.

These studies (along with Bankert et al. 2015; Nixon & Lubow 2015) have
additionally either verified or motivated their models with particle- and
grid-based numerical studies of accreting retrograde systems. These
simulations generally substantiated the postulate that the orbital evolution
was dominated by the physical capture of gas at binary apocenter, and that
gravitational forces were a subdominant effect. We have tested this directly
in our full-domain solutions and found that the opposite is true: the binary
orbital evolution is almost entirely explained by gravitational forces on the
binary, and the physical accretion of mass and momentum is a secondary effect.
Moreover, the dominant phase for driving eccentricity occurs during orbital
waxing, between pericenter and apocenter, and the effect at apocenter is
small. Similar to prograde scenarios, the evolution of the binary is
controlled by time-dependent, non-axisymmetric features in the disk morphology
that exert strong gravitational forces, and we have demonstrated that the
strongest forces come from the intra-orbital material that passes within the
semi-major axis of the binary $r<a$ . We hypothesize that this had previously
been missed because older simulations did not include the inner-most regions
of the accretion flow either by construction (Bankert et al., 2015) or because
of the specifics of their numerical scheme and setup (Nixon et al., 2011a;
Roedig & Sesana, 2014; Nixon & Lubow, 2015). This is akin to using a sink-
prescription that is too large and aggressive (see Appendix) such that it
washes away the complex interplay between intra-obit material and the
circumbinary flow as well as the non-trivial exchanges of energy and angular
momentum between this material and the binary (as was noted for small mass-
ratio binaries $q<0.1$ by Amaro-Seoane et al. 2016).

###  5.2 Observational implications

The presence of circumbinary disks can leave observable imprints in both
electromagnetic and gravitational wave radiation from accreting binaries, and
the effects of such accretion phases leave a lasting impact on populations of
binaries that have undergone gas-mediated phases. Such effects have been
discussed in detail for prograde accretion scenarios (e.g., Farris et al.,
2015; Ryan & MacFadyen, 2017; Bortolas et al., 2021; Major Krauth et al.,
2023), and many retrograde effects have been discussed in previous works (c.f.
Schnittman & Krolik, 2015). We expand on these by discussing the implications
of this study in the context of recent prograde results.

_Disk-mediated decay and GWs_ — In prograde solutions, it is possible for
accretion-mediated phases of binary evolution to both expand the binary orbit
and to facilitate binary inspiral. Namely, near-circular-orbit binaries tend
to circularize and expand their orbits due to accretion and interaction with a
CBD for disk scale-heights $h/r\gtrsim 0.05$; however, prograde binaries with
initial eccentricity $e\gtrsim 0.08$ evolve towards an equilibrium
eccentricity of $e_{\rm eq}\sim 0.4$ (or some value $0.25\lesssim e_{\rm
eq}\lesssim 0.5$ for $q<1$ Siwek et al. 2023b) where they shrink their semi-
major axis at a rate $\mathcal{A}\equiv d\log a/d\log M$ that is order unity
(Zrake et al., 2021; D’Orazio & Duffell, 2021).

For retrograde accretion, we can describe the orbital evolution with a simple
fitting formula,

$\displaystyle\dot{e}$ $\displaystyle=$
$\displaystyle\mathcal{E}\frac{\dot{M}_{0}}{M}\begin{cases}e&e<e_{*}\\\
e_{*}&e\geq e_{*}\end{cases}$ $\displaystyle\frac{\dot{a}}{a}$
$\displaystyle=$ $\displaystyle\mathcal{A}\frac{\dot{M}_{0}}{M}.$ (14)

For the remaining discussion we choose $\mathcal{E}=30$, $e_{*}=0.1$, and
$\mathcal{A}=-10$. Note, however, that a range of $e_{*}\sim 0.1-0.2$ and
$\mathcal{E}e_{*}=2.5-3.5$ are consistent with the Disco and Sailfish results
in Figure 1. Hence, in contrast to the prograde results, and in accordance
with previous retrograde accretion studies, retrograde accretion scenarios
facilitate binary decay at all eccentricities, and do so at a rate
$\mathcal{A}\approx-10$; a factor of $2-3$ faster than prograde disks at
$e_{\rm eq}$. For black holes accreting at their Eddington rate, this implies
a retrograde gas-mediated inspiral timescale
$\tau_{a}=a/\dot{a}=-\mathcal{A}^{-1}M/\dot{M}\sim 4.5\,{\rm Myr}$ where we’ve
used $M/\dot{M}\approx 4.5\times 10^{7}\,{\rm yr}$ as the Eddington-limited
mass doubling time of an accreting black hole—or Salpeter time—with an
accretion efficiency 0.1. This is shorter than the expected $10-100\,{\rm
Myr}$ lifetimes of quasars (Martini, 2004).

In contrast with prograde solutions, binaries accreting from retrograde disks
have their eccentricity pumped at all eccentricities. Our $\dot{e}(e)$
solution is nearly linear for $e\lesssim 0.1$, with an estimated form
$30e\,\dot{M}/M$, and approximately constant for $e\geq 0.1$. Therefore, for
initial eccentricities $e_{0}\lesssim 0.1$, the binary will grow its
eccentricity exponentially, with an e-folding timescale of $(\sim
30\dot{M}/M)^{-1}$, quickly driving initially small eccentricities into the
constant $\dot{e}$ regime, where eccentricity grows at approximately twice the
mass doubling rate. Hence, for any initial binary eccentricity, retrograde
circumbinary disks will drive eccentricities $e\rightarrow 1$ in approximately
half of a Salpeter time, $\sim 20$ Myr. Comparing the mutual evolution of
$\dot{a}$ and $\dot{e}$, this implies that $a$ will decrease by 5 e-foldings
$a_{0}\mathrm{e}^{-5}$ in the time required for $e\rightarrow 1$.

A retrograde disk will simultaneously pump binary eccentricity and shrink the
semi-major axis until the effects of energy and angular momentum loss by GWs
begin to dominate at either high-$e$ or small-$a$. In the large eccentricity
limit, there exists some eccentricity $e_{\star}$ at which GWs will begin to
damp binary eccentricity at a rate faster than the disk can pump it. We can
estimate a _lower-bound_ for the eccentricity achieved from retrograde
accretion before GWs take over by finding $e_{\star}$ such that

$\dot{e}_{GW}(M,a_{0}\mathrm{e}^{-5},e_{\star})=-\dot{e}_{CBD}(e_{\star})$
(15)

as a function of binary mass $M$ and initial separation $a_{0}$ (assuming
$q=1$). This lower bound is likely the more realistic scenario, but we also
compute an _upper-bound_ for binary eccentricity at the onset of a GW
dominated regime by using $a_{0}$ instead of $a_{0}e^{-5}$. This ignores the
commensurate shrinking of semi-major axis during eccentricity pumping444We
note that the quantities $a(1+e)$ and $a(1-e)$ are not conserved and are
always decaying.. These upper- and lower-bounds on the disk-driven
eccentricity are shown as contours of $\log_{10}(1-e_{\star})$ in $a_{0}-M$
space in Figure 12 (the _left_ and _right_ panels respectively). The light
shaded region in the upper right corner denotes where $a_{0}$ would not fit
into a gravitationally stable, steady-state thin disk (using Equation 16 in
Haiman et al. 2009), and the dark shaded region in the bottom right corner
illustrates where a binary of given mass is within the component ISCOs. The
black-dashed lines show initial semi-major axes associated to the initial
binary period $P_{0}$.

Massive binaries that undergo periods of prograde accretion before merging in
the LISA band are predicted to be driven to an equilibrium eccentricity
$e_{\rm eq}=0.4$. Binaries that enter a GW-dominated regime at $e_{\rm eq}$
will retain only a sub-percent eccentricity very near the LISA detection
threshold (Cuadra et al., 2009; Zrake et al., 2021). Massive binaries that
have their eccentricity pumped to near unity by retrograde accretion, however,
may retain significantly more eccentricity upon entering the LISA band and are
comparatively more likely to lie above the detector’s eccentricity detection
threshold (Garg et al., 2023). Assuming the lower-bound as the more realistic
scenario in Figure 12, we still observe that phases of retrograde accretion
could drive very large eccentricities $e_{\star}\gtrsim 0.9$ into massive
binaries that are eventually detected by LISA, if accretion commences at large
enough separations $a_{0}\sim 10^{-2}$pc. Moreover, because we observe
retrograde minidisks, binaries that have undergone retrograde accretion will
acquire spins that are counter-aligned with the binary orbital angular
momentum. Thus, a signpost of such retrograde accretion phases may be the
measurement of notable eccentricity in conjunction with negative effective
spin parameters $\chi_{\rm ef}<0$ in the GW signal analysis.

Gravitational wave emission at such large eccentricities will also shift GW
power from the standard, circular $2f_{b}$—with $f_{b}$ the binary
frequency—to higher frequencies $f=f_{b}(1+e)^{1/2}/(1-e)^{3/2}$. Rapid
environment driven coalescence will additionally diminish GW power at
frequencies where such effects are active. For the most massive binaries
$M>10^{8}\,M_{\odot}$, like those posited to source the low-frequency
gravitational wave background (Agazie et al., 2023b; Antoniadis et al., 2023a;
Reardon et al., 2023; Xu et al., 2023), the effect of a retrograde accretion
phase that both pumps eccentricity and rapidly shrinks the binary orbit would
be to reduce GW power at the low-frequency end of the background, as was
tentatively observed (Agazie et al., 2023a; Antoniadis et al., 2023b).
However, Figure 12 suggests that eccentricity pumping for such massive
binaries may not be as effective as for smaller systems, and that self-
gravitating disks, where our results break down (c.f. Franchini et al., 2021),
likely become relevant. Further modeling of the effects of retrograde
accretion on low-frequency gravitational wave backgrounds is warranted.

It is worth mentioning that varying the disk scale-height $h/r$ can
dramatically alter the standard circular prograde picture presented above, as
disks with $h/r\lesssim 0.05$ recover binary inspiral (Tiede et al., 2020;
Penzlin et al., 2022) and show tentative evidence that those with values
approaching the expected theoretical limit of $h/r\sim 0.01-0.001$ for
geometrically thin disks around massive black holes possess orders of
magnitude faster inspiral rates $\mathcal{A}\approx\mathcal{O}(10^{2})$
(Dittmann & Ryan, 2022). The dependence of eccentricity evolution for such
thin disks, however, is yet to be explored, and it remains unclear how
changing the scale-height in retrograde scenarios changes the solutions.

_Electromagnetic observations_ — While 2D isothermal hydrodynamics simulations
cannot self-consistently produce lightcurves of the accretion flow, the
accretion rate can serve as an approximation for the variability in the system
luminosity, and recent magnetohydrodynamics simulations have shown that the
variability of Poynting fluxes from component jets track the binary accretion
rate (e.g. Combi et al., 2022). One of the most prominent features of prograde
accretion solutions are the presence of a 5 orbit accretion variability
associated with an $m=1$ density mode that orbits the inner edge of the
prograde CBD for near-circular binaries. By contrast, the accretion rate for
circular retrograde binaries shows almost no variability, except for small
fluctuations at the orbital frequency. As documented in §4.2 and §4.5,
retrograde solutions also retain a high-degree of axisymmetry for $e<0.55$,
and as such are strongly modulated only on the orbital frequency of the
binary; and the periodic variability for eccentricities in this range are
almost indistinguishable, with the exception of a possible double peak
emerging for $e=0.4-0.5$, making them hard to distinguish in time-domain
observational data. At large eccentricities $e>0.55$, however, retrograde
binary accretion manifests a two-orbit double peaked accretion signature
characterized by a large flare, followed by a flare of approximately half the
magnitude. This kind of behavior has not been observed in any prograde
configurations and could serve as a unique observational signature of
retrograde binary accretion. At the highest eccentricities this “double
flaring” feature could manifest as quasi-periodic (or periodic if temporally
and flux resolved) eruptions on timescales of super-massive black hole binary
orbits (days to years).

Beyond variability in the accretion rate and system luminosity, the existence
of axisymmetric density waves propagating through the CBD for $0<e\lesssim
0.55$ may also have observational implications. The higher density rings, with
possibly higher temperature when considering non-isothermal equations of
state, would cause a time dependent variation of the disk spectra compared to
that of a steady disk, causing the opposite of a spectral notch proposed for
prograde disks (Gültekin & Miller, 2012; Roedig et al., 2014), but also
varying periodically in time. In addition, these retrograde ring waves may
manifest in images of circumbinary disks around stellar binaries, resolvable
with instruments such as the Atacama Large Millimeter Array (ALMA, e.g., Alves
et al., 2019). Synthetic observations of retrograde systems should explore
this possibility (e.g., Ragusa et al., 2021), as recent theoretical work
predicts $10\%$ of stellar circumbinary disk systems to form retrograde
(Elsender et al., 2023).

###  5.3 Numerical limitations

In order to densely sample the binary eccentricity and accurately calculate
the orbital evolution, we have made a number of simplifying assumptions.
Foremost, these simulations have employed an isothermal treatment of the gas
thermodynamics where any heat generated from viscosity or shocks is assumed to
instantly leave the system as radiation. In future work it would be
informative to include an equation of state that accounts for such heating
processes and subsequently cools the disk on an appropriate timescale. We have
also restricted ourselves to only two dimensional, co-planar configurations,
but possible out of plane effects have been previously reported (Nixon et al.,
2011b; Roedig & Sesana, 2014). We have additionally ignored magnetic fields
and considered only a simple constant-$\nu$ model for the disk viscosity. A
more detailed treatment would include electromagnetic effects and resolve the
magneto-rotational turbulence that would self-consistently govern angular
momentum transport in the disk.

The other major simplification of this study is that it focused exclusively on
equal mass-ratio binaries. Because of the complexity and importance of the
inner-most regions of the accretion flow to our results, we posit that varying
$q$ may produce substantially different results and plan to explore this in
future work.

##  6 Conclusions

We have conducted a number of simulations of equal mass binaries of varying
eccentricity accreting from retrograde circumbinary disks using the grid based
hydrodynamics codes Disco and Sailfish. We have continuously characterized the
effects of such retrograde accretion on the binary’s orbital elements,
documented the morphology and phase-dependent behavior of retrograde accretion
flows at numerous eccentricities, and determined observational signatures
associated with retrograde solutions.

In accordance with most theoretical expectation and previous analytic work, we
find that retrograde accretion scenarios shrink the binary semi-major axis at
all eccentricities with a near-constant rate $d\log a/d\log M=-10$, and
simultaneously pump eccentricity for any $0<e<0.8$ (Figure 1). The latter
point, however, is in contrast with previous SPH simulations and impulse-
approximation based models that predicted near-circular binaries would have
their eccentricity damped, remaining effectively circular. Moreover, we find
that the dominant contribution to the binary orbital evolution is
gravitational forces from the retrograde binary minidisks and retrograde-
bridge that form in the inner-most $r<a$ of the accretion flow (Figures 6 and
7). Asymmetries in these features yield gravitational forces that act at all
phases of the binary orbit, and the primary contributions to orbital
eccentricity pumping occur during binary waxing (Figures 4, 5, and 8).

We found that the morphology of the retrograde CBD’s posses three separate
regimes as we vary binary eccentricity (see Figure 3). First, near-circular
binaries $e<0.01$ display a predominantly time-invariant structure in the co-
rotating frame of the binary, with the exception of small wiggles in the
retrograde-bridge (Figure 2). In regime (_ii_), binaries with eccentricities
$0.025<e<0.55$ yield phase-dependent disk behavior that repeats every binary
orbit and is characterized by an axisymmetric density wave that is driven into
the disk once per orbit during binary waning (Figure 4). Lastly, binaries with
$e>0.55$ in regime (_iii_) manifest two-orbit periodic disk oscillations
characterized by the forcing of non-axisymmetric $m=2$ spiral density waves in
the “first” orbit and a predominantly axisymmetric density ring (with
comparatively weak spiral arms) in the “second” orbit (Figure 5). The
emergence of even-azimuthal-mode spiral density waves at $e>0.55$ is
consistent with previous predictions that disk resonances can manifest in
viscous, retrograde disks at large eccentricities (see also § 4.5 and Figure
9).

We have additionally analyzed the characteristic variability of the accretion
rates for binaries of many different eccentricities in retrograde CBD’s
(Figure 11). We observed that circular, retrograde solutions are dramatically
different from their prograde counterparts as they exhibit almost no accretion
variability due to their functionally steady-state configurations. Binary
eccentricities in retrograde regime (_ii_) all exhibit strong accretion rate
modulation at the orbital frequency with the possible emergence of a kinked or
double peaked structure at $e=0.4-0.5$. For eccentricities in regime (_iii_),
$e>0.55$, accretion rate modulations start varying at both the binary period
and twice the binary period as the accretion oscillates between minor- and
major-spikes.

In light of modern numerical studies that have begun characterizing the
behavior of prograde circumbinary accretion across system parameters like
eccentricity, mass ratio, and disk thickness, we have revisited the
complimentary retrograde scenario. Occurrences of retrograde accretion can be
an important contributor to binary orbital decay and eccentricity growth, and
they exhibit unique observational features for binary searches in large-scale
surveys. We have demonstrated that high resolution, full-domain treatments are
required to accurately quantify these effects and how they contrast with
prograde solutions. As such, retrograde investigations should continue
alongside their prograde counterparts as the community continues to
characterize these systems.

## Acknowledgements

The authors sincerely thank Jonathan Zrake and Paul Duffell for making their
codes Sailfish and Disco available for use in this work. C.T. and D.J.D extend
their gratitude to all members of the GW-Astro group at the NBIA and Jeff J.
Andrews for useful discussions. We also thank the anonymous referee for a
constructive report. D.J.D. received funding from the European Union’s Horizon
2020 research and innovation programme under Marie Sklodowska-Curie grant
agreement No. 101029157. D.J.D. and C.T. received support from the Danish
Independent Research Fund through Sapere Aude Starting Grant No. 121587.

## Data Availability

The data used for this study can be shared upon reasonable request to the
corresponding author. The codes used for the presented simulations are also
publicly available on Github at https://github.com/duffell/Disco and
https://github.com/clemson-cal/sailfish.

Figure A1: The effects of sink radius, sink rate, and gravitational softening
on the full orbital evolution of the binary (_left_) and the component of
evolution due to accretion alone (_right_) at fixed binary eccentricity
$e=0.3$. The pink curve shows the fiducial sink rate at varying sink radius
for both Disco (_dashed-dotted + circles_ ; also marked D in the supplementary
runs that vary the softening radius) and Sailfish (_solid + crosses_ ; SF in
the supplementary runs). The cyan curve shows the same experiment for a fast-
sink. The black diamond shows the use of a “standard” acceleration-free sink,
and the triangles illustrate the effects of varying the gravitational
softening radius. Note the change in normalization to $\dot{M}_{0}$ vs. the
measured $\dot{M}$. We observe that our results remain largely insensitive to
the sink radius and rate so long as the sink region is sufficiently resolved;
with the exception of large-fast sinks.

## References

  * Agazie et al. (2023a) Agazie G., et al., 2023a, arXiv e-prints, p. arXiv:2306.16220
  * Agazie et al. (2023b) Agazie G., et al., 2023b, ApJ, 951, L8
  * Alves et al. (2019) Alves F. O., Caselli P., Girart J. M., Segura-Cox D., Franco G. A. P., Schmiedeke A., Zhao B., 2019, Science, 366, 90
  * Amaro-Seoane et al. (2016) Amaro-Seoane P., Maureira-Fredes C., Dotti M., Colpi M., 2016, A&A, 591, A114
  * Antoniadis et al. (2023a) Antoniadis J., et al., 2023a, arXiv e-prints, p. arXiv:2306.16214
  * Antoniadis et al. (2023b) Antoniadis J., et al., 2023b, arXiv e-prints, p. arXiv:2306.16227
  * Bankert et al. (2015) Bankert J., Krolik J. H., Shi J., 2015, The Astrophysical Journal, 801, 114
  * Barnes & Hernquist (1996) Barnes J. E., Hernquist L., 1996, ApJ, 471, 115
  * Bate et al. (2010) Bate M. R., Lodato G., Pringle J. E., 2010, MNRAS, 401, 1505
  * Bortolas et al. (2021) Bortolas E., Franchini A., Bonetti M., Sesana A., 2021, ApJ, 918, L15
  * Combi et al. (2022) Combi L., Lopez Armengol F. G., Campanelli M., Noble S. C., Avara M., Krolik J. H., Bowen D., 2022, ApJ, 928, 187
  * Cuadra et al. (2009) Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, MNRAS, 393, 1423
  * D’Orazio & Duffell (2021) D’Orazio D. J., Duffell P. C., 2021, ApJ, 914, L21
  * D’Orazio et al. (2013) D’Orazio D. J., Haiman Z., MacFadyen A., 2013, MNRAS, 436, 2997
  * Dempsey et al. (2020) Dempsey A. M., Muñoz D., Lithwick Y., 2020, ApJ, 892, L29
  * Dittmann & Ryan (2021) Dittmann A. J., Ryan G., 2021, ApJ, 921, 71
  * Dittmann & Ryan (2022) Dittmann A. J., Ryan G., 2022, MNRAS, 513, 6158
  * Duffel et al. (2023) Duffel P., et al., 2023, In Prep
  * Duffell (2016) Duffell P. C., 2016, preprint,  (arXiv:1605.03577)
  * 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
  * Dunhill et al. (2015) Dunhill A. C., Cuadra J., Dougados C., 2015, MNRAS, 448, 3545
  * Elsender et al. (2023) Elsender D., Bate M. R., Lakeland B. S., Jensen E. L. N., Lubow S. H., 2023, MNRAS, 523, 4353
  * Farris et al. (2014) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ, 783, 134
  * Farris et al. (2015) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015, MNRAS, 447, L80
  * Franchini et al. (2021) Franchini A., Sesana A., Dotti M., 2021, MNRAS, 507, 1458
  * Franchini et al. (2022) Franchini A., Lupi A., Sesana A., 2022, ApJ, 929, L13
  * Garg et al. (2023) Garg M., Tiwari S., Derdzinski A., Baker J. G., Marsat S., Mayer L., 2023, MNRAS,
  * Goldreich & Tremaine (1979) Goldreich P., Tremaine S., 1979, ApJ, 233, 857
  * Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  * Gültekin & Miller (2012) Gültekin K., Miller J. M., 2012, ApJ, 761, 90
  * Haiman et al. (2009) Haiman Z., Kocsis B., Menou K., Lippai Z., Frei Z., 2009, Classical and Quantum Gravity, 26, 094032
  * Hayasaki et al. (2007) Hayasaki K., Mineshige S., Sudou H., 2007, Publications of the Astronomical Society of Japan, 59, 427
  * Heath & Nixon (2020) Heath R. M., Nixon C. J., 2020, A&A, 641, A64
  * Hobbs et al. (2011) Hobbs A., Nayakshin S., Power C., King A., 2011, MNRAS, 413, 2633
  * King & Pringle (2006) King A. R., Pringle J. E., 2006, MNRAS, 373, L90
  * Kley & Nelson (2012) Kley W., Nelson R. P., 2012, ARA&A, 50, 211
  * Lai & Muñoz (2022) Lai D., Muñoz D. J., 2022, arXiv e-prints, p. arXiv:2211.00028
  * Lodato et al. (2009) Lodato G., Nayakshin S., King A. R., Pringle J. E., 2009, MNRAS, 398, 1392
  * MacFadyen & Milosavljević (2008) MacFadyen A. I., Milosavljević M., 2008, ApJ, 672, 83
  * Mahesh et al. (2023) Mahesh S., McWilliams S. T., Pirog M., 2023, arXiv e-prints, p. arXiv:2305.01533
  * Major Krauth et al. (2023) Major Krauth L., Davelaar J., Haiman Z., Westernacher-Schneider J. R., Zrake J., MacFadyen A., 2023, arXiv e-prints, p. arXiv:2304.02575
  * Martini (2004) Martini P., 2004, Coevolution of Black Holes and Galaxies, p. 169
  * Mayer (2013) Mayer L., 2013, Classical and Quantum Gravity, 30, 244008
  * Milosavljević & Merritt (2003) Milosavljević M., Merritt D., 2003, ApJ, 596, 860
  * 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 et al. (2019) Muñoz D. J., Miranda R., Lai D., 2019, ApJ, 871, 84
  * Nixon & Lubow (2015) Nixon C., Lubow S. H., 2015, Monthly Notices of the Royal Astronomical Society, 448, 3472
  * Nixon et al. (2011a) Nixon C. J., Cossins P. J., King A. R., Pringle J. E., 2011a, MNRAS, 412, 1591
  * Nixon et al. (2011b) Nixon C. J., King A. R., Pringle J. E., 2011b, Monthly Notices of the Royal Astronomical Society: Letters, 417, L66
  * Orosz et al. (2012) Orosz J. A., Welsh W. F., et. al 2012, Science, 337, 1511
  * Penzlin et al. (2022) Penzlin A. B. T., Kley W., Audiffren H., Schäfer C. M., 2022, A&A, 660, A101
  * Pringle (1991) Pringle J. E., 1991, MNRAS, 248, 754
  * Ragusa et al. (2021) Ragusa E., et al., 2021, MNRAS, 507, 1157
  * Reardon et al. (2023) Reardon D. J., et al., 2023, ApJ, 951, L6
  * Roedig & Sesana (2014) Roedig C., Sesana A., 2014, MNRAS, 439, 3476
  * Roedig et al. (2014) Roedig C., Krolik J. H., Miller M. C., 2014, ApJ, 785, 115
  * Ryan & MacFadyen (2017) Ryan G., MacFadyen A., 2017, ApJ, 835, 199
  * Schnittman & Krolik (2015) Schnittman J. D., Krolik J. H., 2015, ApJ, 806, 88
  * 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. (2023a) Siwek M., Weinberger R., Muñoz D. J., Hernquist L., 2023a, MNRAS, 518, 5059
  * Siwek et al. (2023b) Siwek M., Weinberger R., Hernquist L., 2023b, MNRAS, 522, 2707
  * Sudarshan et al. (2022) Sudarshan P., Penzlin A. B. T., Ziampras A., Kley W., Nelson R. P., 2022, A&A, 664, A157
  * 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 et al. (2020) Tiede C., Zrake J., MacFadyen A., Haiman Z., 2020, ApJ, 900, 43
  * Tiede et al. (2022) Tiede C., Zrake J., MacFadyen A., Haiman Z., 2022, ApJ, 932, 24
  * Wang et al. (2022) Wang H.-Y., Bai X.-N., Lai D., Lin D. N. C., 2022, arXiv e-prints, p. arXiv:2212.07416
  * Wang et al. (2023) Wang H.-Y., Bai X.-N., Lai D., 2023, ApJ, 943, 175
  * Webbink (1976) Webbink R. F., 1976, ApJ, 209, 829
  * Westernacher-Schneider et al. (2022) Westernacher-Schneider J. R., Zrake J., MacFadyen A., Haiman Z., 2022, Phys. Rev. D, 106, 103010
  * Westernacher-Schneider et al. (2023) Westernacher-Schneider J. R., Zrake J., MacFadyen A., Haiman Z., 2023, arXiv e-prints, p. arXiv:2307.01154
  * Xu et al. (2023) Xu H., et al., 2023, Research in Astronomy and Astrophysics, 23, 075024
  * Zrake et al. (2021) Zrake J., Tiede C., MacFadyen A., Haiman Z., 2021, ApJ, 909, L13

$\begin{array}[]{ccc}\includegraphics[scale={0.35}]{Figures/rs0p025_soft0p025_ts1_Dens2D_CTspec_checkpoint_0019.h5}&&\\\
\includegraphics[scale={0.35}]{Figures/rs0p05_soft0p025_ts1_Dens2D_CTspec_checkpoint_0019.h5}&\includegraphics[scale={0.35}]{Figures/rs0p05_soft0p05_ts1_Dens2D_CTspec_checkpoint_0019.h5}&\includegraphics[scale={0.35}]{Figures/rs0p05_soft0p1_ts1_Dens2D_CTspec_checkpoint_0019.h5}\\\
&&\includegraphics[scale={0.35}]{Figures/rs0p1_soft0p1_ts1_Dens2D_CTspec_checkpoint_0019.h5}\\\
\end{array}$

Figure A2:  Density snapshots for different sink radii and gravitational
softening. The diagonal displays density snapshots corresponding to the sink
tests in Figure A1, with $t_{s}=1\Omega_{b}$ using the DISCO code (pink dash-
dotted lines in A1). From top left to bottom right, the diagonal changes sink
radius and softening together with values ($r_{s},r_{\rm soft}$) = ($0.025a$,
$0.025a$), ($0.05a$, $0.05a$), ($0.01a$, $0.1a$). The middle column fixes the
sink size to $0.05a$, but varies the softening from left to right as $r_{\rm
soft}=0.025a$, $0.05a$, $0.1a$.

## Appendix A: Dependence on Numerical Parameters

### Mass and Momentum Sinks

Unlike in prograde solutions, binaries embedded in a retrograde accretion disk
interact directly with the CBD every orbit. Because of this, as noted by
Amaro-Seoane et al. (2016), retrograde solutions can be sensitive to the sink
prescriptions that remove mass and momentum at each binary component location.
In order to examine the effect of varying the sink implementation on our
results, we performed an additional suite of simulations at fixed eccentricity
$e=0.3$ in both Disco (_dash-dotted_ lines, _circle_ markers) and Sailfish
(_solid_ lines, _cross_ markers) where we primarily varied the sink radius
$r_{s}$ at two different sink rates $t_{s}^{-1}=\gamma\Omega_{b}$ (see
definitions in §3). We chose to keep the gravitational softening the same size
as the sink radius, however below we relax this assumption as well. The
results of this study are shown in Appendix Figure A1 where the top-panels
show the measured change in binary semi-major axis $\dot{a}$, the bottom
panels show the measured change in eccentricity $\dot{e}$, the left column
shows values derived from both accretion effects and gravitational forces, and
the right column shows the effect due only to accreted mass and momentum
$\dot{a}_{\rm acc}$, $\dot{e}_{\rm acc}$ (note the change in normalization
from the measured binary accretion rate to the rate imposed at infinity
$\dot{M}_{0}$).

The pink lines show the effect of varying $r_{s}$ at the fiducial sink rate
$t_{s}^{-1}=1\,\Omega_{b}$., We find that at the fiducial sink-rate our
results are largely insensitive to changes in the sink radius, with $\sim
10\%$ change to orbital element rates of changes. At the fiducial sink rate,
there is a small decreasing trend in the magnitude of the total rates with
smaller sink radius. The effects due to accreted mass and momentum
$\dot{a}_{\rm acc}$, $\dot{e}_{\rm acc}$ do not vary significantly with
changing $r_{s}$.

The cyan lines illustrate the same experiment in the fast-sink regime with
$t_{s}^{-1}=10\Omega_{b}$. For small enough sink radius $r_{s}\leq 0.05\,a$,
these results for the full $\dot{a}$ and $\dot{e}$ appear converged to the
fiducial system, and only excessively large sinks $r_{s}\sim 0.1\,a$ begin to
show pathological behavior in the disk morphology and measured orbital
evolution. This can also be seen in the emergence of an outsized contribution
to the orbital evolution from the sink terms in the right column of Figure A1.
We do note, however, that in Disco, while the full results are largely
unaffected by the fast sink, the disks with fast sinks developed asymmetric
and kinked retrograde-bridges, likely leading to the sign change in
$\dot{e}_{\rm acc}$. Disco runs with fast sinks and $r_{s}=0.1\,a$ also
developed numerical instabilities that prevented the completion of those runs.

For completeness, we also include one Sailfish run with a “standard”
acceleration-free sink where $\mathbf{S_{J}}=S_{\Sigma}\mathbf{v}$ marked with
a black diamond. Changing from spinless to acceleration-free sinks similarly
does not meaningfully alter our results, and only serves to change the sign of
the eccentricity evolution component from the accretion terms alone.

### Gravitational Softening

Because the changes in $\dot{a}_{\rm acc}$, $\dot{e}_{\rm acc}$ are much
smaller than the changes to the total rates (except for the fastest and
largest sinks), the changes to the total rates must be primarily from the
changing gravitational softening length. To test this, we run a smaller and a
larger softening length for the fiducial sink, sampling $r_{\rm
soft}=r_{s}/2$, $r_{s}$, $2r_{s}$. These are denoted by the triangles in
Figure A1. The primary result, as expected, is that $\dot{a}_{\rm acc}$,
$\dot{e}_{\rm acc}$ are very weakly, or not affected, while the total rates
(left panel) vary on a level similar to that from changing the sink radius
above. This is due to the gravitational softening affecting the density
structure in the vicinity of the minidisks, and thus, affecting the primary
source of gravitational force on the binary (see §4.3).

Figure A2 shows log-density snapshots for the $e=0.3$ binary sink and
softening test, from the Disco runs at fiducial sink rate. From left to right,
the gravitational softening is increased from small to large. From top to
bottom the sink radius is increased from small to large. Hence, the diagonal
from top left to bottom right shows the results of the sink test denoted by
the dot-dashed pink lines in Figure A1, where softening and sink radius are
increased together. The middle row shows the effect of changing the softening
while fixing the sink properties (also plotted in Figure A1, see legend). It
is clear from the density snapshots that the wake structure is affected by
softening, thus resulting in different gravitational forces which are the
primary reason for differing values in Figure A1. While the general behavior
remains unaltered, the magnitude of semi-major axis shrinking and eccentricity
driving may be sensitive to the chosen gravitational softening. We note that
the strongest effect appears in $\dot{a}$ from Sailfish, but the increased
non-linear features from a smaller softening radius destabilized this run, and
it would require higher resolution to characterize accurately.

Note that we do not expect reliable results for smaller softening radii as the
fiducial resolution employed (in the Disco runs) covers the diameter of the
core of the softened potential by only four cells, or by approximately 16
cells in area. Higher resolution is needed to probe the effects of gas
interacting with stronger regions of the point-mass potential. We leave a more
detailed exploration of this effect to future work.