3D MHD Simulations of Circumbinary Disks

Three Dimensional MHD Simulation of Circumbinary Accretion Disks –2. Net Accretion Rate

Ji-Ming Shi11affiliation: Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544 and Julian H. Krolik22affiliation: Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218

When an accretion disk surrounds a binary rotating in the same sense, the binary exerts strong torques on the gas. Analytic work in the 1D approximation indicated that these torques sharply diminish or even eliminate accretion from the disk onto the binary. However, recent 2D and 3D simulational work has shown at most modest diminution. We present new MHD simulations demonstrating that for binaries with mass ratios of 1 and 0.1 there is essentially no difference between the accretion rate at large radius in the disk and the accretion rate onto the binary. To resolve the discrepancy with earlier analytic estimates, we identify the small subset of gas trajectories traveling from the inner edge of the disk to the binary and show how the full accretion rate is concentrated onto them as a result of stream-disk shocks driven by the binary torques.

accretion, accretion disks — binaries: general — MHD — methods: numerical

1 Introduction

Label Type of Simulation ResolutionaaCell counts are listed as for 3D MHD simulations and for 2D hydrodynamic simulations. Radial ExtentbbIn units of . The azimuthal extent is for all MHD simulations except S3DEQ, which is .
S2DE Hydrodynamics 0.0
S3DE MHD 0.0
B3DE MHD 1.0
B3DEq MHD 0.1
Table 1: Properties of Accretion Disk Simulations

Binary systems often form within a gaseous disk environment. Because the tidal torques exerted by the binary on the disk are repulsive when the disk and the binary orbit in the same sense, it has long been thought that these torques severely limit, or perhaps entirely prevent, accretion from the disk onto the binary Pringle (1991). As a result, a very low-density cavity forms within of the binary center of mass, where is the binary’s semi-major axis. Despite this prediction, observations indicate accretion onto such binary systems at a level comparable to their single star counterparts. This is clearly detected for low-mass binary stars in nearby star-forming regions (e.g., White & Ghez, 2001). High velocity gas flows are observed bridging gaps that are believed to be cleared out by protoplanetary companions (e.g., Casassus et al., 2013; Rosenfeld et al., 2014). There are even detections of accretion onto planetary mass companions (e.g., Bowler et al., 2011; Zhou et al., 2014). A growing number of dual AGN candidates have also been reported (e.g., Komossa et al., 2003; Rodriguez et al., 2006; Comerford et al., 2011); if systems like these become mutually bound, they could become accreting binaries. It is therefore important to understand how gas is able to accrete despite these torques, both to be able to use electromagnetic signatures (Rödig et al., 2014) as diagnostics and also to learn how accretion influences the evolution of these systems.

Much work has already been expended investigating how the binary torque reshapes the surrounding circumbinary disk (e.g., Artymowicz & Lubow, 1994; Crida et al., 2006; MacFadyen & Milosavljević, 2008; Shi et al., 2012; Noble et al., 2012). However, how the torque regulates the gas accretion is a less well-developed subject. The 1D analysis of Pringle (1991) has been extended, but adopting that paper’s assumption of zero accretion (e.g., Milosavjević & Phinney, 2005; Lodato et al., 2009; Liu & Shapiro, 2010). More recently, there have been a number of 2D and 3D simulations of such systems (e.g., Artymowicz & Lubow, 1996; Bryden et al., 1999; MacFadyen & Milosavljević, 2008; Dotti et al., 2009; Shi et al., 2012; Noble et al., 2012; D’Orazio et al., 2013; Farris et al., 2014). These consistently find the cavity predicted analytically, but also accretion at a rate comparable to or somewhat smaller than in the outer disk. Typically, the accretion takes place in narrow, high velocity streams emanating from the edge of the cavity and spiraling inward toward the central binary. Moreover, these studies have shown significant consequences of the accretion for the development of the system. In protoplanetary disks, the amount of accretion through the gap controls the mass growth of a gap-clearing planet (Bryden et al., 1999; Lubow et al., 1999; D’Angelo et al., 2006) and the surface density within the gap (Zhu et al., 2011; Fung et al., 2014). The angular momentum associated with the advecting gas also strongly influence the orbital evolution of the binary (Shi et al., 2012; Roedig et al., 2012).

It is therefore crucial to quantify the net accretion fraction () of a prograde circumbinary disk, i.e. the ratio between the time-averaged net accretion rate and the rate that would take place in the absence of binary torque :


