Hydrodynamics of Circumbinary Accretion: Angular Momentum Transfer and Binary Orbital Evolution

The orbital angular momentum of the binary is subject to torques due to accretion and gravitational interactions:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dL}}_{{\rm{b}}}}{{dt}} & = & {\left.\displaystyle \frac{{{dL}}_{{\rm{b}}}}{{dt}}\right|}_{\mathrm{acc}}+{\left.\displaystyle \frac{{{dL}}_{{\rm{b}}}}{{dt}}\right|}_{\mathrm{grav}}\\ & = & {\dot{\mu }}_{{\rm{b}}}{l}_{{\rm{b}}}+{\mu }_{{\rm{b}}}\displaystyle \frac{{{dl}}_{{\rm{b}}}}{{dt}}{\left|\mathop{}\limits_{\mathrm{acc}}+{\mu }_{{\rm{b}}}\displaystyle \frac{{{dl}}_{{\rm{b}}}}{{dt}}\right|}_{\mathrm{grav}},\end{array}\end{eqnarray} \tag{ 2 }$$

where ${l}_{{\rm{b}}}={a}_{{\rm{b}}}^{2}{{\rm{\Omega }}}_{{\rm{b}}}{(1-{e}_{{\rm{b}}}^{2})}^{1/2}$ is the binary's specific angular momentum and is defined by ${L}_{{\rm{b}}}={\mu }_{{\rm{b}}}{l}_{{\rm{b}}}$, where ${\mu }_{{\rm{b}}}\,={M}_{1}{M}_{2}/{M}_{{\rm{b}}}={q}_{{\rm{b}}}{M}_{{\rm{b}}}/{(1+{q}_{{\rm{b}}})}^{2}$ is the binary's reduced mass. For the low-mass disks of interest in this paper, the quantities ${M}_{{\rm{b}}}$, ${q}_{{\rm{b}}}$, and ${l}_{{\rm{b}}}$ in Equation (2) change a negligible amount over the timescales of the simulations (i.e., these quantities evolve over a long, secular timescale, $\sim {M}_{{\rm{b}}}/{\dot{M}}_{{\rm{b}}}$), and we determine their time derivatives from simulation output in postprocessing. Hence, the binary evolution is not "live" (see Section 2.1 above).

The rate of change of total angular momentum in the binary ${\dot{J}}_{{\rm{b}}}$ should also include the angular momentum transfer to the spins of the "stars." Thus

$$\begin{eqnarray}&&{\dot{J}}_{{\rm{b}}}={\dot{L}}_{{\rm{b}}}+{\dot{S}}_{1}+{\dot{S}}_{2}.\end{eqnarray} \tag{ 3 }$$

Below, we describe in detail how the different contributions to Equations (2) and (3) are computed from numerical simulations.

Mass change: Even in the absence of specific torques in Equation (2), the angular momentum of the binary changes via changes in the reduced mass ${\mu }_{{\rm{b}}}$. Theses changes are simply related to the changes in mass ratio and total mass:

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{\mu }}_{{\rm{b}}}}{{\mu }_{{\rm{b}}}}=\displaystyle \frac{{\dot{M}}_{{\rm{b}}}}{{M}_{{\rm{b}}}}+\displaystyle \frac{1-{q}_{{\rm{b}}}}{1+{q}_{{\rm{b}}}}\displaystyle \frac{{\dot{q}}_{{\rm{b}}}}{{q}_{{\rm{b}}}},\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{b}}}={\dot{M}}_{1}+{\dot{M}}_{2}\,\,\,\mathrm{and}\,\,\,\,\,\,{\dot{q}}_{{\rm{b}}}=\displaystyle \frac{1+{q}_{{\rm{b}}}}{{M}_{{\rm{b}}}}\left({\dot{M}}_{2}-{q}_{{\rm{b}}}{\dot{M}}_{1}\right).\end{eqnarray} \tag{ 5 }$$

As in ML16, we compute the accretion rate by keeping track of ${\rm{\Delta }}{m}_{\mathrm{acc},i}(t)$, the cumulatively accreted mass at time t, over the course of the simulation. We compute ${\dot{M}}_{i}$ in postprocessing by taking the time derivative using a centered-differences approximation:

$$\begin{eqnarray}&&{\dot{M}}_{i}(t)\equiv \displaystyle \frac{d}{{dt}}{\rm{\Delta }}{m}_{\mathrm{acc},i}(t).\end{eqnarray} \tag{ 6 }$$

Specific gravitational torque: The specific gravitational torque acting on the binary has the general form

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{d{{\boldsymbol{l}}}_{{\rm{b}}}}{{dt}}\right|}_{\mathrm{grav}} & = & \displaystyle \frac{{M}_{1}}{{\mu }_{{\rm{b}}}}{{\boldsymbol{r}}}_{1}\times {{\boldsymbol{f}}}_{\mathrm{grav},1}+\displaystyle \frac{{M}_{2}}{{\mu }_{{\rm{b}}}}{{\boldsymbol{r}}}_{2}\times {{\boldsymbol{f}}}_{\mathrm{grav},2}\\ & = & {{\boldsymbol{r}}}_{{\rm{b}}}\times ({{\boldsymbol{f}}}_{\mathrm{grav},1}-{{\boldsymbol{f}}}_{\mathrm{grav},2}),\end{array}\end{eqnarray} \tag{ 7 }$$

where ${{\boldsymbol{r}}}_{{\rm{b}}}={{\boldsymbol{r}}}_{1}-{{\boldsymbol{r}}}_{2}$, and ${{\boldsymbol{f}}}_{\mathrm{grav},i}$ is the force per unit mass acting on the primary/secondary:

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{grav},i}=-{ \mathcal G }\displaystyle \sum _{k}{m}_{k}\displaystyle \frac{{{\boldsymbol{r}}}_{1}-{{\boldsymbol{r}}}_{k}}{| {{\boldsymbol{r}}}_{1}-{{\boldsymbol{r}}}_{k}{| }^{3}},\end{eqnarray} \tag{ 8 }$$

where m_k and ${{\boldsymbol{r}}}_{k}$ are the mass and position vector of the kth computational cell, and the sum extends over all cells in the simulation domain. A proper account of the acceleration history of the "stars" requires a densely sampled force evaluation in time. In practice, we follow the gravitational "kicks" that act on the primary and secondary, evaluating them twice per time step, mimicking the kick–drift–kick operation that one would implement for a truly "live" binary orbit integration. As with the accreted mass, at each time t, we store the cumulative velocity change ${\rm{\Delta }}{{\boldsymbol{\upsilon }}}_{\mathrm{grav},i}(t)$. In postprocessing, the gravitational force per unit mass is computed as

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{grav},{\rm{i}}}=\displaystyle \frac{d}{{dt}}{\rm{\Delta }}{{\boldsymbol{\upsilon }}}_{\mathrm{grav},i}(t).\end{eqnarray} \tag{ 9 }$$

This approach is based on momentum conservation; as such, it takes into account all time steps in the simulation and does not rely on the instantaneous forces computed from sparse simulation snapshots.

