THREE-DIMENSIONAL MHD SIMULATION OF CIRCUMBINARY ACCRETION DISKS. II. NET ACCRETION RATE

The model setup is very similar to the one described in Shi et al. (2012). Our account here is therefore very brief and emphasizes the few differences between our new simulations and the one presented in our earlier paper. We also summarize the main properties of all runs in this paper in Table 1 for reference.

Table 1.  Properties of Accretion Disk Simulations

Label	Type of Simulation	q	Resolutiona	Radial Extentb	$\langle \dot{M}({r}_{\mathrm{in}}){\rangle }_{t}$ c $(\sqrt{{GMa}}{{\rm{\Sigma }}}_{0})$
S2DE	Hydrodynamics	0.0	256 × 64	$(0.8,40)$	...
S3DEQ	MHD	0.0	$512\times 400\times 96$	$(0.8,40)$	0.0081
S3DE	MHD	0.0	$512\times 400\times 384$	$(0.8,40)$	0.0085
B3DE	MHD	1.0	$512\times 400\times 384$	$(0.8,40)$	0.011–0.014
B3DEq	MHD	0.1	$480\times 400\times 384$	$(1.02,40)$	0.013–0.017

Notes.

^aCell counts are listed as $r\times \theta \times \phi$ for 3D MHD simulations and $r\times \theta$ for 2D hydrodynamic simulations. ^bIn units of a. The azimuthal extent is $(0,2\pi )$ for all MHD simulations except S3DEQ, which is $(0,\pi /2)$. ^cTime1averaged accretion rate measured at the inner boundary. The accretion rate in the quarter-disk simulation S3DEQ is multiplied by 4. For binary runs, time averages at both early and late stages are presented and separated by a dash.

Download table as:  ASCIITypeset image

We construct our disks in an inertial frame with origin at the center of mass. We set the gravitational constant G and the central mass M to be unity, whether it is a single object or a binary. The binary separation a = 1 is the simulation lengthscale.³ The time unit is therefore $\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}=\sqrt{{a}^{3}/({GM})}$. The density is normalized to the initial midplane value ${\rho }_{0}=1$ and the surface density unit is defined to be ${{\rm{\Sigma }}}_{0}={\rho }_{0}a$. Similar to Shi et al. (2012), we adopt a globally isothermal equation of state with a fixed sound speed, but we set ${c}_{{\rm{s}}}=0.1$, twice that of the previous paper. A larger sound speed allows us to achieve better resolution and also perform longer duration simulations. Again, we assume the disk mass is considerably smaller than the mass of the central object and therefore neglect the self-gravity of the disk.

