THREE-DIMENSIONAL MAGNETOHYDRODYNAMIC SIMULATIONS OF CIRCUMBINARY ACCRETION DISKS: DISK STRUCTURES AND ANGULAR MOMENTUM TRANSPORT

We construct our physical disk model in the inertial frame in which the center of mass of the binary is at rest at the origin. Treating the simplest case first, we assume an equal-mass circular restricted binary system, i.e., the binary orbits circularly in the disk plane, and we neglect binary evolution as the disk inflow timescale is usually much smaller than the shrinkage timescale of the binary, and certainly more than the duration of the simulation (see estimate in Section 5). We set the gravitational constant G, total binary mass M, and binary separation a to be unity, and therefore the binary frequency $\Omega _{\rm bin}=\sqrt{GM/a^3}$ is unity as well. We also assume the circumbinary disk to be cold and thin. As we are mainly concerned with the orbital dynamics of the flow, we choose a simple global isothermal sound speed c_s = 0.05. The disk flares at larger radii because the ratio of height to radius H/R = c_s/RΩ_K = 0.05(R/a)^1/2, where H denotes the scale height of the disk, R is the disk radius in cylindrical coordinates, and $\Omega _{\rm K}=\sqrt{GM/R^3}$ is the Keplerian frequency. As the fluid is well coupled to the magnetic field even in a cold disk where the ionization fraction is far below unity (see, e.g., Stone et al. 2000), we assume ideal MHD. For this reason, we include no explicit diffusivity except the von Neumann–Richtmyer bulk viscosity in compressive regions that ensures the right jump conditions for shock waves. The effect on damping the angular momentum is negligible compared to other transport mechanisms. Because we are most concerned with the inner part of the disk, not accretion onto the binary, we excise a central region. This cutout must be well inside the inner edge of the disk so that it does not affect fully resolving the inner disk and any gas leakage from the inner edge. We choose to cut out the area within 0.8a. This region is beyond the main extent of the interior disks that surround each binary member because these disks are each tidally truncated at ∼0.3a from their central objects (Paczyński 1977). Once the disk reaches the quasi-steady stage, we find that the inner edge of the disk is located at r ≃ 2a, which is about a factor of 2.5 outside the cutout. We approximate the time-dependent binary potential in the disk region using Newtonian dynamics because the disk region we are interested in is far from the gravitational radii of the individual binary components. The disk mass is assumed to be much smaller than the binary mass, so that the Toomre parameter satisfies Q ≈ (H/R)(M/M_d) ≫ 1, where M_d denotes the disk mass. We therefore neglect the contribution of the disk self-gravity to the potential.

The properties of circumbinary disks require us to resolve three different length scales when setting up the grid. The first is the half-disk thickness H. The computational domain has to contain at least several (≲ 4) scale heights on each side of the midplane, and for each H, several tens of cells are needed. The second is the maximum growth-rate wavelength of the MRI, $\lambda _{\rm MRI}\equiv 8\pi /\sqrt{15}v_A/\Omega (R)$, where v_A is the Alfvén speed and Ω(R) is the disk rotational frequency. We require λ_MRI to be resolved by at least six grid elements. The last one is the spiral density wavelength, λ_d ∼ 2πc_s/Ω_bin. There should be many cells across this wavelength (MM08); we require at least six.

In order to satisfy the above requirements, we follow the scheme proposed in Noble et al. (2010) to construct our grid in spherical coordinates (r, θ, ϕ). We adopt a logarithmic grid in the radial direction, which provides a constant Δr/r. The vertical grid is derived by mapping a simple linear function y(x) = x for x ∈ [0, 1] through a polynomial transformation (see Equation (6) in Noble et al. 2010):

$$\begin{equation} \theta (y) = \frac{\pi }{2}\left[1+(1-\xi)(2y-1)+\left(\xi -\frac{2\theta _c}{\pi }\right)(2y-1)^n\right],\vspace*{3pt} \end{equation} \tag{ 1 }$$

where ξ, θ_c, and n are parameters that define the shape of the polynomial. Note that y = 0.5 is exactly mapped to the midplane. The merit of this θ-grid is that for n > 1 it ensures dense and nearly uniform cell elements close to the disk midplane. The azimuthal grid is evenly spaced and covers all 2π. Theoretically speaking, the aspect ratios of the cell shape should be as isotropic as possible, but the azimuthal cell widths can be a factor of a few longer than the other two, as the shearing tends to draw out features in this direction.

The grid resolution used in the present simulation is 400 × 160 × 540 in (r, θ, ϕ), with a computational domain covering [r_in, r_out] radially, [θ_c, π − θ_c] meridionally, and [0, 2π] azimuthally (see also in Table 1), where r_in = 0.8a, r_out = 16a, and θ_c = 0.2 is the width of the cutout around the polar axis. The other parameters used in Equation (1) are ξ = 0.9 and n = 9. Using this grid, we are able to resolve one disk scale height with 20 cells at the inner boundary of the computational domain. The resolution grows to ∼40 cells per scale height at radius ∼3a. Within two scale heights of the midplane, λ_MRI is resolved by more than six cells if the plasma β is not greater than ∼100. The grid also resolves λ_d with at least six cells in the r–ϕ plane for r ≲ 7a.

We choose outflow boundary conditions in the radial and meridional directions for the gas, which permit only flows going outward; any inward velocities are set to zero. As we cover all 2π of ϕ, the boundary conditions for ϕ-grid are simply periodic. For the magnetic field in the radial and meridional directions, we set the transverse components of the field to be zero in the ghost zones. The components normal to the boundaries are calculated by imposing the divergence-free constraint.

We begin with a prograde disk orbiting in the binary plane. Between r_min = 3a and r_max = 6a, its density is constant in the midplane. The initial disk is axisymmetric, and the polar angle density distribution is $\rho =\rho _0 \exp {[-(\theta -\pi /2)^2/(\sqrt{2}H/r)^2]}$, in which ρ₀ = 1 is the unit of disk density. This form provides initial hydrostatic balance vertically for a point-mass potential and zero radial pressure gradient along the midplane. For a first-order approximation, the difference between a binary potential and a point mass in the midplane can be described by the temporally and azimuthally averaged quadrupole moment of the binary potential. We therefore modify the angular frequency of the initial disk to account for the quadrupole contribution:

$$\begin{equation} \Omega (r)^2 \approx \Omega _{\rm K}^2\left[1+\frac{3}{4}\left(\frac{a}{r}\right)^2\frac{q}{(1+q)^2}\right], \end{equation} \tag{ 2 }$$

where q = 1 is the mass ratio of the binary. Here we replace R with r as r = Rsin θ ≈ R for regions near the midplane. We also add 1/(rΣ)dP/dr to the right-hand side of Equation (2) to compensate for the small radial gradient of the vertically integrated pressure P = Σc²_s.

The initial magnetic field is a single poloidal loop within the main body of the disk. It is subthermal with plasma β = 100 on average. The field is constructed from the vector potential $\mathbf {A} = (0,0,A_{\phi })$, and we define A_ϕ by $A_{\phi }= A_0\sqrt{\rho }\sin (2\pi r/kH)(r/r_{\rm min}-1)(1-r/r_{\rm max})-\sqrt{\rho _{\rm cut}}$ if A_ϕ > 0 and zero otherwise, where k = 2Ω_bina/c_s, ρ_cut = 10⁻³ρ₀, and A₀ is a constant determined by the constraint on the averaged β.

Three-dimensional MHD simulations usually produce a large amount of data, which pose challenges for data storage, transport, and post-simulation analysis. In order to facilitate the study of the spatial and temporal properties of the disk, we write out spatially averaged history data at a frequency of one dump per time unit and write out 3D snapshots every five units of time.

We use two different types of spatially averaged history data. The first is defined by either an integration or an average over radial shells. Shell averaging for variable X is defined by

$$\begin{equation} \langle X\rangle (r,t) \equiv \frac{\int X r^2\sin \theta d\theta d\phi }{\int r^2\sin \theta d\theta d\phi }, \end{equation} \tag{ 3 }$$

where ∫r²sin θdθdϕ ≡ A_x(r) is the shell surface area. With this definition, the density-weighted shell average is 〈X〉_ρ ≡ 〈ρX〉/〈ρ〉. For example, the net disk accretion rate is

$$\begin{equation} \dot{M}(r,t)\equiv \int \rho v_{\rm r} r^2 \sin \theta d\theta d\phi =A_{\rm x}\langle \rho v_r\rangle, \end{equation} \tag{ 4 }$$

and the average specific angular momentum l = 〈v_ϕr〉_ρ = 〈ρv_ϕr〉/〈ρ〉. The surface density is a vertically integrated and azimuthally averaged quantity:

$$\begin{equation} \Sigma (r,t) \equiv \frac{1}{2\pi }\int \rho r \sin \theta d\theta d\phi. \end{equation} \tag{ 5 }$$

The second type of average is 2D, either an azimuthal average of poloidal slices or a vertical average referred to the equatorial plane. We define the vertical average by

$$\begin{equation} \langle Y\rangle _z (r,\phi,t) \equiv \frac{\int Y r\sin \theta d\theta }{\int r\sin \theta d\theta } \end{equation} \tag{ 6 }$$

and the density-weighted vertical average by

$$\begin{equation} \langle Y\rangle _{z,\rho }(r,\phi,t)\equiv \frac{\int \rho Y r\sin \theta d\theta }{\int \rho r\sin \theta d\theta } = \frac{\langle \rho Y\rangle _z}{\langle \rho \rangle _z}. \end{equation} \tag{ 7 }$$

Similarly we have the azimuthal average as 〈Y〉_ϕ(r, θ, t) ≡ 1/2π∫Ydϕ.

We need to be very careful about the definition of the r–ϕ component of the internal stress in the present simulation. We follow Hawley & Krolik (2001) and define the stress as

$$\begin{equation} w_{r\phi } (r,\theta,\phi,t)= t_{r\phi } + r_{r\phi } = -\frac{B_rB_\phi }{4\pi } + \rho \delta v_r\delta v_{\phi }, \end{equation} \tag{ 8 }$$

where t_rϕ is the Maxwell stress and r_rϕ is the Reynolds stress. The perturbed velocities δv_r and δv_ϕ are calculated from

$$\begin{eqnarray} \delta v_r (r,\theta,\phi,t) &=& v_r - \langle v_r\rangle _{\rho } \nonumber \\ \delta v_{\phi }(r,\theta,\phi,t) &=& v_{\phi } - \langle v_{\phi }\rangle _{\rho }. \end{eqnarray} \tag{ 9 }$$

The vertically integrated and azimuthally averaged total stress can then be described by

$$\begin{eqnarray} W_{r\phi }&\equiv & \frac{1}{2\pi }\int w_{r\phi }r\sin \theta d\theta d\phi = L_z\langle w_{r\phi }\rangle = T_{r\phi } + R_{r\phi },\nonumber \\ T_{r\phi }& =& L_z\left\langle \frac{-B_rB_{\phi }}{4\pi }\right\rangle,\nonumber \\ R_{r\phi }& =& L_z\left[\langle \rho v_r v_{\phi }\rangle - \frac{\langle \rho v_r\rangle \langle \rho v_{\phi }\rangle }{\langle \rho \rangle }\right], \end{eqnarray} \tag{ 10 }$$