For the extreme small mass-ratio case (we define the mass-ratio ), previous studies using viscous hydrodynamics simulations have found that there is a threshold mass-ratio , such that for smaller the accretion rate is nearly unchanged when compared to that for a single point-mass case; i.e. when . For values of closer to unity, the results so far are somewhat mixed. Most simulations have found that is reduced by a factor of a few (Bryden et al., 1999; Lubow et al., 1999; Shi et al., 2012; Noble et al., 2012; D’Orazio et al., 2013), but the dependence on remains poorly defined. For example, D’Orazio et al. (2013) saw diminished accretion for binaries, while Farris et al. (2014) saw no signs of suppression for . The situation is further clouded by the fact that previous work (except for Shi et al. (2012) and Noble et al. (2012)) has used the approximation of 2D hydrodynamics, in which internal accretion stresses, both within the accretion disk itself and in the streams traversing the cavity, are described by the phenomenological viscosity model, rather than the actual physical mechanism, correlated MHD turbulence.

In this paper, we carry out 3D MHD simulations to investigate accretion in circumbinary disks around equal mass () and unequal mass () binaries. We assume the disk and binary are coplanar, orbit in the same sense, and the binary orbit is a circle with fixed semi-major axis . For a control experiment, we also simulate their single point-mass counterpart. Comparison between the single mass run and the circumbinary runs provides the answers to two key questions: (1) To what extent does the accretion rate of a binary differ from that of a single object, i.e., what is the value of ? (2) How does the accretion rate from circumbinary disks depend on the mass ratio of the binary, i.e., how does vary with ? As we will show, the answer to the first question is also the answer to the second: for all . This answer leads to a further question that we will also attempt to answer: what was missing from the early 1D analytic studies that led them to conclude ?

We organize this paper as follows: In § 2, we describe the physical model and numerical setup of our single point-mass and circumbinary disk simulations. In § 3, we present our simulation results. We analyze these results physically in § 4. Finally, we summarize our conclusions and discuss the implications of our findings in § 5.

Figure 1: (a) Density distribution (color) at the end of the 2D hydrodynamical simulation S2DE, which is used to provide initial conditions for the 3D MHD quarter disk run S3DEQ. Note that the vertical scale is stretched relative to the horizontal. The black contour lines mark the initial magnetic field loops. (b)  Density isosurfaces (color) of the bottom half disk of S3DEQ at . This data is used as the initial conditions of S3DE after combining four identical copies to make the whole disk. (c)   Density isosurfaces (color) of the bottom half disk of S3DE at . This data becomes the initial conditions for our circumbinary simulations.

2 Numerical Simulations

In this section, we discuss details about the numerical experiments we carried out in order to explore the effects of the binary torques on the accretion rate. We used the same code as Shi et al. (2012), which is a modern version of the 3D, time-explicit Eulerian finite-differencing ZEUS code for MHD (Stone & Norman, 1992a, b; Hawley & Stone, 1995).

2.1 Model Setup

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.

We construct our disks in an inertial frame with origin at the center of mass. We set the gravitational constant and the central mass to be unity, whether it is a single object or a binary. The binary separation is the simulation lengthscale111In the simulation, has no physical interpretation other than the code-unit of length.. The time unit is therefore . The density is normalized to the initial midplane value and the surface density unit is defined to be . Similar to Shi et al. (2012), we adopt a globally isothermal equation of state with a fixed sound speed, but we set , 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 , the maximum growth-rate wavelength of the MRI , and the spiral density wavelength , where is the Alfven speed and 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 , and ). There were cells in , covering a computational domain spanning radially, meridionally and azimuthally, where , , is the radius of the circular central cut-out ( for and , increased to for to confine the secondary inside the excision), and for and , diminished to for run as a result of the larger . In terms of minimum cell size, our polar angle resolution is a factor better than in Shi et al. (2012). With a doubled disk scaleheight ( vs. ), 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 , where is the Alfven speed associated with the vertical () component of the magnetic field or azimuthal () component (Hawley et al., 2011, 2013). We typically have  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 Numerical Experiments

Figure 2: Accretion history of the single mass quarter-circle simulation S3DEQ (black) and the full-circle S3DE MHD (red) disks. Top: Inner edge accretion rate as a function of time. Middle: Time-averaged accretion rate as a function of radius. Bottom: Mass interior to a given radius as a function of time. The accretion rates and enclosed mass for the quarter-disk are multiplied by for comparison with whole disk values. Note the very similar behaviors between the quarter disk and the complete disk, in agreement with previous simulations: the time averaged accretion of MHD turbulent disks is nearly independent of the extent of the azimuthal domain if it is more than (Hawley, 2000; Papaloizou & Nelson, 2003).

2.2.1 Single Mass Run:

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 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 in . 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 to with constant midplane density and is vertically hydrostatic. We evolve the disk for , when the disk reaches a steady state (see Figure 1(a)).

