Circumbinary Accretion: From Binary Stars to Massive Binary Black Holes

INTRODUCTION

Circumbinary disk (CBD) accretion plays an important role in the evolution of many types of astrophysical systems, ranging from young binary stars, main-sequence and post-main-sequence binaries to supermassive binary black holes. Figure 1 illustrates the basic features of such accretion: Gas from large distances gradually spirals toward the binary in a CBD driven by viscous dissipation; the disk is truncated at a few binary separations by the varying gravitational force from the binary (the “egg beater”), forming a cavity; the gas at the inner edge falls inward through accretion streams towards individual stars (or black holes), forming circum-single disks (CSDs, or “mini-disks”) and eventually accreting onto each stars.

Circumbinary accretion has long been suggested to exist around binary massive black holes (MBHs) following galaxy mergers (e.g., [22, 202, 87, 205, 77, 53, 44]). This has been demonstrated in many numerical simulations over the years. An example can be found in Mayer et al. [195]: two galaxies (each containing a MBH), initially separated at $\sim100$ kpc, collide with each other, and eventually end up with two MBHs separated by $\sim10$ pc and surrounded by an extended ($\sim100$ pc) disk/torus at the center of the merged galaxy.

One of the key questions concerning circumbinary accretion is: Does the binary lose or gain angular momentum and how does the binary orbit evolve? The first discussion of this issue appeared in Begelman et al. [22]:

… infall of gas onto the binary can also lead to some orbital evolution. Gas may be flung out of the system, acquiring energy (and angular momentum) at the expense of the binary; alternatively, gas may accrete onto the larger hole, causing orbital contraction as the product $Mr$ is adiabatically invariant. In either case, the evolution time scale is

As we discuss in this review (see § 3), this issue has been controversial and the prevailing view has been challenged by recent studies. The evolution of SMBH binaries undergoing gas accretion may directly impact the low-frequency gravitational wave signatures probed by space interferometers such as LISA (Amaro-Seoane et al. [2]) and Pulsar Timing Arrays (Burke-Spolaor et al. 2019) (see § 4.1).

Circumbinary accretion disks are a natural byproduct of binary star formation via disk fragmentation (e.g., [34, 142, 228]). A number of these disks have been observed around Class I/II young stellar binaries – well-known examples include GG Tau, DQ Tau, and UZ Tau E (e.g., [83, 191, 242]) – and recently even around much younger Class 0 objects like L1448 IRS3B (Tobin et al. [294]) and IRAS 16293-2422 A (Maureira et al. [194]). With ALMA, many CBDs have been discovered through direct imaging (e.g., Czekala et al. [57]). Recent observations have revealed that the disk can be highly misaligned with the central binary (Kennedy et al. [131]) – an issue we will address in § 5. The observed properties of “mature” stellar binaries (such as the mass ratio distribution) may have been shaped by circumbinary accretion at the earlier (proto-stellar) phase (see § 4.2).

Circumbinary disks have been also found around many post-AGN binaries (e.g., Van Winckel [302]). These ”second-generation” disks, likely formed as a result of binary interactions during the AGB phase, can impact the evolution of the systems (see § 4.4).

Starting from NASA’s Kepler mission, planets have been found around stellar binaries using the transit method (e.g., Doyle et al. [78]). So far, more than a dozen of such systems are known (see Welsh & Orosz [310], Kostov et al. [140], and references therein). An interesting features of these circumbinary systems is that many of the planets are found very close (within a factor of 1.5) to the stability limit, i.e., if the semi-major axis of the planet is a bit smaller (e.g. by 10% in the case of Kepler-16b), the planet would be ejected from the system. These planets are unlikely to have formed in-situ, but must have migrated from far away in the protoplanetary disk. The dynamics of CBDs can strongly influence the formation and migration of the planets (see § 4.3).

Recently, a special type of circumbinary accretion has gained interest, in connection with the gravitational wave sources detected by LIGO/VIRGO. It has been suggested that merging stellar black-hole binaries can be produced in AGN disks (e.g., [282, 197, 283] (Bart 2017)). The hydrodynamical flows generated by binaries embedded in AGN disks have several distinct features compared to normal circumbinary accretion disks (see § 6).

THEORY AND SIMULATION OF CIRCUMBINARY DISK ACCRETION: BASIC CONCEPTS AND OVERVIEW

We first review some of the key theoretical concepts related to circumbinary accretion. The analytic aspect of binary-disk interaction through gravitational forcing is relatively straightforward, at least in the linear regime (i.e., when the perturbation of the disk by the binary is weak). The gravitational potential produced by the binary (with total mass $M_{\mathrm{b}}$, semi-major axis $a_{\mathrm{b}}$ and eccentricity $e_{\mathrm{b}}$) on the disk (assumed aligned with the binary plane) at the position $\mathbf{r}=(r,\phi)$ (measured from the binary’s center of mass) can be written as Goldreich & Tremaine [101]

where $\Omega_{\mathrm{b}}=(GM_{\mathrm{b}}/a_{\mathrm{b}}^3)^{1/2}$ is the mean angular frequency (or “mean motion”) of the binary and $\kappa_{\mathrm{b}}$ is the radial epicyclic frequency. For binaries in Keplerian orbits, $\kappa_{\mathrm{b}}=\Omega_{\mathrm{b}}$. The potential component $\Phi_{mn}$ depends on $e_{\mathrm{b}}$ and $a_{\mathrm{b}}/r$; to the leading order in $e_{\mathrm{b}}$, we have $\Phi_{mn}\sim e_{\mathrm{b}}^{|n|}\Phi_{mm}$. For $m>0$, the $(mn)$-potential rotates with the pattern frequency

(where $N=m+n$ and the second equality assumes $\kappa_{\mathrm{b}}=\Omega_{\mathrm{b}}$), and excites spiral density waves at the Lindblad resonances (LRs), where

with the upper (lower) sign corresponding to the outer (inner) LR. As the CBD is approximately Keplerian, $\kappa\simeq\Omega\simeq(GM_{\mathrm{b}}/r^3)^{1/2}$, the LRs are located at

The torque on the disk at a LR is Goldreich & Tremaine [100]

where $\Sigma$ is the disk surface density, $D=\kappa^2-m^2(\Omega-\omega_{mn})^2$, and

At the outer LRs (which are most relevant for CBDs), $(\mathrm{d}D/\mathrm{d}\ln r)=-3N^2\Omega_{\mathrm{b}}^2/(m+1)$, we find $T_{mn}^{\mathrm{OLR}}>0$, i.e., the disk particles gain angular momentum from the binary through resonant gravitational torques. This is a general result: A rotating potential always transports angular momentum from higher to lower angular velocity [181, 99]. The associated energy transfer rate to the disk through the LR is given by

These expressions are useful for determining how gravitational binary-disk interaction affects the binary orbit (see § 4.4).

Ωb: $(GM_{\mathrm{b}}/a_{\mathrm{b}}^3)^{1/2}$, mean angular frequency (mean motion) of the binary; $P_{\mathrm{b}}=2\pi/\Omega_{\mathrm{b}}$ is the binary orbital period

Ω(r): $(GM_{\mathrm{b}}/r^3)^{1/2}$, angular frequency of the CBD

LR: Lindblad resonance; OLR: outer LR

On the other hand, the disk particles lose angular momentum through viscous torque. Assuming the $\alpha$-ansatz for the kinematic viscosity coefficient, $\nu= \alpha c_{\mathrm{s}}^2/\Omega$ (where $c_{\mathrm{s}}$ is the disk sound speed), the viscous torque is given by [e.g. 250]

where $h = H/r$ is the disk aspect ratio. A gap is opened at the $(mn)$-LR if $T_{mn} \geq T_{\nu}$ [7]. In this picture, the radius $r_{\mathrm{cav}}$ of the inner cavity of a CBD is determined by largest radius at which a gap can be cleared.

As an example, for a circular binary, the dominant potential has $m = 2$, $n = 0$, with $\Phi_{mn} \simeq-3G\mu_{\mathrm{b}}a_{\mathrm{b}}^2/(4r^3)$ (assuming $r \gg a_{\mathrm{b}}$, where $\mu_{\mathrm{b}}$ is the reduced mass of the binary). The OLR is located at $r_{\mathrm{LR}} \simeq(3/2)^{2/3}a_{\mathrm{b}}$, and the LR torque is $T_{20}^{\mathrm{LR}} \simeq(49\pi^2/2)\Sigma(\Phi_{20}/\Omega_{\mathrm{b}})^2$. Gap opening at the OLR requires $T_{20}^{\mathrm{LR}} \gtrsim T_{\nu}$, i.e.

a condition easily satisfied for binaries with comparable component masses. For a small but finite $e_{\mathrm{b}}$, gap is cleared by the $m = 2$, $n = -1$ potential (with $\Phi_{mn} \propto e_{\mathrm{b}}$), and the corresponding OLR is at $r_{\mathrm{LR}} \simeq3^{2/3}a_{\mathrm{b}}$.

Miranda & Lai [207] considered eccentric binaries with general mass ratios and binary-disk inclination angles. For typical disk $c_{\mathrm{s}}$ and viscosity parameter, the inner radius $r_{\mathrm{cav}}$ of the CBD is found to be $(2$-$3)a_{\mathrm{b}}$ and depends on $e_{\mathrm{b}}$ in a discrete manner. Misaligned disks generally have smaller inner radii than aligned disks. In any case, it is important to recognize that such theoretical calculation of the inner disk radius has obvious limitations, as real disks are expected to have a fuzzy and dynamical inner boundary, with gas streaming into the cavity. Numerical simulations are needed to capture the whole complexities of binary-disk interactions, especially near the inner truncation radius.

More subtle resonant binary-disk interactions can also play a role in the dynamics of CBDs, such as the parametric instability associated with Lindblad resonances (see § 3.3) that can excite disk eccentricity [113, 172 b]. In addition, for an eccentric binary, secular (orbital-averaged) interaction may affect the disk eccentricity evolution [208, 174 see Section 3.3].

Because of the importance of circumbinary accretion in various astrophysical contexts, many numerical simulations have been carried out over the years. Some works were based on Smoothed Particle Hydrodynamics (SPH) [e.g. 20, 8, 87, 53, 262, 239, 258]. Others used grid-based Eulerian hydrodynamical methods: these include, among others, Günther & Kley [107] (hybrid grid), MacFadyen & Milosavljević [182] (FLASH, polar grid with inner cavity excised), Hanawa et al. [111] (nested Cartesian grid), de Val-Borro et al. [59] (cartesian grid), D’Orazio et al. [75] (FLASH, polar grid with inner cavity excised), Lines et al. [161] (FARGO, polar grid with inner cavity excised), Miranda et al. [208] (PLUTO, polar grid with inner cavity excised). Most simulations solved hydrodynamical equations with parameterized viscosity, but some short-duration (essentially) Newtonian MHD simulations have also been carried out [226, 271, 270, 37] and even GRMHD simulations, in which a time-dependent, sewn-together metric [193 e.g.,] can be prescribed in order to to take into account the Kerr metric in the vicinity of each black hole [50, 51].

Numerical simulations of circumbinary accretion are challenging because of the wide spatial range involved, and the multiples timescales on which variability takes place (Figure 2). The accreting gas flows from a large, viscously-evolving disk surrounding the binary, transitions into plunging accretion streams, which then feed the circum-single disks around each individual binary component, which accretes mass at spatial scales typically much smaller than the binary separation. In addition, since the flow in the disk and near the binary is highly dynamical, to determine the long-term effect of the flow on the binary, sufficiently long simulations with careful averaging must be carried out.

In recent years, several finite-volume moving mesh codes have been used to study circumbinary accretion. The first is DISCO that utilizes a moving ring grid [90, 74, 286] Duffell et al. 2020). Another is AREPO ([217, 220, 218]). The general-purpose grid-based Godunov code ATHENA++ has also been used for long-term simulations of CBDs ([213, 307 b]).

SIMULATIONS OF CIRCUMBINARY ACCRETION: KEY RESULTS