As discussed in Shi et al. (2012, see Section 2.2), the simulation grid needs to resolve three different length scales: the disk scale height H, the maximum growth-rate wavelength of the MRI ${\lambda }_{\mathrm{MRI}}=8\pi /\sqrt{15}{v}_{{\rm{A}}}/{\rm{\Omega }}(R)$, and the spiral density wavelength ${\lambda }_{{\rm{d}}}\sim 2\pi {c}_{{\rm{s}}}/{{\rm{\Omega }}}_{\mathrm{bin}}$, where v_A is the Alfvén speed and ${\rm{\Omega }}(R)$ is the disk rotational frequency. We constructed our grid in spherical coordinates following the same scheme as in Shi et al. (2012), which was originally proposed by Noble et al. (2010): logarithmic in the radial direction, uniformly spaced in azimuthal angle, and spaced according to a polynomial function in polar angle in order to concentrate cells near the orbital plane (Equation (1) of Shi et al. (2012), with $\xi =0.825$, ${\theta }_{c}=0.1$ and n = 9). There were $[{N}_{r},400,384]$ cells in $(r,\theta ,\phi )$, covering a computational domain spanning $[{r}_{\mathrm{in}},{r}_{\mathrm{out}}]$ radially, $[{\theta }_{c},\pi -{\theta }_{c}]$ meridionally, and $[0,2\pi ]$ azimuthally, where ${r}_{\mathrm{out}}=40a$, ${\theta }_{c}=0.1$, ${r}_{\mathrm{in}}$ is the radius of the circular central cut out ($0.8a$ for q = 0 and q = 1, increased to $1.02a$ for $q=0.1$ to confine the secondary inside the excision), and ${N}_{r}=512$ for q = 0 and q = 1, diminished to 480 for $q=0.1$ run as a result of the larger ${r}_{\mathrm{in}}$. In terms of minimum cell size, our polar angle resolution is a factor ∼1.4 better than in Shi et al. (2012). With a doubled disk scaleheight (${c}_{{\rm{s}}}=0.1$ versus ${c}_{{\rm{s}}}=0.05$), we therefore have about twice as many cells per scaleheight as in Shi et al. (2012). In recent years, standards have been developed for achieving resolution adequate to describe the principal features of MRI-driven MHD turbulence: at least 10–20 cells per $2\pi {v}_{{\rm{A}}z,\phi }/{\rm{\Omega }}$, where ${v}_{{\rm{A}}z,\phi }$ is the Alfvén speed associated with the vertical (z) component of the magnetic field or azimuthal (ϕ) component (Hawley et al. 2011, 2013). We typically have $\simeq 20$ cells per characteristic vertical wavelength and at least 10 cells per characteristic azimuthal wavelength.

We choose strict outflow boundary conditions for the inner and outer radial surfaces and also at the edges of the polar axis cutout; all inward velocities are set to zero on these boundaries. Periodic boundary conditions are used for the azimuthal direction. 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.

2.2.1. Single Mass Run: q = 0

The single mass simulation both serves as a reference point and provides the initial conditions for the circumbinary simulations discussed later. In order to speed up these lengthy simulations, we built a three-step ladder.

Step 1. 2D Hydrodynamic Disk: We start from a 2D axisymmetric, inviscid, hydrodynamic simulation (S2DE) in the $(r,\theta )$ plane with the same numerical and physical parameters as described above, except that the parameter ξ in the θ-grid definition is 0.92 and the number of cells is only lower $256\times 64$ in $r\times \theta$. The goal of this 2D simulation is to get rid of initial transients and to find a quasi-equilibrium disk solution. The initial configuration follows that of Shi et al. (2012, Section 2.3): the disk stretches from $r=3a$ to $r=6a$ with constant midplane density ${\rho }_{0}$ and is vertically hydrostatic. We evolve the disk for $2000\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$, when the disk reaches a steady state (see Figure 1(a)).

Zoom In Zoom Out Reset image size

Step 2. 3D MHD Quarter Disk: We then use the final dump from S2DE as the initial conditions for a $\pi /2$ wedge of a 3D MHD disk (S3DEQ). The initial disk for S3DEQ is therefore axisymmetric. Its initial density and angular velocity are interpolated from the results of S2DE data. Motions in other directions are neglected as they are very small. The initial contours of magnetic vector potential ${A}_{\phi }$ are a set of nested poloidal loops following the contours of the density within the main body of the disk (the contour lines in Figure 1(a)). The values of the ${A}_{\phi }$ contours are ${A}_{0}(\rho -0.1{\rho }_{0})$, where A₀ is a constant determined by requiring the average plasma $\beta =100$. The magnetic field ${\boldsymbol{B}}$ is computed by taking the curl of the vector potential. Run S3DEQ begins at t = 0 and lasts until $t=1330\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$. A snapshot of the quarter-disk at $t=800\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ is shown in Figure 1(b). We find the quarter disk approaches steady accretion for $r\lt 5a$ between $t=800\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ and $1300\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ (see Figure 2).

Zoom In Zoom Out Reset image size

