3D MHD Simulation of Circumbinary Disk

Three Dimensional MHD Simulation of Circumbinary Accretion Disks: Disk Structures and Angular Momentum Transport

Ji-Ming Shi11affiliation: Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218 22affiliation: Center for Integrative Planetary Sciences, Astronomy Department, University of California at Berkeley, Berkeley, CA94720; jmshi@astro.berkeley.edu and Julian H. Krolik11affiliation: Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218 Stephen H. Lubow33affiliation: Space Telescope Science Institute, Baltimore, MD 21218 John F. Hawley44affiliation: Department of Astronomy, University of Virginia, Charlottesville, VA 22903

We present the first three-dimensional magnetohydrodynamic (MHD) simulations of a circumbinary disk surrounding an equal mass binary. The binary maintains a fixed circular orbit of separation . As in previous hydrodynamical simulations, strong torques by the binary can maintain a gap of radius . Streams curve inward from toward the binary; some of their mass passes through the inner boundary, while the remainder swings back out to the disk. However, we also find that near its inner edge the disk develops both a strong asymmetry and growing orbital eccentricity. Because the MHD stresses introduce more matter into the gap, the total torque per unit disk mass is times larger than found previously. The inner boundary accretion rate per unit disk mass is times greater than found from previous hydrodynamical calculations. The implied binary shrinkage rate is determined by a balance between the rate at which the binary gains angular momentum by accretion and loses it by gravitational torque. The large accretion rate brings these two rates nearly into balance, but in net, we find that , and its magnitude is about times larger than predicted by the earlier hydrodynamic simulations. If the binary comprises two massive black holes, the accretion rate may be great enough for one or both to be AGN, with consequences for the physical state of the gas both in the disk body and in its inner gap.

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

1 Introduction

Various types of astronomical binary systems can be embedded in gaseous disks, from young binary stars to stars with growing planets to binary black holes. Such disks have been observed directly in nearby star-forming regions. One of the best resolved (around a young binary) is in GG Tau (e.g., Dutrey et al., 1994; Krist et al., 2005). To date, however, only a few possible disk-planet systems have been directly imaged (e.g., Greaves et al., 2008; Kalas et al., 2008; Hashimoto et al., 2011). Although there is no direct evidence for the existence of circumbinary disks involving binary black hole systems, it is generally believed that such configurations should exist near the centers of galaxies after a galaxy merger (e.g., Begelman et al., 1980; Ivanov et al., 1999; Merritt & Milosavjević, 2005; Escala et al., 2004, 2005; Mayer et al., 2007; Dotti et al., 2009).

Tidal forces exerted by the binary can sometimes clear a gap in the disk. When the binary mass ratio , where denotes the mass of the star and the mass of the planet, the gap that is formed is an annular ring around the primary through which the secondary travels. Whether such a gap opens depends on whether the secondary’s mass is sufficiently large to overcome the gap closing effects of internal stresses. Independent of whether such a gap exists, the secondary can exchange angular momentum with the gas via gravitational torques. Inward orbital migration of the secondary may occur on the disk inflow timescale if the disk mass is large compared to the secondary mass (Lin & Papaloizou, 1986, 1993). Otherwise the migration will be slower (Ivanov et al., 1999; Armitage & Natarajan, 2002). There has been extensive theoretical study of this situation, using both analytic and numerical methods (e.g., Goldreich & Tremaine, 1980; Lin & Papaloizou, 1986, 1993; Bryden et al., 1999; Ivanov et al., 1999; Bate et al., 2003; Nelson & Papaloizou, 2003; Winters et al., 2003).

When is closer to unity, the gap can include the entire binary itself. In this case, the resulting configuration can contain as many as three disks: one around each member of the binary and one that orbits outside the binary, called the circumbinary disk. Observational evidence of such large gap clearing and circumbinary disks has been found in several young stellar binaries (Mathieu, 1994). Numerical simulations have also been applied to study binary systems with a circumbinary disk (e.g., Artymowicz & Lubow, 1994; Bate & Bonnel, 1997; Günther & Kley, 2002; Escala et al., 2005; MacFadyen & Milosavljević, 2008; Hayasaki et al., 2007; Cuadra et al., 2009; Hanawa et al., 2010). These studies have found that the radius of the disk inner edge depends on several factors, including the binary separation, mass ratio, eccentricity, and the strength of the disk turbulence (Artymowicz & Lubow, 1994).

