The Evolution of Accreting Binaries: from Brown Dwarfs to Supermassive Black Holes

Circumbinary accretion occurs throughout the universe, from the formation of stars and planets to the aftermath of major galactic mergers. We present an extensive investigation of circumbinary accretion disks, studying circular binaries with mass ratios ($q\equiv M_2/M_1$) from 0.01 to 1 and at each mass ratio probing the effects of disk thickness and viscosity. We study disks with aspect ratios $H/r\in\{0.1, 0.05, 0.033\}$, and vary both the magnitude and spatial dependance of viscosity. Although thin accretion disks have previously been found to promote rapid inspirals of equal-mass binaries, we find that gravitational torques become weaker at lower mass ratios and most binaries with $0.01\leq q\leq0.04$ outspiral, which may delay the coalescence of black hole binaries formed from minor mergers and cause high-mass exoplanets to migrate outwards. However, in a number of cases, the disks accreting onto binaries with mass ratios $\sim 0.07$ fail to develop eccentric modes, leading to extremely rapid inspirals. Variability in black hole accretion correlates with disk eccentricity, and we observe variability above the $\sim10\%$ level even for mass ratios of $0.01$. We demonstrate that the spatial dependence of the viscosity (e.g. $\alpha$ vs constant-$\nu$) significantly affects the degree of preferential accretion onto the secondary, resolving discrepancies between previous studies. Colder circumbinary disks remain eccentric even at $q\sim0.01$ and sustain deep, asymmetric cavities.