Step 3. 3D MHD full Disk: We then patch together four identical copies of the $\pi /2$ disk from S3DEQ at $t=800\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ to build the initial conditions for a full $2\pi$ disk (S3DE). This run continues from $t=800\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ through $1340\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$. The evolution history (see Figure 2), nearly the same as the quarter disk run, indicates that by $\sim 800\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ the full disk undergoes quasi-steady accretion.⁴ We show a clipped isosurfaces plot of this fully turbulent disk at $t=1000\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ in Figure 1(c). The disk at this moment serves as the initial conditions for the binary runs that are discussed next.

2.2.2. Binary Runs: $q\ne 0$

The initial disk adjusts to the new potential within 1–2 binary periods after the equal-mass binary replaces the point-mass. During this initial transient phase, the binary clears out a cavity around itself in which the density is ${10}^{-3}$–${10}^{-2}$ of what it was when there was a point-mass at the center. A small amount of mass initially found at $r\lesssim 2a$ passes through the inner boundary (about $2\%$ of the total disk mass), but most of the disk mass within $r\lt 2a$ is pushed outward by the binary torque. After $\gtrsim 50\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$, the circumbinary disk gradually reaches a quasi-steady state. In its inner portion ($r\lesssim 4a$), the mass interior to r becomes nearly constant in time (first panel in Figure 3) so that the surface densityʼs radial profile also changes only very slowly (Figure 4). We also note that the late-time average surface density of the outer disk follows the ${\rm{\Sigma }}\propto {r}^{-2}$ scaling observed in Shi et al. (2012) quite well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After the first $\sim 50\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$, the accretion histories show that the q = 1 disk has also reached a quasi-steady state with respect to this property (Figure 5). Clearly, the accretion is not hampered by the central binary (compare the red dashed curve for the q = 0 disk with the green curve for q = 1). Although the strong binary torque clears out a low density cavity near the center, the cavity is not empty of gas. Streams of gas still flow through the potential maxima at the L2 and L3 points (see Figure 10) and maintain accretion at a rate comparable to the single mass case. The overall time-averaged accretion rate of the q = 1 disk is in fact a bit greater than that of the q = 0 case: $\dot{M}{(r={r}_{\mathrm{in}})\simeq 0.011({GMa})}^{1/2}{{\rm{\Sigma }}}_{0}$ over t = 1050–$1250\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ ("early-time") and $\simeq 0.014{({GMa})}^{1/2}{{\rm{\Sigma }}}_{0}$ over 1250–$1400\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ ("late-time"), compared to $0.0085{({GMa})}^{1/2}{{\rm{\Sigma }}}_{0}$ of the q = 0 disk (from $t=1000\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ until the end of the simulation). These numbers translate to an accretion efficiency $\epsilon \simeq 1.3$–1.6. This result is in approximate agreement with the viscous hydrodynamics simulations of D'Orazio et al. (2013), in which they found the accretion rate of an equal mass binary is ∼0.93 times the single-mass case (see their Table 3). It is in even better agreement with the results of Farris et al. (2014), who found a ratio of 1.55.

Zoom In Zoom Out Reset image size

We show the radial dependence of the time-averaged accretion rate $\dot{M}(r)$ in Figure 5. For both early (t = 1050–$1250\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$) and late time averages (t = 1250–$1400\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$), $\dot{M}(r)$ for both q = 1 and $q=0.1$ exceeds that of the q = 0 disk at all radii. The q = 1 disk achieves a fairly steady accretion within $r\lesssim 4a$, quite similar to a single-mass disk. At late times, the accretion rate at smaller radii increases by $\sim 20\%$, and the zone of nearly constant accretion radius as a function of radius extends outward from $r\simeq 4a$, characteristic of early times, to $r\simeq 25a$.

The different levels of accretion observed in Figure 5 are directly connected to variations in the internal stresses. As shown in Figure 6, although the Maxwell stress changes little with time, the Reynolds stress rises by a factor of ∼2 after $t=1250\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ at all radii outside $r\simeq 2a$. This increase in internal stress then drives a larger accretion rate at late time.

Zoom In Zoom Out Reset image size

