Three Dimensional MHD Simulation of Circumbinary Accretion Disks: Disk Structures and Angular Momentum Transport
Ji-Ming Shi, Julian H. Krolik, Stephen H. Lubow, John F. Hawley
TThree Dimensional MHD Simulation of Circumbinary AccretionDisks: Disk Structures and Angular Momentum Transport
Ji-Ming Shi , and Julian H. Krolik andStephen H. Lubow andJohn F. Hawley ABSTRACT
We present the first three-dimensional magnetohydrodynamic (MHD) simu-lations of a circumbinary disk surrounding an equal mass binary. The binarymaintains a fixed circular orbit of separation a . As in previous hydrodynamicalsimulations, strong torques by the binary can maintain a gap of radius (cid:39) a .Streams curve inward from r (cid:39) a toward the binary; some of their mass passesthrough 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 m = 1 asymmetry and growing orbital eccentricity. Because the MHD stressesintroduce more matter into the gap, the total torque per unit disk mass is (cid:39) (cid:39)
40 times greater than found from previous hydrodynamical calcu-lations. The implied binary shrinkage rate is determined by a balance betweenthe rate at which the binary gains angular momentum by accretion and loses itby gravitational torque. The large accretion rate brings these two rates nearlyinto balance, but in net, we find that ˙ a/a <
0, and its magnitude is about 2 . Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218 Center for Integrative Planetary Sciences, Astronomy Department, University of California at Berkeley,Berkeley, CA94720; [email protected] Space Telescope Science Institute, Baltimore, MD 21218 Department of Astronomy, University of Virginia, Charlottesville, VA 22903 a r X i v : . [ a s t r o - ph . H E ] F e b Subject headings: accretion, accretion disks — binaries: general — MHD —methods: numerical
1. INTRODUCTION
Various types of astronomical binary systems can be embedded in gaseous disks, fromyoung binary stars to stars with growing planets to binary black holes. Such disks have beenobserved directly in nearby star-forming regions. One of the best resolved (around a youngbinary) is in GG Tau (e.g., Dutrey et al. 1994; Krist et al. 2005). To date, however, only afew possible disk-planet systems have been directly imaged (e.g., Greaves et al. 2008; Kalaset al. 2008; Hashimoto et al. 2011). Although there is no direct evidence for the existenceof circumbinary disks involving binary black hole systems, it is generally believed that suchconfigurations should exist near the centers of galaxies after a galaxy merger (e.g., Begelmanet al. 1980; Ivanov et al. 1999; Merritt & Milosavjevi´c 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 thebinary mass ratio q ≡ M /M (cid:28)
1, where M denotes the mass of the star and M themass of the planet, the gap that is formed is an annular ring around the primary throughwhich the secondary travels. Whether such a gap opens depends on whether the secondary’smass is sufficiently large to overcome the gap closing effects of internal stresses. Independentof whether such a gap exists, the secondary can exchange angular momentum with the gasvia gravitational torques. Inward orbital migration of the secondary may occur on the diskinflow timescale if the disk mass is large compared to the secondary mass (Lin & Papaloizou1986, 1993). Otherwise the migration will be slower (Ivanov et al. 1999; Armitage &Natarajan 2002). There has been extensive theoretical study of this situation, using bothanalytic 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 q is closer to unity, the gap can include the entire binary itself. In this case, theresulting configuration can contain as many as three disks: one around each member of thebinary and one that orbits outside the binary, called the circumbinary disk. Observationalevidence of such large gap clearing and circumbinary disks has been found in several youngstellar binaries (Mathieu 1994). Numerical simulations have also been applied to study 3 – q (cid:46) α -prescription todescribe 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 internalstress to pressure ratio α was taken to be constant everywhere and at all times. Although thatmight be a reasonable approximation for vertically-integrated and time-averaged conditionsin the main body of a disk, it becomes unrealistic for highly time-dependent turbulentaccretion 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 forboth the circumbinary disk and the binary, we need more a realistic description of theunderlying internal stresses. It is now generally recognized that whenever the material ofthe 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 bythe magnetorotational instability (MRI). It is therefore necessary to study circumbinarydisks using MHD simulations in which internal stress arises self-consistently from turbulencegenerated by the MRI. To date, this has been done only for extreme mass ratio star-planetsystems, using either unstratified (Nelson & Papaloizou 2003; Winters et al. 2003; Baruteauet 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 simula-tion of a circumbinary disk around an equal mass binary. To simplify our model, we assumethe disk and binary to be coplanar. The binary orbits on a fixed circular orbit. On the basisof this simulation, we will try to answer the following questions: (1) What is the inner diskstructure? Is the disk truncated or not? (2) How is angular momentum transported withinthe disk? Can the internal stress balance the binary torque and therefore allow the disk toachieve a quasi-steady state? (3) Is there any eccentricity growth of the disk? If so, what isthe cause? (4) How does the accretion rate onto a binary compare with the rate onto singlepoint mass?We organize this paper as follows: In §
2, we describe the physical model and numericalprocedures of our circumbinary disk simulations. In §
3, we present our simulation results.We then discuss the binary torque and binary contraction in § § § §
2. NUMERICAL SIMULATION
In this section, we discuss in detail the numerical procedure of this work. The code usedis 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 thetime-dependent binary potential.
We construct our physical disk model in the inertial frame in which the center of massof the binary is at rest at the origin. Treating the simplest case first, we assume an equalmass circular restricted binary system, i.e. the binary orbits circularly in the disk plane, andwe neglect binary evolution as the disk inflow time scale is usually much smaller than theshrinkage timescale of the binary, and certainly more than the duration of the simulation(see estimate in § G , total binary mass M and thebinary separation a to be unity, and therefore the binary frequency Ω bin = (cid:112) GM/a isunity as well. We also assume the circumbinary disk to be cold and thin. As we aremainly concerned with the orbital dynamics of the flow, we choose a simple global isothermalsound speed c s = 0 .
05. The disk flares at larger radii because the ratio of height to radius
H/R = c s /R Ω K = 0 . R/a ) / , where H denotes the scale height of the disk, R is thedisk radius in cylindrical coordinates, and Ω K = (cid:112) GM/R is the Keplerian frequency.As the fluid is well coupled to the magnetic field even in a cold disk where the ionizationfraction is far below unity(see, e.g., Stone et al. 2000), we assume ideal MHD. For thisreason, we include no explicit diffusivity except the von Neumann-Richtmyer bulk viscosityin compressive regions that ensures the right jump conditions for shock waves. The effecton 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 thebinary, we excise a central region. This cut-out must be well inside the inner edge of thedisk so that it does not affect fully resolving the inner disk and any gas leakage from theinner edge. We choose to cut out the area within 0 . a . This region is beyond the mainextent of the interior disks that surround each binary member because these disks are eachtidally truncated at ∼ . a from their central objects (Paczy´nski 1977). Once the diskreaches the quasi-steady stage, we find that the inner edge of the disk is located at r (cid:39) a ,which is about a factor of 2 . Q ≈ ( H/R )( M/M d ) (cid:29)
1, where M d denotes the disk mass. We therefore neglectthe contribution of the disk self-gravity to the potential. The properties of circumbinary disks require us to resolve three different length scaleswhen setting up the grid. The first is the half disk thickness H . The computational domainhas to contain at least several ( (cid:46)
4) scale heights on each side of the midplane, and for each H , several tens of cells are needed. The second is the maximum growth-rate wavelengthof the MRI, λ MRI ≡ π/ √ v A / Ω( R ), where v A is the Alfven speed, and Ω( R ) is the diskrotational frequency. We require λ MRI to be resolved by at least six grid elements. The lastone is the spiral density wavelength, λ d ∼ πc s / Ω bin . There should be many cells across thiswavelength (MM08); we require at least six.In order to satisfy the above requirements, we follow the scheme proposed in Noble etal. (2010) to construct our grid in spherical coordinates ( r , θ , φ ). We adopt a logarithmicgrid in the radial direction, which provides a constant ∆ r/r . The vertical grid is derived bymapping a simple linear function y ( x ) = x for x ∈ [0 ,
1] through a polynomial transformation(see equation (6) in Noble et al. (2010)): θ ( y ) = π (cid:20) − ξ )(2 y −
1) + (cid:18) ξ − θ c π (cid:19) (2 y − n (cid:21) , (1)where ξ , θ c and n are parameters that define the shape of the polynomial. Note that y = 0 . θ -grid is that for n > π . Theoretically speaking, the aspect ratios of the cell shapeshould be as isotropic as possible, but the azimuthal cell widths can be a factor of a fewlonger than the other two, as the shearing tends to draw out features in this direction.The grid resolution used in the present simulation is 400 × ×
540 in ( r, θ, φ ), witha computational domain covering [ r in , r out ] radially, [ θ c , π − θ c ] meridionally and [0 , π ] az-imuthally (see also in Table 1), where r in = 0 . a , r out = 16 a , and θ c = 0 . ξ = 0 . n = 9. Using this grid, we are able to resolve one disk scale height with 20 cells at theinner boundary of the computational domain. The resolution grows to ∼
40 cells per scale 7 –height at radius ∼ a . Within two scale heights of the midplane, λ MRI is resolved by morethan six cells if the plasma β is no greater than ∼ λ d with atleast six cells in the r – φ plane for r (cid:46) a .We choose outflow boundary conditions in the radial and meridional directions for thegas, which permit only flows going outward; any inward velocities are set to zero. As wecover all 2 π of φ , the boundary conditions for φ –grid are simply periodic. For the magneticfield in the radial and meridional directions, we set the transverse components of the fieldto be zero in the ghost zones. The components normal to the boundaries are calculated byimposing the divergence-free constraint. We begin with a prograde disk orbiting in the binary plane. Between r min = 3 a and r max = 6 a , its density is constant in the midplane. The initial disk is axisymmetric, andthe polar angle density distribution is ρ = ρ exp [ − ( θ − π/ / ( √ H/r ) ], in which ρ = 1is the unit of disk density. This form provides initial hydrostatic balance vertically for apoint mass potential and zero radial pressure gradient along the midplane. For a first orderapproximation, the difference between a binary potential and a point-mass in the midplanecan be described by the temporally and azimuthally averaged quadrupole moment of thebinary potential. We therefore modify the angular frequency of the initial disk to accountfor the quadrupole contribution:Ω( r ) ≈ Ω (cid:20) (cid:16) ar (cid:17) q (1 + q ) (cid:21) , (2)where q = 1 is the mass ratio of the binary. Here we replace R with r as r = R sin θ ≈ R forregions near the midplane. We also add 1 / ( r Σ) dP/dr to the right hand side of equation (2)to compensate for the small radial gradient of the vertically integrated pressure P = Σ c s .The initial magnetic field is a single poloidal loop within the main body of the disk.It is subthermal with plasma β = 100 on average. The field is constructed from the vectorpotential A = (0 , , A φ ) and we define A φ by A φ = A √ ρ sin(2 πr/kH )( r/r min − − r/r max ) − √ ρ cut if A φ > k = 2Ω bin a/c s , ρ cut = 10 − ρ , and A isa constant determined by the constraint on the averaged β . 8 – Three-dimensional MHD simulations usually produce a large amount of data, whichposes challenges for data storage, transport and post-simulation analysis. In order to fa-cilitate the study of the spatial and temporal properties of the disk, we write out spatiallyaveraged history data at a frequency of one dump per time unit and write out 3D snapshotsevery five units of time.We use two different types of spatially averaged history data. The first is defined byeither an integration or an average over radial shells. Shell-averaging for variable X is definedby (cid:104) X (cid:105) ( r, t ) ≡ (cid:82) Xr sin θdθdφ (cid:82) r sin θdθdφ , (3)where (cid:82) r sin θdθdφ ≡ A x ( r ) is the shell surface area. With this definition, the density-weighted shell average is (cid:104) X (cid:105) ρ ≡ (cid:104) ρX (cid:105) / (cid:104) ρ (cid:105) . For example, the net disk accretion rate is˙ M ( r, t ) ≡ (cid:90) ρv r r sin θdθdφ = A x (cid:104) ρv r (cid:105) , (4)and the average specific angular momentum l = (cid:104) v φ r (cid:105) ρ = (cid:104) ρv φ r (cid:105) / (cid:104) ρ (cid:105) . The surface density isvertically integrated and azimuthally averaged quantity:Σ( r, t ) ≡ π (cid:90) ρr sin θdθdφ. (5)The second type of average is two-dimensional, either an azimuthal average of poloidalslices or a vertical average referred to the equatorial plane. We define the vertical averageby (cid:104) Y (cid:105) z ( r, φ, t ) ≡ (cid:82) Y r sin θdθ (cid:82) r sin θdθ , (6)and the density weighted vertical average by (cid:104) Y (cid:105) z,ρ ( r, φ, t ) ≡ (cid:82) ρY r sin θdθ (cid:82) ρr sin θdθ = (cid:104) ρY (cid:105) z (cid:104) ρ (cid:105) z . (7)Similarly we have the azimuthal average is (cid:104) Y (cid:105) φ ( r, θ, t ) ≡ π (cid:82) Y dφ .We need to be very careful about the definition of the r – φ component of the internalstress in the present simulation. We follow Hawley & Krolik (2001) and define the stress as w rφ ( r, θ, φ, t ) = t rφ + r rφ = − B r B φ π + ρδv r δv φ , (8) 9 –where t rφ is the Maxwell stress and r rφ is the Reynolds stress. The perturbed velocities δv r and δv φ are calculated from δv r ( r, θ, φ, t ) = v r − (cid:104) v r (cid:105) ρ δv φ ( r, θ, φ, t ) = v φ − (cid:104) v φ (cid:105) ρ . (9)The vertically integrated and azimuthally averaged total stress can then be described by W rφ ≡ π (cid:90) w rφ r sin θdθdφ = L z (cid:104) w rφ (cid:105) = T rφ + R rφ ,T rφ = L z (cid:28) − B r B φ π (cid:29) ,R rφ = L z (cid:20) (cid:104) ρv r v φ (cid:105) − (cid:104) ρv r (cid:105)(cid:104) ρv φ (cid:105)(cid:104) ρ (cid:105) (cid:21) , (10)where T rφ and R rφ are the average Maxwell stress and Reynolds stress respectively, and L z ≡ A x / (2 πr ) is the vertical integral length, which in our case equals the height of thecomputational domain at a given radius r . In a similar way, we can obtain the verticallyaveraged stress (cid:104) w rφ (cid:105) z and azimuthally averaged stress (cid:104) w rφ (cid:105) φ as well. To study the mechanism for disk eccentricity growth, we also carried out a set of numer-ical experiments using two-dimensional viscous hydrodynamic simulations. Their purposewas both to distinguish hydrodynamic from MHD effects and to test the effect of where theinner boundary is placed.The hydrodynamic simulations in this paper used the α -disk prescription. We chose α = 0 . α = 0 . r, φ ) coordinates in the inertial frame. Following the grid scheme of theMHD simulation, we set the radial grid to be logarithmically spaced, while the azimuthalgrid is spaced uniformly.Among these simulations, B2D.rin=0.8 serves as the control run. The parameters whichdescribe the physical properties of the hydrodynamic disk and the binary are kept the sameas in the MHD case. The resolution of the control run is 512 × r × φ . Its gridcovers the same physical extent in radial and azimuthal directions as the MHD one. Theinitial surface density of the two-dimensional disk is simply taken from a vertical integrationof the initial condition of the MHD disk. Other hydrodynamic simulations are reruns of 10 –B2D.rin=0.8 at t = 500 with various locations of the inner boundary. The initial disks ofthe reruns are obtained by truncating the restart data of B2D.rin=0.8 to the desired radiuswhile keeping ∆ r/r and ∆ φ fixed. The properties of both the reruns and the control runcan be found in Table 1.
3. RESULTS
We present one 3D MHD simulation, called B3D, in this paper. This simulation wasterminated after the inner disk ( r < a ) reached a quasi-steady stage for several hundredtime units. Longer simulations require better resolution to resolve the MRI in the growingdensity concentration region (see § (cid:39)
480 code units, which corresponds to (cid:39)
77 binary orbits. The gas quicklyfills in the initially empty region within several binary orbital periods, and after ∼
100 unitsof time, the disk becomes fully turbulent. The binary torque then is able to maintain a lowdensity gap, and the inner part of the disk finally reaches a quasi-steady state after t ∼ ∼ § § § After the first 100 time units, the binary torque starts to clear out a low density gapbetween the inner boundary and r ∼ a . In the top-left panel of Figure 1, we show thevertically integrated surface density of the circumbinary disk at t = 120 in the x − y 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. Fromthe other three panels in Figure 1, we find that as the simulation continues the gap persists, We also performed a short duration rerun for t = 300–322 with higher dump rates: ten historydumps 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.
11 –the surface density gradually increases outside the gap, and finally an incomplete ring ofdense gas forms due to the combined effects of mass accretion and gravitational torque. Thestream also persists, but only one arm at a time. We will study the properties and effects ofthe transient stream in § t = 250 the disk appears to be elliptical and slightlyoff center, and at t = 350 and 450 the disk obviously becomes eccentric. We also observe agrowing azimuthal density asymmetry. We will discuss the eccentricity growth and densityasymmetry in § r = a ,1 . a , 2 a , 3 a and 4 a as a function of time in Figure 2. We find that the mass at r < a (i.e.,inside the initial disk’s inner edge) undergoes dramatic growth during the first ∼
100 timeunits as the initial disk fills in regions with < r min . The radial pressure gradient at the edgeof our initial disk leads to an inflow during the first several orbits, and once toroidal fielddevelops near the disk inner edge, the Maxwell stress quickly removes the angular momentumof the low density flow and drives inflow. As this transient phase passes, the mass in thegap levels off and stays quasi-steady (falling very slowly) after t ∼ t = 200 and the end of the simulation. On theother hand, the mass inside r = 3 a and 4 a rises by a factor of two after t = 100, indicatingthat this region does not achieve inflow equilibrium. The mass-interior profiles also possesshigh frequency fluctuations (typical time scale ∼ ∼
350 time units, there are sloweroscillations (time scale (cid:39) t = 250–350 (∆ T ) and t = 350–450 (∆ T ). We find there is only minor change inthe averaged surface density and accretion rate of the disk region < . a over ∆ T and ∆ T .If we define the unit of surface density Σ ≡ ρ a , the peak density grows in time, but onlyslowly, rising from (cid:39) . to (cid:39) . . Its position r p also changes, but similarly slowly,gradually moving from (cid:39) . a to (cid:39) a . The shape of the surface density profile likewisehardly changes in time. On the inside, the disk is truncated exponentially, (cid:104) Σ( r < a ) (cid:105) t ∼ Σ exp(3 . r/a − .
4) for both intervals, where (cid:104) (cid:105) t represents an average over time. Towardlarger radii, the surface density is ∝ r − , the predicted profile of a ‘decretion disk’ (Pringle1991). All that changes is that the region where Σ( r ) ∝ r − extends farther outward at latertimes. 12 –Averaged over the quasi-steady period, the accretion rate at the inner boundary is˙ M (cid:39) . GM a ) / Σ . However, the disk evolves throughout the simulation at r > a . Wesee gas inflow rates as large as ∼ . GM a ) / Σ outside the region ∼ a during both timeintervals. 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 surfacedensity peak. It is therefore convenient to define a ‘disk mass’ M d ≡ πr (cid:104) Σ p (cid:105) t , where r p = 3 a and (cid:104) Σ p (cid:105) t = 0 . . 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 (cid:104) ˙ M (cid:105) = 1 . × − M d Ω bin . (11)The surface density of such disks is very uncertain. A sense of scale, however, can be gleanedfrom translating eqn. 11 into the luminosity that would be produced if the accretion wereconverted into radiation at customary black hole efficiency (10%); in Eddington units, it is L/L E = 0 . M − / a / . τ p , where we have scaled to a total black hole mass of 10 M (cid:12) , abinary separation of 0.1 pc, and a disk column density whose Thomson optical depth is τ p . Unlike a steady accretion disk in a point-mass potential, where the internal stress simplytransports angular momentum outward and thereby drives an inflow, the binary consistentlyinteracts gravitationally with the circumbinary disk by torquing the surrounding gas. Theangular momentum delivered through these torques is transported by two mechanisms: MHDstresses, mostly due to MRI-driven turbulence; and Reynolds stresses associated with coher-ent gas motions in the gap region. We have checked that the code’s numerical shear viscosityis negligible compared to the transport mechanisms discussed here. In this section, we firstinvestigate the mechanisms of angular momentum transport and the properties of the binarytorque, and then discuss the balance between binary torques and the stresses.
In the first panel of Figure 4, we show the vertically-integrated and time- and azimuthally-averaged stresses and their sum as a function of radius. Within r < a , the character ofthese stresses is very similar in our two averaging periods, ∆ T and ∆ T . The Reynoldsstress shows big fluctuations radially. Its first peak, at r (cid:39) . a , is ∼ . c s . It dips ataround 3 a to nearly zero and then rises up again slowly. The ratio of the first peak to thesecond peak is ∼
10. However, the Maxwell stress has a flatter radial profile. It is roughly 13 – ∼ . c s on average and varies only slowly over a range of at most a factor of two. Itdiminishes within the gap and decreases smoothly towards larger distance. Comparing thetwo, we find the Reynolds stress exceeds the Maxwell stress by a factor ∼ r ∼ a .On the other hand, beyond the gap region, the MHD stress always dominates the Reynoldsstress by a factor (cid:38)
3. Within one binary separation, the Maxwell stress also exceeds theReynolds stress, mainly because the Reynolds stress falls with the lower gas density near theboundary. The total internal stress therefore follows the Reynolds stress at r (cid:46) a , whileit is close to the magnetic stress outside that region. Similar results have been reported inprevious studies on gap formation in protoplanetary disks with an embedded planet usingMHD simulations (Nelson & Papaloizou 2003; Winters et al. 2003). It is not a surprisingoutcome as the disk with a planet is just a special case of a circumbinary with an extrememass ratio.In the language of the α SS parameter of Shakura & Sunyaev (1973), we can also measurethe stresses in units of the pressure. In the right-hand panel of Figure 4, we present the timeaveraged stresses in these units. The most significant distinction compared with the absolutestresses is that the stress ratios increase steeply inward due to the low surface density in thegap. The total stress ratio peaks at α SS (cid:39)
11 at r = 1 . a . Its value in the disk body isalmost two orders of magnitude smaller.Our simulation shows that the circumbinary disk possesses highly time-dependent struc-tures in the horizontal plane, making it important to look at the instantaneous spatial distri-bution of the stresses in addition to the time-averaged values. Vertically-integrated snapshotsof both Reynolds stress (top left) and Maxwell stress (top right) taken at t = 305 time unitsare shown in Figure 5. On top of each plot, we also draw the surface density with 15 contourlines, logarithmically spaced from the density floor value to the density maximum. Nearlyall the Reynolds stress is confined within r < a . The most interesting point is that we findrelatively large stress in a single gas stream extending from ∼ a (see the red stripe inthe top left panel), where the gas gains angular momentum ( δv φ >
0) and is being pushedout by the binary. The gas being kicked thus goes on an eccentric orbit, creating negativestress at r ∼ a and − π/ < φ < π/ δv r > δv φ < π/ < φ < π/ δv r and δv φ < Having discussed the mechanism for angular momentum transport, let us now considerthe torque exerted on the disk by the binary. We plot the time-averaged torque density andthe integrated torque in Figure 6, where the torque is calculated approximately using thesurface density instead of the local density itself. The definition follows MM08: dTdr ≡ − (cid:90) Σ( r, φ ) ∂ Φ ∂φ rdφ (12)is the local torque, and the integrated torque is T ( r ) ≡ (cid:90) rr in dTdr (cid:48) dr (cid:48) , (13)where Φ = Φ( r, φ ) is the binary potential in the midplane.The first panel of Figure 6 shows the averages of the local torques for the two intervals∆ T and ∆ T are very similar, which again indicates that the inner part of the disk reaches aquasi-steady state. The local torque density is positive at r (cid:39) a and peaks at r (cid:39) . . a with dTdr (cid:39) . . GM Σ as the peak value. It then goes negative at r (cid:39) . a , reachingits first negative maximum at r (cid:39) . . a , but with a smaller magnitude (cid:39) − . GM Σ .Similar to what was previously found in MM08 and Cuadra et al. (2009), we find the torquedensity oscillates around zero but damps quickly toward larger distance. In addition, thetorque density for r (cid:46) a is negative because the gas within that region is advanced in phase(greater angular frequency) with respect to the binary, an effect that would appear in purelyhydrodynamic treatments as well as MHD.The integrated binary torque is displayed in the top right panel of Figure 6. Thetotal binary torque exerted on the disk is T ( ∞ ) (cid:39) . . GM a Σ . In order to makea comparison with the torque found in MM08, some normalization is needed. First, wenormalize the torque to the surface density at the radius where the local binary torque 15 –peaks. We denote that radius as r torq . In our MHD simulation, the first positive peakof the torque takes place at r torq (cid:39) . a while it is ∼ . a in MM08. In terms of thisnormalization, T ( ∞ ) / Σ( r torq ) ∼ . GM a in the hydrodynamic simulation of MM08. Inour MHD simulation, we find that it is (cid:39) . GM a (averaged in the quasi-steady period), afactor of 1 . α -viscosity treatment.However, as the peak of the torque is located in the gap region, where the surface densityrises rapidly with increasing radius, this means of normalization is unreliable. We thereforeturn to a different method, normalizing the binary torque to the disk mass. With thismethod, the normalized torque is T ( ∞ ) /M d (cid:39) . × − GM a − in our MHD simulation.However it is only ∼ . × − GM a − in MM08, ∼
14 times smaller. Just as for the massaccretion rate, there is a natural set of units for the total torque: it is convenient to describeit in terms of how rapidly the torque removes the binary’s orbital angular momentum. Inthat language, T ( ∞ ) = 2 . × − j bin M d Ω bin , (14)where j bin = ( GM a ) / / r p is roughly around (cid:46) a for both MHD and hydrodynamic sim-ulations, we can write the binary torque as some factor times GM a Σ p , and then we have T ( ∞ ) (cid:39) . × − GM a Σ p in our present simulation, while MM08 gives ∼ . × − GM a Σ p , ∼
14 times smaller. The largest part of this contrast can be explained by the fact that inthe MHD treatment Σ( r torq ) / Σ p (cid:39) .
15, a factor of ∼ α . In our case, α ∼ α SS (cid:38) . r (cid:39) a , a factor of (cid:38)
20 greater than the constant α = 0 .
01 assumed by MM08.In order to show how the local binary torque evolves with time, we also plot in Figure 6the specific torque density (bottom left), i.e., the absolute value of the ratio of torque densityto surface density. By dividing by the surface density, we remove the effects due to redistri-bution of the surface density. We find the specific torque slightly shifts to larger radius anddiminishes its amplitude from ∆ T to ∆ T . In the bottom right panel, we plot the historyof the total torque smoothed over a short time span ∼ T bin , where T bin = 2 π/ Ω bin is the bi-nary period. Before smoothing, the total torque fluctuates very rapidly ( ∼ . bin ) between ∼ − . GM a Σ and 0 . GM a Σ in the quasi-steady period. The power spectrum of thisfluctuation is dominated by a peak at ∼ . bin , approximately the beat frequency betweenthe orbital frequency of the matter being torqued most strongly (at the peak of dT /dr , i.e. r ∼ . . a ) and 2Ω bin . After smoothing, we find the binary torque has a sudden rise 16 –around t = 60 which is due to the initial disk filling the small radius region. The torquethen levels out between t = 100–200. Once in the quasi-steady stage, the torque shows aslowly decreasing trend, (cid:46)
30% in fractional terms over the final ∼
300 time units. It alsoappears to oscillate at the disk orbital frequency ( r ∼ a ) during the late time, a fluctu-ation caused by the eccentric movement of the density lump. As both the disk eccentricityand the density of the lump grow with time, the torque normalized by the time-dependentΣ p r simply reflects the increasing Σ p and r p . For a steady circumbinary disk, the angular momentum at a given radius should notfluctuate too much over time. Therefore the angular momentum deposition from the binarymust either be transported outward by Maxwell and Reynolds stresses or advected towardsthe hole. In order to test this idea, we perform an estimate of the differential angularmomentum flux averaged over a short time interval between t = 300–320 using the hightime-resolution rerun. Based on the conservation of angular momentum, we write the angularmomentum budget for a ring of disk between ( r − ∆ r , r + ∆ r ) as ∂J∂t + π ∆ r [( r + ∆ r ) F ( r + ∆ r ) − ( r − ∆ r ) F ( r − ∆ r )]= π ∆ r [( r − ∆ r ) W rφ ( r − ∆ r ) − ( r + ∆ r ) W rφ ( r + ∆ r )] + dTdr , (15)where J ( r, t ) = 2 πr (cid:104) ρv φ r (cid:105) L z is the shell-integrated angular momentum, F ( r, t ) = (cid:104) ρv r (cid:105)(cid:104) ρv φ (cid:105) L z / (cid:104) ρ (cid:105) is the mean angular momentum flux due to mass motion, and W rφ ( r, t ) = T rφ + R rφ is thetotal internal stress. We write the conservation law this way so that a negative differentialflux due to advection (the second term on the left hand side of equation (15)) indicates a netangular momentum inflow, while a negative flux due to internal stresses (the first term onthe 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 ofradius. In its y –legend, we use dF J /dr to represent all the terms in the equation, which areall radial derivatives of one or another variety of shell-integrated angular momentum flux( F J ).During this time interval, the local angular momentum stays steady, with only verysmall temporal variations from the averaged ∂J/∂t (solid black curve), as expected for aquasi-steady accretion disk. The binary torque (blue dashed curve) at r ∼ a is almostbalanced by the differential flux due to the sum of Reynolds stress and Maxwell stress (cyandashed curve). In other words, the circumbinary disk within that region reaches a balance 17 –between binary torque and internal stresses. As a result, the angular momentum flux carriedby accretion flows oscillates about zero (magenta curves within 1–3 a ). Clearly, the internaltorque due to the Reynolds stress (green dash-dotted curve) dominates the Maxwell stress(red dash-dotted curve) at r ∼ a , which suggests that most of the angular momentumdumped by the binary at r ∼ . a is transported outward via the Reynolds stress, whilethe angular momentum drawn from the binary at r ∼ . . a is compensated by the torqueof the Reynolds stress. Finally, we comment that torque balance is not reached at r (cid:46) a (where negative binary torque dominates) or r (cid:38) a (where internal stress dominates). Bysumming all the other torques, we find the differential angular momentum flux that mustbe advected inward in those two regions (magenta curve). Directly computing the secondterm on the left hand side of equation (15) yields the (very similar) black dashed curve inFigure 7.To sum up, we have found three principal conclusions about the angular momentumbudget. First, the Reynolds stress dominates the Maxwell stress in the gap region, while thelatter exceeds the former outside that. Secondly, the distribution of both stresses followsthe gas streams and are vertically confined: within two scale heights for Reynolds stressand four scale heights for Maxwell stress. Third, we find the radial profile of the binarytorque 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 torquein 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 greaterthan that found in MM08. Lastly, we find that in the quasi-steady state, the binary torqueis roughly balanced by the torque due to the Reynolds and Maxwell stresses in the innerdisk ( r < a ). Most of the angular momentum dumped by the binary in the gap region istransported outward by the torque of the Reynolds stress associated with the gas streams. Although we begin our simulation with an axisymmetric disk configuration, which meansthere is zero eccentricity at the start, the disk eccentricity grows throughout our simulation(e.g., compare the disk surface density snapshots in Figure 1). 18 –To quantify the eccentricity, we define the local eccentricity as e ( r, t ) = 1 T bin (cid:90) t + T bin t (cid:12)(cid:12)(cid:10) ρv r e iφ (cid:11)(cid:12)(cid:12) (cid:104) ρv φ (cid:105) dt (cid:48) , (16)where (cid:104) (cid:105) is the shell average defined as in equation (3). We show this local eccentricity inthe space-time diagram in Figure 8 (top panel). The local eccentricity of the gas in the gap( r (cid:46) a ) rises significantly during the period t ∼ e ( r ) in the gapsimply reflects the fact that the dynamics there are strongly non-Keplerian. The eccentricityin the main body of the disk ( r > a ) grows rapidly over the period t = 100–300. Thetypical value increases from much less than 0 .
01 throughout the disk body at t = 100 to ∼ .
08 between r = 2 –3 a at t = 300. It keeps growing slowly after 300 time units, andat the end of the simulation a ring of eccentric disk forms at 2 a < r < . a with typical e ( r ) ∼ .
1. A more detailed picture of time-variation in the eccentricity is shown in thebottom panel of Figure 8, where we show its evolution from t = 300 to 320. For clarity,we use the data from the high time resolution simulation, and do not take the temporalsmoothing as in equation (16). We find strong radial extending structures which connect thegap with the disk body. These features bend at ∼ . a 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. Wewill discuss the stream effects on the disk eccentricity in § t = 150–250(solid black), ∆ T (red dotted) and ∆ T (blue dashed).The distribution curves in the disk body ( r > a ) shift toward larger radius and greatereccentricity as the disk evolves in time. Outside the peaks around 2 a , the slopes of thedistribution curves are constant. The radial dependence of e ( r ) in the disk body can beroughly described as ∝ exp [ − . r/a )].Interestingly, the eccentricity distribution found by MM08 was quite similar to ours inboth amplitude and shape, despite that calculation’s very different internal stresses, initialsurface density profile, and duration. In our own 2-d hydrodynamical simulations, we havelikewise found qualitatively similar eccentricity growth. However, the eccentricity profilein those simulations matches the MHD simulation (and MM08) only at a specific time( t ∼ − ). From these comparisons, it is clear that the primary mechanism for drivingeccentricity must be hydrodynamic response to gravitational forcing by the binary, but itsquantitative development can be affected both by disk physics (e.g., MHD stresses) and initialconditions (e.g., the initial surface density profile). We will discuss specific mechanisms atgreater length in § r = 2 a and r = 4 a as a measure of the disk eccentricity. InFigure 10, we show the averaged disk eccentricity e disk constructed by taking volume averagesfrom r = 2 a to r = 4 a and time averaging over one binary period: e disk = 1 T bin (cid:90) t + T bin t dt (cid:48) (cid:12)(cid:12)(cid:12)(cid:82) a a dr (cid:82) ρv r e iφ r sin θdθdφ (cid:12)(cid:12)(cid:12)(cid:82) a a dr (cid:82) ρv φ r sin θdθdφ . (17)At early times ( t (cid:46) t = 100to t = 250. By fitting the eccentricity curve during that interval, we find the growth rateis γ e (cid:39) . bin , where the subscript e denotes eccentricity. The growth rate of diskeccentricity slows by a factor of 4 after ∼ To follow changes in orientation of the eccentricity, we define the longitude of the apoapse (cid:36) in terms of the complex phase angle of the disk eccentricity for 2 a ≤ r ≤ a . Itstime evolution is shown in Figure 11. Here the angle (cid:36) measures the angular location ofthe apoapse with respect to the x –direction in the inertial frame. The disk eccentricity isvigorously 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 thistime interval is very small. The angle (cid:36) then gradually rotates in a prograde manner. Byfitting the curve between t = 200–480, we find ˙ (cid:36) (cid:39) . × − Ω bin . Linear perturbationtheory predicts the precession rate for particles around the binary potential is ˙ (cid:36) (cid:39) Ω − κ for a cold disk, where Ω is the angular frequency and κ is the epicyclic frequency of the disk.To a first order approximation, Ω( r ) ≈ Ω K (cid:104) (cid:0) ar (cid:1) q (1+ q ) (cid:105) (also see equation (2)), and κ ( r ) ≈ Ω K (cid:104) − (cid:0) ar (cid:1) q (1+ q ) (cid:105) . Thus we find the prediction is ˙ (cid:36) (cid:39) . a/r ) / Ω bin for q = 1,or (cid:39) . × − Ω bin evaluated at r (cid:39) a . The precession rates measured in both our MHDsimulation and MM08 are close to this prediction. We therefore suggest that the precessionis mainly due to the quadrupolar potential. 20 – There is another disk characteristic closely connected to the eccentricity: the densityconcentration 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 theeccentric disk. To quantitatively study the asymmetric feature we introduce the m = 1 modestrength component S m=1 ( r, t ) ≡ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ρe imφ rdθdφ (cid:12)(cid:12)(cid:12)(cid:12) ( m = 1) , (18)which is basically the m = 1 component of the Fourier transform of the surface density. Weplot the space-time diagram of S m=1 in Figure 12. The orbiting lump appears as a zig-zagpattern whose peak moves in and out between r ∼ a and 4 a after about 300 time units.The pattern period, which is also the orbital period of the eccentric disk, is ∼
30 time units.We will discuss the lump in § As shown in Figure 1, gas is not entirely absent from the gap. Consistent gas flowslaunch from the inner edge of the disk and stream toward the binary. The streams not onlyaffect the accretion rate of the circumbinary disk, but also the structure of the inner part ofthe disk.In a frame that co-rotates with the binary, we find the streams follow roughly thestatic bi-symmetric potential of the binary and pass through the inner boundary about 0 . . L L r (cid:46) a in § t = 122 and125 in the first two panels of Figure 13, and draw the velocity field vectors measured in theco-moving frame on top of it. The snapshots at t = 122 show the inflowing gas is regulatedby 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 tothe next, a direct signal that the binary torque changes its sign from one quadrant to theother. In the co-rotating frame of the binary, the gas in the second (near φ = π ) and fourth(near φ = 0) quadrants is pulled toward the closer component of the binary (see t = 122snapshots). The tangential velocity results in a Coriolis force, which tends to move matteraway from the center. It takes about half the binary period for the gas streams in those two 21 –quadrants to be kicked out and strongly interact with the disk edge. As shown in the latersnapshots at t = 125, the velocity is mostly outward at r ∼ . a in those two quadrants, andthere are strong density enhancements near the regions where the velocity gradient is large.There is also an equivalent way to think about the transition between t = 122 and t = 125.In the inertial frame, the binary almost always rotates faster than the disk; therefore, thesecond and fourth quadrants have positive binary torque. As a result, gas in these regionsgains 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 diskbecomes significant. Instead of a pair of streams with comparable size and strength, we find,for instance at t = 301, the stream on the right dominates the left. For an eccentric disk,the apocenter side of the disk is far from both members of the binary. Thus, the gas streamis strongly torqued, and therefore carries more mass, only when one of the members getsclose to the pericenter of the disk. After half a binary period, the binary members switchtheir phases (in the inertial frame) and a strong stream forms around the other member. Weillustrate this process with a time sequence of snapshots at t = 301,302, 303 and 304 . t = 301 is gradually pushed away at t = 302and 303 and joins the apocenter of the eccentric disk. Meanwhile, the stream associatedwith the other member of the binary strengthens as that mass approaches the pericenter at t = 304 . | f g + f L | / |∇ Φ | is always (cid:46) /
3, where f g = −∇ P/ρ is the force density of the gas pressure gradient and f L = ( ∇ × B ) × ∇ B / πρ is the Lorentzforce density. Because gas pressure and magnetic forces are smaller than the gravitationalforce, the trajectory of the gas is essentially ballistic. Consequently, we neglect the fluid andmagnetic effects when calculating the interaction between the streams and the disk inneredge in § We found in the previous subsection that the Maxwell stress closely follows the gasstreams in the x – y plane and is well confined near to the midplane. Now we show the fieldstructure in the gap region, especially in the stream, has similar features. In Figure 14, wedraw the vertical averaged magnetic field as vectors on top of the logarithmic scaled surfacedensity at t = 300 .
5. The location and length of the arrows show the position and relativestrength 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 22 –field in the upper stream (top half plane). The ordered field in the gap therefore producesgreater Maxwell stress only within the two streams.
We present the accretion rate through the inner boundary as a function of time in thetop panel of Figure 15. The sampling rate of the time sequence is determined by the outputrate 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 dashedcurve, which is derived by smoothing the time sequence using a boxcar average with width15 time units). These two time variations are the most significant modes in the temporalprofile of ˙ M . By constructing the power spectral density of the accretion rate using Fouriertransforms (middle panel, Figure 15), we are able to identify them: The higher frequencymode is at about twice the binary rotation rate, and it is the same rate at which the streamsare pulled inward by the binary potential; the lower frequency mode is ∼ . bin , whichis the same as the orbital frequency of the lump during the later stages of the simulation(as shown in § bin if q = 1; MM08;Hayasaki et al. 2008; Cuadra et al. 2009; R¨odig et al. 2011; Sesana et al. 2011), a lowfrequency component due to the motion of the dense part of the disk (MM08; R¨odig et al.2011), and another component (in our case at 1 . bin ) created by a beat between the binaryorbital frequency and the disk orbital frequency (R¨odig et al. 2011). The principal contrastbetween our results and this earlier work is that the dominant frequency is 2Ω bin rather thanΩ bin because the members of the binary have exactly equal mass.To further verify the cause of these two modes, we pick out two times ( t = 400 and 441as marked in the top panel by black arrows) and plot their two-dimensional local accretionrates in the inertial frame in the bottom panel of the same graph. There were accretionoutbursts 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 verticallyintegrated accretion rate. The black contour lines show the surface density in a logarithmicscale from 10 − to 10 − . . These lines pick out the location of the streams, while the whitecontour lines show the surface density between 1 . . M are indeed coming from the infallingstreams located in the fourth quadrants of both snapshots. However, the local accretion ratein 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 23 –higher density than other times. Thus we conclude that the time dependent accretion showsboth high frequency outbursts corresponding to the periodic inflow from the streams and lowfrequency modulation due to the orbiting lump. The low frequency mode becomes importantonly when the lump forms and maintains itself after about t = 350.The time averaged accretion rate at the inner edge, normalized with the peak density,is ∼ . √ GM a Σ p , a factor ∼
40 greater than found in MM08. Interestingly, during thequasi-steady period, the accretion rate in the gap is nearly twice the rate measured at r p ,where the surface density peaks. This contrast explains the slow mass depletion inside thegap shown in Figure 2. One of the most important questions about accretion onto a binaryis the ratio of the accretion rate accepted by the binary to the rate at which mass is fed intothe disk from the outside. If the latter rate can be estimated by the maximum accretionrate in the disk, we find that accretion through the inner boundary of the simulation, andtherefore presumably onto the binary, is about 1/3 of the accretion rate supplied. This ratiois a factor of two greater than previously found by MM08. However, in evaluating thesenumbers one should be aware that the disk as a whole never reaches equilibrium over thecourse of the simulation. That the time-averaged accretion rate as a function of radius isfar 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 estimateof 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 MHDsimulation. However, our hydrodynamic simulations with different cut-outs suggest that therate hardly changes as long as the cut-out size is smaller than (cid:46) a .
4. DISCUSSION4.1. Binary Torque and Linear Resonance Theory
Using the results we have on the disk stresses and the binary torque, we want to furtherelucidate the angular momentum transport mechanisms within a MHD turbulent circumbi-nary disk. The binary strongly torques the surrounding gas in the gap as the gas formsstreams; the outward-moving portions of the streams deliver that angular momentum to theinner part of the circumbinary disk. By contrast, the torque directly exerted on the diskbody is much weaker and diminishes rapidly toward larger distance. In that sense, we cansay that the binary mainly torques the gas in the low density gap; moreover, the magnitudeof that torque is proportional to the mass of gas in the gap. This fact suggests that thebinary torque is very sensitive to the basic physics coupling the streams to the disk such asthe stresses and the equation of state. For instance, scaled with the peak surface density, 24 –our MHD turbulent disk provides a factor of 14 times stronger binary torque than that ofMM08.This fact also leads to a strong contrast with the commonly-used linear resonance theoryof interaction between disks and planets or other orbiting satelites (Goldreich & Tremaine1982). To illustrate that contrast, we calculate both the torque density and total torqueas functions of radius, assuming most of the torque comes from the Ω bin : Ω = 3 : 2 outerLindblad resonance at r /a = (3 / / (cid:39) . / d ln Σ /dr ) − (cid:46) . a ,which is comparable to the first wavelength 5[( m + 1) / (3 m )] / ( H/r ) / a ∼ . a , where H/r is evaluated at the resonance radius, and m = 2 is used here (see Figure 1 and the pressure-dominated case in Table 1 of Meyer-Vernet & Sicardy (1987)). Because the surface densityat r is so much smaller than at r ∼ a , one might expect the maximum torque to moveaway from the resonance point toward the inner disk. Secondly, linear theory assumes nearlycircular orbits because it also assumes q (cid:28)
1. However, in our simulation q = 1, and, largelyfor this reason, the orbits of gas within the inner edge of the disk body are significantlynon-Keplerian: they form streams. Because the torque depends so strongly on whatever gasmass is in the gap, this distinction is qualitatively important. Lastly, as the mass ratio ofthe binary of our simulation is unity, the forcing term, which is roughly controlled by theratio of the perturbing gravitational potential to the sound speed, becomes large enough todrive nonlinear density disturbances. The longer wavelength oscillation found in our torquedensity might be partly explained as a result of these nonlinear effects (Yuan & Cassen1994). The time-averaged rate at which the binary loses angular momentum to the disk is T ( ∞ ) (cid:39) . GM a Σ . Meanwhile, the time-averaged angular momentum flux across the 25 –inner boundary due to accretion is ˙ M j in (cid:39) . GM a Σ , where j in (cid:39) . √ GM a is theaveraged specific angular momentum at r in . On average, therefore, the binary actually gainsangular momentum in net.However, even though the binary gains angular momentum, it does not automaticallyfollow that the binary separation grows. The evolution of the binary also depends on itsgrowth in mass. If the orbit remains circular, the shrinkage rate for an equal mass binary is˙ aa = 2 ˙ JJ − MM , (19)where J = M j bin = 0 . M √ GM a denotes the angular momentum of the binary, and itsrate of change is ˙ J = ˙ M j in − T ( ∞ ).After separating the accretion term of ˙ J and combining it with the second term inequation (19), we have ˙ aa = − T ( ∞ ) J + 2 ˙ M j in J (cid:18) − j bin j in (cid:19) . (20)If j in /j bin > / (cid:39) . − . (cid:112) Ga/M Σ and 0 . (cid:112) Ga/M Σ respectively.The net result is orbital compression, but at a rate much smaller than either of the two terms,˙ a/a (cid:39) − . (cid:112) Ga/M Σ .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˙ aa = 2 ˙ JJ − MM + (cid:18) e − e (cid:19) ˙ e bin e bin , (21)where e bin denotes the binary’s orbital eccentricity. In addition, the angular momentum ofthe binary should also be adjusted to account for the eccentricity: J = M j bin (cid:112) − e . InAppendix A, we show that the rate of eccentricity growth depends on how the accretionflow interacts with the interior disk of each individual binary member, an interaction oursimulation did not treat. However, we also show that the best estimate we can make fromthe data we do have is that there is little evidence for any significant ˙ e bin /e bin . Appendix Aalso 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, andlinear theory may not apply in this regime. If, as we did in § M d = π Σ p (3 a ) , and further assume that the binary orbit remains circular, then 26 –˙ a/a (cid:39) − . × − ( M d /M )Ω bin . For comparison, MM08 found a rate (cid:39) − . (cid:112) Ga/M Σ ,which can be scaled used the definition of M d to (cid:39) − × − ( M d /M )Ω bin , somewhat slowerthan ours. In other words, although we find a torque an order of magnitude larger than thatof an α –viscosity hydrodynamic model, we also find an accretion rate so much larger thatthe net effect on the binary orbit is reduced to an increase by only a factor of 2.7.Another way to understand this result is to note that our T ( ∞ ) / ˙ M (cid:39) . GM a ) / ,a factor of 3 smaller than in MM08. That is, the greater internal stresses produced byMHD effects, especially in the gap region, both introduce more mass to that region (therebymagnifying the torque) and lead to even greater accretion, which returns angular momentumto the binary. Unfortunately, because the cancellation is so close, our calculation cannot bedefinitive in regard to the shrinkage rate of comparable-mass binaries in Nature. Furtheruncertainty comes from the fact that we cannot track the gas flowing into the cutoff regionwith the present model, so we do not know exactly how much energy and angular momentumaccretion might bring to the binary, yet our sizable accretion rate makes these contributionssignificant. The net outcome in any particular case may therefore depend on the specificdetails (gas equation of state, binary mass ratio, etc.), as all of these can influence j in , ˙ M , T ( ∞ ), as well as the detailed mechanics of how the stream joins the interior disks associatedwith the members of the binary. The distribution of eccentricity in radius plays an important role in the evolution of thedisk. For the eccentricity to grow substantially over time in a large circumbinary disk, it needsto be confined in radius. Otherwise, the energy and angular momentum (more precisely,angular momentum deficit) associated with an eccentric disturbance are spread over a largeradius, resulting in a small eccentricity in any one location. We consider here a linear modelfor the eccentricity distribution and compare it with the results of the simulations.The linear equation governing the eccentricity distribution and evolution is given byGoodchild & Ogilvie (2006, hereafter, GO06). This also allows generation and damping ofthe 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 theeccentricity was small and growing exponentially, i.e., the linear regime. We define the latestage of the exponential growth phase as the period t = 200 −
250 time units (see Figure 10).Surface density Σ( r ) is taken directly from the simulation by averaging over this time interval. 27 –The disk inner radius is taken to be r i = 2 a and the outer radius r o = 10 a . We assume thesame equation of state as in the simulation: isothermal, with sound speed 0 . bin a . Theadiabatic index is taken to be γ = 1 − . i . The small imaginary part (which introduces slowdamping) assures mathematical consistency. We model the forces driving the initiation ofeccentricty as a Gaussian in radius (see equation (B11)), centered on r c = 2 . a and havinga width w = 0 . a . These parameters are not well constrained by the simulation. It followsthat the the results are also insensitive to their values within certain limits. Parameter s c , thepeak value of the eccentricity injection distribution in equation (B11), is determined by thecondition that the growth rate equal that in the simulation, γ e = − Im( ω ) (cid:39) . bin , where ω is the eigenvalue defined in equation (B14). The value is found to be s c = 0 . a Ω bin .The physically relevant solution on long timescales corresponds to the fastest growingmode. This solution was determined by the methods described in Lubow (2010, hereafterL10). Figure 16 plots the eccentricity distribution for this mode, along with the eccentricitydistribution obtained from the simulation. The middle dashed line corresponds to the eccen-tricity distribution based on linear theory with the set of parameters described above. Thedistribution matches well the one obtained from the simulation (the solid line). This distri-bution closely resembles the distribution for the fundamental free precession mode in the disk(the uppermost dashed line). That mode has an eccentricity injection of zero, s c = 0. Theconfinement of eccentricity is due to the effects of differential precession of the binary ∂ r ˙ (cid:36) g (see equation (B8) in Appendix B), while pressure effects act to spread the mode. Differen-tial precession effectively creates a potential well that leads to trapping of the fundamentalmode. Further localization is due to the injection of eccentricity, which is concentrated nearthe disk inner edge (solid curve on lower left). The faster eccentricity is injected into a smallregion, the more the eccentricity distribution narrows. For somewhat higher eccentricityinjection rates, the distribution is further confined (lowest dashed curve). Therefore, we seethat the eccentricity distribution is a consequence of the excitation of the fundamental freeeccentric mode in the circumbinary disk and that the eccentricity is confined within a radialwidth ∼ a .The excitation of the fundamental mode occurs because that mode overlaps well withthe eccentricity injection distribution. As discussed above, the location and width of theeccentricity injection are not well determined by the simulations. On the other hand, thelinear equation has solutions that match the simulation results only if r c (cid:46) . a , whichindicates that the source of the eccentricity should be located near the inner edge of ourMHD disk. For larger values of r c , there is no value of s c that could provide the growth ratefound in the simulation for an eccentricity distribution that resembles the simulated one.This radius limit has a weak dependence on the width of the injected eccentricity w . If thewidth is doubled to 0 . a , then the limiting radius increases by about 1%. 28 –The outer regions of the simulation exhibit an exponential falloff of eccentricity withan e-folding length of about 0 . a (see Figure 9 in § r (cid:29) a ) using linear theory. It isdominated by an exponential falloff in r with an e-folding length also of about 0.8a.To sum up the analysis in this subsection, we find the linear model can explain theeccentricity distribution of our simulation. It shows the confinement of the disk eccentricityis mainly due to the disk precession and eccentricity growth. The linear model also putsconstraints on the location of the source for eccentricity growth. Finally, it explains theexponential falloff of eccentricity profile. One possible mechanism for eccentricity growth is a tidal instability. This instability isbelieved 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 thatsurrounds a star (Lubow 1991, hereafter L91). As pointed out in L91, this instability canalso be at work in circumbinary disks at radial locations that satisfy the conditionΩ( r ) = m Ω bin m + 2 , (22)where m is an integer that corresponds to the azimuthal wavenumber of the tidal componentof the potential involved in the instability. In the case analyzed here of an equal mass binary, m = 1 is excluded because that component of the tidal potential is absent. For m = 2, theinstability is associated with the Ω bin : Ω = 2 : 1 resonance that occurs at r (cid:39) . a . Inthis model, the m = 2 tidal field interacts with the m = 1 eccentric motions of the disk toproduce disturbances of the form exp (3 iφ − i Ω bin t ). These disturbances launch waves ofthat form at the 2 : 1 resonance. The waves in turn interact with the tidal field to producea stress that increases the disk eccentricity exponentially in time.However, the 2 : 1 resonance lies at r (cid:39) . a , inside the circumbinary disk gap. Further-more, the gas motions are far from circular there. Even if the gas followed circular orbits inthis region, the model of L91 predicts in this case an eccentricity growth rate ∼ a Ω bin /w e ,where w e is the width of the eccentricity distribution and is ∼ a . This very rapid growthrate is a consequence of the powerful m = 2 tidal potential for an equal mass binary, thedominant component of the tidal field. However, this rate is more than an order of magni-tude faster than the growth rate measured in the simulation. It may be possible that theeffects of the resonance are still felt, but at a reduced level, near the disk inner edge as a 29 –consequence of finite width of the resonance. Thus, tidal instability may play a role in theeccentricity growth, but the evidence suggests it is at most a limited role. As discussed in § ∼ a toward larger radius. In this sectionwe examine the possibility that the impact of streams in the gap striking the inner edge ofthe disk is the mechanism of eccentricity injection into the disk.During the exponential growth phase of eccentricity, a pair of gas streams are pulledin from the disk inner edge by the binary, and then flung back to the disk after half binaryperiod. It is easy to show that if the disk is on a circular orbit and free of eccentricity, the twostreams are not capable of breaking the bisymmetry that is intrinsic for an equal mass binarypotential. We speculate that a small eccentricity breaks the symmetry by inducing unequalstrength stream pairs (in terms of both the density and the velocity) and/or asynchronousstream-disk interaction. Once the symmetry of stream impact is broken, stream impactcould amplify the initially small disk eccentricity, in turn increasing the asymmetry of thestreams themselves. This process could then lead to sustained disk eccentricity growth.We find two sorts of evidence that support this mechanism of eccentricity growth. Thefirst is based on a special set of simulations we performed. In these, we eliminated part or allof the returning streams by enlarging the central cut-off. If the removal of outward streamshalts eccentricity growth, we may take that as evidence in support of a stream impact originfor eccentricity excitation. We do not consider cutoff radii large enough to affect the circularorbit region of the disk because otherwise, the result might be equally well interpreted interms of a resonance model ( § § a , similar to B3D.In the gap region, there is a pair of streams at t = 300 −
700 while the disk eccentricitygrows exponentially. A single strong stream appears when the disk eccentricity becomessignificant. The simulation is terminated at t = 1500, and at that time the disk eccentricityis saturated. We find the stream dynamics and the eccentricity growth in the control run arevery similar to what were found in B3D, demonstrating that hydrodynamic effects alone canyield eccentricity qualitatively similar to MHD. Thus, the results based on hydrodynamicsimulations can provide a clear indication of the role of the streams in eccentricity generation. 30 –We display the growth history of the disk eccentricity in Figure 17 (black solid). Here thedisk eccentricity is calculated using a two dimensional definition of equation (17). Similarto the MHD simulation, we find exponential growth from t ∼ (cid:39) . bin , very close to that of the MHD simulation ( (cid:39) . bin ). Wethen restart from t = 500 with five different values of r in while keeping all other parametersfixed. We terminate the reruns at t = 1500, the stopping time for B2D.ri=0.8.In Figure 17, we plot the evolution of disk eccentricities for cases with different r in .The case with r in = 1 . a (blue dash-dotted) loses only the tips of the streams, and itseccentricity evolution closely follows the control run, with exponential growth continuingwithout any interruptions until it becomes saturated. The disks with larger cut-offs ( r in =1 . a : green dashed curve; 1 . a : red dotted curve), however, cease exponential growth oftheir eccentricities right after the restart because most of the streams forming in the gapare lost in the hole and never get a chance to interact with the disk edge. There are alsocases in between. When r in = 1 . a (cyan dash dot dot curve) and r in = 1 . a (magenta longdash curve), the growth rate of the eccentricity diminishes after the restart. It takes longertimes for both disks to reach the eccentricity of the control run. We attribute these resultsto partial loss of outward stream motion.When r in = 1 . a , the simulation domain includes the region well inside the gap wheregas motions are noncircular and are dominated by the streams. In this case, the simulationregion also extends well inside the 2:1 resonance at r (cid:39) . a whose presence is required forthe mechanism of § t = 300, when the eccentricity growth is slow, the inner portion of thedisk forms an eccentric ring between r ∼ a with e ∼ .
08. Its pericenter is at r ∼ a and φ ∼ π/
2, and the apocenter is at ∼ a and φ ∼ π/ dedt = √ − e a d n (cid:20) R sin f + S f + 3 e + e cos(2 f )2(1 + e cos f ) (cid:21) , (23)where a d is the semi-major axis of the eccentric ring of disk and n = (cid:112) GM/a is the meanmotion of the disk. Quantity f is the true anomaly where the stream impact occurs. R and S are the radial and angular components of the disturbance force density, and are definedas ( R, S ) = ˙ m s ( v s − u e ) m e , (24)where ˙ m s denotes the mass injection rate of the stream, m e is the total mass of the eccentricring, v s and u e are the stream and disk velocities near the impact, respectively.Since the impact takes place over a short time, we must use the higher time resolutionsimulation ( t = 300–322) in order to estimate the disk perturbation. We carry out ouranalysis using the vertically-integrated two-dimensional data. For the Gauss equation toapply, the properties of the eccentric ring, such as e , a d , n , and the longitude of pericenter (cid:36) , need to be nearly constant over the orbital period of the ring. This assumption holds wellbecause the timescale for the ring properties to change is much longer. In addition, the masstransferred to the ring by the stream impact over an orbital period is considerably smallerthan the mass of the eccentric ring. We use time-averaged values for those parameters. Thevalues are: e = 0 . a d = 2 . a and (cid:36) = 1 .
48 radians. The stream-disk impact is localizedin space. We analyze the stream-disk impact over an arc that is defined by r = 2 a and φ ∈ [ − π/ , π/ f and v s are directly measured on the arc by taking thedensity weighted averages at each time. The mass injection rate is calculated by integratingthe mass flux along the arc. The instantaneous velocity of the eccentric edge u e is obtainedfrom the unaffected matter at a slightly greater radius, which we take to be an arc at r = 2 . a .In fact, the disk velocity is not very sensitive to the location we choose, so long as it is takenwithin the eccentric ring.Using these parameters, we calculate the rate of change of the eccentricity from equa-tion (23). The result is presented in Figure 18. We find that each impact lasts for less thanone time unit. The impacts induce peak values of de/dt that range between 4 . × − Ω bin and 0 . bin at each half binary period. Most of the time, de/dt remains close to zerobecause stream-disk impacts are so brief. Taking time averages of the curve, we obtain anestimate of the impact-induced rate of change of eccentricity (cid:104) de/dt (cid:105) t (cid:39) . × − Ω bin . Thecorresponding growth rate γ e (cid:39) . bin is consistent with the growth rate at earlier times 32 –in the simulation ( (cid:39) . bin ). The growth rate contribution predicted by the stream im-pact at this later time may not be quite the same as the impact contribution at the earlierstage. Nonetheless, the approximate correspondence is encouraging. Indeed, one questionthat 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, maylimit eccentricity growth. The disk after t = 400 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. Asdiscussed in § stream reabsorption . When the lump approaches pericen-ter, the binary peels off a gas stream that is denser than streams drawn from lower densityregions. This relatively dense stream is then kicked back out by the binary torque. Althoughthe azimuthal location reached by the returning stream may not be exactly where its start-ing 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 compressesit, restoring the density loss incurred by shearing during the time the stream passed nearthe binary. The left panel in Figure 19 shows the moment an episode of stream reabsorptiontaking place.The second mechanism makes use of streams drawn from regions of the inner disk otherthan the lump. We call it lump feeding , as this channel of development not only enhancesthe concentration of the lump, but also increases its mass. After t = 400, the densityenhancement orbits eccentrically, and as it follows its orbit toward apocenter, both v r and v φ diminish. Meanwhile, streams traversing the gap continue to strike the disk inner edge, 33 –creating a new, smaller density enhancement at orbital phase angles behind the lump. Theorbital speed of this new density concentration is much greater than that of the main lumpbecause it is at rather smaller radius. As a result, the newly-created and smaller lumpcan catch up with the main lump and join it before the lump returns from apocenter andaccelerates. In the right panel of Figure 19, we show an example of this process, picturedat a time ( t = 422) shortly before an outgoing stream reaches the lump. Three time unitslater, the small enhancement (the green region close to r = 2 a at φ (cid:46) π/
2) joins the mainlump (green and red colored area in the second quadrant). The density-weighted velocity inthe enhancement is ( v r , v φ ) ∼ (0 . , . bin a (we average the stream gas where the surfacedensity is greater than 1 . ). However, the lump has a slower velocity ∼ (0 . , . bin a (averaging over locations where Σ > ). As a result of both the lump feeding and streamreabsorption, the density lump grows both in mass and density contrast.Another way to illustrate the connection between streams and the lump is to look atthe space-time diagram in Figure 12. In the diagram, the stream clearly stands out as theradially stretched features between r (cid:46) . a and r (cid:38) a which propagate away from thebinary with a pattern speed of about 0 . bin a . This pattern repeats itself at a rate abouttwice the binary frequency. Meanwhile, the zigzag oscillation at late time ( t > ∼ a and 4 a shows the movement of the density lump as it orbits around the binaryon an eccentric orbit with a period of ∼
30 time units. Clearly, when the lump passesthe pericenter, the stream reabsorption process promotes growth in the density contrast,driving an increase of the mode strength at r = 2–2 . a 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 nowmuch further away from the stream in azimuthal angle, and thus we find the descending partat r (cid:38) . a 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 (lowmode strength regions) around 2 . a between the zig-zags actually show the moments thatthe newly kicked out stream reaches the radial coordinate (but is distant in azimuth) of thelump, so that it smooths the m = 1 mode at that radius. The zig-zag feature’s growth inboth 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. 34 –
5. CONCLUSIONS5.1. Specific Results
In this work, we have performed the first three-dimensional MHD simulation of a cir-cumbinary disk around an equal-mass binary, which in this case follows a circular orbit. Themain results are:1. The disk exhibits a number of nearly steady features. There is a low-density gap within r ≤ a , an eccentric inner disk at (2–3) a . At early times, there is a pair of gas streamsthat flow into the gap from disk inner edge. They are nearly steady in the corotatingframe of the binary. At later times, there is a single dominant stream whose mass fluxis time varying with a period of half the binary orbital period. Parts of these streamsare torqued so strongly by the binary that they return and impact on the disk inneredge.2. Some aspects of the disk evolve secularly. At late time, the disk inner edge developsan asymmetric density concentration (‘the lump’) whose mass, density contrast, andorbital eccentricity grow steadily. We find the lump is due to a combination of streamimpact and disk eccentricity.3. The disk eccentricity grows exponentially during the time t = 100–250 (16–40 binaryorbits) with a growth rate of ∼ . bin . The growth rate then slows down signifi-cantly. By the end of the simulation, the disk body at 2 a (cid:46) r (cid:46) a has reached aneccentricity ∼ .
1. Its pericenter precesses slowly due to the quadrupolar componentof the binary potential. Stream impact largely accounts for the eccentricity growth.4. Reynolds stress associated with the streams is the leading transport mechanism in thegap region, while Maxwell stress dominates in the disk body.5. Although the profile of the binary torque is similar to that of previous hydrodynamicsimulations, its magnitude depends strongly on the magnitude of internal disk stresses.Normalized to the disk mass near the density peak, we find a measured total torque T ( ∞ ) (cid:39)
14 times greater than found by MM08. Because the torque at any particularradius is proportional to the gas mass there, the integrated torque is directly propor-tional to the gas mass in the gap region, where the binary torques are strongest. Thismass, always a small fraction of the disk mass, reaches the gap through the action ofinternal stresses in the accretion flow. Relative to the gas pressure, MHD stresses inthe disk body are about an order of magnitude larger than the “viscous” stresses esti-mated phenomenologically in previous work; in the gap, this ratio increases by another 35 –order of magnitude. As a result, the density of matter in the gap is also about an orderof magnitude larger than found in hydrodynamic calculations, and this contrast leadsdirectly to the larger torque.6. After time averaging, the accretion rate at the inner boundary is ∼
30% of the peak ratein the disk body. Compared with MM08 and normalized to the peak surface density,the accretion rate at the inner boundary in this simulation is ∼
40 times greater. Thiscontrast, like the contrast in the total torque, can be largely attributed to the strongerinternal stresses due to MHD effects. The time-dependent accretion rate is stronglymodulated on both the binary orbital frequency (by the stream) and the inner-diskorbital frequency (by the lump).7. Previous work on the angular momentum budget of the binary has focussed on thetorque it exerts on the disk. Because we find a substantially larger accretion rate thanprevious calculations, we also find that the binary’s gain in angular momentum due toaccretion can substantially offset its loss by torque. As a result, the estimated binarycontraction rate ˙ a/a ∼ − × − ( M d /M )( GM ) / a − / is only slightly larger than therate estimated by MM08. M d < M vs. M d > M These detailed results have a number of more general implications for the evolution ofcircumbinary disks, particularly in the context of supermassive black hole binaries. It isgenerally believed that stellar dynamical effects become relatively ineffective at driving theevolution of this sort of binary when its separation falls much below ∼ a/a .Defining the orbital shrinkage time t shrink ≡ | a/ ˙ a | , the time taken by torque alone to change 36 –the binary orbital angular momentum t torque ≡ M j bin /T ( ∞ ), and the disk accretion time t acc ≡ M d / ˙ M , we can rewrite this rate as t acc t shrink = 2 (cid:12)(cid:12)(cid:12)(cid:12) − t acc t torque + M d M (cid:18) j in j bin − (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) . (25)Using the values found in the simulation makes t acc /t shrink (cid:39) . M d /M , i.e., the ratio of theaccretion time to the orbital shrinkage time is roughly the same as the ratio of the disk massto the binary mass. Indeed, Lodato et al. (2009) argued that, unless the disk mass wereat least comparable to the binary’s secondary mass, it would be ineffective in driving binaryevolution. Although it is true that a mass comparable to the secondary’s mass must passthrough the inner region of the disk over an orbital evolution time t shrink , it is not necessarilytrue that this mass must be there all at once.Suppose, for example, that M d (cid:28) M . Because t acc (cid:28) t shrink in this case, if the diskmass is put in place once and for all, it would certainly be drained long before the binaryhas significantly evolved. On the other hand, one could also imagine situations in which thedisk is continuously replenished. In that case, the binary orbit could change substantiallydue to this interaction even though the instantaneous disk mass is always much smaller thanthe binary mass; it’s just that this process takes longer. Before determining how much timeis required, it is important to recall that the close subtraction between the torque term andthe accretion term may make the numerical value of ( t acc /t shrink ) / ( M d /M ) rather sensitive tospecific assumptions of our simulation, e.g., the binary mass ratio and the disk gas’s equationof state. If the mass accretion rate were reduced by a factor of two while the torque washeld fixed, the shrinkage rate would increase by a factor (cid:39)
4; conversely, if it were increasedrelative to a constant torque by a ratio (cid:38) .
2, the binary orbit would actually expand overtime.These processes could also be affected by the fact that the accretion rate onto thebinary is a substantial fraction of the accretion rate in the outer disk. The specific numberwe found ( (cid:39) z Ω ) wouldalter saturation of the magneto-rotational instability, but more rapid inflow is one plausibleconsequence.The opposite case is also worth considering, in which M d (cid:29) M , so that t acc (cid:29) t shrink . Ifthis is so, then even a one-time deposit of mass in the disk might drive strong evolution inthe binary orbit. However, when the binary separation decreases, the position of the strongtorques moves to smaller distance from the center of mass as well, raising the question ofhow much mass might be there to receive those torques. We can estimate that time withinthe disk body of our simulation via t inflow ( r ) ≡ (cid:82) r dr (cid:48) / (cid:104) v r (cid:105) ρ,t , where (cid:104) v r (cid:105) ρ,t represents thetime- and density-weighted shell-averaged inflow velocity. In a steady-state disk around apoint-mass, t acc would match t inflow at r = r p . However, to the degree that binary torquesretard accretion, t acc exceeds the inflow time measured at the surface density peak, which is t inflow ( r p ) ∼ × Ω − . At sufficiently large radius, t inflow ( r ) will nearly always be greaterthan t acc because t inflow ( r ) ∝ r / ( r/H ) ; only a rapid outward flare in the disk thicknesscould alter this conclusion. It has long been thought that if t inflow ( r p ) ∼ t acc (cid:29) t shrink , thedisk and the binary would decouple once the binary shrinks by a factor of a few because toolittle gas would be found close enough to the binary to feel the torques; if orbital compressiondepends solely upon interaction with a circumbinary disk, the shrinkage time could thereforenever be much shorter than t inflow ( r p ) (e.g., Gould & Rix 2000; Armitage & Natarajan 2002).However, the dynamical behavior seen both in our simulation and in others hints thatthis assumption might not be valid. The streams occurring in our simulation travel inwardon a timescale ∼ Ω − , which is (cid:28) t inflow . Although a smaller binary separation would makethe trajectories of the streams at radii not far inside the disk edge more nearly angularmomentum-conserving, and therefore make it more difficult for them to penetrate to muchsmaller radii, the internal magnetic stresses within the streams might be strong enough toreplace the binary torques. Indeed, we found that the nominal “ α ” in the gap region inour simulation was one or two orders of magnitude larger than in the disk body. There isalso another mechanism, commonly seen in global MHD simulations of disks around point-masses, that might contribute in this way. As emphasized by Guan, Hawley & Krolik (2011), 38 –it is easy for disks far from inflow equilibrium to drive mass-inflow rates much greater thanthose expected on the basis of inflow equilibrium estimates. For example, when an annulardisk with a seed poloidal magnetic field is created, orbital shear rapidly creates toroidal fieldfrom radial field. The − B r B φ stress acting on the low density matter near the disk’s inneredge rapidly removes its angular momentum, propelling it inward. As the binary separationbecomes small compared to r p , the gravitational potential comes to resemble a point-mass,and similar dynamics might be expected. Because the exertion of a sizable torque requiresonly a small fraction of the disk mass to orbit in the gap region, where the torque per unitmass is greatest, a strong torque might continue even though the bulk of the disk massremains far from the binary.A massive disk is also likely to be subject to self-gravity, a mechanism not treated in oursimulation. If M d /M (cid:38) ( c s /v orb ), self-gravitational instability can be triggered (Goldreich& Lynden-Bell 1965). The evolution of this instability depends on the ratio of the coolingtime to the dynamical timescale. If the ratio is greater than unity, the disk can maintainmarginal stability, and the self-gravity would add extra stress to facilitate the accretion(Gammie 2001; Lodato & Rice 2004, 2005). If the cooling time is short, which is likely tobe the case for disks around massive black holes, the disk is likely to fragment and muchof its mass may be transformed into stars (Shlosman & Begelman 1989; Gammie 2001;Nayakshin 2006), reducing the gas surface density (Lodato et al. 2009). The consequencesfor interaction with the binary remain somewhat unclear. In sufficiently massive disks, starsmight form with orbits taking them sufficiently close to the binary that the summed torquesfrom many individual star-binary encounters could continue to shrink the binary; in effect,the stellar loss-cone is repopulated by local star formation (Cuadra et al. 2009). In lessmassive disks, stars might form only in the disk body. If most of the disk mass is convertedinto stars, the gas mass available to fill the gap would be reduced, possibly leading to botha smaller torque and a smaller accretion rate. On the other hand, gravitational torques dueto non-axisymmetric fluctuations in the stellar density might create stresses strong enoughto overcome the reduction in gaseous disk mass. (Binney & Tremaine 2008). Given all thesecomplicated possibilities, it is clear that further study of self-gravitating circumbinary diskswill be necessary before any strong conclusions can be drawn on the fate of massive disks. The work done by the binary on the streams can heat the inner edge of the disk throughstream impact. The total energy per unit time delivered by the binary torque can be esti- 39 –mated as L torq ≡ Ω bin T ( ∞ ) (cid:39) . × − Ω a M d . (26)For disks around massive binary black holes, L torq ∼ . × M / a − / . N ergs s − , (27)where M = M/ M (cid:12) is the binary mass, a . = a/ . N = N H / . × cm − is the column density of the disk. Thus, the expected luminosity is wellbelow a typical AGN luminosity, and would therefore be difficult to detect unless its powerwere emitted in a small number of lines. Similarly, we can get the heating rate for disksaround stellar binaries: L torq ∼ . × M / a − / M d , − ergs s − , (28)where M = M/ M (cid:12) is the mass of the binary in solar masses, a = a/
10 AU is the binaryseparation in units of 10 AU, and M d , − = M d / . (cid:12) is the disk mass in units of 0 . (cid:12) .We note that, with our ideal MHD model, this luminosity only sets the upper limit of thepossible heating rate due to the binary torque. In reality, disks outside 10 AU might possesslow-ionization dead zones (Gammie 1996) which suppress the MHD turbulence and diminishthe mean effective α SS . Based on the accretion rate observed in disks around young stars,which is typically ∼ − M (cid:12) yr − (e.g., Hartigan et al. 1995; Gullbring et al. 1998) for TTauri stars, we expect the actual luminosity will be two to three orders of magnitude smallerthan estimated in equation (28). Because stream impact is periodic with the period of thebinary, the resulting luminosity might be modulated with a period (cid:39) M − / a / yrs. In the black hole binary context, circumbinary disk studies have been in part motivatedby a search for mechanisms to solve the “last parsec problem”, the expected slow-down inorbital evolution by stellar encounters when the binary separation becomes (cid:46)
APPENDIXA. Binary Eccentricity Evolution
In this section, we estimate the growth rate of the binary’s eccentricity based on asimple model of accretion and energy exchange between the disk and the binary. The orbitaleccentricity e bin can be written in terms of the total orbital energy E , angular momentum J , and total mass M : e bin = (cid:115) EJ µG M , (A1)where µ = m m / ( m + m ) is the binary’s reduced mass, m and m are the masses of thebinary components, and M = m + m . Taking the time derivative of equation (A1), wefind ˙ e bin : ˙ e bin e bin = e − e (cid:32) ˙ EE + 2 ˙ JJ − MM (cid:33) . (A2)Note that if the mass ratio m /m = 1 (as in our simulation), ˙ e bin is independent of anychange in m /m .In order to estimate the rate of change of the orbital energy E , we must understandhow the binary and circumbinary disk exchange energy. The binary can lose orbital energyby doing work on the surrounding disk through its time-averaged torque. In addition, theaccretion flow penetrating the gap ultimately becomes attached to the binary, bringing itsenergy to the binary. However, when the accreting matter begins to orbit around an individ-ual member of the binary, we must distinguish how much of its energy becomes associatedwith that motion as opposed to the binary orbit. In addition, there may also be radiationlosses. Both of these effects could be computed if the binary members were on the grid ofthe simulation, but they were not in the calculation we performed. Consequently, the bestwe can do here is to bound the range of possibilities.For the greatest physical insight, it is convenient to group the contributions accordingto whether they are associated with the binary-disk torque or with accretion, i.e., writing theterms in the parentheses of equation (A2) as ( ˙ E/E + 2 ˙
J /J ) torq + ( ˙ E/E + 2 ˙
J /J ) acc − M /M ,where the subscript ‘torq’ denotes torque related contributions, and ‘acc’ the accretion relatedquantities.Consider the binary-disk interaction first. When the binary’s orbit is exactly circular,the work done by the disk on the binary can be written in terms of the azimuthal component 42 –of the binary’s gravity acting on the disk: W = (cid:90) Σ( r, φ ) (cid:18) ∂ Φ r∂φ (cid:19) Ω bin rdrdφ = − T ( ∞ )Ω bin . (A3)In our simulation, this contribution (in terms of the change in binary orbital energy) is (cid:39) − . GM a Σ Ω bin (also see section 5.3). Because the orbital energy of the binary is − GM / (8 a ), the ratio W/E = 0 . a Σ /M )Ω bin . The rate at which the orbital angularmomentum changes due to torque on the disk is − . GM a Σ ; in ratio to the binaryangular momentum, (1 / M √ GM a , this contribution to ˙
J /J is − . a Σ /M )Ω bin . Thecombination ( ˙ E/E + 2 ˙
J /J ) torq for these two pieces is therefore zero to within the accuracywith which we can do the calculation.To describe the energy and angular momentum brought to the binary orbit by accretion,we employ the following toy model. Around each of the binary member, the accreting gaswould presumably form an interior disk: a disk extending out to a distance r d from itscentral mass. We imagine each binary component together with its disk as a ‘hard sphere’of that radius orbiting around the center of mass of the binary system. When the stream ofaccreting matter reaches the sphere’s surface, it sticks to the sphere, adding its momentumto the sphere’s. It also delivers its potential energy. For simplicity, we will deal only inquantities averaged over times much longer than the binary’s orbital period, but shorterthan the timescale of the binary’s orbital evolution.In this model, the rate of change of kinetic energy in the binary orbit due to accretionis ˙ E k = ˙ M(cid:126)v s · (cid:126)v bin , (A4)where E k is the binary’s kinetic energy, (cid:126)v s the velocity of the stream, and (cid:126)v bin the velocityof the sphere, i.e., the orbital velocity of the individual binary component. For an equalmass binary, the time averaged | (cid:126)v bin | = (1 / (cid:112) GM/a with respect to the center of mass ofthe binary. From energy conservation, we see that | (cid:126)v s | = (cid:112) (cid:15) (cid:48) in + GM/a + GM/r d , where (cid:15) (cid:48) in is the specific energy of the stream when it hits the ‘hard sphere’. This energy could bedifferent from (cid:15) in , the one measured when the stream crosses the simulation boundary. Weapproximate r d = 0 . a as the radius of the tidally truncated interior disk for an equal-massbinary (Paczy´nski 1977). Mass addition from the stream diminishes the binary’s potentialenergy by ˙ E p = − ˙ M ( GM/ a ). This can be viewed alternatively as the potential energy ofthe mass stream with respect to the other binary member or as twice the potential energyper unit mass of the binary as a whole because there is a contribution both from the arrivingmatter and from the deepening of the potential due to the increase in total mass associatedwith accretion. The potential energy of the stream with respect to the binary member it isjoining is associated entirely with the orbit of the stream around that member and does not 43 –contribute to the orbital energy of the binary. Summing the kinetic and potential energybrought to the binary by accretion and forming its ratio to the orbital energy, we find (cid:32) ˙ EE (cid:33) acc = 4 ˙ MM (cid:34) − (cid:104) cos θ (cid:105) (cid:115) (cid:15) (cid:48) in GM/a + 1 + ar d (cid:35) , (A5)where (cid:104) cos θ (cid:105) is the time averaged angle between (cid:126)v s and (cid:126)v bin . Assuming (cid:15) (cid:48) in (cid:39) (cid:15) in = − . GM/a (its value measured at the simulation inner boundary r = r in by taking massflux weighted shell averages of the kinetic and potential energy), the accretion contributionto the change in binary energy is (cid:32) ˙ EE (cid:33) acc = 4 ˙ MM (1 − . (cid:104) cos θ (cid:105) ) . (A6)Angular momentum is brought to the binary at the fractional rate (cid:32) ˙ JJ (cid:33) acc = 3 .
76 ˙ MM (A7)because the specific angular momentum of the binary is (1 / √ GM a and the specific angularmomentum of the stream when it crosses the simulation boundary is j in = 0 . √ GM a .The total rate of change of eccentricity (see equation (A2))is then the sum of the accre-tion contributions due to energy, angular momentum, and mass acquisition˙ e bin e bin = 1 − e e (6 . (cid:104) cos θ (cid:105) − .
52) ˙
MM . (A8)We have previously found that ˙
M /M = 0 . a Σ /M )Ω bin , so that˙ e bin e bin = 1 − e e (6 . (cid:104) cos θ (cid:105) − . × − (cid:18) M d M (cid:19) Ω bin . (A9)In other words, the eccentricity decreases unless (cid:104) cos θ (cid:105) is very nearly unity. Because oursimulation assumes a circular orbit, this result seems surprising because there is no obviousreason why (cid:104) cos θ (cid:105) cannot be significantly smaller than one. In fact, however, the orbitalmechanics require (cid:104) cos θ (cid:105) (cid:39)
1, as can be seen from the following argument.When measuring the specific orbital energy of the stream as it crosses the inner problemboundary, we find its azimuthal velocity is (cid:39) . (cid:112) GM/a , almost twice the magnitude ofits radial velocity ( (cid:39) − . (cid:112) GM/a ). Thus, the stream’s velocity is nearly in the azimuthaldirection. Moreover, this azimuthal speed is rather greater than the binary orbital speed 44 –0 . (cid:112) GM/a , and, as shown by the six panels of Figure 13, most matter crosses the innerboundary at an orbital phase slightly ahead of the nearest member of the binary. Therefore,when the stream crosses the boundary, it does so at a small angle and then rapidly catchesup with the other member of the binary, traveling at almost constant radius. In other words,when it strikes that disk, it does so with (cid:104) cos θ (cid:105) (cid:39) e bin /e bin , as well as the disk interac-tion contribution, almost exactly cancels. To within the accuracy of this simulation, theseeffects appear to have at most a weak effect on the binary eccentricity.There is, however, a resonant interaction between binaries and their surrounding disksthat has the potential to drive a linear instability in the binary eccentricity. The resonance in-volves forcing by a binary potential component of the form Φ m,l ( r, φ, t ) = φ m,l ( r ) cos ( mφ − Ω bin t ).It is strongest at the 1:3 resonance, where m = 2 and l = 1 (Artymowicz et al. 1991). Usingthe same formalism as previously, the rate of growth of eccentricity can be written as d ln edt = − ˙ E − Ω bin ˙ J e | E | (A10)when e (cid:28)
1. Using ˙ E = Ω p ˙ J with Ω p = Ω bin l/m = Ω bin / d ln edt = − ˙ J e µ Ω bin a , (A11)where µ is the reduced mass of the binary. To evaluate this eccentricity growth rate for theresonance, we determine the torque ˙ J = − T , as defined in equation 11 of Artymowicz &Lubow (1994). This torque is determined through the use of equations (12)-(14) and (21)in that paper. We obtain d ln edt = 49 π m m M a Σ M Ω b , (A12)where surface density Σ is evaluated near the resonance at r (cid:39) / a and for an equal massbinary m = m = M/
2. If the density varies near this radius, then it should be suitablyaveraged over the resonance width of order ∼ ( H/r ) / r . There is a caveat regarding thisresult, however: it is derived in the context of linear theory and assumes an axisymmetricdisk whose fluid follows circular orbits. It is uncertain how much the growth rate mightchange in a highly non-axisymmetric, eccentric disk like the one found in our simulation.Moreover, once the binary eccentricity begins to grow, new resonances appear whose effecttends to counteract eccentricity growth (Lubow & Artymowicz 1992). Thus, it is hard toevaluate the ultimate impact of this possible instability. 45 – B. Linear Equation for Eccentricity Distribution
We apply the linear equation for eccentricity evolution given by GO06 that can bewritten as i∂ r ( f ( r ) ∂ r E ( r, t )) + if ( r ) E ( r, t ) + J ( r ) s ( r ) E ( r, t )= J ( r ) ∂ t E ( r, t ) , (B1)where E ( r, t ) = e ( r, t ) exp ( i(cid:36) ( r, t )) is the complex eccentricity, for real eccentricity e andperiapse angle (cid:36) . E is related to the linear perturbations from the axisymmetric circular velocity. For thevelocity expressed cylindrical coordinates as ( u (cid:48) ( r, t ) exp ( − iφ ) , v (cid:48) ( r, t ) exp ( − iφ )) , we havethat u (cid:48) ( r, t ) = ir Ω( r ) E ( r, t ) (B2)and v (cid:48) ( r, t ) = 12 r Ω( r ) E ( r, t ) , (B3)where Ω( r ) is the Keplerian orbital frequency about the binary of mass M given byΩ( r ) = (cid:114) GMr . (B4)Quantity J ( r ) is the disk angular momentum per unit radius divided by π and is given by J ( r ) = 2 r Ω( r )Σ( r ) . (B5)Functions f ( r ) and f ( r ) are given by f ( r ) = γP ( r ) r , (B6) f ( r ) = dPdr r + J ( r ) ˙ (cid:36) g ( r ) , (B7)where P ( r ) is the two-dimensional (vertically integrated) disk pressure, Σ( r ) is the disksurface density. Quantity ˙ (cid:36) g is the gravitational precession rate of a free particle on aneccentric orbit about the equal mass binary, which is given by˙ (cid:36) g ( r ) = − r Ω ∂ r ( r ∂ r Φ ( r )) (B8)where Φ ( r ) = − GMπr χK ( χ ) (B9) 46 –with χ = 2 ra (0 . a + r ) (B10)and K is the complete elliptic integral of the first kind. The terms involving quantities f and f describe the eccentricity propagation and precession, respectively. Quantity γ is thegas adiabatic index. Real function s ( r ) is the eccentricity injection rate due to some sourcethat is assumed to be distributed as a gaussian of width w centered at radius r c and centralvalue s c / ( √ πw ) s ( r ) = s c exp ( − ( r − r c ) /w ) √ πw . (B11)Following along the lines G06 and L10, we adopt the inner boundary condition ∂ r E ( r i , t ) = 0 . (B12)For a Keplerian disk, the divergence of the velocity is proportional to ∂ r E . Consequently,this boundary condition is equivalent to requiring that the Lagrangian density perturbationnear the disk inner edge vanishes. The outer boundary condition is taken to be E ( r o , t ) = 0 . (B13)To ensure that the trapping of the eccentricity is not an artifact of the outer boundarycondition, we test that the resulting mode structure is independent of the outer radius r o .To obtain modes, we write the eccentricity as E ( r, t ) = E ( r ) exp ( iωt ) (B14)where ω is a complex eigenfrequency. The eccentricity distribution equation (B1) has in-finitely many modes and eigenfrequencies. C. Eccentricity Distribution Far From Disk Inner Edge
We determine here the analytic form of the eccentricity distribution far from the inneredge of the circumbinary disk, based on the linear theory of Appendix B. In that limit, weassume that the local precession rate is small compared with the eigenvalue ω defined inequation (B14). This assumption is expected to hold, since the gravitational precession rate˙ (cid:36) g declines rapidly with r . We also expect the pressure contribution to the precession rate todecline, since the pressure is expected to decline with r . We assume power law dependencesfor the gas sound speed and surface density c s ( r ) = c s x − b , (C1)Σ( r ) = Σ x − b , (C2) 47 –where x = r/a and c s , Σ , b , and b are constants. Equation (B1) then simplifies to x − b +3 / E (cid:48)(cid:48) ( x ) + ( − b − b + 3) x − b +1 / E (cid:48) ( x ) = ˜ ωE ( x ) , (C3)where we ignored the eccentricity source term that is assumed to be localized near the diskinner edge. Dimensionless frequency ˜ ω is a constant given by˜ ω = 2 γ (cid:18) Ω bin ac s (cid:19) ω Ω bin . (C4)For x (cid:29)
1, equation (C3) has an asymptotic solution in lowest order of the form E ( x ) ∼ exp ( − c x c ) x c , (C5)where c = 4 √ ˜ ω b , (C6) c = 14 (1 + 2 b ) , (C7) c = 18 (9 − b − b ) . (C8)In this relation we ignore the overall scale factor for the eccentricity distribution, since we areinterested in its function form. The appropriate root for √ ˜ ω is taken such that Re ( √ ˜ ω ) > γ = 1 , b = 0 (isothermal), c s = 0 . bin a ,and ω (cid:39) − . bin i . Quantity b = 1 for an isothermal constant α decretion disk, although b is closer to -1 in the outer simulated region during much of the simulation. In any case,quantity b has no influence on the dominant factor, the exponential. We adopt b = 1 andobtain e ( x ) = | E ( x ) | ∼ exp ( − . x / ) x / . (C9)We approximate the exponent of equation (C9) by means of a Taylor series about x = 3(typical of the simulated outer disk region) and obtain in leading order e ( x ) ∼ exp ( − . x ) x / (C10)and the length scale for the exponential decay is then ∼ . a . 48 – REFERENCES
Antonucci, R.R.J. 1993, ARA&A, 31, 473Armitage, P., & Natarajan, P. 2002, ApJ, 567, L9Armitage, P., & Natarajan, P. 2005, ApJ, 634, 921Artymowicz, P. et al. 1991, ApJ, 370, L35Artymowicz, P. & Lubow S.H. 1994, ApJ, 421, 651Artymowicz, P. & Lubow S.H. 1996, ApJ, 467, 77Balsara D.S. & Krolik, J.H. 1993, ApJ, 402, 109Bate, M.R. et al. 2003, MNRAS, 341, 213Baruteau, C. et al. 2011, arXiv:astro-ph/1107.2774Bate, M.R., & Bonnel, I.A. 1997, MNRAS, 285,33Begelman, M.C., Blandford, R.D., & Rees, M.J. 1980, Nature, 287, 307Berczik, P. et al. 2006, ApJ, 642, 21Binney, J., & Tremaine, S. 2008, Galactic Dynamics(2nd ed.;Princeton, NJ: Princeton Univ.Press)Brouwer,D., & Clemence, G.M.1961, Methods of Celestial Mechanics (New York: Academic)Bryden, G. et al. 1999, ApJ, 514, 344Corrales, Lia R. et al. 2010, MNRAS, 404, 947Cuadra, J. et al. 2009, MNRAS, 393, 1423D’Angelo, G., Lubow, S.H. & Bate, M.R. 2006, ApJ, 652, 1698de Val-Borro, M. et al. 2011, MNRAS, 423, 2679Dotti, M. et al. 2009, MNRAS, 396, 1640Dutrey, A., Guilloteau, S.,& Simon, M. 1994, A&A, 286, 149Escala, A. et al, 2004, ApJ, 607, 765 49 –Escala, A. et al, 2005, ApJ, 630, 152Farris, B.D. et al. 2011, Phys. Rev. D, 84, 024032Gammie, C.F. 1996, ApJ, 457, 355Gammie, C.F. 2001, ApJ, 553, 174Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125Goldreich, P. & Tremaine, S. 1980, ApJ, 241, 425Goldreich, P. & Tremaine, S. 1982, ARA&A, 20, 249Goodchild, S. & Ogilvie, G. 2006, MNRAS, 368, 1123Gould, A.,& Rix, H.-W. 2000, ApJ, 532, L29Greaves, J.S. et al. 2008, MNRAS, 391, L74Guan, X., Hawley, J.F. & Krolik, J.H. 2011, ApJ, 738, 84Gullbring, E. et al 1998, ApJ, 492, 323G¨unther, R. & Kley, W. 2002, A&A, 387, 550Hanawa, T., Ochi, Y., & Ando, K. 2010, ApJ, 708, 485Hartigan,P. et al. 1995, ApJ, 452, 736Hashimooto, J. et al. 2011, ApJ, 729, L17Hawley, J.F. & Krolik, J.H. 2001, ApJ, 548, 348Hawley, J.F. & Stone, J.M. 1995, CoPhC, 89, 127Hayasaki, K. et al. 2007, PASJ, 59, 427Hayasaki, K. et al. 2008, ApJ, 682, 1134Ivanov, P.B. et al. 1999, MNRAS, 307, 79Kalas, P. et al. 2008, Science, 322,1345Key, J.S. & Cornish, N.J. 2011, Phys. Rev. D, 83, 3001Khan, F. et al. 2011, ApJ, 732, 89 50 –Kley, W. & Dirksen, G. 2006, A&A, 447, 369Kley, W., Papaloizou, J.C.B.,& Ogilvie, G.I. 2008, A&A, 487, 671Kocsis, B., & Loeb, A. 2008, Phys. Rev. Lett., 101, 041101Krist, J.E., Stapelfeldt, K.R., & Watson, A.M., 2002, ApJ, 570, 785Krolik, J.H. 2007, ApJ, 661, 52Krolik, J.H. 2010, ApJ, 709, 774Krolik, J.H. & Begelman, M.C. 1986, ApJ, 308, L55Lin, D.N.C., & Papaloizou, J.C.B. 1986, ApJ, 309, 846Lin, D.N.C., & Papaloizou, J.C.B. 1993, in Protostars and Planets III, ed, E. Levy & J.Lunine (Tucson: Univ. Arizona Press), 749Lodato, G.,& Rice, W.K.M.2004, MNRAS, 351, 630Lodato, G.,& Rice, W.K.M.2005, MNRAS, 358, 1489Lodato, G. et al. 2009, ApJ, 398, 1392Lubow,S.H. 1991,ApJ, 381, 268Lubow, S.H. & Artymowicz, P. 1992, in Binaries as Tracers of Star Formation, A. Duquennoy& M. Mayor, eds., p.145Lubow, S.H. 1994, ApJ, 432,224Lubow, S.H. et al. 1999, ApJ, 526, 1001Lubow, S.H. & D’Angelo, G. 2006, ApJ, 641, 526Lubow, S.H. 2010, MNRAS, 406,2777MacFadyen, A.I. & Milosavljevi´c, M. 2008, ApJ, 672, 83Mathieu, R.D. 1994, ARA&A, 32,465Mayer, L. et al. 2007, Science, 316, 1874Merritt, D. & Milosavjevi´c, M. 2005, Living Rev. Relative., 8, 8Meyer-Vernet, N.& Sicardy, B., 1987, Icarus, 69, 157 51 –Nayakshin S. 2006, MNRAS, 372, 143Nelson, R.P. & Papaloizou, J.C.B. 2003, MNRAS, 339, 993Nixon, C.J., Cossins, P.J, King, A.R. & Pringle, J.E. 2011, MNRAS, 412, 1591Noble, S.C., Krolik, J.H., & Hawley, J.F. 2010, ApJ, 711, 959Osaki, Y. 1996, PASP, 108,39Paczy´nski, B. 1977, ApJ, 216, 822Papaloizou, J.C.B., Nelson, R.P.,& Masset, F. 2001, A&A,366,263Perets, H., & Alexander, T. 2008, ApJ, 677, 146Pier, E.A. & Krolik, J.H. 1992, ApJ, 399, L23Preto, M. et al. 2011, ApJ, 732, 26Pringle, J.E. 1991, MNRAS, 248, 754R¨odig, C. et al. 2011, arXiv:astro-ph/1104.3868RR´ozyczka, M., & Laughlin, G. 1997, in ASP Conf.Ser. 121, IAU Colloq. 163, AccretionPhenomena and Related Outflows, 792Sesana,A. et al. 2011,arXiv:astro-ph/1107.2927Schnittman, J.D., & Krolik, J.H.2008, ApJ, 684, 835Shakura, N.I. & Sunyaev, R.A. 1973, A&A, 24, 337Shi,J.,& Krolik, J.H. 2008,ApJ, 679, 1018Shlosman, I.,& Begelman, M.1989, ApJ, 341, 685Stone, J.M., Gammie, C.F., Hawley, J.F. & Balbus, S.A. 2000, in Protostars and PlanetsIV, ed. V. Mannings & Boss(Tuscon: Univ. Arizona Press), 589Stone, J.M. & Norman, M.L. 1992,ApJS, 80, 753Stone, J.M. & Norman, M.L. 1992, ApJS, 80, 791Tanaka, T. et al. 2011, arXiv:1107.2937Uribe, A. et al. 2011, ApJ, 736, 85 52 –Winters, W.F. et al. 2003,ApJ, 589, 543Yuan, C. & Cassen, P. 1994, ApJ, 437, 338
This preprint was prepared with the AAS L A TEX macros v5.2.
53 –Fig. 1.—
A series of snapshots of the disk surface density at different times. Density contours areon a linear scale. The color scale encoding the density (see the color bar for each panel) has twicethe range in the bottom two panels as in the top two; the white regions are density peaks whichare off the scale. White dots show the position of the binary; the faint white solid circle shows theboundary of the central cut-out; the white dash-dotted circles represent the radii r = 1, 2 and 3 a .
54 –Fig. 2.—
Disk mass profile history. The mass interior to r = 1 . a is given by the solid curve, 1 . a the dotted curve, 2 . a the dashed curve, 3 a the dash-dotted curve, and 4 a the dash-triple-dottedcurve.
55 –Fig. 3.—
Time- and shell-integrated surface density (top) and accretion rate (bottom) averagedover t = 250–350 (∆ T , solid) and t = 350–450 (∆ T , dashed). The dotted line in the surfacedensity plot displays the ∝ r − density profile of a ‘decretion’ disk.
56 –Fig. 4.—
Top: Time-averaged Reynolds stress (green), Maxwell stress (red) and the total stress(black) as functions of radius, solid for ∆ T and dashed for ∆ T . Bottom: Same stresses butnormalized with the local pressure.
57 –Fig. 5.—
Snapshots of Reynolds stress (left column) and Maxwell stress (right column) at t = 305averaged either vertically (top row) or azimuthally (bottom row). The z = ± H and ± H levelsare plotted as white dash-dotted lines.
58 –Fig. 6.—
Time averages (dashed for ∆ T , solid for ∆ T ) of local (top left), total(top right), andspecific (bottom left) binary torque as functions of radius. The linear theory prediction using 1/4the time-averaged surface density at the 3 : 2 resonance is shown as a red-dotted curve in theupper two panels. The history of the total torque is shown in the bottom right panel. The curveis smoothed by boxcar average whose width is twice the binary period.
59 –Fig. 7.—
Radial derivatives of the different sorts of time-averaged and shell-integrated angularmomentum flux ( dF J /dr ) appearing in equation (15) as functions of radius. The averaging period is t = 300–320. Net rate of change of the local angular momentum (solid black curve); binary torquedensity (blue dashed curve); Reynolds stress (green dash-dotted curve); Maxwell stress (red dash-dotted curve); the sum of Reynolds and Maxwell stresses (cyan dashed curve); angular momentumflux due to gas advection (black dashed curve; the inferred flux according to conservation law isshown with the magenta dotted curve for comparison).
60 –Fig. 8.—
Top: Space-time diagram of the local eccentricity (smoothed over T bin ); color contoursare logarithmic (see color bar). Bottom: an enlarged section of the space-time diagram (withouttime smoothing) to show short timescale behavior in the inner disk. White regions are eccentricitypeaks off the scale.
61 –Fig. 9.—
Time-averaged disk eccentricity as a function of radius at three epochs: t = 150–250(black solid curve), t = 250–350 (red dotted curve) and t = 350–450 (blue dashed curve). The graydash-dotted line shows a radial dependence ∝ exp [2 − . r/a )]. Table 1. Properties of Simulations of Circumbinary Accretion DisksLabel Type of Simulation Resolution a Radial Extent b B3D MHD 400 × ×
512 (0 . , × . , × . , × . , × . , × . , × . , a Resolution is listed as r × θ × φ in MHD simulation and r × φ in Hydro-dynamic simulations. b The azimuthal extent is (0 , π ) for all simulations.
62 –Fig. 10.—
Evolution of the disk eccentricity. The dash-dotted lines show linear fits for theexponential growth phase and the later saturation phase.
63 –Fig. 11.—
Phase angle of the disk eccentricity vector as a function of time.
Fig. 12.—
Space-time diagram of the m=1 Fourier component of the surface density. Colorcontours are logarithmic (see color bar).
64 –Fig. 13.—
Three pairs of snapshots of the surface density in the binary’s co-rotating frame inlogarithmic scale (white regions are off the scale). The density-weighted vertically averaged velocityis shown by black arrows. To show how phase-dependent stream effects change secularly over time,the time separation between the right and left panels in each pair is a fraction of a binary orbit,while the intervals from the first pair to the rest are many orbits. The Roman numerals show thequadrants referred in section 3.4.
65 –Fig. 14.—
Vertically-averaged horizontal magnetic field superposed on surface density contoursat t = 300 .
5. Color surface density contours are logarithmic (see color bar), white are peak valuesthat is off the scale.
66 –Fig. 15.—
Top: Accretion rate as a function of time (in units of the binary orbital period dividedby 2 π ) during the later stages of the simulation. Middle: Power spectral density of the accretionrate as a function of frequency in units of the binary orbital angular frequency Ω bin , or 2 π timesthe binary frequency. Bottom: Snapshots of radial mass flux (cid:104) ρv r (cid:105) z L z at the times shown by thearrows in the top panel. Color contours are in a linear scale (see color bar). Contour lines representthe surface density in two groups: black contours show low surface density (10 − < Σ < − . )on a logarithmic scale in order to highlight the streams; white contours show high surface density(1 . < Σ < .
67 –Fig. 16.—
Disk eccentricity as a function of radius from r = 2 a to r = 5 a , normalized by its valueat the inner edge. The heavy solid line is the eccentricity distribution obtained from the simulationswhen the eccentricity evolution is in the linear regime. The middle dashed line is the eccentricitydistribution determined by linear theory (equation (B1), with the input parameters appropriate forthe simulation (see parameter details in § . s ( r ) /s ( r c ). The upper dashed curve is the eccentricity distributionfor the fundamental free eccentric mode. The lower dashed curve shows the result of a highereccentricity growth rate, about twice as high as in the middle dashed curve.
68 –Fig. 17.—
The evolution of the disk eccentricity with varying cut-off sizes: r in = 0 . a (blacksolid curve), 1 . a (blue dash-dotted curve), 1 . a (cyan dash-dot-dotted curve), 1 . a (magenta longdashed curve), 1 . a (green dashed curve) and 1 . a (red dotted curve). There is a qualitative changewhen the cut-off becomes larger than (cid:39) . a
69 –Fig. 18.—
Rate of change of the disk eccentricity due to the stream impact. Solid line shows theinstantaneous rate; dashed line shows the time-averaged rate.
Fig. 19.—
Snapshots of the lump concentration and feeding mechanisms at two specific times, t = 443 (left) and tt