Step 2. 3D MHD Quarter Disk: We then use the final dump from S2DE as the initial conditions for a 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 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 contours are , where is a constant determined by requiring the average plasma . The magnetic field is computed by taking the curl of the vector potential. Run S3DEQ begins at and lasts until . A snapshot of the quarter-disk at is shown in Figure 1(b). We find the quarter disk approaches steady accretion for between and (see Figure 2).

Step 3. 3D MHD full Disk: We then patch together four identical copies of the disk from S3DEQ at to build the initial conditions for a full disk (S3DE). This run continues from through . The evolution history (see Figure 2), nearly the same as the quarter disk run, indicates that by the full disk undergoes quasi-steady accretion222The time averaged accretion rate at different radii suggests our disk is approaching inflow equilibrium for - (see the second panel in Figure 2). Outside that radius, the accretion rate gradually shrinks, changing sign to outflow at . A small amount of mass outflow at large radius carries the angular momentum transported from inside.. We show a clipped isosurfaces plot of this fully turbulent disk at 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:

We performed two binary simulations, one with mass ratio (B3DE) and one with (B3DEq). For both runs, we use S3DE data at as the initial conditions. Both also retain the domain size, numerical resolution, and physical parameters of the run; the only change is to replace the point-mass at the center with either an equal-mass binary or a binary with the same total mass. In both binary cases, the orbit is constant and circular, and rotates counter-clockwise, prograde with respect to the disk.

Everything about the case is the same except for the extent of the radial grid. Because the secondary is from the center of mass, we increase the inner excision size from to . This change results in a truncation of the initial disk of the single mass run, but we keep all other aspects, such as the resolution and outer boundary of the domain, intact.

3 Results

Figure 3:  History of disk mass interior to a given radius for (top) and circumbinary disks.

3.1 Run

The initial disk adjusts to the new potential within 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 of what it was when there was a point-mass at the center. A small amount of mass initially found at passes through the inner boundary (about of the total disk mass), but most of the disk mass within is pushed outward by the binary torque. After , the circumbinary disk gradually reaches a quasi-steady state. In its inner portion (), the mass interior to becomes nearly constant in time (first panel in Figure 3) so that the surface density’s radial profile also changes only very slowly (Fig. 4). We also note that the late-time average surface density of the outer disk follows the scaling observed in Shi et al. (2012) quite well.

Figure 4: Surface densities for the initial condition (red dashed) and the time-averaged (black) and (green) cases. Solid curves are early-time, dashed curves are late-time averages. Solids are for the early time averages and dashed for the late time averages. For reference, a scaling is shown by the black dotted curve.
Figure 5: Top: Accretion rates through the inner boundary of (red dashed), (black), and (green) disks. Bottom: Time-averaged accretion rates for (red dashed), (black), and (green) as a function of radius. For the cases, both early (solid) and late (dashed) time averages are presented.

After the first , the accretion histories show that the disk has also reached a quasi-steady state with respect to this property (Fig. 5). Clearly, the accretion is not hampered by the central binary (contrast the red dashed curve for the disk with the green curve for ). 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 and points (see Fig. 10) and maintain accretion at a rate comparable to the single mass case. The overall time-averaged accretion rate of the disk is in fact a bit greater than that of the case: over (‘early-time’) and over (‘late-time’), compared to of the disk (from until the end of the simulation). These numbers translate to an accretion efficiency . 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 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.

We show the radial dependence of the time-averaged accretion rate in Figure 5. For both early (-) and late time averages (), for both and exceeds that of the disk at all radii. The disk achieves a fairly steady accretion within , quite similar to a single-mass disk. At late times, the accretion rate at smaller radii increases by , and the zone of nearly constant accretion radius as a function of radius extends outward from , characteristic of early times, to .

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 after at all radii outside . This increase in internal stress then drives a larger accretion rate at late time.

