Eccentric Minidisks in Accreting Binaries

We show that gas disks around the components of an orbiting binary system (so-called minidisks) may be susceptible to a resonant instability that causes the minidisks to become significantly eccentric. Eccentricity is injected by, and also induces, regular impacts between the minidisks at roughly the orbital period of the binary. Such eccentric minidisks are seen in vertically integrated, two-dimensional simulations of a circular, equal-mass binary accreting from a circumbinary gas disk with a Γ-law equation of state. Minidisk eccentricity is suppressed by the use of an isothermal equation of state. However, the instability still operates and can be revealed in a minimal disk-binary simulation by removing the circumbinary disk and feeding the minidisks from the component positions. Minidisk eccentricity is also suppressed when the gravitational softening length is large (≳4% of the binary semimajor axis), suggesting that its absence could be an artifact of widely adopted numerical approximations; a follow-up study in three dimensions with well-resolved, geometrically thin minidisks (aspect ratios ≲0.02) may be needed to assess whether eccentric minidisks can occur in real astrophysical environments. If they can, the electromagnetic signature may be important for discriminating between binary and single black hole scenarios for quasiperiodic oscillations in active galactic nuclei; in turn, this might aid in targeted searches with pulsar timing arrays for individual supermassive black hole binary sources of low-frequency gravitational waves.