In a one-dimensional model (one that considers only radial dependences), a circumbinary disk would in principle behave as a “decretion” disk’ (Pringle, 1991): the binary loses angular momentum to the surrounding disk so that the disk is repelled. However, the one-dimensional model neglects any non-axisymmetric properties of the disk. For instance, there could be some directions where the binary torque becomes weak enough that matter could leak through the disk inner edge. The gas then penetrates the gap and may accrete onto the binary. Indeed, this penetration effect has been witnessed in numerous two– or three–dimensional simulations of disks with various binary mass ratios using either smoothed particle hydrodynamics (SPH; Artymowicz & Lubow, 1996; Escala et al., 2005; Cuadra et al., 2009; Dotti et al., 2009) or grid-based methods (Bryden et al., 1999; Günther & Kley, 2002; MacFadyen & Milosavljević, 2008; de Val-Borro et al., 2011). Such work has found that in the low-density gap there are gas flows from the disk inner edge to the binary in form of narrow, high-velocity spiral streams. The flow rate through the gap depends on the binary and the disk properties. Compared to the accretion rate near the disk center that would be expected in the absence of binary torque, the accretion rate appears to be reduced (e.g., Lubow et al., 1999; Lubow & D’Angelo, 2006; MacFadyen & Milosavljević, 2008, hereafter MM08). However, Rózyczka & Laughlin (1997) found no reduction at all.

Disk eccentricity is another potential non-axisymmetric property. A disk can become eccentric as a result of its interaction with the binary. Circumbinary disks can, of course, gain eccentricity through direct driving by an eccentric binary (Artymowicz et al., 1991; Papaloizou et al., 2001; Rödig et al., 2011). However, there is also evidence that disk eccentricity can arise even when the orbit of the binary is circular (e.g., Papaloizou et al., 2001, MM08). Simulations of protoplanetary disks have shown that these disks are subject to a resonant mode coupling instability (Kley & Dirksen, 2006; D’Angelo et al., 2006). Through this instability, the disk’s eccentricity can grow although the planet is on a fixed circular orbit. This instability follows the same tidal resonance mechanism found for eccentrically unstable circumstellar disks in superhump binaries (Lubow, 1991; Kley et al., 2008). For circumstellar disks, Lubow (1994) found that orbiting secondaries can drive eccentricity by stream impacts if they are strongly modulated in time. Recently, MM08 found that the disk around an equal mass binary on a circular orbit also became eccentric after a large number of binary orbits. They suggested an eccentricity generation mechanism that involved the action of the binary torque on the gas within the low density gap. Conversely, disk eccentricity might also excite eccentricity in the binary (e.g., Papaloizou et al., 2001). For coalescing massive black holes, the residue eccentricity in the emitted gravitational waves might be detected by proposed instruments like the Laser Interferometer Space Antenna (LISA), perhaps signaling gas-driven evolution (Armitage, & Natarajan, 2005; Key & Cornish, 2011).

Almost all previous studies of circumbinary disks have adopted the -prescription to describe internal stresses and treat either as evolving through diffusion (in one dimension) or as a result of “viscous” stresses (in two or three dimensions). In these efforts, the internal stress to pressure ratio was taken to be constant everywhere and at all times. Although that might be a reasonable approximation for vertically-integrated and time-averaged conditions in the main body of a disk, it becomes unrealistic for highly time-dependent turbulent accretion flows and low density regions outside the disk body(Hawley & Krolik, 2001). Since the exchange of angular momentum between the binary and the disk is crucial for both the circumbinary disk and the binary, we need more a realistic description of the underlying internal stresses. It is now generally recognized that whenever the material of the disk is sufficiently ionized so as to be well-coupled to any embedded magnetic field, the principal mechanism of angular momentum transport is MHD turbulence induced by the magnetorotational instability (MRI). It is therefore necessary to study circumbinary disks using MHD simulations in which internal stress arises self-consistently from turbulence generated by the MRI. To date, this has been done only for extreme mass ratio star-planet systems, using either unstratified (Nelson & Papaloizou, 2003; Winters et al., 2003; Baruteau et al., 2011) or stratified (Uribe et al., 2011) MHD simulations.

It is the goal of this project to construct the first three-dimensional (3D) MHD simulation of a circumbinary disk around an equal mass binary. To simplify our model, we assume the disk and binary to be coplanar. The binary orbits on a fixed circular orbit. On the basis of this simulation, we will try to answer the following questions: (1) What is the inner disk structure? Is the disk truncated or not? (2) How is angular momentum transported within the disk? Can the internal stress balance the binary torque and therefore allow the disk to achieve a quasi-steady state? (3) Is there any eccentricity growth of the disk? If so, what is the cause? (4) How does the accretion rate onto a binary compare with the rate onto single point mass?