Specific accretion torque: Accretion itself can induce a specific torque on the binary (e.g., see Roedig et al. 2012), separate from the torque inherent to the change in binary mass. Similar to Equation (7), we have

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{l}}}_{{\rm{b}}}}{{dt}}\left|\Space{0ex}{3ex}{0ex},| ,\mathop{}\limits_{\mathrm{acc}}={{\boldsymbol{r}}}_{{\rm{b}}}\times ({{\boldsymbol{f}}}_{\mathrm{acc},1}-{{\boldsymbol{f}}}_{\mathrm{acc},2})\right..\end{eqnarray} \tag{ 10 }$$

We compute the accretion forces by keeping track of the velocity kicks as mass is incorporated into the "stars." As one would do in "live" accretion integrations, conservation of linear momentum dictates that the velocity of an accreting body of mass M_i is updated after each time step according to

$$\begin{eqnarray}&&\delta {{\boldsymbol{v}}}_{\mathrm{acc},i}=\displaystyle \frac{{{\boldsymbol{\upsilon }}}_{i}(t){M}_{i}+\delta {{\boldsymbol{p}}}_{\mathrm{acc},i}}{{M}_{i}+\delta {m}_{\mathrm{acc},i}}-{{\boldsymbol{\upsilon }}}_{i}(t)\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&\approx \,\displaystyle \frac{\delta {{\boldsymbol{p}}}_{\mathrm{acc},{\rm{i}}}}{{M}_{i}}-\displaystyle \frac{\delta {m}_{\mathrm{acc},{\rm{i}}}}{{M}_{i}}{{\boldsymbol{\upsilon }}}_{i}(t),\end{eqnarray} \tag{ 11b }$$

where $\delta {m}_{\mathrm{acc},i}$ and $\delta {{\boldsymbol{p}}}_{\mathrm{acc},i}$ are the changes in mass and momentum during a single time step. The accumulation of $\delta {{\boldsymbol{\upsilon }}}_{\mathrm{acc},i}$ over time then defines ${\rm{\Delta }}{{\boldsymbol{\upsilon }}}_{\mathrm{acc},i}(t)$. As in Equation (9), the accretion force is computed in postprocessing as

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{acc},i}=\displaystyle \frac{d}{{dt}}{\rm{\Delta }}{{\boldsymbol{\upsilon }}}_{\mathrm{acc},i}(t),\end{eqnarray} \tag{ 12 }$$

from which we compute the associated torques (Equation (10)). This specific torque vanishes if accretion onto the "star" is isotropic, in which case, only the mass—but not the velocity—changes. Therefore, we typically refer to the torque resulting from Equation (12) as the anisotropic accretion torque.

Note that, had we stored the cumulative evolution of accreted linear momentum ${\rm{\Delta }}{{\boldsymbol{p}}}_{\mathrm{acc},i}(t)$ instead of ${\rm{\Delta }}{{\boldsymbol{\upsilon }}}_{\mathrm{acc},i}(t)$, the total accretion torque could be written as

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{d{{\boldsymbol{L}}}_{{\rm{b}}}}{{dt}}\right|}_{\mathrm{acc}} & = & {{\boldsymbol{r}}}_{1}\times \displaystyle \frac{d}{{dt}}{\rm{\Delta }}{{\boldsymbol{p}}}_{\mathrm{acc},1}(t)+{{\boldsymbol{r}}}_{2}\times \displaystyle \frac{d}{{dt}}{\rm{\Delta }}{{\boldsymbol{p}}}_{\mathrm{acc},2}(t)\\ & \approx & {\mu }_{{\rm{b}}}{{\boldsymbol{r}}}_{{\rm{b}}}\times ({{\boldsymbol{f}}}_{\mathrm{acc},1}-{{\boldsymbol{f}}}_{\mathrm{acc},2})+({{\boldsymbol{r}}}_{{\rm{b}}}\times {{\boldsymbol{\upsilon }}}_{{\rm{b}}})\displaystyle \frac{d{\mu }_{{\rm{b}}}}{{dt}},\end{array}\end{eqnarray} \tag{ 13 }$$

where ${{\boldsymbol{\upsilon }}}_{{\rm{b}}}={\dot{{\boldsymbol{r}}}}_{{\rm{b}}}={{\boldsymbol{\upsilon }}}_{1}-{{\boldsymbol{\upsilon }}}_{2}$. Note that the second term in Equation (13) is already accounted for by the first term in Equation (2), which describes the accretion-induced torque at constant specific angular momentum.

Spin torque: Gas cells typically orbit around the "stars" before being drained, therefore carrying a small amount of angular momentum relative to a moving center. This residual angular momentum is irrevocably lost from the system due to accretion (e.g., Bate et al. 1995). The amount of angular momentum—or "spin"— ${\rm{\Delta }}{{\boldsymbol{S}}}_{\mathrm{acc},i}$ that is transferred to the "star" can be tracked as follows: each time some mass and momentum are drained from the kth gas cell, the spin angular momentum $({{\boldsymbol{r}}}_{k}-{{\boldsymbol{r}}}_{i})\times {f}_{k}{{\boldsymbol{p}}}_{k}$ (where f_k is the fraction of mass and linear momentum drained from the kth cell) is added to ${\rm{\Delta }}{{\boldsymbol{S}}}_{\mathrm{acc},i}$.⁵ As before, the spin torque is computed in postprocessing as

$$\begin{eqnarray}&&{\dot{{\boldsymbol{S}}}}_{\mathrm{acc},i}=\displaystyle \frac{d}{{dt}}{\rm{\Delta }}{{\boldsymbol{S}}}_{\mathrm{acc},i}(t).\end{eqnarray} \tag{ 14 }$$

In summary, in order to compute ${\dot{L}}_{{\rm{b}}}$ and ${\dot{J}}_{{\rm{b}}}$ (Equations (2) and (3)) from simulation output, we need ${\dot{M}}_{{\rm{b}}}$ and ${\dot{q}}_{{\rm{b}}}$ from Equations (5) and using Equation (6), $({{dl}}_{{\rm{b}}}/{dt})\left|{}_{\mathrm{grav}}\right.$ from Equation (7), $({{dl}}_{{\rm{b}}}/{dt})\left|{}_{\mathrm{acc}}\right.$ from Equation (10), and ${\dot{S}}_{\mathrm{1,2}}$ from Equation (14).

To determine the secular orbital evolution of an eccentric binary, we also need to compute the rate of change of the orbital energy, which is fully determined from ${\dot{M}}_{{\rm{b}}}$, ${{\boldsymbol{f}}}_{\mathrm{grav},{\rm{i}}}$, and ${{\boldsymbol{f}}}_{\mathrm{acc},{\rm{i}}}$. This is discussed in Section 4.

In quasi-steady-state, the time-averaged angular momentum transfer rate across the CBD, $\langle {\dot{J}}_{{\rm{d}}}(R,t)\rangle$, is constant in R and equal to the angular momentum transfer rate onto the binary $\langle {\dot{J}}_{{\rm{b}}}\rangle$. The balance of angular momentum currents throughout the CBD is (see Appendix A of MML17)