The increase in Reynolds stress occurs at the same time ($\sim 1250\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$) as the diskʼs spiral density waves change from a tightly wrapped m = 2 pattern to a single-armed wave (Figure 7). Although the spiral waves have little effect on the Maxwell stress, the vertically integrated Reynolds stress is enhanced along the density crest of the m = 1 wave. For instance, the region of large Reynolds stress between $r=4a$ and $8a$ at $t=1300\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ (the yellow spiral arm winding from 9 o'clock to 6 o'clock) is absent in the $t=1200\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$ snapshot. It appears that the single-armed spiral wave is more effective at conveying angular momentum outward than the two-armed wave.

Zoom In Zoom Out Reset image size

The accretion flow in the $q=0.1$ disk is quite similar to that seen in the q = 1 case. The disk reaches its quasi-steady state after the first $\sim 50\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$. The history of enclosed disk mass within a given radius in Figure 3 indicates a slowly evolving steady state in the inner part of the circumbinary disk. Further evidence that an approximate state of inflow equilibrium has been reached for $r\lt 4a$ is given by the time-averaged radial profiles of $\dot{M}(r)$ (Figure 5). Also like the q = 1 case, the accretion rate for $q=0.1$ gradually increases over the course of the simulation, rising by $\simeq 30\%$.

These two phases of accretion are also closely connected to the phase transition of the disk structure for the $q=0.1$ run. In Figure 8, we show snapshots of disk properties before and after the transition. It is clear that the disk evolves from a relatively compact two-armed spiral structure into a large scale single-armed structure. We notice that this transition happens faster than in the q = 1 run, possibly due to the asymmetry of the binary system itself. A common feature shared by these two phases is the gas stream attached to the secondary, which indicates stronger accretion onto the secondary than onto the primary.

Zoom In Zoom Out Reset image size

Like the q = 1 case, the accretion rate through the inner boundary in the $q=0.1$ case (black solid curve in upper panel of Figure 5) appears to have two phases as well. At early times (t = 1050–$1250\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$), the time-averaged rate $\dot{M}{(r={r}_{\mathrm{in}})\simeq 0.013({GMa})}^{1/2}{{\rm{\Sigma }}}_{0}$, while the accretion rate gradually increases so that the average at late times (t = 1250–1500) is $\simeq 0.017{({GMa})}^{1/2}{{\rm{\Sigma }}}_{0}$. In both phases, the time-averaged $\dot{M}(r={r}_{\mathrm{in}})$ exceeds the q = 1 case by $\sim 20\%$, and is ∼1.5–2 times the single mass case.⁵ By comparison, in the 2D hydrodynamic simulations of D'Orazio et al. (2013), the time-averaged accretion rate onto a $q=0.1$ binary was $\simeq 0.7\times$ that of a single-mass system, while those of Farris et al. (2014) found a ratio $\simeq 1.7$. Just as for the q = 1 case, our results show a qualitatively similar but quantitatively somewhat greater accretion efficiency than found by D'Orazio et al. (2013), and quite good quantitative agreement with Farris et al. (2014).

In both the q = 1 and $q=0.1$ disks, the fractional amplitude of fluctuations in the mass accretion rate is a few times greater than for a point-mass (q = 0). As we will argue below, the larger fluctuations are driven by the binary itself. However, the accretion rate fluctuations in the binary cases also appear to depend significantly on mass-ratio, both in amplitude and in frequency. The typical peak-to-trough amplitude contrast for $q=0.1$ is $\simeq 0.01$–0.015 ${{\rm{\Sigma }}}_{0}{({GMa})}^{1/2}$, about as large as the mean accretion rate. This is about twice the amplitude seen when q = 1.