We organize this paper as follows: In  2, we describe the physical model and numerical procedures of our circumbinary disk simulations. In  3, we present our simulation results. We then discuss the binary torque and binary contraction in  4.1 and 4.2. In  4.3, we discuss possible mechanisms for disk eccentricity growth. We explain the formation of an asymmetric density concentration near the disk inner edge in  4.4. Finally, in  5 we summarize our conclusions.

2 Numerical Simulation

In this section, we discuss in detail the numerical procedure of this work. The code used is a modern version of the 3D, time-explicit Eulerian finite-differencing ZEUS code for MHD (Stone & Norman, 1992a, b; Hawley & Stone, 1995). We modified the code to cope with the time-dependent binary potential.

2.1 Physical Model

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 time scale is usually much smaller than the shrinkage timescale of the binary, and certainly more than the duration of the simulation (see estimate in  5). We set the gravitational constant , total binary mass and the binary separation to be unity, and therefore the binary frequency 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 . The disk flares at larger radii because the ratio of height to radius , where denotes the scale height of the disk, is the disk radius in cylindrical coordinates, and 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 cut-out 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 . This region is beyond the main extent of the interior disks that surround each binary member because these disks are each tidally truncated at 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 , which is about a factor of outside the cut-out. 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 , where denotes the disk mass. We therefore neglect the contribution of the disk self-gravity to the potential.

2.2 Grid Scheme and Boundary Conditions

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 . The computational domain has to contain at least several () scale heights on each side of the midplane, and for each , several tens of cells are needed. The second is the maximum growth-rate wavelength of the MRI, , where is the Alfven speed, and is the disk rotational frequency. We require to be resolved by at least six grid elements. The last one is the spiral density wavelength, . 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 (,,). We adopt a logarithmic grid in the radial direction, which provides a constant . The vertical grid is derived by mapping a simple linear function for through a polynomial transformation (see equation (6) in Noble et al. (2010)):


where , and are parameters that define the shape of the polynomial. Note that is exactly mapped to the midplane. The merit of this -grid is that for it ensures dense and nearly uniform cell elements close to the disk midplane. The azimuthal grid is evenly spaced and covers all . 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 in , with a computational domain covering radially, meridionally and azimuthally (see also in Table 1), where , , and is the width of the cut-out around the polar axis. The other parameters used in equation (1) are and . Using this grid, we are able to resolve one disk scale height with cells at the inner boundary of the computational domain. The resolution grows to cells per scale height at radius . Within two scale heights of the midplane, is resolved by more than six cells if the plasma is no greater than . The grid also resolves with at least six cells in the plane for .

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 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.

2.3 Initial Setup

We begin with a prograde disk orbiting in the binary plane. Between and , its density is constant in the midplane. The initial disk is axisymmetric, and the polar angle density distribution is , in which 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:


where is the mass ratio of the binary. Here we replace with as for regions near the midplane. We also add to the right hand side of equation (2) to compensate for the small radial gradient of the vertically integrated pressure .

The initial magnetic field is a single poloidal loop within the main body of the disk. It is subthermal with plasma on average. The field is constructed from the vector potential and we define by if and zero otherwise, where , , and is a constant determined by the constraint on the averaged .

2.4 Diagnostics

Three-dimensional MHD simulations usually produce a large amount of data, which poses 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 is defined by


where is the shell surface area. With this definition, the density-weighted shell average is . For example, the net disk accretion rate is


and the average specific angular momentum . The surface density is vertically integrated and azimuthally averaged quantity:


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


and the density weighted vertical average by


Similarly we have the azimuthal average is .

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


where is the Maxwell stress and is the Reynolds stress. The perturbed velocities and are calculated from


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


where and are the average Maxwell stress and Reynolds stress respectively, and is the vertical integral length, which in our case equals the height of the computational domain at a given radius . In a similar way, we can obtain the vertically averaged stress and azimuthally averaged stress as well.

2.5 Hydrodynamic Simulations

To study the mechanism for disk eccentricity growth, we also carried out a set of numerical experiments using two-dimensional 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 throughout the disk because that is roughly the ratio of stress to pressure in the disk body of our MHD simulation (note MM08 used ). The hydrodynamic disks are evolved in polar 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 which 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 for . Its grid covers the same physical extent in radial and azimuthal directions as the MHD one. The initial surface density of the two-dimensional 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 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 and fixed. The properties of both the reruns and the control run can be found in Table 1.

