Eccentric Binaries in Retrograde Disks

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 [46] (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; 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 spiral 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 (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 $\sim5\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 $\sim5$ 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 \lesssim0.1$ have their eccentricity damped towards circular orbits—which they are driven apart by the gas—but all other initial binary eccentricities $e \gtrsim0.1$ are driven towards an equilibrium eccentricity $e \sim0.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} \lesssim0.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 isotropically, 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., $\sim10\%$ 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 material³. 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.

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, \ldots$ as $\Omega(r) = m\omega_b/(m \pm1)$. 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 enhance 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 $\propto1+q$ for when $q \ll1$ 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 \gtrsim0.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 § 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 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 \to1$ 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 a 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 \sim0.6$ where prograde 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).

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 Navier-Stokes equations for viscous, locally isothermal hydrodynamics in 2-dimensions. We refer the reader to an upcoming code-comparison paper Duffell et al. (2023) for implementation-specific details as well as a direct comparison of prograde CBD solutions.

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/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$:

As both sets of initial conditions capture steady state disk solutions far from the binary, and because we allow the disk to relax for $\gtrsim500$ 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 $\tilde{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

and $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 disk's self-gravity (e.g., the disk Toomre parameter $Q \sim\mathcal{M}^{-1}(M/M_d) \gg1$).

Source terms and boundary conditions

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

where $v_i$ is the mid-plane fluid velocity, $P = c_s^2\Sigma$ is the vertically averaged gas pressure, $\boldsymbol{\tau}$ is the viscous stress tensor, $\mathbf{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

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

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$:

Different sink behaviors are achieved through calculation of the adjusted gas-velocity $\mathbf{v}_i^*$ 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

where $\mathbf{v}$ is the gas velocity in the inertial frame, $\mathbf{v}_i$ is the binary component velocity, and $\hat{\mathbf{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^* = \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).

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

We take $n = 5\times10^3$ 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.5a$, slightly higher inside $r < 0.5a$, 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 \gg$ few $\times\,a$, with possible exception at large eccentricity to be discussed later.

Numerical Results

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

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 \approx3$ at $e \approx0.15$ and settle to a near constant growth rate $de/d\log M \approx2.25$ by $e \gtrsim0.4$. The shrinking of the binary $a$ shows almost no dependence on binary eccentricity and always retains a value $d\log a/d\log M \approx-10$.⁴

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 \sim0.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).

Disk morphology

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 orbiting 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 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 § 6.

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 \lesssim0.55$, and (iii) the forcing of non-axisymmetric spiral density waves at large eccentricities $e \gtrsim0.55$. Circular binaries (i) with $e \lesssim0.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 \lesssim0.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 \approx0.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 \approx0.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 Figure 10 and Figure 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.

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 $\dot{a}_{\mathrm{grav}}$, $\dot{e}_{\mathrm{grav}}$ (illustrated as light-purple and light-orange), and the components due 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}_{\mathrm{grav}}$, $\dot{e}_{\mathrm{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 $\sim10\%$ 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 $\lvert(\mathbf{r}_{b} \times\mathbf{f}_{g})/\sqrt{GMa(1-e^{2})}\rvert$ and power $\lvert(\mathbf{f}_{g} \cdot\mathbf{v}_{b})/(GM/2a)\rvert$, with $\mathbf{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 \lesssim0.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.

Eccentricity driving

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—or 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

such that $\dot{e}_{\mathrm{grav}}=\dot{e}_{\mathrm{power}}+\dot{e}_{\mathrm{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 $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}_{\mathrm{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 Figure 4 & Figure 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 \approx0.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 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$.

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 \ge0.6$. The simulations presented in this study also suggest evidence for the emergence of retrograde resonances at eccentricities $e \gtrsim0.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 CBD.⁵ To check, albeit indirectly, for the presence of such resonances, we calculate timeseries of the azimuthal density mode $m = \{1, 2, 3, 4\}$ as

Figure 9 shows the time averaged strength of these azimuthal density modes $\overline{\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 \ge0.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.

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

has been suggested as evidence for non-resonant effects in the cavity-carving process for equal mass binaries (Mahesh et al. 2023).

Fourier vector with angular frequency $\omega$ in frequency space,

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 \times300$ 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_{\mathrm{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 \approx0.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 \sim0.2$. From $e \sim0.2 - 0.5$, the amplitude of accretion rate modulations saturates to $\sim50\%$ of the mean, and the time of peak accretion rate drifts from just after apocenter towards just before pericenter with increasing $e$. At $e \gtrsim0.3$, a second peak in the accretion rate develops, occurring just after apocenter and before the other peak, which, for $0.3 \lesssim e \lesssim0.5$, occurs just before pericenter. For $e \gtrsim0.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 $\sim2$ 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 \gtrsim0.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 $\sim0.7$ compared to $e \sim0.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).

Discussion

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 \sim0.15$ after which it stabilizes to $\mathrm{d}e/\mathrm{d}\log M \sim2$.

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-orbit 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)).

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 \gtrsim0.05$; however, prograde binaries with initial eccentricity $e \gtrsim0.08$ evolve towards an equilibrium eccentricity of $e_{\rm eq} \sim0.4$ (or some value $0.25 \lesssim e_{\rm eq} \lesssim0.5$ for $q < 1$ Siwek et al. (2023b)) where they shrink their semi-major axis at a rate $\mathcal{A} \equiv\mathrm{d}\log a/\mathrm{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,

For the remaining discussion we choose $\mathcal{E}=30$, $e_*=0.1$, and $\mathcal{A}=-10$. Note, however, that a range of $e_* \sim0.1$–$0.2$ and $\mathcal{E}e_* \sim2.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} \sim4.5\ \mathrm{Myr}$ where we’ve used $M/\dot{M} \approx4.5 \times10^7\ \mathrm{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 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 \lesssim0.1$, with an estimated form $30e\ \dot{M}/M$, and approximately constant for $e \gtrsim0.1$. Therefore, for initial eccentricities $e_0 \lesssim0.1$, the binary will grow its eccentricity exponentially, with an e-folding timescale of $(\sim30\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 \to1$ in approximately half of a Salpeter time, $\sim20$ Myr. Comparing the mutual evolution of $\dot{a}$ and $\dot{e}$, this implies that $a$ will decrease by 5 e-foldings $a_0e^{-5}$ in the time required for $e \to1$.

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

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_0e^{-5}$. This ignores the commensurate shrinking of semi-major axis during eccentricity pumping.⁶

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\gtrsim0.9$ into massive binaries that are eventually detected by LISA, if accretion commences at large enough separations $a_0 \sim10^{-2}\,\mathrm{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_{\mathrm{eff}} < 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 \lesssim0.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 \sim0.01-0.001$ for geometrically thin disks around massive black holes possess orders of magnitude faster inspiral rates $\mathcal{A} \sim 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 \lesssim0.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).

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.

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 (Figure 6 and Figure 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 (Figure 4, Figure 5, and Figure 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 Skłodowska-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.

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 $\sim10\%$ 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 \leq0.05a$, these results for the full $\dot{a}$ and $\dot{e}$ appear converged to the fiducial system, and only excessively large sinks $r_s \sim0.1a$ 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.1a$ 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 $S_J = S_{\Sigma}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 ran 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.

This paper has been typeset from a TEX/LATEX file prepared by the author.