At early times in the q = 1 run, the accretion rate fluctuates periodically at a frequency of $\sim 1.5\;{{\rm{\Omega }}}_{\mathrm{bin}}$, as shown in Figure 9. This is the stage, called the "transient state" by D'Orazio et al. (2013), in which the disk is beginning to become elliptical and two streams of nearly equal strength run inward from the diskʼs inner edge; D'Orazio et al. (2013) similarly found a periodic modulation with this frequency. Later in the simulation, the $1.5\;{{\rm{\Omega }}}_{\mathrm{bin}}$ peak in the power spectrum splits into two lower magnitude spikes, their frequencies centered on the original peak frequency but separated by $\simeq 0.1\;{{\rm{\Omega }}}_{\mathrm{bin}}$. This split indicates a change in the disk structure from point-symmetric toward more eccentric shape.

Zoom In Zoom Out Reset image size

Unlike the q = 1 case, the $q=0.1$ case can induce m = 1 disk asymmetry immediately. Consequently, right from the start we see strong peaks in the accretion rate power spectrum at $\simeq 0.7\;{{\rm{\Omega }}}_{\mathrm{bin}}$ and $\simeq 1.3\;{{\rm{\Omega }}}_{\mathrm{bin}}$. These modes are beats between the binary frequency ${{\rm{\Omega }}}_{\mathrm{bin}}$ and the orbital frequency of the small density enhancement near the diskʼs inner edge, $\simeq 0.3\;{{\rm{\Omega }}}_{\mathrm{bin}}$. Another peak at $\simeq 1.4\;{{\rm{\Omega }}}_{\mathrm{bin}}$ might be the first harmonic of the $0.7\;{{\rm{\Omega }}}_{\mathrm{bin}}$ beat. At later time, we still see the $0.7\;{{\rm{\Omega }}}_{\mathrm{bin}}$ and $1.4\;{{\rm{\Omega }}}_{\mathrm{bin}}$ beat frequencies, but they shift slightly toward higher frequency as the disk gap expands by a small amount. Interestingly, we do not see the strong single peak at $1.0\;{{\rm{\Omega }}}_{\mathrm{bin}}$ observed by D'Orazio et al. (2013) for this mass ratio.

In the body of an accretion disk, the inflow speed is generically much slower than the orbital speed, ∼α ${(H/r)}^{2}{v}_{\mathrm{orb}}$, where α is the usual ratio of vertically integrated stress to vertically integrated pressure, and H is the local scale height. On the other hand, this flow, although only $\sim H$ thick in the vertical direction, takes place, on average, around the entire circumference of the disk, through an area $2\pi r$ wide.

By contrast, the flow across the cavity (see Figure 10) is restricted to very narrow streams. Along the central density maximum of the streams, they are typically ∼2–$3H$ wide if measured sideways from the maximum to where the density drops by $90\%$. Moreover, the density in the streams as they approach the inner boundary is ∼3–10 times lower than the density in the disk body. Thus, in order to carry the same mass inflow, the inward velocity in a stream must be $\sim (r/H)({\rho }_{\mathrm{disk}}/{\rho }_{\mathrm{stream}})$ larger than the typical disk inflow speed, or $\sim ({\rho }_{\mathrm{disk}}/{\rho }_{\mathrm{stream}})\alpha (H/r){v}_{\mathrm{orb}}$. This condition can be easily achieved because the absence of stable closed orbits within $r\sim 2a$ when q is not too small leads to a characteristic inward speed $\sim 0.3{v}_{\mathrm{orb}}$, whereas $({\rho }_{\mathrm{disk}}/{\rho }_{\mathrm{stream}})\alpha (H/r)\sim 0.03$–0.1 in the conditions of our simulation, and often smaller in real disks.

Zoom In Zoom Out Reset image size

Thus, one way of looking at the contrast between the 1D and the 2D/3D results is to observe that the 1D picture, in which there is no inward motion within $r\sim 2a$, is correct at almost all, but not quite all, positions. We explore this idea further in the next subsection.

The reason gas accretion occurs through such narrow channels is that only a small fraction of the gas possesses the proper initial conditions (in position–velocity phase space) for "infall" trajectories, i.e., trajectories that begin from the diskʼs inner edge and reach the inner cutoff ${r}_{\mathrm{in}}$. Most of the gas near the diskʼs inner edge that begins moving inward is turned around and flung back out to the disk by binary torques. This fate of the majority of the initially inflowing gas is what led the 1D analysis to predict no net accretion.