3 Results

We present one 3D MHD simulation, called B3D, in this paper.111We also performed a short duration rerun for with higher dump rates: ten history dumps and one 3D dump per time unit. They are used when high time resolution is required, e.g. when we try to investigate the angular momentum budget and the stream dynamics. This simulation was terminated after the inner disk () reached a quasi-steady stage for several hundred time units. Longer simulations require better resolution to resolve the MRI in the growing density concentration region (see  3.3.3).

B3D ran for code units, which corresponds to binary orbits. The gas quickly fills in the initially empty region within several binary orbital periods, and after units of time, the disk becomes fully turbulent. The binary torque then is able to maintain a low density gap, and the inner part of the disk finally reaches a quasi-steady state after . We note that this single simulation consumed K CPU hours on the Kraken Cray XT5 system.

In the following subsections, we first describe the overall evolution of the circumbinary disk and its quasi-steady state. We then try to discuss how the angular momentum is transported in the inner disk in 3.2. We will mainly discuss the characteristic disk structures in 3.3 and 3.4. The field structure will be considered in  3.5 and the temporal properties of the accretion in the last subsection.

3.1 A Secularly Evolving Quasi-steady State

After the first time units, the binary torque starts to clear out a low density gap between the inner boundary and . In the top-left panel of Figure 1, we show the vertically integrated surface density of the circumbinary disk at in the plane. In that panel, we can clearly see 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 § 3.4 and 4. At the disk appears to be elliptical and slightly off center, and at and the disk obviously becomes eccentric. We also observe a growing azimuthal density asymmetry. We will discuss the eccentricity growth and density asymmetry in § 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.

We show the enclosed mass (integrated disk mass interior to a given radius) inside , , , and as a function of time in Figure 2. We find that the mass at (i.e., inside the initial disk’s inner edge) undergoes dramatic growth during the first time units as the initial disk fills in regions with . 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 . We therefore define the quasi-steady stage as the period between and the end of the simulation. On the other hand, the mass inside and rises by a factor of two after , indicating that this region does not achieve inflow equilibrium. The mass-interior profiles also possess high frequency fluctuations (typical time scale time units) due to the orbital modulation of mass contained in the streams. In addition, after time units, there are slower oscillations (time scale time units) of the enclosed mass caused by the eccentrically orbiting density lump.

In Figure 3, we display the time-averaged surface density and accretion rate over two intervals: () and (). We find there is only minor change in the averaged surface density and accretion rate of the disk region over and . If we define the unit of surface density , the peak density grows in time, but only slowly, rising from to . Its position also changes, but similarly slowly, gradually moving from to . The shape of the surface density profile likewise hardly changes in time. On the inside, the disk is truncated exponentially, for both intervals, where represents an average over time. Toward larger radii, the surface density is , the predicted profile of a ‘decretion disk’ (Pringle, 1991). All that changes is that the region where extends farther outward at later times.

Averaged over the quasi-steady period, the accretion rate at the inner boundary is . However, the disk evolves throughout the simulation at . We see gas inflow rates as large as outside the region during both time intervals. As a result, the surface density continuously build 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’ , where and . 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


The surface density of such disks is very uncertain. A sense of scale, however, can be gleaned from translating eqn. 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 , where we have scaled to a total black hole mass of , a binary separation of 0.1 pc, and a disk column density whose Thomson optical depth is .

3.2 Angular Momentum Budget

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 binary torques and the 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 , the character of these stresses is very similar in our two averaging periods, and . The Reynolds stress shows big fluctuations radially. Its first peak, at , is . It dips at around to nearly zero and then rises up again slowly. The ratio of the first peak to the second peak is . However, the Maxwell stress has a flatter radial profile. It is roughly on average and varies only slowly over a range of at most a factor of two. It diminishes within the gap and decreases smoothly towards larger distance. Comparing the two, we find the Reynolds stress exceeds the Maxwell stress by a factor at . On the other hand, beyond the gap region, the MHD stress always dominates the Reynolds stress by a factor . 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 , 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.