$$\begin{eqnarray}&&{\dot{J}}_{{\rm{d}}}={\dot{J}}_{{\rm{d}},\mathrm{adv}}-{\dot{J}}_{\mathrm{dvisc}}-{\dot{J}}_{{\rm{d}},\mathrm{grav}},\end{eqnarray} \tag{ 15 }$$

where ${\dot{J}}_{{\rm{d}},\mathrm{adv}}$, ${\dot{J}}_{{\rm{d}},\mathrm{visc}}$, and ${\dot{J}}_{{\rm{d}},\mathrm{grav}}$ are the advective, viscous, and gravitational contributions, respectively (in the notation of MML17, ${\dot{J}}_{{\rm{d}},\mathrm{grav}}$ is ${T}_{\mathrm{grav}}^{\gt r}$). The computation of the angular momentum current in the disk ${\dot{J}}_{{\rm{d}}}$ from the AREPO output is more involved than in the PLUTO case, since there is no polar grid. In order to obtain ${\dot{J}}_{{\rm{d}},\mathrm{adv}}$, ${\dot{J}}_{{\rm{d}},\mathrm{visc}}$, and ${\dot{J}}_{{\rm{d}},\mathrm{grav}}$, we first compute 2D maps of angular momentum fluxes ${F}_{J}(R,\phi )$ by remapping or interpolating the primitive variables of AREPO output onto a structured polar grid. The currents $\dot{J}(R)$ are then calculated from

$$\begin{eqnarray}&&\dot{J}(R)={\int }_{0}^{2\pi }{F}_{J}(R,\phi ){Rd}\phi =2\pi \langle {F}_{J}(R,\phi )R{\rangle }_{\phi }.\end{eqnarray} \tag{ 16 }$$

The advective, viscous, and gravitational fluxes are, respectively,

$$\begin{eqnarray}&&{F}_{J,\mathrm{adv}}=-R{\rm{\Sigma }}{\upsilon }_{R}{\upsilon }_{\phi }=-{\rm{\Sigma }}\displaystyle \frac{\left[{xy}({\upsilon }_{y}^{2}-{\upsilon }_{x}^{2})+{\upsilon }_{x}{\upsilon }_{y}({x}^{2}-{y}^{2})\right]}{\sqrt{{x}^{2}+{y}^{2}}}\,,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\begin{array}{rcl}{F}_{J,\mathrm{visc}} & = & -{R}^{2}\nu {\rm{\Sigma }}\left[\displaystyle \frac{\partial }{\partial R}\left(\displaystyle \frac{{\upsilon }_{\phi }}{R}\right)+\displaystyle \frac{1}{{R}^{2}}\displaystyle \frac{\partial {\upsilon }_{R}}{\partial \phi }\right]\\ & = & -\displaystyle \frac{\nu {\rm{\Sigma }}}{\sqrt{{x}^{2}+{y}^{2}}}\left[2{xy}\left(\displaystyle \frac{\partial {\upsilon }_{y}}{\partial y}-\displaystyle \frac{\partial {\upsilon }_{x}}{\partial x}\right)\right.\\ & & \left.+\ ({x}^{2}-{y}^{2})\left(\displaystyle \frac{\partial {\upsilon }_{x}}{\partial y}+\displaystyle \frac{\partial {\upsilon }_{y}}{\partial x}\right)\right],\end{array}\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}&&{F}_{J,\mathrm{grav}}=-{\int }_{R}^{{R}_{\mathrm{out}}}{\rm{\Sigma }}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{b}}}}{\partial \phi }{dR}^{\prime} =-{\int }_{R}^{{R}_{\mathrm{out}}}{\rm{\Sigma }}\left(-{{xa}}_{y}+{{ya}}_{x}\right){dR}^{\prime} ,\end{eqnarray} \tag{ 19 }$$

where we have made the dependence on Cartesian coordinates (the ones used by AREPO) explicit. In Equation (19), ${\boldsymbol{a}}=({a}_{x},{a}_{y})=-{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{b}}}$ is the cell-centered acceleration of the gas due to the binary's potential ${{\rm{\Phi }}}_{{\rm{b}}}$. These maps are computed from simulation snapshots, stacked to produce a time average and then integrated over azimuth (Equation (16)) to produce the time-averaged angular momentum currents of Equation (15).

In order to test how close our initial condition is to a true steady state, we carry out three test integrations with ${e}_{{\rm{b}}}=0$: (1) we run an AREPO simulation with an open, diode-like boundary condition at $R={a}_{{\rm{b}}}$ for $3000{P}_{{\rm{b}}}$ and then "release" the boundary to allow for direct accretion onto the binary, integrating for an additional 500 orbits (see ML16); (2) we allow for direct accretion onto the "stars" immediately at the start-up of the simulation and integrate for 500 orbits; (3) we use the results of the PLUTO simulation from MML17 at $t=3000{P}_{{\rm{b}}}$, remapping it onto an AREPO grid and rescaling the units such that the measured accretion rate at $R=15{a}_{{\rm{b}}}$ equals unity; then we proceed to evolve this setup by integrating, as in cases (1) and (2), for another 500 orbits. The accretion rates corresponding to these three tests are depicted in Figure 1. In all cases, after gas crosses the $R={a}_{{\rm{b}}}$ boundary, only 5–20 orbits are needed for the CSDs to form. Another 20–30 orbits are needed for the binary to start accreting with an averaged rate equal to ${\dot{M}}_{0}$ to within a few percent. Eventually, ${\dot{M}}_{0}$ is matched to within 1% in all cases, although this can take from 100 orbits (top panel) to 300 orbits (bottom panel). Overall, these three simulations are indistinguishable from each other, and we conclude that once quasi-steady state is reached, the initial conditions become irrelevant.

Zoom In Zoom Out Reset image size

From here on, we define quasi-steady state as the state of the simulation for which the averaged accretion rate onto the binary $\langle {\dot{M}}_{{\rm{b}}}\rangle$ over an interval $T=100{P}_{{\rm{b}}}$ is equal to the mass supply rate ${\dot{M}}_{0}$ to within 2% and where the time-averaged mass accretion profile in the CBD $\langle {\dot{M}}_{{\rm{d}}}\rangle (R)$ is also $\approx {\dot{M}}_{0}$ (i.e., flat to within an rms tolerance of 5%). Unless stated otherwise, all simulations presented henceforth include a preliminary phase ($3000{P}_{{\rm{b}}}$) in which the disk is evolved using a diode boundary condition at ${R}_{\mathrm{in}}={a}_{{\rm{b}}}(1+{e}_{{\rm{b}}})$. If a flat $\langle {\dot{M}}_{{\rm{d}}}\rangle$ profile is reached within that time period, the boundary is removed and the simulation is resumed in the self-consistent fashion described in the preceding paragraphs.

The simulation output is processed to produce time series for ${\dot{M}}_{{\rm{b}}}$ and ${\dot{q}}_{{\rm{b}}}$, the specific torques ${{dl}}_{{\rm{b}}}/{dt}{| }_{\mathrm{grav}}$ and ${{dl}}_{{\rm{b}}}/{dt}{| }_{\mathrm{acc}}$, and the spin change rate ${\dot{S}}_{{\rm{b}}}\equiv {\dot{S}}_{1}+{\dot{S}}_{2}$ (see Section 2.2.1). These quantities and the total ${\dot{J}}_{{\rm{b}}}$ are shown in Figure 2 for a time interval of $200{P}_{{\rm{b}}}$. We find that in quasi-steady state, the time-averaged ${\dot{J}}_{{\rm{b}}}$ is given by