The increase in Reynolds stress occurs at the same time () as the disk’s spiral density waves change from a tightly-wrapped pattern to a single-armed wave (Fig. 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 wave. For instance, the region of large Reynolds stress between and at (the yellow spiral arm winding from 9 o’clock to 6 o’clock) is absent in the snapshot. It appears that the single-armed spiral wave is more effective at conveying angular momentum outward than the two-armed wave.

Figure 6: Time averaged stress to pressure ratios (top) and their associated angular momentum fluxes for the disk. Different line styles denote different time spans: - (solid), - (dashed), and - of the run (dotted); different colors distinguish Reynolds (green) and Maxwell (red) stresses. There is an increase of Reynolds stress and its corresponding angular momentum flux at later times, owing to the spiral waves.

3.2 Run

The accretion flow in the disk is quite similar to that seen in the case. The disk reaches its quasi-steady state after the first . 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 is given by the time-averaged radial profiles of (Figure 5). Also like the case, the accretion rate for gradually increases over the course of the simulation, rising by .

These two phases of accretion are also closely connected to the phase transition of the disk structure for the 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 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.

Figure 7: Four disk diagnostics before and after the spiral wave phase transition at in the simulation. Upper row: , lower row: . From left to right, the variables shown are: surface density, density weighted radial velocity (normalized to the sound speed), vertically integrated Maxwell and Reynolds stresses in units of . The two white dots in the central cut-out region represent the binary members. White indicates off-scale high values.

Figure 8: The same diagnostics as in Figure 7, but for the disk at and . In this case, the spiral wave transition took place at , slightly earlier than in the equal mass case.

Like the case, the accretion rate through the inner boundary in the case (black solid curve in upper panel of Figure 5) appears to have two phases as well. At early times (), the time-averaged rate , while the accretion rate gradually increases so that the average at late times () is . In both phases, the time-averaged exceeds the case by , and is times the single mass case 333Here we have ignored the effects from the different as they are negligibly small. The accretion rate of the equal mass binary hardly changes between and . . By comparison, in the 2D hydrodynamic simulations of D’Orazio et al. (2013), the time-averaged accretion rate onto a binary was that of a single-mass system, while those of Farris et al. (2014) found a ratio . Just as for the 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).

3.3 Accretion rate fluctuations

In both the and disks, the fractional amplitude of fluctuations in the mass accretion rate is a few times greater than for a point-mass (). 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 is , about as large as the mean accretion rate. This is about twice the amplitude seen when .

At early times in the run, the accretion rate fluctuates periodically at a frequency of , 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 peak in the power spectrum splits into two lower magnitude spikes, their frequencies centered on the original peak frequency but separated by . This split indicates a change in the disk structure from point-symmetric toward more eccentric shape.

Figure 9: Fourier decomposed power spectrum of the accretion through the inner boundary of a binary (top panel) and a binary (bottom panel).

Unlike the case, the case can induce disk asymmetry immediately. Consequently, right from the start we see strong peaks in the accretion rate power spectrum at and . These modes are beats between the binary frequency and the orbital frequency of the small density enhancement near the disk’s inner edge, . Another peak at might be the first harmonic of the beat. At later time, we still see the and 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 observed by D’Orazio et al. (2013) for this mass ratio.

4 Analysis

Having seen that is actually slightly greater than unity, we now turn to an effort to understand this result. The first point to raise is that it is unlikely can persist for a long period of time. If it were to do so, the inner region of the circumbinary disk () would be drained of mass, inevitably leading to a reduction in the accretion rate onto the binary. Thus, the better way to think about the values of seen in our simulations, a few tens of percent greater than unity, is that the spiral waves excited in a circumbinary disk by the members of the binary create a sufficient enhancement of the Reynolds stress to raise the accretion rate per unit mass in the inner disk by a few tens of percent. By this means, an accretion rate equal to that injected at large radius can be sustained by a surface density somewhat smaller than required when the potential is due to a point-mass. Over longer times than we can follow with this kind of simulation, we expect that the surface density in the inner disk will decline to this level, leaving the disk in true inflow equilibrium.

With that clarification, it is time to consider the question of why 2D and 3D simulations consistently see substantial accretion from circumbinary disks onto the central binary despite the contrary prediction made by 1D studies. One clue to the answer comes from the structure of the accretion flow through the cavity: narrow streams.

4.1 Stream Structure

Figure 10: Left column: Time averaged midplane density (top), midplane (middle) and inner boundary (bottom) accretion rate , both for the binary over the last of the simulation. All figures in a frame comoving with the binary. Right column: same as left, but for binary. Here negative means inflow. The plus symbols in the midplane plots mark the L2 and L3 points. Summed separately, regions of inward and outward mass flux have comparable magnitude; their net, although smaller in magnitude, is consistently inward.

In the body of an accretion disk, the inflow speed is generically much slower than the orbital speed, , where is the usual ratio of vertically-integrated stress to vertically-integrated pressure, and is the local scale height. On the other hand, this flow, although only thick in the vertical direction, takes place, on average, around the entire circumference of the disk, through an area wide.