INTRODUCTION
There is strong evidence that some galaxies host supermassive black hole binaries at their centers (SMB-HBs; Begelman et al. 1980).These objects are powerful sources of low-frequency gravitational wave (GW) radiation, and their population has long been theorized to generate a stochastic GW background (Carr 1980;Rajagopal & Romani 1995).The recent detection of such a background with pulsar timing arrays (Agazie et al. 2023;Antoniadis et al. 2023a;Xu et al. 2023;Reardon et al. 2023a), and indications that the SMBHB population is indeed its likely source (Afzal et al. 2023; Anto-Corresponding author: John Ryan Westernacher-Schneider john.westernacher.schneider@gmail.comniadis et al. 2023b;Reardon et al. 2023b), creates a new imperative to identify individual SMBHB systems.
Galaxies that host SMBHBs are likely to be active, due to the presence of copious circumnuclear gas associated with the galaxy merger that gave rise to the black hole pair (Barnes & Hernquist 1992).Such binary AGNs may exhibit quasi-periodic oscillations (QPOs) connected in some way to the binary's orbital motion, and more than 200 binary AGN candidates have been proposed based on QPO detections in electromagnetic (EM) surveys (e.g.Graham et al. 2015;Charisi et al. 2016;Liu et al. 2019Liu et al. , 2020;;Chen et al. 2020;Peñil et al. 2022;Chen et al. 2022;Li et al. 2023).Those SMBHB candidates, and others yet to be discovered, could become joint GW-EM sources as the pulsar timing arrays continue to accumulate sensitivity.
The identification of an individual low-frequency GW source with a particular time-varying AGN will be chal- Σ from two of our high-and very high-resolution runs, on independent color scales.Left: run labeled VH3 with resolution ∆x = 0.00125 a (see Table 1).Right: the corresponding depletion experiment with resolution ∆x = 0.0025 a, where the circumbinary disk has been removed (run labeled H2).
lenging, firstly, because of theoretical uncertainties in the relationship between the binary's GW and EM temporal signatures, and secondly, because AGN periodicity can also signify processes that have nothing to do with a binary, such as jet precession or accretion disk instabilities (see e.g.Hawkins 2002, for a review).It means that many of the binary AGN candidates so far identified could be single accreting SMBHBs, exhibiting real periodicities.To separate the single black hole AGNs from the binaries will thus require a theoretical understanding of the unique temporal characteristics of the binary accretion process.
In a black hole binary accretion system, gas comes to the binary through a circumbinary disk (CBD), which feeds gas into minidisks around each of the black holes.The electromagnetic emission from the minidisks and CBD could be significantly affected by their dynamical interaction with the binary, leading to an imprint on many bands across the electromagnetic spectrum (e.g.Sesana et al. 2012;Farris et al. 2015;Panessa et al. 2019).
Variability in γ-rays is also expected in binary blazars, due to modulation of the rates of mass delivery to the black holes (Sobacchi et al. 2017;Peñil et al. 2020).Computational studies have revealed a variety of mechanisms that could lead to detectable, periodically varying EM output, associated specifically with the dynamics of binary accretion (e.g.Komossa 2006;MacFadyen & Milosavljević 2008;Noble et al. 2012;Roedig et al. 2012;Shi & Krolik 2015;D'Orazio et al. 2016;Muñoz et al. 2019;Duffell et al. 2020;Zrake et al. 2021;Franchini et al. 2023).These range from periodic variations at or near the orbital period, to those over tens of binary or-bits, and they originate from distinct spatial regions of the accretion system.
In a previous study (Westernacher-Schneider et al. 2022, henceforth W22), we reported multi-wavelength light curves of thermal emission from accreting black hole binaries, computed from vertically integrated 2dimensional viscous hydrodynamics simulations with a detailed treatment of the radiative cooling.In that paper, we identified a strong QPO in the disk thermal luminosity at roughly the binary orbital period.However, we did not give a detailed discussion of the physical origin of the QPO, nor whether the same one had been identified previously by other authors.In this paper we report on the physical mechanism of the QPO we found in W22, and present very high resolution simulations revealing that it arises from an instability in which the minidisks become eccentric and exchange mass at a regular interval -see Fig. 1 for sample snapshots.The instability is sensitive to the thermodynamic treatment, being generally suppressed when locally isothermal or targettemperature "β-cooling" approximations are used.
Aside from the need to understand the temporal structure of EM output from accreting binaries to identify SMBHBs among candidates, there is also a new imperative to understand CBD morphology as radio imaging (at especially high resolution by ALMA) has revealed dramatic examples of well-resolved substructures in proto-planetary disks including those around young stellar binaries (e.g.Andrews et al. 2018;Huang et al. 2020;Sierra et al. 2021).
Our paper is organized as follows.In §2 we review theoretical results on periodic light curves from accreting binaries.In §3 we describe the simulation setup and diagnostics used to quantify the minidisk morphology.We describe the numerical approach in §3.2.Our simulation results are reported in §4, including subsections on the irrelevance of infall from the CBD for driving minidisk eccentricity ( §4.1), the central role of the minidiskminidisk interaction ( §4.2), the periodic mass trade between minidisks and precession effects ( §4.3), dependence on various parameters and prescriptions ( §4.4- §4.7), numerical convergence of the minidisk eccentricity growth rate ( §4.8), and a summary of our key numerical findings ( §4.9).We provide discussions in §5 on the mechanism of the eccentric instability ( §5.1), a comparison to other eccentricity-driving mechanisms ( §5.2), the role of gravitational softening ( §5.3 & §5.4), guidance for future 3-dimensional studies ( §5.5), and the observable consequences of the minidisk mass trade ( §5.6).We conclude in §6, and our suite of simulations is tabulated in the appendix.

BACKGROUND
We briefly review some of the mechanisms known to cause periodic oscillations in the light curves of accreting binary systems.Several of these are related to the formation of an m = 1 over-density in the CBD, referred to here as the "lump" (e.g.MacFadyen & Milosavljević 2008). 1 The lump forms near the inner edge of the CBD at r ≳ 3a, and orbits the binary at roughly the Kepler frequency (for reference that is ∼ 5 binary orbits if the lump orbits at a radius of 3a; in general we denote the lump orbital frequency as f lump ).Its presence leads to a modulation in the rate of mass delivered from the CBD to the minidisks on the time scale 1/f lump , because gas orbits in the CBD are generally eccentric, and feeding to the binary is enhanced when the lump orbits through its closest approach to the binary.Variations in the rate of feeding to the minidisks are not necessarily transmitted to the black holes.Indeed, the time scale for mass to accrete through the minidisks is typically in the range of 10's of binary orbits.However, simulations by a number of authors (e.g.Farris et al. 2014;Zrake et al. 2021) indicate that in spite of possible buffering by the minidisks, the 5 − 10 orbit lump period can still be detectable at the ∼ 10% level in the time series of the accretion power.Also, gas injection to a minidisk involves the impacts of gas streams from the CBD, and the resulting disturbances may propagate to the black hole in a fraction of an or-bit, much faster than the viscous rate.Furthermore, in radiatively efficient environments, the EM light curves could reflect the stream-minidisk impacts even if the black holes themselves accrete steadily.QPOs associated with lump-induced variation of gas delivery from the CBD may thus be a detectable feature of binary AGN light curves.
For a circular, equal-mass binary, the presence of the lump can lead to a second kind of periodic oscillation, at the frequency 2(f bin − f lump ) (e.g.Noble et al. 2012).This "binary-lump beat frequency" is the frequency at which one black hole or the other overtakes the orbiting lump.Simulations indicate that mass delivery to the binary is modulated at this frequency when the CBD extends inwards far enough for tidal stripping of the CBD to operate at any orbital phase, and thus be enhanced any time a binary component passes the lump.When the low-density cavity around the binary is very large, gas is only tidally stripped from the near side of the eccentric CBD, and the enhanced feeding rate at the frequency 2(f bin − f lump ) is suppressed (see examples in Farris et al. 2015, and W22) (one could say the duty cycle of this mode is reduced to a fraction of the lump orbit around its closest approach to the binary).For reference, when the lump period is 5 binary orbits, the binary-lump beating operates at a frequency of 1.6 times the binary orbital frequency.
In W22, we computed light-curves of thermal disk emission from an equal-mass black hole binary, and reported a QPO operating at between 1 and 2 times the binary orbital frequency.The similarity in frequency made it easy to confuse this feature with the binarylump beating effect, however there was an important difference to suggest it had a distinct physical origin: the QPO from W22 showed sensitivity to the length scale parameter r soft used in the code to soften the gravitational potential.In particular, the W22 frequency seemed to approach the orbital frequency as r soft was decreased.The binary-lump beating involves a coupling between the outer edges of the minidisks and the inner edge of the CBD, so it should not be sensitive to how gravity is numerically modeled very near the black holes.
We have carried out a detailed analysis since the publication of W22, and confirmed that indeed the QPO we saw there had nothing to do with the lump, nor the CBD in any direct way.Instead, we found the effect arises due to the minidisks developing a significant eccentricity, and experiencing regular collisions with one another as a result.The minidisks have opposing eccentricity vectors, and the disks collide to produce an EM flare when the long ends of the disks strike one another.The minidisk eccentricity vectors undergo retro-grade precession, and the collisions occur at the beat frequency between the minidisk precession and the binary orbit.As we show below, the rate of the minidisk precession increases with r soft , and this accounts for the observation from W22 that the eccentric minidisk beat frequency approaches the orbital frequency as softening is decreased.In the subsequent sections we present a detailed characterization of the eccentric minidisk beating effect, and an investigation of the conditions that lead to the growth of minidisk eccentricity.

NUMERICAL SETUP
Following W22, we study a binary with total mass M = M 1 + M 2 = 8 × 10 6 M ⊙ and semi-major axis a ≃ 10 −3 pc, yielding an orbital period T bin ≃ 1 yr, but in this paper we consider only an equal mass binary (mass ratio q = 1) on a circular orbit (e = 0).The vertically integrated pressure and density are P and Σ, respectively, and ϵ is the specific internal energy.We use an adiabatic equation of state P = Σϵ(Γ − 1) with Γ = 5/3, neglecting radiation pressure.In order to obtain numerically tractable Mach numbers ≃ O(10) while neglecting radiation pressure, the accretion rate must often be extremely super-Eddington (see e.g.D 'Orazio et al. 2013;Westernacher-Schneider et al. 2022).In this work, gas Mach numbers are chosen in the initial conditions, and are subsequently determined self-consistently, and are typically in the range ∼ 7 − 25.
In addition to cooling models, we also compare our results with locally isothermal runs, where the prescribed sound speed profile corresponds to a uniform orbital Mach number of 10.This is consistent with the target ϵ 0 we use in our target-temperature cooling runs.
Disk initial conditions correspond to near-equilibrium configurations about a single gravitating object.The gas configuration is allowed to settle around the orbiting binary over several viscous times before analysis begins.This corresponds to ∼ 1000 − 3000 binary orbits, depending on the model -see Table 1.The constantα models exclusively use self-consistent radiative cooling, and their initial conditions are Σ ∝ r −3/5 and P ∝ r −3/2 .The constant-ν models initially have Σ = constant, corresponding to a spatially uniform accretion rate.The subset of constant-ν models with radiative cooling initially have P ∝ r −3/4 , corresponding to local balance of viscous heating and radiative cooling, and yielding a Mach number profile M ∝ r −1/8 .The subset of constant-ν models with target-temperature cooling instead initially have uniform M = 10, which is a popular Mach profile used in both isothermal and targettemperature cooled models, and the target-temperature cooling term drives towards this initial Mach number.
As in W22, black holes are represented by torque-free sinks (see Dempsey et al. 2020;Dittmann & Ryan 2021) with radius r sink and sink rate s.Our gravity model derives from a Plummer potential Φ P ∝ (r 2 +r 2 soft ) −1/2 .As recognized in the literature (e.g.Huré & Pierens 2009;Müller et al. 2012), although some type of softening is numerically necessary to regulate the divergence at a Newtonian point mass, in two-dimensional calculations it physically represents the vertical integration of the force of gravity when the disk has finite thickness (we discuss this in §5.4).In our run with r soft = 0, we use a purely Newtonian force ⃗ F N outside the sinks, and transition to a Plummer force ⃗ F P (softened using r sink ) inside the sinks using a functional form ⃗ F = θ ⃗ F P + (1 − θ) ⃗ F N , where θ = 1 − (r/r sink ) 2 2 for distances r < r sink from the black hole, θ = 0 otherwise.This regulates the singular behavior at the black hole location while achieving zero softening in the regions of interest (i.e.regions outside the sink).
Lastly, we perform a set of "decretion" runs, where the binary is initialized in near-vacuum, and the sinks are replaced by sources, non-zero only within distances r s from each point mass, given by where ⃗ U are the conserved variables and ⃗ U 0 are their target values, given by a rigid circular rotation at speed (1/2)GM/r soft , a uniform surface density 0.1M/a 2 , and a chosen uniform Mach number.This source term strongly drives ⃗ U to ⃗ U 0 inside the source.A circumbinary decretion disk is prevented from forming by allowing ejected material to flow off the grid.Our suite of simulations is summarized in Table 1.

Minidisk diagnostics
To quantify minidisks individually, we integrate hydrodynamic quantities over the spatial region within a distance of 0.5a from their host black hole.These diagnostics are the minidisk mass, and the center-of-mass (COM) vector measured relative to the host black hole's location.Visual inspection confirms that the COM vector points in the direction of the farthest edge of the minidisk.
Finite minidisk eccentricity is found to be indicated by persistent non-zero COM amplitude over tens of binary orbits.A comparison between integrating over distances of 0.5a and 0.25a from black holes is provided in Fig. 2, and indicates that trends are robust to the size of the integration region.Crucially, the coherence of minidisk eccentricity shows up as an orderly precession, and this is indicated by a steady linear trend in the COM phase over tens of binary orbits.A tell-tale sign of a lack of persistent eccentricity is a jagged COM phase over time.This usually indicates that the COM vector is reflecting smaller scale or more transient features in the minidisks, rather than the coherent lopsidedness characteristic of the eccentric minidisk instability.A prototypical case of steady, coherent minidisk eccentricity is shown in Fig. 3 (top and bottom panels).We use two diagnostics to characterize the relationship between minidisks: relative orientation and mass flux.Their relative orientation is quantified by their relative COM phases.A prototypical case of a steady relative orientation of π radians is shown in Fig. 3 (middle panel).The mass-trading between minidisks Ṁtrade (t) is quantified by the root-mean-square (RMS) flux of mass across a line of length a through the origin, orthogonally bisecting the black hole separation.

Solution scheme
We use the Sailfish code, which is same code we used in W22, and we refer the reader to that work for most numerical details (Sailfish is a GPU implementation of Mara3 which was used in Tiede et al. 2020, Zrake et al. 2021and Tiede et al. 2022).A summary of parameters for our suite of runs is provided in Table 1.The computational domain usually extends to 12 a; exceptions are the high-resolution runs labeled H1−H4 (which extend to 15 a), the very high-resolution zoom-in run VH3 (which is evolved for a short time and whose grid extends to 7.5 a), and the decretion runs D0−D1 (which . Time series of the COM vectors and minidisk masses in a run (H2 in Table 1) where the CBD is abruptly removed (vertical dashed line).Although the minidisks are steadily depleted (bottom panel), the retrograde precession (top panel), anti-alignment (top-middle), and COM amplitude (bottom-middle) eventually restore the same qualitative behavior as when the CBD was present (Fig. 3).
extend to 1.75 a because there is no circumbinary disk).Combi et al. 2022), we assess the effect such an obstruction by placing a third sink at the origin of the binary system of radius 0.05 a, in addition to the two orbiting sinks representing black holes.Note our Cartesian grid has no singular behavior at the origin.

RESULTS
In Fig. 1 we show snapshots of the surface density Σ (raised to the power of 1/2 to improve visual contrast) from two of our high-resolution runs, focusing on the minidisks.Both runs use self-consistent thermal cooling with a nominal Mach number of M ∼ 11 and α-viscosity Σ in the CBD depletion experiment from the moment the CBD is removed until about 55 binary orbits later (run H2 in Table 1).The structure of the minidisks suffers a transient disruption which reestablishes over time via mass trading (see Fig. 4).
with α = 0.1.The runs shown in the left and right panels of Fig. 1 are models VH3 and H2 respectively (see Table 1).VH3 is a zoomed-in version of model H3, with double resolution (∆x = 0.00125 a).The H2 and H3 models have the same parameters, except that the one on the right (model H2) has had the CBD removed, to demonstrate that the minidisks develop eccentricity even if they do not interact with gas infall from the environment.In both cases, the minidisk eccentricity is persistent, i.e. the images in Fig. 1 are a good representation of how the disks would look at a randomly selected time in a well-evolved simulation.Fig. 1 also reveals that the minidisks settle into a configuration with their apsides oriented 180 • away from each other.In the sections below we examine these effects in detail.

Infall from the CBD is not required to drive minidisk eccentricity
In order to determine whether the minidisk eccentricity is driven by gas infall from the CBD, we restarted a run from a well-developed state, with the CBD subtracted and replaced with a near-vacuum (model H2).After the restart, the minidisks continue to evolve, but can no longer acquire gas from the environment.The right panel of Fig. 1 shows that the eventual relaxed state of the minidisks is again eccentric, and still has the disk apsides anti-parallel to one another.Fig. 4 shows Figure 6.An animated version of this figure is available in the HTML version of the final article (also available with higher quality at https://youtu.be/om15kZRhC18).A single frame from the animation is shown here.The animation is a 3 minute and 20 second video showing the evolution of √ Σ for all 50 binary orbits in the decretion experiment (run D1 in Table 1).The minidisks form from the inside out, fill their Roche lobes, and begin a synchronized mass trade, resulting in eccentricity growth (see Fig. 7).
the time series of the minidisk diagnostics following the depletion.Without feeding from the CBD, the mass in the minidisks diminishes over time as they accrete into the sinks (4th panel).The minidisk eccentricity (indicated by COM amplitude, 3rd panel), opposing orientation (2nd panel), and precession (1st panel) undergo a short disruption at roughly 20 binary orbits after the CBD is depleted, but minidisk eccentricity then restores over the subsequent 40 binary orbits.The maintenance of eccentric minidisks in the absence of gas infall from the CBD indicates that eccentricity is being injected by interactions between the minidisks.A video showing this run is provided in Fig. 5.
We have also confirmed that the eccentric minidisks can be established if mass is supplied from the sinks in a "decretion" run (model D1).For this model, a circular equal-mass binary was initialized in near-vacuum, with mass being steadily added to the system from the sinks (rather than subtracted, see §3).Animations of the decretion run (see Fig. 6) illustrate how the eccentric minidisks are established.First, gas flows out from the particle positions and forms minidisks around the binary components.Then as the minidisks grow in size and overflow their Roche lobes, they "collide" and exchange mass across the inner Lagrange point.After the mass-trading event, the minidisks recede inside their Roche lobes, but develop a small amount of eccentric- . Time series of the COM vectors and minidisk masses in a run where the CBD is removed and mass is supplied from the component locations ("decretion" run, D1 in Table 1).The behavior is largely the same as the CBD removal experiment shown in Fig. 4. Note the top-middle panel is zoomed in on π, compared to Fig. 4. The relative orientation of the minidisks begins at π due to the mirror symmetry of the system.
ity.Subsequently, the minidisks collide preferentially at their "long end," leading to further eccentricity injection and then stronger collisions.Fig. 7 shows the time series of minidisk diagnostics in the decretion experiment, and indicates that the eccentric minidisks, retrograde precession, and anti-parallel orientations become fully established.Fig. 8 shows the results of one final experiment, in which the minidisks are subtracted but the CBD gas is retained (model H3).The results are similar to the decretion run: the minidisks refill, this time from gas infalling from the CBD, and over the course of about 30 orbits they settle into the characteristic eccentric, antialigned configuration.In the 3rd panel we also show the model VH3, which is a zoomed-in version of H3 with double resolution (∆x = 0.00125 a), and which shows that the minidisk eccentricity settles to the same level.A video of this run is provided in Fig. 9.
Figure 8. Similar to Fig. 3, but for the minidisk refilling experiment run with a smaller softening length of r soft = 0.01 a (run H3 in Table 1).The VH3 run with double the resolution of H3 is also shown in the middle-bottom panel.

Role of the minidisk-minidisk interaction
The visual impression given by animations of the decretion and minidisk refilling experiments (models D1 and H3) is that interaction between the minidisks, and the associated mass exchanges, mediate the eccentricity growth.To see how things would be changed without the minidisk-minidisk interaction, we performed a run (model H4) where one of the minidisks is replaced by a large absorber of radius 0.45 a.In this configuration, one minidisk refills from the circumbinary disk, and can lose mass to the companion absorber, but does not receive stream impacts from a companion minidisk.Time series of the minidisk diagnostics in Fig. 10 show that the COM amplitude of the lone minidisk (2nd panel) exhibits large oscillations around roughly 0.05a.A video of this run is provided in Fig. 11.In contrast, with no absorber present (model H3; Fig. 8), the COM amplitude is about 0.1a with relatively little variation.Also, the minidisks undergo a steady rate of retrograde precession in the "normal" run H3, and that precession is not seen when the absorber is present (Fig. 8  Σ over roughly 60 binary orbits in the minidisk refilling experiment (run H3 in Table 1).Mass is drawn from the CBD to reform the minidisks, which fill their Roche lobes and begin a synchronized mass trade, resulting in eccentricity growth (see Fig. 8).

top panel
vs. Fig. 10 top panel).These observations indicate that some eccentricity must be injected by gas infall to the minidisks, but not persistently enough to account for the minidisks observed in the "normal" run H3.The persistent eccentricity, seen in runs that include both minidisks, indicates the minidisk-minidisk coupling is a likely cause of the directionally coherent eccentricity injection.

Periodic mass trade and apsidal precession
The eccentric minidisks collide periodically and exchange mass.This effect is shown quantitatively in the top panel of Fig. 12, where we plot the RMS mass flux across the midline between the binary components (the midline rotates at the binary orbital frequency).The spikes in the RMS mass flux correspond to a rate of mass exchange that exceeds the average mass flow to the binary by factors of 10 − 20.The mass transferred per collision exceeds 20% of the disk mass when the instability is most aggressive (Fig. 2), so the mass-trading events are dynamically significant.The pulses are very regular, and the interval is on the order of the binary orbital period.
The frequency of the mass trading events is accurately predicted by the beat frequency f bin − f prec associated with the binary orbital frequency f bin and the apsidal precession frequency f prec of the minidisks.This is con- .Similar to Fig. 8, but with one black hole replaced with an absorber of radius 0.45a (run H4 in Table 1).Without a companion, the minidisk develops some eccentricity, but the magnitude is smaller and more variable in time.
firmed in Fig. 13, which shows periodograms of the RMS mass flux between minidisks for runs S0-S5 (same runs as the first row of Fig. 14).Apsidal precession of eccentric disks in binary systems is generally governed by a combination of pressure gradients, viscous stresses, and the tidal field of the companion (see e.g.Goodchild & Ogilvie 2006;Kley et al. 2008).The precession rate seen in our simulations is additionally found to be sensitive to the gravitational softening length, r soft .For example, the run shown in the left-most panel of the top row of Fig. 14 has a precession rate that is consistent with zero, and that run (model S0) also has a zero gravitational softening length.We also performed a decretion run with zero softening (model D0 in Table 1) where the initial condition is a low-density atmosphere and gas is injected from the component locations, and with a small viscosity α = 0.001.This experiment reveals that kinematic shear viscosity tends to drive retrograde disk precession; in contrast to model S0 (which has α = 0.1), in model D0, which has much lower viscosity, we found the minidisk precession becomes prograde with a period of ∼ 47 binary orbits.A video showing this run is provided in Fig. 15.When gravitational softening is zero and the viscosity is negligible, the precession rate can still be positive or negative, and is then determined by a competition between tidal interaction with the companion (which drives prograde precession) and pressure gradients, which generally drive retrograde precession provided that radial derivatives of pressure are not too positive; Goodchild & Ogilvie (2006).1).Mass is drawn from the CBD to reform the lone minidisk, which fills its Roche lobe.Without mass trading from a companion minidisk, it does not develop a persistent value of eccentricity (see Fig. 10).

Dependence on gravitational softening
To determine the necessary conditions for the instability to operate, we have systematically "switched off" different pieces of physics.A summary of the visual results is presented in the panels of Fig. 14.The top row shows a sequence of representative images (again, color showing Σ 1/2 ) from runs where the gravitational softening length r soft is increased from 0.0 to 0.05a.It is visually evident that the instability gets weaker with larger r soft , and becomes too weak to see when r soft ≳ 0.03a.This trend is corroborated in Fig. 2, where we plot the distance of the minidisk COM from the respective binary component (as described in 3.1) as a function of r soft .Although the minidisk eccentricity is not visually obvious for r soft ≳ 0.03a, precession is nonetheless quantifiable (see e.g. the high-resolution model H1 with r soft = 0.04 a in Fig. 3), and Fig. 13 shows the cadence of minidisk mass exchange is still well-predicted by f bin − f prec .

Dependence on the orbital Mach number
The second row of images in Fig. 14 shows representative minidisk morphologies for a range of nominal Mach numbers in the range 7 − 25, and significant and persistent minidisk eccentricity is seen for all cases.The top panel of Fig. 16 shows the minidisk COM diagnostic as a function of the Mach number, and confirms that there is no clear dependence of the minidisk eccentricity on the disk temperature in the range we have simulated here.

Dependence on gas viscosity and suppression by target temperature profiles
The third row of Fig. 14 shows how the minidisk morphology depends on the gas viscosity, and reveals that high enough viscosity, ν ≳ 10 −4 reliably suppresses the instability, resulting in roughly circular minidisks.This is also corroborated in the bottom panel of Fig. 16.The fourth row of Fig. 14 (other than the right-most image) shows that the instability is significantly suppressed by the use of target temperature profiles.The degree of suppression is not markedly affected by the rate of driving towards the target temperature profile (compare 2nd and 3rd panels).Crucially, use of the isothermal equation of state (4th row, 4th panel), which is widely used in studies of binary accretion, strongly suppresses minidisk eccentricity in circumbinary accretion.However, in §4.8 we show that suppression by the isothermal equation of state can be overcome in some cases.

Effect of a "hole" near the origin
Many studies of binary accretion use a grid code with cylindrical polar coordinates, and such geometries could induce anomolous flow patterns in the vicinity of the coordinate origin.Given that mass transferred between the minidisks generally passes through the origin, we found it germane to examine how a source of systematic numerical error, such as arising from a coordinate singularity or inner boundary, might affect how the instability behaves.We modeled the source of error using a "hole" placed at the origin (runs labeled HOLE and S2 in Table 1), which is included as a third sink term as described in §3.2.
The minidisk morphology when the hole is present is shown at the bottom right panel of Fig. 14.The time se- ries of the minidisk COM, shown in Fig. 17, shows that the hole diminishes the average minidisk eccentricity by about half, and also reveals a large amplitude, slow oscillation about the mean eccentricity.This suggests the instability could be mischaracterized in simulations that use a cylindrical polar coordinate grid, unless particular care is taken to avoid numerical errors near the origin.This experiment also yielded serendipitous insight into the dynamics of the eccentricity driving.The slow oscillation of the minidisk COM is seen in both disks, but these oscillations are 180 • out of phase with one another; one disk gets more eccentric while the other circularizes.The oscillation period for the run shown in Fig. 17 is roughly 9 orbits, which is also the minidisk apsidal precession period for that run.We now understand that the radializing minidisk is hogging the gas falling in from the CBD, and that the circularizing minidisk is relatively starved.The explanation for this may be as follows: (a) the CBD cavity is eccentric, (b) the minidisks are eccentric, and (c) one of the disks extends more in the direction of the near side of the cavity wall, thereby receiving more of the infalling gas.These conditions apply too when no hole is present, but then gas flows relatively unimpeded from the disk which catches more of the infall to the one which catches less, and both disks remain fed.By inhibiting the mass and momentum transfer between disks, the hole leads to the starvation of one disk at a time, and also allows that disk to circularize.This is a further indication that the minidisk-minidisk interaction ( §4.2) is instrumental in driving the instability.

Numerical convergence of eccentricity growth rates
We performed a resolution study of the exponential growth rate of the minidisk COM amplitude using decretion experiments.We used the decretion scenario in order to remove the complicating effects of the mass in- fall from the CBD.Indeed, this scenario most cleanly isolates the role of the minidisk-minidisk interaction in driving the instability.When the CBD is removed, the instability occurs even with an isothermal EOS. 2 .Since isothermal runs are significantly less computationally expensive than radiatively cooled runs, we used the 2 Note that we have not observed robust development of the instability in any runs that both use a target temperature profile and include the CBD.The CBD introduces an asymmetry in mass feeding to the minidisks by feeding them one at a time, and the degree of its influence on the minidisk-minidisk mass exchange is evidently sensitive to the cooling model.The runs shown in Fig. 18 have nearly perfect symmetry in the minidiskminidisk mass exchange, and therefore indicate that a high degree of such symmetry plays an important role in this eccentric instability.Why CBDs suppress minidisk eccentricity when using target temperature profiles, but not otherwise, is not yet clearsee the discussion at the end of §5.1 Figure 15.An animated version of this figure is available in the HTML version of the final article (also available with higher quality at https://youtu.be/rFEvJFRePDA).A single frame from the animation is shown here.The animation is a 3 minute video showing the evolution of √ Σ over 50 binary orbits in the low-viscosity (α = 0.001), zero-softening decretion experiment (run D0 in Table 1).The minidisks form from the inside out, fill their Roche lobes, and begin a synchronized mass trade, resulting in eccentricity growth, but in this case the minidisks precess slowly in a prograde sense.
isothermal EOS with no CBD to illustrate the numerical convergence of the growth rate over a wider range of resolutions.Fig. 18 shows the exponential growth of the minidisk COM amplitude, in runs where the grid resolution is 200, 400, 800, 1600, and 3200 zones per semimajor axis.All of these runs develop eccentric minidisks, but the lowest resolution runs displayed in each panel, 200 (400) zones per a for an isothermal (radiatively cooled Γ-law) equation of state, show a spuriously large growth rate and early saturation.The growth rate is consistently measured to be approximately 0.07f bin at each subsequent doubling of the grid resolution.Saturation occurs around a consistent value.Note that the series of convergence runs displayed in Fig. 18 is not enumerated in Table 1.

Summary of key numerical findings
The results of our numerical investigation strongly point to the minidisk-minidisk interaction as a necessary and sufficient condition for the growth of persistent minidisk eccentricity.Namely, in runs which only have minidisks and their unimpeded interaction (models H2 and D1, §4.1), minidisk eccentricity is seen to grow and saturate at a persistent level, showing that the minidisk interaction is sufficient to generate persistent eccentricity; whereas shutting off the minidisk-minidisk  1).Bottom: Time-averaged COM amplitude of minidisks versus kinematic viscosity ν, for runs using a Γ-law equation of state (runs V1 -V10 in Table 1).Error bars correspond to the standard deviation of the COM amplitude.
interaction (model H2, §4.2), or impeding their interaction (model HOLE, §4.7), prevents persistent minidisk eccentricity.The minidisk eccentricity growth rate numerically converges ( §4.8), showing it is a physical effect.The minidisk-minidisk interaction manifests as the regular trading of mass between the minidisks across the inner Lagrange point, at roughly the orbital period of the binary.Departure of the observed mass trade interval from the orbital period is due to apsidal precession of the minidisks ( §4.3).Precession can in general be prograde or retrograde, but when r soft ≳ 0.01a it is always retrograde and gets faster with increased r soft ( §4.4).This effect is consistent with known retrograde precession of ballistic particles in a softened gravitational potential, which we also checked with numerical integrations of eccentric particle orbits in softened potentials.
The instability likely exists in a formal sense regardless of how thermodynamics is modeled, however it seems to be suppressed in scenarios where a CBD is present, and a target temperature profile is used, as with β-cooling or the locally isothermal EOS.

Mechanism of the instability
Our numerical calculations in §4 clearly point to a mechanism for minidisk eccentricity growth: we propose that minidisk eccentricity is injected by regular impacts between the minidisks satisfying a certain resonant phase condition.In this section, we present a heuristic for runs with and without a "hole" at r < 0.05a (runs labeled HOLE and S2 in Table 1).
picture of this whereby the outer edge of a minidisk is pictured as an eccentric ring of test particles in Keplerian orbit around a binary component, as illustrated in Fig. 19.The particles in the ring are subject to an external forcing term ⃗ f e (ν, θ), which depends on the ring eccentricity e, and the orbital phases, ν and θ, of the ring particle and of the binary orbit respectively.The zero of the binary orbital phase is chosen so that θ = 0 means the eccentricity vector ⃗ e 1 of the BH1 minidisk points horizontally to the right.
The forcing term needs to capture the dynamical effects of head-on impacts between gas parcels in opposing minidisks.Since more mass is exchanged per impact as the eccentricity grows (Fig. 2), the forcing amplitude must increase with e. Impacts occur when θ ≃ 0, and for small eccentricities they mainly affect the particles at the long ends of the minidisks around ν ≃ π.The impact force is directed opposite the particle's velocity ⃗ v orb (ν), and is proportional in magnitude to its speed.A possible forcing term is then where δ is a Dirac delta function and the constant in Eqn. 2 is positive.The ring evolves "rigidly", in the sense that, only that part of the forcing which determines the total torque and power applied to the ring is included in the equation of motion.This also means the ring does not precess, although one could estimate the rotation rates of ⃗ e 1,2 by averaging the local rate of apsidal rotation over the particle phases ν; non-elliptical distortions obviously cannot be captured in the ring approximation.Integration of ⃗ r × ⃗ f e and ⃗ f e • ⃗ v orb over dν yields respectively the ring specific torque l and specific power Ė, and in turn, an expression for ė(t) via e = 1 + 2Eℓ 2 /(GM ) 2 .A detailed solution of the e(t) equation is not needed to appreciate that circular rings subject to the forcing term in Eqn. 2 are unstable to small-amplitude perturbations.When 0 < e ≪ 1, the ring particles near the far turnaround points (overlapping ellipses in Fig. 19) experience a weak retrograde impulse, corresponding to the minidisk-minidisk impact.Backwards forcing near apocenter drives an angular momentum deficit, i.e. it increases the particle eccentricity, and that effect is not compensated around the minidisk pericenter because of the factor δ(ν − π).The larger e leads to a stronger retrograde impulse via Eqn.2, completing a feedback loop in which e(t) grows exponentially.In §4.8 we determined the growth rate to be ≃ 0.07f bin ; this rate empirically fixes the constant in Eqn. 2.
A stochastic forcing term could be added as a model of gas falling in from the CBD and impacting the minid-sks.The result should be the appearance of a non-zero but random-walking minidisk eccentricity, like what we observed in the simulations from §4.2 and §4.7,where the minidisk-minidisk interaction was suppressed by use of a large absorber, or a "hole", respectively.
The mechanism proposed here for the eccentric minidisk instability does not deal with the hydrodynamical energy budget, nor the complicating effects of feeding from a CBD, and therefore cannot account for the instability's apparent sensitivity to the thermodynamical treatment when a CBD is present.Our numerical results are consistent with two possible interpretations.The first, is that the regularity of minidisk-minidisk impacts is compromised by the appearance of spiral arms, and that spiral arm formation is directly sensitive to the equation of state.The second, is that the equation of state directly affects the CBD morphology, which in turn sets the cadence and regularity of minidisk feeding, to which the eccentricity evolution is sensitive.In the second scenario, the absence of spiral arms (Fig. 14) could be a red herring, or it could be a consequence of the disks already being eccentric.This issue needs to be investigated further.

Comparison to other eccentricity mechanisms
The physical picture just proposed is adapted from one that was described in Lubow (1994, henceforth L94) to explain the growth of eccentric disks around the white dwarf accretors of SU Ursae Majoris (SU UMa) binary systems.Those systems are seen to exhibit so-called "superhump" oscillations during periods of enhanced mass transfer from the donor star.The oscillations are widely interpreted as signifying an eccentric disk around the white dwarf, which precesses and causes the observed superhump mode to occur at the beat frequency with the binary orbital period.Eccentricity is known to be excited by the 3:1 Lindblad resonance (e.g.Whitehurst & King 1991;Lubow 1991;Franchini & Martin 2019) operating in the outer edge of the disk, however L94 was exploring an alternative in which the eccentricity was driven instead by the gas stream from the donor star impacting the disk around the white dwarf.
The ballistic particle-ring approximation with external forcing was used in L94 to analyze the eccentricity injection by the impacting gas stream, however with a different forcing term from the one in Eqn. 2. In L94, the forcing strength was set proportional to the rate of mass flow from the donor star, which would not be in resonance with any waves excited in the disk.L94 showed that the ring eccentricity is excited during periods of increasing mass flow, but then is dissipated after the mass flow rate stabilizes to a new level.It was concluded in L94 that stream impacts were not a viable scenario for eccentricity injection in SU UMa systems, and that the 3:1 resonance was the more likely culprit.
It is relevant to note that we considered the 3:1 Lindblad resonance as a possible mechanism for the eccentric minidisk instability.However, that has been shown to succeed only when the binary mass ratio is q ≲ 1/3 (see e.g.Whitehurst & King 1991;Murray et al. 2000), whereas we see the eccentric minidisk instability operating when q = 1.Besides, the Lindblad resonance is a tidal interaction, and our results from §4.2 and §4.7 indicate the eccentric minidisk instability is being driven by resonant mass exchange.
In §4.2 we established that the resonant interaction could be destroyed by replacing one minidisk with a large absorber, but we also saw that some eccentricity was nonetheless developing in the extant minidisk, albeit without the coherent directionality.We interpreted this as arising from stochastic eccentricity injection by gas infall from the CBD, however we have considered the possibility that minidisks might also be susceptible to some kind of secular instability.For example, Kozai-Lidov oscillations can cause an initially circular minidisk to develop eccentricity (Martin et al. 2014;Franchini et al. 2019), but this mechanism does not apply in our case because it requires a large disk inclination with respect to the binary orbital plane.In another example, isolated α-disks were found by Lyubarskij et al. (1994) to be unstable to eccentricity growth by a viscous overstability.Later work by Ogilvie (2001) pointed out that viscous overstability could be an unphysical aspect of the α-disk model, because it is suppressed when accounting for the finite relaxation time of magnetohydrodynamic turbulence expected in accretion disks.Possible scenarios where viscous overstability may be physical were further elucidated in Latter & Ogilvie (2006).Ogilvie (2001) also showed that bulk viscosity suppresses viscous overstability, and this fact was used by Kley et al. (2008) to test whether viscous overstability was important in the development of eccentric disks.
A simple way to assess the importance of viscous overstability is to perform runs with non-zero kinematic bulk viscosity λ.Indeed, it was argued in Kley et al. (2008) that, if viscous overstability were important, then using λ = 2ν would suppress disk eccentricity.We checked this case (see the panel labeled "λ = 2ν = 10 −4 " in Fig. 14), but we found no significant suppression of eccentricity (minidisk COM amplitude is still ≃ 0.09 a).Ogilvie ( 2001) also derived a quenching condition for viscous overstability when α = 0.1, namely that the bulk α-viscosity parameter is > 0.35.Thus, we also checked a case with λ/ν > 0.35 (see the panel labeled . The geometry of the minidisk-minidisk impacts.Dotted lines depict the outer edges of the minidisks, with the first black hole (BH1) in blue and the second (BH2) in black.The coordinate system is chosen so the minidisk apsides are horizontal, and the binary has orbital phase θ = 0 when BH1's minidisk's eccentricity vector points to the right (top panel).The red circle depicts the binary orbit.The true anomaly ν of a ring particle orbiting BH1 is shown in the bottom panel.The minidisks do not precess here, i.e. their orientations are fixed in the inertial frame; the middle and bottom panels show the configuration at orbital phases θ = π/2 and θ = π respectively.The minidisks collide and eccentricity is injected to the particle rings around the orbital phase θ = 0 (top panel).
"λ = 4ν = 10 −4 " in Fig. 14), and we again found no suppression of minidisk eccentricity (minidisk COM amplitude is still ≃ 0.09 a).In both bulk viscosity tests, we reduced ν to below 10 −4 √ GM a to ensure the tests were not affected by the viscous suppression demonstrated in the bottom panel of Fig. 16.We conclude that viscous overstability is not a likely explanation for the appearance of eccentric minidisks.

Why gravitational softening produces less eccentric minidisks
We found (top row of Fig. 14) that minidisk eccentricity is suppressed by gravitational softening.This can be understood in terms of ballistic particle trajectories in softened gravitational potentials.Consider the effective potential for a gas parcel of specific angular momentum ℓ orbiting in the softened potential of an object with mass M , The turning points in this potential are fixed by the specific orbital energy E of the gas parcel.Orbital eccentricity is not defined in the usual sense when r soft > 0, however, if ℓ and E are both fixed, then the radial distance between the turning points can be easily seen to decrease with increasing r soft .The result is that a given forcing amplitude (i.e. a fixed value of the constant in Eqn. 2) results in a smaller geometrical distortion of the disk when r soft is larger.This effect could be accounted for by replacing e(t) in Eqn. 2 with a different function that reflects the degree to which a ring with parameters ℓ and E is non-circular.Doing so would predict a slower growth rate and could account for the observed reduction of minidisk eccentricity with larger r soft .

Softening-driven apsidal precession
In the vertically integrated thin disk setting, gravity is softened at second order in h/r, where h is the vertical disk height measured from the midplane.This can be seen by introducing an ansatz for the vertical density profile, say ρ = ρ 0 (1 − (z/h) 2 ) for |z| < h, ρ = 0 otherwise, where ρ 0 has no dependence on z.The amplitude of the horizontal component of the gravitational force density is ρ(GM/R 3 )r, where R 2 = z 2 + r 2 and r is the cylindrical radial coordinate.Taking advantage of the fact that z/r ≪ 1, we can write Integrating over z ∈ [−h, h] and defining Σ ≡ The factor in square brackets is ≤ 1, and therefore weakens ("softens") the gravitational force.Gravitational softening can thus be understood as modeling the finite thickness of the disk.If the factor in square brackets also decreases for decreasing r (a condition which depends upon h(r)), it will soften more at smaller r, similar to the Plummer potential we use to model gravity in our simulations.
In practice, a commonly used model for softened gravity in thin disks derives from the Plummer potential (Plummer 1911), where r soft is the softening length.Based on comparisons to three-dimensional disks, the softening length in the Plummer potential ought to be on the order of the disk scale height (Müller et al. 2012). 3 In this section, based on the Plummer model of gravitational softening, we describe conditions under which one might expect softening-driven retrograde apsidal precession of planar eccentric disks.To understand the role of softening, we consider a single gravitating mass, and neglect the disk self-gravity and hydrodynamic effects such as pressure gradients and effective viscosity (see a discussion of these other effects in Goodchild & Ogilvie 2006).We therefore operate in the ballistic approximation around a single gravitating mass, whereby fluid elements are treated as test masses moving freely under gravity.We take the softening length to be linear in h with a constant of proportionality that is of order unity (Müller et al. 2012).
Consider first the case of a razor thin disk, such that r soft ∝ h = 0.In this case, the gravitational force is Newtonian, thus we expect eccentric orbits to be closed ellipses (i.e.zero precession).The same expectation holds whenever the disk has a constant aspect ratio, h ∝ r, because then r soft ∝ r and the Plummer potential becomes proportional to the Newtonian potential.In this case, eccentric orbits are still closed ellipses, but it is as though the central gravitating object has a suppressed mass.This condition is representative of gas-dominated α-disks, as they have relatively constant aspect ratios (since h ∝ r 21/20 or h ∝ r 9/8 when electron or freefree scattering opacity dominate, respectively; see e.g.Haiman et al. 2009).
On the other hand, radiation-dominated disks have constant disk scale height (see e.g.Haiman et al. 2009), i.e. h/r ∼ 1/r, so that the Plummer force is approximately Newtonian for r ≫ h but weaker than Newtonian for r ≃ h.In this case, the deficit of (vertically integrated) gravity near the central object causes eccentric orbits to precess in the retrograde direction.This can be understood intuitively as follows.At large distances, the eccentric orbit is approximately Newtonian, i.e. an ellipse.But close to the central object, gravity becomes increasingly weaker than Newtonian, and unable to close the particle trajectory to an ellipse.This causes its next apocenter to be rotated in the direction opposite of the orbital motion.We verified this picture numerically by evolving test particle trajectories in a Plummer potential.

Eccentric precessing minidisks in 2D versus 3D
As guidance for three-dimensional studies, in this section we point to regions of parameter space where the effects of minidisk eccentricity and retrograde precession may reveal themselves.
Since minidisk eccentricity is triggered by masstrading activity between minidisks, it is important that such activity not be disrupted by, e.g.artificial obstructions between them.Thus, the entire region between the binary should be resolved.Since minidisk eccentricity is suppressed by viscosity, three-dimensional studies seeking to reveal eccentric minidisks should have weaker effective viscosity.The value of α = 0.01 achieved in Oyang et al. (2021), for example, should be amply low, since we find eccentric minidisks with viscosity as high as α = 0.1.Since minidisk eccentricity is suppressed by gravitational softening, and softening represents the finite thickness of disks, a three-dimensional investigation should strive to make the disk thinner.Simulating thinner disks in three dimensions increases computational cost due to the higher resolution required in the z-direction.Although this does not increase the number of cells required in the vertical direction, the cells are smaller, which tightens the time step constraint.If one strives to simulate the r soft = 0.01 a case (which yields obvious minidisk eccentricity in 2 dimensions, e.g. e ∼ 0.5), and assuming the softening length is ≃ 0.5 h where h is the disk half-thickness, then one requires that h does not exceed ≃ 0.02 a within the minidisks.Note that the insensitivity of minidisk eccentricity to Mach number in our study is not expected to be repro-duced in three-dimensional studies, because the effective softening length is intimately tied to Mach number via M ∼ (h/r) −1 ; whereas in our study, the softening length is instead an independent ad hoc parameter, artificially decoupled from Mach number.
On the other hand, testing softening-driven retrograde precession in three-dimensional studies requires disks that are sufficiently thin (such that minidisk eccentricity is appreciable), but still sufficiently thick that the effect of the implied gravitational softening on precession dominates over other hydrodynamical and tidal effects (see §4.3).Our results suggest that the r soft = 0.01 a case should be sufficient (i.e.h ≃ 0.02 a).However, the functional form of the Plummer potential also suggests that disks with nearly constant aspect ratios will not undergo retrograde precession (see §5.4); instead, three-dimensional studies seeking to reveal retrograde minidisk precession should focus on flatter disk profiles (such as constant disk heights expected in radiationdominated disks).Note that retrograde precession should not require a binary, so a targeted simulation of a disk with flat height profile and eccentricity initialized to e.g. e ≃ 0.5 around a single gravitating object ought to be a sufficient test of the effect.Cylindrical polar coordinates, rather than spherical, would be efficient for simulations of flat disks.
In their relativistic simulations, Avara et al. ( 2023) reported time-varying tilts of the minidisks out of the equatorial plane, with tilt angles comparable to the aspect ratio of the disk.If severe enough, such tilts could cause eccentric minidisks to miss each other at the phase θ = 0 shown in Fig. 19, thereby inhibiting eccentricity growth.Thus, it is conceivable that three-dimensional effects not captured in the vertically integrated approach could prevent minidisk eccentricity growth in some scenarios.
There are subtleties about our two-dimensional models which may be important to understand when comparing with three-dimensional simulations.Firstly, the Plummer potential is an ad hoc model of the gravitational softening that occurs when integrating out the vertical degree of freedom in thin disks.In particular, it is not derived in a controlled perturbative procedure in powers of the local disk aspect ratio.To do so would require greater knowledge of the local vertical density profile.Thus, a lack of softening-driven retrograde precession in constant aspect ratio disks is only predicted to the extent that the Plummer potential is a reasonable model of softened gravity in that regime (see Müller et al. 2012, for some validations of the Plummer potential against three-dimensional calculations).However, we expect that retrograde precession in flat disks (i.e.h = constant) should be a generic consequence of gravitational softening, independent of the applicability of the Plummer model.
Secondly, when performing a vertical integration of the hydrodynamic equations of a thin disk, the gravitational force softens beginning at second order in the disk aspect ratio.Instead, if a polar integration is performed, the magnitude of the gravitational force per unit area can be calculated at fully nonlinear order as GM Σ polar /R 2 since the coordinates conform with the spherical symmetry of the point mass gravitational potential.Here, we defined Σ polar ≡ π/2+θ h π/2−θ h ρRdθ to be the exact surface density, where θ h defines the disk's local polar extent measured from the midplane.In other words, with polar integration, gravity does not appear to have a softened functional form at fully nonlinear order.Thus, it is more cautious and nuanced to say that retrograde precession is possibly a finite-thickness effect, which manifests via softened gravity under vertical integration, but may have alternative physical interpretations in other two-dimensional reductions.Ambiguity in the physical interpretation of thin disks was recognized in Abramowicz et al. (1997), e.g. the use of spherical versus cylindrical coordinates trades between a vertical gravitational force and a vertical centrifugal force.On this note, it is worth pointing out that pressure and finite disk thickness (which are not mutually exclusive) are known to influence the apsidal precession rates of disks (see e.g.Kato 1983;Goodchild & Ogilvie 2006), with pressure effects in particular giving rise to retrograde forcing (as long as radial derivatives of pressure are not too positive) which becomes stronger for thicker disks (Goodchild & Ogilvie 2006).This increased retrograde driving with thicker disks (also reported in Kley et al. 2008) is at least directionally consistent with the softening-driven precession we describe in this work (i.e.thicker disks imply larger softening, which implies stronger retrograde driving).
Lastly, two-dimensional models of disks with explicit viscosity are, strictly speaking, turbulence closure models.Hence, the evolved variables must be understood as suitable averages (e.g.Favre-averaged quantities, see Favre 1969).A strict comparison with three-dimensional simulations would therefore necessitate computing such averaged quantities.If eccentric minidisks and softening-driven retrograde precession are manifested in three dimensions in an average sense, but not in any instantaneous sense, or if these effects are complete artifacts of the two-dimensional models, this would be a peculiar aspect of turbulence models of thin disks that theorists ought to be aware of.

Observable consequences
In our runs varying the softening length (labeled S0 -S5 in Table 1, and see the first row of Fig. 14), ∼ 2−13% of the minidisk mass is exchanged between minidisks per trading event (see the red curve in Fig. 2).Averaged over time, this corresponds to a mass exchange rate of ∼ 0.6 − 2.1 times the total mass accretion rate (see Fig. 12).Mass trading events can therefore be significant hydrodynamic events that cause observable EM flares.In our previous work W22, we reported simulated light curves 4 from accreting equal-mass circular binaries which exhibit periodic flaring at near-orbital frequency.In this work, we explain this periodicity as corresponding to mass trading events between minidisks.Flares are caused by shock-heating during minidisk-minidisk collisions, and we have checked that the flares are coincident 4 Note that the simulations presented here, as well as in W22 and in many similar previous studies, predict only the (mainly UV) thermal emission from the optically thick plasma, heated by viscosity and also shocks.At higher photon energies, additional non-thermal radiation is expected to be produced, similar to the hot coronal X-ray emission observed in AGN powered by solitary SMBHs (see e.g.Sesana et al. 2012).
in time with these collisions; we leave a more detailed study of each collision to future work.The frequency of mass trading is a beat frequency f bin − f prec between the binary orbit f bin and the signed minidisk precession f prec (negative values meaning retrograde precession), indicated with the dotted vertical line in the bottom panel of Fig. 20.The system's optical periodogram (obtained in the same way as W22) exhibits a peak at this beat frequency.The well-established m = 1 overdensity in the circumbinary disk (called the "lump") has a pattern speed which imprints on the optical emission at a frequency of f lump = 0.1 per binary orbit (see the leftmost peak in Fig. 20).We note that f bin − f prec is a distinct physical phenomenon from the beat frequency between the binary and the lump, 2(f bin − f lump ) (shown as the dashdotted vertical line in Fig. 20), which forms when all sides of the cavity wall are sufficiently close to the binary that the binary can strip material from it at almost all lump phases (e.g.cavities as reported in Noble et al. 2012;Bowen et al. 2019;Gutiérrez et al. 2021;Combi et al. 2022, which are smaller, less offset, and less eccentric than the run we show in Fig. 20).The top panel of Fig. 20 shows a snapshot of Σ (raised to the 1/3rd power to improve contrast), showing that the cavity is far too large and offset for the binary-lump beat frequency to form.
Also note that a similar minidisk mass-exchange phenomenon was observed in relativistic simulations (Bowen et al. 2017;Avara et al. 2023).That phenomenon was reported as being an effect enhanced to a significant level by relativistic gravity, characterized by a sloshing flow behavior resulting mostly in alternating mass transfer between minidisks, and associated with large morphological transitions of the minidisks (referred to as "disk-dominated" and "stream-dominated" states in Avara et al. 2023).This contrasts with our finding in many ways, since we find a strong phenomenon occurring even at the Newtonian level, characterized by eccentric, precessing minidisks, and sychronized mass trading events.Assuming the phenomenon we report in this work is indeed distinct from that reported in Bowen et al. (2017) and Avara et al. (2023), there is potential to confuse their observational signatures.The cadence of the sloshing effect between minidisks reported in the more recent work by Avara et al. (2023) is roughly 1.4× per orbit, quite similar to our f bin − f prec when the softening length is ∼ 0.04 a (see Fig. 13).
The rate of mass exchanged in the sloshing mechanism (as reported in Avara et al. 2023), is at the level of 0.1 × ṀBHs , whereas the eccentric minidisk instability can lead to a much larger mass exchange rate of ∼ 0.6 − 2.1 × ṀBHs (e.g.see our run S1 in Fig. 12).
If eccentric minidisks do form in accreting SMBHBs, we showed in W22 that they could produce a detectable QPO at or near the binary orbital period, likely in the UV or optical bands.This knowledge could aid in the identification of EM counterparts to future individualsource detections by the pulsar timing arrays, or assist in targeted searches by placing a prior on the binary orbital period given the QPO periodicity (e.g.Arzoumanian et al. 2020).

CONCLUSIONS & OUTLOOK
In this work, we showed that accreting, circular, equal mass binaries are prone to an instability which grows a significant eccentricity in the minidisks around the binary components.The mechanism originates in mass trading between the minidisks, which tends to synchronize and become periodic, driving up eccentricity, and causing the eccentric minidisks to maintain opposing orientations.This process is especially strong in the limit of small gravitational softening.Gas impacts from a circumbinary disk are neither necessary nor sufficient to explain this effect.
We investigated the dependence of minidisk eccentricity on model details.We found that many model choices, such as the use of artificial target temperature profiles (e.g. the use of β-cooling or locally isothermal equations of state), large gravitational softening (e.g.r = 0.05 a), large viscosity (e.g.ν = 10 −3 √ GM a), and a grid obstruction between the minidisks, all suppress minidisk eccentricity.This may partly explain why significant minidisk eccentricity in circular equal-mass binaries has not been previously reported in the literature.We found that minidisk eccentricity is robust to large bulk-to-shear viscosity ratios, which suggests this phenomenon is robust to a finite relaxation time of magnetohydrodynamic turbulence (Ogilvie 2001).
We also showed that eccentric minidisks tend to precess steadily in the retrograde direction when gravity is softened.In the limit of zero softening, the precession can in general be prograde, zero, or retrograde, depending on the balance of driving from hydrodynamic and tidal forces.The minidisks trade mass at a beat frequency, f bin − f prec , between the binary orbital frequency f bin and the minidisk precession frequency f prec ; note the minidisk precession frequency is negative when the precession is retrograde.This "eccentric minidisk beat frequency" imprints on light curves from thermal disk emission, as we reported in Westernacher-Schneider et al. (2022); in this work we clarified that the physical origin of such periodicity is the minidisk mass trade.Although the frequency can be similar, this effect is distinct from the "binary-lump beat frequency" 2(f bin −f lump ) formed between the lump and the binary, which occurs when the cavity wall is sufficiently close to the binary that the binary can draw material from the lump at most lump phases.
In a careful interpretation of the two-dimensional thin disk setting, we argued that softening-driven retrograde precession is a finite-thickness effect, even though a precise physical interpretation is not obvious.We believe that future three-dimensional simulations could observe the eccentric minidisk instability, but acknowledge that high resolution may be required due to the possible need for a rather thin disk, with h/r ≲ 0.02.Threedimensional simulations could also help clarify the physical meaning of gravitational softening commonly used in vertically integrated hydrodynamical settings.Seeing as we restricted this work to circular, equal-mass binaries, future work should also determine the range of binary parameters where this eccentric minidisk instability operates.
All simulations were performed on Clemson University's Palmetto cluster, and we gratefully acknowledge the Palmetto HPC support team.We acknowledge support from National Science Foundation grants AST-2006176 (to ZH) and AST-1715661 (to ZH and AM), NASA ATP grant 80NSSC22K0822 (to ZH and AM), and use of the software NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), SciPy (Jones et al. 2001-), CuPy (Okuta et al. 2017).We thank the KITP at UC Santa Barbara for their hospitality during the Bi-nary22 workshop, where some of the early work for this project was performed.We also acknowledge Julian Krolik, Mark Avara, Alessia Franchini, Matthew Bate, and Steve Lubow for valuable discussions at that workshop and since.KITP is supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

Figure 1 .
Figure1.Snapshots of √ Σ from two of our high-and very high-resolution runs, on independent color scales.Left: run labeled VH3 with resolution ∆x = 0.00125 a (see Table1).Right: the corresponding depletion experiment with resolution ∆x = 0.0025 a, where the circumbinary disk has been removed (run labeled H2).

Figure 5 .
Figure5.An animated version of this figure is available in the HTML version of the final article (also available with higher quality at https://youtu.be/9pltm6oOHhE).A single frame from the animation is shown here.The animation is a 3 minute and 39 second video showing the evolution of √ Σ in the CBD depletion experiment from the moment the CBD is removed until about 55 binary orbits later (run H2 in Table1).The structure of the minidisks suffers a transient disruption which reestablishes over time via mass trading (see Fig.4).

Figure 9 .
Figure9.An animated version of this figure is available in the HTML version of the final article (also available with higher quality at https://youtu.be/GCh7yW-QuY8).A single frame from the animation is shown here.The animation is a 3 minute and 46 second video showing the evolution of √ Σ over roughly 60 binary orbits in the minidisk refilling experiment (run H3 in Table1).Mass is drawn from the CBD to reform the minidisks, which fill their Roche lobes and begin a synchronized mass trade, resulting in eccentricity growth (see Fig.8).
Figure10.Similar to Fig.8, but with one black hole replaced with an absorber of radius 0.45a (run H4 in Table1).Without a companion, the minidisk develops some eccentricity, but the magnitude is smaller and more variable in time.

Figure 11 .
Figure11.An animated version of this figure is available in the HTML version of the final article (also available with higher quality at https://youtu.be/Th13XvxKsxA).A single frame from the animation is shown here.The animation is a 3 minute and 46 second video showing the evolution of √ Σ over roughly 60 binary orbits in the single minidisk refilling experiment (run H4 in Table1).Mass is drawn from the CBD to reform the lone minidisk, which fills its Roche lobe.Without mass trading from a companion minidisk, it does not develop a persistent value of eccentricity (see Fig.10).

Figure 12 .
Figure 12.Top: Time series of the inter-minidisk mass trading rate Ṁtrade .Bottom: Time series of the instantaneous mass accretion rate ṀBHs to the black holes.Both signals are normalized by the time-averaged mass accretion rate ⟨ ṀBHs⟩, and are computed from our run S1.

Figure 13 .
Figure 13.Normalized PSDs of the RMS mass flux between minidisks for different values of r soft .Data is from the fiducial runs with Γ-law EOS shown in Fig. 2. Vertical solid lines indicate the "eccentric minidisk beat frequency" f bin − fprec, computed from the minidisk COM phase evolution, accurately predicting the peak frequencies.

Figure 16 .
Figure16.Top: Time-averaged COM amplitude of minidisks versus Mach number M(r = a, t = 0), for runs using a Γ-law equation of state (runs M7 -M25 in Table1).Bottom: Time-averaged COM amplitude of minidisks versus kinematic viscosity ν, for runs using a Γ-law equation of state (runs V1 -V10 in Table1).Error bars correspond to the standard deviation of the COM amplitude.

Figure 17 .
Figure17.Time series of the COM amplitude of minidisks, for runs with and without a "hole" at r < 0.05a (runs labeled HOLE and S2 in Table1).

Figure 18 .
Figure18.Resolution study of the exponential growth phase of minidisk COM amplitudes for decretion, with the locally isothermal (top panel) and radiatively cooled Γ-law (bottom panel) equations of state.The fiducial resolution is ∆x = 0.0025 a.The higher resolution Γ-law runs were too expensive to evolve to saturation.
Figure20.Top: Images of the scaled surface density Σ 1/3 , showing the extent of the eccentric cavity forming in radiatively cooled Γ-law runs.Bottom: Normalized PSD of the optical light curve for the same system as in Fig.3.The peak around 0.1f bin is the lump frequency and the peak at 1.34f bin is the "eccentric minidisk beat frequency" f bin − fprec.The large cavity leads to suppression of the "binary-lump beat frequency" 2(f bin − f lump ) introduced in §2.