$$\begin{eqnarray}&&\langle {\dot{J}}_{{\rm{b}}}\rangle \approx 0.676{\dot{M}}_{0}{{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}.\end{eqnarray} \tag{ 27 }$$

Of the five different contributions to ${\dot{J}}_{{\rm{b}}}$, two quantities dominate over the rest: the mass change $\langle {\dot{M}}_{{\rm{b}}}\rangle \approx {\dot{M}}_{0}$ (first panel) and the specific gravitational torque $\langle {\dot{l}}_{{\rm{b}},\mathrm{grav}}\rangle \,\approx 1.567({\dot{M}}_{0}/{M}_{{\rm{b}}}){{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}$ (third panel).

Zoom In Zoom Out Reset image size

The stationary nature of ${\dot{M}}_{{\rm{b}}}$ and ${\dot{J}}_{{\rm{b}}}$ means that, in an averaged sense, a constant amount of angular momentum is transferred to the binary by the disk per unit accreted mass. This constant is the eigenvalue of the accretion flow (e.g., Paczynski 1991; Popham & Narayan 1991), that is,

$$\begin{eqnarray}&&\displaystyle \frac{\langle {\dot{J}}_{{\rm{b}}}\rangle }{\langle {\dot{M}}_{{\rm{b}}}\rangle }\equiv {l}_{0}.\end{eqnarray} \tag{ 28 }$$

Our result in Figure 2 implies

$$\begin{eqnarray}&&{l}_{0}\left|{}_{{e}_{{\rm{b}}}=0}\approx 0.68{{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}\right..\end{eqnarray} \tag{ 29 }$$

The eigenvalue l₀ can also be obtained by computing $\langle {\dot{J}}_{{\rm{d}}}\rangle (R)$, the angular momentum current as a function of radius in the CBD (Equation (15)). We compute ${\dot{J}}_{{\rm{d}}}(R)$ from the combination of ${\dot{J}}_{{\rm{d}},\mathrm{adv}}$, ${\dot{J}}_{{\rm{d}},\mathrm{visc}}$, and ${\dot{J}}_{{\rm{d}},\mathrm{grav}}$, in turn computed from the time-averaged angular momentum flux maps F_J as described in Section 2.2.2. We compute these angular momentum flux maps from simulation snapshots produced at a high output rate (15 snapshots per orbit) followed by time-averaging. Figure 3 shows the result. Importantly, we see that the net angular momentum current $\langle {\dot{J}}_{{\rm{d}}}\rangle (R)$ (dark red) is nearly flat for R between ${a}_{{\rm{b}}}$ and $10{a}_{{\rm{b}}}$. Furthermore, $\langle {\dot{J}}_{{\rm{d}}}\rangle (R)$ is very close to $\langle {\dot{J}}_{{\rm{b}}}\rangle$. In other words, all of the angular momentum crossing the CBD is received by the binary. This provides convincing evidence that the exchange of mass and angular momentum between the CBD and the binary are in a true quasi-steady state, and that our computed eigenvalue (Equation (29)) is reliable.

Zoom In Zoom Out Reset image size

Dependence on sink radius: In our simulations, the "stars" behave as sink particles. In most of our simulations, we adopt the "stellar radius" of ${r}_{\mathrm{acc}}=0.02{a}_{{\rm{b}}}$, but it is natural to evaluate how sensitive our results are to the choice of ${r}_{\mathrm{acc}}$. Recently, Tang et al. (2017) posited that the sign of the net angular momentum transfer rate $\langle {\dot{J}}_{{\rm{b}}}\rangle$ (note that Tang et al. 2017 only computed $\langle {\dot{L}}_{{\rm{b}}}\rangle$ and did not include spin torques) can depend critically on the way gas is drained from the computational domain by the stars (see Section 5.2.1). We test the influence of the accretion algorithm on our results by carrying out two additional simulations with ${r}_{\mathrm{acc}}=0.04{a}_{{\rm{b}}}$ and ${r}_{\mathrm{acc}}=0.06{a}_{{\rm{b}}}$ (the softening lengths are updated accordingly, but all other parameters remain unchanged). The results of this test are depicted in Figure 4, where we contrast the gas density field (to illustrate the effect of a larger sink radius) with the evolution of ${\dot{M}}_{{\rm{b}}}$, ${\dot{S}}_{{\rm{b}}}$, and ${\dot{J}}_{{\rm{b}}}$. The mean values $\langle {\dot{J}}_{{\rm{b}}}\rangle$ and $\langle {\dot{M}}_{{\rm{b}}}\rangle$ are robust against ${r}_{\mathrm{acc}}$, but $\langle {\dot{S}}_{{\rm{b}}}\rangle$ grows almost proportionally to ${r}_{\mathrm{acc}}$. Since $\langle {\dot{J}}_{{\rm{b}}}\rangle$ appears unmodified, the increase in $\langle {\dot{S}}_{{\rm{b}}}\rangle$ must be accompanied by a decrease in the torque ${\mu }_{{\rm{b}}}\langle {\dot{l}}_{{\rm{b}},\mathrm{grav}}\rangle$ (not shown). Indeed, over the interval $[3500,3700]{P}_{{\rm{b}}}$, the specific gravitational torque $\langle {\dot{l}}_{{\rm{b}},\mathrm{grav}}\rangle$ is $1.534({\dot{M}}_{0}/{M}_{{\rm{b}}}){{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}$ when ${r}_{\mathrm{acc}}=0.04{a}_{{\rm{b}}}$ and $1.482({\dot{M}}_{0}/{M}_{{\rm{b}}}){{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}$ when ${r}_{\mathrm{acc}}=0.06{a}_{{\rm{b}}}$, both smaller than the value of $1.567({\dot{M}}_{0}/{M}_{{\rm{b}}}){{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}$ in our fiducial simulation (Figure 2). The anisotropic accretion torque ${\mu }_{{\rm{b}}}\langle {\dot{l}}_{{\rm{b}},\mathrm{acc}}\rangle$ remains negligible in all cases. For all the values of ${r}_{\mathrm{acc}}$ explored, $\langle {\dot{S}}_{{\rm{b}}}\rangle$ is a small contribution to the total transfer of angular momentum to the binary $\langle {\dot{J}}_{{\rm{b}}}\rangle$. Presumably, a much larger value of ${r}_{\mathrm{acc}}$ than explored here—one that forbids the formation of CSDs—could make $\langle {\dot{S}}_{{\rm{b}}}\rangle$ an important contributor to $\langle {\dot{J}}_{{\rm{b}}}\rangle$, one that competes with ${\mu }_{{\rm{b}}}{\dot{l}}_{{\rm{b}},\mathrm{grav}}$, perhaps to the point of erasing the positive torque contribution of the CSDs, turning all of the positive gain in orbital angular momentum into a growth in spin.

Zoom In Zoom Out Reset image size

A major result of this paper (see Figures 2–3) is that $\langle {\dot{J}}_{{\rm{b}}}\rangle \simeq \langle {\dot{J}}_{{\rm{d}}}\rangle \gt 0$; that is, an accreting binary gains angular momentum from a steady-supply disk. From Figure 2, we see that a major contribution to $\langle {\dot{J}}_{{\rm{b}}}\rangle$ comes from the total gravitational torque ${\mu }_{{\rm{b}}}\langle {\dot{l}}_{{\rm{b}},\mathrm{grav}}\rangle$ (third panel), which is positive, in apparent contradiction with what one might expect from theoretical work (Goldreich & Tremaine 1980; Artymowicz & Lubow 1994; Lubow & Artymowicz 1996), since the linear theory of Lindblad torques of Goldreich & Tremaine (1979) predicts that an outer disk exerts negative torques on the inner binary. Note that our result for the angular momentum current ${\dot{J}}_{{\rm{d}},\mathrm{grav}}$ in the CBD (Figure 3, yellow curve) is indeed consistent with this theoretical expectation. However, a self-consistent calculation of $\langle {\dot{J}}_{{\rm{b}}}\rangle$ needs to include the entire distribution of gas and not only that of the CBD.

Figure 5 depicts the dissection of the total gravitational torque ${T}_{{\rm{b}},\mathrm{grav}}={\mu }_{{\rm{b}}}{\dot{l}}_{{\rm{b}},\mathrm{grav}}$ into three components (from top to bottom): (1) the outer disk torque ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{out})}$ ($R\gt {R}_{\mathrm{cav}};$ with ${R}_{\mathrm{cav}}$ as defined in Equation (22)), (2) the "cavity+streams" torque ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{cav})}$ (${a}_{{\rm{b}}}\leqslant R\leqslant {R}_{\mathrm{cav}}$), and (3) the "inner" torque ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{in})}$ ($R\lt {a}_{{\rm{b}}}$). The short-term evolution of these torques is depicted on the left panels of the figure, and the long-term evolution is shown on the right. On average, the sum of these three different torques must equal $\langle {T}_{{\rm{b}},\mathrm{grav}}\rangle ={\mu }_{{\rm{b}}}\langle {\dot{l}}_{{\rm{b}},\mathrm{grav}}\rangle =0.392{\dot{M}}_{0}{a}_{{\rm{b}}}^{2}{{\rm{\Omega }}}_{{\rm{b}}}$.