INTRODUCTION
Circumbinary accretion disks, in a general sense, occur whenever a binary of sufficient mass interacts with a surrounding gaseous disk.When the binary members are of comparable mass, such systems arise, for example, during star formation (e.g., Boss 1986;Bate & Bonnell 1997;Park et al. 2022), where interactions with between the disk and binary may shape the architecture of observed stellar systems (El-Badry et al. 2019;Hwang et al. 2022); and after the merger of comparable-mass galaxies, where interactions with gaseous disks may expedite supermassive black hole coalescence (e.g., Begelman et al. 1980;Tiede et al. 2020;Dittmann & Ryan 2022).Circumbinary accretion systems containing binaries of more disparate masses can result from, for example, the forma-tion of high-mass planets in protoplanetary disks (e.g., Chauvin et al. 2004;Bryan et al. 2018), and the accretion of gas onto the binaries resulting from the merger of massive and dwarf galaxies.
Binary black holes are likely to form in the centers of galaxies following mergers, and gravitational interactions between a binary and the surrounding gas and stars may help these black holes reach sufficiently small separation that the binary can be driven to merge by gravitational radiation (e.g., Begelman et al. 1980;Haiman et al. 2009;Khan et al. 2013).While major mergers typically result in the formation of binaries with comparable masses, the sheer abundance of dwarf galaxies across cosmic time suggests that mergers between galaxies of greatly disparate masses should be much more common (e.g., Haehnelt & Kauffmann 2000;Reed et al. 2007).Although the orbital evolution of binaries formed from mergers involving dwarf galaxies may be complicated by factors such as the dark matter distributions in the merging galaxies (Tamfal et al. 2018), tidal stripping (Taffoni et al. 2003), and threebody dynamics (Amaro-Seoane et al. 2010), the binaries formed from such galaxy mergers may dominate the population of black hole mergers detected by LISA (e.g., Bellovary et al. 2019;Amaro-Seoane et al. 2022, which suggest a peak at q ∼ 0.05).Even if gas does not play a dynamically significant role in the orbital evolution of these binary systems, binary-driven variability may be detectable in these systems, as well as presently studied binary AGN candidates (e.g.Graham et al. (2015); Charisi et al. (2016); Chen et al. (2020), but see Liu et al. (2016); Vaughan et al. (2016); Liu et al. (2018) for cautionary remarks).Additionally, appreciable alterations of binary black hole masses and orbital parameters can appreicably modify the spectrum of the stochastic gravitational wave background (e.g., Rajagopal & Romani 1995;Phinney 2001;Siwek et al. 2020;Franchini et al. 2021;Agazie et al. 2023).
On more local scales, as the population of detected exoplanets with high mass ratios grows we can begin to glean insights into different populations thereof.For example, although gravitational instabilities might produce large-mass large-separation directly imaged exoplanets (e.g., Bowler & Nielsen 2018), their mass distribution may also be consistent with the core accretion scenario (Nielsen et al. 2019;Wagner et al. 2019).Additionally, the factors that limit accretion onto high-mass planets are also uncertain but of great interest (e.g., Rosenthal et al. 2020;Li et al. 2021), as these affect the upper limits on planetary masses.
The problem of accretion disks around unequal-mass binaries has been investigated numerous times (e.g., Artymowicz & Lubow 1996;Bate & Bonnell 1997;Hayasaki et al. 2007;Shi & Krolik 2015;D'Orazio et al. 2016;Muñoz et al. 2020;Ragusa et al. 2020;Penzlin et al. 2022;Mahesh et al. 2023;Cimerman & Rafikov 2023).Although recent studies of thin disks have found nearly universally that the lower-mass binary member accretes preferentially (e.g., Farris et al. 2014;Duffell et al. 2020;Dittmann & Ryan 2021;Siwek et al. 2023), the precise degree of preferential accretion (measured as the ratio of the accretion rate onto the secondary to that onto the primary) at a given mass ratio has been found to vary by factors of ∼ 3 even between studies utilizing different versions of the same code (e.g., Duffell et al. 2020;Dittmann & Ryan 2021), due at least in part to different treatments of accretion and viscosity.Additionally, studies of accretion onto binaries from thick disks have observed preferential accretion instead onto the primary, behavior which has been attributed to the effects of gas pressure on the dynamics of gas within the Hill sphere of the secondary (c.f.Ochi et al. 2005;Hanawa et al. 2010;Young et al. 2015).
Herein we present a systematic survey of accretion onto binaries of unequal masses.We study binaries with mass ratios q ∈ [0.01, 0.02, 0.03, 0.04, 0.05, 0.07, 0.1, 0.14, 0.2, 0.3, 0.5, 1].For each mass ratio, we simulate disks with viscosities ν = 0.001Ω b a 2 b and M = {10, 20, 30}, where Ω b is the orbital angular frequency of the binary, a b its semi-major axis, and M is the characteristic azimuthal Mach number of the disk.At M = 10 we additionally simulate disks with viscosities of ν = 0.0005Ω b a 2 b and α−viscosities, setting α = 0.1.We describe our numerical methodology and simulation diagnostics in Section 2. We present the main results of our study in Section 3, discuss their relevance to various astrophysical scenarios in Section 4, and summarize our results in Section 5. We study the convergence of our simulations in Appendix A, demonstrating not only that our results are generally converged with respect to both numerical resolution and mass removal rate, but also the pernicious effects of overly simple sink prescriptions.Appendix B makes direct comparisons between our work and two prior investigations, Duffell et al. (2020) and Dempsey et al. (2021), in the low-mass ratio regime.Appendix C collects the quantitative results of our main suite of simulations.

METHODS
This study primarily utilizes Athena++ (Stone et al. 2020),1 an Eulerian code, to solve the compressible Navier-Stokes equations, given by (1) Σ is the disk surface density, Π is the verticallyintegrated pressure, v is the fluid velocity vector, Φ is the gravitational potential, ↔ I is the identity matrix, S Σ is a mass sink term, S p is a momentum sink term, ↔ σ is the velocity shear tensor, ν is the kinematic viscosity, and additional damping source terms are given by D Σ and D p , the details of which are presented in Section 2.2.
We implement a 'locally isothermal' equation of state similarly to Moody et al. (2019) by setting the adiabatic index of the fluid to Γ = 1.0001 and implementing a source term which instantaneously sets Π in each cell consitently with a prescribed local sound speed c s : The squared sound speed is given by where M is a fiducial Mach number which we hold fixed in each simulation.We use a softened gravitational potential, where the total potential is given by Φ = i Φ i and Φ i is the gravitational potential of a point mass of mass M i at position x i given by where ϵ g,i is a gravitational softening length.We hold the binary orbit fixed, at a semi-major axis a b and mass M = M 1 +M 2 , over the course of each simulation, calculating positions analytically each time step, which is justified as long as evolutionary timescales are many times longer than the binary orbital period and simulation duration.
We consider multiple kinematic viscosity (ν) prescriptions, including a globally constant kinematic viscosity and an α-viscosity ansatz (Shakura & Sunyaev 1973) such that ν = αc s H, where H is the disk scale height which we calculate using which is appropriate to a thin disk in vertical hydrostatic equilibrium about multiple point masses with the gravitational potentials given by Equation (5).Because the degree of preferential accretion onto the secondary depends on the functional form of the viscosity, we test an additional viscosity prescription with an even stronger distance dependence than that of an α-viscosity in Section 3.3.1.The surface density sink term is given by where s i is a function specifying the sink profile for each particle, Ω b is the angular frequency of the binary, and γ i is the sink rate for each particle.Throughout this work, we use the sink profile where the sink radius for each particle r s,i is a purely numerical parameter defining a characteristic radius over which the sink term acts.
We use a torque-controlled momentum sink term, given by where v i is the velocity of the sink particle, ri and φi are the unit basis vectors at position x of a polar coordinate system centred on the sink particle, and δ is a dimensionless control parameter.Throughout most of this work, we set δ = 0, which results in a sink term that is torque-free in the frame of a given sink particle, leading to no change in the spin angular momentum of this sinks.Such sink methods make simulation results more robust to numerical parameters such as γ (Dittmann & Ryan 2021), and are able to reproduce analytic steadystate disk profiles in single-object disks (Dempsey et al. 2020;Dittmann & Ryan 2021).However, in Appendix A, we present the results of a few simulations which used δ = 1, which sets v * = v, to illustrate how this choice intrinsically leads to spurious results that depend strongly on the sink rate (γ).

Diagnostics
Over the course of each simulation, we record the matter and angular momentum accreted by each binary component, as well as the angular momentum exchanged between each member of the binary and the surrounding gas through gravitational interactions.We integrate these quantities on every time step over the entire simulation, and output these volume-and time-integrals at a cadence of 100 outputs each binary orbital period.
The accretion rate onto each particle is given by where We define the momentum sink term for each particle analogously as S p,i = S Σi v * i via Equation ( 9).Although these integrals formally extend over the entire domain, the sink profile function s i (x) sharply truncates contributions from beyond more than a few r s from a given sink particle.The torque on each particle due to accretion Ja is then given by The change in angular momentum of each particle is computed every time step using The total torque on the system of sink particles is the sum of the torques on each component: which is also the rate of change of the binary orbital angular momentum because we employ (spin) torquefree sink particles.We calculate the orbital evolution of the binary from the torques measured in our simulations using the equation for the binary orbital angular momentum where M is the total mass of the binary and we have specialized to circular binaries.Previous studies have found that small eccentricities e b ≲ 0.08 are damped through interaction with circumbinary disks (Muñoz et al. 2019;Zrake et al. 2021), so this approximation is appropriate for binaries which form with low initial eccentricities.
The torque on the binary can be expressed as from which an expression for the evolution of the binary semi-major axis follows: where q = M 2 /M 1 is the mass ratio, λ ≡ Ṁ2 / Ṁ1 is the degree of preferential accretion onto the secondary, l 0 = Jorb / Ṁ is the orbital angular momentum change in the binary per unit mass, and l b is the specific angular momentum of the binary Ω b a 2 b q/(1 + q) 2 .If gas is able to circularize before accreting, then the contribution to the overall torque due to accretion is given by Ja = Ṁ1 l 1 + Ṁ2 l 2 , where l 1 = q 2 /(1 + q) 2 and l 2 = 1/(1 + q) 2 are the specific angular momenta of the binary components.When this is the case (only when, in addition to the accretion flow becoming circular near the accreting masses, the entire angular momentum accreted contributes to the orbital angular momentum of the system rather than the spin), Equation ( 16) reduces to Notably, if the aforementioned assumptions hold, the evolution of the binary semi-major axis only depends on the gravitational torque, overall accretion rate, and the specific angular momentum of the binary.The evolution of the mass ratio through accretion can be calculated similarly, using (20) We have also computed a number of one-dimensional radial histograms: binning in radius from both the binary barycenter and in coordinate systems centered on each objects, and also binning quantities as functions of fluid element semi-major axis about the barycenter.The semi-major axis of each fluid element is calculated from its radial location and the magnitude of its velocity (v) as a = (2GM/r − v 2 ) −1 , and the eccentricity vector of each fluid element is calculated as e = (v 2 x − (x • v)v)/GM − x/r, and where r ≡ |x| is the distance from the binary barycenter.Every time step, we include time-and either volume-or mass-weighted samples of fluid surface density, eccentricity, and gravitational torque density.We use the lattermost to calculate radially integrated gravitational torque profiles We used these histograms to measure various quantities describing the disk.For example, we use histograms of the surface density and mass-weighted eccentricity, as functions of binary semi-major axis, to determine quantities such as the characteristic cavity semi-major axis (a cav ) and eccentricity (e cav ).We define the cavity location as that of the minimum of the time-averaged |∂ 2 Σ(a)/∂a 2 | for a b < a < 6a b , which is valid for all of our simulations, as the surface density profile Σ(a) always has an inflection point at a > a b .We define the cavity eccentricity as the average eccentricity of the disk at a cav .We derived radial profiles of the disk argument of periapsis over time using the mass-weighted histograms of e, defining the argument of periapsis as ω = tan −1 e y /e x .We measured the average disk precession rate by measuring the time derivative of ω between a = a cav and a = a cav + a b .

Initial Conditions and Discretisation
The domain of our simulations was centered on the binary barycenter.We initialized each simulation with an axisymmetric surface density distribution of the form (21) where r c sets the initial radius of a cavity around the binary, ξ sets the steepness of this cavity and η is a constant which sets the minimum density in our initial condition, p sets the surface density profile far from the binary and is given by p = 0 for simulations using a constant kinematic viscosity and p = −1/2 for simulations using an α-viscosity prescription, ℓ 0 is an initial guess for the eventual values of l 0 in each simulation, and Σ 0 is a fiducial surface density.Setting ℓ 0 near l 0 has been shown to significantly expedite the process of simulations settling into steady states of globally constant (in a time-averaged sense) accretion rate (Miranda et al. 2017) and prevent temporary suppression of the accretion rate onto the binary at high Mach numbers (Dittmann & Ryan 2022).We set r c = 2.5a b , ξ = 10, and η = 0.0001 in all of our simulations.The initial radial velocity in our simulations was given by where R = r 2 + ϵ 2 g .The initial azimuthal velocity (v ϕ = rΩ) is specified through approximate hydrostatic equilibrium taking into account both pressure gradients and the binary quadrupole moment, Our simulations are carried out in Cartesian coordinates, their domains extending from −10a b to 10a b in both x and y.Although we used outflow boundary conditions, we also set the source terms D Σ and D p to damp fluid variable towards their initial values near the outer boundary.For each conserved variable U i , the corresponding source term D i takes the form where and we set β = 10, r d = 9.5a b , and l d = 0.5a b .These damping terms serve to suppress m = 4 instabilities related to our use of a Cartesian grid, and have been employed previously in Dittmann et al. (2023).
Our simulations at high mass ratios (q ≥ 0.1) utilized static mesh refinement, refining the grid by a factor of two within [−5a b , 5a b ], [−2.5a b , 2.5a b ], and [−1.25a b , 1.25a b ].We typically resolve our base grid using 384 cells along each dimension, reaching a resolution of ∆x = ∆y ≈ 0.0065a b at the highest refinement level.These simulations held ϵ g,i = 0.035a b and r s,i = 0.035a b .For our ν = 0.001Ω b a 2 b simulation we set γ i = 2.67, and for our ν = 0.0005Ω b a 2 b simulations and those using an α−viscosity prescription we set γ i = 1.33.
For our simulations at low mass ratios (q < 0.1), we utilize the adaptive mesh refinement framework (although in a static manner) to keep the largest cell size within each mesh block as close to ∆x ≈ max[(r−r 2 ), r c ]ϵ as possible, where ϵ is a target cell aspect ratio, r c is a cutoff length, and (r − r 2 ) is the difference between the radial location of a given cell and the radius of the orbit of the secondary.This choice of refinement focuses resolution near the orbit of the secondary while not requiring very small cells near the much more massive primary, allowing us to resolve the flow of gas within the Hill sphere of the secondary (of radius r H = a b (q/3) 1/3 ) while reducing the time step penalties incurred by uniformly resolving the inner regions of the accretion flow.We typically set ϵ = 4/300 and r c = 0.15a b over a grid with a base resolution of ∆x ≈, resulting in cell sizes of ∆x ≈ 0.0016a b near the secondary and ∆x ≈ 0.0065a b near the primary.
For our low-mass ratio runs, we scale the softening and sink radii of each point mass according to ϵ g,i /ϵ g,0 = r s,i /r s,0 = g(m i ), where g(m i ) = (2m i ) 1/3 so that the sink and softening lengths are scaled according to the tidal radius about each particle, and set r s,0 = ϵ g,0 = 0.035a b .We scale the sink rate according to γ i = γ 0 g(m i ) 2 in our constant-ν simulations, and according to γ i = γ 0 g(m i ) 3/2 for constant-α simulations to keep a roughly constant ratio of the viscous timescale at the sink radius ∼ r s,i 2/ν to the accretion timescale (∼ Ω −1 b /γ).We set γ i = 2.67 in our ν = 0.001Ω b a 2 b simulations, and for our ν = 0.0005Ω b a 2 b simulations and those using an α−viscosity prescription we set γ i = 1.33.
All simulations used 3rd−order piecewise-parabolic reconstruction (Felker & Stone 2018), the second-order van Leer integrator (Stone & Gardiner 2009) with a CFL factor of 0.15, and the Roe approximate Riemann solver (Roe 1981).We have found that, particularly at low mass ratios, high-order spatial reconstruction schemes significantly aid convergence.We chose initial values of ℓ 0 for each Mach number based on the simulations presented in Dittmann & Ryan (2022) and the observation that the angular momentum current through the disk varies quite weakly with binary mass ratio, at least at M = 10 (e.g., Muñoz et al. 2020;Dittmann & Ryan 2021), using ℓ 0 = {0.7,0.1, −0.5} for M = {10, 20, 30} respectively.

RESULTS
To begin a sweeping overview of our simulations, we present snapshots of the gas surface densities from a subset thereof in Figure 1; a number of trends are imme-Figure 1. Surface density profiles from a subset of our simulations, increasing in binary mass ratio from left to right.Profiles from simulations with q < 0.1 are presented at t = 1000π, and profiles from simulations with q ≥ 0.1 are presented at t = 4000π.Snapshots from simulations with q ≥ 0.1.Results from simulations of M = 10 disks with α = 0.1 and ν = 0.0005a2 b Ω b are presented in the first and second rows respectively.The third, fourth, and fifth rows display results from simulations of ν = 0.001a 2 b Ω b disks at M = {10, 20, 30} respectively.In all cases deeper and larger cavities form at higher mass ratios.At higher Mach numbers the cavity remains visibly eccentric even at very low mass ratios, although in thicker disks the cavity eccentricities driven by lower mass ratio binaries become more subtle.Tick markings are spaced 2a b apart.
diately apparent.Across all disk viscosities and Mach numbers, higher mass ratio binaries open deeper and more pronounced cavities, 2 which follow directly from the weaker shocks driven by commensurately weaker torques from lower mass ratio binaries.Higher mass ratio binaries often drive more eccentric cavities, but this is not always the case: for example, out of our M = 20 simulations, the disk orbiting the q = 0.03 binary is more eccentric than that around the q = 0.05 binary.Previous studies of equal-mass binaries have noted that inner edges of higher Mach number circumbinary disks have higher surface densities than those around lower Mach number disks (e.g., Tiede et al. 2020), which relates to the angular momentum current through those disks as demonstrated in Dittmann & Ryan (2022); we observe here the same trend at lower mass ratios, at least for q ≳ 0.1.Additionally, as noted by previous studies of protoplanetary and circumbinary disks, lower viscosity disks are characterized by deeper and wider cavities at a Values of Ṁ2/ Ṁ1, Jg/ Ṁ , and J/ Ṁ from q = 0.01 to q = 1 at M ∈ {10, 20, 30} disks with ν = 0.001, averaged in time over the last 500 orbits of each simulation.Although all quantities are non-monotonic, as soon as symmetry is broken and M2 < M1, Ṁ2 > Ṁ1.Although, when moving from high to low mass ratios, gravitational torques initially grow weaker, the gravitational torque suddenly (in q) becomes negative, at mass ratios which depend on Mach number.These negative torques correspond to decreases in Ṁ2/ Ṁ1 in the M = 20 and M = 30 cases, which are discussed in Section 3.2.1.Numerical values are provided in Table 4.
We provide visual summaries of the quantities related to binary orbital evolution Only a few trends appear to be universal.The total torque on binaries in M = 10 disks is always positive (although in some cases the gravitational torque on the binary is negative).Additionally, unequal-mass binaries always accrete such that their mass ratio increases over time, and that disks with α viscosities accrete more pref- .The evolution of binary mass ratios and semimajor axes in the same simulations shown in Figure 2. At larger mass ratios, we find results in general agreement with previous studies.However, at intermediate mass ratios, which have not been previously probed at M ̸ = 10, we find extremely rapid inspirals, following from both the appreciably negative torque, small specific angular momentum of the binary, and that λ > q.At very low mass ratios (q ≲ 0.03), the rate of change of the binary semi-major axis can be quite large unless the torque is very nearly balanced by the change in the mass of the binary.Numerical values are provided in Table 4.
erentially onto the secondary than those with globally constant kinematic viscositites.
The semi-major axis and eccentricity of the inner edge of the circumbinary disk tend to increase as the binary mass ratio increases, although not monotonically.At higher mass ratios the circumbinary disk always precesses in a prograde sense, and a rate which decreases (non-monotonically) at lower mass ratios.For the M = 10 disks considered here, the precession rate becomes negative, or retrograde with respect to the binary, at many (but not all) mass ratios below q ≲ 0.05.
At mass ratios between q ∼ 0.1 and q ∼ 1, a few other trends trends emerge across disk parameters.For example, material tends to accrete more and more preferentially onto the secondary as mass ratios decrease from q = 1, Ṁ2 / Ṁ1 peaking, at least locally in q, between q = 0.1 and q = 0.2.The gravitational torque on the binary also becomes gradually weaker, lessening in magnitude -this can be qualitatively understood as the accretion flow near the more massive primary approaching axisymmetry at smaller mass ratios, while the secondary shrinks in mass and gravitational influence.We examine these trends more quantitatively in the following subsections, before moving on to features more localized in mass ratio such as the sudden decreases in the gravitational torque at in M = 30 disks at mass ratios q = 0.05 and q = 0.07.We find appreciable changes in Values of Ṁ2/ Ṁ1, Jg/ Ṁ , and J/ Ṁ from q = 0.01 to q = 1 with M = 10 and disks with viscosities ν = 0.001, ν = 0.0005, and α = 0.1, averaged in time over the last 500 orbits of each simulation.At all mass ratios below 1, we find that binaries accreting from α−viscous disks accrete more preferentially onto the secondary than those accreting from constant−ν disks.Depending on the viscosity prescription, binaries at intermediate mass ratios, 0.05 ≲ q ≲ 0.2 experience negative gravitational torques, which we discuss further in Section 3.2.2.Numerical values are provided in Table 4.
neither binary orbital evolution nor disk morphology related to the stability of particle orbits about the fourth and firth Lagrange points of the circular restricted three body problem, contrary to the advocation of D'Orazio et al. ( 2016) -evidently viscous and hydrodynamical effects are more significant in this regime.
3.1.Accretion morphology around q ≳ 0.1 binaries We begin by reviewing simulations which display characteristics most similar to those of equal-mass binaries, which have been studied in great detail previously.Higher mass ratio binaries open pronounced eccentric cavities in their circumbinary disks.Whenever one of the binary members approaches the pericenter of the cavity, some material is captured by the binary while other gas remains unbound.The latter receives a gravitational 'kick' and is launched across the cavity, as illustrated in Figure 1.As gas from the edge of the cavity is perturbed by the binary, material which leads the binary orbit causes a positive gravitational torque on the binary, and material which lags behind the binary orbit exerts a negative torque.The total torque on the binary is then derived from these gravitational torques and the .The evolution of binary mass ratios and semimajor axes in the same simulations shown in Figure 4. We find that all unequal-mass binaries are driven towards q = 1, although binaries accreting from α disks accrete towards equal mass at a faster rate.At intermediate mass ratios, where the gravitational torque is negative or small in magnitude (see Equation 19), binaries are driven to inspiral, while at other mass ratios the binaries can outspiral rather quickly.Notably, the rate of outspiral at low mass ratios varies strongly with viscosity, disks with constant kinematic viscosity leading to significantly faster outspirals.Numerical values are provided in Table 4.
angular momentum carried by gas as it is captured by the binary, as virtually all gas remains within the sphere of influence of the binary after entering (e.g.Tiede et al. 2022).This advected angular momentum then manifests as gravitational torques proximate to the binary and finally as an accretion torque.In M = 10 disks, the overall gravitational torque on the binary is usually positive, but in higher Mach number disks, proportionally more gas is launched through the cavity each pericenter passage than is accreted, leading to stronger negative torques on the binary in colder, thinner disks (Dittmann & Ryan 2022).
We illustrate in Figure 8 azimuthally averaged surface density and radially integrated torque profiles from our higher-mass ratio M = 30 simulations.To assist with interpretation, we also present time-averaged surface density profiles for a subset of these simulations, along with streamlines of the time-averaged momentum profiles, in the frame rotating with the binary, in Figure 9. Considering first the equal-mass binary, we note the net negative gravitational torque from r ≲ a b /2 shown in the bottom panel of Figure 8.This corresponds to the streams of gas which flow from one member of the binary to the other, illustrated in the top left panel of Figure 9, which lag behind the binary and exert a negative torque.From r ∼ a b /2 to r ∼ a b , a strong positive torque is exerted on the binary by gravitational interac- .Orbital properties of the inner edge of the circumbinary disk cavity/gap -the semi-major axis, scalar eccentricity, and argument of periapsis -from q = 0.01 to q = 1 at M ∈ {10, 20, 30} disks with ν = 0.001, averaged in time over the last 500 orbits of each simulation.Although all three quantities tend to decline as mass ratios decrease from q = 1, they do not so monotonically.Notably, at certain intermediate mass ratios at M = 20 and M = 30, the cavity becomes much smaller, more circular, and ceases to precess in a meaningful way; however, even in these cases the typical magnitude of the disk eccentricity is larger than for some lower-q binaries which still have well-defined, precessing, eccentric disks.Interestingly, at low mass ratios, some M = 10 disks precess in a retrograde rather than prograde sense.Numerical values are provided in Table 4.
tions between each binary component and the opposite minidisk, as well as the portions of accreting gas streams which lead the binary.The final appreciable component of the gravitational torque then arises between r ∼ a b and r ∼ 2a b , from gas which on average lags the binary, that which is stripped from the cavity walls only to be flung through the cavity rather than captured by the binary.
As binary mass ratios decrease, the accretion morphology becomes more axisymmetric, many of the torque components weaken, and the circumbinary disk gradually shifts inward.As the secondary moves outward, so does the radius at which mass typically flows between primary and secondary minidisks, as illustrated in Figure 9.This shift corresponds to the low-radius region of negative torques extending to larger radii, but becoming less prominent, as illustrated in Figure 8. Pronounced positive torques arise both in the minidisk about the secondary and from gas accreting from the circumbinary  6, but for our simulations of M = 10disks with viscosities ν = 0.001, ν = 0.0005, and α = 0.1.The lowest constantviscosity disks, at least for q ≳ 0.1, tend to have larger cavities and precess more slowly.However, at those same higher mass ratios disks with α = 0.1 tend precess appreciably faster than disks with ν = 0.001, even when they have smaller cavities, illustrating the important dynamical differences which arise from different spatial viscosity profiles.Generally, as binary mass ratio decreases, so does the disk precession rate, albeit non-monotonically.Similarly, at q ≲ 0.05, many disks exhibit retrograde precession, although curiously the α = 0.1 disk around a q = 0.01 binary precesses in a prograde sense.Numerical values are provided in Table 4.
disk at some radii.Material which accretes directly onto the primary tends to lead it, exerting a positive torque, while material which accretes onto the secondary leads the binary at smaller radii, but lags the binary at larger radii.At lower mass ratios, less material is perturbed by the secondary, and lower average surface densities within the cavity correspond to smaller torques.Similar trends occur for all disks around higher mass ratio binaries.

Thin disks and the absence of eccentric modes
At intermediate mass ratios, binaries accreting from M = 20 and M = 30 disks experience strikingly nonmonotonic gravitational torques at the same mass ratios where the eccentricity of their circumbinary disks also drops dramatically.These transitions occur at M = 20 between q = 0.04 and q = 0.14 (inclusive), and in M = 30 disks only at q = 0.05 and q = 0.07.As 10 0 a/a b 0 1 2 Σ q = 0.10 q = 0.14 q = 0.20 q = 0.30 q = 0.50 q = 1.00 . Average surface density (top panel) and cumulative gravitational torque acting on the binary (bottom panel, the gravitational torque density integrated from r ′ = 0 to r ′ = r, normalized by the accretion rate onto the binary) profiles corresponding to binaries with 0.1 ≤ q ≤ 1 accreting from M = 30 disks.Each profile was averaged over the final 500 binary orbits of each simulation.
Figure 9. Time-averaged, in the rotating frame of the binary, surface density profiles with momenta streamlines, for binaries with M = 30 disks at q ∈ {0.1, 0.2, 0.5, 1.0}.The equipotential contours of the first, second, and third Lagrange points are illustrated in each case using white curves.
In a time-averaged sense, the flow onto the primary from the circumbinary disk transitions from both lagging and leading the primary at high mass ratios, to primarily trailing the binary at lower mass ratios, contributing to a less negative overall gravitational torque as q decreases.illustrated in Figure 10 and the precession rates and cavity eccentricities plotted in Figures 6 and 7, the circumbinary disks corresponding to these extreme negative torques are not simply near-circular, but lack a quasi-global, precessing, eccentric mode throughout.In particular, a number of the disks around q = 0.01 binaries have even lower average eccentricities, but sill possess large-scale uniformly precessing eccentric modes.
The low disk eccentricities of and strong negative gravitational torques upon binaries at these mass ratios are illustrated in Figures 6 and 2 respectively.The absence of coherent eccentric modes in these systems, and the presence of such modes at both higher and lower mass ratios, is illustrated in Figure 10, which plots the azimuthally averaged disk argument of periapsis as a function of space and time for a set of binaries.The physics governing the growth, saturation, and precession rate of eccentric modes is beyond the scope of this work, but we plan to pursue a deeper understanding thereof in a future study.However, we observe steady precession, on periods of hundreds to thousands of binary orbits, in disks which support eccentric modes, and a lack of such precession in those which do not.Notably, as cavity semi-major axes and eccentricities drop as eccentric modes disappear at certain mass ratios, the pericenter of the cavity stays nearly constant in most cases (4).
Some of the resulting changes to the accretion flow are illustrated in Figure 11, which plots time-averaged surface density profiles for a pair of binaries with eccentric and non-eccentric disks at both M = 20 and M = 30.A few of the relevant binaries are also featured in Figure 1, particularly those at q = 0.03 and q = 0.05.Although not obvious from the time-averaged surface density profiles, circumbinary disks without eccentric modes at these mass ratios accrete much more onto the primary than those with eccentric modes.However, the morphology of this accretion changes dramatically, almost entirely lagging behind the primary rather than leading it and thus exerting a negative torque rather than positive.The streams of gas which accrete onto the secondary from disks without eccentric modes also appear to be higher-density, and more asymmetric azimuthally, lagging the binary more pronouncedly and further contributing to the strong negative torques in these systems.

Thick disks and mild torques
Thicker, M = 10 disks also exhibit torques which become slightly negative at intermediate mass ratios, which are often accompanied by changes in cavity size and eccentricity (c.f.Figures 4 and 7  .Space-time diagrams of the disk argument of periapsis over the final 500 orbits of each simulation, averaged azimuthally in bins of constant fluid element semi-major axis, with the cavity semi-major axis marked using a dashed white line.Rarely is the entire disk precisely in phase.However, for each Mach number shown here, M = 20 in the top row and M = 30 in the bottom row, disks around both higher-and lower-mass binaries precess at similar rates throughout.However, disks around binaries at the intermediate mass ratios shown here lack a uniformly precessing eccentric mode, and their eccentricity at small semi-major axes is dominated by orbital-timescale variability.P orb is the binary orbital period 2πΩ −1 b .We used t0 = 1500 P orb for simulations of binaries with q ≥ 0.1 and t0 = 500 P orb for simulations with q < 0.1.in α = 0.1 and ν = 0.001 disks respectively, smaller torques on binaries with q = 0.2 than those with q = 0.1 of q = 0.3.We have reproduced these results herein, and found similar behavior in disks with ν = 0.0005 at q = 0.1.While it is tantalizing to guess that the same mechanisms at play in the previously discussed M = 20 and M = 30 disks might be at play here, these M = 10 disks still exhibit uniformly precessing eccentric modes even when the magnitude of their eccentricity drops at these mass ratios.Additionally, while the cavity size and eccentricity change, they do not do so while keeping nearly constant preicenter as did the higher Mach number disks (Table 4), and the secondary actually accretes more preferentially as the torque drops whereas more accretion onto the primary occurred from more circular disks in the higher Mach number cases.
To illuminate some of the causes of these trends in gravitational torque, which occur at different mass ratios for different disk viscosities, we plot in Figure 12 a few profiles of surface density and integrated gravitational torque density for pertinent mass ratios.For the ν = 0.001 disk, the change in cavity size between the q = 0.3 and = 0.2 binaries is similar to that between the q = 0.14 and q = 0.1 binaries at ν = 0.0005, suggesting that the same mechanism governs both cases.
The integrated torque profiles in all cases largely follow each other with nearly constant offsets beyond r ≳ a b , suggesting that the torque contributions from larger radii, including the negative torques dominant from a b ≲ r ≲ 2a b , are largely similar for each mass ratio and viscosity.However, changes to the morphology at small radii occur according to the angular momentum carried by the infalling gas.Upon examination of timeaveraged surface density profiles, shown in Figure 13, it appears that as the mass ratio drops, an appreciable amount of gas which flows about the primary finds its way into the Roche lobe of the secondary, not remaining bound to the primary.Additionally, at the lower mass ratios in these cases, less gas seems to flow from the secondary to the primary, in accordance with the larger values of Ṁ2 / Ṁ1 observed at the lower mass ratios of the binaries considered here.

Trends in preferential accretion
At all non-unity mass ratios, we observe that the secondary accretes at a higher rate than the primary.While many previous studies have identified this trend at M = 10, recent studies have found various degrees of (dis-)agreement in the magnitude of Ṁ2 / Ṁ1 .Despite prior studies using constant-α viscosity prescriptions  10 in which eccentric modes appear at intermittent mass ratios.As binaries accrete from disks with prominent eccentric modes (q = 0.2 M = 20 and q = 0.04 M = 30 here), gas that flows onto the primary more often leads it, contributing to a more positive torque.Additionally, disks with eccentric modes accrete in a more spatially uniform manner onto the secondary, while the secondary-bound accretion streams from disks without eccentric modes tend to lag behind the binary.The aforementioned effects conspire to create a much more negative gravitational torque on binaries accreting from disks; the difference in torque follows from a fundamental change in the accretion flow morphology, and the primaries actually accrete (proportionally) more from disks without eccentric modes than from disks with prominent eccentric modes.
(e.g., Muñoz et al. 2020;Siwek et al. 2023), constantν viscosity prescriptions (e.g., Dittmann & Ryan 2021), and heterodox formulations of the equations of viscous hydrodynamics (Duffell et al. 2020), the role of viscosity as been heretofore unappreciated: α-viscosity, at least in M = 10 disks, enhances the degree of preferential accretion onto the secondary, as shown in Figure 4. Additionall, at very low mass ratios, we find a higher degree of preferential accretion in colder disks, as shown in Figure 2.

Viscosity
When equal partition of material between the primary and secondary is not guaranteed by symmetry, α viscosity causes higher rates of accretion onto the secondary relative to those onto the primary than in comparable simulations of disks with constant kinematic viscosity.Qualitatively, from examination and comparison of the Figure 13.Time-averaged profiles of the disk surface density, with streamlines of the time-averaged momentum field, in the frame of the binary for a few binaries accreting from M = 10 disks.At lower mass ratios more gas which is stripped from the cavity by the primary eventually becomes bound to the secondary, and less gas flows from the minidisk of the secondary onto the primary.These morphological changes coincide with decreases in the gravitational torques at q = 0.2 for these systems. .Accretion rate time series, over fifteen orbital periods, for a representative sample of binaries of mass ratios.In each panel, the accretion rate onto the primary is plotted in blue ( Ṁ1/ Ṁ0), and the accretion rate onto the secondary is plotted in orange ( Ṁ2/ Ṁ0), in both cases normalized by the equilibrium accretion rate for each simulation.At lower mass ratios (q ≲ 0.1), the accretion rate onto the primary is fairly constant, and well below the accretion rate onto the secondary.
constant-α and constant-ν binaries presented in Figure 13, it appears that less matter flows between the primary and secondary in α−viscous disks than in constant-ν disks, while the amount of gas flowing directly onto the primary through the cavity does not appear appreciably different.It may be that the relatively flat surface density profiles of constant-ν disks in Newtonian potentials allow more material to flow between the binary members, while the more centrally-peaked surface density profiles that arise from α-viscosity prescriptions (Σ ∝ r −1/2 in a Newtonian potential) allow less gas to flow from the minidisk of the secondary to that of the primary.
In an attempt to test this, we have carried out an additional simulation at q = 0.1, M = 10 using an ad hoc viscosity prescription, 3 defining ν = ℵc s ‫,ח‬ setting ℵ = 0.1 and defining (26) 3 Although a different radial viscosity profile could be achieved within an α-viscosity framework by varying the sound speed profile of the disk, such changes also alter gas dynamics even for constant viscosity profiles, as illustrated in Dittmann & Ryan (2022), so we have restricted ourselves to a simpler scenario.
in analogy to Equation 6 but with a viscosity profile ν ∝ r far from the binary and ν ∝ r i within each minidisk, resulting in a more sharply peaked surface density profile within the minidisks.We find that the average degree of preferential accretion in this case, Ṁ2 / Ṁ1 over the final 500 orbits of the 2000-orbit simulation, is 7.06, compared to ∼ 5.9 for an α−viscosity and ∼ 4.4 for a globally constant viscosity at the same Mach number and mass ratio.While this can not rule out that the viscosity profile within the cavity plays some additional role, it further demonstrates that the spatial viscosity profile governs accretion onto the binary.

The rate of accretion onto low-mass objects
As shown in Figure 2, we find that at q ≲ 0.03 the secondaries accrete more preferentially in colder disks.The rate at which mass accretes through some surface can often be parameterized by Ṁ ∼ ρ c A c v c , where ρ c , A c , and v c are characteristic densities, areas, and velocities respectively.When the flow of gas around an object is spherically symmetric, the characteristic speed is typically taken to be the sound speed c s , and the area to be the surface area of a sphere with a radius given by R B = 2GM/c 2 s , roughly balancing thermal and gravitational energy, resulting in the Bondi accretion rate ṀB ∝ ρM 2 c −3 s (Bondi 1952).However, when the Hill radius of an embedded object becomes smaller than the Bondi radius, the surface area of the Hill sphere becomes Figure 15.The relative standard deviation of the accretion rate onto the primary (δ Ṁ1) and onto the secondary (δ Ṁ2).Accretion is never steady in an absolute sense.However, in cases where the cavity eccentricity is below e ≲ 10 −3 or so, the level of normalized variability is below ∼ 1%.At most mass ratios, the degree of variability is greater for binaries in colder disks.
a natural characteristic area.Similarly, when the Hill sphere becomes larger than the disk scale height, the characteristic area becomes 2πHr H .The latter regime characterizes all of the binaries studied herein.
Thus, we might expect an accretion rate onto the secondary along the lines of Ṁ ∝ ρHr H v c ∝ Σr H v c .Whether the characteristic velocity is taken to simply be the sound speed, or related to the Hill radius and Keplerian shear as in Rosenthal et al. (2020), the accretion rate onto the secondary increasing as disks cool contradicts these expectations.We note that such predictions appear to be fairly accurate when measuring the rate at which material enters the sphere of influence of objects in disks (defined directionally, only measuring matter which is moving towards the secondary rather than that moving away from it) (Choksi et al. 2023).However, we have measured not just the rate at which gas enters the sphere of influence of the secondary, but the rate at which it accretes onto the secondary.Indeed, it appears that the accretion streams onto the secondary are more pronounced in higher Mach number disks, and that the corresponding cavities are larger and deeper, preventing accretion onto the primary, most likely because of weak shocks which become stronger in cooler disks.

Variability
In the high mass ratio (q ≳ 0.2) regime the accretion rate onto the binary is often governed by interactions with the pronouncedly eccentric cavities of their circumbinary disks and the overdense regions on the inner cavity edge, leading to modulation of the accretion rate on the orbital period of the binary (or half, in the case of equal-mass binaries) as binary members pass by the cavity pericenter, on the orbital period of the inner edge of the cavity and its overdensites, and the associated harmonics and beat frequencies (e.g., MacFadyen & Milosavljević 2008;Shi et al. 2012;Noble et al. 2012;D'Orazio et al. 2013;Shi & Krolik 2015;Miranda et al. 2017;Bowen et al. 2019;Muñoz et al. 2020;Dittmann & Ryan 2021;Noble et al. 2021;Westernacher-Schneider et al. 2022;Combi et al. 2022;Krauth et al. 2023).
Some studies (e.g., Duffell et al. 2020) have argued that below q ≲ 0.2 accretion variability on supraorbital timescales disappears, or that below q ≲ 0.04 the accretion is characterized by steady-state behavior rather than fluctuations (e.g., D'Orazio et al. 2016): based on our higher-fidelity simulations, 4 we find that significant variability persists at least to q = 0.01, and that in some cases this can occur on the orbital period of the cavity rather than that of the binary.To illustrate this, we show in Figure 14 a sequence of accretion rate time series, at q ∈ {0.01, 0.03, 0.1, 0.3, 1.0} and M ∈ {10, 20, 30}.Evidently, below q ≲ 0.1, the accretion rate onto the primary, and its variability, become negligible compared to that onto the secondary.However, in most cases the accretion rate onto the secondary is still highly variable.Noting that at M = 10, the secondary of the q = 0.01 binary exhibits accretion variability on the orbital period of the cavity, while the q = 0.03 binary exhibits no appreciable variability, it appears, based on Figure 6, that the critical disk eccentricity below which the accretion flow becomes steady, if one exists, is e ≲ 10 −3 .Although the orbital period of the cavity dominates the accretion rate periodicity of the the q = 0.01, M = 10 binary, oftentimes the most prominent period at q ≲ 0.1 is that of the binary orbit, although this is not to say that other periods are not present.
In an absolute sense, the accretion flow onto neither object is ever steady.However, in a number of lowmass ratio cases exhibit only very minute variability, to the point of being nearly imperceptible in visualizations such as Figure 14.To more quantitatively assess variability, we have calculated for each binary the relative standard deviation (also known as the coefficient of vari-  (2020), and in addition to third-rather than second-order spatial reconstruction, we have also used more robust sink methods and time-averaged relevant quantities between outputs rather than outputting point values in time.Notably, Dittmann & Ryan (2021) illustrated that the sink methods used in prior studies led to spurious accretion variability.
ation) for the accretion rate onto each object (δM i ), the standard deviation divided by the mean of each time series.Figure 15 displays this metric for each accreting object.Thinner disks typically exhibit fluctuations on the order of ∼ 10% or more for at least the secondary, if not both objects.Thicker, M = 10 disks, on the other hand, can exhibit appreciably less accretion variability, particularly at q ≤ 0.05, although the q = 0.01 binary accreting from a ν = 0.001 disk still exhibits accretion rate variability at the ∼ 10% level.While these low-variability cases often correspond to binaries with low-eccentricity disks (e ≲ 10 −2 , Figure 7), other disksboth thick and thin -with similarly small eccentricities lead to orders of magnitude greater variability at the same mass ratios.

DISCUSSION
Our simulations provide insight into the evolution high-mass protoplanets, supermassive black hole binaries, and stellar binaries, among other systems.While the survey was performed at high resolution and covered a range of Mach number, disk viscosity, and mass ratio, many limitations remain.In particular we neglected binary eccentricity, used the overly simplistic locally isothermal equation of state, and, most importantly, performed the simulations using only two dimensional hydrodynamics.That being said, these calculations can still provide general insight into a number of astrophysical systems and their dependence on physical parameters.

Stars and Planets
The Gaia mission has unveiled a substantial population of wide binaries stars with separations reaching ∼ 1000 au and mass ratios near unity (q ≳ 0.95) (e.g., El-Badry et al. 2019), many of which are thought to have formed in and been shaped by circumbinary disks (Hwang et al. 2022).While accretion from all but the hottest, thickest disks (such as those studied by Ochi et al. 2005, for example) should cause the binary mass ratio to approach unity, the widening of binaries appears to have more sensitive dependencies.For thicker, M ≲ 10 disks, it appears that, as discussed in Section 3.2.2,binaries that form with mass ratios q ≳ 0.1 − 0.3 may spiral outward due to interactions with their circumbinary disk, depending on the effective disk viscosity.
Concerning the interactions between planets and circumstellar disks, gravitational torques in the linear, very low mass regime have been thoroughly characterized, oftentimes leading to negative torques on the planet (e.g., Goldreich & Tremaine 1979;Tanaka et al. 2002;Paardekooper & Papaloizou 2009).Further complications arise as planets become more massive, strongly perturbing the gas around them and potentially exciting eccentricity in the disk (e.g., Kley & Dirksen 2006;Kanagawa et al. 2018;Dempsey et al. 2021).In the fairly high-mass regime studied here, relevant to the evolution of brown dwarfs and super-Jupiters (0.01 ≲ q ≲ 0.1), we have identified both positive torques, associated with the presence of long-lived precessing eccentric disks (correlating with but not depending on the amplitude of the disk eccentricity), and negative torques which transpire in the absence of these modes.The positive torques identified at q ∼ 0.01 may help planets formed through core accretion to move outwards to large separations, and outward migration may facilitate giant planet formation from otherwise under-massive cores (e.g., Piso & Youdin 2014).

The Evolution and Detectability of SMBH binaries
Studies of equal-mass binaries have suggested that thin accretion disks, as often invoked for active galactic nuclei as predicted by simple models of optically thick disks accreting at below the Eddington rate (e.g., Shakura & Sunyaev 1973), 5 have identified very strong negative torques, which could potentially drive binaries into the gravitaitonal wave-dominated regime of orbital evolution within the duration of a typical accretion episode (Dittmann & Ryan 2022).However, as binary mass ratios decrease, these torques become appreciably weaker: at M = 20, d log a b /d log M increases from ∼ −1.5 at q = 1 to 1.2 at q = 0.2; and at M = 30, d log a b /d log M increases from ∼ −6.5 at q = 1 to ∼ −2.2 at q = 0.1.In some cases, for binaries around which eccentric disk modes are suppressed may inspiral at even faster rates, although we have not determined the reason behind this behavior.While near-equal-mass binaries formed from major mergers may inspiral quickly and those with mass ratios q ∼ 0.3 more slowly, there is a chance that those with mass ratios q ∼ 0.05, like those predicted to result most frequently from galaxy mergers (e.g., Amaro-Seoane et al. ( 2022); Bellovary et al. (2019), but see also Tamfal et al. (2018)), inspiral extremely quickly (d log a b /d log M ∼ −35).However, we also find that at lower mass ratios, q ≲ 0.04, binaries often outspiral due to interactions with circumbinary disks, which may cause such binaries to stall.
5 Notably, near-Eddington or very below-Eddington accretion rates are expected to beget fairly thick accretion disks (e.g., Narayan & Yi 1995;Sądowski et al. 2011).Disks with sufficient heating to avoid gravitational instability at large radii are also oftentimes predicted to be relatively thick (e.g., Sirko & Goodman 2003;Dittmann & Miller 2020;Gilbaum & Stone 2022) Active galactic nuclei are often inferred to accrete at rates between ∼ 1 and ∼ 1/100 times the Eddington limit (e.g., Shen et al. 2008).Assuming that this holds true for binary systems as well, preferential accretion onto the secondary by the factors of ∼ 2 − 7 typically observed herein suggests that the more massive of the SMBHs typically accretes at an extremely sub-Eddington rate, due to both the tenuous supply of gas and higher mass.The resulting circumprimary accretion flows may be advection-dominated (e.g., Narayan & Yi 1995), or emit much more subtle radiation, making detection more challenging.We also note that because unequal-mass binaries often exhibit prominent accretion rate variability on their orbital period, such variations may be easier to robustly detect than the longertimescale cavity-dominated variability that dominates in near-equal-mass binaries.

The importance of viscosity
Interestingly, the functional form of the Navier-Stokes viscosity strongly affects the rate of preferential accretion onto the secondary (Figure 4), and at some mass ratios can affect the rate and direction of the evolution of the binary semi-major axis (Figure 5).Although this clarifies some of the differences between prior investigations into the problem of unequal-mass binaries (e.g., Muñoz et al. 2020;Dittmann & Ryan 2021), it is somewhat unfortunate, because neither α−viscosity nor a constant kinematic viscosity (of the magnitude considered herein) is truly present in astrophysical accretion disks: because the results of these simplified simulations are dependent on the viscous model employed, simulations with realistically driven angular momentum transport, via processes such as magnetohydrodynamic turbulence (e.g., Velikhov 1959) or gravitational instability (e.g., Gammie 2001;Rafikov 2015) will be necessary to assuredly understand preferential accretion, although in a broad sense the evolution of the binary semi-major axis appears less sensitive to viscous model assumptions.Recent advancements in numerical techniques may facilitate such studies, particularly for studies of the late inspirals of binary black holes (e.g., Avara et al. 2023), and novel stochastic viscosity prescriptions may also hold promise (Turner & Reynolds 2023).

Disk Eccentricity and Orbital Stability
In a number of cases, examined in Section 3.2.1, the circumbinary disk fails to develop a well-defined eccentric mode throughout the disk, an absence that is correlated strongly with changes in accretion morphology, leading to pronouncedly negative torques on the binary.Similar observations have been made previously, albeit with less precision (e.g., Dempsey et al. 2021), the key difference being that disks may have smallermagnitude eccentricities but still possess quasi-global eccentric modes and exhibit positive gravitational torques.Additionally, although disks with low eccentricities oftentimes coincide with small degrees of accretion variability, as shown in Section 3.4, this is by no means universal.
Orbits about the triangular Lagrange points in the circular restricted three-body problem are linearly unstable for binary mass ratios q > (25 − √ 621)/2 ≈ 0.04 (e.g., Danby 1988), and for this reason, we focused disproportionately on proximate mass ratios.Examining Figures 2 to 7, and 15, the only changes we see associated specifically with q = 0.04 occur in colder (M = 20, M = 30) disks, where eccentric modes are present at q = 0.03 but not at q = 0.05.
A previous paper, D 'Orazio et al. (2016) suggested that the above instability threshold has drastic consequences, causing a transition from a steady flow at q ≲ 0.04 to a strongly fluctuating flow above that mass ratio.D'Orazio et al. ( 2016) also suggested that this critical mass ratio marked a transition from a nearlycircular circumbinary disk with a narrow annular gap below q ≈ 0.04 to an eccentric disk with a large hollow cavity above.Our M = 10 simulations bear some resemblance to this picture, although they seem appear to change gradually below q ≲ 0.1.However, our simulations of thinner disks display qualitatively different behaviour, maintainting eccentric disks and accretion rate fluctuations greater than 10% even down to q = 0.01.However, we do confirm that in a broad sense, disk eccentricity, cavity size, and disk precession rate decrease along with the binary mass ratio.

SUMMARY
We have carried out a broad survey of circumbinary accretion, studying circular binaries with mass ratios between 0.01 and 1; disk aspect ratios between 0.1 and 0.0 33 (azimuthal Mach numbers between 10 and 30), using both a constant-ν viscosity prescription (with ν = 0.0005a 2 b Ω b and ν = 0.001a 2 b Ω b ) and an α viscosity prescription.Although two-dimensional, our simulations have used very high resolutions, down to δx ∼ 0.00165a b , and third-order spatial reconstruction in addition to robust sink methodology -in fact, our tests in Appendix A indicate that without using torque-free sinks (Dempsey et al. 2020;Dittmann & Ryan 2021), it is impossible to achieve converged results.Our main results are: • The spatial dependance of viscosity strongly determines the relative distribution of accretion be-tween primary and secondary, likely because the more centrally concentrated surface density profiles reduces the amount of gas which flows from the Roche lobe of the secondary to that of the primary.
• Cold (M ≥ 20) disks that do not develop quasiglobal precessing eccentric modes deliver very strong negative gravitational torques to the binary.
• Cavity eccentricity and the variability of accretion onto the binary broadly correlate with one another, and cavity eccentricities below e ∼ 10 −3 are often needed to preclude variability.
• The strong negative gravitational torques observed for circular binaries and thin disks become weaker as the binary mass ratio decreases.
• We find little evidence at these Mach numbers and viscosities that the stability of orbits about the fourth and fifth Lagrange points in the circular restricted three body problem plays a significant role in the dynamics of circumbinary accretion.
These results should further our understanding of the orbital evolution of giant planets, binary stars, and supermassive black holes.

CONVERGENCE
We have performed a small set of convergence tests at mass ratios q = 0.01 and q = 0.1, the results of which suggest that our fiducial simulation parameters are sufficient to derive sufficiently converged results.First, at q = 0.1, M = 10, ν = 0.001Ω b a 2 b we performed and additional simulation at a base resolution of 512 cells in both the x and y directions, resulting in a resolution near the binary of ∆x ≈ 0.0049a b .As a cross-check between our different refinement schemes, we also performed a simulation with the aforementioned physical parameters but the grid, sink, and softening configuration employed in our lower mass ratio simulations.The results of these tests are displayed in Table 1.The total torque on the binary in very well converged, varying by < 1% between simulations.At constant gravitational softening length, the degree of preferential accretion is also very well converged with resolution, varying by < 1%.However, Ṁ2 / Ṁ1 appears to depend fairly strongly on the gravitational softening length, decreasing as the softening length scale decreases.One possible explanation for this trend is that, similarly to α−viscosity but in this case because of the decreased shear in the inner disk, gravitational softening leads to a steady-state surface density profile which decreases with radius which may reduce the rate of mass transfer between the Roche lobe of the primary and secondary.
Because certain mass removal algorithms and parameter choices can lead to unphysical biases (Dempsey et al. 2020;Dittmann & Ryan 2021), we test the extent to which our choices of sink rate γ i affect our results.Although our primary suite of simulations all used torque free sinks, setting δ = 0 in Equation (10), we carried out a limited number of ν = 0.001, q = 0.1, M = 10 simulations to test the extent to which use of standard (δ = 1) sinks may have contributed to discrepancies in the values of Ṁ2 / Ṁ1 previously reported in the litera-  1.The results of our resolution tests.The torque on the binary is converged regardless of numerical parameters.However, the degree of preferential accretion depends fairly strongly on the gravitational softening used in each simulation.
ture (e.g., Farris et al. 2014;Muñoz et al. 2020;Siwek et al. 2023).Time-averaged values of Ṁ2 / Ṁ1 and Jg / Ṁ are shown in Figure 16.Although the values measured in simulations which used torque-free sinks are broadly consistent with one another, results derived using standard sinks depend strongly on the sink rate.At low sink rates, torque-free and standard sinks are in better agreement, but if the sink rate is too low (in this case γ ≲ 4/3) the results become erroneous, as the sink is unable to remove mass sufficiently quickly and artificial density gradients develop.Figure 17 illustrates time-averaged density profiles of the minidisks in these simulations.Notably, at slower sink rates, the surface density profile from the standardsink simulations approaches that of the simulations which used torque-free sinks.Although the simulations using standard sinks display no semblance of convergence with sink rate, the surface density profiles in simulations using torque-free sinks and sink rates γ ≥ 8/3 approximately converge far from the sinks.Furthermore, the surface density profile around the primary in the γ = 4/3 simulation, which deviates from the others well beyond the sink region, suggests that mass is being removed very slightly too slowly, allowing the disk structure outside of the sink to be affected.Thus, we consider both our default algorithmic choices and sink rates to be well-justified.
We performed a similar test at q = 0.01, for a disk with ν = 0.001 and M = 10.Using out fiducial resolution, sink size, and softening parameters, we tested the sink rates γ 0 = {4/3, 8/3, 16/3}.Taking our fiducial sink rate, γ 0 = 8/3 as a baseline, we found that the accretion rate ratio Ṁ2 / Ṁ1 and gravitational torque Jg / Ṁ agreed to better than 3% between the γ 0 = 1.33 and γ 0 = 2.67 runs, and agreed to better than 1% between the γ 0 = Figure 17.Time-and azimuth-averaged surface density profiles in the frame of the primary (top panel) and secondary (bottom panel).Different line colors, from dark to light, indicate different sink rates from slow to fast.Solid lines plot profiles derived from simulations using standard sinks and dashed lines plot profiles derived from simulations using torque-free sinks, which are the default for the simulations in this work.We observe that standard sinks significantly deplete the surface density profile well beyond the sink radius (0.035a b ), and the profiles do not converge at any value of the sink rate.Torque-free sinks, on the other hand, appear to produce converged surface density profiles throughout the bulk of the minidisks at γ ≥ 8/3.

B. COMPARISONS WITH PREVIOUS WORK
Over the last few years, a few studies have investigated the orbital evolution of binaries with q ∼ 10 −2 and found results in disagreement with what we have reported in this work.Namely, when nominally studying binaries with 10 −2 ≤ q ≤ 1, in M = 10, ν = 0.001 disks, Duffell et al. (2020) found not only that binaries inspiral at q < 0.1 -where we find that binaries should outspiral -but also values of Ṁ2 / Ṁ1 appreciably higher than in our simulations.Additionally, Dempsey et al. (2021) reported negative gravitational torques in a simulation of a q = 0.01 binary with a M = 10 disk.Below, we identify some of the modeling choices responsible for these discrepancies.

B.1. Comparison with Duffell et al. 2020
The present study differs from Duffell et al. (2020) in a number of ways: the code utilized (Athena++ vs Disco), the simulation initial conditions, the simulation diagnostics, and the fact that Duffell et al. (2020) solved a fundamentally different set of equations describing viscous hydrodynamics.To investigate this, we have carried out a small suite of simulations using different versions of the Disco code.We find results which are consisten with our main, Athena++-derived, results when solving the same set of equations and reporting similar diagnostic quantities.
The most important factor contributing to the differences in the values of Ṁ2 / Ṁ1 measured in our simulations and those in Duffell et al. (2020) is that we solved the Navier-Stokes (NS) equations, given by Equation (2), with the velocity shear tensor defined by Table 2.The results from a series of simulations of accretion onto q = 0.02 binaries from M = 10, ν = 0.001 disks, using different versions of Disco, solving both the Navier-Stokes equations and the set of equations solved in Duffell et al. (2020), and at different levels of numerical dissipation.Gravitational torques are presented in units of 3πνΣ0a 4 b Ω b , as in Duffell et al. (2020), rather than units of Ṁ a 2 b Ω b as in every other section of this paper.For our simulations using Disco, we report per-timestep time averages of each quantity over the last 500 binary orbits; because Disco v1 lacks this functionality, for those simulation we report the mean of point values which were output 100 times per orbit, over the last 500 orbits.When using Disco to solve the Navier-Stokes equations, we find consistent results with our Athena++ simulations.When solving a modified equation set and reproducing errors in torque calculations, we are able to reproduce the results of Duffell et al. (2020).
simulations solving the NS * equations find systematically higher values of Ṁ2 / Ṁ1 , which also depend more strongly on numerical dissipation than results from simulations solving the NS equations.
Some more recent studies (Dittmann & Ryan 2021, 2022) have used an updated version of Disco, which solves the NS equations (as well as the NS * equations, for the sole purpose of making comparisons to Disco v1), and includes a number of other upgrades such as using cell centroids rather than centres where appropriate and time-averaging outputs every time step properly throughout the Runge-Kutta time integration.We denote this version of Disco simply as Disco.7 Our simulations using both Disco and Disco v1 indicate that the largest differences between the gravita-tional torques measured in our Athena++ simulations and those reported in Duffell et al. (2020) are due to errors in the simulation diagnostics used in Duffell et al. (2020), which we have identified in the Disco v1 source code. 8Specifically, Disco v1 uses the lever arm of the secondary, in the barycentric frame, to calculate the torque on both the secondary and the primary.Our code comparison and source code examination also suggests that Duffell et al. (2020) mislabeled the gravitiatonal torques on the primary and secondary, matching the variable names in the source code.We denote a verison of Disco v1 that corrects these typos as Disco v1 * .As shown in Table 2, we find results in good agreement with Duffell et al. (2020) and the base version of Disco v1.However, if we correct the typos in the torque calculation and reporting, we find good agreement between Disco v1 * , Athena++, and Disco.
When attempting to compare our results to those of Duffell et al. (2020), we have only focused on a single mass ratio, q = 0.02 rather than adiabatically varying the mass ratio.We use the same initial conditions described in Section 2.2, albeit using a different mesh: rather than a Cartesian domain extending from −10a b to 10a b in x and y with static mesh refinement as in our Athena++ simulations, we use a cylindrical domain which extends from r = 0 to r = 30, with approximately uniform resolution within r < a b and approximately log-uniformly spaced cells from a b ≤ r ≤ 30a b , setting the number of azimuthal cells at each radius so that r∆ϕ ≈ ∆r.All simulations presented here used the HLLC Riemann solver and set the motion of each grid annulus to match the average azimuthal velocity of the fluid elements therein.To test the effects of numerical dissipation, particularly in the NS * formalism, we tested multiple resolutions (N r = 512, as in Duffell et al. (2020), and N r = 768) as well as different slope limiting coefficients θ PLM , where θ PLM = 1 corresponds to a more numerically dissipative approximation of the slope in a piecewise linear reconstruction and θ PLM = 1.5 leads to a slightly less dissipative, but still total variation diminishing, reconstruction using a generalized minmod formulation (see Kurganov & Tadmor 2000;Duffell 2016).For the sake of comparison, all of these simulations used standard (δ = 1) sinks with γ = 3 and r s = 0.05a b , as in Duffell et al. (2020).Thus, these 8 See https://github.com/duffell/Disco/blob/833bb65ab937bc8e18d109cd87afa95ddda5880c/report.c#L149 and https://github.com/duffell/Disco/blob/833bb65ab937bc8e18d109cd87afa95ddda5880c/Planets/bin_ varyq.c.We note that Disco v1 names the torque on the secondary "Torque," and the torque on the primary "Torque2."3. The results from a series of simulations attempting to reconcile the differences between our results and those of Dempsey et al. (2021) for q = 0.01 binaries.Both a large softening length and disallowance of accretion onto the secondary are necessary to reproduce the results of Dempsey et al. (2021).
results can not be converged, but may at best appear so by chance (see Appendix A).

B.2. Comparison with Dempsey et al. 2021
One of the primary results of Dempsey et al. ( 2021) was the transition from inspiral at lower mass ratios to outspirals at high mass ratios.More precisely, Dempsey et al. (2021) found that, in terms of the gap-width scaling parameter K ′ = q 2 α −1 h −1 (the ratio of the onesided gravitational torque in the absence of a gap to the viscous torque, scaled by the square of the disk aspect ratio h = H/r, e.g., Lin & Papaloizou 1986;Dempsey et al. 2021), binaries with K ′ ≲ 10 experience negative gravitational torques whereas those with K ′ ≳ 10 experience positive gravitational torques (See Figure 2 in Dempsey et al. (2021)).However, as an example, we have found positive torques in all of our simulations at q = 0.01, which have K ′ ≤ 1.
One factor which may have contributed to Dempsey et al. (2021) finding negative gravitational torques in all simulations with K ′ ≲ 10 is their omission, as detailed in their first Appendix, of terms in the gravitational torque related to performing simulations in the stellocentric frame.The neglected terms can be sufficiently large to lead to positive torques in some simulations for which they reported negative torques.
Two other factors contributing to the difference in results between our study and Dempsey et al. (2021) are that while we allowed both point masses to accrete and used gravitational softening values that were much smaller than the Hill radius of a given point mass, Dempsey et al. (2021) did not allow accretion onto the secondary and also used much larger values of the softening, 0.06 for M = 10 disks, following Müller et al. (2012).We have conducted an additional set of simulations of q = 0.01, M = 10, ν = 0.001 systems, the results of which are presented in Table 3, and determined that both larger gravitational softening values and disallowance of accretion onto the secondary are necessary to produce negative gravitational torques in such systems.

C. SIMULATION SUMMARIES
For general convenience, we have collected a number of results concerning disk morphology and binary evolution from our main suite of simulations in Table 4.
a b /d log M , and dq/d log M ) as functions of mass ratio in Figures 2 and 3, for simulations which varied M ∈ {10, 20, 30} at constant kinematic viscosity ν = 0.001, and in Figures 4 and 5 for simulations of M = 10 disks that employed Navier-Stokes viscosities of different magnitudes and spatial variations.Figures 6 and 7 plot the similarly partitioned cavity semi-major axis, eccentricity, and precession rate.Numerical values are collected in Appendix C.
Figure3.The evolution of binary mass ratios and semimajor axes in the same simulations shown in Figure2.At larger mass ratios, we find results in general agreement with previous studies.However, at intermediate mass ratios, which have not been previously probed at M ̸ = 10, we find extremely rapid inspirals, following from both the appreciably negative torque, small specific angular momentum of the binary, and that λ > q.At very low mass ratios (q ≲ 0.03), the rate of change of the binary semi-major axis can be quite large unless the torque is very nearly balanced by the change in the mass of the binary.Numerical values are provided in Table4.
Figure5.The evolution of binary mass ratios and semimajor axes in the same simulations shown in Figure4.We find that all unequal-mass binaries are driven towards q = 1, although binaries accreting from α disks accrete towards equal mass at a faster rate.At intermediate mass ratios, where the gravitational torque is negative or small in magnitude (see Equation19), binaries are driven to inspiral, while at other mass ratios the binaries can outspiral rather quickly.Notably, the rate of outspiral at low mass ratios varies strongly with viscosity, disks with constant kinematic viscosity leading to significantly faster outspirals.Numerical values are provided in Table4.
Figure6.Orbital properties of the inner edge of the circumbinary disk cavity/gap -the semi-major axis, scalar eccentricity, and argument of periapsis -from q = 0.01 to q = 1 at M ∈ {10, 20, 30} disks with ν = 0.001, averaged in time over the last 500 orbits of each simulation.Although all three quantities tend to decline as mass ratios decrease from q = 1, they do not so monotonically.Notably, at certain intermediate mass ratios at M = 20 and M = 30, the cavity becomes much smaller, more circular, and ceases to precess in a meaningful way; however, even in these cases the typical magnitude of the disk eccentricity is larger than for some lower-q binaries which still have well-defined, precessing, eccentric disks.Interestingly, at low mass ratios, some M = 10 disks precess in a retrograde rather than prograde sense.Numerical values are provided in Table4.

Figure 7 .
Figure7.The same disk properties plotted in Figure6, but for our simulations of M = 10disks with viscosities ν = 0.001, ν = 0.0005, and α = 0.1.The lowest constantviscosity disks, at least for q ≳ 0.1, tend to have larger cavities and precess more slowly.However, at those same higher mass ratios disks with α = 0.1 tend precess appreciably faster than disks with ν = 0.001, even when they have smaller cavities, illustrating the important dynamical differences which arise from different spatial viscosity profiles.Generally, as binary mass ratio decreases, so does the disk precession rate, albeit non-monotonically.Similarly, at q ≲ 0.05, many disks exhibit retrograde precession, although curiously the α = 0.1 disk around a q = 0.01 binary precesses in a prograde sense.Numerical values are provided in Table4.

Figure 10
Figure10.Space-time diagrams of the disk argument of periapsis over the final 500 orbits of each simulation, averaged azimuthally in bins of constant fluid element semi-major axis, with the cavity semi-major axis marked using a dashed white line.Rarely is the entire disk precisely in phase.However, for each Mach number shown here, M = 20 in the top row and M = 30 in the bottom row, disks around both higher-and lower-mass binaries precess at similar rates throughout.However, disks around binaries at the intermediate mass ratios shown here lack a uniformly precessing eccentric mode, and their eccentricity at small semi-major axes is dominated by orbital-timescale variability.P orb is the binary orbital period 2πΩ −1 b .We used t0 = 1500 P orb for simulations of binaries with q ≥ 0.1 and t0 = 500 P orb for simulations with q < 0.1.

Figure 11 .
Figure11.Time-averaged profiles of the disk surface density, with streamlines of the time-averaged momentum field, in the frame of the binary for a few of the simulations depicted in Figure10in which eccentric modes appear at intermittent mass ratios.As binaries accrete from disks with prominent eccentric modes (q = 0.2 M = 20 and q = 0.04 M = 30 here), gas that flows onto the primary more often leads it, contributing to a more positive torque.Additionally, disks with eccentric modes accrete in a more spatially uniform manner onto the secondary, while the secondary-bound accretion streams from disks without eccentric modes tend to lag behind the binary.The aforementioned effects conspire to create a much more negative gravitational torque on binaries accreting from disks; the difference in torque follows from a fundamental change in the accretion flow morphology, and the primaries actually accrete (proportionally) more from disks without eccentric modes than from disks with prominent eccentric modes.

Figure 12 .
Figure12.Average surface density (top panel) and cumulative gravitational torque (bottom panel, the gravitational torque density integrated from r ′ = 0 to r ′ = r, normalized by the accretion rate onto the binary) profiles for binaries accreting from M = 10 disks.Similar changes in both the surface density and cumulative torque profiles occur for both viscosities, the disk moving inwards and the overall torque becoming negative, although these changes Each profile was averaged over the final 500 binary orbits of each simulation.
Figure14.Accretion rate time series, over fifteen orbital periods, for a representative sample of binaries of mass ratios.In each panel, the accretion rate onto the primary is plotted in blue ( Ṁ1/ Ṁ0), and the accretion rate onto the secondary is plotted in orange ( Ṁ2/ Ṁ0), in both cases normalized by the equilibrium accretion rate for each simulation.At lower mass ratios (q ≲ 0.1), the accretion rate onto the primary is fairly constant, and well below the accretion rate onto the secondary.

4
Besides using cells near r ∼ a b about ∼ 12 times smaller than D'Orazio et al. (2016) and ∼ 6 times smaller than Duffell et al.

Figure 16 .
Figure16.The effect of sink rate on the degree of preferential accretion.At sink rates that are too slow, mass builds up regardless of the sink method used, leading to erroneous torques and preferential accretion measurements.Once sufficiently fast sinks are employed, torque-free (δ = 0) sinks converge, while standard (δ = 1) sinks, if sufficiently fast, can lead to spurious results of any magnitude.
, Disco (Duffell 2016; Dittmann & Ryan 2021) ACKNOWLEDGMENTS We wish to thank Adam Dempsey, Dan D'Orazio, Paul Duffell, and Cole Miller for discussions during the course of this project and feedback on a draft of this paper.AJD is also grateful for many useful discussions that took place at the Summer 2022 Center for Computational Astrophysics N-Body Workshop, especially with Jillian Bellovary.The authors acknowledge the University of Maryland supercomputing resources (http://hpcc.umd.edu) that were made available for conducting the research reported in this paper, and the ASTRA cluster administered by the Center for Theory and Computation within the University of Maryland Department of Astronomy.Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Colleges and Universities.AJD was supported in part by NASA