In the language of the parameter of Shakura & Sunyaev (1973), we can also measure the stresses in units of the pressure. In the right-hand 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 at . 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 Reynolds stress (top left) and Maxwell stress (top right) taken at time units are shown in Figure 5. On top of each plot, we also draw the surface density with contour lines, logarithmically spaced from the density floor value to the density maximum. Nearly all the Reynolds stress is confined within . The most interesting point is that we find relatively large stress in a single gas stream extending from (see the red stripe in the top left panel), where the gas gains angular momentum () and is being pushed out by the binary. The gas being kicked thus goes on an eccentric orbit, creating negative stress at and ( and ). In the same range of radii but a different azimuthal location (), the gas falls back toward the binary, creating positive stress (both and ). 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 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.

To illustrate how the stresses depend on the disk height, we present azimuthally-averaged snapshots at the same time in Figure 5. We find 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:


is the local torque, and the integrated torque is


where is the binary potential in the midplane.

The first panel of Figure 6 shows the averages of the local torques for the two intervals and 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 and peaks at with as the peak value. It then goes negative at , reaching its first negative maximum at , but with a smaller magnitude . Similar to what was previously found in MM08 and Cuadra et al. (2009), we find the torque density oscillates around zero but damps quickly toward larger distance. In addition, the torque density for 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 . 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 . In our MHD simulation, the first positive peak of the torque takes place at while it is in MM08. In terms of this normalization, in the hydrodynamic simulation of MM08. In our MHD simulation, we find that it is (averaged in the quasi-steady period), a factor of 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 in our MHD simulation. However it is only in MM08, 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,


where is the specific angular momentum of the binary.

Alternatively, as is roughly around for both MHD and hydrodynamic simulations, we can write the binary torque as some factor times , and then we have in our present simulation, while MM08 gives , times smaller. The largest part of this contrast can be explained by the fact that in the MHD treatment , a factor of 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 internal stress near the disk inner edge. Therefore, the binary torque should increase linearly with the effective . In our case, in the disk body at , a factor of greater than the constant 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. By dividing by the surface density, we remove the effects due to redistribution of the surface density. We find the specific torque slightly shifts to larger radius and diminishes its amplitude from to . In the bottom right panel, we plot the history of the total torque smoothed over a short time span , where is the binary period. Before smoothing, the total torque fluctuates very rapidly () between and in the quasi-steady period. The power spectrum of this fluctuation is dominated by a peak at , approximately the beat frequency between the orbital frequency of the matter being torqued most strongly (at the peak of , i.e. ) and . After smoothing, we find the binary torque has a sudden rise around which is due to the initial disk filling the small radius region. The torque then levels out between . Once in the quasi-steady stage, the torque shows a slowly decreasing trend, in fractional terms over the final time units. It also appears to oscillate at the disk orbital frequency () 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 simply reflects the increasing and .

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 either be transported outward by Maxwell and Reynolds stresses or advected towards 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 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 (, ) as


where is the shell-integrated angular momentum, is the mean angular momentum flux due to mass motion, and 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 –legend, we use to represent all the terms in the equation, which are all radial derivatives of one or another variety of shell-integrated angular momentum flux ().