Zoom In Zoom Out Reset image size

From the torque dissection, we confirm that the gravitational torque acting on the binary due to the material outside $R={a}_{{\rm{b}}}$ is indeed negative, but that the material inside $R={a}_{{\rm{b}}}$ exerts a positive torque that nearly doubles in magnitude the negative outer torque. We also find that the oscillation amplitudes of these three torque contributions can be much larger than their mean values. In order to capture both the rapid sign changes in ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{out})}$, ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{cav})}$, and ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{in})}$ and their long-term averages, a fine sampling in time as well as a long-term integration are required. We find that a sample rate of 15 snapshots per binary orbit captures ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{out})}$ and ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{cav})}$ adequately because they vary on timescales of $\sim {P}_{{\rm{b}}};$ however, because ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{in})}$ can vary on timescales comparable to the orbital period within the CSDs, we need a sampling rate of ∼100 snapshots per orbit to capture such rapid variability. This rapid variability is difficult to capture from a postprocessing analysis of simulation snapshots, but is automatically accounted for by the strategy described in Section 2.2.1.

The contribution of the cavity region to the gravitational torque can be better understood if we explore the spatial distribution of the torque in two dimensions. For this, we construct gravitational torque density maps. Evidently, since $-{{\boldsymbol{T}}}_{{\rm{b}},\mathrm{grav}}={{\boldsymbol{T}}}_{{\rm{d}},\mathrm{grav}}$ (the torque exerted on the disk by the binary), the torque per unit area can be computed directly from the cell-centered values of density and gravitational acceleration:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\boldsymbol{T}}}_{{\rm{b}},\mathrm{grav}}}{d{\rm{A}}}({\boldsymbol{r}}) & = & -\displaystyle \frac{d{{\boldsymbol{T}}}_{{\rm{d}},\mathrm{grav}}}{d{\rm{A}}}({\boldsymbol{r}})={\rm{\Sigma }}({\boldsymbol{r}}){\boldsymbol{r}}\times {\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{b}}}({\boldsymbol{r}})\\ & \approx & -{{\rm{\Sigma }}}_{k}{{\boldsymbol{r}}}_{k}\times {{\boldsymbol{a}}}_{k},\end{array}\end{eqnarray} \tag{ 30 }$$

where ${{\boldsymbol{a}}}_{k}$ is the gravitational acceleration of the gas and the subscript k denotes evaluation at the center of a cell.

To explore the long-term effect of the nonaxisymmetric features in the circumbinary gas, we can take time averages of both the density and torque density fields. This is shown in Figure 6, where we have stacked the density field Σ and the torque per unit area (Equation (30); 15 snapshots per orbit) rotated by an angle that matches the rotation of the binary. The coordinate axes in this figure are thus ξ and η, related to the original coordinates by the area-preserving transformation $x={a}_{{\rm{b}}}(\xi \cos {f}_{{\rm{b}}}-\eta \sin {f}_{{\rm{b}}})$ and $y={a}_{{\rm{b}}}(\xi \sin {f}_{{\rm{b}}}+\eta \cos {f}_{{\rm{b}}})$.⁶ We carry out this average over intervals of ∼30 orbits in duration, at different times in the simulation. In addition, we perform the average over 200 orbits (rightmost panels of Figure 6). All of the averaging sets produce essentially the same result, implying that most of the variability takes place on timescales shorter than 30 orbits. Figure 6 reveals that the streamers are the dominant nonaxisymmetric feature that persists in time. From the torque density maps (bottom panels), it can be seen that different portions of the gas can "torque up" (regions in red) or "torque down" (regions in blue) the binary, but only the persistent nonaxisymmetric features contribute to the nonzero net torque.

Zoom In Zoom Out Reset image size

