Eccentric Minidisks in Accreting Binaries

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 nonzero 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 Figure 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 telltale 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 Figure 3 (top and bottom panels).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 π rad is shown in Figure 3 (middle panel). The mass trading between minidisks ${\dot{M}}_{\mathrm{trade}}(t)$ is quantified by the rms flux of mass across a line of length a through the origin, orthogonally bisecting the black hole separation.

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 that was used in Tiede et al. 2020, 2022; Zrake et al. 2021). A summary of parameters for our suite of runs is provided in Table 1 in the Appendix. The computational domain usually extends to 12a; exceptions are the high-resolution runs labeled H1−H4 (which extend to 15a), the very high resolution zoom-in run VH3 (which is evolved for a short time and whose grid extends to 7.5a), and the decretion runs D0−D1 (which extend to 1.75a because there is no CBD). We use Courant–Friedrichs–Lewy numbers in the range of 0.01–0.1. Radiative cooling requires delicate numerical treatment (we follow the method of Ryan & MacFadyen 2017). Motivated by studies that have a coordinate singularity or inner boundary between the minidisks (e.g., d'Ascoli et al. 2018; Bowen et al. 2019; Combi et al. 2022), we assess the effect of such an obstruction by placing a third sink at the origin of the binary system of radius 0.05a, in addition to the two orbiting sinks representing black holes. Note that our Cartesian grid has no singular behavior at the origin.

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 they can no longer acquire gas from the environment. The right panel of Figure 1 shows that the eventual relaxed state of the minidisks is again eccentric and still has the disk apsides antiparallel to one another. Figure 4 shows 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 (fourth panel). The minidisk eccentricity (indicated by COM amplitude, third panel), opposing orientation (second panel), and precession (first 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 Figure 5.

Zoom In Zoom Out Reset image size

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 Section 3). Animations of the decretion run (see Figure 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 eccentricity. Subsequently, the minidisks collide preferentially at their "long end," leading to further eccentricity injection and then stronger collisions. Figure 7 shows the time series of minidisk diagnostics in the decretion experiment and indicates that the eccentric minidisks, retrograde precession, and antiparallel orientations become fully established.

Zoom In Zoom Out Reset image size

Figure 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, anti-aligned configuration. In the third panel we also show the model VH3, which is a zoomed-in version of H3 with double resolution (Δx = 0.00125a), and which shows that the minidisk eccentricity settles to the same level. A video of this run is provided in Figure 9.

Zoom In Zoom Out Reset image size

The visual impression given by animations of the decretion and minidisk refilling experiments (models D1 and H3) is that the 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.45a. In this configuration, one minidisk refills from the CBD and can lose mass to the companion absorber, but it does not receive stream impacts from a companion minidisk. Time series of the minidisk diagnostics in Figure 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 Figure 11. In contrast, with no absorber present (model H3; Figure 8), the COM amplitude is about 0.1a with relatively little variation. In addition, 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 (top panel of Figure 8 vs. top panel of Figure 10). 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 that the minidisk–minidisk coupling is a likely cause of the directionally coherent eccentricity injection.

Zoom In Zoom Out Reset image size

The eccentric minidisks collide periodically and exchange mass. This effect is shown quantitatively in the top panel of Figure 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 (Figure 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.

Zoom In Zoom Out Reset image size

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 confirmed in Figure 13, which shows periodograms of the rms mass flux between minidisks for runs S0–S5 (same runs as the first row of Figure 14).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 leftmost panel of the top row of Figure 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 in the Appendix) 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 that the minidisk precession becomes prograde with a period of ∼47 binary orbits. A video showing this run is provided in Figure 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).

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 Figure 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 Figure 2, where we plot the distance of the minidisk COM from the respective binary component (as described in Section 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.04a in Figure 3), and Figure 13 shows that the cadence of minidisk mass exchange is still well predicted by f_bin − f_prec.

The second row of images in Figure 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 Figure 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.

Zoom In Zoom Out Reset image size

The third row of Figure 14 shows how the minidisk morphology depends on the gas viscosity and reveals that high enough viscosity, ν ≳ 10⁻⁴, reliably suppresses the instability, resulting in roughly circular minidisks. This is also corroborated in the bottom panel of Figure 16. The fourth row of Figure 14 (other than the rightmost panel) 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 toward the target temperature profile (compare second and third panels). Crucially, use of the isothermal EOS (fourth row, fourth panel), which is widely used in studies of binary accretion, strongly suppresses minidisk eccentricity in circumbinary accretion. However, in Section 4.8 we show that suppression by the isothermal EOS can be overcome in some cases.

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 in the Appendix), which is included as a third sink term as described in Section 3.2.