To test this explanation and identify exactly what the special conditions for successful inward flow are, we begin by solving a large number of restricted three-body problems further restricted to orbits entirely in the midplane. The binary for these trajectory integrations has a circular orbit and q = 1. The test-particles have initial conditions evenly distributed within a phase-space volume defined by $1.5a\leqslant r\leqslant 2.5a$, $0\leqslant \phi \leqslant 2\pi$, $-0.4\sqrt{{GM}/a}\leqslant {v}_{r}\leqslant 0$, and $-0.4\leqslant {\omega }^{\prime }-[{{\rm{\Omega }}}_{K}(r)-{{\rm{\Omega }}}_{\mathrm{bin}}]\leqslant +0.4$. Here ${\omega }^{\prime }$ is the angular frequency in a frame co-rotating with the binary. The initial specific angular momentum is therefore $j\equiv ({\omega }^{\prime }+{{\rm{\Omega }}}_{\mathrm{bin}}){r}^{2}$.

Our results are shown in the left panel of Figure 11 in the form of a parallel coordinates plot.⁶ When the initial radius is $r=2a$ (at which a circular orbit requires $j\simeq 1.45\sqrt{{GMa}}$ when q = 1), particles whose initial j is $\gtrsim 1.5\sqrt{{GMa}}$ (${\omega }^{\prime }\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}\gtrsim -0.63$) cannot reach the inner boundary; particles with $1.3\lesssim j/\sqrt{{GMa}}\lesssim 1.5$ ($-0.68\lesssim {\omega }^{\prime }\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}\lesssim -0.63$) can, but only if they also have ${v}_{r}\lesssim -(0.1$–$0.3)\sqrt{{GM}/a}$; particles with $j\lesssim 1.3\sqrt{{GMa}}$ are able to travel to ${r}_{\mathrm{in}}$ even with v_r only slightly negative. The behavior of both these classes of particles can be understood by reference to the approximate (i.e., ignoring the quadrupolar contribution) effective potential ${V}_{\mathrm{eff}}(r)\simeq -{GM}/r+{j}^{2}/(2{r}^{2})$ at $r={r}_{\mathrm{in}}$. If ${v}_{r}^{2}\ll {V}_{\mathrm{eff}}({r}_{\mathrm{in}})$, the particle radial kinetic energy is negligible, and infall can happen only when ${V}_{\mathrm{eff}}({r}_{\mathrm{in}})\leqslant 0$, i.e., $j\simeq 1.3\sqrt{{GMa}}$ or less, or equivalently ${\omega }^{\prime }$ is no more than $\simeq -0.68\;{{\rm{\Omega }}}_{\mathrm{bin}}$. Alternatively, if v_r is sufficiently negative, the particle radial kinetic energy can be large enough to overcome a positive effective potential barrier.

Zoom In Zoom Out Reset image size

These conditions for infall are not easily met in a real disk because the mean j is close to the circular orbit value, $\simeq 1.45\sqrt{{GMa}}$, and radial infall speeds $\lesssim -(0.1$–$0.3)\sqrt{{GM}/a}$ are rare in the disk body. Consequently, only those few fluid elements with j well below the mean can contribute to the inflow, and they are found in a tightly constrained portion of phase space. In the right panel of Figure 11, we show a selection of fluid elements taken from a snapshot of our simulation at $t=1208\;{{\rm{\Omega }}}_{\mathrm{bin}}^{-1}$. The samples are drawn from a ring of disk around $r=2a$ with a width of $0.5a$ in order to capture the inner edge of the disk. Comparing the orange regions (all cells) with the blue ones (those falling in), one can conclude that indeed only a tiny fraction of the disk proper falls to the binary, while the rest of its fluid elements are pushed back out even if they initially move in a small distance. The criterion governing which fluid elements are able to move inward is identical to that identified for the test-particles: $j\lesssim 1.3\sqrt{{GMa}}$. In contrast, most of the orange fluid elements have either too large a j (or ${\omega }^{\prime }$), i.e. large ${V}_{\mathrm{eff}}$, or too small an inflow velocity $-{v}_{r}$, and are therefore flung out. Particularly low j fluid elements move inward so quickly there is too little time for binary torques to substantially raise their angular momentum; because they can pass the effective potential barrier at ${r}_{\mathrm{in}}$, they cross the inner boundary and are permanently removed from the circumbinary disk.