During this time interval, the local angular momentum stays steady, with only very small temporal variations from the averaged (solid black curve), as expected for a quasi-steady accretion disk. The binary torque (blue dashed curve) at is almost balanced by the differential flux due to the sum of Reynolds stress and Maxwell stress (cyan dashed curve). In other words, the circumbinary disk within that region reaches a balance between binary torque and internal stresses. As a result, the angular momentum flux carried by accretion flows oscillates about zero (magenta curves within ). Clearly, the internal torque due to the Reynolds stress (green dash-dotted curve) dominates the Maxwell stress (red dash-dotted curve) at , which suggests that most of the angular momentum dumped by the binary at is transported outward via the Reynolds stress, while the angular momentum drawn from the binary at is compensated by the torque of the Reynolds stress. Finally, we comment that torque balance is not reached at (where negative binary torque dominates) or (where 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. Secondly, the distribution of both stresses follows the gas streams and are vertically confined: within two scale heights for Reynolds stress and four scale heights for Maxwell stress. Third, we find 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 (). 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 Disk Eccentricity

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


where is the shell average defined as in equation (3). We show this local eccentricity in the space-time diagram in Figure 8 (top panel). The local eccentricity of the gas in the gap () rises significantly during the period . The increasing in the gap simply reflects the fact that the dynamics there are strongly non-Keplerian. The eccentricity in the main body of the disk () grows rapidly over the period . The typical value increases from much less than throughout the disk body at to between at . It keeps growing slowly after time units, and at the end of the simulation a ring of eccentric disk forms at with typical . 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 to . 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 which connect the gap with the disk body. These features bend at pointing to later times at both ends, which suggests 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  4.3.3.

In Figure 9, we present the radial distribution of the disk eccentricity averaged over three time intervals: (solid black), (red dotted) and (blue dashed). The distribution curves in the disk body () shift toward larger radius and greater eccentricity as the disk evolves in time. Outside the peaks around , the slopes of the distribution curves are constant. The radial dependence of in the disk body can be roughly described as .

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 2-d hydrodynamical 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 (). 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 both by 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 4.3.1.

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


At early times (), 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 to . By fitting the eccentricity curve during that interval, we find the growth rate is , where the subscript denotes eccentricity. The growth rate of disk eccentricity slows by a factor of after . This suggests some nonlinear mechanisms come into play to limit and perhaps eventually saturate eccentricity growth.

3.3.2 Precession of the Eccentric Disk

To follow changes in orientation of the eccentricity, we define the longitude of the apoapse in terms of the complex phase angle of the disk eccentricity for . Its time evolution is shown in Figure 11. Here the angle measures the angular location of the apoapse with respect to the –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 , we find . Linear perturbation theory predicts the precession rate for particles around the binary potential is for a cold disk, where is the angular frequency and is the epicyclic frequency of the disk. To a first order approximation, (also see equation (2)), and . Thus we find the prediction is for , or evaluated at . 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.

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 mode strength component


which is basically the component of the Fourier transform of the surface density. We plot the space-time diagram of in Figure 12. The orbiting lump appears as a zig-zag pattern whose peak moves in and out between and after about time units. The pattern period, which is also the orbital period of the eccentric disk, is time units. We will discuss the lump in  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.

3.4 Streams in the Gap

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 not only affect the accretion rate of the circumbinary disk, but also the structure of the inner part of the disk.

In a frame that co-rotates with the binary, we find the streams follow roughly the static bi-symmetric potential of the binary and pass through the inner boundary about radians beyond the Lagrange and points (see Figure 13). This phase offset is consistent with the negative average torque at in § 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 and in the first two panels of Figure 13, and draw the velocity field vectors measured in the co-moving frame on top of it. The snapshots at show the inflowing gas is regulated by the binary torque as it leaves the inner disk edge and streams towards 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 co-rotating frame of the binary, the gas in the second (near ) and fourth (near ) quadrants is pulled toward the closer component of the binary (see 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 , the velocity is mostly outward at 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 and . In the inertial frame, the binary almost always rotates faster than the disk; therefore, the second and fourth quadrants have positive binary torque. As a result, gas in these regions gains angular momentum from the binary and then drifts away from the center.

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, for instance at , 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 ,, and in the corotating frame. The stream on the right at is gradually pushed away at and 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 .

We also measure the fluid and magnetic effects on the dynamics of the streams. In the inner disk and the gap, the force ratio is always , where is the force density of the gas pressure gradient and 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 § 4.3.3.

3.5 Field Structure

We found in the previous subsection that the Maxwell stress closely follows the gas streams in the plane and is well confined near to the midplane. Now we show 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 . 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.

3.6 Time Dependent Accretion

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 . 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 time units). These two time variations are the most significant modes in the temporal profile of . 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 , which is the same as the orbital frequency of the lump during the later stages of the simulation (as shown in § 3.3). Our result in Figure 15 is largely consistent with previous findings: a dominant frequency associated with the binary orbital frequency ( if ; MM08; Hayasaki et al., 2008; Cuadra et al., 2009; Rödig et al., 2011; Sesana et al., 2011), 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 ) 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 rather than because the members of the binary have exactly equal mass.

To further verify the cause of these two modes, we pick out two times ( and as marked in the top panel by black arrows) and plot their two-dimensional 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 to . These lines pick out the location of the streams, while the white contour lines show the surface density between to in a linear scale; these contours delineate the lump region. We find the outbursts of 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 .

The time averaged accretion rate at the inner edge, normalized with the peak density, is , a factor greater than found in MM08. Interestingly, during the quasi-steady period, the accretion rate in the gap is nearly twice the rate measured at , 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 cut-outs suggest that the rate hardly changes as long as the cut-out size is smaller than .

4 Discussion

4.1 Binary Torque and Linear Resonance Theory