The minidisk morphology when the hole is present is shown in the bottom right panel of Figure 14. The time series of the minidisk COM, shown in Figure 17, shows that the hole diminishes the average minidisk eccentricity by about half, and it also reveals a large-amplitude, slow oscillation about the mean eccentricity. This suggests that 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.

Zoom In Zoom Out Reset image size

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 Figure 17 is roughly nine 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 that catches more of the infall to the one that 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 it also allows that disk to circularize. This is a further indication that the minidisk–minidisk interaction (Section 4.2) is instrumental in driving the instability.

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 infall 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; ⁶ see Section 5.1. Since isothermal runs are significantly less computationally expensive than radiatively cooled runs, we used the isothermal EOS with no CBD to illustrate the numerical convergence of the growth rate over a wider range of resolutions. Figure 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) EOS, 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 Figure 18 is not enumerated in Table 1 in the Appendix.

Zoom In Zoom Out Reset image size

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 that only have minidisks and their unimpeded interaction (models H2 and D1; Section 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 interaction (model H2; Section 4.2), or impeding their interaction (model HOLE; Section 4.7), prevents persistent minidisk eccentricity. The minidisk eccentricity growth rate numerically converges (Section 4.8), showing that 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 (Section 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 (Section 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.

Our numerical calculations in Section 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 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 Figure 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 the binary orbit, respectively. The zero of the binary orbital phase is chosen so that θ = 0 means that the eccentricity vector e ₁ of the BH1 minidisk points horizontally to the right.

Zoom In Zoom Out Reset image size

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 (Figure 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

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{e}(\nu ,\theta )=-\mathrm{const}\times e(t)\times \ \delta (\theta )\ \delta (\nu -\pi )\ {{\boldsymbol{v}}}_{\mathrm{orb}}(\nu ),\end{eqnarray} \tag{ 2 }$$

where δ is a Dirac delta function and the constant in Equation (2) is positive. The ring evolves "rigidly," in the sense that only that part of the forcing that determines the total torque and power applied to the ring is included in the equation of motion. This also means that 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 ν; nonelliptical distortions obviously cannot be captured in the ring approximation. Integration of r × f _e and f _e · v _orb over d ν yields the ring specific torque $\dot{{\ell }}$ and specific power $\dot{{ \mathcal E }}$, respectively, and, in turn, an expression for $\dot{e}(t)$ via $e=\sqrt{1+2{ \mathcal E }{{\ell }}^{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 Equation (2) are unstable to small-amplitude perturbations. When 0 < e ≪ 1, the ring particles near the far turnaround points (overlapping ellipses in Figure 19) experience a weak retrograde impulse, corresponding to the minidisk–minidisk impact. Backward 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 Equation (2), completing a feedback loop in which e(t) grows exponentially. In Section 4.8 we determined the growth rate to be ≃0.07f_bin; this rate empirically fixes the constant in Equation (2).

A stochastic forcing term could be added as a model of gas falling in from the CBD and impacting the minidsks. The result should be the appearance of a nonzero but random-walking minidisk eccentricity, like what we observed in the simulations from Sections 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 EOS. The second is that the EOS 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 (Figure 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.

The physical picture just proposed is adapted from one that was described in Lubow (1994, hereafter 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., Lubow 1991; Whitehurst & King 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, but with a different forcing term from the one in Equation (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 Sections 4.2 and 4.7 indicate that the eccentric minidisk instability is being driven by resonant mass exchange.

In Section 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 nonzero 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⁻⁴" in Figure 14), but we found no significant suppression of eccentricity (minidisk COM amplitude is still ≃ 0.09a). 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 "λ = 4ν = 10⁻⁴" in Figure 14), and we again found no suppression of minidisk eccentricity (minidisk COM amplitude is still ≃0.09a). In both bulk viscosity tests, we reduced ν to below ${10}^{-4}\sqrt{{GMa}}$ to ensure that the tests were not affected by the viscous suppression demonstrated in the bottom panel of Figure 16. We conclude that viscous overstability is not a likely explanation for the appearance of eccentric minidisks.

We found (top row of Figure 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,

$$\begin{eqnarray*}&&{u}_{\mathrm{eff}}(r)=\displaystyle \frac{{{\ell }}^{2}}{2{r}^{2}}-\displaystyle \frac{{GM}}{\sqrt{{r}^{2}+{r}_{\mathrm{soft}}^{2}}}.\end{eqnarray*}$$

The turning points in this potential are fixed by the specific orbital energy ${ \mathcal E }$ of the gas parcel. Orbital eccentricity is not defined in the usual sense when r_soft > 0; however, if ℓ and ${ \mathcal 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 Equation (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 Equation (2) with a different function that reflects the degree to which a ring with parameters ℓ and ${ \mathcal E }$ is noncircular. Doing so would predict a slower growth rate and could account for the observed reduction of minidisk eccentricity with larger r_soft.

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, ρ = ρ₀(1 − (z/h)²) for ∣z∣ < h and ρ = 0 otherwise, where ρ₀ has no dependence on z. The amplitude of the horizontal component of the gravitational force density is ρ(GM/R³)r, where R² = z² + r² and r is the cylindrical radial coordinate. Taking advantage of the fact that z/r ≪ 1, we can write

$$\begin{eqnarray}&&\rho \displaystyle \frac{{GM}}{{R}^{3}}r=\rho \displaystyle \frac{{GM}}{{r}^{2}}\left[1-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{z}{r}\right)}^{2}\right]+{ \mathcal O }{\left(\displaystyle \frac{z}{r}\right)}^{4}.\end{eqnarray} \tag{ 3 }$$

Integrating over z ∈ [ − h, h] and defining ${\rm{\Sigma }}\equiv {\int }_{-h}^{h}\rho \,{dz}$ yields

$$\begin{eqnarray}&&{\int }_{-h}^{h}{dz}\,\rho \displaystyle \frac{{GM}}{{R}^{3}}r\simeq {\rm{\Sigma }}\displaystyle \frac{{GM}}{{r}^{2}}\left[1-\displaystyle \frac{3}{10}{\left(\displaystyle \frac{h}{r}\right)}^{2}\right].\end{eqnarray} \tag{ 4 }$$

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 that depends on 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),

$$\begin{eqnarray}&&{\rm{\Phi }}=-\displaystyle \frac{{GM}}{\sqrt{{r}^{2}+{r}_{\mathrm{soft}}^{2}}},\end{eqnarray} \tag{ 5 }$$

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

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 free–free scattering opacity dominates, 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 to the orbital motion. We verified this picture numerically by evolving test particle trajectories in a Plummer potential.

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 mass trading 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 owing 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.01a case (which yields obvious minidisk eccentricity in two dimensions, e.g., e ∼ 0.5), and assuming that the softening length is ≃0.5h, where h is the disk half-thickness, then one requires that h does not exceed ≃0.02a within the minidisks. Note that the insensitivity of minidisk eccentricity to Mach number in our study is not expected to be reproduced in three-dimensional studies, because the effective softening length is intimately tied to Mach number via ${ \mathcal M }\sim {(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 Section 4.3). Our results suggest that the r_soft = 0.01a case should be sufficient (i.e., h ≃ 0.02a). However, the functional form of the Plummer potential also suggests that disks with nearly constant aspect ratios will not undergo retrograde precession (see Section 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 radiation-dominated 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 Figure 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 that may be important to understand when comparing with three-dimensional simulations. First, 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.

Second, 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² since the coordinates conform with the spherical symmetry of the point mass gravitational potential. Here we define ${{\rm{\Sigma }}}_{\mathrm{polar}}\equiv {\int }_{\pi /2-{\theta }_{h}}^{\pi /2+{\theta }_{h}}\rho {Rd}\theta$ 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.

In our runs varying the softening length (labeled S0–S5 in Table 1 in the Appendix; see also the first row of Figure 14), ∼2%–13% of the minidisk mass is exchanged between minidisks per trading event (see the red curve in Figure 2). Averaged over time, this corresponds to a mass exchange rate of ∼0.6–2.1 times the total mass accretion rate (see Figure 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 ⁸ from accreting equal-mass circular binaries that 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 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 Figure 20. The system's optical periodogram (obtained in the same way as in W22) exhibits a peak at this beat frequency.

Zoom In Zoom Out Reset image size

The well-established m = 1 overdensity in the CBD (called the "lump") has a pattern speed that imprints on the optical emission at a frequency of f_lump = 0.1 per binary orbit (see the leftmost peak in Figure 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 dashed–dotted vertical line in Figure 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; Combi et al. 2022; Gutiérrez et al. 2022, which are smaller, less offset, and less eccentric than the run we show in Figure 20). The top panel of Figure 20 shows a snapshot of Σ (raised to the 1/3 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 that 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.04a (see Figure 13).

The rate of mass exchanged in the sloshing mechanism (as reported in Avara et al. 2023) is at the level of $0.1\times {\dot{M}}_{\mathrm{BHs}}$, whereas the eccentric minidisk instability can lead to a much larger mass exchange rate of $\sim 0.6\mbox{--}2.1\times {\dot{M}}_{\mathrm{BHs}}$ (see, e.g., our run S1 in Figure 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 individual-source 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).