By contrast, the flow across the cavity (see Fig 10) is restricted to very narrow streams. Along the central density maximum of the streams, they are typically wide if measured sideways from the maximum to where the density drops by . Moreover, the density in the streams as they approach the inner boundary is 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 larger than the typical disk inflow speed, or . This condition can be easily achieved because the absence of stable closed orbits within when is not too small leads to a characteristic inward speed , whereas in the conditions of our simulation, and often smaller in real disks.

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 , is correct at almost all, but not quite all, positions. We explore this idea further in the next subsection.

4.2 Orbital properties of inward trajectories

Figure 11: Left: A parallel coordinate plot for test particles at distributed uniformly in , (shown in units of , and . The last coordinate is the initial specific angular momentum in units of . The green lines show the initial locations of all infalling particles; note that all have , the specific angular momentum that supports a circular orbit in the region sampled. The red lines show particles with enough angular momentum to trace nearly-circular orbits, but also fall inward; all have especially large inward speed as well. The red particles concentrate at two opposite azimuth angles where the binary potential is relatively weaker than at other angles, making it easier for particles to fall in.
                                                                                                                                                  Right: Similar to the left but now the lines represent individual cells drawn from the actual simulation of the disk at . The orange lines show all cells within annulus having ; the blue lines show all the infalling cells; the green lines show infalling cells with initial radius exactly . Comparing the orange and blue regions, only a small fraction of the disk cells actually fall in. The green cells found here are consistent in their properties with the green test particles shown on the left: for .

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 . 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 . The test-particles have initial conditions evenly distributed within a phase-space volume defined by , , , and . Here is the angular frequency in a frame co-rotating with the binary. The initial specific angular momentum is therefore .

Our results are shown in the left panel of Figure 11 in the form of a Parallel Coordinates plot444Parallel Coordinates is a common method for multidimensional visualization. For our case, each particle is represented by a line connecting multiple attributes shown as vertical coordinates in parallel.. When the initial radius is (at which a circular orbit requires when ), particles whose initial is () cannot reach the inner boundary; particles with () can, but only if they also have ; particles with are able to travel to even with 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 at . If , the particle radial kinetic energy is negligible, and infall can happen only when , i.e., or less, or equivalently is no more than . Alternatively, if is sufficiently negative, the particle radial kinetic energy can be large enough to overcome a positive effective potential barrier.

These conditions for infall are not easily met in a real disk because the mean is close to the circular orbit value, , and radial infall speeds are rare in the disk body. Consequently, only those few fluid elements with 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 . The samples are drawn from a ring of disk around with a width of 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: . In contrast, most of the orange fluid elements have either too large a (or ), i.e. large , or too small an inflow velocity , and are therefore flung out. Particularly low 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 , 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 of the annulus between and . For a typical surface density and an infall timescale , the inferred inflow rate would be , 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 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 Fig. 6, the Reynolds stress in the case is only of the pressure, and the pressure is only 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 ( as opposed to ). The work done by the torques also increases the streams’ energy, giving them an outward radial speed .

Figure 12: Surface density (color contours) and mass-weighted vertically-averaged velocity in the orbital plane (arrows) at .

An example of such a rapidly-moving outward stream can be seen near , in Figure 12. Part of the stream (–2.4, –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 and . 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 (, ) to (, ). This shock deflection is the origin of the low angular momentum material which can now travel all the way to the inner boundary.

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 , provided only that is small enough that the stream-disk interaction is supersonic. In our simulation, in which the sound speed is , the Mach numbers of these shocks are –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.

4.3 One-armed Spiral Wave

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.

\pdfmark[]pdfmark=/ANN,Subtype=/FileAttachment,Raw=/F 0/T (f13.mp4)/Contents (Media File (video/mp4))/AP ¡¡/N¡¡¿¿/R¡¡¿¿/D¡¡¿¿¿¿/FS filespec1\pdfmark[]pdfmark=/ANN,Subtype=/Screen,Border=0 0 0,Raw=/_objdef screenannot1/F 5/T (f13.mp4)/Contents (Media File (video/mp4))/BS ¡¡/S/S/W 0¿¿/P ThisPage/AA aadict1/AP ¡¡/N¡¡¿¿/R¡¡¿¿/D¡¡¿¿¿¿\pdfmarkpdfmark=/PUT,Raw=screenannot1 ¡¡/A ¡¡/R mediarendition1/S/Rendition/OP 0/JS ( app.focusRect=true;if(focusonplayer==undefined)var focusonplayer=0; var settings=privateData: paused: false , autoPlay: false, visible: false, volume: 100 , showUI: true, startAt: 0; var events=new app.media.Events(onBlur: function (e) if(focusonplayer ¿ 0)focusonplayer=0; , afterBlur: function (e) if(focusonplayer==0)try e.target.settings.privateData.paused=false; e.target.play(); catch(e) , onFocus: function (e) focusonplayer=1; , afterFocus: function (e) if(!e.target.isPlaying)try e.target.settings.privateData.paused=false; e.target.play(); if(!e.target.isPlaying)if( e.target.settings.startAt.time —— e.target.settings.startAt.frame —— e.target.settings.startAt.marker —— e.target.id == ’vnd.adobe.swname:AAPL_QuickTime’ ) e.target.seek(e.target.settings.startAt); else e.target.stop(); e.target.play(); catch (e) e.target.visible=true;, onPlay: function (e) e.target.settings.privateData.paused=false; , onPause: function (e) e.target.settings.privateData.paused=true; , afterReady: function (e) try if( e.target.settings.startAt.time —— e.target.settings.startAt.frame —— e.target.settings.startAt.marker ) e.target.play(); e.target.pause(); e.target.stop(); e.target.settings.privateData.paused=false; e.target.seek(e.target.settings.startAt); e.target.visible=true; e.target.settings.privateData.paused=false; e.target.play(); catch (e) ); var player1=app.media.openPlayer(settings: settings, events: events ); ) /AN screenannot1¿¿¿¿\pdfmarkpdfmark=/OBJ,Raw=/type/dict/_objdef mediarendition1\pdfmarkpdfmark=/PUT,Raw=mediarendition1 ¡¡/C mediaclipdata1/S/MR/SP ¡¡/BE¡¡/O 0.0¿¿¿¿/P ¡¡/BE ¡¡/F 2/C true/D ¡¡/S /F¿¿/A false¿¿¿¿¿¿\pdfmarkpdfmark=/OBJ,Raw=/_objdef mediaclipdata1/type/dict\pdfmarkpdfmark=/PUT,Raw=mediaclipdata1 ¡¡/D filespec1/P ¡¡/TF(TEMPACCESS)¿¿/S/MCD/CT (video/mp4)¿¿\pdfmarkpdfmark=/OBJ,Raw=/_objdef filespec1/type/dict\pdfmarkpdfmark=/PUT,Raw=filespec1 ¡¡/F(f13.mp4)/Type/Filespec¿¿\pdfmarkpdfmark=/OBJ,Raw=/_objdef fstream1/type/stream\pdfmarkpdfmark=/PUT,Raw=fstream1(f13.mp4) (r) file\pdfmarkpdfmark=/PUT,Raw=fstream1 ¡¡/Type/EmbeddedFile/Subtype(video/mp4)¿¿\pdfmarkpdfmark=/PUT,Raw=filespec1 ¡¡/EF ¡¡ /F fstream1 ¿¿¿¿\pdfmarkpdfmark=/OBJ,Raw=/_objdef pageopenaction1/type/dict\pdfmarkpdfmark=/PUT,Raw=pageopenaction1 ¡¡/R mediarendition1/S/Rendition/OP 2/JS ( tryif(player1.isOpen)player1.page=this.pageNum; player1.visible=true; elsethrow ’isClosed’; catch(e)) /AN screenannot1¿¿\pdfmarkpdfmark=/OBJ,Raw=/_objdef pagecloseaction1/type/dict\pdfmarkpdfmark=/PUT,Raw=pagecloseaction1 ¡¡/R mediarendition1/S/Rendition/OP 1/JS ( tryplayer1.settings.privateData.paused=false; if(!player1.isPlaying) player1.play(); player1.stop(); if( player1.settings.startAt.time —— player1.settings.startAt.frame —— player1.settings.startAt.marker —— player1.id == ’vnd.adobe.swname:AAPL_QuickTime’ ) player1.seek(player1.settings.startAt); focusonplayer=-1; player1.visible=false; catch(e) )/AN screenannot1¿¿\pdfmarkpdfmark=/OBJ,Raw=/type/dict/_objdef aadict1\pdfmarkpdfmark=/PUT,Raw=aadict1 ¡¡/PO pageopenaction1/PC pagecloseaction1¿¿

Figure 13: Streams returning to the circumbinary disk exciting spiral waves. Click the figure to play a short animation of surface density in the and cases 666Also available at http://www.astro.princeton.edu/~jmshi/leakage.htm.. Spiral waves can clearly be seen to begin when a stream pushed outward by the binary torques strikes the inner edge of the disk.

The most likely cause of this 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 Figure 7 and 8) with a standard density wave pattern (e.g., Equation (36) of Rafikov, 2002). The resulting pattern speed is () for the () case, which suggests wave excitation occurs at when and when . As shown in Figure 7 and 8, the density is strongly enhanced near . This connection is most clearly seen in an animation (Fig. 6).

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 density waves are excited (see 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 () at late time (bottom row) in the () disk, which could explain the slightly greater accretion rate observed in disks than in the disk.

Figure 14: Radial derivatives of time-averaged and shell-integrated angular momentum flux (AMF) as functions of radius for the (left) and (right) cases before and after the phase transition, where red shows the differentiated Maxwell AMF and green represents the differentiated Reynolds AMF. Negative means outward transport of angular momentum. Time averages are long enough that the net change of local angular momentum are close to zero. The late time averages of the () cases show that the Reynolds stress dominates the angular momentum budget at (), a location coinciding with the stream impact regions shown in Figure 7 for the case and Figure 8 for the case.

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 () 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)), at all radii, so that the isothermal sound speed drops . 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.