Using the results we have on the disk stresses and the binary torque, we want to further elucidate the angular momentum transport mechanisms within a 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 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 satelites (Goldreich & Tremaine, 1982). To illustrate that contrast, we calculate both the torque density and total torque as functions of radius, assuming most of the torque comes from the outer Lindblad resonance at (e.g., Meyer-Vernet & Sicardy, 1987). The estimate is shown in the top panel of Figure 6, with the red lines representing 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 , which is comparable to the first wavelength , where is evaluated at the resonance radius, and is used here (see Figure 1 and the pressure-dominated case in Table  of Meyer-Vernet & Sicardy (1987)). Because the surface density at is so much smaller than at , one might expect the maximum torque to move away from the resonance point toward the inner disk. Secondly, linear theory assumes nearly circular orbits because it also assumes . However, in our simulation , 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).

4.2 Binary Orbital Evolution

The time-averaged rate at which the binary loses angular momentum to the disk is . Meanwhile, the time-averaged angular momentum flux across the inner boundary due to accretion is , where is the averaged specific angular momentum at . 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


where denotes the angular momentum of the binary, and its rate of change is .

After separating the accretion term of and combining it with the second term in equation (19), we have


If (here it is ), the shrinkage rate is controlled by the competition between the first and second terms. Plugging in our numbers, we find the two terms are surprisingly comparable in magnitude: and respectively. The net result is orbital compression, but at a rate much smaller than either of the two terms, .

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


where denotes the binary’s orbital eccentricity. In addition, the angular momentum of the binary should also be adjusted to account for the eccentricity: . 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 . 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  3.2.2, we define a characteristic disk mass , and further assume that the binary orbit remains circular, then . For comparison, MM08 found a rate , which can be scaled used the definition of to , 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 , a factor of 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 , , , as well as the detailed mechanics of how the stream joins the interior disks associated with the members of the binary.