where T_rϕ and R_rϕ are the average Maxwell stress and Reynolds stress, respectively, and L_z ≡ A_x/(2πr) is the vertical integral length, which in our case equals the height of the computational domain at a given radius r. In a similar way, we can obtain the vertically averaged stress 〈w_rϕ〉_z and azimuthally averaged stress 〈w_rϕ〉_ϕ as well.

To study the mechanism for disk eccentricity growth, we also carried out a set of numerical experiments using 2D viscous hydrodynamic simulations. Their purpose was both to distinguish hydrodynamic from MHD effects and to test the effect of where the inner boundary is placed.

The hydrodynamic simulations in this paper used the α-disk prescription. We chose α = 0.1 throughout the disk because that is roughly the ratio of stress to pressure in the disk body of our MHD simulation (note that MM08 used α = 0.01). The hydrodynamic disks are evolved in polar (r, ϕ) coordinates in the inertial frame. Following the grid scheme of the MHD simulation, we set the radial grid to be logarithmically spaced, while the azimuthal grid is spaced uniformly.

Among these simulations, B2D.rin=0.8 serves as the control run. The parameters that describe the physical properties of the hydrodynamic disk and the binary are kept the same as in the MHD case. The resolution of the control run is 512 × 1024 for r × ϕ. Its grid covers the same physical extent in the radial and azimuthal directions as the MHD one. The initial surface density of the 2D disk is simply taken from a vertical integration of the initial condition of the MHD disk. Other hydrodynamic simulations are reruns of B2D.rin=0.8 at t = 500 with various locations of the inner boundary. The initial disks of the reruns are obtained by truncating the restart data of B2D.rin=0.8 to the desired radius while keeping Δr/r and Δϕ fixed. The properties of both the reruns and the control run can be found in Table 1.

After the first 100 time units, the binary torque starts to clear out a low-density gap between the inner boundary and r ∼ 2a. In the top left panel of Figure 1, we show the vertically integrated surface density of the circumbinary disk at t = 120 in the x–y plane. In that panel, we can clearly see that the disk is truncated at around twice the binary separation. We also find two streams emanating from the disk edge toward the binary components. From the other three panels in Figure 1, we find that as the simulation continues the gap persists, the surface density gradually increases outside the gap, and finally an incomplete ring of dense gas forms due to the combined effects of mass accretion and gravitational torque. The stream also persists, but only one arm at a time. We will study the properties and effects of the transient stream in Sections 3.4 and 4. At t = 250 the disk appears to be elliptical and slightly off-center, and at t = 350 and 450 the disk obviously becomes eccentric. We also observe a growing azimuthal density asymmetry. We will discuss the eccentricity growth and density asymmetry in Sections 3.3 and 4. Despite the slow growth of both disk eccentricity and asymmetry, the principal disk structures, the gap and stream, are sufficiently time steady to suggest that in a qualitative sense the simulation has achieved a statistically stationary state. An analysis based on angle-averaged quantities will further demonstrate that.

Zoom In Zoom Out Reset image size