In Figure 7 we show the total angular momentum transfer rate ${\dot{J}}_{{\rm{b}}}$ and its various contributions for a ${q}_{{\rm{b}}}=1$, ${e}_{{\rm{b}}}=0.6$ binary. As expected, we find that ${\dot{M}}_{{\rm{b}}}$ is no longer modulated at a frequency of $\sim \tfrac{1}{5}{{\rm{\Omega }}}_{{\rm{b}}}$, but that the dominant frequency $\sim {{\rm{\Omega }}}_{{\rm{b}}}$ (ML16). There are variations over much longer timescales: the mass ratio change ${\dot{q}}_{{\rm{b}}}$ flips signs every $\sim 415\,{P}_{{\rm{b}}}$. The gravitational specific torque ${\dot{l}}_{{\rm{b}},\mathrm{grav}}$ is again positive and larger than in the circular case. The contribution of the anisotropic accretion torque ${\mu }_{{\rm{b}}}\langle {\dot{l}}_{{\rm{b}},\mathrm{acc}}\rangle \approx -0.21{\dot{M}}_{0}{{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}$ is no longer negligible. As in the circular case, the contribution of the spin torque $\langle {\dot{S}}_{{\rm{b}}}\rangle$ to the total torque $\langle {\dot{J}}_{{\rm{b}}}\rangle$ is small ($\sim 2 \%$).

Zoom In Zoom Out Reset image size

As in Section 3.1, we can compute the accretion eigenvalue l₀ directly from the result shown in Figure 7, giving

$$\begin{eqnarray}&&{l}_{0}\left|{}_{{e}_{{\rm{b}}}=0.6}\approx 0.81{{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}\right..\end{eqnarray} \tag{ 31 }$$

We also compute the net angular momentum current in the CBD ${\dot{J}}_{{\rm{d}}}(R)$ (Section 2.2.2) averaged over 300 orbits, as shown in Figure 8. We see that $\langle {\dot{J}}_{{\rm{d}}}\rangle (R)$ is approximately constant (independent of R) and equal to $\langle {\dot{J}}_{{\rm{b}}}\rangle$, indicating that the global angular momentum balance is achieved. In this case, the match between $\langle {\dot{J}}_{{\rm{d}}}\rangle /{\dot{M}}_{0}$ and ${l}_{0}=0.81{{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}$ is remarkable for $R\gt 4{a}_{{\rm{b}}}$ and $R\lt 2{a}_{{\rm{b}}}$.

Zoom In Zoom Out Reset image size

As in the case of the circular binary, the eccentric binary gains angular momentum, which, as before, is largely a result of the positive net gravitational torque $\langle {T}_{{\rm{b}},\mathrm{grav}}\rangle ={\mu }_{{\rm{b}}}\langle {\dot{l}}_{{\rm{b}},\mathrm{grav}}\rangle \,\approx 0.796{\dot{M}}_{0}{{\rm{\Omega }}}_{{\rm{b}}}{a}_{{\rm{b}}}^{2}$. As before, we study the spatial distribution of gravitational torques in the simulation by splitting ${T}_{{\rm{b}},\mathrm{grav}}$ into three components: ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{out})}$, ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{cav})}$, and ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{in})}$ (Figure 9). We see a behavior similar to that of a circular binary (Figure 5): the sum ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{out})}+{T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{cav})}$ is negative (dominated by ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{cav})}$), while the positive sign of the total torque $\langle {T}_{{\rm{b}},\mathrm{grav}}\rangle$ is entirely due to the positive sign of ${T}_{{\rm{b}},\mathrm{grav}}^{(\mathrm{in})}$. The torque time series shown in Figure 9 exhibit some degree of regularity, but they differ from the nearly exact periodicity seen in Figure 5. As noted above, the accretion flow onto eccentric binaries is modulated over long secular timescales, and it is not surprising that a few tens of orbits cannot capture the entire range of time variability.

Zoom In Zoom Out Reset image size

The fact that tens of orbits is too short a time span to reveal the underlying stationary behavior of the accretion flow can be illustrated by repeating the map-stacking analysis for the ${e}_{{\rm{b}}}=0.6$ case (see Figure 10 and Figure 6). The stacking of these images is trickier than in the circular case: in order for the two "stars" to line up at all snapshots, we need to introduce a rotating–pulsating coordinate system, a transformation that is well known in the study of the elliptical restricted three-body problem (ER3BP; e.g., Nechvíle 1926; Kopal & Lyttleton 1963; Szebehely & Giacaglia 1964; Szebehely 1967; Musielak & Quarles 2014). The scaled, rotated coordinates ξ and η are related to the barycentric coordinates x and y by $x={r}_{{\rm{b}}}(\xi \cos {f}_{{\rm{b}}}-\eta \sin {f}_{{\rm{b}}})$ and $y={r}_{{\rm{b}}}(\xi \sin {f}_{{\rm{b}}}+\eta \cos {f}_{{\rm{b}}})$, where ${r}_{{\rm{b}}}={a}_{{\rm{b}}}(1-{e}_{{\rm{b}}}^{2})/(1+{e}_{{\rm{b}}}\cos {f}_{{\rm{b}}})$ and ${f}_{{\rm{b}}}$ is the true anomaly of the binary. This transformation is not area preserving, so scalars must be transformed by multiplying by the Jacobian of the transformation ($={r}_{{\rm{b}}}^{2}$). In Figure 10 (top row, left to right), we show Σ stacked over the time intervals $[3470\,{P}_{{\rm{b}}},3500\,{P}_{{\rm{b}}}]$, $[3530\,{P}_{{\rm{b}}},3560\,{P}_{{\rm{b}}}]$, $[3590\,{P}_{{\rm{b}}},3620\,{P}_{{\rm{b}}}]$, and $[3470\,{P}_{{\rm{b}}},3799\,{P}_{{\rm{b}}}]$. These reveal a striking inherent property of accretion onto eccentric binaries: the "loading" of the CSDs alternates in time; that is, the gas surface density is higher in one CSD than in the other. This disk loading can be directly tied to the sign of ${\dot{q}}_{{\rm{b}}}$ in Figure 7. When we average over $370{P}_{{\rm{b}}}$—almost the entire precessional period identified from the time series of ${\dot{q}}_{{\rm{b}}}$—we see that most of the asymmetries disappear (the last column of Figure 10), and that the distribution of gas starts to resemble that of the circular case (Figure 6). Note how the shape of the cavity also changes as preferential accretion switches recipients, as is expected if alternation is caused by the precession of the cavity edge.

Zoom In Zoom Out Reset image size

We use this rotating system to study the spatial distribution of torques, albeit only in a qualitative manner (Figure 10). Some similarities can be found with the analogous plot for the circular case (Figure 6, bottom row), although the spatial distribution of torques in the eccentric case is more diffuse and variable. A direct comparison to the spatial torque dissection of Figure 9 is not possible, since the barycentric radial coordinate R changes from snapshot to snapshot in the rotating–pulsating coordinate system.

The most important result of our study is that in quasi-steady state the binary gains angular momentum. This result was also strongly suggested in our recent study (MML17) using the polar-grid code PLUTO, in which the innermost region of the circumbinary cavity is excised from the simulation domain. MML17 showed that a quasi-steady state can be reached in the CBD, with global balance of both $\langle {\dot{M}}_{{\rm{d}}}\rangle (R)$ and $\langle {\dot{J}}_{{\rm{d}}}\rangle (R)$. In fact, for ${e}_{{\rm{b}}}=0$, the accretion eigenvalue obtained by MML17 (${l}_{0}\approx 0.8{a}_{{\rm{b}}}^{2}{{\rm{\Omega }}}_{{\rm{b}}}$ ) is close to the value obtained in this paper (${l}_{0}\simeq 0.68{a}_{{\rm{b}}}^{2}{{\rm{\Omega }}}_{{\rm{b}}};$ Equation (29)). At high eccentricities, however, the mismatch between our new results and those of MML17 grows (e.g., they found ${l}_{0}\gtrsim {a}_{{\rm{b}}}^{2}{{\rm{\Omega }}}_{{\rm{b}}}$ for ${e}_{{\rm{b}}}=0.6$). This discrepancy may be understood from the fact that, for finite eccentricities, the mass flow ${\dot{M}}_{{\rm{d}}}(R)$ at $R={a}_{{\rm{b}}}(1+{e}_{{\rm{b}}})$ (the inner boundary of the computation domain adopted by MML17) can be both inward and outward (see Figure 2 of ML16), a behavior that a diode-like open boundary (adopted by MML17) cannot capture. The mismatch also underscores the importance of resolving the flows inside the cavity. Because the simulations of MML17 did not fully include the inner cavity, one should take those results as suggestive. In this paper, calculating $\langle {\dot{J}}_{{\rm{b}}}\rangle$ in two different ways and obtaining agreement (direct torque and angular momentum current across the CBD), we have obtained robust evidence that the binary can gain angular momentum while accreting ($\langle {\dot{J}}_{{\rm{b}}}\rangle \gt 0$). This result is insensitive to the accretion radius (the size of the "stars") as long as ${r}_{\mathrm{acc}}\lesssim 0.1{a}_{{\rm{b}}}$ (see Table 1).