4.3 Disk Eccentricity Growth

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 time units (see Figure 10). Surface density is taken directly from the simulation by averaging over this time interval. The disk inner radius is taken to be and the outer radius . We assume the same equation of state as in the simulation: isothermal, with sound speed . The adiabatic index is taken to be . The small imaginary part (which introduces slow damping) assures mathematical consistency. We model the forces driving the initiation of eccentricty as a Gaussian in radius (see equation (B11)), centered on and having a width . These parameters are not well constrained by the simulation. It follows that the the results are also insensitive to their values within certain limits. Parameter , 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, , where is the eigenvalue defined in equation (B14). The value is found to be .

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, . The confinement of eccentricity is due to the effects of differential precession of the binary (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 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 (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 .

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 , which indicates that the source of the eccentricity should be located near the inner edge of our MHD disk. For larger values of , there is no value of 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 . If the width is doubled to , 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 (see Figure 9 in  3.3.1). In Appendix C, we derive an asymptotic form for the eccentricity distribution (valid for ) using linear theory. It is dominated by an exponential falloff in with an e-folding length also of about 0.8a.

To sum up the analysis in this subsection, we find the linear model can explain the eccentricity distribution of our simulation. It shows 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 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


where 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, is excluded because that component of the tidal potential is absent. For , the instability is associated with the resonance that occurs at . In this model, the tidal field interacts with the eccentric motions of the disk to produce disturbances of the form . These disturbances launch waves of that form at the resonance. The waves in turn interact with the tidal field to produce a stress that increases the disk eccentricity exponentially in time.

However, the resonance lies at , 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 , where is the width of the eccentricity distribution and is . This very rapid growth rate is a consequence of the powerful 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 it is at most a limited role.

4.3.3 Eccentricity Growth through Stream Impact

As discussed in 3.3.1, the space-time diagram reveals a special feature of the growth: the propagation of eccentricity from radii close to 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 cut-off. 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 (§ 4.3.2). Because these are essentially hydrodynamic effects and involve no vertical dynamics, these special simulations were performed in 2-d with no magnetic field (see  2.5 and Table 1 for details of the simulations).

In the control run (B2D.rin=0.8), the disk is truncated around , similar to B3D. In the gap region, there is a pair of streams at while the disk eccentricity grows exponentially. A single strong stream appears when the disk eccentricity becomes significant. The simulation is terminated at , and at that time the disk eccentricity is saturated. We find 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 two dimensional definition of equation (17). Similar to the MHD simulation, we find exponential growth from . The growth rate in that period is , very close to that of the MHD simulation (). We then restart from with five different values of while keeping all other parameters fixed. We terminate the reruns at , the stopping time for B2D.ri=0.8.

In Figure 17, we plot the evolution of disk eccentricities for cases with different . The case with (blue dash-dotted) 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 cut-offs (: green dashed curve; : red dotted curve), however, cease 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 (cyan dash dot dot curve) and (magenta long dash 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 , 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 whose presence is required for the mechanism of § 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 , when the eccentricity growth is slow, the inner portion of the disk forms an eccentric ring between with . Its pericenter is at and , and the apocenter is at and in 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


where is the semi-major axis of the eccentric ring of disk and is the mean motion of the disk. Quantity is the true anomaly where the stream impact occurs. and are the radial and angular components of the disturbance force density, and are defined as


where denotes the mass injection rate of the stream, is the total mass of the eccentric ring, and 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 () in order to estimate the disk perturbation. We carry out our analysis using the vertically-integrated two-dimensional data. For the Gauss equation to apply, the properties of the eccentric ring, such as , , , 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: , and  radians. The stream-disk impact is localized in space. We analyze the stream-disk impact over an arc that is defined by and . The quantities and 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 is obtained from the unaffected matter at a slightly greater radius, which we take to be an arc at . 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 that range between and at each half binary period. Most of the time, 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 . The corresponding growth rate is consistent with the growth rate at earlier times in the simulation (). 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.

4.4 Interpretation of the Lump

The disk after shows a significant asymmetry which we call the “density lump”. We believe the lump is due to the combined effect of stream impact and disk eccentricity. As discussed in § 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 an episode of stream reabsorption taking place.

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 , the density enhancement orbits eccentrically, and as it follows its orbit toward apocenter, both and 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 () shortly before an outgoing stream reaches the lump. Three time units later, the small enhancement (the green region close to at ) joins the main lump (green and red colored area in the second quadrant). The density-weighted velocity in the enhancement is (we average the stream gas where the surface density is greater than ). However, the lump has a slower velocity (averaging over locations where ). As a result of both the lump feeding and stream reabsorption, the density lump grows both in mass and density contrast.

Another way to illustrate the connection between streams and the lump is to look at the space-time diagram in Figure 12. In the diagram, the stream clearly stands out as the radially stretched features between and which propagate away from the binary with a pattern speed of about . This pattern repeats itself at a rate about twice the binary frequency. Meanwhile, the zigzag oscillation at late time () between radius and shows the movement of the density lump as it orbits around the binary on an eccentric orbit with a period of 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 in the ascending part of the zig-zag. 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 further away from the stream in azimuthal angle, and thus we find the descending part at is less affected by the stream features. As the lump approaches pericenter again, another cycle of stream reabsoprtion and lump feeding begins. Note that the blue dots (low mode strength regions) around between the zig-zags 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 mode at that radius. The zig-zag 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.

5 Conclusions

5.1 Specific Results

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

  1. The disk exhibits a number of nearly steady features. There is a low-density gap within , an eccentric inner disk at . At early times, there is a pair of gas streams that flow into the gap from 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 the lump is due to a combination of stream impact and disk eccentricity.

  3. The disk eccentricity grows exponentially during the time (16–40 binary orbits) with a growth rate of . The growth rate then slows down significantly. By the end of the simulation, the disk body at has reached an eccentricity . Its pericenter precesses slowly due to the quadrupolar component of the binary potential. Stream impact largely accounts for the eccentricity growth.

  4. Reynolds stress associated with the streams is the leading transport mechanism in the gap region, while 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 times greater than 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 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 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 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 focussed 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 is only slightly larger than the rate estimated by MM08.

5.2 Consequences for Orbital Evolution: vs.

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  pc (Begelman et al., 1980). Recently, there have been several attempts to solve this ‘final–parsec problem’ with stellar dynamics. Although they 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., 2011; Farris et al., 2011).

We have previously described this process in terms of the orbital shrinkage rate . Defining the orbital shrinkage time , the time taken by torque alone to change the binary orbital angular momentum , and the disk accretion time , we can rewrite this rate as


Using the values found in the simulation makes , 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 , it is not necessarily true that this mass must be there all at once.

Suppose, for example, that . Because 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’s 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 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 ; conversely, if it were increased relative to a constant torque by a ratio , 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 () 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 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 AGN 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 ) would alter saturation of the magneto-rotational instability, but more rapid inflow is one plausible consequence.

The opposite case is also worth considering, in which , so that . 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 , where represents the time- and density-weighted shell-averaged inflow velocity. In a steady-state disk around a point-mass, would match at . However, to the degree that binary torques retard accretion, exceeds the inflow time measured at the surface density peak, which is . At sufficiently large radius,