We show the enclosed mass (integrated disk mass interior to a given radius) inside r = a, 1.5a, 2a, 3a, and 4a as a function of time in Figure 2. We find that the mass at r < 3a (i.e., inside the initial disk's inner edge) undergoes dramatic growth during the first ∼100 time units as the initial disk fills in regions with <r_min. The radial pressure gradient at the edge of our initial disk leads to an inflow during the first several orbits, and once toroidal field develops near the disk inner edge, the Maxwell stress quickly removes the angular momentum of the low-density flow and drives inflow. As this transient phase passes, the mass in the gap levels off and stays quasi-steady (falling very slowly) after t ∼ 200. We therefore define the quasi-steady stage as the period between t = 200 and the end of the simulation. On the other hand, the mass inside r = 3a and 4a rises by a factor of two after t = 100, indicating that this region does not achieve inflow equilibrium. The mass-interior profiles also possess high-frequency fluctuations (typical timescale ∼3 time units) due to the orbital modulation of mass contained in the streams. In addition, after ∼350 time units, there are slower oscillations (timescale ≃ 20–30 time units) of the enclosed mass caused by the eccentrically orbiting density lump.

Zoom In Zoom Out Reset image size

In Figure 3, we display the time-averaged surface density and accretion rate over two intervals: t = 250–350 (ΔT₁) and t = 350–450 (ΔT₂). We find that there is only a minor change in the averaged surface density and accretion rate of the disk region <2.5a over ΔT₁ and ΔT₂. If we define the unit of surface density Σ₀ ≡ ρ₀a, the peak density grows in time, but only slowly, rising from ≃ 0.6Σ₀ to ≃ 0.7Σ₀. Its position r_p also changes, but similarly slowly, gradually moving from ≃ 2.5a to ≃ 3a. The shape of the surface density profile likewise hardly changes in time. On the inside, the disk is truncated exponentially, 〈Σ(r < 2a)〉_t ∼ Σ₀exp (3.8r/a − 8.4) for both intervals, where 〈〉_t represents an average over time. Toward larger radii, the surface density is ∝r⁻², the predicted profile of a "decretion disk" (Pringle 1991). All that changes is the region where Σ(r)∝r⁻² extends farther outward at later times.

Zoom In Zoom Out Reset image size

Averaged over the quasi-steady period, the accretion rate at the inner boundary is $\dot{M}\simeq 0.018(GMa)^{1/2}\Sigma _0$. However, the disk evolves throughout the simulation at r > 3a. We see gas inflow rates as large as ∼0.06(GMa)^1/2Σ₀ outside the region ∼3a during both time intervals. As a result, the surface density continuously builds up around this radius.

Many properties of these disks can be expected to scale with the mass near the surface density peak. It is therefore convenient to define a "disk mass" M_d ≡ πr²_p〈Σ_p〉_t, where r_p = 3a and 〈Σ_p〉_t = 0.65Σ₀. In terms of this mass (and the natural unit of time for the simulation), the mean accretion rate through the inner boundary can be rewritten as

$$\begin{equation} \langle \dot{M} \rangle = 1.0 \times 10^{-3}M_{\rm d}\Omega _{\rm bin}. \end{equation} \tag{ 11 }$$

The surface density of such disks is very uncertain. A sense of scale, however, can be gleaned from translating Equation (11) into the luminosity that would be produced if the accretion were converted into radiation at customary black hole efficiency (10%); in Eddington units, it is L/L_E = 0.036M^−1/2₈a^1/2_0.1τ_p, where we have scaled to a total black hole mass of 10⁸ M_☉, a binary separation of 0.1 pc, and a disk column density whose Thomson optical depth is τ_p.

Unlike a steady accretion disk in a point-mass potential, where the internal stress simply transports angular momentum outward and thereby drives an inflow, the binary consistently interacts gravitationally with the circumbinary disk by torquing the surrounding gas. The angular momentum delivered through these torques is transported by two mechanisms: MHD stresses, mostly due to MRI-driven turbulence, and Reynolds stresses associated with coherent gas motions in the gap region. We have checked that the code's numerical shear viscosity is negligible compared to the transport mechanisms discussed here. In this section, we first investigate the mechanisms of angular momentum transport and the properties of the binary torque and then discuss the balance between the binary torques and stresses.

3.2.1. Maxwell and Reynolds Stresses

In the first panel of Figure 4, we show the vertically integrated and time- and azimuthally averaged stresses and their sum as a function of radius. Within r < 3a, the character of these stresses is very similar in our two averaging periods, ΔT₁ and ΔT₂. The Reynolds stress shows big fluctuations radially. Its first peak, at r ≃ 1.7a, is ∼0.35Σ₀c²_s. It dips at around 3a to nearly zero and then rises up again slowly. The ratio of the first peak to the second peak is ∼10. However, the Maxwell stress has a flatter radial profile. It is roughly ∼0.1Σ₀c²_s on average and varies only slowly over a range of at most a factor of two. It diminishes within the gap and decreases smoothly toward larger distance. Comparing the two, we find that the Reynolds stress exceeds the Maxwell stress by a factor ∼4 at r ∼ 1a–2a. On the other hand, beyond the gap region, the MHD stress always dominates the Reynolds stress by a factor ≳ 3. Within one binary separation, the Maxwell stress also exceeds the Reynolds stress, mainly because the Reynolds stress falls with the lower gas density near the boundary. The total internal stress therefore follows the Reynolds stress at r ≲ 2a, while it is close to the magnetic stress outside that region. Similar results have been reported in previous studies on gap formation in protoplanetary disks with an embedded planet using MHD simulations (Nelson & Papaloizou 2003; Winters et al. 2003). It is not a surprising outcome as the disk with a planet is just a special case of a circumbinary with an extreme mass ratio.

Zoom In Zoom Out Reset image size

In the language of the α_SS parameter of Shakura & Sunyaev (1973), we can also measure the stresses in units of the pressure. In the bottom panel of Figure 4, we present the time-averaged stresses in these units. The most significant distinction compared with the absolute stresses is that the stress ratios increase steeply inward due to the low surface density in the gap. The total stress ratio peaks at α_SS ≃ 11 at r = 1.1a. Its value in the disk body is almost two orders of magnitude smaller.

Our simulation shows that the circumbinary disk possesses highly time-dependent structures in the horizontal plane, making it important to look at the instantaneous spatial distribution of the stresses in addition to the time-averaged values. Vertically integrated snapshots of both the Reynolds stress (top left) and the Maxwell stress (top right) taken at t = 305 time units are shown in Figure 5. On top of each plot, we also draw the surface density with 15 contour lines, logarithmically spaced from the density floor value to the density maximum. Nearly all the Reynolds stress is confined within r < 3a. The most interesting point is that we find relatively large stress in a single gas stream extending from ∼1a to 2a (see the red stripe in the top left panel), where the gas gains angular momentum (δv_ϕ > 0) and is being pushed out by the binary. The gas being kicked thus goes on an eccentric orbit, creating a negative stress at r ∼ 2a–3a and −π/4 < ϕ < π/2 (δv_r > 0 and δv_ϕ < 0). In the same range of radii but a different azimuthal location (π/2 < ϕ < 5π/4), the gas falls back toward the binary, creating positive stress (both δv_r and δv_ϕ < 0). This bi-symmetric property provides a near cancellation so that only the stress within the stream contributes significantly to the radial profile of the Reynolds stress. On the other hand, the Maxwell stress is more evenly distributed in the disk body, although the typical magnitude is a factor of four less than the typical magnitude of the Reynolds stress in the inner disk. In the gap, the Maxwell stress is strongly positive in the streams, where the field lines are collimated by the gas flow, but negative elsewhere.

Zoom In Zoom Out Reset image size

To illustrate how the stresses depend on the disk height, we present azimuthally averaged snapshots at the same time in Figure 5. We find that the Reynolds stress is well concentrated in the midplane. Its strength decreases to about half the midplane value at only one scale height away from the midplane, and about one-tenth at two scale heights. However, the Maxwell stress is less sensitive to the altitude. It appears to be mainly confined within four scale heights but also shows positive contributions at higher altitude near the pole due to field buoyancy.

3.2.2. Binary Torque

Having discussed the mechanism for angular momentum transport, let us now consider the torque exerted on the disk by the binary. We plot the time-averaged torque density and the integrated torque in Figure 6, where the torque is calculated approximately using the surface density instead of the local density itself. The definition follows MM08:

$$\begin{equation} \frac{dT}{dr} \equiv -\int \Sigma (r,\phi)\frac{\partial \Phi }{\partial \phi } rd\phi \end{equation} \tag{ 12 }$$

is the local torque, and the integrated torque is

$$\begin{equation} T(r) \equiv \int ^{r}_{r_{\rm in}} \frac{dT}{dr^{\prime }}dr^{\prime }, \end{equation} \tag{ 13 }$$

where Φ = Φ(r, ϕ) is the binary potential in the midplane.

Zoom In Zoom Out Reset image size

The first panel of Figure 6 shows that the averages of the local torques for the two intervals ΔT₁ and ΔT₂ are very similar, which again indicates that the inner part of the disk reaches a quasi-steady state. The local torque density is positive at r ≃ 1a–2a and peaks at r ≃ 1.5a–1.6a with dT/dr ≃ 0.03–0.04GMΣ₀ as the peak value. It then goes negative at r ≃ 2a–2.4a, reaching its first negative maximum at r ≃ 2.0a–2.2a, but with a smaller magnitude ≃ − 0.02GMΣ₀. Similar to what was previously found in MM08 and Cuadra et al. (2009), we find that the torque density oscillates around zero but damps quickly toward larger distance. In addition, the torque density for r ≲ a is negative because the gas within that region is advanced in phase (greater angular frequency) with respect to the binary, an effect that would appear in purely hydrodynamic treatments as well as MHD.

The integrated binary torque is displayed in the top right panel of Figure 6. The total binary torque exerted on the disk is T(∞) ≃ 0.011–0.013GMaΣ₀. In order to make a comparison with the torque found in MM08, some normalization is needed. First, we normalize the torque to the surface density at the radius where the local binary torque peaks. We denote that radius as r_torq. In our MHD simulation, the first positive peak of the torque takes place at r_torq ≃ 1.6a while it is ∼1.7a in MM08. In terms of this normalization, T(∞)/Σ(r_torq) ∼ 0.07GMa in the hydrodynamic simulation of MM08. In our MHD simulation, we find that it is ≃ 0.12GMa (averaged in the quasi-steady period), a factor of 1.7 greater than in the hydrodynamic and α-viscosity treatment.

However, as the peak of the torque is located in the gap region, where the surface density rises rapidly with increasing radius, this means of normalization is unreliable. We therefore turn to a different method, normalizing the binary torque to the disk mass. With this method, the normalized torque is T(∞)/M_d ≃ 6.5 × 10⁻⁴GMa⁻¹ in our MHD simulation. However, it is only ∼4.5 × 10⁻⁵GMa⁻¹ in MM08, ∼ 14 times smaller. Just as for the mass accretion rate, there is a natural set of units for the total torque: it is convenient to describe it in terms of how rapidly the torque removes the binary's orbital angular momentum. In that language,

$$\begin{equation} T({\infty }) = 2.6 \times 10^{-3} j_{\rm bin}M_{\rm d}\Omega _{\rm bin}, \end{equation} \tag{ 14 }$$

where j_bin = (GMa)^1/2/4 is the specific angular momentum of the binary.

Alternatively, as r_p is roughly around ≲ 3a for both MHD and hydrodynamic simulations, we can write the binary torque as some factor times GMaΣ_p, and then we have T(∞) ≃ 1.9 × 10⁻²GMaΣ_p in our present simulation, while MM08 gives ∼1.4 × 10⁻³GMaΣ_p, ∼14 times smaller. The largest part of this contrast can be explained by the fact that in the MHD treatment Σ(r_torq)/Σ_p ≃ 0.15, a factor of ∼8 greater than that of MM08.

This contrast can also be understood in terms of the angular momentum budget. The gravitational torque per unit mass approximately balances the torque per unit mass due to the internal stress near the disk inner edge. Therefore, the binary torque should increase linearly with the effective α. In our case, α ∼ α_SS ≳ 0.2 in the disk body at r ≃ 2a–5a, a factor of ≳ 20 greater than the constant α = 0.01 assumed by MM08.

In order to show how the local binary torque evolves with time, we also plot in Figure 6 the specific torque density (bottom left), i.e., the absolute value of the ratio of torque density to surface density. On dividing by the surface density, we remove the effects due to redistribution of the surface density. We find that the specific torque slightly shifts to larger radius and diminishes its amplitude from ΔT₁ to ΔT₂. In the bottom right panel, we plot the history of the total torque smoothed over a short time span ∼2T_bin, where T_bin = 2π/Ω_bin is the binary period. Before smoothing, the total torque fluctuates very rapidly (∼1.5Ω_bin) between ∼ − 0.05GMaΣ₀ and 0.06GMaΣ₀ in the quasi-steady period. The power spectrum of this fluctuation is dominated by a peak at ∼1.5Ω_bin, approximately the beat frequency between the orbital frequency of the matter being torqued most strongly (at the peak of dT/dr, i.e., r ∼ 1.6a–1.7a) and 2Ω_bin. After smoothing, we find that the binary torque has a sudden rise around t = 60, which is due to the initial disk filling the small radius region. The torque then levels out between t = 100 and 200. Once in the quasi-steady stage, the torque shows a slowly decreasing trend, ≲ 30% in fractional terms over the final ∼300 time units. It also appears to oscillate at the disk orbital frequency (r ∼ 2a–3a) during the late time, a fluctuation caused by the eccentric movement of the density lump. As both the disk eccentricity and the density of the lump grow with time, the torque normalized by the time-dependent Σ_pr²_p simply reflects the increasing Σ_p and r_p.

3.2.3. Torque Balance

For a steady circumbinary disk, the angular momentum at a given radius should not fluctuate too much over time. Therefore, the angular momentum deposition from the binary must be either transported outward by Maxwell and Reynolds stresses or advected toward the hole. In order to test this idea, we perform an estimate of the differential angular momentum flux averaged over a short time interval between t = 300 and 320 using the high time-resolution rerun. Based on the conservation of angular momentum, we write the angular momentum budget for a ring of disk between (r − Δr, r + Δr) as

$$\begin{eqnarray} \frac{\partial J}{\partial t}&+& \frac{\pi }{\Delta r}[(r+\Delta r)^2F(r+\Delta r)-(r-\Delta r)^2F(r-\Delta r)] \nonumber \\ &=& \frac{\pi }{\Delta r}[(r-\Delta r)^2W_{r\phi }(r-\Delta r)\nonumber\\ &&-(r+\Delta r)^2 W_{r\phi }(r+\Delta r)] + \frac{dT}{dr}, \end{eqnarray} \tag{ 15 }$$

where J(r, t) = 2πr〈ρv_ϕr〉L_z is the shell-integrated angular momentum, F(r, t) = 〈ρv_r〉〈ρv_ϕ〉L_z/〈ρ〉 is the mean angular momentum flux due to mass motion, and W_rϕ(r, t) = T_rϕ + R_rϕ is the total internal stress. We write the conservation law this way so that a negative differential flux due to advection (the second term on the left-hand side of Equation (15)) indicates a net angular momentum inflow, while a negative flux due to internal stresses (the first term on the right-hand side) means angular momentum is removed and transported to larger radii. We calculate the time averages of all these terms and plot them in Figure 7 as functions of radius. In its y-legend, we use dF_J/dr to represent all the terms in the equation, which are all radial derivatives of one or another variety of shell-integrated angular momentum flux (F_J).

Zoom In Zoom Out Reset image size

During this time interval, the local angular momentum stays steady, with only very small temporal variations from the averaged ∂J/∂t (solid black curve), as expected for a quasi-steady accretion disk. The binary torque (blue dashed curve) at r ∼ 1a–3a is almost balanced by the differential flux due to the sum of the Reynolds stress and the Maxwell stress (cyan dashed curve). In other words, the circumbinary disk within that region reaches a balance between the binary torque and internal stresses. As a result, the angular momentum flux carried by accretion flows oscillates about zero (magenta curves within 1a–3a). Clearly, the internal torque due to the Reynolds stress (green dash-dotted curve) dominates the Maxwell stress (red dash-dotted curve) at r ∼ 1a–3a, which suggests that most of the angular momentum dumped by the binary at r ∼ 1a–1.8a is transported outward via the Reynolds stress, while the angular momentum drawn from the binary at r ∼ 1.8a–2.4a is compensated by the torque of the Reynolds stress. Finally, we comment that torque balance is not reached at r ≲ a (where the negative binary torque dominates) or r ≳ 3a (where the internal stress dominates). By summing all the other torques, we find the differential angular momentum flux that must be advected inward in those two regions (magenta curve). Directly computing the second term on the left-hand side of Equation (15) yields the (very similar) black dashed curve in Figure 7.

To sum up, we have found three principal conclusions about the angular momentum budget. First, the Reynolds stress dominates the Maxwell stress in the gap region, while the latter exceeds the former outside that. Second, the distribution of both stresses follows the gas streams and is vertically confined: within two scale heights for the Reynolds stress and four scale heights for the Maxwell stress. Third, we find that the radial profile of the binary torque is similar to that of previous hydrodynamic results, but the amplitude is different: normalized by the surface density at the location where local torque peaks, the binary torque in the present MHD simulation is twice as great as in previous hydrodynamic calculations; however, when normalized by the disk mass, our torque is an order of magnitude greater than that found in MM08. Lastly, we find that in the quasi-steady state, the binary torque is roughly balanced by the torque due to the Reynolds and Maxwell stresses in the inner disk (r < 3a). Most of the angular momentum dumped by the binary in the gap region is transported outward by the torque of the Reynolds stress associated with the gas streams.

3.3.1. Evolution of the Disk Eccentricity

Although we begin our simulation with an axisymmetric disk configuration, which means there is zero eccentricity at the start, the disk eccentricity grows throughout our simulation (e.g., compare the disk surface density snapshots in Figure 1).

To quantify the eccentricity, we define the local eccentricity as

$$\begin{equation} e(r,t) = \frac{1}{T_{\rm bin}}\int ^{t+T_{\rm bin}}_{t}\frac{\big|\big\langle \rho v_r e^{i\phi }\big\rangle \big|}{\langle \rho v_{\phi }\rangle } dt^{\prime }, \end{equation} \tag{ 16 }$$

where 〈〉 is the shell average defined as in Equation (3). We show this local eccentricity in the spacetime diagram in Figure 8 (top panel). The local eccentricity of the gas in the gap (r ≲ 2a) rises significantly during the period t ∼ 100–200. The increasing e(r) in the gap simply reflects the fact that the dynamics there are strongly non-Keplerian. The eccentricity in the main body of the disk (r > 2a) grows rapidly over the period t = 100–300. The typical value increases from much less than 0.01 throughout the disk body at t = 100 to ∼0.08 between r = 2a and r = 3a at t = 300. It keeps growing slowly after 300 time units, and at the end of the simulation a ring of eccentric disk forms at 2a < r < 3.3a with typical e(r)  ∼  0.1. A more detailed picture of time variation in the eccentricity is shown in the bottom panel of Figure 8, where we show its evolution from t = 300 to 320. For clarity, we use the data from the high time-resolution simulation and do not take the temporal smoothing as in Equation (16). We find strong radial extending structures that connect the gap with the disk body. These features bend at ∼1.5a pointing to later times at both ends, which suggests the propagation of eccentricity from the gap to both the disk and the binary. The pattern repeats once every half binary period, indicating a stream-related origin. We will discuss the stream effects on the disk eccentricity in Section 4.3.3.

Zoom In Zoom Out Reset image size

In Figure 9, we present the radial distribution of the disk eccentricity averaged over three time intervals: t = 150–250 (solid black), ΔT₁ (red dotted), and ΔT₂ (blue dashed). The distribution curves in the disk body (r > 2a) shift toward larger radius and greater eccentricity as the disk evolves in time. Outside the peaks around 2a, the slopes of the distribution curves are constant. The radial dependence of e(r) in the disk body can be roughly described as ∝exp [ − 1.3(r/a)].

Zoom In Zoom Out Reset image size

Interestingly, the eccentricity distribution found by MM08 was quite similar to ours in both amplitude and shape, despite that calculation's very different internal stresses, initial surface density profile, and duration. In our own 2D hydrodynamic simulations, we have likewise found qualitatively similar eccentricity growth. However, the eccentricity profile in those simulations matches the MHD simulation (and MM08) only at a specific time (t ∼ 800Ω⁻¹_bin). From these comparisons, it is clear that the primary mechanism for driving eccentricity must be hydrodynamic response to gravitational forcing by the binary, but its quantitative development can be affected by both disk physics (e.g., MHD stresses) and initial conditions (e.g., the initial surface density profile). We will discuss specific mechanisms at greater length in Section 4.3.1.

The spacetime diagram suggests we may take the radial average of the local eccentricity for the disk body between r = 2a and r = 4a as a measure of the disk eccentricity. In Figure 10, we show the averaged disk eccentricity e_disk constructed by taking volume averages from r = 2a to r = 4a and time averaging over one binary period:

$$\begin{equation} e_{\rm disk} =\frac{1}{T_{\rm bin}}\int ^{t+T_{\rm bin}}_{t} dt^{\prime } \frac{\big| \int ^{4a}_{2a}dr \int \rho v_{r} e^{i\phi } r^2\sin \theta d\theta d\phi \big|}{\int ^{4a}_{2a} dr\int \rho v_{\phi } r^2\sin \theta d\theta d\phi }. \end{equation} \tag{ 17 }$$

At early times (t ≲ 100), Figure 10 shows the transient phase of a non-equilibrium disk in a non-Keplerian potential. During this phase, the noise develops significantly, especially in the low-density region of the initial configuration, and all modes undergo dramatic growth. After the transient, the eccentricity of the disk undergoes exponential growth from t = 100 to t = 250. By fitting the eccentricity curve during that interval, we find that the growth rate is γ_e ≃ 0.018Ω_bin, where the subscript e denotes eccentricity. The growth rate of disk eccentricity slows by a factor of four after ∼250. This suggests that some nonlinear mechanisms come into play to limit and perhaps eventually saturate eccentricity growth.

Zoom In Zoom Out Reset image size

3.3.2. Precession of the Eccentric Disk

To follow changes in the orientation of the eccentricity, we define the longitude of the apoapse ϖ in terms of the complex phase angle of the disk eccentricity for 2a ⩽ r ⩽ 4a. Its time evolution is shown in Figure 11. Here the angle ϖ measures the angular location of the apoapse with respect to the x-direction in the inertial frame. The disk eccentricity is vigorously perturbed by the motion of the gas in the gap during the first one hundred units, and the arbitrary orientation in that period is partly because the eccentricity during this time interval is very small. The angle ϖ then gradually rotates in a prograde manner. By fitting the curve between t = 200 and 480, we find $\dot{\varpi } \simeq 3.2\times 10^{-3}\Omega _{\rm bin}$. Linear perturbation theory predicts that the precession rate for particles around the binary potential is $\dot{\varpi }\simeq \Omega - \kappa$ for a cold disk, where Ω is the angular frequency and κ is the epicyclic frequency of the disk. To a first-order approximation, Ω(r) ≈ Ω_K[1 + (3/8)(a/r)²q/(1 + q)²] (also see Equation (2)), and κ(r) ≈ Ω_K[1 − (3/8)(a/r)²q/(1 + q)²]. Thus, we find that the prediction is $\dot{\varpi }\simeq 0.19 (a/r)^{7/2}\Omega _{\rm bin}$ for q = 1, or ≃ 4.1 × 10⁻³Ω_bin evaluated at r ≃ 3a. The precession rates measured in both our MHD simulation and MM08 are close to this prediction. We therefore suggest that the precession is mainly due to the quadrupolar potential.

Zoom In Zoom Out Reset image size

3.3.3. Late-time Lump

There is another disk characteristic closely connected to the eccentricity: the density concentration that appears near the inner disk edge at late time. We call it the "lump" (see Figure 1). The lump orbits around the binary at nearly the same orbital speed as the eccentric disk. To quantitatively study the asymmetric feature, we introduce the m = 1 mode strength component

$$\begin{equation} S_{{m=1}}(r,t) \equiv \left|\int \rho e^{im\phi } rd\theta d\phi \right| (m=1), \end{equation} \tag{ 18 }$$

which is basically the m = 1 component of the Fourier transform of the surface density. We plot the spacetime diagram of S_{m = 1} in Figure 12. The orbiting lump appears as a zigzag pattern whose peak moves in and out between r ∼ 2a and 4a after about 300 time units. The pattern period, which is also the orbital period of the eccentric disk, is ∼30 time units. We will discuss the lump in Section 4.4 and explain it as due to a combination of disk eccentricity and the action of periodic streams impacting upon the edge of the disk.

Zoom In Zoom Out Reset image size

As shown in Figure 1, gas is not entirely absent from the gap. Consistent gas flows launch from the inner edge of the disk and stream toward the binary. The streams affect not only the accretion rate of the circumbinary disk but also the structure of the inner part of the disk.

In a frame that corotates with the binary, we find that the streams follow roughly the static bi-symmetric potential of the binary and pass through the inner boundary about 0.3–0.4 rad beyond the Lagrange L2 and L3 points (see Figure 13). This phase offset is consistent with the negative average torque at r ≲ a in Section 3.2.2. Early in the simulation, the two streams have roughly equal strength, i.e., similar density contrast between the stream and the nearby gap region. We show the surface density near the inner edge at t = 122 and 125 in the first two panels of Figure 13 and draw the velocity field vectors measured in the comoving frame on top of it. The snapshots at t = 122 show that the inflowing gas is regulated by the binary torque as it leaves the inner disk edge and streams toward the saddle points. The sign of the angular speed near the inner boundary alternates from one quadrant to the next, a direct signal that the binary torque changes its sign from one quadrant to the other. In the corotating frame of the binary, the gas in the second (near ϕ = π) and fourth (near ϕ = 0) quadrants is pulled toward the closer component of the binary (see t = 122 snapshots). The tangential velocity results in a Coriolis force, which tends to move matter away from the center. It takes about half the binary period for the gas streams in those two quadrants to be kicked out and strongly interact with the disk edge. As shown in the later snapshots at t = 125, the velocity is mostly outward at r ∼ 1.5a in those two quadrants, and there are strong density enhancements near the regions where the velocity gradient is large. There is also an equivalent way to think about the transition between t = 122 and t = 125. In the inertial frame, the binary almost always rotates faster than the disk; therefore, the second and fourth quadrants have a positive binary torque. As a result, gas in these regions gains angular momentum from the binary and then drifts away from the center.

Zoom In Zoom Out Reset image size

However, the disk inner edge changes its morphology once the eccentricity of the disk becomes significant. Instead of a pair of streams with comparable size and strength, we find that, for instance at t = 301, the stream on the right dominates the left. For an eccentric disk, the apocenter side of the disk is far from both members of the binary. Thus, the gas stream is strongly torqued, and therefore carries more mass, only when one of the members gets close to the pericenter of the disk. After half a binary period, the binary members switch their phases (in the inertial frame) and a strong stream forms around the other member. We illustrate this process with a time sequence of snapshots at t = 301, 302, 303, and 304.1 in the corotating frame. The stream on the right at t = 301 is gradually pushed away at t = 302 and 303 and joins the apocenter of the eccentric disk. Meanwhile, the stream associated with the other member of the binary strengthens as that mass approaches the pericenter at t = 304.1.

We also measure the fluid and magnetic effects on the dynamics of the streams. In the inner disk and the gap, the force ratio $|\mathbf {f_{g}}+\mathbf {f_{L}}|/|\nabla \Phi |$ is always ≲ 1/3, where $\mathbf {f_{g}}=-\nabla P/\rho$ is the force density of the gas pressure gradient and $\mathbf {f_{L}}=(\nabla \times \mathbf {B})\times \nabla \mathbf {B}/4\pi \rho$ is the Lorentz force density. Because gas pressure and magnetic forces are smaller than the gravitational force, the trajectory of the gas is essentially ballistic. Consequently, we neglect the fluid and magnetic effects when calculating the interaction between the streams and the disk inner edge in Section 4.3.3.

We found in the previous subsection that the Maxwell stress closely follows the gas streams in the x–y plane and is well confined near to the midplane. Now we show that the field structure in the gap region, especially in the stream, has similar features. In Figure 14, we draw the vertical averaged magnetic field as vectors on top of the logarithmic scaled surface density at t = 300.5. The location and length of the arrows show the position and relative strength of the field. Clearly, the field is well collimated by the motion of the gas streams, with stronger inward field lines in the lower stream (bottom-half plane) and weaker outward field in the upper stream (top-half plane). The ordered field in the gap therefore produces greater Maxwell stress only within the two streams.

Zoom In Zoom Out Reset image size

We present the accretion rate through the inner boundary as a function of time in the top panel of Figure 15. The sampling rate of the time sequence is determined by the output rate into the history files, which was once per Ω⁻¹_bin. There are high-frequency outbursts (the narrow spikes in the time sequence) and also a low-frequency modulation (the dashed curve, which is derived by smoothing the time sequence using a boxcar average with width 15 time units). These two time variations are the most significant modes in the temporal profile of $\dot{M}$. By constructing the power spectral density of the accretion rate using Fourier transforms (middle panel, Figure 15), we are able to identify them: the higher frequency mode is at about twice the binary rotation rate, and it is the same rate at which the streams are pulled inward by the binary potential; the lower frequency mode is ∼0.2Ω_bin, which is the same as the orbital frequency of the lump during the later stages of the simulation (as shown in Section 3.3). Our result in Figure 15 is largely consistent with previous findings: a dominant frequency associated with the binary orbital frequency (2Ω_bin if q = 1; MM08; Hayasaki et al. 2008; Cuadra et al. 2009; Rödig et al. 2011; Sesana et al. 2012), a low-frequency component due to the motion of the dense part of the disk (MM08; Rödig et al. 2011), and another component (in our case at 1.8Ω_bin) created by a beat between the binary orbital frequency and the disk orbital frequency (Rödig et al. 2011). The principal contrast between our results and this earlier work is that the dominant frequency is 2Ω_bin rather than Ω_bin because the members of the binary have exactly equal mass.

Zoom In Zoom Out Reset image size

To further verify the cause of these two modes, we pick out two times (t = 400 and 441 as marked in the top panel by black arrows) and plot their 2D local accretion rates in the inertial frame in the bottom panel of the same graph. There were accretion outbursts at both times, but the former is in the valley of a low-frequency oscillation, while the latter is at the crest. In the snapshots, the color contours represent the vertically integrated accretion rate. The black contour lines show the surface density in a logarithmic scale from 10⁻⁴ to 10^−0.5. These lines pick out the location of the streams, while the white contour lines show the surface density between 1.4 and 3.0 in a linear scale; these contours delineate the lump region. We find that the outbursts of $\dot{M}$ are indeed coming from the infalling streams located in the fourth quadrants of both snapshots. However, the local accretion rate in the stream from the left panel is much weaker than that from the right. In the right plot, the lump is passing the pericenter region. The binary then pulls a gas stream with relatively higher density than other times. Thus, we conclude that the time-dependent accretion shows both high-frequency outbursts corresponding to the periodic inflow from the streams and low-frequency modulation due to the orbiting lump. The low-frequency mode becomes important only when the lump forms and maintains itself after about t = 350.

The time-averaged accretion rate at the inner edge, normalized with the peak density, is ${\sim} 0.028\sqrt{GMa}\Sigma _{\rm p}$, a factor of ∼40 greater than that found in MM08. Interestingly, during the quasi-steady period, the accretion rate in the gap is nearly twice the rate measured at r_p, where the surface density peaks. This contrast explains the slow mass depletion inside the gap shown in Figure 2. One of the most important questions about accretion onto a binary is the ratio of the accretion rate accepted by the binary to the rate at which mass is fed into the disk from the outside. If the latter rate can be estimated by the maximum accretion rate in the disk, we find that accretion through the inner boundary of the simulation, and therefore presumably onto the binary, is about 1/3 of the accretion rate supplied. This ratio is a factor of two greater than previously found by MM08. However, in evaluating these numbers one should be aware that the disk as a whole never reaches equilibrium over the course of the simulation. That the time-averaged accretion rate as a function of radius is far from constant is a symptom of this fact. Thus, our estimate (and for the same reasons, that of MM08 as well) should be taken as no better than an order-of-magnitude estimate of this fraction. We plan a more rigorous approach to this problem in future work. Lastly, we note that the effects of the excision on the accretion rate are not explored in this MHD simulation. However, our hydrodynamic simulations with different cutouts suggest that the rate hardly changes as long as the cutout size is ≲ a.

Using the results we have on the disk stresses and the binary torque, we want to further elucidate the angular momentum transport mechanisms within an MHD turbulent circumbinary disk. The binary strongly torques the surrounding gas in the gap as the gas forms streams; the outward-moving portions of the streams deliver that angular momentum to the inner part of the circumbinary disk. By contrast, the torque directly exerted on the disk body is much weaker and diminishes rapidly toward larger distance. In that sense, we can say that the binary mainly torques the gas in the low-density gap; moreover, the magnitude of that torque is proportional to the mass of gas in the gap. This fact suggests that the binary torque is very sensitive to the basic physics coupling the streams to the disk such as the stresses and the equation of state. For instance, scaled with the peak surface density, our MHD turbulent disk provides a factor of 14 times stronger binary torque than that of MM08.

This fact also leads to a strong contrast with the commonly used linear resonance theory of interaction between disks and planets or other orbiting satellites (Goldreich & Tremaine 1982). To illustrate that contrast, we calculate both the torque density and total torque as functions of radius, assuming that most of the torque comes from the Ω_bin: Ω = 3: 2 outer Lindblad resonance at r₂/a = (3/2)^2/3 ≃ 1.3 (e.g., Meyer-Vernet & Sicardy 1987). The estimate is shown in the top panel of Figure 6, with the red lines representing 1/4 of the predicted values. The linear theory prediction must be reduced by roughly that amount to match the measured total torque. In addition to having a larger amplitude than the simulation found, the linear theory also predicts a much shorter wavelength for its radial oscillations.

Three distinctions between the circumstances of the MHD circumbinary disk we simulate and the assumptions made in linear theory may explain these contrasts. The first is the surface density profile. Standard linear theory assumes that any surface density gradients have length scales much longer than the wavelength of the density wave excited by the resonance. However, the surface density profile in the gap is steep enough that (dln Σ/dr)⁻¹ ≲ 0.3a, which is comparable to the first wavelength 5[(m + 1)/(3m²)]^1/3(H/r)^2/3a ∼ 0.5a, where H/r is evaluated at the resonance radius, and m = 2 is used here (see Figure 1 and the pressure-dominated case in Table 1 of Meyer-Vernet & Sicardy 1987). Because the surface density at r₂ is so much smaller than at r ∼ 2a–3a, one might expect the maximum torque to move away from the resonance point toward the inner disk. Second, linear theory assumes nearly circular orbits because it also assumes q ≪ 1. However, in our simulation q = 1, and, largely for this reason, the orbits of gas within the inner edge of the disk body are significantly non-Keplerian: they form streams. Because the torque depends so strongly on whatever gas mass is in the gap, this distinction is qualitatively important. Lastly, as the mass ratio of the binary of our simulation is unity, the forcing term, which is roughly controlled by the ratio of the perturbing gravitational potential to the sound speed, becomes large enough to drive nonlinear density disturbances. The longer wavelength oscillation found in our torque density might be partly explained as a result of these nonlinear effects (Yuan & Cassen 1994).

The time-averaged rate at which the binary loses angular momentum to the disk is T(∞) ≃ 0.012GMaΣ₀. Meanwhile, the time-averaged angular momentum flux across the inner boundary due to accretion is $\dot{M} j_{\rm in} \simeq 0.017 GMa\Sigma _0$, where $j_{\rm in} \simeq 0.94 \sqrt{GMa}$ is the averaged specific angular momentum at r_in. On average, therefore, the binary actually gains angular momentum in net.

However, even though the binary gains angular momentum, it does not automatically follow that the binary separation grows. The evolution of the binary also depends on its growth in mass. If the orbit remains circular, the shrinkage rate for an equal-mass binary is

$$\begin{equation} \frac{\dot{a}}{a} = 2\frac{\dot{J}}{J} - 3 \frac{\dot{M}}{M}, \end{equation} \tag{ 19 }$$

where $J\,{=}\,M j_{\rm bin} \,{=}\, 0.25 M\sqrt{GMa}$ denotes the angular momentum of the binary and its rate of change is $\dot{J}=\dot{M} j_{\rm in} - T(\infty)$.

After separating the accretion term of $\dot{J}$ and combining it with the second term in Equation (19), we have

$$\begin{equation} \frac{\dot{a}}{a} = -2\frac{T(\infty)}{J} + 2 \frac{\dot{M}j_{\rm in}}{J}\left(1-\frac{3}{2}\frac{j_{\rm bin}}{j_{\rm in}}\right). \end{equation} \tag{ 20 }$$

If j_in/j_bin > 3/2 (here it is ≃ 3.75), the shrinkage rate is controlled by the competition between the first and second terms. Plugging in our numbers, we find that the two terms are surprisingly comparable in magnitude: $-0.096\sqrt{Ga/M}\Sigma _0$ and $0.082\sqrt{Ga/M}\Sigma _0$, respectively. The net result is orbital compression, but at a rate much smaller than either of the two terms, $\dot{a}/a\simeq -0.014 \sqrt{Ga/M}\Sigma _0$.

However, it is also possible that interaction with the disk may induce eccentricity. If so, an extra term on the right-hand side of Equation (19) should be included, making it

$$\begin{equation} \frac{\dot{a}}{a} = 2\frac{\dot{J}}{J} - 3 \frac{\dot{M}}{M} + \left(\frac{2e_{\rm bin}^2}{1-e_{\rm bin}^2}\right)\frac{\dot{e}_{\rm bin}}{e_{\rm bin}}, \end{equation} \tag{ 21 }$$

where e_bin denotes the binary's orbital eccentricity. In addition, the angular momentum of the binary should also be adjusted to account for the eccentricity: $J = M j_{\rm bin} \ssty\sqrt{1-e_{\rm bin}^2}$. In Appendix A, we show that the rate of eccentricity growth depends on how the accretion flow interacts with the interior disk of each individual binary member, an interaction our simulation did not treat. However, we also show that the best estimate we can make from the data we do have is that there is little evidence for any significant ${\dot{e}}_{\rm bin}/e_{\rm bin}$. Appendix A also briefly discusses the possibility of unstable eccentricity growth raised by linear theory; unfortunately, we can say little about it because our simulation data do not bear on it, and linear theory may not apply in this regime. If, as we did in Section 3.2.2, we define a characteristic disk mass M_d = πΣ_p(3a)², and further assume that the binary orbit remains circular, then $\dot{a}/a\simeq - 8.0 \times 10^{-4} ({M_{\rm d}}/{M})\Omega _{\rm bin}$. For comparison, MM08 found a rate ${\simeq} {-} 0.003 \sqrt{Ga/M}{\Sigma _0}$, which can be scaled using the definition of M_d to ≃ − 3 × 10⁻⁴(M_d/M)Ω_bin, somewhat slower than ours. In other words, although we find a torque an order of magnitude larger than that of an α-viscosity hydrodynamic model, we also find an accretion rate so much larger that the net effect on the binary orbit is reduced to an increase by only a factor of 2.7.

Another way to understand this result is to note that our $T(\infty)/\dot{M} \simeq 0.67 (GMa)^{1/2}$, a factor of three smaller than in MM08. That is, the greater internal stresses produced by MHD effects, especially in the gap region, both introduce more mass to that region (thereby magnifying the torque) and lead to even greater accretion, which returns angular momentum to the binary. Unfortunately, because the cancellation is so close, our calculation cannot be definitive in regard to the shrinkage rate of comparable-mass binaries in nature. Further uncertainty comes from the fact that we cannot track the gas flowing into the cutoff region with the present model, so we do not know exactly how much energy and angular momentum accretion might bring to the binary, yet our sizable accretion rate makes these contributions significant. The net outcome in any particular case may therefore depend on the specific details (gas equation of state, binary mass ratio, etc.), as all of these can influence j_in, $\dot{M}$, T(∞), as well as the detailed mechanics of how the stream joins the interior disks associated with the members of the binary.

4.3.1. Eccentricity Distribution

The distribution of eccentricity in radius plays an important role in the evolution of the disk. For the eccentricity to grow substantially over time in a large circumbinary disk, it needs to be confined in radius. Otherwise, the energy and angular momentum (more precisely, angular momentum deficit) associated with an eccentric disturbance are spread over a large radius, resulting in a small eccentricity in any one location. We consider here a linear model for the eccentricity distribution and compare it with the results of the simulations.

The linear equation governing the eccentricity distribution and evolution is given by Goodchild & Ogilvie (2006, hereafter GO06). This also allows generation and damping of the disk eccentricity. The basic equation and boundary conditions are given in Appendix B. For our circumbinary disk, we adopt input parameters taken from the simulations while the eccentricity was small and growing exponentially, i.e., the linear regime. We define the late stage of the exponential growth phase as the period t = 200–250 time units (see Figure 10). Surface density Σ(r) is taken directly from the simulation by averaging over this time interval. The disk inner radius is taken to be r_i = 2a and the outer radius r_o = 10a. We assume the same equation of state as in the simulation: isothermal, with sound speed 0.05Ω_bina. The adiabatic index is taken to be γ = 1 − 0.1i. The small imaginary part (which introduces slow damping) assures mathematical consistency. We model the forces driving the initiation of eccentricity as a Gaussian in radius (see Equation (B11)), centered on r_c = 2.1a and having a width w = 0.05a. These parameters are not well constrained by the simulation. It follows that the results are also insensitive to their values within certain limits. Parameter s_c, the peak value of the eccentricity injection distribution in Equation (B11), is determined by the condition that the growth rate equal that in the simulation, γ_e = −Im(ω) ≃ 0.018Ω_bin, where ω is the eigenvalue defined in Equation (B14). The value is found to be s_c = 0.0145aΩ_bin.

The physically relevant solution on long timescales corresponds to the fastest growing mode. This solution was determined by the methods described in Lubow (2010, hereafter L10). Figure 16 plots the eccentricity distribution for this mode, along with the eccentricity distribution obtained from the simulation. The middle dashed line corresponds to the eccentricity distribution based on linear theory with the set of parameters described above. The distribution matches well the one obtained from the simulation (the solid line). This distribution closely resembles the distribution for the fundamental free precession mode in the disk (the uppermost dashed line). That mode has an eccentricity injection of zero, s_c = 0. The confinement of eccentricity is due to the effects of differential precession of the binary $\partial _r \dot{\varpi }_{g}$ (see Equation (B8) in Appendix B), while pressure effects act to spread the mode. Differential precession effectively creates a potential well that leads to trapping of the fundamental mode. Further localization is due to the injection of eccentricity, which is concentrated near the disk inner edge (solid curve on the lower left). The faster eccentricity is injected into a small region, the more the eccentricity distribution narrows. For somewhat higher eccentricity injection rates, the distribution is further confined (the lowest dashed curve). Therefore, we see that the eccentricity distribution is a consequence of the excitation of the fundamental free eccentric mode in the circumbinary disk and that the eccentricity is confined within a radial width ∼a.

Zoom In Zoom Out Reset image size

The excitation of the fundamental mode occurs because that mode overlaps well with the eccentricity injection distribution. As discussed above, the location and width of the eccentricity injection are not well determined by the simulations. On the other hand, the linear equation has solutions that match the simulation results only if r_c ≲ 2.6a, which indicates that the source of the eccentricity should be located near the inner edge of our MHD disk. For larger values of r_c, there is no value of s_c that could provide the growth rate found in the simulation for an eccentricity distribution that resembles the simulated one. This radius limit has a weak dependence on the width of the injected eccentricity w. If the width is doubled to 0.1a, then the limiting radius increases by about 1%.

The outer regions of the simulation exhibit an exponential falloff of eccentricity with an e-folding length of about 0.8a (see Figure 9 in Section 3.3.1). In Appendix C, we derive an asymptotic form for the eccentricity distribution (valid for r ≫ a) using linear theory. It is dominated by an exponential falloff in r with an e-folding length also of about 0.8a.

To sum up the analysis in this subsection, we find that the linear model can explain the eccentricity distribution of our simulation. It shows that the confinement of the disk eccentricity is mainly due to the disk precession and eccentricity growth. The linear model also puts constraints on the location of the source for eccentricity growth. Finally, it explains the exponential falloff of the eccentricity profile.

4.3.2. Eccentricity Growth through Tidal Instability

One possible mechanism for eccentricity growth is a tidal instability. This instability is believed to explain the superhump phenomenon in cataclysmic binaries (see Osaki 1996). The instability involves the action of the 3:1 eccentric Lindblad resonance in a disk that surrounds a star (Lubow 1991, hereafter L91). As pointed out in L91, this instability can also be at work in circumbinary disks at radial locations that satisfy the condition

$$\begin{equation} \Omega (r) = \frac{m \Omega _{\rm bin}}{m+2}, \end{equation} \tag{ 22 }$$

where m is an integer that corresponds to the azimuthal wavenumber of the tidal component of the potential involved in the instability. In the case analyzed here of an equal-mass binary, m = 1 is excluded because that component of the tidal potential is absent. For m = 2, the instability is associated with the Ω_bin: Ω = 2: 1 resonance that occurs at r ≃ 1.6a. In this model, the m = 2 tidal field interacts with the m = 1 eccentric motions of the disk to produce disturbances of the form exp (3iϕ − 2iΩ_bint). These disturbances launch waves of that form at the 2:1 resonance. The waves in turn interact with the tidal field to produce a stress that increases the disk eccentricity exponentially in time.

However, the 2: 1 resonance lies at r ≃ 1.6a, inside the circumbinary disk gap. Furthermore, the gas motions are far from circular there. Even if the gas followed circular orbits in this region, the model of L91 predicts in this case an eccentricity growth rate ∼aΩ_bin/w_e, where w_e is the width of the eccentricity distribution and is ∼a. This very rapid growth rate is a consequence of the powerful m = 2 tidal potential for an equal-mass binary, the dominant component of the tidal field. However, this rate is more than an order of magnitude faster than the growth rate measured in the simulation. It may be possible that the effects of the resonance are still felt, but at a reduced level, near the disk inner edge as a consequence of finite width of the resonance. Thus, tidal instability may play a role in the eccentricity growth, but the evidence suggests that it is at most a limited role.

4.3.3. Eccentricity Growth through Stream Impact

As discussed in Section 3.3.1, the spacetime diagram reveals a special feature of the growth: the propagation of eccentricity from radii close to ∼2a toward larger radius. In this section, we examine the possibility that the impact of streams in the gap striking the inner edge of the disk is the mechanism of eccentricity injection into the disk.

During the exponential growth phase of eccentricity, a pair of gas streams are pulled in from the disk inner edge by the binary and then flung back to the disk after half binary period. It is easy to show that if the disk is on a circular orbit and free of eccentricity, the two streams are not capable of breaking the bisymmetry that is intrinsic for an equal-mass binary potential. We speculate that a small eccentricity breaks the symmetry by inducing unequal strength stream pairs (in terms of both the density and the velocity) and/or asynchronous stream–disk interaction. Once the symmetry of stream impact is broken, stream impact could amplify the initially small disk eccentricity, in turn increasing the asymmetry of the streams themselves. This process could then lead to sustained disk eccentricity growth.

We find two sorts of evidence that support this mechanism of eccentricity growth. The first is based on a special set of simulations we performed. In these, we eliminated part or all of the returning streams by enlarging the central cutoff. If the removal of outward streams halts eccentricity growth, we may take that as evidence in support of a stream impact origin for eccentricity excitation. We do not consider cutoff radii large enough to affect the circular orbit region of the disk because otherwise, the result might be equally well interpreted in terms of a resonance model (Section 4.3.2). Because these are essentially hydrodynamic effects and involve no vertical dynamics, these special simulations were performed in 2D with no magnetic field (see Section 2.5 and Table 1 for details of the simulations).

In the control run (B2D.rin=0.8), the disk is truncated around 2a, similar to B3D. In the gap region, there is a pair of streams at t = 300–700 while the disk eccentricity grows exponentially. A single strong stream appears when the disk eccentricity becomes significant. The simulation is terminated at t = 1500, and at that time the disk eccentricity is saturated. We find that the stream dynamics and the eccentricity growth in the control run are very similar to what were found in B3D, demonstrating that hydrodynamic effects alone can yield eccentricity qualitatively similar to MHD. Thus, the results based on hydrodynamic simulations can provide a clear indication of the role of the streams in eccentricity generation. We display the growth history of the disk eccentricity in Figure 17 (black solid). Here the disk eccentricity is calculated using a 2D definition of Equation (17). Similar to the MHD simulation, we find exponential growth from t ∼ 300 and 700. The growth rate in that period is ≃ 0.017Ω_bin, very close to that of the MHD simulation (≃ 0.018Ω_bin). We then restart from t = 500 with five different values of r_in while keeping all other parameters fixed. We terminate the reruns at t = 1500, the stopping time for B2D.ri=0.8.

Zoom In Zoom Out Reset image size

In Figure 17, we plot the evolution of disk eccentricities for cases with different r_in. The case with r_in = 1.0a (blue dash-dotted curve) loses only the tips of the streams, and its eccentricity evolution closely follows the control run, with exponential growth continuing without any interruptions until it becomes saturated. The disks with larger cutoffs (r_in = 1.3a: green dashed curve; 1.7a: red dotted curve), however, cease the exponential growth of their eccentricities right after the restart because most of the streams forming in the gap are lost in the hole and never get a chance to interact with the disk edge. There are also cases in between. When r_in = 1.1a (cyan triple-dot-dashed curve) and r_in = 1.2a (magenta long dashed curve), the growth rate of the eccentricity diminishes after the restart. It takes longer times for both disks to reach the eccentricity of the control run. We attribute these results to partial loss of outward stream motion.

When r_in = 1.3a, the simulation domain includes the region well inside the gap where gas motions are noncircular and are dominated by the streams. In this case, the simulation region also extends well inside the 2:1 resonance at r ≃ 1.6a, whose presence is required for the mechanism of Section 4.3.2. Consequently, any resonant instability should not be seriously affected by this inner boundary change. However, the eccentricity growth rate did change. This evidence then favors the stream impact mechanism.

The other evidence favoring a stream impact origin of the eccentricity comes from directly measuring the rate of change of eccentricity due to the stream–disk interaction in the MHD simulation. While the eccentricity grows exponentially, any asymmetry in stream impacts is necessarily small, so the effects of the gas stream are mild and difficult to directly measure. By the time t = 300, when the eccentricity growth is slow, the inner portion of the disk forms an eccentric ring between r ∼ 2a and 3a with e ∼ 0.08. Its pericenter is at r ∼ 2a and ϕ ∼ 3π/2, and the apocenter is at ∼3a and ϕ ∼ π/2 in the inertial frame. The eccentric ring slowly precesses around the binary. In this state there is only a single dominant stream. In this configuration, the binary interacts strongly only with the periastron side of the disk. At this stage, the effects of the stream are strong enough for us to reliably measure its effects on the disk.

The Gauss equations of celestial mechanics (Brouwer & Clemence 1961) provide a way to relate the perturbing effects of the stream to the evolution of disk eccentricity. This method has been previously used to calculate the effects of an external perturbing force on a ring in the context of dwarf-nova accretion disks (Lubow 1994). The evolution of the disk eccentricity is described by

$$\begin{equation} \frac{de}{dt} = \frac{\sqrt{1-e^2}}{a_{\rm d} n}\left[R\sin f + S\frac{4\cos f + 3e +e \cos (2f)}{2(1+e\cos f)}\right],\vspace*{3pt} \end{equation} \tag{ 23 }$$

where a_d is the semi-major axis of the eccentric ring of disk and $n=\ssty\sqrt{GM/a_{\rm d}^3}$ is the mean motion of the disk. Quantity f is the true anomaly where the stream impact occurs. R and S are the radial and angular components of the disturbance force density and are defined as

$$\begin{equation} (R,S) = \frac{\dot{m_s}(\mathbf {v_s} - \mathbf {u_e})}{m_e}, \end{equation} \tag{ 24 }$$

where $\dot{m_s}$ denotes the mass injection rate of the stream, m_e is the total mass of the eccentric ring, and $\mathbf {v_s}$ and $\mathbf {u_e}$ are the stream and disk velocities near the impact, respectively.

Since the impact takes place over a short time, we must use the higher time-resolution simulation (t = 300–322) in order to estimate the disk perturbation. We carry out our analysis using the vertically integrated 2D data. For the Gauss equation to apply, the properties of the eccentric ring, such as e, a_d, n, and the longitude of pericenter ϖ, need to be nearly constant over the orbital period of the ring. This assumption holds well because the timescale for the ring properties to change is much longer. In addition, the mass transferred to the ring by the stream impact over an orbital period is considerably smaller than the mass of the eccentric ring. We use time-averaged values for those parameters. The values are e = 0.08, a_d = 2.5a, and ϖ = 1.48 rad. The stream–disk impact is localized in space. We analyze the stream–disk impact over an arc that is defined by r = 2a and ϕ ∈ [ − π/4, π/4]. The quantities f and $\mathbf {v_s}$ are directly measured on the arc by taking the density-weighted averages at each time. The mass injection rate is calculated by integrating the mass flux along the arc. The instantaneous velocity of the eccentric edge $\mathbf {u_e}$ is obtained from the unaffected matter at a slightly greater radius, which we take to be an arc at r = 2.1a. In fact, the disk velocity is not very sensitive to the location we choose, so long as it is taken within the eccentric ring.

Using these parameters, we calculate the rate of change of the eccentricity from Equation (23). The result is presented in Figure 18. We find that each impact lasts for less than one time unit. The impacts induce peak values of de/dt that range between 4.0 × 10⁻³Ω_bin and 0.013Ω_bin at each half binary period. Most of the time, de/dt remains close to zero because stream–disk impacts are so brief. Taking time averages of the curve, we obtain an estimate of the impact-induced rate of change of eccentricity 〈de/dt〉_t ≃ 1.5 × 10⁻³Ω_bin. The corresponding growth rate γ_e ≃ 0.019Ω_bin is consistent with the growth rate at earlier times in the simulation (≃ 0.018Ω_bin). The growth rate contribution predicted by the stream impact at this later time may not be quite the same as the impact contribution at the earlier stage. Nonetheless, the approximate correspondence is encouraging. Indeed, one question that arises is why, given the estimate just made, the growth rate is smaller at this later time. We can only speculate that nonlinear effects, significant at this time but not earlier, may limit eccentricity growth.

Zoom In Zoom Out Reset image size

The disk after t = 400 shows a significant asymmetry that we call the "density lump." We believe the lump is due to the combined effect of stream impact and disk eccentricity. As discussed in Section 3.4, only one stream strongly interacts with the binary at any given time once the disk is noticeably eccentric. The stream moves outward after gaining angular momentum from the binary. Once it reaches the disk inner edge, the stream is shocked and its fluid compressed. Its mass is added to the nearby gas mass, increasing the local density. This region of enhanced density then moves around the binary at about the local orbital speed, decaying over roughly one orbital period due to pressure forces and shearing. One might not therefore expect such a small-scale density enhancement to become a large-scale, long-lasting lump.

We have, however, identified two mechanisms that sustain the lump and foster its growth. First, the lump can grow more concentrated via streams coming from the lump itself. We name this mechanism stream reabsorption. When the lump approaches pericenter, the binary peels off a gas stream that is denser than streams drawn from lower density regions. This relatively dense stream is then kicked back out by the binary torque. Although the azimuthal location reached by the returning stream may not be exactly where its starting point has arrived at this time, given the relatively large azimuthal extent of the lump, the chances that the returning stream strikes somewhere in the lump are relatively good. Moreover, the forward shock driven by the returning stream into the lump gas compresses it, restoring the density loss incurred by shearing during the time the stream passed near the binary. The left panel in Figure 19 shows the moment of an episode of stream reabsorption taking place.

Zoom In Zoom Out Reset image size

The second mechanism makes use of streams drawn from regions of the inner disk other than the lump. We call it lump feeding, as this channel of development not only enhances the concentration of the lump but also increases its mass. After t = 400, the density enhancement orbits eccentrically, and as it follows its orbit toward apocenter, both v_r and v_ϕ diminish. Meanwhile, streams traversing the gap continue to strike the disk inner edge, creating a new, smaller density enhancement at orbital phase angles behind the lump. The orbital speed of this new density concentration is much greater than that of the main lump because it is at rather smaller radius. As a result, the newly created and smaller lump can catch up with the main lump and join it before the lump returns from apocenter and accelerates. In the right panel of Figure 19, we show an example of this process, pictured at a time (t = 422) shortly before an outgoing stream reaches the lump. Three time units later, the small enhancement (the green region close to r = 2a at ϕ ≲ π/2) joins the main lump (green and red colored area in the second quadrant). The density-weighted velocity in the enhancement is (v_r, v_ϕ) ∼ (0.17, 0.66)Ω_bina (we average the stream gas where the surface density is greater than 1.5Σ₀). However, the lump has a slower velocity ∼(0.07, 0.52)Ω_bina (averaging over locations where Σ > 2Σ₀). As a result of both the lump feeding and stream reabsorption, the density lump grows in both mass and density contrast.

Another way to illustrate the connection between streams and the lump is to look at the spacetime diagram in Figure 12. In the diagram, the stream clearly stands out as the radially stretched features between r ≲ 1.5a and r ≳ 2a that propagate away from the binary with a pattern speed of about 0.2Ω_bina. This pattern repeats itself at a rate of about twice the binary frequency. Meanwhile, the zigzag oscillation at late time (t > 400) between radius ∼2a and 4a shows the movement of the density lump as it orbits around the binary on an eccentric orbit with a period of ∼30 time units. Clearly, when the lump passes the pericenter, the stream reabsorption process promotes growth in the density contrast, driving an increase of the mode strength at r = 2a–2.5a in the ascending part of the zigzag. The ascending part at even larger radius is able to maintain its strength via lump feeding. However, once the lump passes apocenter, the lump feeding is limited as the lump is now much farther away from the stream in azimuthal angle, and thus we find that the descending part at r ≳ 2.5a is less affected by the stream features. As the lump approaches pericenter again, another cycle of stream reabsorption and lump feeding begins. Note that the blue dots (low mode strength regions) around 2.5a between the zigzags actually show the moments that the newly kicked out stream reaches the radial coordinate (but is distant in azimuth) of the lump, so that it smooths the m = 1 mode at that radius. The zigzag feature's growth in both strength and radial range demonstrates that the density of the lump increases gradually, and as the disk becomes more eccentric, the semi-major axis of the lump's orbit grows.

In this work, we have performed the first 3D MHD simulation of a circumbinary disk around an equal-mass binary, which in this case follows a circular orbit. The main results are as follows:

• 1.  
  The disk exhibits a number of nearly steady features. There is a low-density gap within r ⩽ 2a and an eccentric inner disk at 2a–3a. At early times, there is a pair of gas streams that flow into the gap from the disk inner edge. They are nearly steady in the corotating frame of the binary. At later times, there is a single dominant stream whose mass flux is time varying with a period of half the binary orbital period. Parts of these streams are torqued so strongly by the binary that they return and impact on the disk inner edge.
• 2.  
  Some aspects of the disk evolve secularly. At late time, the disk inner edge develops an asymmetric density concentration ("the lump") whose mass, density contrast, and orbital eccentricity grow steadily. We find that the lump is due to a combination of the stream impact and disk eccentricity.
• 3.  
  The disk eccentricity grows exponentially during the time t = 100–250 (16–40 binary orbits) with a growth rate of ∼0.018Ω_bin. The growth rate then slows down significantly. By the end of the simulation, the disk body at 2a ≲ r ≲ 3a has reached an eccentricity ∼0.1. Its pericenter precesses slowly due to the quadrupolar component of the binary potential. Stream impact largely accounts for the eccentricity growth.
• 4.  
  The Reynolds stress associated with the streams is the leading transport mechanism in the gap region, while the Maxwell stress dominates in the disk body.
• 5.  
  Although the profile of the binary torque is similar to that of previous hydrodynamic simulations, its magnitude depends strongly on the magnitude of internal disk stresses. Normalized to the disk mass near the density peak, we find a measured total torque T(∞) ≃ 14 times greater than that found by MM08. Because the torque at any particular radius is proportional to the gas mass there, the integrated torque is directly proportional to the gas mass in the gap region, where the binary torques are strongest. This mass, always a small fraction of the disk mass, reaches the gap through the action of internal stresses in the accretion flow. Relative to the gas pressure, MHD stresses in the disk body are about an order of magnitude larger than the "viscous" stresses estimated phenomenologically in previous work; in the gap, this ratio increases by another order of magnitude. As a result, the density of matter in the gap is also about an order of magnitude larger than that found in hydrodynamic calculations, and this contrast leads directly to the larger torque.
• 6.  
  After time averaging, the accretion rate at the inner boundary is ∼30% of the peak rate in the disk body. Compared with MM08 and normalized to the peak surface density, the accretion rate at the inner boundary in this simulation is ∼40 times greater. This contrast, like the contrast in the total torque, can be largely attributed to the stronger internal stresses due to MHD effects. The time-dependent accretion rate is strongly modulated on both the binary orbital frequency (by the stream) and the inner-disk orbital frequency (by the lump).
• 7.  
  Previous work on the angular momentum budget of the binary has focused on the torque it exerts on the disk. Because we find a substantially larger accretion rate than previous calculations, we also find that the binary's gain in angular momentum due to accretion can substantially offset its loss by torque. As a result, the estimated binary contraction rate $\dot{a}/a \sim -8\times 10^{-4}(M_{\rm d}/M)(GM)^{1/2}a^{-3/2}$ is only slightly larger than the rate estimated by MM08.

These detailed results have a number of more general implications for the evolution of circumbinary disks, particularly in the context of supermassive black hole binaries. It is generally believed that stellar dynamical effects become relatively ineffective at driving the evolution of this sort of binary when its separation falls much below ∼1 pc (Begelman et al. 1980). Recently, there have been several attempts to solve this "final-parsec problem" with stellar dynamics. Although there are potential solutions to this problem, they require either special non-axisymmetric stellar distributions (Berczik et al. 2006; Khan et al. 2011; Preto et al. 2011) or extra perturbers such as giant molecular clouds (Perets & Alexander 2008). Given the uncertainty about whether these candidate mechanisms suffice to solve the problem, the prospect that angular momentum loss to a surrounding disk may push the binary through this barrier is an attractive one (e.g., Ivanov et al. 1999; Gould & Rix 2000; Armitage & Natarajan 2002; Escala et al. 2005; Kocsis & Loeb 2008; Schnittman & Krolik 2008; Cuadra et al. 2009; Dotti et al. 2009; Corrales et al. 2010; Tanaka et al. 2012; Farris et al. 2011).

We have previously described this process in terms of the orbital shrinkage rate $\dot{a}/a$. Defining the orbital shrinkage time $t_{\rm shrink}\equiv |a/\dot{a}|$, the time taken by torque alone to change the binary orbital angular momentum t_torque ≡ Mj_bin/T(∞), and the disk accretion time $t_{\rm acc} \equiv M_{\rm d}/\dot{M}$, we can rewrite this rate as

$$\begin{equation} \frac{t_{\rm acc}}{t_{\rm shrink}} = 2\left|-\frac{t_{\rm acc}}{t_{\rm torque}} + \frac{M_{\rm d}}{M}\left(\frac{j_{\rm in}}{j_{\rm bin}} - \frac{3}{2}\right)\right|. \end{equation} \tag{ 25 }$$

Using the values found in the simulation makes t_acc/t_shrink ≃ 0.6M_d/M, i.e., the ratio of the accretion time to the orbital shrinkage time is roughly the same as the ratio of the disk mass to the binary mass. Indeed, Lodato et al. (2009) argued that, unless the disk mass were at least comparable to the binary's secondary mass, it would be ineffective in driving binary evolution. Although it is true that a mass comparable to the secondary's mass must pass through the inner region of the disk over an orbital evolution time t_shrink, it is not necessarily true that this mass must be there all at once.

Suppose, for example, that M_d ≪ M. Because t_acc ≪ t_shrink in this case, if the disk mass is put in place once and for all, it would certainly be drained long before the binary has significantly evolved. On the other hand, one could also imagine situations in which the disk is continuously replenished. In that case, the binary orbit could change substantially due to this interaction even though the instantaneous disk mass is always much smaller than the binary mass; it is just that this process takes longer. Before determining how much time is required, it is important to recall that the close subtraction between the torque term and the accretion term may make the numerical value of (t_acc/t_shrink)/(M_d/M) rather sensitive to specific assumptions of our simulation, e.g., the binary mass ratio and the disk gas's equation of state. If the mass accretion rate were reduced by a factor of two while the torque was held fixed, the shrinkage rate would increase by a factor ≃ 4; conversely, if it were increased relative to a constant torque by a ratio ≳ 1.2, the binary orbit would actually expand over time.

These processes could also be affected by the fact that the accretion rate onto the binary is a substantial fraction of the accretion rate in the outer disk. The specific number we found (≃ 30%) is not terribly well defined because the lack of inflow equilibrium in the outer disk makes the denominator in that ratio very uncertain. However, a reinterpretation of that figure as, more vaguely, tens of percent nonetheless has significant implications. If the accretion rate in the outer disk were sufficient to supply an active galactic nucleus (AGN), even the fraction leaking through the disk's inner edge would still be large enough to fuel AGN activity, even if somewhat weaker. As a result, the inner edge of the disk would be illuminated in much the same way AGNs are known to illuminate the inner edge of their "obscuring tori" (Antonucci 1993). Low-density gas in the gap region would then be strongly heated by absorption of the AGN continuum, and much of it driven outward in a wind (Krolik & Begelman 1986; Balsara & Krolik 1993). Such a wind would lead to a reduction of the rate at which mass is captured by the members of the binary. As a result, the AGN luminosity would be smaller, and an equilibrium in which only a fraction of the accretion rate through the disk inner edge is captured by one or the other of the black holes might be achieved. An immediate consequence would be an increase in the ratio of torque to angular momentum captured by the binary, and therefore a factor of several increase in the shrinkage rate. Of course, a qualitative change in the equation of state of gas in the gap might well lead to an order-unity change in the torque, as well. At the same time, radiation forces associated with illumination of the disk body can thicken it (Pier & Krolik 1992; Krolik 2007; Shi & Krolik 2008). It is unclear how such a situation (effective vertical gravity weaker than zΩ²) would alter saturation of the MRI, but more rapid inflow is one plausible consequence.

The opposite case is also worth considering, in which M_d ≫ M, so that t_acc ≫ t_shrink. If this is so, then even a one-time deposit of mass in the disk might drive strong evolution in the binary orbit. However, when the binary separation decreases, the position of the strong torques moves to smaller distance from the center of mass as well, raising the question of how much mass might be there to receive those torques. We can estimate that time within the disk body of our simulation via t_inflow(r) ≡ ∫^rdr'/〈v_r〉_{ρ, t}, where 〈v_r〉_{ρ, t} represents the time- and density-weighted shell-averaged inflow velocity. In a steady-state disk around a point mass, t_acc would match t_inflow at r = r_p. However, to the degree that binary torques retard accretion, t_acc exceeds the inflow time measured at the surface density peak, which is t_inflow(r_p) ∼ 5 × 10²Ω⁻¹_bin. At sufficiently large radius, t_inflow(r) will nearly always be greater than t_acc because t_inflow(r)∝r^3/2(r/H)²; only a rapid outward flare in the disk thickness could alter this conclusion. It has long been thought that if t_inflow(r_p) ∼ t_acc ≫ t_shrink, the disk and the binary would decouple once the binary shrinks by a factor of a few because too little gas would be found close enough to the binary to feel the torques; if orbital compression depends solely on interaction with a circumbinary disk, the shrinkage time could therefore never be much shorter than t_inflow(r_p) (e.g., Gould & Rix 2000; Armitage & Natarajan 2002).

However, the dynamical behavior seen both in our simulation and in others hints that this assumption might not be valid. The streams occurring in our simulation travel inward on a timescale ∼Ω⁻¹_bin, which is ≪t_inflow. Although a smaller binary separation would make the trajectories of the streams at radii not far inside the disk edge more nearly angular momentum-conserving, and therefore make it more difficult for them to penetrate to much smaller radii, the internal magnetic stresses within the streams might be strong enough to replace the binary torques. Indeed, we found that the nominal "α" in the gap region in our simulation was one or two orders of magnitude larger than in the disk body. There is also another mechanism, commonly seen in global MHD simulations of disks around point masses, that might contribute in this way. As emphasized by Guan et al. (2011), it is easy for disks far from inflow equilibrium to drive mass-inflow rates much greater than those expected on the basis of inflow equilibrium estimates. For example, when an annular disk with a seed poloidal magnetic field is created, orbital shear rapidly creates toroidal field from radial field. The −B_rB_ϕ stress acting on the low-density matter near the disk's inner edge rapidly removes its angular momentum, propelling it inward. As the binary separation becomes small compared to r_p, the gravitational potential comes to resemble a point mass, and similar dynamics might be expected. Because the exertion of a sizable torque requires only a small fraction of the disk mass to orbit in the gap region, where the torque per unit mass is greatest, a strong torque might continue even though the bulk of the disk mass remains far from the binary.

A massive disk is also likely to be subject to self-gravity, a mechanism not treated in our simulation. If M_d/M ≳ (c_s/v_orb), self-gravitational instability can be triggered (Goldreich & Lynden-Bell 1965). The evolution of this instability depends on the ratio of the cooling time to the dynamical timescale. If the ratio is greater than unity, the disk can maintain marginal stability, and the self-gravity would add extra stress to facilitate the accretion (Gammie 2001; Lodato & Rice 2004, 2005). If the cooling time is short, which is likely to be the case for disks around massive black holes, the disk is likely to fragment and much of its mass may be transformed into stars (Shlosman & Begelman 1989; Gammie 2001; Nayakshin 2006), reducing the gas surface density (Lodato et al. 2009). The consequences for interaction with the binary remain somewhat unclear. In sufficiently massive disks, stars might form with orbits taking them sufficiently close to the binary that the summed torques from many individual star–binary encounters could continue to shrink the binary; in effect, the stellar loss-cone is repopulated by local star formation (Cuadra et al. 2009). In less massive disks, stars might form only in the disk body. If most of the disk mass is converted into stars, the gas mass available to fill the gap would be reduced, possibly leading to both a smaller torque and a smaller accretion rate. On the other hand, gravitational torques due to non-axisymmetric fluctuations in the stellar density might create stresses strong enough to overcome the reduction in gaseous disk mass (Binney & Tremaine 2008). Given all these complicated possibilities, it is clear that further study of self-gravitating circumbinary disks will be necessary before any strong conclusions can be drawn on the fate of massive disks.

The work done by the binary on the streams can heat the inner edge of the disk through stream impact. The total energy per unit time delivered by the binary torque can be estimated as

$$\begin{equation} L_{\rm torq} \equiv \Omega _{\rm bin}T(\infty)\simeq 6.5\times 10^{-4} \Omega _{\rm bin}^3 a^2 M_{\rm d}. \end{equation} \tag{ 26 }$$

For disks around massive binary black holes,

$$\begin{equation} L_{\rm torq} \sim 1.9 \times 10^{38} M_6^{3/2} a_{0.1}^{-1/2} N_{24}\, \hbox{ erg\, s$^{-1}$}, \end{equation} \tag{ 27 }$$

where M₆ = M/10⁶ M_☉ is the binary mass, a_0.1 = a/0.1 pc is the separation, and N₂₄ = N_H/1.7 × 10²⁴ cm⁻² is the column density of the disk. Thus, the expected luminosity is well below a typical AGN luminosity and would therefore be difficult to detect unless its power were emitted in a small number of lines. Similarly, we can get the heating rate for disks around stellar binaries:

$$\begin{equation} L_{\rm torq} \sim 7.9 \times 10^{33} M_0^{3/2} a_{10}^{-5/2} M_{\rm d,-1} \,{\rm erg} \,{\rm s}^{-1}, \end{equation} \tag{ 28 }$$

where M₀ = M/M_☉ is the mass of the binary in solar masses, a₁₀ = a/10 AU is the binary separation in units of 10 AU, and M_{d, -1} = M_d/0.1 M_☉ is the disk mass in units of 0.1 M_☉. We note that, with our ideal MHD model, this luminosity only sets the upper limit of the possible heating rate due to the binary torque. In reality, disks outside 10 AU might possess low-ionization dead zones (Gammie 1996) that suppress the MHD turbulence and diminish the mean effective α_SS. Based on the accretion rate observed in disks around young stars, which is typically ∼10⁻⁸ M_☉ yr⁻¹(e.g., Hartigan et al. 1995; Gullbring et al. 1998) for T Tauri stars, we expect that the actual luminosity will be two to three orders of magnitude smaller than estimated in Equation (28). Because stream impact is periodic with the period of the binary, the resulting luminosity might be modulated with a period ≃ 16M^−1/2₀a₁₀^3/2 yr.

In the black hole binary context, circumbinary disk studies have been in part motivated by a search for mechanisms to solve the "last parsec problem," the expected slow-down in orbital evolution by stellar encounters when the binary separation becomes ≲ 1 pc. Our results have, in one sense, quantitatively weakened the constraints on this proposed solution: per unit disk mass, our nominal result is an orbital shrinkage rate a few times faster than previously thought. However, in a number of other ways, we have raised potential complications. One is that the nominal shrinkage rate is very sensitive to parameters because the larger accretion rate per unit disk mass we find leads to a near cancellation between the binary's loss of angular momentum through torques on the circumbinary disk and acquisition of angular momentum through accretion. Another stems from the likely back-reaction on both the disk proper and the matter in the gap caused by AGN illumination as a result of accretion onto the binary. Still another is the question of how the asymmetry in the disk (the "lump") affects dynamics driven by disk self-gravity when the disk is sufficiently massive. For all these reasons, our work has enriched, rather than settled, the question of whether, or to what degree, circumbinary disks alleviate the "last parsec problem."

Finally, we acknowledge that, as the first MHD simulation of a circumbinary disk with order-unity mass ratio, this work is limited in several aspects.

Most of all, it would be interesting to know the long-term behavior of the lump. The reason why we ended the simulation when we did was that matter accumulates more rapidly than the magnetic field in the lump. Consequently, its internal Alfvén speed diminished and the ability of the simulation code to resolve adequately the MHD turbulence along with it. Such a situation leads to an artificial weakening of the magnetic field. To follow properly what actually happens in the lump will require significantly better spatial resolution.

There are also numerous different parameters of the disks that need to be explored, such as the mass ratio and the disk thickness. For example, a mass ratio different from unity might result in noticeable changes in both the structure of the disk and the leakage rate because the secondary would come closer to the disk inner edge. A smaller thickness would likely reduce the inflow rate in the disk body. Such a change might diminish the ratio of gas mass in the gap to disk mass, but the considerations discussed above about MHD dynamics specific to the gap might counteract this effect.

We have also assumed that the disk orbits precisely in the binary orbital plane, and in a prograde sense. Inclined orbits might well occur, and their interaction with the binary's quadrupole potential will lead to precession of the inner disk's orbital plane. Coplanar, but retrograde, gas orbits eliminate the resonance structure characterizing the linear theory of binary–disk interaction and could lead to other qualitative effects (Nixon et al. 2011).