Our result ($\langle {\dot{J}}_{{\rm{b}}}\rangle \gt 0$) contradicts the widely held view that a binary loses angular momentum to its surrounding disk and experiences orbital decay. In fact, while there have been many numerical simulations of circumbinary accretion (see references in Section 1), only a few address the issue of angular momentum transfer in a quantitative way. Three such studies use simulations that excise the inner cavity (MacFadyen & Milosavljević 2008; Shi et al. 2012; MML17). The work of MacFadyen & Milosavljević (2008) was the most similar to MML17: these authors 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}_{{\rm{b}}}$ and ${R}_{\mathrm{out}}=100{a}_{{\rm{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 viscosity parameter (by contrast, MML17 used $\alpha =0.1$ and 0.05), the "viscous relaxation" radius $t=4000{P}_{{\rm{b}}}$ (the typical duration of their runs) is only about $3{a}_{{\rm{b}}}$, and their surface density profile is far from steady state even for $R\gg {R}_{\mathrm{in}}$. We suggest that the result of MacFadyen & Milosavljević (2008) reflects a "transient" phase of their simulations. Shi et al. (2012) obtained a positive value of ${\dot{J}}_{{\rm{b}}}$ in their 3D MHD simulations of CBD (truncated at ${R}_{\mathrm{in}}=0.8{a}_{{\rm{b}}}$). However, the duration of their simulations is only $\sim 100{P}_{{\rm{b}}}$, and it is unlikely that a quasi-steady state is reached. Their value of l₀, which is too small to cause orbital expansion, may not properly characterize the long-term evolution of the binary. As noted above, our previous study (MML17) using the polar-grid code PLUTO (with ${R}_{\mathrm{in}}={a}_{{\rm{b}}}(1+{e}_{{\rm{b}}})$ for eccentric binaries) did achieve quasi-steady state in terms of mass and angular momentum transport in the CBD. For ${e}_{{\rm{b}}}=0$, the value of l₀ measured by MML17 was close to the value presented in this paper (Table 1), while the l₀ values for ${e}_{{\rm{b}}}\gt 0$ were unreliable for understandable reasons (see above).

In a recent paper, Tang et al. (2017here TMH17) carried out hydrodynamical simulations of accretion onto circular binaries using the DISCO code (Duffell & MacFadyen 2012; Duffell 2016). Their simulations are similar to ours in many respects (i.e., locally isothermal equation of state with $H/R=0.1$, $\alpha =0.1$). Their accretion prescription is different from ours: inside a "sink" radius (measured from each "star"), they assumed that the gas is depleted at the rate $d{\rm{\Sigma }}/{dt}=-{\rm{\Sigma }}/{t}_{\mathrm{sink}}$, with ${r}_{\mathrm{sink}}=0.1{a}_{{\rm{b}}}$ and ${t}_{\mathrm{sink}}$ a free parameter. By contrast, the AREPO code in its Navier–Stokes formulation (Muñoz et al. 2013) self-consistently incorporates viscous transport in all of our computation domain in a unified manner, and gas accretion onto the "stars" is treated using an absorbing sink algorithm with a single parameter ${r}_{\mathrm{acc}}=0.02{a}_{{\rm{b}}}$ (see Section 2.1). Using a direct force computation, TMH17 found that the net torque on the binary is negative unless ${t}_{\mathrm{sink}}$ is much smaller than ${P}_{{\rm{b}}}$. The discrepancy between our results and those of TMH17 is important and deserves further study and comparison. For now, we offer the following comments: (1) In our simulations, we have achieved global mass transfer balance ($\langle {\dot{M}}_{1}\rangle +\langle {\dot{M}}_{2}\rangle \,\simeq \langle {\dot{M}}_{{\rm{d}}}\rangle (R)\simeq {\dot{M}}_{0};$ see Figures 1 and 2) and global angular momentum transfer balance ($\langle {\dot{J}}_{{\rm{b}}}\rangle \simeq \langle {\dot{J}}_{{\rm{d}}}\rangle (R)\simeq \mathrm{constant};$ see Figures 3 and 8); TMH17 did not demonstrate that they have achieved these. (2) The mass depletion prescription adopted by TMH17 may be problematic in terms of mass conservation. TMH17 found that the more "filled" the CSDs are (corresponding to larger ${t}_{\mathrm{sink}}$), the more negative the net torque on the binary is (see their Figure 2); this is counterintuitive, as on physical grounds one expects that CSDs contribute a positive gravitational torque on the binary (Savonije et al. 1994). (3) It is unclear whether the definition and calculation of the accretion torque in TMH17 are consistent with angular momentum conservation (see their Equation (8)).

Turning to eccentric binaries, our paper represents the first work to compute the secular evolution rates of ${a}_{{\rm{b}}}$ and ${e}_{{\rm{b}}}$, consistently taking into account accretion via direct, long-term simulations. Most previous (analytic) works (such as Artymowicz & Lubow 1996; Lubow & Artymowicz 2000) assumed that accretion was negligible (${\dot{M}}_{{\rm{b}}}=0$) and that ${\dot{a}}_{{\rm{b}}}$ and ${\dot{e}}_{{\rm{b}}}$ were solely due to gravitational Lindblad torques. Hayasaki (2009) carried out an analytic calculation including the effect of accretion coupled to the Lindblad torques from the CBD; in contrast to our new findings, he concluded that, for $\alpha ={h}_{0}=0.1$, binaries migrate inward if they accrete at the same rate as the CBD.

In this paper, we have "solved" the problem of circumbinary accretion in the "cleanest" scenario and in quasi-steady state, taking advantage of AREPO's capability to resolve a wide dynamical range in space (from $70{a}_{{\rm{b}}}$ down to $0.02{a}_{{\rm{b}}}$). We should, however, bear in mind the limitations of our work when confronting real astrophysical applications.

Our simulations are two-dimensional and adopt a simple equation of state as well as the α-viscosity prescription. Three-dimensional MHD simulations (Noble et al. 2012; Shi et al. 2012; Shi & Krolik 2015) would represent the transport of angular momentum more realistically. In addition, general relativistic effects are bound to be important at close separations in the case of supermassive binary black holes (SMBBHs; Farris et al. 2012).

The value of α chosen (0.1) in our simulations is in part due to computational expediency. A larger viscosity results in a shorter viscous time—still much longer than all other relevant timescales—which allows us to reach the quasi-steady state in the inner region ($R\lesssim 10{a}_{{\rm{b}}}$) of the CBD. Very small values of α would prevent us from attaining such a steady state. MML17 considered $\alpha =0.1$ and 0.05 in their PLUTO simulations and found the same l₀ value for circular binaries. At very low viscosities, one can expect more pronounced spiral arms (in the CBD and CSDs), which would increase the tidal torques; at the same time, the tidal cavity surrounding the binary becomes larger and the CSDs are truncated at smaller radii (Artymowicz & Lubow 1994; Miranda & Lai 2015), thus removing the gas that would exert the strongest torques. These opposing effects may cancel each other out (at least for $\alpha =0.05-0.1$ and ${e}_{{\rm{b}}}=0$), as suggested by the results of MML17. On the other hand, the smoothed particle hydrodynamics simulations by Ragusa et al. (2016) suggest that accretion may be suppressed when the viscosity is sufficiently small (small α or small h₀), although it is not clear if those simulations were run for a sufficiently long integration time. Exploring the dependence of our results on α and h₀ will be the subject of future work.

We have assumed a steady mass supply at a large distance, which may not exist under realistic astrophysical conditions. We expect that our general results should hold even if the external supply rate ${\dot{M}}_{0}$ varies in time, provided it does so smoothly and slowly, on a timescale longer than the accretion time at the outer edge of the CBD, ${\tau }_{\mathrm{acc}}\sim {R}_{\mathrm{out}}^{2}/\nu$. As several of the quantities measured in this work are scaled by ${\dot{M}}_{0}$, such as ${\dot{J}}_{{\rm{b}}}$, ${\dot{a}}_{{\rm{b}}}$, and ${\dot{e}}_{{\rm{b}}}$, a slow change in the supply rate will translate into a proportional amplification or reduction of these quantities, without affecting their qualitative behavior. On the other hand, clear violations of the steady-supply condition—such as accretion via gas clumps (e.g., Goicovic et al. 2016) or from tori of small radial extent—might change the response of the binary.

We have focused on ${q}_{{\rm{b}}}=1$ in this paper, but the behavior of ${q}_{{\rm{b}}}\lt 1$ binaries may be different. In particular, the phenomenon of alternating accretion described in ML16, attributed to apsidal precession in the CBD, has not been thoroughly explored for mass ratios different from unity. The value of $\langle {\dot{q}}_{{\rm{b}}}\rangle$—and whether it averages to zero over secular timescales—is very important for properly determining $\langle {\dot{J}}_{{\rm{b}}}\rangle$ when ${q}_{{\rm{b}}}\ne 1$ (Equations (2) and (4)). The sign of $\langle {\dot{q}}_{{\rm{b}}}\rangle$ from CBD accretion is still an unresolved question, and recent work by Young & Clarke (2015) suggests that it depends on the aspect ratio of the disk (see also Clarke 2008). No grid-based simulation has explored the sign of $\langle {\dot{q}}_{{\rm{b}}}\rangle$ while systematically varying both ${q}_{{\rm{b}}}$ and ${e}_{{\rm{b}}}$, and this will be the subject of future work.

Finally, we do not include self-gravity of the disk in our simulations, as we have implicitly assumed that the local disk mass, ${R}^{2}{\rm{\Sigma }}(R)$ (with $R\sim 5{a}_{{\rm{b}}}$), is much less than the binary mass. A massive, self-gravitating disk may qualitatively change the dynamics of circumbinary accretion studied in this paper (e.g., Cuadra et al. 2009; Roedig et al. 2012), although some disk fragmentation studies have already reported outward migration of protobinaries (Kratter et al. 2010a); see below.

Supermassive binary BHs: It is widely assumed that CBDs around SMBBHs extract angular momentum from such binaries, shrinking their orbital radii from ∼1 pc to ∼0.01 pc, before gravitational radiation can push the binary toward merger within a Hubble time (e.g., Begelman et al. 1980; Armitage & Natarajan 2002; MacFadyen & Milosavljević 2008; Rafikov 2016; Tang et al. 2017). Our findings contradict this assumption, thus putting into question the relevance of circumbinary gas as a solution to the "final parsec problem." At least for equal-mass binary black holes and non-self-gravitating disks, our results suggest that circumbinary accretion may prevent rather than assist the merger of such objects.

Formation of stellar binaries: The finding that binaries embedded in accretion disks are able to gain angular momentum poses problems to binary formation theory. Of the proposed mechanisms for binary formation at intermediate separations ($\lesssim 10\,\mathrm{au};$ e.g., Tohline 2002) the still-qualitative "fragmentation then migration" scenario has emerged as a strong contender over the past decade (e.g., Krumholz et al. 2009; Kratter et al. 2010a; Zhu et al. 2012). This mechanism is grounded in the difficulty for massive disks to fragment at distances smaller than 50–100 au from the central object (Matzner & Levin 2005; Rafikov 2005; Boley et al. 2006; Whitworth & Stamatellos 2006; Stamatellos & Whitworth 2008; Cossins et al. 2009; Kratter et al. 2010b) and thus requires binaries to reduce their separations by one or two orders of magnitude after they have formed (see Moe & Kratter 2018 for a recent analysis of binaries at small separations and their possible formation mechanisms). If binaries embedded in disks tend to expand, the viability of this mechanism would be questionable (e.g., Satsuka et al. 2017).

How can binaries still shrink?: Given the idealized nature of our simulations, we cannot rule out the possibility of binary shrinkage. We envision several possibilities where binary decay may occur:

• (i)  
  Unequal-mass binaries and locally massive disks. We have only explored binaries with ${q}_{{\rm{b}}}=1$ and small disk masses, but systems with ${q}_{{\rm{b}}}\lt 1$ and disks with local mass $\sim {\rm{\Sigma }}{a}_{{\rm{b}}}^{2}\gtrsim {M}_{2}$ could behave very differently (Cuadra et al. 2009; Roedig et al. 2012). In particular, when ${q}_{{\rm{b}}}\ll 1$, the binary would evolve in a fashion similar to Type II migration (e.g., Lin & Papaloizou 1986; Nelson et al. 2000) or a modified version of it (e.g., Duffell et al. 2014; Dürmann & Kley 2015). Such massive disks may exist in the early stage immediately following galaxy merger (producing a massive gas torus surrounding an SMBH binary) or the early (class zero) stage of stellar binary formation.
• (ii)  
  Three-dimensional angular momentum loss effects. Our 2D calculations cannot capture effects that are intrinsically 3D, such as winds and outflows. MHD winds can take away angular momentum in accretion disks, and such a mechanism may operate in circumbinary accretion as well. The 3D ideal MHD simulations of Kuruwita et al. (2017) suggest that newly formed binaries could contract significantly within hundreds of orbits as they grow in mass by a factor of ∼100, again implying that binary shrinkage may take place at a much earlier phase of binary formation than that explored in the present work.