5 Conclusions

In this paper, we carried out 3D global MHD simulations of circumbinary disks with binary mass ratios and and contrasted them with a disk orbiting a solitary point-mass. We found two major results:

  1. The time-averaged accretion rate from a circumbinary disk with either or is indistinguishable from that of a circum-solitary disk whose central mass is the same. In other words, essentially all the mass supply given the disk at large radius ultimately leaves its inner edge and travels to the binary. The similarity of the and cases in this regard suggests that this result depends at most weakly on binary mass-ratio. This result confirms, with physical internal fluid stresses, the conclusion reached by Farris et al. (2014) on the basis of 2-d hydrodynamics and a phenomenological viscous stress.

  2. The key reason why initial 1D analyses suggested that the accretion efficiency (the parameter we call ) is rather than is that they omitted consideration of the small volume in orbital phase space from which trajectories can travel inward from the disk’s inner edge all the way to the binary, avoiding the strong torques that, for most of phase space, push matter back outward. Only those fluid elements with specific angular momentum less than the circular orbit value at the disk’s inner edge can cross the gap and reach the binary. Ironically, fluid elements with exactly this property are created as a consequence of the binary torques themselves: these torques add enough angular momentum to other gas traveling through the gap to propel it back out to the disk; it shocks upon reaching the disk, and a portion of its mass is deflected onto accretion orbits.

