Modeling Magnetorotational Turbulence in Protoplanetary Disks with Dead Zones
aa r X i v : . [ a s t r o - ph . E P ] A ug A P J A
CCEPTED
Preprint typeset using L A TEX style emulateapj v. 5/2/11
MODELING MAGNETOROTATIONAL TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES S ATOSHI O KUZUMI
Department of Physics, Nagoya University, Nagoya, Aichi 464-8602, Japan; [email protected]
AND S HIGENOBU H IROSE
Institute for Research on Earth Evolution, JAMSTEC, Yokohama, Kanagawa 236-0001, Japan
ApJ Accepted
ABSTRACTTurbulence driven by magnetorotational instability (MRI) crucially affects the evolution of solid bodies inprotoplanetary disks. On the other hand, small dust particles stabilize MRI by capturing ionized gas particlesneeded for the coupling of the gas and magnetic fields. To provide an empirical basis for modeling the co-evolution of dust and MRI, we perform three-dimensional, ohmic-resistive MHD simulations of a verticallystratified shearing box with an MRI-inactive “dead zone” of various sizes and with a net vertical magnetic fluxof various strengths. We find that the vertical structure of turbulence is well characterized by the vertical mag-netic flux and three critical heights derived from the linear analysis of MRI in a stratified disk. In particular,the turbulent structure depends on the resistivity profile only through the critical heights and is insensitive tothe details of the resistivity profile. We discover scaling relations between the amplitudes of various turbulentquantities (velocity dispersion, density fluctuation, vertical diffusion coefficient, and outflow mass flux) andvertically integrated accretion stresses. We also obtain empirical formulae for the integrated accretion stressesas a function of the vertical magnetic flux and the critical heights. These empirical relations allow to predictthe vertical turbulent structure of a protoplanetary disk for a given strength of the magnetic flux and a givenresistivity profile.
Subject headings: dust, extinction — planets and satellites: formation — protoplanetary disks INTRODUCTIONPlanets are believed to form in protoplanetary gasdisks. The standard scenario for planet formation con-sists of the following steps. Initially, submicron-sized dustgrains grow into kilometer-sized planetesimals by collisionalsticking and/or gravitational instability (Safronov 1969;Goldreich & Ward 1973; Weidenschilling & Cuzzi 1993).Planetesimals undergo further growth toward Moon-sizedprotoplanets through mutual collision assisted by gravita-tional interaction (Wetherill & Stewart 1989). Accretion ofthe disk gas onto protoplanets leads to the formation of gasgiants (Mizuno 1980; Pollack et al. 1996), while terrestrialplanets form through the giant impacts of protoplanets afterthe gas disk disperses by viscous accretion onto the centralstar and other effects (Chambers & Wetherill 1998).Turbulence in protoplanetary disks plays a decisive roleon planet formation as well as on disk dispersal. Theimpact of turbulence is particularly strong on the forma-tion of planetesimals since the frictional coupling of gasand dust particles governs the process. Classically, plan-etesimal formation has been attributed to the collapse ofa dust sedimentary layer by self-gravity (Safronov 1969;Goldreich & Ward 1973) and/or the collisional growth ofdust grains (Weidenschilling & Cuzzi 1993). The pres-ence of strong turbulence is preferable for dust growthwhen the dust particles is so small that Coulomb repulsionis effective (Okuzumi 2009; Okuzumi et al. 2011). How-ever, strong turbulence acts against the growth of macro-scopic dust aggregates since it makes their collision disrup-tive (Weidenschilling 1984; Johansen et al. 2008). Further-more, turbulence causes the diffusion of a dust sedimen-tary layer, making planetesimal formation via gravitationalinstability difficult as well (Weidenschilling 1984). Turbu- lence is also known to concentrate dust particles of particularsizes, but its relevance to planetesimal formation via grav-itational instability is still under debate (Cuzzi et al. 2001,2008, 2010; Pan et al. 2011). More recently, it has been sug-gested that two-fluid instability of gas and dust can producedust clumps with density high enough for gravitational col-lapse, but successful dust coagulation to macroscopic sizesseems to be still required for this mechanism to becomeviable (Youdin & Goodman 2005; Johansen & Youdin 2007;Johansen et al. 2007; Bai & Stone 2010). Besides, turbulencealso affects planetesimal growth as turbulent density fluctua-tions gravitationally interact with planetesimals and can raisetheir random velocities above the escape velocity (Ida et al.2008; Nelson & Gressel 2010). The fluctuating gravitationalfield can even cause random orbital migration of protoplanets(Laughlin et al. 2004; Nelson & Papaloizou 2004). Thus, tounderstand the growth of solid bodies in various stages, it isessential to know the strength and spatial distribution of diskturbulence.Interestingly, the evolution of solid bodies is not only af-fected but also affects disk turbulence. The most viable mech-anism for generating disk turbulence is the magnetorotationalinstability (MRI; Balbus & Hawley 1991). This instabilityhas its origin in the interaction between the gas disk and mag-netic fields, and therefore requires a sufficiently high ioniza-tion degree to operate. Importantly, whether the MRI oper-ates or not in each location of the disk is strongly depen-dent on the amount of small dust grains because they ef-ficiently capture ionized gas particles and thus reduce theionization degree (Sano et al. 2000; Ilgner & Nelson 2006;Okuzumi 2009). This implies that dust and MRI-driven tur-bulence affect each other and thus evolve simultaneously.The purpose of this study is to present an empirical ba- OKUZUMI & HIROSEsis for studying the coevolution of solid particles and MRI-driven turbulence. It is computationally intensive to simulatethe evolution of dust and MRI-driven turbulence simultane-ously, since the evolutionary timescale of solid bodies is gen-erally much longer than the dynamical timescale of the tur-bulence. For example, turbulent eddies grow and decay on atimescale of one orbital period (e.g., Fromang & Papaloizou2006), while dust particles grow to macroscopic sizes andsettle to the midplane spending 100–1000 orbital periods(e.g., Nakagawa et al. 1981; Dullemond & Dominik 2005;Brauer et al. 2008). However, this also means that MRI-driven turbulence can be regarded as quasi-steady in each evo-lutionary stage of dust evolution. Motivated by this fact, werestrict ourselves to time-independent ohmic resistivity, butinstead focus on how the quasi-steady structure of turbulencedepends on the vertical profile of the resistivity.To characterize the vertical structure of MRI-driven tur-bulence, we perform a number of three-dimensional MHDsimulations of local stratified disks including resistivity andnonzero net vertical magnetic flux. Inclusion of a nonzeronet vertical flux is important as it determines the satura-tion level of turbulence (Hawley et al. 1995; Sano et al. 2004;Suzuki & Inutsuka 2009, see also our Section 4). Simi-lar simulations have been done in a number of previousstudies (e.g., Miller & Stone 2000; Suzuki & Inutsuka 2009;Oishi & Mac Low 2009; Suzuki et al. 2010; Turner et al.2010; Gressel et al. 2011; Simon et al. 2011; Hirose & Turner2011). One important difference between our study and pre-vious ones is that we focus on general dependence of the sat-urated turbulent state on the model parameters such as the re-sistivity and net magnetic vertical flux.Our modeling of MRI-driven turbulence follows two steps.In the first step, we seek scaling laws giving the rela-tions among turbulent quantities. We express the relationsas a function of the vertically integrated accretion stress,which is the quantity that determines the rate at whichturbulent energy is extracted from the differential rotation(Balbus & Papaloizou 1999). As we will see, excellent scal-ing relations are obtained if we divide the integrated stressinto two components that characterize the contributions fromdifferent regions in the stratified disk (which we will call the“ disk core ” and “ atmosphere ”) In the second step, we findout empirical formulae that predict the vertically integratedstresses as a function of the resistivity profile and verticalmagnetic flux.The plan of this paper is as follows. In Section 2, we de-scribe the method and setup used in our MHD simulations.In Section 3, we introduce “critical heights” derived from thelinear analysis of MRI in stratified disks. As we will see later,these critical heights are useful to characterize the turbulentstructure observed in our simulations. We present our sim-ulation results in Section 4, and obtain scaling relations andpredictor functions for the quasi-steady state of turbulence inSection 5. In Section 6, we simulate dust settling in a deadzone to model the diffusion coefficient for small particles asa function of height. Effects of numerical resolutions on oursimulation results are discussed in Section 7. Our findings aresummarized in Section 8. SIMULATION SETUP AND METHODIn this section, we describe the setup and method adoptedin our stratified resistive MHD simulations.2.1.
Setup
Our MHD simulations adopt the shearing box approxima-tion (Hawley et al. 1995). We consider a small patch of diskcentered on the midplane of an arbitrary distance from thecentral star, and model it as a stratified shearing box corotat-ing with the angular speed Ω at the domain center. We use theCartesian coordinate system ( x , y , z ), where x , y , and z standfor the radial, azimuthal, and vertical distance from the do-main center. In addition, we assume that the gas is isothermalthroughout the box; thus, the sound velocity c s of the gas isconstant in both time and space.2.1.1. Initial Conditions
For the initial condition, we assume that the gas disk is ini-tially in hydrostatic equilibrium and is threaded by uniformvertical magnetic field B z . The assumption of the hydrostaticequilibrium leas to the initial gas density profile ρ = ρ exp (cid:18) - z h (cid:19) , (1)where ρ is the initial gas density at the midplane and h = c s Ω (2)is the pressure scale height. The ratio between the initialmidplane gas pressure P = ρ c s and the initial magnetic pres-sure B z / π defines the initial plasma beta β z ≡ πρ c s B z . (3)In this paper, the strength of the initial magnetic flux will bereferred by β - z rather than B z .2.1.2. Resistivity Profile
The main purpose of this study is to see how the turbulencedepends on the vertical profile of the ohmic resistivity. Weadopt a simple analytic resistivity profile based on the follow-ing consideration. For fixed temperature, the resistivity is in-versely proportional to the ionization degree (Blaes & Balbus1994). In protoplanetary disks, the ionization degree at eachlocation is determined by the balance of ionization (by, e.g.,cosmic rays and X-rays) and recombination (in the gas phaseand on grain surfaces). Detailed structure of the resistivityprofile depends on what processes dominate the ionizationand recombination. However, a general tendency is that theionization degree decreases toward the midplane of the disk,because the ionization rate is lower as the column depth isgreater and because the recombination rate is higher as the gasdensity is higher (see, e.g., Sano et al. 2000). Based on thisfact, we give the resistivity profile η ( z ) such that η increasesas z decreases. To be more specific, we adopt the followingresistivity profile η = η mid exp (cid:18) - z h η (cid:19) , (4)where η mid is the resistivity at the midplane and h η is the scaleheight of η .Equation (4) satisfies the important property of realistic re-sistivity profiles mentioned above. Furthermore, as shown in Note that the “gas scale height” is often defined as H = √ h = √ c s / Ω in the literature on stratified MHD simulations. ODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES 3
Table 1
Model Parameters and Initial Critical HeightsModel η mid c s h h h η β z log( Λ ) h ideal , h h Λ , h h res , h Ideal 0 · · · × ∞ . . . .
02 0 . × - . . . . .
02 1 . × - . . . . .
02 4 . × - . . . . .
02 20 . × - . . . . . . × - . . . . . . × - . . . . . . × - . . . . . . × - . . . . . . × - . . . . . . × - . . . . . . × - . . . . .
02 1 . × - . . . . .
02 1 . × - . . . . .
02 1 . × - . . . . .
02 1 . × - . . . . . · · · × - . . . . Appendix, Equation (4) exactly reproduces the vertical resis-tivity profile of a disk in some limited cases. However, it willbe useful to examine possible influences of limiting η to Equa-tion (4). In order to do that, we also consider a resistive profileused in Fleming & Stone (2003), η FS03 = η exp (cid:18) - z h (cid:19) exp Σ Σ CR √ π Z ∞| z | / h e - z ′ / dz ′ ! , (5)which is characterized by two parameters η and Σ / Σ CR (seeEquation (10) of Fleming & Stone 2003). Physically, Equa-tion (5) corresponds to the resistivity profile when the ioniza-tion degree is determined by the balance between cosmic-rayionization and gas-phase recombination (see also Appendix).We construct 17 simulation models using Equations (4) and(5). 16 models are constructed from Equation (4) with var-ious sets of the parameters ( η mid , h η , β z ). Table 1 lists theparameters adopted in each run. The model “Ideal” assumeszero resistivity throughout the simulation box. Models X0–X3 are defined by the same value of η mid but different valuesof h η . The difference between X, Y, and W models is in thevalue of η mid . The initial plasma beta is taken to be 3 × in all models except X1a–X1d. In addition, we construct onemodel from Equation (5) with η mid = η exp( Σ / Σ CR ) = 0 . Σ / Σ CR = 34 .
7. This set of parameters corresponds tothe “larger dead zone” model of Fleming & Stone (2003) de-fined by Re M , mid = 100 and Re M , z =2 √ h = 5 . × , whereRe M ≡ c s h /η is the magnetic Reynolds number with the typi-cal length and velocity set to be h and c s , respectively. We re-fer to the model with these parameters as model FS03L. Theinitial plasma beta in the FS03L model is taken to be 3 × .Note that our FS03L run is not exactly the same as the largerdead zone run of Fleming & Stone (2003) since they assumedzero net vertical magnetic flux.2.2. Method
We solve the equations of the isothermal resistive MHD us-ing the ZEUS code (Stone & Norman 1992). The domain isa box of size 2 √ h × √ h × √ h along the radial, az-imuthal, and vertical directions, divided into 40 × ×
200 grid cells. For all runs, the radial and azimuthal boundary con-ditions are taken to be shearing-periodic and periodic, respec-tively. Therefore, the vertical magnetic flux is a conservedquantity of our simulations. The vertical boundary conditionfor all runs except X1d is the standard ZEUS outflow condi-tion, where the fields on the boundaries are computed fromthe extrapolated electromotive forces; only for run X1d, weassume for numerical stability that the magnetic fields are ver-tical on the top and bottom boundaries as done by Flaig et al.(2010). We have checked that the results hardly depend on thechoice of the two types of boundary conditions. Note that theradial and azimuthal components of the mean magnetic fieldsare not conserved due to the outflow boundary condition. Wealso added a small artificial resistivity near the boundaries fornumerical stability (Hirose et al. 2009). A density floor of10 - ρ is applied to prevent high Alfvén speeds from haltingthe calculations. CHARACTERISTIC WAVELENGTHS AND CRITICALHEIGHTSAs we will see in the following sections, it is useful toanalyze the simulation results using the knowledge obtainedfrom the linear analysis of MRI. In this subsection, we intro-duce several quantities that characterize the linear evolutionof MRI. 3.1.
Characteristic Wavelengths of MRI
According to the local linear analysis of MRI includingohmic resistivity (Sano & Miyama 1999), the wavelength ofthe most unstable MRI mode can be approximately expressedas λ local ≈ max { λ ideal , λ res } , (6)where λ ideal ≡ π v Az Ω (7)and λ res ≡ π η v Az (8)are the characteristic wavelengths of MRI modes in the idealand resistive MHD limits, respectively, and v Az = B z / √ πρ isthe vertical component of the Alfvén velocity. Equation (6)can be written as λ local ≈ λ ideal max { , Λ - } , where Λ is theElsasser number defined by Λ ≡ v Az η Ω . (9)The Elsasser number determines the growth rate of the MRI.If Λ >
1, ohmic diffusion does not affect the most unstablemode lying at λ = λ ideal , and the local instability occurs rapidlyat a wavelength λ local ≈ λ ideal and at a rate ≈ Ω . If Λ < λ local ≈ λ res = Λ - λ ideal and at a slower rate ≈ Λ - Ω . We will refer to theformer case as the “ideal MRI,” and to the latter case as the“resistive MRI.”Figure 1 schematically illustrates how λ ideal and λ res varywith height | z | . In general, λ ideal grows toward higher | z | be-cause the Alfvén speed v Az increases as the density decreases.By contrast, λ res grows toward lower | z | because λ res is in-versely proportional to v Az and because η increases with de-creasing | z | (see the discussion in Section 2.1.2). OKUZUMI & HIROSE Figure 1.
Schematic illustration showing the vertical structure of a stratifiedprotoplanetary disk with vertical magnetic fields. The horizontal axis showsthe distance z from the midplane, while the vertical axis shows the charac-teristic wavelengths λ ideal (solid curve) and λ res (dashed curve) of MRI ateach z as well as the gas scale height h (dotted line). The set of three ra-tios λ ideal / h , λ res / h , and λ ideal / λ res ≡ Λ defines four layers. At | z | > h ideal ( Λ > λ ideal > h ), MRI is stabilized due to the weak gas pressure com-pared to the magnetic tension. At h Λ < | z | < h ideal ( Λ > λ ideal < h ;dark gray region), MRI operates without affected by ohmic dissipation. At h res < | z | < h Λ ( Λ < λ res < h ; light gray region), ohmic dissipation iseffective and MRI operates only weakly. At | z | < h res ( Λ < λ res > h ),ohmic dissipation perfectly stabilizes MRI. We refer to the regions | z | < h ideal and | z | > h ideal as the “disk core” and “atmosphere,” respectively. The global instability of a stratified disk can be describedin terms of the local analysis. As shown by Sano & Miyama(1999), the gas motion at height z is unstable if the local un-stable wavelength λ local is shorter than the scale height of thedisk, i.e., max { λ ideal , λ res } . h . (10)3.2. Critical Heights
With the global instability criterion (Equation (10)) togetherwith the vertical dependence of λ ideal and λ res , we can definethree different critical heights for a stratified disk.1. The first one is h ideal defined by λ ideal ( z = h ideal ) = h , (11)or equivalently, β z ( z = h ideal ) = 8 π , where β z ( z ) =8 πρ ( z ) c s / B z ( z ). At | z | & h ideal ( β z . π ), MRI does notoperate because the wavelengths of the unstable modesexceed the disk thickness ∼ h (Sano & Miyama 1999).We refer to the region | z | > h ideal as the “atmosphere”and to the region | z | h ideal as the “disk core.”2. The second one is h Λ defined by Λ ( z = h Λ ) = 1 , (12)or equivalently, λ ideal ( z = h Λ ) = λ res ( z = h Λ ). The layer h Λ | z | h ideal is the so-called “active layer,” whereMRI operates without affected by ohmic diffusion norgas stratification. The region | z | h Λ is what we callthe “dead zone,” where ohmic diffusion stabilizes themost unstable ideal MRI mode. For convenience, we Figure 2.
Vertical profiles of the ohmic resistivity η (upper panel) and theinitial Elsasser number Λ t =0 (lower panel) for models X0 (blue dot-dashedcurve), X1 (blue dotted curve), X2 (blue dashed curve), X3 (blue solid curve),Y3 (red curve), W3 (green curve), and FS03L (black curve). regard h Λ as zero when a dead zone is absent. This isthe case for model Ideal.3. The third one is h res defined by λ res ( z = h res ) = h . (13)Ohmic diffusion allows the resistive MRI to operate at h res . | z | . h Λ . At | z | . h res , ohmic diffusion stabilizesall the unstable MRI modes. Note that some previousstudies (e.g., Gammie 1996; Sano et al. 2000) used theterminology “dead zone” for the region | z | h res ratherthan | z | h Λ . In fact, as we will see in Section 4, theset of h Λ and h res best characterizes our dead zone. Weregard h res as zero when λ res is less than h at all heights.This is the case for Y models.The critical heights in the initial state ( h ideal , , h Λ , , and h res , ) are shown in Table 1 for all of our 17 simulations. Us-ing Equations (1) and (4), one can analytically calculate theinitial critical heights for models except FL03L as h ideal , = (cid:20) (cid:18) β z π (cid:19)(cid:21) / h , (14) h Λ , = (cid:18) Λ - + ( h / h η ) (cid:19) / h , (15) h res , = " (cid:0) π β - z Λ - (cid:1) + h / h η ) / h , (16)where Λ = 2 c s h η mid β z (17)is the initial Elsasser number at the midplane.ODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES 5 Figure 3.
Horizontally averaged Maxwell stress w M (upper panel) anddensity-weighted velocity dispersion ρδ v (bottom panel) as a function oftime t and height z for model X1. The solid, dotted, and dashed lines showthe critical heights z = h ideal , h Λ , and h res , respectively. Figure 2 shows the vertical profiles of the resistivity η andthe initial Elsasser number Λ t =0 for some of our models. Theinitial midplane Elsasser number Λ and initial critical heights( h ideal , , h Λ , , h res , ) are listed in Table 1 for all models. As onecan see from Table 1 and the lower panel of Figure 2, modelslabeled by the same number are arranged so that they havesimilar values of h Λ , .For turbulent states, we evaluate v Az in the Elsasser numberand the characteristic wavelengths as ( B z / πρ ) / , where theoverbars denote the horizontal averages. SIMULATION RESULTS4.1.
The Fiducial Model
We select model X1 as the fiducial model to describe indetail. Figure 3 shows how MRI-driven turbulence reachesa quasi-steady state in run X1. The upper and lower pan-els plot the horizontal averages of the Maxwell stress w M = - δ B x δ B y / π and the density-weighted velocity dispersion ρδ v , respectively, as a function of time t and height z . Thesolid, dotted, and dashed lines are the loci of the criticalheights h ideal , h Λ , and h res , respectively. As seen in the fig-ure, a quasi-steady state is reached within the first 40 orbits.The critical height h ideal measured in the quasi-steady state isslightly lower than that in the initial nonturbulent state. This isbecause the ideal MRI wavelength λ ideal ∝ v Az is increased bythe fluctuation in the vertical magnetic field, δ B z . By contrast, h Λ and h res are almost unchanged, because the fluctuation ofthe magnetic field is suppressed in the dead zone.Figure 4 shows the vertical structure of the disk averagedover a time interval 175 < Ω t / (2 π ) < h· · ·i denote the averages over time and horizontaldirections.In Figure 4(a), we compare the averaged gas density h ρ i with the initial density given by Equation (1). We see thatthe density is almost unchanged in the disk core ( | z | < h ideal )but is considerably increased in the atmosphere ( | z | > h ideal ). This is because the magnetic pressure is negligibly small inthe disk core but dominates over the gas pressure in the atmo-sphere. Figure 4(a) also shows the amplitude of the densityfluctuation, h δρ i / . As one can see, the density fluctuationis small ( h δρ i / ≪ h ρ i ) except at | z | ≫ h ideal .The magnetic activity in the disk can be seen in Fig-ure 4(b), where the vertical profiles of the magnetic energies( h δ B i / π , h B y i / π , and h B z i / π ) and Maxwell stress h w M i are plotted. One can see that these quantities peak near theouter boundaries of the active layers, | z | ≈ h ideal . This is be-cause the largest channel flows develop at locations where λ ideal ≈ h (see, e.g., Suzuki & Inutsuka 2009). In the deadzone, ohmic dissipation suppresses the fluctuation in the mag-netic fields, h δ B i , leaving the initial vertical field ( h B z i ≈ B z )and coherent toroidal fields ( h B y i ≈ h B y i ) generated by thedifferential rotation. Figure 4(c) shows the density-weighted velocity disper-sions h ρ ih δ v i and h ρ ih δ v z i and the Reynolds stress h w R i = h ρδ v x δ v y i . These quantities characterize the kinetic energy inthe random motion of the gas. Comparing Figures 4(b) and(c), we find that the drop in these quantities in the disk coreis not as significant as the drop in h δ B i and h w M i . This isan indication that sound waves generated in the active layerspenetrate deep inside the dead zone (Fleming & Stone 2003).Furthermore, we find that h ρ ih δ v i is approximately constant,i.e., the velocity dispersion h δ v i is inversely proportional tothe mean density h ρ i , in the disk core. This means that thekinetic energy density of fluctuation is nearly constant in thedisk core. This is another indication of sound waves, becausethe amplitude of the velocity fluctuation δ v is generally pro-portional to 1 / √ ρ for freely propagating sound waves. Theroot-mean-squared random velocity h δ v i / is shown in Fig-ure 4(d). The random velocity is subsonic in the disk core andexceeds the sound speed only in the atmosphere.An indication of freely sound waves can be also found inthe density fluctuation. Shown in Figure 4(e) is the meansquared density fluctuation h δρ i divided by the mean den-sity h ρ i . Since h δρ i / ≪ h ρ i , the quantity h δρ i / h ρ i ≈h δρ /ρ i is approximately proportional to the thermal energydensity of fluctuation, h c s δρ / ρ i . In the disk core, we seethat h δρ i / h ρ i is roughly constant along the vertical direc-tion, meaning that the amplitude of the density fluctuation, h δρ i / , is proportional to the square root of the mean den-sity h ρ i . The similarity between h ρ ih δ v i and h δρ i / h ρ i ispeculiar to sound waves, for which δ v / c s ∼ δρ/ρ .Figure 4(f) displays the profile of the vertical mass flux h ρ v z i . In the atmosphere, the vertical flux is outward, i.e., h ρ v z i > z > h ideal and h ρ v z i < z < - h ideal . This outflowresults from the breakup of large channel flows at the outerboundaries of the active layers, | z | ≈ h ideal (Suzuki & Inutsuka2009; Suzuki et al. 2010). The vertical mass flux reachesa constant value at | z | & h . This fact allows us to mea-sure a well-defined outflow flux for each simulation (see Sec-tion 5.1.3). 4.2. Model Comparison In our simulations, ohmic resistivity is not high enough to remove shear-generated, coherent toroidal fields. In this sense, our dead zone is an “undeadzone” in the terminology of Turner & Sano (2008). In the disk core ( | z | . h ideal ), h ρ ih δ v i is approximately equal to h ρδ v i since the density fluctuation is small (see Figure 4(a)). OKUZUMI & HIROSE
Figure 4.
Vertical profiles of various quantities averaged over x – y planes and over a time interval 175 orbits < t <
350 orbits for model X1. (a) Initial (dashedcurve) and time-averaged (solid curve) gas densities, and amplitude of the density fluctuation (dotted curve). (b) Magnetic energies h δ B i / π (solid curve), h B y i / π (dot-dashed curve), and h B z i / π (dashed curve), and Maxwell stress h w M i (dotted curve) normalized by the initial midplane gas pressure P = ρ c s . (c)Density-weighted velocity dispersions h ρ ih δ v i (thick solid curve) and h ρ ih δ v z i (dashed curve), and Reynolds stress h w R i (dotted curve). (d) Root-mean-squaredrandom velocity h δ v i / . (e) h δρ i / h ρ i , which represents the thermal energy of fluctuation. (f) Vertical mass flux h ρ v z i , plotted for positive (solid curve) andnegative (dashed curve) values. The regions shaded in dark and light gray indicate where ideal and resistive MRIs operate, respectively (see also Figure 1). Figure 5.
Vertical profiles of temporally and horizontally averaged Maxwell stress h w M i (red curves) and density-weighted velocity dispersion h ρ ih δ v i (blackcurves) normalized by the initial midplane gas pressure P = ρ c s for β z = 3 × models. The dark gray bars indicate where MRI operates without affected byohmic resistivity ( h Λ < | z | < h ideal ), while the light gray bars show where MRI operates but is weakened by ohmic resistivity ( h res < | z | < h Λ ). ODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES 7We now investigate how the vertical structure of turbulencedepends on the resistivity profile and vertical magnetic flux.Figure 5 displays the temporal and horizontal averages ofthe Maxwell stress h w M i and the density-weighted velocitydispersion h ρ ih δ v i as a function of z for various β z = 3 × models. As in Figures 1 and 4, the dark and light gray bars ineach panel indicate the heights where ideal and resistive MRIsoperate, respectively.The effect of changing the size of the dead zone can be seenin Figure 5(a), where models X1, X2, and X3 are compared.These models are characterized by the same values of β z and Λ but different values of h η . For all the models, h w M i sharplyfalls at z ∼ h res , meaning that h res well predicts where the resis-tivity shuts off the magnetic activity. By contrast, h ρ ih δ v i ex-hibits a flat profile at | z | . h ideal with no distinct change across z = h res nor z = h Λ . The only clear difference is the value of h ρ ih δ v i in the disk core, i.e., the value is lower when the deadzone is wider. Note that h ρ ih δ v i decreases more slowly thanthe column density of the active layers ( h Λ < | z | < h ideal ). Forexample, the active column density in model X1 is 20 timessmaller than that in model X3. However, the midplane valueof h ρ ih δ v i in the former is only five times smaller than thatin the latter. This suggests that even a very thin active layercan provide a large velocity dispersion near the midplane.Interestingly, the vertical structure of turbulence depends onthe critical heights ( h ideal , h Λ , and h res ) but are very insensitiveto the details of the resistivity profile. This can be seen in Fig-ure 5(b), where we compare runs with similar critical heights(runs X3, W3, and FS03L). We see that these models pro-duce very similar vertical profiles of h w M i and h ρ ih δ v i eventhough they assume quite different resistivity profiles (see theupper panel of Figure 2). This suggests that the vertical struc-ture of turbulence is determined by the values of the criticalheights.Importance of distinguishing h Λ and h res is illustrated inFigures 5(c) and (d). These panels compare five models (Y1–4 and Ideal) in which the resistive MRI is active at the mid-plane, i.e., h res = 0. Figure 5(c) shows models with h Λ > . h .We see that the profiles of h ρ ih δ v i and h w M i are very similarfor the three models. This implies that the vertical structure isdetermined by the value of h res when h Λ & . h . Figure 5(d)shows what happens when h Λ falls below 0 . h . Model Ideal isclearly different from the other models. Model Y4 ( h Λ = 0 . both features. In the lowerhalf of the disk ( z < z > h w M i is closer to that in model Ideal.Next, we see how the saturation level of turbulence dependson the vertical magnetic flux. Figure 6 shows the vertical pro-files of h w M i and h ρ ih δ v i for three runs with different valuesof β z (X1b, X1, and X1d). We see that these values increasewith decreasing β z . The peak value of h w M i is approximately10 - P , 10 - P , and 10 - P for runs X1b, X1, and X1d, re-spectively ( P = ρ c s is the initial midplane gas pressure). Thisindicates a linear scaling between the turbulent stress and β - z . SCALING RELATIONS AND PREDICTORFUNCTIONSNow we seek how the amplitudes of turbulent quantities de-pend on the vertical magnetic flux and the resistivity profile.We do this in two steps. First, we derive relations between theamplitudes of turbulent quantities and the vertically integratedturbulent stress. We then obtain empirical formulae that pre-
Figure 6.
Vertical profiles of temporally and horizontally averaged Maxwellstress h w M i (red curves) and density-weighted velocity dispersion h ρ ih δ v i (black curves) for models with different β z . The dotted, solid, and dashedcurves correspond to models X1b, X1, and X1d, respectively. dict the integrated stress as a function of the vertical magneticflux and the resistivity profile.5.1. Scaling Relations between Turbulent Quantities andVertically Integrated Accretion Stresses
The ultimate source of the energy of turbulence is theshear motion of the background flow. The accretion stress w xy ≡ w R + w M determines the rate at which the free energyis extracted. Therefore, we expect that the accretion stress isrelated to the amplitudes of turbulent quantities, such as thegas velocity dispersion and outflow mass flux.To quantify the rate of the energy input in the simulationbox, we introduce the effective α parameter α ≡ R h w xy i dz Σ c s , (18)where Σ = R h ρ i dz is the gas surface density. In the classical,one-dimensional viscous disk theory (Lynden-Bell & Pringle1974), the parameter α is related to the turbulent viscosity ν trub as ν trub = (3 / α c s / Ω , where the prefactor 3 / α also charac-terizes the vertically integrated mass accretion rate.As we will see below, it is useful to decompose α as α = α core + α atm , where α core ≡ R | z | < h ideal h w xy i dz Σ c s (19)and α atm ≡ R | z | > h ideal h w xy i dz Σ c s (20)are the contributions from the disk core ( | z | < h ideal ) and atmo-sphere ( | z | > h ideal ), respectively. Table 2 shows the values of α , α core , and α atm as well as the time-averaged critical heights( h ideal , h Λ , h res ) for all our simulations.5.1.1. Velocity Dispersion
Random motion of the gas crucially affects the growthof dust particles as it enhances the collision velocity be-tween the particles via friction forces (Weidenschilling 1984;Johansen et al. 2008). Here, we seek how the velocity disper-sion h δ v i is related to the integrated accretion stress.First, we focus on the velocity dispersion at the midplane, h δ v i mid . Figure 7(a) shows h δ v i mid versus the total accretion OKUZUMI & HIROSE Table 2
Time-averaged Properties of MHD SimulationsModel h ideal h h Λ h h res h α - α core - α atm - h δ v i mid - c s h δ v z i mid h δ v i mid h δρ i mid - h ρ i ˙ m w - h ρ i mid c s Ideal 2 . . . . .
22 6 . . . . . . .
30 1 . .
13 0 .
34 0 .
12 2 . . . . . .
42 1 . .
24 0 .
32 0 .
16 2 . . . . . .
68 1 . .
45 0 .
35 0 .
27 2 . . . . . . . . .
32 0 .
66 2 . . . . . . . . .
26 1 . . . . . . . . . .
31 0 .
84 3 . . . . . . . . .
23 1 . . . . . . . . . .
26 2 . . . . . . .
37 1 . .
24 0 .
36 0 .
16 3 . . . . . .
59 1 . .
38 0 .
37 0 .
23 2 . . . . . . . . .
39 0 .
51 2 . . . . .
061 0 .
014 0 .
047 0 . .
17 0 . . . . . .
19 0 .
056 0 .
13 0 .
031 0 .
25 0 .
023 0 . . . . . . . .
93 0 .
40 0 .
50 9 . . . . . . .
47 1 . . . . . . . . .
35 0 .
70 3 . Figure 7.
Gas velocity dispersion at the midplane, h δ v i mid , for all runspresented in this study. Panels (a) and (b) plot the data versus α and α core ,respectively. The symbols correspond to models Ideal (circle), X0–X3 (opensquares), Y1–Y4 (triangles), W1–W3 (crosses), X1a–X1d (filled squares),and FS03L (plus sign). The lines show the best linear fits (Equation (21) forpanel (b)). stress α for all our runs. The value of h δ v i mid for each run islisted in Table 2. One can see a rough linear correlation be-tween the velocity dispersion and the accretion stress (for ref-erence, a linear fit h δ v i mid = 0 . α c s is shown by the dashedline). However, detailed inspection shows that h δ v i mid de-creases more rapidly than α as the dead zone increases in size.As found from Table 2, this is because the contribution fromthe atmosphere, α atm , is insensitive to the size of the deadzone in the disk core. In Figure 7(b), we replot the data byreplacing α with the accretion stress in the disk core, α core .Comparison between Figures 7(a) and (b) shows that h δ v i mid more tightly correlates with α core rather than with α . We findthat the data can be well fit by a simple linear relation h δ v i mid = 0 . α core c s , (21)which is shown by the solid line in Figure 7(b). This resultindicates that the accretion stress in the atmosphere does notcontribute to the velocity fluctuation near the midplane.Once h δ v i mid is known, it is also possible to reproduce thevertical profile of the velocity dispersion. For the disk core( | z | < h ideal ), we already know that h δ v i is inversely propor-tional to the mean gas density h ρ i and that h ρ i hardly deviatesfrom the initial Gaussian profile. From these facts, we canpredict the vertical distribution of h δ v i as h δ v i ≈ h δ v i mid h ρ i mid h ρ i ≈ h δ v i mid exp (cid:18) z h (cid:19) ≈ . α core c s exp (cid:18) z h (cid:19) , (22)where Equation (21) has been used in the final equality. InFigure 8, we compare the vertical profiles of the random ve-locity h δ v i / directly obtained from runs Ideal, X1, and X1awith the predictions from Equation (22), where the valuesof α core are taken from Table 2. We see that Equation (22)successfully reproduces the vertical profiles of h δ v i / in thedisk core. We remark that Equation (22) greatly overestimatesthe velocity dispersion at | z | ≫ h ideal , where the gas densitycan no longer be approximated by the initial Gaussian profile(see Figure 4(a)). 5.1.2. Density Fluctuation
Density fluctuations generated by MRI-driven turbulencegravitationally interact with planetesimals and larger solidbodies, affecting their collisional and orbital evolution in pro-toplanetary disks (Laughlin et al. 2004; Nelson & Papaloizou2004; Nelson & Gressel 2010; Gressel et al. 2011). Here, weexamine how the amplitude of the density fluctuations is de-termined the vertically integrated accretion stress.As in Section 5.1.1, we begin with the analysis of the den-sity fluctuations at the midplane, h δρ i / . We find from Ta-ble 2 that h δρ i mid more tightly correlates with α core than withODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES 9 Figure 8.
Vertical profiles of the random velocity h δ v i / directly obtainedfrom MHD simulations (black curves), compared with the predictions fromEquation (22) with α core taken from Table 2 (gray curves). The solid, dashed,and dotted curves correspond to runs Ideal, X1, and X1a, respectively. Thepredicted profiles are plotted only at | z | < h ideal , where Equation (22) is valid. Figure 9.
Gas density fluctuation at the midplane, h δρ i mid / h ρ i , versus α core for all runs presented in this study. The symbols correspond to modelsIdeal (circle), X0–X3 (open squares), Y1–Y4 (triangles), W1–Y3 (crosses),X1a–X1d (filled squares), and FS03L (plus sign). The line shows the bestlinear fit (Equation (23)). α . Figure 9 shows h δρ i mid / h ρ i mid versus α core for all runs.The best linear fit is given by h δρ i mid = 0 . α core h ρ i , (23)which is shown by the solid line in Figure 9. If we use thisequation with Equation (21), we can also obtain the rela-tion between the velocity dispersion and density fluctuation, h δ v i mid / c s = 1 . h δρ i mid / h ρ i . This is consistent with theidea that the fluctuations near the midplane are created bysound waves, for which δ v / c s ∼ δρ/ρ (see also Section 4.1).As shown in Section 4.1, h δρ i is roughly proportional to h ρ i along the vertical direction in the disk core. Hence, if h δρ i / is given, one can reconstruct the vertical profile ofthe density fluctuation in the disk core according to h δρ i ≈ h δρ i mid h ρ ih ρ i mid ≈ h δρ i mid exp (cid:18) - z h (cid:19) ≈ . α core h ρ i exp (cid:18) - z h (cid:19) , (24)where we have used h ρ i ≈ h ρ i mid exp( - z / h ) and Equa- Figure 10.
Outflow mass flux ˙ m w normalized by h ρ i mid c s for all runs pre-sented in this study. Panels (a) and (b) plot the data versus α and α atm , re-spectively. The symbols correspond to models Ideal (circle), X0–X3 (opensquares), Y1–Y4 (triangles), W1–W3 (crosses), X1a–X1d (filled squares),and FS03L (plus sign). The lines show the best linear fits (Equation (26) forpanel (b)). tion (23) in the second and third equalities, respectively.5.1.3. Outflow Flux
We have seen in Section 4.1 and Figure 4(f) that MRIdrives outgoing gas flow at the outer boundaries of the ac-tive layers. The MRI-driven outflow has been first observedby Suzuki & Inutsuka (2009) in shearing-box simulations andbeen recently demonstrated by Flock et al. (2011) in globalsimulations. Suzuki & Inutsuka (2009) and Suzuki et al.(2010) point out that this outflow might contribute to thedispersal of protoplanetary disks, although it is still unclearwhether the outflow can really escape from the disks (see be-low). Meanwhile, MRI also contributes to the accretion of thegas in the radial direction. For consistent modeling of thesetwo effects, we seek how the accretion stress and outflow fluxare correlated with each other.We evaluate the outflow mass flux in the following way. Asseen in Section 4.1, the temporally and horizontally averagedvertical mass flux h ρ v z i is nearly constant at heights | z | & h .Using this fact, we define the outflow mass flux ˙ m w as the sumof |h ρ v z i| averaged near upper and lower boundaries, ˙ m w = R h bnd h bnd - h |h ρ v z i| dzh + R - h bnd + h - h bnd |h ρ v z i| dzh , (25)where h bnd = 5 √ h ≈ h is the height of the upper and lowerboundaries of the simulation box.Figure 10(a) shows ˙ m w normalized by h ρ i mid c s versus α forall simulations. The dimensionless quantity ˙ m w / ( h ρ i mid c s ) isequivalent to C w used in Suzuki et al. (2010). For model Ideal,the value ˙ m w = 5 × - h ρ i mid c s is consistent with the β z = 10 ideal run of Suzuki & Inutsuka (2009). The dashed line showsthe best linear fit ˙ m w / ( h ρ i mid c s ) = 0 . α . It can be seen that0 OKUZUMI & HIROSEthe linear fit captures a rough trend but still considerably over-estimates the outflow flux for models Ideal and Y4. As seenin Table 2, these are the models in which α core dominates over α atm . This implies that the turbulence in the disk core (whichis the source of α core ) does not contribute to the outflow. InFigure 10(b), we replot the data by replacing α with α atm . Wefind that ˙ m w more tightly correlates with α atm than with α .The best linear fit is found to be ˙ m w = 0 . α atm h ρ i mid c s . (26)This result is consistent with the idea that the outflowis driven at the outer boundaries of the active layers(Suzuki & Inutsuka 2009), because the dominant contributionto α atm comes from heights very close to h ideal .Although outflow from the simulation box is a general phe-nomenon in our simulations, it is unclear whether the outflowleaves or returns to the disk. In fact, the outflow velocity ob-served in our simulations does not exceed the sound speedeven at the vertical boundaries. Since the escape velocityis higher than the sound speed, this means that the outflowdoes not have an outward velocity enough to escape out ofthe disk. Acceleration of the outflow beyond the escape ve-locity has not been directly demonstrated by previous simu-lations as well (Suzuki et al. 2010; Flock et al. 2011). How-ever, Suzuki et al. (2010) point out a possibility that magne-tocentrifugal forces and/or stellar winds could accelerate theoutflow to the escape velocity. If the escape of the outflowwill be confirmed in the future, our scaling formula for ˙ m w will certainly become a useful tool to discuss the dispersal ofprotoplanetary disks.5.2. Saturation Predictors for the Accretion Stresses
In the previous subsection, we have shown that the ampli-tudes of various turbulent quantities scale with the verticallyintegrated stresses α core and α atm . The next step is to find outhow to predict α core and α atm in the saturated state from thevertical magnetic flux B z (or equivalently β z ) and the resis-tivity profile η . As shown in Section 4.2, the turbulent stateof a disk depends on the resistivity only through the criticalheights of the dead zone, h Λ and h res . Furthermore, the val-ues of h Λ and h res are only weakly affected by the nonlinearevolution of MRI since the fluctuations in B z and ρ are smallinside the dead zone (see Section 4.1). Therefore, we expectthat the effect of the resistivity can be well predicted by thevalues of h Λ and h res in the initial state, i.e., h Λ , and h res , .With this expectation, we try to derive saturation predictorsfor α core and α atm as a function of β z , h Λ , , and h res , .First, we focus on α core . Figure 11(a) plots α core versus β - z for all our simulations. We see that α core scales roughlylinearly with β - z . The deviation from the linear scaling isexpected to come from the difference in the dead zone size,i.e., h Λ , and h res , . In Figure 12, we plot the product α core β z as a function of h res , . For models except Ideal and Y4, wefind that α core β z is well predicted by a simple formula α core β z = 510 exp( - . h res , / h ) . (27)Figure 11(b) replot the data in Figure 11(a) by replacing β - z with exp( - . h res , / h ) β - z . For models Ideal and Y4, Equa-tion (27) underestimates α core . As explained in Section 4.2,these models exhibit higher magnetic activity near the mid-plane than the other models because of no or a thin dead zone(2 h Λ < h ). We expect that the higher magnetic activity gives Figure 11.
Disk core accretion stress α core for all our simulations. (a) Vs. theinverse initial plasma beta β - z . (b) Vs. β - z multiplied by exp( - . h res , / h )(see also Figure 12). (c) Vs. the final predictor function, Equation (28)(solid line). The symbols correspond to models Ideal (circle), X0–X3 (opensquares), Y1–Y4 (triangles), W1–Y3 (crosses), X1a–X1d (filled squares),and FS03L (plus sign). The dashed lines in panels (a) and (b) are linearfits, shown only for reference. Figure 12. α core β z versus h res , for all runs. The symbols correspond tomodels Ideal (circle), X0–X3 (open squares), Y1–Y4 (triangles), W1–W3(crosses), X1a–X1d (filled squares), and FS03L (plus sign). The solid lineshows an exponential fit for runs except Ideal and Y4 (Equation (27)). ODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES 11
Figure 13.
Atmosphere accretion stress α atm for all our simulations. (a) Vs.the inverse initial plasma beta β - z . (b) Vs. the final predictor function, Equa-tion (29) (solid line). The symbols correspond to models Ideal (circle), X0–X3 (open squares), Y1–Y4 (triangles), W1–W3 (crosses), X1a–X1d (filledsquares), and FS03L (plus sign). The dashed line in panel (a) is a linear fit,only shown for reference. additional contribution to α core . Taking into account this ef-fect, we arrive at the final predictor function, α core = 510 exp( - . h res , / h ) β - z + .
011 exp( - . h Λ , / h ) . (28)Here, the numerical factors 0 .
011 and 3 . α core for all our models. Note thatthe second term of the predictor function is assumed to haveno explicit linear dependence on β - z unlike the first term. Infact, it is possible to reproduce our data by multiplying thesecond term by a prefactor 3 × /β z . However, as we willsee below, the absence of the prefactor makes the predictorfunction consistent with the results of ideal MHD simulationsin the literature.The predictor function for α atm can be obtained in a similarway. Figure 13(a) shows α atm versus β - z for all our runs. Wefind that a simple linear relation α atm = 530 β - z well fits to thedata except for models Ideal and Y4. This means that α atm is characterized only by β - z as long as the dead zone is thick(2 h Λ > h ). To take into account the cases of thin dead zones,we add a term proportional to exp( - . h Λ , / h ) as has beendone for α core , and obtain α atm = 530 β - z + . - . h Λ , / h ) , (29)where the prefactor 0 . α atm for all our runs.Our predictor functions indicate that the vertically inte-grated accretion stress is inversely proportional to β z when Figure 14.
Accretion stress h w xy i versus gas density h ρ i at the upper bound-ary of the active layer ( z = h ideal ) for all runs. The symbols correspond tomodels Ideal (circle), X0–X3 (open squares), Y1–Y4 (triangles), W1–Y3(crosses), X1a–X1d (filled squares), and FS03L (plus sign). The solid lineshows a linear fit (Equation (30)). Figure 15.
Vertically integrated accretion stress α versus the column density Σ active of the active layer for β z = 3 × models. The symbols correspondto models Ideal (circle), X0–X3 (open squares), Y1–Y4 (triangles), W1–Y3(crosses), and FS03L (plus sign). The dashed line is a linear function, onlyshown for reference. a large dead zone is present ( h Λ & h ). As we show below, thisdependence originates from the magnitude of the accretionstress at the outer boundaries of the active layers, | z | ≈ h ideal .When a dead zone exists, the dominant contribution to α comes from the accretion stress at that location (see Figures 5and 6). As shown Figure 14, our simulations suggest that theaccretion stress at | z | = h ideal obeys a simple relation h w xy i ( h ideal ) ≈ . h ρ i ( h ideal ) c s . (30)This means that the averaged accretion stress at | z | = h ideal is18% of the averaged gas pressure h ρ c s i at the same height. Bythe definition of h ideal , the gas density at | z | = h ideal is related to h B z i at the same height as h ρ i ( h ideal ) = (2 π ) h B z i ( h ideal ) / π c s .Since our simulations suggest h B z i ( h ideal ) ∼ B z , the rela-tion means h ρ i ( h ideal ) ∼ π ) B z / π c s ∼ β - z ρ . Usingthis fact, Equation (30) can be rewritten into the linear relationbetween h w xy i ( h ideal ) and β - z : h w xy i ( h ideal ) ∼ β - z ρ c s . (31)When a dead zone is present, the level of α is determined by h w xy i ( h ideal ) (see above), so we have α ∝ β - z .We remark that the vertically integrated stress does notscale linearly with the column density of active layers. Thisis shown in Figure 15, where we compare α with the columndensity Σ active of the active region h Λ < | z | < h ideal . We seethat α decreases much more slowly than Σ active when Σ active isless than 10% of the total gas surface density. This reflects the2 OKUZUMI & HIROSEfact that the dominant contribution to α comes from the outerboundaries of the active zones, | z | ≈ h ideal .It is useful to see how the predictor functions work whena dead zone is absent. If h Λ , = h res , = 0, Equations (28) and(29) predict the total accretion stress α = 1 . × β - z + . α is constant ( α ≈ - ) for β z & andincreases linearly with β - z ( α ≈ - (10 /β z )) for β z . .Strikingly, this prediction is consistent with the finding bySuzuki et al. (2010, see their Figure 2). The existence of thefloor value α ≈ - at low net vertical magnetic fluxes (i.e.,at high β z ) is also supported by recent stratified MHD sim-ulations with zero net flux (Davis et al. 2010). These factssuggest that our predictor functions are applicable even whena dead zone is absent. VERTICAL DIFFUSION COEFFICIENTAs seen in Section 4, sound waves excited in the upperlayers create fluctuations in the gas velocity near the mid-plane. It has been well known that fully developed MRI-driven turbulence causes the diffusion of small dust particles(Johansen & Klahr 2005; Turner et al. 2010). However, it hasnot been fully understood how the sound waves propagatinginside a dead zone affect the dynamics of dust particles there.For example, Suzuki & Inutsuka (2009) speculated that thesound waves might promote dust sedimentation by transfer-ring the downward momentum to dust particles. On the otherhand, Turner et al. (2010) reported that the waves excite ver-tical oscillation of dust particles deep inside the dead zoneand thus prevent the formation of a thin dust layer. Since dustsedimentation is crucial to planetesimal formation via gravita-tional instability, it is worth addressing here how it is affectedby the velocity dispersion created by sound waves.Here, we focus on the dynamics of small dust particles, andmodel the swarm of the particles as a passive scalar as waspreviously done by Johansen & Klahr (2005) and Turner et al.(2010). We assume that dust particles are so small and theirstopping time τ s is much shorter than the turnover time ofturbulence ( ∼ Ω - ). We also assume that the dust density islower than the gas density and hence the dust has no effecton the gas motion. Under these assumptions, the velocity ofdust particles relative to the gas can be approximated by theterminal velocity V T = - Ω τ s z ˆ z , where ˆ z is the unit vector forthe z –direction. Then, the equation of continuity for dust isgiven by ∂ρ d ∂ t + ∇ · [ ρ d ( v + V T )] = 0 , (32)where ρ d is the dust density. Equation (32) has an advantagethat the time step can be taken longer than τ s in numericalcalculation.We have solved Equation (32) for four models (Ideal, X1,X3, and Y1) with the initial condition that the dust-to-gasmass ratio f ≡ ρ d /ρ is constant throughout the simulationbox. To extract the effect of the quasi-stationary turbulence,we insert the dust 100 (for models X1, X3, and Y1) or 250 (formodel Ideal) orbits after the beginning of the MHD calcula-tions. The stopping time τ s is set to τ s = 0 . Ω - for modelIdeal and τ s = 0 . Ω - for the other models. The longer τ s has been adopted for model Ideal to allow the dust to settleappreciably in the stronger turbulence. In reality, the stoppingtime of a dust particle depends on the gas density and henceon z , but we ignore this dependency for simplicity.Figure 16 shows the temporal evolution of the dust density at the midplane, ρ d , mid , for the four MHD runs. The solidcurves show the horizontally averaged ρ d , mid observed in theMHD runs. For comparison, the evolution of ρ d , mid in a hy-drostatic, laminar disk is also shown by the dotted curves. Ir-respectively of the presence or absence of a dead zone, thedust density observed in the MHD runs is higher than that ina laminar disk at all moments. This means that sound wavespropagating in a dead zone do not promote but prevent dustsettling as turbulence does in an active zone.To illustrate more clearly the diffusive nature of the ve-locity dispersion in a dead zone, we try to compare theabove results with a simple advection-diffusion theory. Here,we consider a one-dimensional advection-diffusion equation(Dubrulle et al. 1995) ∂ρ d ∂ t = - ∂ ( ρ d V T ) ∂ z + ∂∂ z (cid:18) D z ρ ∂∂ z ρ d ρ (cid:19) , (33)where V T is the terminal velocity given above and D z is thevertical diffusion coefficient for dust. If the dust particles aresufficiently small ( τ s Ω ≪ D z is equal to the diffusion co-efficient for gaseous contaminants (e.g., Youdin & Lithwick2007). The first term in the right-hand side of Equation (33)represents the downward advection of dust due to settling,while the second term represents the diffusion of dust inthe stratified gas. For disks with no or a small ( | z | . h )dead zone, it is known that Equation (33) well describes theevolution of ρ d if the diffusion coefficient is assumed to be(Fromang & Papaloizou 2006) D z ≈ h δ v z i / Ω . (34)We here examine whether Equations (33) and (34) workwell even when a dead zone is present. We solve Equa-tion (33) with D z = b h δ v z i / Ω , where the vertical distribu-tion of h δ v z i is taken from temporally and horizontally av-eraged MHD data and b is a dimensionless fitting parameter.The dashed curves in Figure 16 show the predictions by theadvection-diffusion model, where b is set to be 1 .
0, 0 .
5, 0 . . b ∼ ρ d , mid for all the models. It is striking that a constant b well re-produces the evolution of the dust density at all heights, as isshown in Figure 17. In this figure, the solid and dashed curvesshow the vertical distribution of the dust-to-gas mass ratio f = ρ d /ρ g observed in run X1 and predicted by the advection-diffusion model with b = 0 .
5, respectively. In this run, theboundaries between the active and dead zones are located at | z | = h Λ ≈ . h (see Table 2). However, Equation (33) suc-cessfully predicts the evolution of f even if we do not changethe value of b across the boundaries. This fact supports theidea that sound waves propagating across a dead zone con-tribute to the diffusion of dust particles just as turbulence doesin active zones.It is worth mentioning here that the diffusion coefficient D z increases with | z | as has been pointed out by Turner et al.(2006) and Fromang & Nelson (2009). This effect is partic-ularly significant at high altitudes where the gas density ismuch lower than that at the midplane, because D z ∝ h δ v z i isroughly proportional to the inverse of the gas density. Thedotted curves in Figure (17) show how Equation (33) would Fromang & Nelson (2009) used a diffusion coefficient quadratic in z toexplain the dust distribution in their MHD simulation. Our finding D z ∝ ρ - ODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES 13
Figure 16.
Temporal evolution of the dust density at the midplane in various MHD runs. The four panels show the results for models Ideal (upper left), X1 (upperright), X3 (lower left), and Y1 (lower right). The solid curves show the horizontally averaged dust density ρ d , mid observed in the MHD runs, while the dotted curveshow the evolution of ρ d , mid in a laminar disk. The dashed curves are the predictions from the one-dimensional advection-diffusion equation (Equation (33)). Thediffusion coefficient D z in Equation (33) is assumed to be proportional to h δ v z i / Ω , and the proportionality factor has been chosen to best reproduce the evolutionof ρ d , mid in the MHD runs. Figure 17.
Snapshots of the vertical distribution of the dust-to-gas mass ra-tio f = ρ d / ρ at t = 100, 250, 400, and 550 orbits for model X1. The solidcurves show the horizontally averaged f observed in the MHD simulation.The dashed curves show the solutions to the one-dimensional advection-diffusion equation (Equation (33)) with D z ( z ) = 0 . h δ v z i ( z ) / Ω . The dottedcurves show the solution to Equation (33) with a constant diffusion coeffi-cient D z = 0 . h δ v z i mid / Ω . fail to predict dust evolution if one assumed a constant dif-fusion coefficient D z = 0 . h δ v z i mid / Ω . We see that the con-stant diffusion coefficient model significantly underestimatesthe dust density at | z | ≫ h . This fact will merit considera- does not contradict their assumption because ρ - ∝ exp( z / h ) ≈ + z / h near the midplane. tion when modeling the chemical evolution of protoplanetarydisks, in which the vertical mixing of molecules is of impor-tance (Heinzeller et al. 2011).Finally, we give a simple analytic recipe for the vertical dis-tribution of D z . It is useful to rewrite Equation (34) in termsof h δ v i , for which the scaling relation (Equation (22)) andpredictor function (Equation (28)) are available. Table 2 liststhe ratio of h δ v z i mid to h δ v i mid for all our simulations. It canbe seen that h δ v z i mid ≈ (32% ± × h δ v i mid , indicatingthat h δ v z i mid is roughly equal to a third of h δ v i mid . Fur-thermore, the ratio h δ v z i / h δ v i is approximately constant inthe disk core, as is illustrated in Figure 4(c). Based on thesefacts, we approximate h δ v z i as h δ v i / h δ v i (Equation (22)), we rewrite Equation (34) as D z ≈ h δ v i / Ω ≈ . α core c s / Ω ) exp (cid:18) z h (cid:19) . (35)If one uses this equation together with the predictor functionfor α core (Equation (28)), one can calculate the vertical distri-bution of D z in the disk core for given β z and η . DISCUSSION: EFFECTS OF NUMERICALRESOLUTIONAll MHD simulations presented in the previous sectionswere performed with the numerical resolution of 40 × ×
200 grid cells for the simulation box of size 2 √ h × √ h × √ h . Here, we examine how the numerical resolution af-fects the saturated state of turbulence.4 OKUZUMI & HIROSE Figure 18.
Saturated values of various quantities versus numerical res-olution for model X1. From top to bottom: α , α core , h δ v i mid / c s , h δρ i mid / h ρ i , and ˙ m w / h ρ i mid c s . The horizontal axis shows the numberof grid cells per length √ h in the vertical direction, n z . The value n z = 20corresponds to the resolution adopted in this study. We carry out X1 simulations with changing the numericalresolution to 20 × ×
100 cells and 80 × ×
400 cells.Figure 18 compares the saturated values of various quantitiesobtained from the two runs with the values from the originalX1 run (Table 2). Here, the horizontal axis shows the numberof grid cells per length √ h in the vertical direction, n z (seefootnote 1). The value n z = 20 corresponds to our originalresolution. We see that the change of the resolution hardly af-fects the integrated accretion stresses α and α core and outflowflux ˙ m w , suggesting that the resolution of n z = 20 is sufficientfor these quantities to converge well. By contrast, the veloc-ity and density dispersions h δ v i mid and h δρ i mid increase withimproving the resolution. Since the energy input rate to tur-bulence should be the same if the integrated accretion stressis unchanged, the resolution dependence of the velocity anddensity dispersions is expected to mainly come from artificialdissipation of sound waves in the simulation box. However,we also see that this effect becomes less significant as the res-olution is improved. Detailed inspection shows that the frac-tional increase in h δ v i mid is 81% when going from n z = 10 to n z = 20 but is 40% when going from n z = 20 to n z = 40. Thissuggests that the amplitudes of the velocity and density fluc-tuations should converge to finite values in the limit of highresolutions ( n z → ∞ ). This is to be expected, since soundwaves in a stratified disk physically dissipate through, e.g.,shock formation, particularly at high altitudes where the gasdensity is low and therefore the amplitudes of the waves be-come large. We find that the data for h δ v i mid / c s shown inFigure 18 lie on a curve h δ v i mid / c s = 0 . - . n - . z .This implies a converged value of h δ v i mid / c s ≈ . n z = 20 simula-tion ( h δ v i mid / c s ≈ . h δ v i and h δρ i could be underestimated by a factor of severalin the simulations presented in this study.In summary, we find that the outflow mass flux and verti-cally integrated accretion stress converge well within our nu-merical resolution. This suggests that the predictor functionsfor α core and α atm (Equations (28) and (29)) and the scalingrelation between ˙ m w and α atm are hardly affected by the res-olution. On the other hand, the amplitudes of sound wavescould be underestimated by a factor of several because of thefinite grid size. Future high-resolution simulations will enableto better quantify the scaling relations between h δ v i and α core and between h δρ i and α core (Equations (22) and (24)). SUMMARYGood knowledge about the turbulent structure of protoplan-etary disks is essential for understanding planet formation.To provide an empirical basis for modeling the coevolutionof dust and MRI, we have performed MHD simulations of avertically stratified shearing box with an MRI-inactive “deadzone” of various sizes and with a vertical magnetic flux ofvarious strengths. Our findings are summarized as follows.1. We have introduced the critical heights ( h ideal , h Λ , and h res ) that characterize the MRI in a stratified disk (Sec-tion 3). We have found that the vertical structure ofMRI-driven turbulence depends on the resistivity pro-file only through the critical heights for the dead zone( h Λ and h res ) and is insensitive to the detail of the resis-tivity profile (Section 4.2).2. In the “disk core” ( | z | < h ideal ), the density-weighted ve-locity dispersion h ρ ih δ v i is nearly constant along thevertical direction (Section 4.1). This means that thevelocity dispersion is approximately inversely propor-tional to the gas density. Weak dependence on z is alsofound for h δρ i / h ρ i , meaning that the density fluctua-tion h δρ i / is proportional to the square root of theaveraged density.3. The accretion stresses in the disk core and “atmo-sphere” ( | z | > h ideal ) differently contribute to the tur-bulent structure of a disk (Section 5.1). The velocitydispersion h δ v i and density fluctuation h δρ i in thedisk core depend linearly on the accretion stress inte-grated over the core, α core (Equations (21) and (23)). Bycontrast, the outflow mass flux ˙ m w depends linearly onthe stress integrated over the atmosphere, α atm (Equa-tion (26)).4. We have obtained simple empirical formulae that pre-dict the vertically integrated stresses α core and α atm inthe saturated state (Section 5.2; Equations (28) and(29)). These are written as a function of the strengthof the vertical magnetic flux (or β z ) and the criticalheights of the dead zone measured in the nonturbulentstate ( h res , and h Λ , ). These predictor functions to-gether with the saturation relations described above al-low to calculate various turbulent quantities for a givenresistivity profile and a net vertical flux.5. We have confirmed that the vertical diffusion coeffi-cient D z of contaminants is given by D z ≈ h δ v z i / Ω bothinside and outside a dead zone (Section 6). This impliesthat sound waves propagating across a dead zone con-tribute to the diffusion of dust particles just as turbu-lence does in active zones. We have obtained a simpleanalytic recipe for the vertical distribution of D z as afunction of α core on the basis of our MHD simulationdata (Equation (35)).The empirical formulae obtained in this study enable us topredict the amplitudes of various turbulent quantities in a pro-toplanetary disk with a dead zone. The steps to be performedare as follows.1. Prepare the vertical profile of the ohmic resistivity η ,and find h Λ and h res from Equations (12) and (13).ODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES 15A realistic profile of η in the presence of dust parti-cles can be obtained by solving the ionization state ofthe gas and the charge state of the dust simultaneously(e.g., Sano et al. 2000; Ilgner & Nelson 2006; Okuzumi2009).2. Calculate α core and α atm using the predictor functions,Equations (28) and (29).3. One can now calculate the turbulent viscosity ofthe disk as ν turb = (3 / α core + α atm ) c s / Ω (see Equa-tions (18), (19), and (20)). The vertical distribution ofthe gas velocity dispersion and density fluctuation inthe disk core ( | z | < h ideal ) can be calculated from Equa-tions (22) and (24), respectively. The outflow mass fluxcan be evaluated from Equation (26). For the diffusioncoefficient in the disk core, one can use Equation (35).When using our empirical formulae, it should be kept inmind that our scaling relations for velocity and density fluc-tuations (Equation (22) and (24)) could underestimate theirmean-squared amplitudes relative to the integrated accretionstress by a factor of several because of the numerical dissipa-tion of sound waves (Section 7). Future high-resolution sim-ulations will allow to better quantify the saturation level ofsound wave amplitudes.We are grateful to Neal Turner and Takayoshi Sano for pro-viding us with their MHD simulation data that have motivatedus to start this study. We also thank Shu-ichiro Inutsuka,Takeru Suzuki, Taku Takeuchi, Hidekazu Tanaka, TakayukiTanigawa, Mordecai-Mark Mac Low, and the anonymous ref-eree for useful discussion and fruitful comments. Calculationswere made on the Cray XT4 at the CfCA, National Astronom-ical Observatory of Japan, and SR16000 at the Yukawa Insti-tute for Theoretical Physics, Kyoto University. S.O. is sup-ported by a Grant-in-Aid for JSPS Fellows (22 · z above the midplane.Recombination occurs in the gas phase and on dust sur-faces. The gas-phase recombination dominates if the totalsurface area of dust particles is negligibly small. In this case,the equation for the ionization-recombination equilibrium isgiven by ζ n n = γ ie n i n e = γ ie n e , (A1)where ζ is the ionization rate (the probability per unit time atwhich a molecule is ionized), n n , n i , and n e are the numberdensities of neutrals, ions, and electrons, respectively, and γ ie is the gas-phase recombination rate coefficient. The secondequality in the above equation assumes the charge neutrality in the gas phase, n i = n e . Equation (A1) leads to the electronabundance x e = n e n n = s ζγ ie n n ∝ p ζ exp (cid:18) z h (cid:19) , (A2)which means that the resistivity is proportional to ζ - / exp( - z / h ). Thus, if the vertical dependence of ζ can be neglected, the resistivity profile is given by Equa-tion (4) with h η = √ h . Note that Equation (5) is derivedinstead of Equation (4) if cosmic-ray ionization is assumedand the attenuation of cosmic rays toward the midplane istaken into account (Fleming & Stone 2003).If the total surface area of dust particles is large, recombina-tion occurs mainly on dust surfaces. In this case, the equationfor the ionization-recombination equilibrium is given by ζ n n = γ de n d n e , (A3)where γ de is the sticking rate coefficient for dust–electron col-lision and n d is the number density of dust particles. Thesticking rate coefficient depends on the charge of the dustparticles, and ultimately on n e via the charge neutrality (seeOkuzumi 2009), but we will ignore this dependence in thefollowing. From Equation (A3) , we have x e = ζγ de n d ∝ ζ exp (cid:18) z h d (cid:19) , (A4)where we have assumed that n d ∝ exp( - z / h d ) with h d be-ing the scale height of the dust particles (in fact, one canshow using Equation (33) that n d obeys a Gaussian distribu-tion in sedimentation-diffusion equilibrium if D z ∝ τ s ∝ n - g ∝ exp( z / h ); the condition τ s ∝ n - g is satisfied if the size ofthe dust particles is smaller than the mean free path of thegas). Thus, ignoring the dependence of ζ on z , the resistivity η ∝ x - e is given by Equation (4) with h η = h d . REFERENCESBai, X.-N., & Stone, J. 2010, ApJ, 722, 1437Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163Brauer, F., Dullemond, C. P., & Henning, Th. 2008, A&A, 480, 859Chambers, J. E., & Wetherill, G. W. 1998, Icarus, 136, 304Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ,546, 496Cuzzi, J. N., Hogan, R. C., & Shariff, K. 2008, ApJ, 687, 1432Cuzzi, J. N., Hogan, R. C., Bottke, W. F. & 2010, Icarus, 208, 518Davis, S. W., Stone, J. M., Pessah, M. E. 2010, ApJ, 713, 52Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971Flaig, M., Kley, W., & Kissmann, R., MNRAS, 409, 1297Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, Th. 2011,ApJ, 735, 122Fromang, S., & Nelson, R. P. 2009, A&A, 496, 597Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751Gammie, C. F. 1996, ApJ, 457, 355Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051Gressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742Heinzeller, D., Nomura, H., Walsh, C., & Millar, T. J. 2011, ApJ, 731, 115Hirose, S., Krolik, J. H., & Blaes, O. 2009, ApJ, 691, 16Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292Ilgner, M., & Nelson, R. P. 2006, A&A, 445, 205Johansen, A., Brauer, F., Dullemond, C., Klahr, H., & Henning, T. 2008,A&A, 486, 597Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353Johansen, A., et al. 2007, Nature, 448, 1022Johansen, A., & Youdin, A. 2007, ApJ, 662, 6276 OKUZUMI & HIROSE