In this section, we summarize the key findings from recent simulations. Our discussion will be guided by our own simulation results using AREPO ([217, 220, 218, 276]) and related work using PLUTO ([208]), but we will compare with the results from other simulations when relevant, and discuss more recent progress. AREPO ([280, 234]) is a quasi-Lagrangian Godunov-type moving-mesh code. It has unstructured moving grid with adaptive resolution, and with hydrodynamical equations solved in the moving frame. AREPO has been adapted for viscous accretion disk simulations by Muñoz et al. ((2014, 2015)). In our works on circumbinary accretion, we simulate 2D Newtonian viscous flow, with a locally isothermal equation of state $P = \Sigma c_{\mathrm{s}}^{2}$, where the sound speed $c_{\mathrm{s}}(r)$ is a prescribed function of $r$ – we assume the disk has a constant disk aspect ratio $h \equiv H/r$. We use the Shakura–Sunyaev $\alpha$ prescription for the viscosity $\nu(r)$. Our typical simulations resolve accretion onto individual binary components down to $0.02a_{\mathrm{b}}$.

In the following, we first consider binaries with mass ratio $q_{\mathrm{b}} \sim1$ and disks with $H/r \sim0.1$ and viscosity parameter $\alpha= 0.05 - 0.1$. We examine both “infinite” disks (with a fixed mass supply rate at $r_{\mathrm{out}} \gg a_{\mathrm{b}}$ as well as finite disks. Later in this section, we discuss how the results vary for different parameters as well as various complications.

Short-Term Accretion Variability

For an extended disk, when the gas supply rate $\dot{M}_0$ at large distances is constant, the whole binary-disk system can reach a quasi-steady state, in the sense that the time-averaged accretion rate across the CBD and total mass accretion rate onto the binary, $\langle\dot{M}_{\mathrm{b}}\rangle= \langle\dot{M}_1\rangle+ \langle\dot{M}_2\rangle$, are equal to $\dot{M}_0$. However, the flow rates onto the inner binary cavity and onto individual binary components are highly variable on the binary orbital timescale [217, 208]:

• For nearly circular binaries (with $e_b \lesssim0.05$, the accretion rates $\dot{M}_1$ and $\dot{M}_2$ vary with a dominant period about $5P_b$ (where $P_b$ is the binary orbital period), a result that was already found in earlier simulations (e.g., [182, 75, 271]). This dominant period corresponds to the Kepler period of gas at the inner edge of the disk, $r_{\rm cav} \simeq3a_b$. The reason for this is that for circular binaries, the accretion onto the inner cavity arises from the development of $m = 1$ lumps at the inner disk edge and the lump has a pattern speed equal to the local Kepler velocity (see Figure 3).

• For $e_b \gtrsim0.05$, the dominant variability of $\dot{M}_1$ and $\dot{M}_2$ has a period of $P_b$. The reason is that for eccentric binaries, the mass transfer onto the cavity mainly occurs when the binary is at apocenter,

where the binary component is closest to the disk inner edge and “grabs” the gas from the disk and funnels it inwards (see Figure 4).

The different variability timescales ($P_{\mathrm{b}}$ vs $5P_{\mathrm{b}}$) are important when trying to infer the orbital period of binary massive black holes from the observed variability of dual AGNs (see § 4). These short-term variabilities can be also used to confront observations of T Tauri binaries (e.g. [295, 297]). In some cases, the effective sizes (which can be different from the stellar radii because of the magnetic fields) of the individual accretors can be constrained.

Longer-term Accretion Variability

Binary accretion can vary on a timescale much longer the binary orbital period (Figure 2), even when the mass supply rate $\dot{M}_0$ is constant. For an equal-mass circular binary, the accretion rates onto individual stars are quite similar to each other, following the same variable pattern in time, as illustrated in the left panels of Figure 5. By contrast, for eccentric binaries, one of the binary components can accrete at a rate 10-20 times larger than its companion, even for mass ratio $q_{\rm b}=1$, as illustrated in the right panels of Figure 5.

This “symmetry breaking” between the binary components, however, does not persist. Instead, it alternates over timescales of order $200P_{\rm b}$ (Muñoz & Lai [217]; see also Dunhill et al. [82]). This long-term quasi-periodicity can be attributed to a slowly precessing, eccentric CBD. Indeed, the inner region of the disk generally develops eccentricity and precesses coherently (see below). The apsidal precession rate of an eccentric ring (of radius $r$) around the binary is given by

At $r\sim3a_{\rm b}$ (the cavity radius), this corresponds to a period of 200-300$P_{\rm b}$. Over longer timescales ($\gg200P_{\rm b}$), the net accretion rates onto individual binary components are the same (for $q_{\rm b}=1$).

This “symmetry breaking” behavior (even for $q_{\rm b}=1$) is interesting in connection to AGNs around massive black hole binaries or similar accretion phenomena around binary T Tauri stars. It implies that even when the two binary components are very similar in mass, accretion may predominantly occur in one of components and last for hundreds of orbits. Interestingly, Tofflemire et al. [297] has reported evidence for preferential accretion onto the primary of the T Tauri Binary TWA 3A, in contradiction with expectations. that accretion should be preferential onto the secondary (see Section 3.6 below). Having an eccentricity of $e_{\mathrm{b}} = 0.67$ (Tofflemire et al. [296]), TWA 3A could be an example of this intriguing behavior.

Disk Eccentricity and Precession

Numerical simulations have long shown that eccentricity can develop in the innermost region of the CBD (MacFadyen & Milosavljević [182]; Miranda et al. [208]). While this is most noticeable by the formation of a lopsided tidal cavity (e.g. Kley & Dirksen [135]; Thun et al. [292]; Ragusa et al. [257]), it responds to a phenomenon of much wider extent. Indeed, disk eccentricity $e_{\mathrm{d}}$ can be significant ($\gtrsim0.01$) out to radii of $\sim10$–$15a_{\mathrm{b}}$ (Miranda et al. [208]; Muñoz & Lithwick [219]). Typically, this inner region undergoes coherent apsidal precession, with the rate given by the appropriate spatial average of Eq. (10) (see below). But in some cases, it can also become apsidally locked relative to the binary’s eccentricity vector, as found in the $0.2 \lesssim e_{\mathrm{b}} \lesssim0.4$ equal-mass binary simulations of Miranda et al. [208] and in the unequal-mass eccentric binary simulations of Siwek et al. [276]. An example of the coherent precession of a CBD and its correspondence to a lopsided cavity is illustrated in Figure 6.

The coherence of inner CBD precession requires efficient communication between different disk regions (recall that the apsidal precession rate of a test particle around the binary depends strongly on $r$; see Eq. (10)). In the absence of disk self-gravity, this communication is achieved by gas pressure, which induces disk precession that balances the forced precession (Eq. (10)) (Goodchild & Ogilvie [102]; Teyssandier & Ogilvie [287]). The eccentricity profile of the CBD actually corresponds to the eigenmode of a Schrödinger-like eccentricity equation of Goodchild & Ogilvie [102] (see also Shi et al. 2012, Teyssandier & Ogilvie 2016, Lee et al. 2019a,b, Muñoz & Lithwick 2020) and the global disk precession frequency is the associated eigenfrequency. The mode is trapped between two turning points that result from the sharply truncated CBD; as a result, most of the mode “lives” in the immediate vicinity of the cavity, thus obeying the local precession law (with the mode frequency approximately equal to $\varpi_{\mathrm{d}}$ evaluated at the cavity radius) (Muñoz

Eccentric Lindblad/Corotation Resonances and Mean-Motion Resonances

The term ”ELRs” are often used to describe Lindblad resonances that are associated eccentricities. ELRs and ”eccentric corotation resonances” are related to mean-motion resonances in celestial mechanics [222] as follows. Consider the perturbation with a pattern rotation frequency $\omega_{mn} = (m\Omega_b + n\kappa_b)/m$ (see Eq. (2)). A general Lindblad/corotation resonance occurs when

where $n' = 0$ gives the CR, and $n' = \pm1$ gives the (linear) LRs, and $|n'| > 1$ would involve nonlinear resonances. Using $\kappa_b = \Omega_b - \dot{\varpi}_b$ and $\kappa= \Omega- \dot{\varpi}$, where $\dot{\varpi}_b$ and $\dot{\varpi}$ are the apsidal precession rates of the binary and the disk, respectively, the above resonance condition becomes

This describes the $(m + n):(m + n')$ mean-motion resonance, and the strength of the interaction potential is $\propto e_b^{|n|}e^{|n'|}$.

& Lithwick 2020). Interestingly, the mode is so strongly confined, that pressure-induced precession plays a minor role.

While the eccentricity profile and precession frequency can be understood via a linear analysis, the actual amplitude of the eccentricity eigenmode depends on the details of eccentricity growth, damping, and saturation. The growth of eccentricity likely results from a combination of the direct “hydraulic pumping” from high Mach number streamers colliding directly with the inner edge of the disk as argued by Shi et al. [271] and the tidal excitation via eccentric Lindblad resonances (ELRs) (e.g., Hirose & Osaki [113]; Lubow [172]). ELRs is a parametric instability arising from the modulation of the disk particle’s epicyclic motion by the binary potential. As discussed in § 2, the gravitational potential from the binary on the disk can be decomposed into many harmonic components, each having a pattern rotation frequency $\omega_p = \omega_{mn} = N\Omega_b/m$, where $m, N = 1, 2, \cdots$ (and $n = N - m$). In the presence of this rotating potential, the epicyclic frequency $\kappa$ of the a disk fluid element attains a modulation term, proportional to $\cos m(\omega_p - \Omega)t$ (where $\Omega$ is the angular frequency of the disk). A parametric resonance occurs when $m(\omega_p - \Omega) \simeq2\kappa\simeq2\Omega$, i.e.,

The strength of the forcing depends on $e_b$, with the components $m = 2, N = 1$ (which has $\Phi\propto e_b$), and $m = 2, N = 2$ [with $\Phi\propto(1 - 5e_b^2/2)$] being the most important.

For a binary with finite eccentricity, the secular (orbital-averaged) interaction can also drive the disk eccentricity. If we define the complex eccentricity $\mathcal{E}_d(r,t) = e_d(r,t)\exp[i\varpi_d(r,t)]$ for the disk particle (at radius $r$) and $\mathcal{E}_b = e_b\exp(i\varpi_b)$ for the binary, then the secular contribution to the time evolution of $\mathcal{E}_d$ is given by

The apsidal precession rate $\omega_{db} = \dot{\varpi}_d$ (see Eq. (10)) is driven primarily by the quadrupole potential of the binary, while the eccentricity forcing rate $\nu_{db}$ by the octupole potential and is given by (e.g., Moriwaki & [208] (Nakagawa 2004)

In the presence of an eccentricity-damping force (and neglecting other hydrodynamical effects), this secular forcing tends to drive the disk toward a “forced” eccentricity, given by $E_{\mathrm{d}}=(\nu_{\mathrm{db}}/\omega_{\mathrm{db}})E_{\mathrm{b}}$; this forced eccentricity is apsidally aligned with the binary eccentricity.

Thus, the eccentricity evolution of the CBD around an eccentric binary can be quite complex, driven by the secular and resonant forcings from the binary and hydrodynamical effects. Simulations show that the disk generally exhibits apsidal precession, but the precession rate may not be constant and the line of apses can vary with $r$. In analogy to test particle orbit, the CBD should evolve according to coexisting “free” and “forced” modes (e.g., see Lubow 2022). In some cases, the CBD can lock onto the binary and stop precessing altogether [208, 276].

Angular Momentum Transfer and Orbital Evolution: Circular Binaries

We now discuss the important problem of long-term angular momentum exchange between the disk and the binary. Without accretion, the binary always loses angular momentum to the CBD through gravitational (Lindblad) torque. With accretion, the net torque becomes uncertain. As discussed above, the mass flow rate across the disk, $\dot{M}(r,t)$, and the mass accretion rate onto the binary $\dot{M}_{\mathrm{b}}(t)=\dot{M}_1(t)+\dot{M}_2(t)$, are all highly variable. For an extended disk with a constant supply rate $\dot{M}_0$ at $r=r_{\mathrm{out}}\gg a_{\mathrm{b}}$, a quasi-steady state is eventually reached, where the time-averaged mass flow rate $\langle\dot{M}(r,t)\rangle=\langle\dot{M}_{\mathrm{b}}(t)\rangle=\dot{M}_0$. In general, the net torque on the binary can be obtained in two ways:

• First, it can be computed directly as the sum of the gravitational torque from all the gas plus the accretion torque (due to momentum of the accreting gas onto each binary component), i.e.

When the size of the binary component is much less than $a_{\mathrm{b}}$, the spin torques $\langle\dot{S}_1\rangle$ and $\langle\dot{S}_2\rangle$ are negligible, and the total torque $\langle\dot{J}_{\mathrm{b}}\rangle$ acts on the binary orbit, i.e. $\langle\dot{J}_{\mathrm{b}}\rangle\simeq\langle\dot{L}_{\mathrm{b}}\rangle=\langle\dot{L}_{\mathrm{b}}\rangle_{\mathrm{grav}}+\langle\dot{L}_{\mathrm{b}}\rangle_{\mathrm{acc}}$.

• Second, the angular momentum flow rate in the circumbinary disk (at radius $r$ in the disk) can be computed as (Miranda et al. [208])

where $\dot{J}_{\mathrm{d,adv}}$ is the inward angular momentum advection rate, $\dot{J}_{\mathrm{d,visc}}$ is the outward viscous angular momentum transfer rate, and $\dot{J}_{\mathrm{d,grav}}$ is the torque from the binary acting on the gas exterior to radius $r$.

When the disk reaches a quasi-steady state, the time average $\langle\dot{J}(r,t)\rangle$ is independent of $r$, and the circumbinary disk has two global constants: $\langle\dot{M}(r,t)\rangle=\dot{M}_{\mathrm{b}}=\dot{M}_0$ and $\langle\dot{J}(r,t)\rangle=\langle\dot{J}_{\mathrm{b}}\rangle$.

Using long-term AREPO simulations, Muñoz et al. [220] demonstrated that the global quasi-steady state can be achieved for circumbinary accretion, and $\langle\dot{J}_{\mathrm{b}}\rangle$ computed using both methods are in agreement (see Figure 7). For equal-mass circular binaries (and with $q_{\mathrm{b}}=0.1$, $\alpha=0.1$), the specific angular momentum “eigenvalue” (i.e., the angular momentum transferred to the binary per unit accreted mass) is

where $\Omega_{\mathrm{b}}=(GM_{\mathrm{b}}/a_{\mathrm{b}}^3)^{1/2}$ is the rotation rate of the binary. The result was confirmed independently by Moody et al. [213] using the Athena code. [A similar positive value of $l_0$ was first obtained by Miranda et al. [208] based on PLUTO simulations with an excised cavity.]

For a binary accreting from a finite-sized disk/torus, a global quasi-steady state obviously does not exist. Instead, the accretion proceeds in two phases: an initial transient phase, corresponding to the filling of the binary cavity, followed by a viscous pseudo-stationary phase (achieved after the viscous time $r^{2}/\nu$ at a few disk radii), during which the torus viscously spreads and accretes onto the binary. In the transient phase, the torque on the binary is negative since it is entirely gravitational (see § 2). Muñoz et al. [218] demonstrated that in the viscous phase, the net torque on the binary per unit accreted mass is close to $l_{0}$, the value derived for “infinite” disks. Because no global steady-state is required, this method allows for a more efficient computation of the net torque on the binary in non-steady situations (i.e., when $\langle\dot{M}_{\mathrm{b}}\rangle$ gradually evolves in time).

Using angular momentum conservation, $dJ_{\mathrm{b}}/dt=\langle\dot{M}_{\mathrm{b}}\rangle l_{0}$ with $J_{\mathrm{b}}=\mu_{\mathrm{b}}a_{\mathrm{b}}^{2}\Omega_{\mathrm{b}}$, we have, for a $q_{\mathrm{b}}=1$ binary,

Thus the binary orbit expands at the rate $\dot{a}_{\mathrm{b}}/a_{\mathrm{b}}\simeq2.44\langle\dot{M}_{\mathrm{b}}\rangle/M_{\mathrm{b}}$. Note if we take account of the small amount of spin torque $\dot{S}_{\mathrm{b}}=\dot{S}_{1}+\dot{S}_{2}$, the actual $\dot{L}_{\mathrm{b}}$ is a bit smaller than $\dot{J}_{\mathrm{b}}$, and the resulting $\dot{a}_{\mathrm{b}}/a_{\mathrm{b}}$ is then smaller [Muñoz et al. [220] found $\dot{S}_{\mathrm{b}}\simeq0.028\dot{M}_{\mathrm{b}}a_{\mathrm{b}}^{2}\Omega_{\mathrm{b}}$ when assuming a “stellar” radius of $0.02a_{\mathrm{b}}$, and thus $\dot{a}_{\mathrm{b}}/a_{\mathrm{b}}\simeq2.2\dot{M}_{\mathrm{b}}/M_{\mathrm{b}}$.]

Lessons from Numerical Studies of Long-term Binary Evolution. Until recently, it has always been thought that circumbinary accretion leads to binary orbital decay. Pringle [251] suggested that the strong gravitational torque from the binary on the CBD prevents accretion, and therefore the binary loses angular momentum to the disk (see § 2). Gas accretion changes the picture completely. Only a few previous studies addressed the issue of angular momentum transfer in a quantitative way. Examples include MacFadyen & Milosavljević [182] and Shi et al. [271], whose simulations excise the inner cavity. MacFadyen & Milosavljević [182] considered $H/r = 0.1$ and a disk viscosity with $\alpha= 0.01$, and adopted a polar grid in the domain between $r_{\mathrm{in}} = a_{\mathrm{b}}$ and $r_{\mathrm{out}} = 100a_{\mathrm{b}}$. They found a reduction of mass accretion onto the binary and the dominance of the gravitational torque relative to advective torque (therefore a negative net torque on the binary). However, with their small $\alpha$ parameter the “viscous relaxation” radius at $t = 4000P_{\mathrm{b}}$ (the typical duration of their runs) is only about $3a_{\mathrm{b}}$, and their surface density profile is far from steady state even at $r \sim r_{\mathrm{in}}$. So it is likely that the result of MacFadyen & Milosavljević [182] reflects a “transient” phase of their simulations. Shi et al. [271] obtained a positive value of $\dot{J}_{\mathrm{b}}$ in their 3D MHD simulations of CBDs (truncated at $r_{\mathrm{in}} = 0.8a_{\mathrm{b}}$). However, the duration of their simulations is only $\sim100P_{\mathrm{b}}$, and it is unlikely that a quasi-steady state is reached. Their value of $l_0$, which is too small to cause orbital expansion, may not properly characterize the long-term evolution of the binary. Tang et al. [286] carried out simulations of accretion onto circular binaries using DISCO (Duffell & MacFadyen [81], Duffell [79]) with $h = 0.1$ and $\alpha= 0.1$. For the accretion prescription, they assumed that inside a “sink” radius (measured from each “star”), the gas is depleted at a rate $\mathrm{d}\Sigma/\mathrm{d}t = -\Sigma/t_{\mathrm{sink}}$, with $t_{\mathrm{sink}}$ a free parameter. They claimed that the net torque on the binary is negative unless $t_{\mathrm{sink}}$ is much less than $P_{\mathrm{b}}$. This result is in contradiction with Muñoz et al. [220] and Moody et al. [213], the latter adopted a similar accretion prescription and did not find the same behavior as Tang et al. [286]. More recent works using DISCO (Duffell et al. [80], Dittmann & Ryan [69], D’Orazio & Duffell [73]) and the cartesian-grid code MARA3 (Tiede et al. [293], Zrake et al. [323]) have produced similar results as Muñoz et al. [220, 218] and Moody et al. [213], so it appears that different groups have reached agreement (at least for disks with $h = 0.1$, $\alpha= 0.1$). See below for discussion on the numerical resolution requirement to obtain reliable net torques.

Dependence on Disk Thickness and Viscosity

Most of the recent simulations described thus far in this review are two-dimensional, low-mass, viscous accretion disks with locally isothermal equation of state. In this case, the most important parameters are: $q_{\mathrm{b}}$, $e_{\mathrm{b}}$, $\alpha$ (or $\nu$) and $h$. Of these, $q_{\mathrm{b}}$ and $e_{\mathrm{b}}$ have been studied most extensively, but $h$ and $\alpha$ can be equally important.

An exploration of the dependence of $l_0$ and $\dot{a}_{\mathrm{b}}/a_{\mathrm{b}}$ on $h$ for equal-mass, circular binaries was carried out by Tiede et al. [293]. Using the Godunov code Mara3 (implemented in Cartesian coordinates with static mesh refinement), they varied $h = \{0.02, 0.033, 0.05, 0.1\}$, while fixing the kinematic viscosity $\nu= \sqrt{2} \times10^{-2}a_{\mathrm{b}}^2\Omega_{\mathrm{b}}$ (or $\alpha= 10^{-2}h^{-2}$ at $r = 2a_{\mathrm{b}}$). They found that $l_0$ monotonically decreases with decreasing $h$ and that it fall below $(3/8)a_{\mathrm{b}}^2\Omega_{\mathrm{b}}$ (see Equation (19)) when $h \lesssim0.04$. This occurs in spite of the CBD mean cavity remaining roughly of the same size (since they use constant $\nu$ instead of constant $\alpha$, the truncation estimate of Miranda & Lai [207] would yield identical cavity sizes), albeit the streamers become more erratic, narrow, and denser than in thicker disks. The authors reported a dependence the results on resolution, which is likely due to the poorly resolved CSDs even with mesh refinement. Dittmann & Ryan [70] carried out follow-up work using a version of DISCO, and found that the Tiede et al. [293] results were likely unconverged but qualitatively confirmed the migration transition at $h = 0.04$. Dittmann & Ryan [70] also explored the role of viscosity. While earlier works (Muñoz et al. [218], Duffell et al. [80]) found that the $h = 0.1$ simulations exhibited little dependence on viscosity, when $h$ is varied, lower viscosity appears to further reduces $l_0$ with decreasing $h$ (see Figure 8, right panel). The role played by disk thickness could introduce a crucial distinction between the young stellar binaries accreting from protostellar disks and massive BH binaries accreting from AGN disks. While the former systems can be well represented by disks with $h \sim0.1$, the latter are better characterized by values of $h \sim10^{-2}$ or less (Sirko & Goodman [275], Thompson et al. [289]).

Two difficulties must be addressed when discussing the roles of disk viscosity and aspect ratio. First, a low-viscosity disk simply takes longer to reach a quasi-steady state (Miranda et al. [208], Muñoz et al. [218], Dittmann & Ryan [70]), and this is likely the culprit of simulations reporting suppressed binary accretion at low viscosities (e.g., Ragusa et al. [258]). To obtain reliable results, simulations often need an iterative reassessment of the initial condition (Miranda et al. [208]; Dempsey et al. [61]). Second, the resolution requirements in the circum-single disks (CSDs) have not been rigorously addressed, and these can be very stringent. Insufficient resolution would dilute the strong spiral arms in the CSDs, which are the source of positive torques opposing the negative torque from the CBD and the streams. And the width of these spiral arms can be very sensitive to the choice of $h$. Borrowing a page from the theory of planet-disk interaction (e.g., see Papaloizou et al. [236]), we know that the one-sided Lindblad torque peaks at an azimuthal number $\sim1/(2h)$ (Ward [309]) and the width of a spiral arm is $\sim4\pi r h^2$ (e.g., Masset [189]), setting a minimum resolution requirement for the adequate torque calculation. While this resolution requirement can be easily met by modern simulations in the CBD region, it can become extremely difficult in the CSD region, where the number of resolution elements per radial interval around each accretor can be rather small. Moreover, in the frame of the moving accretor, the computational grid is never of a polar nature, and in many cases, is effectively Cartesian at the CSDs scales. Consequently, numerical diffusion can become very taxing in the CSDs, especially for low values of $h$ where the spiral arms are narrow. Under-resolving these spiral arms has the unwanted consequence of reducing the positive torque stemming from this region, favoring the (well-resolved) negative torques from the CBD.

Recent SPH simulations carried out by Heath & Nixon [112] exhibit a drastically different torque reversal threshold of $h = 0.2$. These simulations however, do not exhibit any substantial CSD formation, and are thus prone to underestimate the positive torque stemming from the CSD region. Franchini et al. [94] showed that the SPH code PHANTOM can only reproduce the results Muñoz et al. [220] if the CSD structure is properly resolved with the use of $10^7$ particles. On the other hand, they confirmed that the meshless particle-based code GIZMO (Hopkins [115]) can produce outward migration if the number of resolution elements is increased within the cavity region.

Long-Term Orbital Evolution: Eccentric Binaries

The secular binary migration rate of Equation (19) assumes that $e_{\rm b} = 0$ at all times. For general eccentric binaries, Muñoz et al. [220] devised a method to compute from simulations both $d l_{\rm b}/dt$ and $d\mathcal{E}_{\rm b}/dt$, where $\mathbf{l}_{\rm b} = \mathbf{r}_{\rm b} \times\mathbf{v}_{\rm b}$ (with $\mathbf{r}_{\rm b} = \mathbf{r}_1 - \mathbf{r}_2$ the binary separation vector, and $\mathbf{v}_{\rm b} = d\mathbf{r}_{\rm b}/dt$) and $\mathcal{E}_{\rm b} = \frac{1}{2}\mathbf{v}_{\rm b} \cdot\mathbf{v}_{\rm b} - GM_{\rm b}/|\mathbf{r}_{\rm b}| = -GM_{\rm b}/(2a_{\rm b})$ are the specific orbital angular momentum and energy of the binary, respectively. In particular, the energy transfer rate can be computed from

where $\mathbf{f}_1 = \mathbf{f}_{1,\mathrm{grav}} + \mathbf{f}_{1,\mathrm{acc}}$ is the force per unit mass on $M_1$ (from gravity and accretion), and similarly for $\mathbf{f}_2$. From these, one can obtain the orbital evolution rates

The results of Muñoz et al. [220] are shown as blue pentagons in Figure 9. They concluded that equal-mass binaries of different eccentricities expand at the rate $\dot{a}_{\rm b}/a_{\rm b} \sim\dot{M}_{\rm b}/M_{\rm b}$ (left panel). But while $\dot{a}_{\rm b}/a_{\rm b}$ is always positive, the non-monotonic feature of $\dot{a}_{\rm b}/a_{\rm b}$ as a function of $e_{\rm b}$ is of interest. Related to this feature is the evolution of $e_{\rm b}$: Muñoz et al. [220] found that $\dot{e}_{\rm b}$ (right panel) is slightly negative for $e_{\rm b} \sim0$, becomes positive at $e_{\rm b} = 0.1$, and becomes negative again at $e_{\rm b} \gtrsim0.5$. Thus it appears that there is an eccentricity “attractor” below $e_{\rm b} \sim0.5$ and above $e_{\rm b} \sim0.3$. This behavior was confirmed by an independent study of Zrake et al. [323] (see Figure 9), who used the MARA3 code to study accretion from finite-sized disks onto eccentric binaries finely sampling a range of values of $e_{\rm b}$, and reported an equilibrium binary eccentricity of $e_{\rm b} \sim0.45$. These results are shown in the right panel of Figure 4 (red circles) alongside an associated fitting function (gray curve). The figure also includes the results of D’Orazio & Duffell [73] for $\dot{a}_{\mathrm{b}}$ and $\dot{e}_{\mathrm{b}}$ (purple squares), who found agreement with Muñoz et al. [220] and Zrake et al. [323] for $\dot{e}_{\mathrm{b}}$, but reported negative $\dot{a}_{\mathrm{b}}$ for intermediate values of $e_{\mathrm{b}}$. These authors argued that eccentric binaries close to $e_{\mathrm{b}} = 0.4$ are able to migrate inward because of the non-axisymmetric distortion of the circumbinary cavity for those parameters.

Accretion onto Unequal-Mass Binaries

Early SPH simulations of young stellar binaries accreting from gaseous environments [e.g. 17, 16, 18] revealed that the accretion flow depends sensitively on the binary mass ratio $q_{\mathrm{b}} = M_2/M_1$, and that even over short timescales, the secondary grows in mass faster than the primary [e.g. 16]. Farris et al. [90] conducted the first systematic study of circumbinary accretion as a function of $q_{\mathrm{b}}$ for circular binaries using an early version of the code DISCO. These simulations confirmed the existence of preferential accretion onto the secondary, and also found that the time-variability of accretion is a function of $q_{\mathrm{b}}$. Using a similar setup, Muñoz et al. [218] carried out AREPO simulations for $q_{\mathrm{b}}$ between 0.1 and 1 until a quasi-steady state was reached. Although they found some quantitative discrepancies with Farris et al. [90], they confirmed that both preferential accretion and the accretion variability depend on $q_{\mathrm{b}}$.

Currently, a consensus has emerged across simulation studies of accreting circular binaries. For fiducial parameters $h = 0.1$ and $\alpha= 0.1$, the key findings are

• Accretion is primarily modulated at two frequencies: $\Omega_{\mathrm{b}}/5$ and $\Omega_{\mathrm{b}}$. For $q_{\mathrm{b}} \ge0.5$, the accretion variability is dominated by the lower frequency mode (a period of about $5P_{\mathrm{b}}$); for $q_{\mathrm{b}} \le0.4$, it is dominated by the higher frequency mode (a period of $P_{\mathrm{b}}$) (see Figure 10). The switch is associated to the disappearance of the orbiting “lump” at low binary mass ratios.

• The secondary accretes more than the primary. The long-term preferential accretion ratio, defined as $\eta= \langle\dot{M}_2\rangle/\langle\dot{M}_{\mathrm{b}}\rangle$, is shown in Figure 11 from a collection of recent studies. While some

discrepancies remain, the monotonically decreasing trend of $\eta$ vs $q_b$ is robust. The result of Muñoz et al. [218], follows an approximately linear relation

• The accretion angular momentum eigenvalue $l_0=\langle\dot{J}_b\rangle/\langle\dot{M}_b\rangle$ stays around $0.68a_b^2\Omega_b$ to $0.8a_b^2\Omega_b$ for $q_b$ in the range 0.1 to 1. Figure 12 collects the results reported by Muñoz et al. [218] and Dittmann & Ryan [69], which are in remarkable agreement. Dittmann & Ryan [69] further showed that $l_0$ depends weakly on the sink prescription, provided that the sink region is small enough. Note that this positive eigenvalue does no guarantee binary expansion, as some accreted angular momentum goes into equalizing the binary’s mass ratio. For $\eta$ given by Eq. (23), the critical eigenvalue above which orbital expansion occurs is

which reduces to $l_{0,\mathrm{crit}}=(3/8)a_b^2\Omega_b$ for $q_b=1$. Equation (24) is shown in Figure 12 as a gray line, hinting at binary contraction for $q_b\approx0.1$.

• The binary orbit expands at a nearly constant rate $\langle\dot{a}_b\rangle/a_b\simeq2\langle\dot{M}_b\rangle/M_b$ for $0.3\leq q_b\leq1$; this rate becomes significantly smaller for $q_b\leq0.2$ and changes sign for $q_b\lesssim0.1$. The

left panel of Figure 13 shows the binary migration rate as computed by Muñoz et al. [218]. As $l_{0,\mathrm{crit}}$ grows for small $q_b$ (Figure 12), $\langle\dot a_b\rangle$ approaches zero, hinting at a reversal of binary migration. Duffell et al. [80] probed this transition by running DISCO simulations that dynamically update the value of $q_b$, covering values down to $q_b=10^{-3}$. They found that the net gravitational torque becomes negative for $q_b\lesssim0.05$, implying inward migration. Recent results by Dempsey et al. [64] for $q_b\ll1$ suggest that this torque reversal transition depends on the dimensionless parameter¹ $K' \equiv q_b^2/(\alpha h^3)$ and is associated with the disk becoming eccentric. For $h=0.05$ and $\alpha=10^{-3}$, Dempsey et al. [64] report that the gravitational torque goes from negative to positive when $q_b\gtrsim1.5\times10^{-3}$, or when $K'$ surpasses 20. For $\alpha>0.03$ and $h=0.1$, $K'\gtrsim20$ is equivalent to $q_b\gtrsim2\times10^{-2}$, in agreement with Duffell et al. [80].

Other Complications

In the previous subsections, we have reviewed the key results from CBD simulations using idealized equation of state (locally isothermal with constant disk aspect ratio $h$) and viscosity prescription ($\alpha$ viscosity or constant $\nu$), and assuming small (local) disk mass (without disk self-gravity). These assumptions may not always apply to realistic disks.

The locally isothermal equation of state is known not to conserve the angular momentum flux carried by density waves ([160, 209]), with a resulting modification of the torque profile in the disk. The effect on the integrated torque, however, is negligible for rapidly cooling disks, and amounts to a correction of $\sim30\%$ in slowly cooling/adiabatic disks ([210]).

Thus, the net contribution of this effect to the migration of binaries is likely minor, except for the borderline cases.

The locally isothermal equation of state assumes an efficient cooling process that instantly brings disk temperature back to equilibrium. Without resorting to radiative transfer, the simplest way to assess the role of gas thermodynamics is to allow the disk temperature to relax towards the equilibrium value over some cooling timescale $t_{\mathrm{cool}}$ (usually parameterized using $\beta= \Omega_{\mathrm{K}}t_{\mathrm{cool}}$, where $\Omega_{\mathrm{K}}$ is the local Kepler frequency). Wang et al. [307] carried out ATHENA++ simulations using this $\beta$-cooling prescription, and found that with a longer cooling time, the accretion variability is gradually suppressed and the morphology of the inner disk becomes more symmetric. They also found that the accretion angular momentum eigenvalue $l_0$ generally decreases with increasing $\beta$, such that an equal-mass, circular binary (with equilibrium $h \sim0.1$) undergoes orbital decay for $\beta\gtrsim0.2$, although the details depend on the viscosity prescription.

A more complex behavior is expected to arise from 3D magnetohydrodynamical (MHD) simulations, which are needed to capture viscous dissipation self-consistently. Such MHD simulations for circumbinary accretion are typically initialized with narrow, finite-supply tori, and are run for a small number of orbits. Thus, they are useful for studying short-term accretion variability, but may not be adequate for probing the long-term flow behavior and the secular evolution of the central binary. Early MHD simulations of CBDs around quasi-circular binaries ([271], using a version of the finite difference code ZEUS) and GRMHD ([226], using HARM3D) were carried out excising the central binary for a run time of 80 binary orbits. These show the similar accretion modulation behavior as 2D viscous hydrodynamics, despite the development of MHD turbulence in the disk. More recent efforts have been focused on the general relativistic aspects of CBD evolution, important for black-hole binaries close to merger (with separations of 10-20 gravitational radii). For example, GRMHD simulations using HARM and its variants have been implemented with an approximate analytic metric of binary black holes in order to include the CSDs in the computational domain, although a circular region at the binary barycenter is still excised ([35, 36, 50, 51]). These simulations (typically run for $\sim20$ orbits) show additional time modulations that are interpreted as purely relativistic effects.

The simulations discussed in the previous subsections all assume that the local disk mass near the binary is much less than $M_{\mathrm{b}}$ (although the mass of an extended disk can be much larger). When this condition is not satisfied, the flow/binary dynamics and evolution can be quite different. Such massive disks (or “envelopes”) naturally occur immediately following galaxy mergers. A massive black-hole binary embedded in such an envelope is expected to undergo orbital decay through dynamical friction (e.g., [87, 53, 95]). Massive CBDs or envelopes may also be important for the formation of close (au-scale) proto-stellar binaries (see § 4.2). If CBDs are massive enough and no substantial CSDs can form, the binary is bound to lose angular momentum to the gas. Numerical simulations show that small-mass proto-binaries embedded in collapsing cores decay in separation quickly (e.g., [16, 18]). Aided by dynamical friction, the pair becomes bound at small separations and continue to grow until it becomes dynamically dominant. Only then a true CBD may form, at which point the binary may have already reached au-scale separations. More recent simulations, including the effects of magnetic fields, find qualitatively similar results, except for a further accelerated orbital decay attributed to magnetic braking of the gas envelope (e.g., [319, 147]). For example, in the simulations of Kuruwita et al. [147], the binary pair decreases in separation from 300 au to 10 au in only $\sim30$ radial incursions, at which point it becomes a Keplerian bound pair and its orbital decay is slowed down. In this scenario, orbital decay occurs early in the star formation process, and circumbinary accretion represents the late stage of binary formation.

APPLICATIONS OF CIRCUMBINARY ACCRETION

Massive Black Hole Binaries and the Final Parsec Problem

In the $\Lambda$CDM paradigm of cosmology, MBH encounters are expected to occur as a natural consequence of hierarchical structure growth. If most galaxies host a MBH in their nuclei [183, 139], then MBH binaries (MBHBs) should be an inevitable consequence of galaxy mergers [22]. Indeed, the discovery of dual AGNs with kiloparsec separations confirm that galactic mergers can form wide MBH pairs [see 58 for a review]. Likewise, since galaxy mergers funnel large amounts of gas toward the center of the merger remnant [e.g. 116], the formation of gaseous CBDs should occur alongside that of MBHBs. It is generally assumed that, if MBHBs form, they evolve toward merger one way or another [for a review, see 49]. Such mergers are of great importance for low-frequency GW astronomy in the LISA band [108, 2] and the Pulsar Timing Array (PTA) band (Burke-Spolaor et al. 2019).

Following galaxy mergers, dynamical friction by gas/stars and stellar “loss-cone” scatterings would serve to reduce the separations of MBHBs to about one parsec, while gravitational radiation would become effective (and make the BHs merge) only if their separations are less than about 0.01 parsec [e.g. 22, 247]. Binary orbital evolution in the intermediate separation regime is a long-standing problem known as the “final parsec problem” [204 a].

CBDs have been posited to be a gas-dynamical solution to this final parsec problem (stellar-dynamical solutions, on the other hand, have been thoroughly studied in the literature; e.g., [316, 304, 106]. In the CBD hypothesis, a MBHB embedded in a viscous accretion disk is forced to migrate inward along with the accretion flow, much like massive planet embedded in a protoplanetary disk would [103]. Most if not all cosmological-scale studies of MBHB coalescence that consider the influence of gas assume CBDs always promote inwardly driven migration [109, 127, 129, 306, 305]. The studies of CBD accretion of recent years put this standard assumption into question, potentially introducing a major caveat to the expected detection rates for LISA [e.g. 267] and PTAs [e.g. 268].

Despite the encouraging discovery of wide dual AGN [e.g. 104], the direct observation of spatially resolved, sub-pc MBHBs will remain all but impossible for the foreseeable future (with the possible exception of very long baseline interferometric observations of nearby galaxies; Burke-Spolaor 2011). Thus, until a direct GW measurement of the inspiral and merger of two MBHs is made, indirect detection methods will be needed to probe the elusive sub-pc regime.

Indirect Signatures of Ongoing Circumbinary Accretion onto Massive Black-Hole Binaries. Spatially resolving MBHBs at sub-parsec separations is difficult. The most compact, directly-imaged MBHB on record is the dual AGN in the radio galaxy 0402+379, with a projected separation 7.3 pc [261]. For more compact sources, photometric and spectroscopic techniques can be used to infer the presence of MBHBs. Circumbinary accretion may prove crucial for uncovering multi-messenger efforts to identify such compact MBHBs in the form of periodically varying electromagnetic (EM) counterparts [e.g. 29, 46]. The sought-after periodic signals in time-variable AGN include (i) time-variable kinematic offsets of the broad emission lines [98], (ii) Doppler boosted broad-band emission from the circum-secondary accretion disk [76], (iii) periodic self-lensing flares [72] (Hu et al. 2019), and (iv) broad-band photometric variability due to pulsed/modulated accretion [e.g. 285]. In practice, the kinematics approach is of limited use at sub-pc because of the confusion and truncation of the distinct broad line regions at these separations [126], making the photometric approach a more viable and data-rich technique [e.g. 46]. Indeed, systematic searches for periodic signals in the Catalina Real-time Transient Survey [105], the Palomar Transient Factory (Charisi et al. 2016), the Panoramic Survey Telescope and Rapid Response System [169] have revealed hundreds of MBHB candidates.

Measuring photometric variability is not exempt from systematic uncertainties [105, 170, 313], but has the potential of being a powerful multi-messenger tool in combination with simultaneous GW detections [46]. At sub-pc separations, MBHBs emit GWs in the frequencies suitable for PTAs, which are forecast to provide individual binary detections within the next decade [206, 129]. In this regime, the binary orbit is still essentially Newtonian, and the associated GW waveforms are well understood. Likewise, the known photometric variability mechanisms (accretion/boosting/lensing) are all modulated on timescales commensurate with the binary orbit. Simultaneous GW/photometric detection of these objects would not only conclusively prove their existence, but could also allow for constrains on the their orbital parameters.

Imprints of Circumbinary Accretion on Massive Black-Hole Binary Orbits. In addition to detecting circumbinary accretion “in action” in the form of the photometric variability of individual sources, CBD physics can also leave an imprint in the population of MBHBs responsible for the stochastic gravitational wave background (GWB) [e.g. 322, 290, 259, 241], which is detectable by the method of pulsar timing [264, 68, 25]. The GWB spectrum can be computed from the GWs emitted from MBHB mergers across cosmic scales [241], although the uncertainties in its amplitude and spectral slope depend on the astrophysics of MBHB assembly [268] and on whether or not binary coalescence is delayed, suppressed, or biased toward a specific range in binary mass $M_{\rm b}$ [e.g. 269, 201]. Moreover, recent calculations also conclude that the GWB can depend on the eccentricity distribution of the binaries [128], and on the preferential accretion ratio [277]. Consequently, the physics of circumbinary accretion can have a major impact in the GWB.

Binary Star formation

The existence of tidally cleared disks [e.g. 83, 191, 123, 43, 120, 57] or rotationally supported structures [294, 194] around young stellar binaries, plus the confirmation that such binaries may accrete in a modulated fashion [14, 122, 295, 297] (Muzerolle et al. 2013), have underscored the importance of binary-disk interaction in early stellar evolution. These stellar CBDs also contain information about circumbinary planet formation (see § 4.3 below) and can serve as laboratories and proxies for the circumbinary physics taking place in galactic nuclei (§ 4.1).

Tidal truncation and pulsed accretion in young binaries (already discussed in § 3) are the most straightforward smoking guns for circumbinary accretion. But other important clues of past circumbinary accretion may be found in the evolved population as well. These include binary mass ratios and the over-abundance of stellar twins, and migration history of compact and wide binaries.

The Over-Abundance of Stellar Twins. Over the last few decades, several studies have consistently shown that stellar binaries exhibit an over-abundance of twins (with mass ratio $q_{\rm b} > 0.95$; [180, 299, 110, 179, 256, 211, 212]). This finding appears to be in agreement with some early hydrodynamical simulations which showed that, in binaries accreting from rotating gas, the secondary grew faster that the primary (Bate 2000, but see also [6, 33]). Indeed, high resolution CBD simulations that systematically compute the ratio of accretion rates in circular binaries find that the secondary accretes more than the primary and the ratio $M_2/(\dot{M}_1+\dot{M}_2)$ decreases monotonically with increasing $q_{\rm b}$ (see Figure 11).

The twin excess was originally reported to be most significant at short orbital periods $P_{\rm b} \lesssim20$ d or $a_{\rm b} \lesssim0.4$ au, (Moe & Di Stefano 2017; see also Tokovinin 2000, Raghavan et al. 2010), which was interpreted as being consistent with the once-held idea that accreting binaries migrate inward (see discussion in [256]. In recent years, however, astrometric measurements (using Gaia data) have revealed that the twin excess fraction is roughly constant for projected binary separations between 0.01 au and 10000 au ([86]). This finding contradicts the assumption that equal-masses and compact orbits are two sides of the same coin. Instead, the large separation of twins may indicate that outward migration due to CBDs has occurred. But even if CBDs promote softer orbits in addition to equal masses, they alone may not produce binary orbits as wide as $\sim10^{4}$ au. El-Badry et al. [86] conjectures that twins form first with separations $a_b \lesssim100$ au while accreting from circumbinary disks, and are subsequently widened by dynamical interactions in their parental clusters. The fact that these binaries also tend to be eccentric ([118]) adds credence to the hypothesis that these orbits have been subject to strong perturbations after gas has dissipated.

Formation of Close Binaries. In light of the emerging paradigm that binaries accreting from warm disks tend to expand (see § 3), the origin of 1-10 au-scale binaries remains a puzzle. Numerical simulations suggest that binary stars initially form with large separations ($\gg10^{2}$ au) (see Offner et al. [228]). The difficulty lies in the impossibility of fragmenting a disk at small separations where cooling is inefficient ([192, 253, 30, 311, 281, 52, 143]), and in the impossibility of directly fragmenting a core when these are stable against non-axisymmetric perturbations (and because of the size of the first hydrostatic core is already a few au in size; Bate [15]). Naturally, a significant inward migration would be needed to produce close binaries.

A likely solution to this conundrum is that the required migration takes place at the early (Class 0-I) stage of star formation, when the newly fragmented binaries are still embedded in massive envelopes (see Tokovinin & Moe [298]). As discussed in § 3.9, some simulations of 3D collapse with embedded binaries do show binary orbital decay at rates faster than the Lindblad torques of the associated CBDs. In some cases, these “binaries” decay because they are not a bound pair initially, and their orbital evolution is dominated by dynamical friction and gas accretion from the envelope.

Planets Around Binaries

One of the most exciting results in exoplanetary science in recent years was the discovery of circumbinary planets (CBPs). The Kepler mission has discovered 13 CBPs around 11 eclipsing binaries [e.g. 78, 310, 279] and the TESS mission has so far (Oct. 2022) detected two CBP systems ([140, 141]). Most of these are large planets (with radius between $0.3R_{\mathrm{J}}$ and $R_{\mathrm{J}}$) in a sub-au, nearly co-planar orbit around the host binary. A handful of CBP systems have been discovered using gravitational microlensing (OGLE-2007-BLG-349, Bennett et al. [24]), and direct imaging (e.g., HD 106906 Bailey et al. [10], Lagrange et al. [149], Rodet et al. [260]). In addition, a number of CBP systems have been inferred using eclipse timing variations (ETVs) of binaries (e.g., NN Serpentis Qian et al. [252], Beuermann et al. [26]), although the validity of these planets remains uncertain.

Recent studies suggest that the occurrence rate of large, Kepler-like CBPs is comparable to that of similar-mass planets in single-star systems ($\sim10\%$ Armstrong et al. [5], Martin & Triaud [184], Martin et al. [185]), indicating that planet formation in circumbinary disks is a robust process.

Close-in planetary orbits around a binary are known to be dynamically unstable [e.g. 84, 246, 71]. For binaries with $e_b \sim1$, the critical planetary semi-major axis $a_{\mathrm{crit}}$ is about a few times $a_b$. Holman & Wiegert [114] provides an approximate expression of $a_{\mathrm{crit}}$ based on simulation of circular, co-planar circumbinary particles:

where $\mu_2 = M_2/M_{\mathrm{b}}$. A striking feature of the CBPs discovered by Kepler is that many of them lie close to the instability boundary (e.g. Kepler-16b has a semi-major axis $a_p = 1.09a_{\mathrm{crit}}$, and Kepler-34b has $a_p = 1.14a_{\mathrm{crit}}$). This feature cannot be explained by the selection bias of transit observations (Li et al. 2016). It has been commonly interpreted as evidence for planetary migration, since the circumbinary cavity acts as a “trap” in which migrating planets can be “parked” [e.g. 190, 244, 245, 136, 137, 240, 48]. The precise stopping location of the planet (for a given binary) depends on various disk parameters (such as viscosity and disk aspect ratio), which affect the intrinsic cavity size, and on the planetary mass – A massive planet opens a gap in the disk and circularizes the inner cavity, and thus tends to migrate closer to the binary. Current hydrodynamical simulations struggle to reproduce systems such as Kepler-34b, an eccentric ($e_{\mathrm{p}} = 0.18$) planet orbiting very close to an eccentric binary ($e_{\mathrm{b}} = 0.52$), because the central cavity is large and has a significant eccentricity, causing the planet to park too far away [e.g. 245, 240]. The eccentricity of the inner CBD (see § 3.3) can also leave an imprint on embedded CBPs. In particular, low-mass planets can inherit the eccentricity of the gas disk and precess with it in a state of apsidal corotation [291, 240]. Clearly, the close-in CBPs provide clues on the intricate interplays of binary-disk and planet-disk interactions.

In-situ formation of CBPs at close-in locations are generally difficult because of the large impact velocities between planetesimals driven by the perturbations on the planetesimal orbits from the binary potential and non-axisymmetric density features within the CBD. Such large impact velocities lead to corrosive collisions and prevent the buildup of large planetary bodies [266, 188, 200, 232]. The motion of planetesimals, however, can depend on the CBD as much as on the central binary. Rafikov [254] and Silsbee & Rafikov [273] found that the (secular) planetesimal dynamics is affected by the gravity and gas drag from a precessing eccentric CBD. These effects (for sufficiently massive CBDs) may suppress the excitation from the binary, and couple to it resonantly, defining specific regions of the CBD where planetesimals can grow into planets. Gas drag is fundamental in permitting planet formation under external eccentricity perturbations because it can lead to coherent eccentric orbits [273] and/or coherent “most circular” closed orbits [243, 163, 315, 39]. Planetesimals in such orbits would still have low relative velocities, meaning that growth via collisions would not be hampered by the global eccentricity. Overall, these studies suggest that close-in CBPs formed at relatively large distances ($\gtrsim10a_{\mathrm{b}}$) and move to their current orbits through disk-driven migration.

Post-Main-Sequence Binaries

The final evolution of 0.8-8$M_{\odot}$ stars involves a rapid transition from the Asymptotic Giant Branch (AGB) over the post-AGB phase towards the planetary nebulae stage. Many post-AGB stars are in binary systems (with a main-sequence companion). It is now well-established that most post-AGB binaries (with periods from 100 to a few thousands days) are surrounded by circumbinary gas-dust disks; the observational evidence comes from the disk-like spectral energy distribution of the systems, IR imaging and CO interferometry that resolves the Keplerian velocity field [e.g. 303, 302, 138]. Many of these binaries also have a bipolar outflow/jet launched from the circumstellar disk around the secondary (main-sequence) companion [e.g. 32, 31]. The circumstellar disk is likely fed from the CBD, although mass transfer from the post-AGB primary star cannot be ruled out.

The origin of the CBDs around post-AGB binaries is somewhat unclear. These disks are formed as a result of binary interactions during the AGB phase. One possibility is that during the common-envelope (CE) stage, not all envelope is ejected, and a fraction of the mass falls back to the remnant binary in the form of a CBD [e.g. 125] (Sandquist et al. 1998). One issue with this scenario is that CE evolution is expected to produce systems with periods ranging from hours to hundreds of days, and yet no CBD has been observed around systems with periods less than 100 days. Another possibility is mass loss through the outer Lagrangian point $L_{2}$ associated with the mass transfer or wind of the AGB star [e.g. 272, 60, 238], although it is unclear if stable (and long-lasting) CBD can be produced in this way.

CBDs may produce some dynamical effects on the post-AGB binaries or other evolved binaries. Many post-AGB binaries with periods 100 to a few thousands days are observed to have significant eccentricities (up to $e_b \simeq0.6$) [302]. This is surprising since such binaries are supposed to have circularized during the AGB phase. Dermine et al. [67] suggested that the eccentricities of the post-AGB binaries could be excited as a result of their gravitational interaction with the CBD. Antoniadis [3] suggested that similar effect could explain the eccentricities of some binary millisecond pulsars with the white-dwarf companions. The eccentricity-growth mechanism can be easily understood from gravitational binary-CBD interaction [e.g. 101, 175]. From § 2, we know that the binary loses angular momentum and energy through the OLR associated with the $(m n)$ potential component at the rates $dJ_b/dt=-T_{mn}$ and $dE_b/dt=-\omega_{mn}T_{mn}$. These imply that the binary eccentricity evolves at the rate

For small but finite $e_b$, the $m=2$, $n=-1$ potential component has $\omega_{mn}=\Omega_b/3$, and $T_{mn}\propto e_b^2$, leading to $\dot{e}_b/e_b>0$.

It is important to note that equation (26) relies on pure gravitational interaction between the binary and the CBD. This could apply in the ”transient” stage before gas accretion onto the binary sets in [see, for example, the simulation of accretion from a finite-sized CBD by Muñoz et al. [218]]. However, as discussed in § 3.5, once the accretion starts (typically within a viscous time of the inner disk) and the system settles into a quasi-steady state, the secular evolution of the binary can be quite different from that predicted by the pure gravitational effect. In addition, when applying to post-AGB binaries, it is not clear that the CBD has a sufficient mass to change the binary eccentricity in an appreciable way even if the accretion effect is ignored [255]. Obviously, more works are needed to understand the origin of the peculiar eccentricities of the post-MS binaries.

MISALIGNED DISKS

In the previous sections, we have focused on aligned disks, i.e., disks that are coplanar with their central binaries. Such alignments may not be realized in many situations. For example, in the current theory of star formation, supersonic turbulence in molecular clouds leads to the creation of clumps, which then condense and fragment into stars and binaries [e.g. 196]. The gas that falls on to the central protostellar core/binary and assembles on to the disk at different times may rotate in different directions [19, 21, 227, 301, 91, 284]. In this scenario, it is reasonable to expect a newly formed binary to be surrounded by a highly misaligned CBD, which forms as a result of continued gas accretion. Similar chaotic accretion may also occur in the formation of massive black holes (MBHs) [e.g. 134]. In particular, MBH binaries at the centers of galaxies may experience accretion of successive low-mass gas “clouds” with uncorrelated angular momenta, which naturally lead to misaligned CBDs.

Observationally, most CBDs around young stellar binaries are found to be aligned with their host binary orbital planes [56]. For example, the gas rich CBDs AK Sco [54] and DQ Tau [55], and the debris disks $\alpha$ CrB and $\beta$ Tri [133], all have mutual disc–binary inclinations less than a few degrees. However, there are some notable exceptions. For example, the pre-MS binary HD 155555 possesses a low inclination ($10$–$20^\circ$) precessing circumbinary ring or disk [47, 312, 42, 248] (see Zhu et al. [320] for two similar systems discovered by ZTF). The disks (circumbinary and two circumstellar) in the system IRS 43 are highly misaligned ($\sim60^\circ$) with each other and with the binary [38]. The gaseous CBD around the eccentric binary HD 98000B is nearly polar (with a misalignment angle $\sim90^\circ$; Kennedy et al. [131]). The 6–10 Gyr old eccentric binary 99 Herculis has a nearly polar ($\sim87^\circ$) debris ring around it [132]. The young ($\sim1$ Myr) hierarchical triple star system GW Ori possesses a misaligned circumtriple disk with three broken dusty rings [144, 28].

Disk Warping, Breaking and Alignment

Consider a CBD surrounding a circular binary (with total mass $m_b$, reduced mass $\mu_b$, semi-major axis $a_b$ and eccentricity $e_b=0$). The orientation of the disk at radius $r$ (from the center of mass of the binary) is specified by the unit normal vector $\hat{\boldsymbol{l}}(r,t)$. Averaging over the binary orbital period and the disk azimuthal direction, the binary imposes a torque per unit mass on the disk element at radius $r$ given, to the leading order in $a_b/r$, by

where $\Omega=(GM_b/r^3)^{1/2}$ is the orbital frequency at radius $r$, $\hat{\boldsymbol{l}}_b$ is the unit vector along the binary angular momentum axis, and

characterizes the (nodal) precession rate of the disk mass element around the binary axis. Since $\omega_{\mathrm{prec}}$ depends strongly on $r$, the disk would lose coherence if different parts of the disk did not “communicate” with each other. In reality, the combination of differential precession and internal fluid stresses can give rise to a coherently warped/twisted disk.

Theoretical studies of warped disks [237, 235] have shown that there are two dynamical regimes for the linear propagation of warps in an accretion disk. For high-viscosity disks with $\alpha\gtrsim h \equiv H/r$ (where $H$ is the half-thickness of the disk and $\alpha$ is the Shakura–Sunyaev parameter such that the viscosity is $\nu=\alpha H^2\Omega$), the warp is communicated via angular momentum advection by the oscillatory internal flow whose amplitude is determined by the viscosity, and satisfies a diffusion-type equation. The corresponding diffusion coefficient $\nu_2$ measures the $r$-$z$ viscous stress and can differ from the usual viscosity $\nu$ (which measures the $r$-$\phi$ stress).² Ogilvie [229] has developed a fully nonlinear theory of diffusive warps that is in agreement with 3D numerical simulations [171]. In such high-viscosity regime, we expect that an inclined CBD at large radii transitions to alignment with the binary at small radii, with the characteristic transition radius (“warp radius”) $r_{\mathrm{warp}}$ determined by $\omega_{\mathrm{prec}}(r_{\mathrm{warp}})\simeq\nu_2/r_{\mathrm{warp}}^2$, giving

(This warp transition behavior is similar to the Bardeen-Petterson effect of accretion disks around spinning black holes; e.g., Bardeen & Petterson [11], Kumar & Pringle [146], Scheuer & Feiler [265]).

Depending on the initial condition, a misaligned viscous CBD may be susceptible to tearing (i.e. breaking up into two or more disconnected “rings”) if the steady warped state cannot be attained. Such disk tearing was observed in the SPH simulations of Nixon et al. [225] when the CBD is initiated with a sufficiently large inclination angle $\theta$ with respect to the binary (see also Nixon et al. 2020). The disk breaking radius can be estimated by comparing the viscous warp torque (per unit mass) $(\nu_2/r^2)r^2\Omega$ to the precessional torque from the binary, $|T_b|\sim\omega_{\mathrm{prec}}r^2\Omega|\sin2\theta|$, giving

Disk breaking requires sufficiently large $\theta$, for which $r_{\mathrm{break}}\sim r_{\mathrm{warp}}$. If the initial $\theta$ is small such that $r_{\mathrm{break}}$ lies inside the inner radius of the CBD, disk breaking would not occur and we expect the disk to evolve into the steady warped state, with a smooth transition between alignment (with the binary) at small radii to misalignment at large radii.

For low-viscosity disks with $\alpha\lesssim h$, a low-amplitude warp propagates as a bending wave at about half the sound speed, $c_{\mathrm{s}}/2$, provided that the disk is sufficiently Keplerian (i.e. the apsidal precession and nodal precession rates are less than $h\Omega$) (Lubow & Ogilvie 2000; see also Lubow et al. 2002, Ogilvie 2006). Protoplanetary disks around young binary stars (with $\alpha\sim10^{-4}$–$10^{-2}$, $h \gtrsim0.05$ and $r/a_{\mathrm{b}} \gtrsim$ a few) generally satisfy these conditions [92]. The nonlinear behavior of low-viscosity warped disks is complicated and poorly understood owing to the resonant excitation of vertical ”breathing” motions [88] and a parametric instability associated with the inertial waves [97, 231, 233, 66, 65].

Because of the efficient communication of warps by the bending waves, an inclined protostellar disk around a binary is approximately flat – this is the case provided that the bending wave crossing time $2r/c_{\mathrm{s}}$ is shorter than the characteristic precession time $\omega_{\mathrm{prec}}^{-1}$, a condition that is well satisfied everywhere in the disk except the inner-most region [317]. The interplay between the small (non-zero) disk twist/warp and viscous dissipation drives the long-term evolution of the disk inclination. Foucart & Lai [92, 93] studied the warp and the dissipative torque that drives the inclination evolution of a CBD around a circular binary. Foucart & Lai [92] considered an infinite disk and included the effect of accretion on to the binary, while Foucart & Lai [93] considered a more realistic disk of finite size and angular momentum, which can precess coherently around the binary. They showed that under typical protoplanetary conditions, both viscous torque associated with disk warp/twist and accretion torque tend to damp the mutual disk–binary inclination on time-scale much shorter than the disk lifetime (a few Myr). In contrast, a circumstellar disk inside a binary can maintain large misalignment with respect to the binary orbital plane over its entire lifetime [177, 93]. Qualitatively, the key difference between CBDs and circumstellar disks is that in former/latter, the binary torque is exerted at the inner/outer region of the disk, which contains small/larger amount of angular momentum, leading to relatively large/small warp and thus faster/slower viscous damping. Overall, these results are consistent with the observations that most CBDs are nearly coplanar with their host binaries [56], while circumstellar disks within young stellar binaries are often misaligned (e.g., [121, 119].

Polar Alignment of Disks Around Eccentric Binaries

Although disks around circular binaries tend to evolve toward alignment, recent works suggested that other outcomes may be possible for disks around eccentric binaries. Aly et al. [1] carried out SPH simulations of disks around eccentric MBHBs (which typically lie in the “viscous” regime of disk warps, with $\alpha\gtrsim h$), and showed that the disk may be driven into polar alignment (i.e. the disk plane is perpendicular to the binary plane). Martin & Lubow [186] found numerically that a circumbinary protoplanetary disk (typically in the bending-wave regime, with $\alpha\lesssim h$) inclined to an eccentric ($e_{\mathrm{b}} = 0.5$) binary by $60^\circ$ will evolve to a polar configuration. This dynamical outcome arises from the combined influences of the gravitational torque on the disk from the eccentric binary and viscous torque from disk warping.

To understand the possibility of polar alignment, it is useful to consider the secular dynamics of a circular test particle around an eccentric binary (e.g., [89, 152]).

where $\mathbf{e}_{\mathrm{b}}$ is the binary eccentricity vector. In the absence of hydrodynamical forces, the time evolution of the test particle’s orbital angular momentum vector $\hat{\mathbf{l}}$ is governed by $d\hat{\mathbf{l}}/dt=\mathbf{T}_{\mathrm{b}}/(r^{2}\Omega)$. It is clear that the evolution of $\hat{\mathbf{l}}$ has four possible fixed points (where $d\hat{\mathbf{l}}/dt=0$): $\hat{\mathbf{l}}=\pm\hat{\mathbf{l}}_{\mathrm{b}}$ and $\hat{\mathbf{l}}=\pm\mathbf{e}_{\mathrm{b}}$. To examine the stability of the fixed points, it is useful to analyse the trajectory of $\hat{\mathbf{l}}(t)$ using the “energy” curves. equation of motion for $\hat{l}(t)$ admits an integral of motion

which is simply related to the quadrupole interaction energy (double-averaged over the binary and test-particle orbits) by (e.g., Tremaine et al. [300]; Liu et al. [168])

Figure 14 shows the test particle trajectories (constant-$\Lambda$ curves) in the $I-\Omega$ (left panel) and $I_e-\Omega_e$ (right panel) planes for $e_b=0.3$. The critical separatrix $\Lambda=0$ is displayed in black in both plots. When $\Lambda>0$, $\hat{l}$ precesses around $\hat{l}_b$ with $I\sim$ constant and $\Omega$ circulating the full range $(0-360^\circ)$. When $\Lambda<0$, $\hat{l}$ precesses around $\boldsymbol{e}_b$ with $I_e\sim$ constant and $\Omega_e$ circulating the full range $(0-360^\circ)$. Thus, the test particle angular momentum axis $\hat{l}$ transitions from precession around $\hat{l}_b$ for $\Lambda>0$ to precession around $\boldsymbol{e}_b$ for $\Lambda<0$. A necessary condition for $\hat{l}$ to precess around $\boldsymbol{e}_b$ is $I_{\mathrm{crit}}<I<180^\circ-I_{\mathrm{crit}}$, where

Zanazzi & Lai [318] carried out a detailed theoretical analysis of the dynamics of inclined, warped disks around eccentric binaries and their long-term evolution driven by viscous torques (see also Lubow & Martin [176]). For disks with $H/r\gtrsim\alpha$ (as appropriate for protoplanetary disks), bending wave propagation effectively couples different regions of the disk, making it precess as a quasi-rigid body. Zanazzi & Lai showed explicitly that the dissipative torque (associated with disk warp/twist) tends to drive the disk to one of two states, depending on the initial sign of $\Lambda$: For $\Lambda>0$, the disk angular momentum axis $\hat{l}_d$ aligns (or anti-aligns) with the binary orbital angular momentum vector $\hat{l}_b$; for $\Lambda<0$, $\hat{l}_d$ aligns with the binary eccentricity vector (polar alignment). They also showed that when the disk has a non-negligible angular momentum compared to the binary, the system’s fixed points are modified and the disk may then evolve to a state of near polar alignment, with the inclination somewhat less than $90^\circ$ (see also Martin & Lubow [187]).

Note that $\Lambda$ depends on both $I$ (the disk-binary inclination) and $\Omega$ (the longitude of ascending node of the disk). Thus for a given $e_b$, the direction of inclination evolution depends not only on the initial $I(0)$, but also on the initial $\Omega(0)$. The time-scale of evolution of the disk–binary inclination angle is approximately given by (assuming the disk surface density profile $\Sigma\propto1/r$)

where $r_{\mathrm{in}},r_{\mathrm{out}}$ are the inner and outer radii of the disk. Thus $\tau_b$ is generally less than a few Myr, the lifetime of proto-planetary disks. This suggests that highly inclined disks may exist around eccentric binaries.

This theoretical expectation was recently confirmed by the observation of the young ($\sim10$ Myr) protoplanetary system HD 98800. Using ALMA observations of dust and CO emissions, Kennedy et al. [131] showed that the inner binary BaBb (with $a_b\simeq1$ AU, $e_b=0.785$) in the system (which is a “2+2” hierarchical quadruple system with two inner binaries orbiting each other with a semi-major axis of 54 AU) is surrounded by a CBD in a polar-aligned configuration. The old circumbinary debris disk nearly perpendicular to the orbital plane of binary 99 Herculis (Kennedy et al. [132]) may also have gone through such polar-alignment process early in its lifetime when the gas was present (e.g., Smallwood et al. [278]).

Figure 15 (from Czekala et al. [56]) shows the inclinations of Kepler circumbinary planets and circumbinary protoplanetary and debris disks, as a function of binary orbital period, semi-major axis and eccentricity. A clear trend emerges: All CBDs orbiting binaries with period less than 30 days (semi-major axis less than 0.4 AU) and/or eccentricity less than 0.2 are consistent with being coplanar, while disks orbiting longer period and/or more eccentric binaries exhibit a wide range of mutual inclinations, from coplanar to polar. The origin of this “critical” orbital period is unclear. But the trend is consistent with the general idea that stellar binaries form/fragment at large separations and migrate inwards to small separations in massive disks/envelopes, during which the binary and disk become aligned (see § 4.2).

The near coplarity of CBDs around short-period binaries implies that planets formed in such disks should be similarly coplanar. This is consistent with the finding that Kepler circumbinary planets, orbiting binaries with $P_{\mathrm{b}}<40$ days, have small mutual inclinations, and indicates that the planet occurrence rate around such binaries is similar to that for single stars (see also [5, 153, 185]). Beyond $P_{\mathrm{b}}>40$ days, however, the existence misaligned CBDs suggests that circumbinary planets around eccentric binaries may be found to have a broad distribution of mutual inclinations, with a possible concentration of polar-aligned systems.

BINARIES EMBEDDED IN “BIG” DISKS

A special type of “circumbinary accretion” has gained significant interest in recent years. This concerns stellar binaries embedded in AGN disks around supermassive black holes (SMBHs).

The detection of gravitational waves from the merging binary black holes (BHs) by the LIGO/Virgo collaboration (The LIGO Scientific Collaboration et al. 2021) has motivated many theoretical studies of the formation channels of the BH binaries. In addition to the isolated binary evolution channel (e.g., Lipunov et al. [162], Belczynski et al. [23], there are several flavors of dynamical channels, including strong gravitational scatterings in dense star clusters (e.g., Portegies Zwart & McMillan [249], Kremer et al. [145]), more gentle “tertiary-induced mergers” (often via Lidov-Kozai mechanism) that take place either in stellar triple/quadrupole systems (e.g., Silsbee & Tremaine [274], Liu & Lai [164, 165]) or in nuclear clusters dominated by a central SMBH (e.g. Antonini & Perets [4], Liu et al. [167], Liu & Lai [166]), and (hydro)dynamical interactions of binaries in AGN disks (e.g., Bartos et al. [12], Stone et al. [282], McKernan et al. [197], Tagawa et al. [283]). It has also been suggested that intermediate-mass BHs (IMBHs) may be formed efficiently in AGN disks via accretions or mergers of stars and compact objects (e.g., McKernan et al. [199, 198]). If these IMBHs are themselves members of binaries, they could be susceptible to “evection resonances” driven by the SMBH’s tidal field (Muñoz et al. [221]), which could then accelerate their GW-driven inspiral and produce an eccentric waveform (see also Bhaskar et al. [27]).

In the AGN disk scenario, an important question concerns how single stellar-mass BHs can be captured into bound binaries and merge. Li et al. (2022b) showed that when the gas effect is negligible, two BHs in tightly-packed orbits around a SMBH become bound to each other in rare, very close encounters due to gravitational wave emission, leading to highly eccentric BH mergers in the LIGO band. At sufficiently high gas densities, the gas drags on the BHs can facilitate the capture process (e.g., Tagawa et al. [283], Li et al. 2022a).

Once a bound BH binary forms in the AGN disk, an important question concerns how the binary evolves in the presence of the surrounding gas. It is tempting to consider this as a “circumbinary accretion” problem (e.g., Stone et al. [282]), with the background AGN disk feeding gas onto the CBD around the binary. For example, one might use the modified Bondi-Hoyle-Lyttleton accretion formula (e.g., Edgar [85]) to estimate the mass supply rate from the background AGN disk (with gas density $\rho_{\rm bg}$ and sound speed $c_{\rm s,bg}$) onto the binary, $\dot{m}_{\rm b} \sim\pi\rho_{\rm bg} v_{\rm eff} r_{\rm acc}\min(r_{\rm acc}, H)$, where $H$ is the disk scale-height and $r_{\rm acc}$ is the accretion radius $r_{\rm acc} \sim Gm_{\rm b}/v_{\rm eff}^{2}$, with $v_{\rm eff}^{2} \sim c_{\rm s,bg}^{2} + v_{\rm H}^{2}$ (the second term accounts for the “Hill velocity”, the velocity shear across the Hill radius). However, such an estimate could be quite misleading, as the strong velocity shear and flow angular momentum can significantly reduce the accretion rate compared to the Bondi-Hoyle-Lyttleton estimate (Li & Lai 2022a). In fact, even when the binary is replaced by a single object, the accretion rate can be much smaller than the Bondi-Hoyle-Lyttleton estimate because of the strong upstream velocity shear. Numerical simulations show that the accretion rate generally depends on the physical size of the accretor, indicating that the physics near the accretor can strongly influence the accretion flow (see Xu & Stone [314]). Thus, it is not clear that the results from circumbinary accretion (§ 3) can be directly adapted and applied to the problem of “Binaries embedded in big disks”.

The hydrodynamical evolution of binaries in AGN disks has been studied numerically by a handful of works so far. Baruteau et al. [13] carried out global simulations in 2D isothermal disks and found that a massive (gap-opening) prograde, equal-mass binary is hardened by dynamical friction from the lagging spiral tails trailing each binary component inside the Hill radius. Li et al. [159] used a similar global setup and found that adequately resolved circum-single disk (CSD) regions in fact lead to expanding binaries. Li et al. [158] further found that the temperature structure of the CSDs plays an important role in the evolution of the binary; when CSDs attain a sufficiently high temperature, binaries contract rather than expand.

Resolving CSDs around each binary component is therefore important, but is computationally demanding in global disk simulations. A useful approach is to use a co-rotating local disk (“shearing-box”) model, where the global cylindrical geometry of the disk is mapped onto local Cartesian coordinates centered at the binary’s center of mass which rotates around the SMBH. This is the approach adopted by Li & Lai (2022a,b) and Dempsey et al. [62]. Note that although there are multiple length scales and velocity scales associated with the problem, only a few dimensionless ratios are important. For example, the relevant length scales are:

• Binary semi-major axis $a_{\mathrm{b}}$;

• Hill radius $R_{\mathrm{H}} \equiv R(M_{\mathrm{b}}/M)^{1/3} \equiv Rq^{1/3}$, where $M_{\mathrm{b}}$ is the total mass of the binary, $M$ is the mass of the SMBH, and $R$ is the orbital radius of the binary around the SMBH (note that this definition of Hill radius differs from the standard one $R_{\mathrm{H}}' = R_{\mathrm{H}}/3^{1/3}$);

• Bondi radius $R_{\mathrm{B}} = GM_{\mathrm{b}}/c_{\mathrm{s,bg}}^2$, where $c_{\mathrm{s,bg}}$ is the sound speed of the background gas (far from the binary);

• Scale height of the background disk $H = c_{\mathrm{s,bg}}/\Omega_{\mathrm{K}}$, with $\Omega_{\mathrm{K}} = (GM/R^3)^{1/2}$.

However, their ratios depend on only two dimensionless parameters:

Similarly, the relevant velocity scales are $c_{\mathrm{s,bg}}$, $v_{\mathrm{b}} = (GM_{\mathrm{b}}/a_{\mathrm{b}})^{1/2}$, the velocity shear across the binary $V_{\mathrm{s}} = (3/2)\Omega_{\mathrm{K}}a_{\mathrm{b}}$, and $|\Delta V_{\mathrm{K}}| \sim h^2 V_{\mathrm{K}}$ (the deviation of the background gas velocity around the SMBH from the Keplerian velocity). The first three are related by the same dimensionless parameters:

For thin discs ($h \ll1$), $|\Delta V_{\mathrm{K}}|$ is very subsonic ($|\Delta V_{\mathrm{K}}|/c_{\mathrm{s,bg}} \sim h \ll1$) and is typically much smaller than $V_{\mathrm{s}}$. Therefore, the flow dynamics and the binary orbital evolution, when appropriately scaled, depend on various physical quantities only through two dimensionless parameters: $q/h^3 = (R_{\mathrm{H}}/H)^3$ and $\lambda$ (in addition to other obvious parameters such as the binary eccentricity and mass ratio). Note that the parameter $q/(3h^3)$ is the ratio of the binary mass $M_{\mathrm{b}}$ to the so-called thermal mass, $3h^3M$. When $R_{\mathrm{H}}/H \gtrsim1$ (which also implies $R_{\mathrm{B}} \gtrsim R_{\mathrm{H}}$, the flow onto the Hill sphere is quasi-2D, we expect a CSD to form around each binary component; 2D shearing-sheet simulations are appropriate in this regime. When $R_{\mathrm{H}}/H \lesssim1$, vertical flows onto the Hill sphere are important, 3D simulations are needed to accurately capture the flow structure [62]. Current models of AGN disks [275, 289] contain regions with a wide range of $R_{\mathrm{H}}/H$ ratio (for typical $M_{\mathrm{b}}/M$), from $\lesssim1$ to $\gg1$ (see Figure 1 of Dempsey et al. [62]). The dynamical stability of the binary requires $\lambda\gtrsim2$; even for $\lambda\gtrsim1$, the flow structure still depends on $R_{\mathrm{H}}/H$ and can be quite different from that around isolated binaries.

Figure 16 illustrates the flow structure of an equal-mass circular binary embedded in the AGN disk, from the 2D shearing-box simulation of Li & Lai (2022a) with $\lambda= 5$, $q/h^{3} = 1$ and the gas satisfying the $\gamma$-law equation of state ($\gamma= 1.6$). Grand spirals originated from the circumbinary flows extend to large distances along the shear flow directions; they can be considered as half bow shocks of the binary accretion with upstream flow gradients due to the shear. Some of these features are also found in the global simulations of [13, 159]. There are also horseshoe flows and the inner/outer shear flows around the binary. Such flow structures are similar to those observed in previous studies of a single accretor (e.g., [96, 321, 148, 9]). Due to the fast orbital velocity $v_{\mathrm{b}}$, the flow close to the binary is much more dynamic. Both CSDs contain two spiral shocks that drive accretion throughout the disk. Each CSD is then encompassed and attached by a small half bow shock, the tail of which is slingshot away along each grand spiral once a binary orbit; the propagation of such waves is visible in Figure 16. More examples of the flow properties (which depend on $\lambda$, $q/h^{3}$, $\gamma$, $M_{1}/M_{2}$ and $e_{\mathrm{b}}$) can be found in Li & Lai (2022a,b) and Dempsey et al. [62] (the latter considered 3D simulations with isothermal equation of state). Overall, the flow structure is qualitatively different from that of isolated “circumbinary accretion” discussed in § 3.

Li & Lai (2022a,b) have used a suite of high-resolution 2D shearing-box simulations of binaries embedded in AGN disks to determine the accretion dynamics and secular evolution of the binary, covering a range of values for the parameters $\lambda$, $q/h^{3}$, $\gamma$, $M_{1}/M_{2}$ and $e_{\mathrm{b}}$. The flows close to the binary (including the CSDs) are generally found to be more massive with decreasing $\gamma$ and with increasing $q/h^{3}$ and $\lambda$; they are also hotter and more turbulent with increasing $\gamma$, $q/h^3$ and $\lambda$. The time-averaged accretion rate $\langle\dot{M}_{\rm b}\rangle$ (in units of $\Sigma_{\rm bg}a_{\rm b}v_{\rm b}$) depends on the physical size of the accretor (see above), and monotonically decreases with $\gamma$ and increases with $q/h^3$. In general, circular comparable-mass binaries contract if the EOS is far from isothermal, with an orbital decay rate of a few times the mass doubling rate. When the EOS is close to isothermal ($\gamma= 1$), the binary orbit expands [however, retrograde binaries always experience orbital decay; see also Li et al. [159]]. Eccentric binaries tend to experience eccentricity damping. Prograde binaries with higher eccentricities or smaller mass ratios tend to have slower orbital decay rates, with some extreme cases exhibiting orbital expansion. Note that some of the quantitative results may be modified by 3D effects when $q/h^3 \lesssim$ a few (Dempsey et al. [62]). The accretion flows are highly variable, and the dominant variability frequency is the apparent binary orbital frequency (in the rotating frame around the central massive BH) for circular binaries but gradually shifts to the radial epicyclic frequency as the binary eccentricity increases. These calculations also suggest that the hardening timescales of the binaries are much shorter than their migration timescales in the AGN disk, for all reasonable binary and disk parameters. Overall, these studies show that the dynamics of binaries embedded in AGN disks is quite different from that of isolated binaries in their own CBDs. Obviously, current simulations are still idealized, and future works will be needed to assess the effects of additional physics (e.g. magnetic fields, radiation and accretion feedbacks) on the evolution of BH binaries in AGN disks.

SUMMARY AND FUTURE PROSPECTS

Accretion disks have long played a central role in many areas of astronomy, from bright X-ray sources associated with accreting compact objects in the Galaxy, to outflows and jets associated with active galactic nuclei, to various complex processes and phenomena associated with star and planet formation. When the central accreting object is replaced by a binary, a new set of dynamical behaviors become important, including Lindblad torques and cavity opening, accretion variabilities on different timescales, eccentricity driving and precession of the inner disk, preferential accretion and long-term binary evolution, and, in the case of misaligned accretion, disk warping and breaking, and secular evolution towards alignment or polar alignment, etc. All these dynamical behaviors have direct or indirect observational manifestations.

Since circumbinary accretion involves intrinsically 2D or 3D phenomena, numerical simulations are crucial to unravel its dynamical behaviors and obtain quantitative answers to some of the key questions, with semi-analytic theory providing complementary insights. These simulations are challenging (compared to “normal” accretion disk simulations) because of the wide range of spatial (from $\gg a_{\rm b}$ for the circumbinary region to $\ll a_{\rm b}$ for the circum-single disks) and time scales (from $\ll P_{\rm b}$ to $\gg P_{\rm b}$) involved. For example, to determine the secular effect of accretion on the binary evolution, sufficiently long simulations must be carried out in order to average out the highly dynamical flow behaviors. As a result, systematic numerical simulations have so far focused on idealized 2D setups, with simple equation of state (usually locally isothermal) and viscosity prescription. Some of key findings include:

• Circumbinary accretion is highly dynamical, and the accretion variability is dominated by the period of $P_{\rm b}$ or $\sim5P_{\rm b}$, depending the binary eccentricity and mass ratio.

• The inner region ($\lesssim10a_{\rm b}$) of the circumbinary disk can develop coherent eccentric structure, which may modulate the accretion and affect the physical processes (such as planet migration) taking place in the disk.

• Over long timescales, a binary undergoing accretion evolves towards equal masses between the components.

• While the gravitational torque between the binary and circumbinary disk tends to drive binary orbital decay, once the inner disk region is sufficiently relaxed and accretion sets in, the net angular momentum transfer between the binary and the disk depends on the competition between the accretion and gravitational torques. Circumbinary accretion does not necessarily lead to binary orbit decay,

as commonly assumed; the secular orbital evolution depends on the binary parameters (such as the mass ratio and eccentricity) and the thermodynamic properties of the accreting gas. In general, the binary orbital decay/expansion rate is of order the mass doubling rate, i.e. $\dot{a}_{b}/a_{b} \sim\pm\dot{M}_{b}/M_{b}$.

Much work remains to go beyond these “idealized” simulations in order to systematically evaluate the roles of gas thermodynamics and magnetic fields (which can generate turbulent viscosity and outflows). Note that the above results apply to the regime where the local disk mass near the binary is much less than $M_{b}$. When this condition is not satisfied, the flow and binary dynamics can be quite different. For example, in the early stages of young binary stars and massive black-hole binaries, the local disk/envelope surrounding the binary can be much more massive than the binary, and the disk self-gravity is important – such massive disks can drive rapid binary orbital decay. Eventually, the local disk ”thins out”, and we are back to the ”proper” circumbinary accretion regime that determines the late-stage binary evolution.

On the observational front, the study of accreting proto-stellar binaries and the detailed characterization of “mature” binary stars (e.g. using GAIA) can shed light on the circumbinary accretion process in the aftermath of star formation. Similarly, the observation of accreting massive black-hole binaries (e.g. through variable lightcurves) and the future detection of low-frequency gravitational waves from such binaries (by the Pulsar Timing Arrays and by the space interferometers such as LISA) will help constrain the assembly history of massive black holes and the role of circumbinary accretion.

Circumbinary disk can often be misaligned with the orbital plane of the central binary, as expected in realistic scenarios of star formation and massive black hole evolution. Such misaligned accretion gives rise several new dynamical features:

• A misaligned disk is warped, and the “degree” of the warp depends on the internal hydrodynamical stresses of the disk (viscosity and bending waves). Under some conditions, a highly warped disk may break up into two or more separate “rings”.

• Over long (viscous) timescales, a misaligned disk around a low-eccentricity binary tends to evolve toward alignment driven by viscous dissipation. When the binary eccentricity is significant, the circumbinary disk can evolve toward “polar alignment”, with the disk plane perpendicular to the binary plane.

Numerical simulations of misaligned disks are challenging because they require capturing the large-scale disk warps while resolving small-scale dissipations. Much remains to be done in the future, e.g. to evaluate the possibility and condition of disk breaking. Observationally, polar-aligned circumbinary disks have already been detected, and future characterization of misaligned disks would help constrain the binary formation/migration process. Will highly misaligned planets be detected in the future?

Binaries embedded in a “big” disk present a special type of circumbinary accretion. The mass supply/accretion onto such a binary depends on the intricate coupling between the small-scale (within the binary) and large-scale flow dynamics. The orbital evolution of the binary depends on several dimensionless parameters as well as the gas thermodynamics. It has been suggested that some merging black-hole binaries detected by LIGO/VIRGO are produced by binaries embedded in AGN disks. Future observations, including more merger events and related electromagnetic counterparts, may provide a more definitive answer.

DISCLOSURE STATEMENT

The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.

Acknowledgments

The writing of this review has been motivated by several seminars/colloquia that DL gave on this topic in the last few years, most recently at KITP, UCB, Caltech, UIUC, KIAA and TDLI. We thanks our collaborators, including Xuening Bai, Francois Foucart, Lars Hernquist, Katlin Kratter, Jiaru Li, Rixin Li, Yoram Lithwick, Adam Dempsey, Ryan Miranda, Volker Springel, Haiyang Wang and J.J. Zanazzi, for their contributions and insights. We also thank Adam Dempsey, Julian Krolik, Gordon Ogilvie, Eliot Quataert and Noam Soker for useful comments on an early version of this article. This work has been supported in part by the NSF grant AST-17152 and NASA grant 80NSSC19K0444, and by Cornell University.

Literature Cited