Because gas flows near the binary are close to ballistic (Lubow & Artymowicz 2000, Shi et al. 2012), we can estimate the accretion rate by counting the number of fluid elements satisfying the inflow criterion. The data shown in the right panel of Figure 11 indicate that the blue fluid elements cover $\sim 1/10$ of the annulus between $r\sim 1.9a$ and $2.2a$. For a typical surface density $\sim 0.2{{\rm{\Sigma }}}_{0}$ and an infall timescale $\sim 2\pi /{{\rm{\Omega }}}_{\mathrm{bin}}$, the inferred inflow rate would be $\dot{M}\sim 0.012\sqrt{{GMa}}{{\rm{\Sigma }}}_{0}$, consistent with our simulation results. Thus, although 90% of inner disk fluid elements cannot go any significant distance inward, the angular momentum distribution is broad enough that its low j tail suffices to carry the full accretion rate.

Having found that the fluid elements able to accrete are defined by their low specific angular momentum, the next question to answer is how their angular momentum is reduced to that level. Ordinary MHD turbulence does not broaden the angular momentum distribution to this degree: as shown in Figure 6, the Reynolds stress in the q = 0 case is only $\simeq 1\%$ of the pressure, and the pressure is only $\sim 1\%$ of the orbital energy per unit mass. Instead, the answer appears to be a consequence of the binary torques themselves. As previously remarked, most of the mass in the streams that move inward from the circumbinary disk returns to the disk after its specific angular momentum is raised by those torques. With that additional angular momentum, its azimuthal velocity is somewhat greater than the local orbital velocity when it strikes the disk ($j\simeq 1.6{({GMa})}^{1/2}$ as opposed to $j\simeq 1.45{({GMa})}^{1/2}$). The work done by the torques also increases the streams' energy, giving them an outward radial speed $\simeq (0.3$–$0.5){({GM}/a)}^{1/2}$.

An example of such a rapidly moving outward stream can be seen near $r\sim 2.3$, $\phi \sim 1.5$ in Figure 12. Part of the stream ($r\simeq 2.3$–2.4, $\phi \simeq 1$–1.5) is still heading outward. However, a part of that stream has already encountered the ridge of high density trending from larger radius to smaller between $\phi \simeq 2$ and $\phi \simeq 3$. When it struck that ridge, a pair of shocks was formed, a forward shock propagating into the disk material and a reverse shock propagating into the material that had been moving outward. In the reverse shock, a portion of the stream mass loses a significant fraction of both its angular momentum and its orbital energy and heads back inward, seen in this snapshot along the line from ($r\simeq 2.3$, $\phi \simeq 2$) to ($r\simeq 1.8$, $\phi \simeq 3$). This shock deflection is the origin of the low angular momentum material which can now travel all the way to the inner boundary.

Zoom In Zoom Out Reset image size