Our results have several observational implications. If all the matter supplied to a circumbinary disk at large radius accretes onto the binary through narrow streams, those streams must shock when they strike the outer edges of accretion disks around the members of the binary. In the context of the formation of stellar binaries, these shocks may be the sites of the H vibrational lines that can often be detected (Beck et al., 2012). Radial inflow streams at velocities approaching free-fall can also potentially explain both the observed low density cavities in transitional disk systems and the relatively normal stellar accretion rates they maintain (Rosenfeld et al., 2014). In the context of supermassive black hole binaries, one can similarly expect strong shocks where the streams strike the outer edges of the individual disks. Hard X-rays from those shocks may be observable when the binary separation is small enough (Roedig et al., 2014; Farris et al., 2015). If the accretion flow fed in at large radius is large enough, still greater luminosity can be generated when the accreted matter approaches the two black holes’ ISCO regions.

We thank an anonymous referee for constructive questions that led to improvement of this paper. We also thank Eugene Chiang and Neal Turner for useful discussions. This research was supported by an allocation of advanced computing resources provided by the National Science Foundation. The computations were performed on Kraken at the National Institute for Computational Sciences (http://www.nics.tennessee.edu/). This work was partially supported by NSF grant AST-1028111 (JHK) and NASA Origins grant NNX13AI57G (J-M S).


  • Artymowicz & Lubow (1994) Artymowicz, P. & Lubow S.H. 1994, ApJ, 421, 651
  • Artymowicz & Lubow (1996) Artymowicz, P. & Lubow S.H. 1996, ApJ, 467, 77
  • Beck et al. (2012) Beck, T.L. et al. 2012,ApJ, 754, 72
  • Bowler et al. (2011) Bowler, B. P., Liu, M. C., Kraus, A. L., Mann, A.W., & Ireland, M. J. 2011, ApJ, 743, 148
  • Bryden et al. (1999) Bryden, G. et al. 1999, ApJ, 514, 344
  • Casassus et al. (2013) Casassus, S., van der Plas, G. et al. 2013, Nature, 493, 191
  • Comerford et al. (2011) Comerford, J.M. et al. 2011, ApJ, 737, 19
  • Crida et al. (2006) Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587
  • Cuadra et al. (2009) Cuadra, J. et al. 2009, MNRAS, 393, 1423
  • D’Angelo et al. (2006) D’Angelo, G., Lubow, S.H., & Bate, M.R. 2006, ApJ, 652, 1698
  • D’Orazio et al. (2013) D’Orazio, D.J., Haiman, Z., & MacFadyen, A. 2013, MNRAS, 436, 2997
  • Dotti et al. (2009) Dotti, M. et al. 2009, MNRAS, 396, 1640
  • Farris et al. (2014) Farris, B.D. et al. 2014, ApJ, 783, 134
  • Farris et al. (2015) Farris, B.D. et al. 2015, MNRAS, 446, L36
  • Fung et al. (2014) Fung, J., Shi, J.-M.,& Chiang, E. 2014, ApJ, 782, 88
  • Hawley (2000) Hawley, J. F. 2000, ApJ, 528, 462
  • Hawley, Gammie & Balbus (1995) Hawley, J.F., Gammie, C.F. & Balbus, S.A. 1995, ApJ, 440, 742
  • Hawley & Stone (1995) Hawley, J.F. & Stone, J.M. 1995, CoPhC, 89, 127
  • Hawley et al. (2011) Hawley, J.F., Guan, X. & Krolik, J.H. 2011, ApJ, 738, 84
  • Hawley et al. (2013) Hawley, J.F., Richers, S.A., Guan, X., & Krolik, J.H. 2013, ApJ, 772, 102
  • Komossa et al. (2003) Komossa, S. 2003, ApJ, 582, 15
  • Liu & Shapiro (2010) Liu, Y.T. & Shapiro, S.L. 2010, Phys. Rev. D, 82, 3011
  • Lodato et al. (2009) Lodato, G. et al. 2009, ApJ, 398, 1392
  • Lubow (1991) Lubow,S.H. 1991,ApJ, 381, 268
  • Lubow et al. (1999) Lubow, S.H. et al. 1999, ApJ526, 1001
  • Lubow & D’Angelo (2006) Lubow, S.H. & D’Angelo, G. 2006, ApJ, 641, 526
  • Lubow & Artymowicz (2000) Lubow, S.H. & Artymowicz, P. 2000, Protostars and Planets IV, 731
  • MacFadyen & Milosavljević (2008) MacFadyen, A.I. & Milosavljević, M. 2008, ApJ, 672, 83
  • Milosavjević & Phinney (2005) Milosavjević, M. & Phinney, E.S. 2005, ApJ, 622,L93
  • Noble et al. (2010) Noble, S.C., Krolik, J.H., & Hawley, J.F. 2010, ApJ, 711, 959
  • Noble et al. (2012) Noble, S.C. et al. 2012, ApJ, 755, 51
  • Papaloizou & Nelson (2003) Papaloizou, J.C.B. & Nelson, R.P. 2003,MNRAS, 339, 983
  • Pringle (1991) Pringle, J.E. 1991, MNRAS, 248, 754
  • Rafikov (2002) Rafikov, R.R. 2002,ApJ, 569,997
  • Roedig et al. (2012) Roedig, C. et al. 2012, A&A, 545, 127
  • Roedig et al. (2014) Roedig, C., Krolik, J.H., & Miller, M.C., ApJ, 785, 115
  • Rodriguez et al. (2006) Rodriguez, C. et al. 2006, ApJ, 646, 49
  • Rosenfeld et al. (2014) Rosenfeld, K. A., Chiang, E., & Andrews, S. M. 2014, ApJ, 782, 62
  • Rózyczka & Spruit (1993) Rózyczka, M., & Spruit, H.C. 1993, ApJ, 417, 677
  • Rödig et al. (2014) Rödig, C., Krolik, J.H. & Miller, M.C. 2014, ApJ785, 115
  • Savonije et al. (1994) Savonije, G.J., Papaloizou, J.C.B., & Lin, D.N.C. 1994, MNRAS, 268, 13
  • Sawada et al. (1986) Sawada, K., Matsuda, T., & Hachisu, I. 1986,MNRAS,219,75
  • Shakura & Sunyaev (1973) Shakura, N.I. & Sunyaev, R.A. 1973, A&A, 24, 337
  • Shi et al. (2012) Shi,J., Krolik, J.H., Lubow,S.H., & Hawley, J.F. 2012,ApJ, 749, 118
  • Stone & Norman (1992a) Stone, J.M. & Norman, M.L. 1992,ApJS, 80, 753
  • Stone & Norman (1992b) Stone, J.M. & Norman, M.L. 1992, ApJS, 80, 791
  • White & Ghez (2001) White, R.J., & Ghez, A.M. 2001, ApJ, 556, 265
  • Zhou et al. (2014) Zhou, Y., Herczeg, G. J., Kraus, A. L., Metchev, S., & Cruz, K. 2014, ApJ, 783, 17
  • Zhu et al. (2011) Zhu, Z. et al. 2011, ApJ, 729, 47
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description