Because this mechanism depends almost entirely on 2D orbital mechanics, its nature should be only weakly dependent on parameters such as the diskʼs aspect ratio $H/r={c}_{{\rm{s}}}/{v}_{\mathrm{orb}}$, provided only that ${c}_{{\rm{s}}}/{v}_{\mathrm{orb}}$ is small enough that the stream–disk interaction is supersonic. In our simulation, in which the sound speed is $0.1{({GM}/a)}^{1/2}$, the Mach numbers of these shocks are ∼3–5; in real disks, they could be considerably larger. Support for the view that the stream–disk interaction is at most weakly affected by the disk thickness also comes from the fact that 2D simulations employing laminar hydrodynamics and an "α viscosity" to mock up the fluidʼs internal shear stress (Farris et al. 2014) find a very similar continuity in the accretion rate. It is possible, however, that the dynamics of the stream–disk shocks may depend upon the equation of state of the shocked gas. In our simulation, we assumed that the gas is isothermal. If, instead, this sort of shock were to occur in a less lossy gas, the immediate post-shock temperature would be much greater than the disk gas temperature. The shocked gas could then swell vertically well beyond the thickness of the disk, permitting it to flow above and below the disk. Investigation of such effects is well beyond the scope of our effort here.

As we have already remarked, a surprising outcome of our simulations is that at late times a strong one-armed spiral wave propagates outward through the circumbinary disk, enhancing the Reynolds stress sufficiently to increase the accretion rate by a few tens of percent. Its origin is worth further attention.

The most likely cause of this m = 1 feature is the complement of the mechanisms studied in the previous subsection. There we focused on the properties of those special fluid elements able to travel all the way from the inner edge of the disk to the binary. Here we focus on all the others, the ones that may initially move inward, but then feel the torques exerted by the binary and are propelled outward again until they strike the diskʼs inner edge. To demonstrate how these streams excite spiral waves, we fit the wave form seen in our simulations (surface density plots in Figures 7 and 8) with a standard density wave pattern (e.g., Equation (36) of Rafikov 2002). The resulting pattern speed is $\simeq 0.17\;{{\rm{\Omega }}}_{\mathrm{bin}}$ ($\simeq 0.21\;{{\rm{\Omega }}}_{\mathrm{bin}}$) for the q = 1 ($q=0.1$) case, which suggests wave excitation occurs at $r\simeq 3.2a$ when q = 1 and $r\simeq 2.8a$ when $q=0.1$. As shown in Figures 7 and 8, the density is strongly enhanced near $r\simeq 3a$. This connection is most clearly seen in an animation (Figure 13).

The animations further show that the point of impact, the azimuthal location of the spiral wave driving point, rotates around the diskʼs inner edge at an angular frequency slightly lower than the binary frequency.

Similar to tidally or mass-transfer induced spiral shocks that transmit angular momentum outward in circumstellar disks in a close binary system (e.g., Sawada et al. 1986; Rózyczka & Spruit 1993), the waves excited by the binary-driven streams can also provide long-range angular momentum transport. Figure 14 shows the net outward angular momentum flux (AMF) associated with the time averaged Reynolds stress before and after large scale m = 1 density waves are excited (see the green curves in the bottom two panels). Compared to the early period (top row), the Reynolds stress contributes a sizable negative torque at disk radius ∼3–$4.5a$ (2.5–$4a$) at late times (bottom row) in the q = 1 ($q=0.1$) disk, which could explain the slightly greater accretion rate observed in $q\ne 0$ disks than in the q = 0 disk.

Zoom In Zoom Out Reset image size

These large scale one-armed spiral density features were not previously reported in Shi et al. (2012), D'Orazio et al. (2013), or Farris et al. (2014). We believe they are due to the greater global isothermal sound speed ($0.1{{\rm{\Omega }}}_{\mathrm{bin}}a$) used here, twice the sound speed in simulation B3D (Shi et al. 2012). The faster sound speed causes spiral waves to stretch farther radially, creating large-scale loosely wrapped density features (Savonije et al. 1994). In D'Orazio et al. (2013; as well as in MacFadyen & Milosavljević 2008), $H/r=0.1$ at all radii, so that the isothermal sound speed drops $\propto \;{r}^{-1/2}$. The wavelength becomes shorter as the wave moves outward, causing the wave to become more and more tightly wrapped as it travels to greater radius. Given this apparent sensitivity to the disk equation of state, it remains to be seen whether such strong spiral waves are generic or rare in real circumbinary disks.