MODELING MRI TURBULENCE IN PROTOPLANETARY DISKS WITH DEAD ZONES

# Modeling Magnetorotational Turbulence in Protoplanetary Disks with Dead Zones

Satoshi Okuzumi Department of Physics, Nagoya University, Nagoya, Aichi 464-8602, Japan; okuzumi@nagoya-u.jp Shigenobu Hirose Institute for Research on Earth Evolution, JAMSTEC, Yokohama, Kanagawa 236-0001, Japan
###### Abstract

Turbulence driven by magnetorotational instability (MRI) crucially affects the evolution of solid bodies in protoplanetary disks. On the other hand, small dust particles stabilize MRI by capturing ionized gas particles needed for the coupling of the gas and magnetic fields. To provide an empirical basis for modeling the coevolution of dust and MRI, we perform three-dimensional, ohmic-resistive MHD simulations of a vertically stratified shearing box with an MRI-inactive “dead zone” of various sizes and with a net vertical magnetic flux of various strengths. We find that the vertical structure of turbulence is well characterized by the vertical magnetic 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 to the details of the resistivity profile. We discover scaling relations between the amplitudes of various turbulent quantities (velocity dispersion, density fluctuation, vertical diffusion coefficient, and outflow mass flux) and vertically integrated accretion stresses. We also obtain empirical formulae for the integrated accretion stresses as a function of the vertical magnetic flux and the critical heights. These empirical relations allow to predict the vertical turbulent structure of a protoplanetary disk for a given strength of the magnetic flux and a given resistivity profile.

###### Subject headings:
dust, extinction — planets and satellites: formation — protoplanetary disks
slugcomment: ApJ Accepted

## 1. Introduction

Planets are believed to form in protoplanetary gas disks. The standard scenario for planet formation consists of the following steps. Initially, submicron-sized dust grains grow into kilometer-sized planetesimals by collisional sticking and/or gravitational instability (Safronov, 1969; Goldreich & Ward, 1973; Weidenschilling & Cuzzi, 1993). Planetesimals undergo further growth toward Moon-sized protoplanets through mutual collision assisted by gravitational interaction (Wetherill & Stewart, 1989). Accretion of the disk gas onto protoplanets leads to the formation of gas giants (Mizuno, 1980; Pollack et al., 1996), while terrestrial planets form through the giant impacts of protoplanets after the gas disk disperses by viscous accretion onto the central star and other effects (Chambers & Wetherill, 1998).

Turbulence in protoplanetary disks plays a decisive role on planet formation as well as on disk dispersal. The impact of turbulence is particularly strong on the formation of planetesimals since the frictional coupling of gas and dust particles governs the process. Classically, planetesimal formation has been attributed to the collapse of a dust sedimentary layer by self-gravity (Safronov, 1969; Goldreich & Ward, 1973) and/or the collisional growth of dust grains (Weidenschilling & Cuzzi, 1993). The presence of strong turbulence is preferable for dust growth when the dust particles is so small that Coulomb repulsion is effective (Okuzumi, 2009; Okuzumi et al., 2011). However, strong turbulence acts against the growth of macroscopic dust aggregates since it makes their collision disruptive (Weidenschilling, 1984; Johansen et al., 2008). Furthermore, turbulence causes the diffusion of a dust sedimentary layer, making planetesimal formation via gravitational instability difficult as well (Weidenschilling, 1984). Turbulence is also known to concentrate dust particles of particular sizes, but its relevance to planetesimal formation via gravitational instability is still under debate (Cuzzi et al., 2001, 2008, 2010; Pan et al., 2011). More recently, it has been suggested that two-fluid instability of gas and dust can produce dust clumps with density high enough for gravitational collapse, but successful dust coagulation to macroscopic sizes seems to be still required for this mechanism to become viable (Youdin & Goodman, 2005; Johansen & Youdin, 2007; Johansen et al., 2007; Bai & Stone, 2010). Besides, turbulence also affects planetesimal growth as turbulent density fluctuations gravitationally interact with planetesimals and can raise their random velocities above the escape velocity (Ida et al., 2008; Nelson & Gressel, 2010). The fluctuating gravitational field can even cause random orbital migration of protoplanets (Laughlin et al., 2004; Nelson & Papaloizou, 2004). Thus, to understand the growth of solid bodies in various stages, it is essential to know the strength and spatial distribution of disk turbulence.

Interestingly, the evolution of solid bodies is not only affected but also affects disk turbulence. The most viable mechanism for generating disk turbulence is the magnetorotational instability (MRI; Balbus & Hawley 1991). This instability has its origin in the interaction between the gas disk and magnetic fields, and therefore requires a sufficiently high ionization degree to operate. Importantly, whether the MRI operates or not in each location of the disk is strongly dependent on the amount of small dust grains because they efficiently capture ionized gas particles and thus reduce the ionization degree (Sano et al., 2000; Ilgner & Nelson, 2006; Okuzumi, 2009). This implies that dust and MRI-driven turbulence affect each other and thus evolve simultaneously.

The purpose of this study is to present an empirical basis for studying the coevolution of solid particles and MRI-driven turbulence. It is computationally intensive to simulate the evolution of dust and MRI-driven turbulence simultaneously, since the evolutionary timescale of solid bodies is generally much longer than the dynamical timescale of the turbulence. For example, turbulent eddies grow and decay on a timescale of one orbital period (e.g., Fromang & Papaloizou, 2006), while dust particles grow to macroscopic sizes and settle 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 evolutionary stage of dust evolution. Motivated by this fact, we restrict ourselves to time-independent ohmic resistivity, but instead focus on how the quasi-steady structure of turbulence depends on the vertical profile of the resistivity.

To characterize the vertical structure of MRI-driven turbulence, we perform a number of three-dimensional MHD simulations of local stratified disks including resistivity and nonzero net vertical magnetic flux. Inclusion of a nonzero net vertical flux is important as it determines the saturation level of turbulence (Hawley et al., 1995; Sano et al., 2004; Suzuki & Inutsuka, 2009, see also our Section 4). Similar simulations have been done in a number of previous studies (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 & Turner, 2011). One important difference between our study and previous ones is that we focus on general dependence of the saturated turbulent state on the model parameters such as the resistivity and net magnetic vertical flux.

Our modeling of MRI-driven turbulence follows two steps. In the first step, we seek scaling laws giving the relations among turbulent quantities. We express the relations as a function of the vertically integrated accretion stress, which is the quantity that determines the rate at which turbulent energy is extracted from the differential rotation (Balbus & Papaloizou, 1999). As we will see, excellent scaling relations are obtained if we divide the integrated stress into two components that characterize the contributions from different regions in the stratified disk (which we will call the “disk core” and “atmosphere”) In the second step, we find out empirical formulae that predict the vertically integrated stresses as a function of the resistivity profile and vertical magnetic flux.

The plan of this paper is as follows. In Section 2, we describe the method and setup used in our MHD simulations. In Section 3, we introduce “critical heights” derived from the linear analysis of MRI in stratified disks. As we will see later, these critical heights are useful to characterize the turbulent structure observed in our simulations. We present our simulation results in Section 4, and obtain scaling relations and predictor functions for the quasi-steady state of turbulence in Section 5. In Section 6, we simulate dust settling in a dead zone to model the diffusion coefficient for small particles as a function of height. Effects of numerical resolutions on our simulation results are discussed in Section 7. Our findings are summarized in Section 8.

## 2. Simulation setup and Method

In this section, we describe the setup and method adopted in our stratified resistive MHD simulations.

### 2.1. Setup

Our MHD simulations adopt the shearing box approximation (Hawley et al., 1995). We consider a small patch of disk centered on the midplane of an arbitrary distance from the central star, and model it as a stratified shearing box corotating with the angular speed at the domain center. We use the Cartesian coordinate system , where , , and stand for the radial, azimuthal, and vertical distance from the domain center. In addition, we assume that the gas is isothermal throughout the box; thus, the sound velocity of the gas is constant in both time and space.

#### 2.1.1 Initial Conditions

For the initial condition, we assume that the gas disk is initially in hydrostatic equilibrium and is threaded by uniform vertical magnetic field . The assumption of the hydrostatic equilibrium leas to the initial gas density profile

 ρ=ρ0exp(−z22h2), (1)

where is the initial gas density at the midplane and

 h=csΩ (2)

is the pressure scale height.111Note that the “gas scale height” is often defined as in the literature on stratified MHD simulations. The ratio between the initial midplane gas pressure and the initial magnetic pressure defines the initial plasma beta

 βz0≡8πρ0c2sB2z0. (3)

In this paper, the strength of the initial magnetic flux will be referred by rather than .

#### 2.1.2 Resistivity Profile

The main purpose of this study is to see how the turbulence depends on the vertical profile of the ohmic resistivity. We adopt a simple analytic resistivity profile based on the following consideration. For fixed temperature, the resistivity is inversely proportional to the ionization degree (Blaes & Balbus, 1994). In protoplanetary disks, the ionization degree at each location is determined by the balance of ionization (by, e.g., cosmic rays and X-rays) and recombination (in the gas phase and on grain surfaces). Detailed structure of the resistivity profile depends on what processes dominate the ionization and recombination. However, a general tendency is that the ionization degree decreases toward the midplane of the disk, because the ionization rate is lower as the column depth is greater and because the recombination rate is higher as the gas density is higher (see, e.g., Sano et al., 2000). Based on this fact, we give the resistivity profile such that increases as decreases. To be more specific, we adopt the following resistivity profile

 η=ηmidexp(−z22h2η), (4)

where is the resistivity at the midplane and is the scale height of .

Equation (4) satisfies the important property of realistic resistivity profiles mentioned above. Furthermore, as shown in Appendix, Equation (4) exactly reproduces the vertical resistivity profile of a disk in some limited cases. However, it will be useful to examine possible influences of limiting to Equation (4). In order to do that, we also consider a resistive profile used in Fleming & Stone (2003),

 ηFS03=η0exp(−z24h2)exp(Σ0ΣCR12√2π∫∞|z|/he−z′2/2dz′), (5)

which is characterized by two parameters and (see Equation (10) of Fleming & Stone 2003). Physically, Equation (5) corresponds to the resistivity profile when the ionization degree is determined by the balance between cosmic-ray ionization 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 various sets of the parameters (, , ). Table 1 lists the parameters adopted in each run. The model “Ideal” assumes zero resistivity throughout the simulation box. Models X0–X3 are defined by the same value of but different values of . The difference between X, Y, and W models is in the value of . The initial plasma beta is taken to be in all models except X1a–X1d. In addition, we construct one model from Equation (5) with and . This set of parameters corresponds to the “larger dead zone” model of Fleming & Stone (2003) defined by and , where is the magnetic Reynolds number with the typical length and velocity set to be and , respectively. We refer to the model with these parameters as model FS03L. The initial plasma beta in the FS03L model is taken to be . Note that our FS03L run is not exactly the same as the larger dead zone run of Fleming & Stone (2003) since they assumed zero net vertical magnetic flux.

### 2.2. Method

We solve the equations of the isothermal resistive MHD using the ZEUS code (Stone & Norman, 1992). The domain is a box of size along the radial, azimuthal, and vertical directions, divided into grid cells. For all runs, the radial and azimuthal boundary conditions are taken to be shearing-periodic and periodic, respectively. Therefore, the vertical magnetic flux is a conserved quantity of our simulations. The vertical boundary condition for all runs except X1d is the standard ZEUS outflow condition, where the fields on the boundaries are computed from the extrapolated electromotive forces; only for run X1d, we assume for numerical stability that the magnetic fields are vertical on the top and bottom boundaries as done by Flaig et al. (2010). We have checked that the results hardly depend on the choice of the two types of boundary conditions. Note that the radial and azimuthal components of the mean magnetic fields are not conserved due to the outflow boundary condition. We also added a small artificial resistivity near the boundaries for numerical stability (Hirose et al., 2009). A density floor of is applied to prevent high Alfvén speeds from halting the calculations.

## 3. Characteristic Wavelengths and Critical Heights

As we will see in the following sections, it is useful to analyze the simulation results using the knowledge obtained from the linear analysis of MRI. In this subsection, we introduce several quantities that characterize the linear evolution of MRI.

### 3.1. Characteristic Wavelengths of MRI

According to the local linear analysis of MRI including ohmic resistivity (Sano & Miyama, 1999), the wavelength of the most unstable MRI mode can be approximately expressed as

 λlocal≈max{λideal,λres}, (6)

where

 λideal≡2πvAzΩ (7)

and

 λres≡2πηvAz (8)

are the characteristic wavelengths of MRI modes in the ideal and resistive MHD limits, respectively, and is the vertical component of the Alfvén velocity. Equation (6) can be written as , where is the Elsasser number defined by

 Λ≡v2AzηΩ. (9)

The Elsasser number determines the growth rate of the MRI. If , ohmic diffusion does not affect the most unstable mode lying at , and the local instability occurs rapidly at a wavelength and at a rate . If , ohmic diffusion stabilizes the most unstable mode, and the local instability occurs at a longer wavelength and at a slower rate . We will refer to the former case as the “ideal MRI,” and to the latter case as the “resistive MRI.”

Figure 1 schematically illustrates how and vary with height . In general, grows toward higher because the Alfvén speed increases as the density decreases. By contrast, grows toward lower because is inversely proportional to and because increases with decreasing (see the discussion in Section 2.1.2).

The global instability of a stratified disk can be described in terms of the local analysis. As shown by Sano & Miyama (1999), the gas motion at height is unstable if the local unstable wavelength is shorter than the scale height of the disk, i.e.,

 max{λideal,λres}≲h. (10)

### 3.2. Critical Heights

With the global instability criterion (Equation (10)) together with the vertical dependence of and , we can define three different critical heights for a stratified disk.

1. The first one is defined by

 λideal(z=hideal)=h, (11)

or equivalently, , where . At (), MRI does not operate because the wavelengths of the unstable modes exceed the disk thickness (Sano & Miyama, 1999). We refer to the region as the “atmosphere” and to the region as the “disk core.”

2. The second one is defined by

 Λ(z=hΛ)=1, (12)

or equivalently, . The layer is the so-called “active layer,” where MRI operates without affected by ohmic diffusion nor gas stratification. The region is what we call the “dead zone,” where ohmic diffusion stabilizes the most unstable ideal MRI mode. For convenience, we regard as zero when a dead zone is absent. This is the case for model Ideal.

3. The third one is defined by

 λres(z=hres)=h. (13)

Ohmic diffusion allows the resistive MRI to operate at . At , ohmic diffusion stabilizes all the unstable MRI modes. Note that some previous studies (e.g., Gammie, 1996; Sano et al., 2000) used the terminology “dead zone” for the region rather than . In fact, as we will see in Section 4, the set of and best characterizes our dead zone. We regard as zero when is less than at all heights. This is the case for Y models.

The critical heights in the initial state (, , and ) are shown in Table 1 for all of our 17 simulations. Using Equations (1) and (4), one can analytically calculate the initial critical heights for models except FL03L as

 hideal,0=[2ln(βz08π2)]1/2h, (14)
 hΛ,0=(2lnΛ−101+(h/hη)2)1/2h, (15)
 hres,0=[2ln(8π2β−1z0Λ−20)1+2(h/hη)2]1/2h, (16)

where

 Λ0=2cshηmidβz0 (17)

is the initial Elsasser number at the midplane.

Figure 2 shows the vertical profiles of the resistivity and the initial Elsasser number for some of our models. The initial midplane Elsasser number and initial critical heights () are listed in Table 1 for all models. As one can see from Table 1 and the lower panel of Figure 2, models labeled by the same number are arranged so that they have similar values of .

For turbulent states, we evaluate in the Elsasser number and the characteristic wavelengths as , where the overbars denote the horizontal averages.

## 4. Simulation Results

### 4.1. The Fiducial Model

We select model X1 as the fiducial model to describe in detail. Figure 3 shows how MRI-driven turbulence reaches a quasi-steady state in run X1. The upper and lower panels plot the horizontal averages of the Maxwell stress and the density-weighted velocity dispersion , respectively, as a function of time and height . The solid, dotted, and dashed lines are the loci of the critical heights , , and , respectively. As seen in the figure, a quasi-steady state is reached within the first 40 orbits. The critical height measured in the quasi-steady state is slightly lower than that in the initial nonturbulent state. This is because the ideal MRI wavelength is increased by the fluctuation in the vertical magnetic field, . By contrast, and are almost unchanged, because the fluctuation of the magnetic field is suppressed in the dead zone.

Figure 4 shows the vertical structure of the disk averaged over a time interval . The dark and light gray bars in each panel indicate the heights where ideal and resistive MRIs operate, respectively (see also Figure 1). The brackets denote the averages over time and horizontal directions.

In Figure 4(a), we compare the averaged gas density with the initial density given by Equation (1). We see that the density is almost unchanged in the disk core () but is considerably increased in the atmosphere (). This is because the magnetic pressure is negligibly small in the disk core but dominates over the gas pressure in the atmosphere. Figure 4(a) also shows the amplitude of the density fluctuation, . As one can see, the density fluctuation is small () except at .

The magnetic activity in the disk can be seen in Figure 4(b), where the vertical profiles of the magnetic energies (, , and ) and Maxwell stress are plotted. One can see that these quantities peak near the outer boundaries of the active layers, . This is because the largest channel flows develop at locations where (see, e.g., Suzuki & Inutsuka, 2009). In the dead zone, ohmic dissipation suppresses the fluctuation in the magnetic fields, , leaving the initial vertical field () and coherent toroidal fields () generated by the differential rotation.222 In our simulations, ohmic resistivity is not high enough to remove shear-generated, coherent toroidal fields. In this sense, our dead zone is an “undead zone” in the terminology of Turner & Sano (2008).

Figure 4(c) shows the density-weighted velocity dispersions and and the Reynolds stress . These quantities characterize the kinetic energy in the random motion of the gas.333 In the disk core (), is approximately equal to since the density fluctuation is small (see Figure 4(a)). Comparing Figures 4(b) and (c), we find that the drop in these quantities in the disk core is not as significant as the drop in and . This is an indication that sound waves generated in the active layers penetrate deep inside the dead zone (Fleming & Stone, 2003). Furthermore, we find that is approximately constant, i.e., the velocity dispersion is inversely proportional to the mean density , in the disk core. This means that the kinetic energy density of fluctuation is nearly constant in the disk core. This is another indication of sound waves, because the amplitude of the velocity fluctuation is generally proportional to for freely propagating sound waves. The root-mean-squared random velocity is shown in Figure 4(d). The random velocity is subsonic in the disk core and exceeds the sound speed only in the atmosphere.

An indication of freely sound waves can be also found in the density fluctuation. Shown in Figure 4(e) is the mean squared density fluctuation divided by the mean density . Since , the quantity is approximately proportional to the thermal energy density of fluctuation, . In the disk core, we see that is roughly constant along the vertical direction, meaning that the amplitude of the density fluctuation, , is proportional to the square root of the mean density . The similarity between and is peculiar to sound waves, for which .

Figure 4(f) displays the profile of the vertical mass flux . In the atmosphere, the vertical flux is outward, i.e., at and at . This outflow results from the breakup of large channel flows at the outer boundaries of the active layers, (Suzuki & Inutsuka, 2009; Suzuki et al., 2010). The vertical mass flux reaches a constant value at . This fact allows us to measure a well-defined outflow flux for each simulation (see Section 5.1.3).

### 4.2. Model Comparison

We now investigate how the vertical structure of turbulence depends on the resistivity profile and vertical magnetic flux.

Figure 5 displays the temporal and horizontal averages of the Maxwell stress and the density-weighted velocity dispersion as a function of for various models. As in Figures 1 and 4, the dark and light gray bars in each panel indicate the heights where ideal and resistive MRIs operate, respectively.

The effect of changing the size of the dead zone can be seen in Figure 5(a), where models X1, X2, and X3 are compared. These models are characterized by the same values of and but different values of . For all the models, sharply falls at , meaning that well predicts where the resistivity shuts off the magnetic activity. By contrast, exhibits a flat profile at with no distinct change across nor . The only clear difference is the value of in the disk core, i.e., the value is lower when the dead zone is wider. Note that decreases more slowly than the column density of the active layers . For example, the active column density in model X1 is 20 times smaller than that in model X3. However, the midplane value of in the former is only five times smaller than that in the latter. This suggests that even a very thin active layer can provide a large velocity dispersion near the midplane.

Interestingly, the vertical structure of turbulence depends on the critical heights (, , and ) but are very insensitive to the details of the resistivity profile. This can be seen in Figure 5(b), where we compare runs with similar critical heights (runs X3, W3, and FS03L). We see that these models produce very similar vertical profiles of and even though they assume quite different resistivity profiles (see the upper panel of Figure 2). This suggests that the vertical structure of turbulence is determined by the values of the critical heights.

Importance of distinguishing and is illustrated in Figures 5(c) and (d). These panels compare five models (Y1–4 and Ideal) in which the resistive MRI is active at the midplane, i.e., . Figure 5(c) shows models with . We see that the profiles of and are very similar for the three models. This implies that the vertical structure is determined by the value of when . Figure 5(d) shows what happens when falls below . Model Ideal is clearly different from the other models. Model Y4 () is interesting because it exhibits both features. In the lower half of the disk (), the Maxwell stress behaves as in the other Y models. In the upper half (), however, the profile of is closer to that in model Ideal.

Next, we see how the saturation level of turbulence depends on the vertical magnetic flux. Figure 6 shows the vertical profiles of and for three runs with different values of (X1b, X1, and X1d). We see that these values increase with decreasing . The peak value of is approximately , , and for runs X1b, X1, and X1d, respectively ( is the initial midplane gas pressure). This indicates a linear scaling between the turbulent stress and .

## 5. Scaling Relations and Predictor Functions

Now we seek how the amplitudes of turbulent quantities depend on the vertical magnetic flux and the resistivity profile. We do this in two steps. First, we derive relations between the amplitudes of turbulent quantities and the vertically integrated turbulent stress. We then obtain empirical formulae that predict the integrated stress as a function of the vertical magnetic flux and the resistivity profile.

### 5.1. Scaling Relations between Turbulent Quantities and Vertically Integrated Accretion Stresses

The ultimate source of the energy of turbulence is the shear motion of the background flow. The accretion stress determines the rate at which the free energy is extracted. Therefore, we expect that the accretion stress is related to the amplitudes of turbulent quantities, such as the gas velocity dispersion and outflow mass flux.

To quantify the rate of the energy input in the simulation box, we introduce the effective parameter

 α≡∫⟨wxy⟩dzΣc2s, (18)

where is the gas surface density. In the classical, one-dimensional viscous disk theory (Lynden-Bell & Pringle, 1974), the parameter is related to the turbulent viscosity as , where the prefactor comes from the slope of the Keplerian rotation. Thus, also characterizes the vertically integrated mass accretion rate.

As we will see below, it is useful to decompose as , where

 αcore≡∫|z|

and

 αatm≡∫|z|>hideal⟨wxy⟩dzΣc2s (20)

are the contributions from the disk core () and atmosphere (), respectively. Table 2 shows the values of , , and as well as the time-averaged critical heights (, , ) for all our simulations.

#### 5.1.1 Velocity Dispersion

Random motion of the gas crucially affects the growth of dust particles as it enhances the collision velocity between the particles via friction forces (Weidenschilling, 1984; Johansen et al., 2008). Here, we seek how the velocity dispersion is related to the integrated accretion stress.

First, we focus on the velocity dispersion at the midplane, . Figure 7(a) shows versus the total accretion stress for all our runs. The value of for each run is listed in Table 2. One can see a rough linear correlation between the velocity dispersion and the accretion stress (for reference, a linear fit is shown by the dashed line). However, detailed inspection shows that decreases more rapidly than as the dead zone increases in size. As found from Table 2, this is because the contribution from the atmosphere, , is insensitive to the size of the dead zone in the disk core. In Figure 7(b), we replot the data by replacing with the accretion stress in the disk core, . Comparison between Figures 7(a) and (b) shows that more tightly correlates with rather than with . We find that the data can be well fit by a simple linear relation

 ⟨δv2⟩mid=0.78αcorec2s, (21)

which is shown by the solid line in Figure 7(b). This result indicates that the accretion stress in the atmosphere does not contribute to the velocity fluctuation near the midplane.

Once is known, it is also possible to reproduce the vertical profile of the velocity dispersion. For the disk core (), we already know that is inversely proportional to the mean gas density and that hardly deviates from the initial Gaussian profile. From these facts, we can predict the vertical distribution of as

 ⟨δv2⟩ ≈ ⟨δv2⟩mid⟨ρ⟩mid⟨ρ⟩≈⟨δv2⟩midexp(z22h2) (22) ≈ 0.78αcorec2sexp(z22h2),

where Equation (21) has been used in the final equality. In Figure 8, we compare the vertical profiles of the random velocity directly obtained from runs Ideal, X1, and X1a with the predictions from Equation (22), where the values of are taken from Table 2. We see that Equation (22) successfully reproduces the vertical profiles of in the disk core. We remark that Equation (22) greatly overestimates the velocity dispersion at , where the gas density can no longer be approximated by the initial Gaussian profile (see Figure 4(a)).

#### 5.1.2 Density Fluctuation

Density fluctuations generated by MRI-driven turbulence gravitationally interact with planetesimals and larger solid bodies, affecting their collisional and orbital evolution in protoplanetary disks (Laughlin et al., 2004; Nelson & Papaloizou, 2004; Nelson & Gressel, 2010; Gressel et al., 2011). Here, we examine how the amplitude of the density fluctuations is determined the vertically integrated accretion stress.

As in Section 5.1.1, we begin with the analysis of the density fluctuations at the midplane, . We find from Table 2 that more tightly correlates with than with . Figure 9 shows versus for all runs. The best linear fit is given by

 ⟨δρ2⟩mid=0.47αcore⟨ρ⟩2mid, (23)

which is shown by the solid line in Figure 9. If we use this equation with Equation (21), we can also obtain the relation between the velocity dispersion and density fluctuation, . This is consistent with the idea that the fluctuations near the midplane are created by sound waves, for which (see also Section 4.1).

As shown in Section 4.1, is roughly proportional to along the vertical direction in the disk core. Hence, if is given, one can reconstruct the vertical profile of the density fluctuation in the disk core according to

 ⟨δρ2⟩ ≈ ⟨δρ2⟩mid⟨ρ⟩⟨ρ⟩mid≈⟨δρ2⟩midexp(−z22h2) (24) ≈ 0.47αcore⟨ρ⟩2midexp(−z22h2),

where we have used and Equation (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 MRI drives outgoing gas flow at the outer boundaries of the active layers. The MRI-driven outflow has been first observed by Suzuki & Inutsuka (2009) in shearing-box simulations and been recently demonstrated by Flock et al. (2011) in global simulations. Suzuki & Inutsuka (2009) and Suzuki et al. (2010) point out that this outflow might contribute to the dispersal of protoplanetary disks, although it is still unclear whether the outflow can really escape from the disks (see below). Meanwhile, MRI also contributes to the accretion of the gas in the radial direction. For consistent modeling of these two effects, we seek how the accretion stress and outflow flux are correlated with each other.

We evaluate the outflow mass flux in the following way. As seen in Section 4.1, the temporally and horizontally averaged vertical mass flux is nearly constant at heights . Using this fact, we define the outflow mass flux as the sum of averaged near upper and lower boundaries,

 ˙mw=∫hbndhbnd−h|⟨ρvz⟩|dzh+∫−hbnd+h−hbnd|⟨ρvz⟩|dzh, (25)

where is the height of the upper and lower boundaries of the simulation box.

Figure 10(a) shows normalized by versus for all simulations. The dimensionless quantity is equivalent to used in Suzuki et al. (2010). For model Ideal, the value is consistent with the ideal run of Suzuki & Inutsuka (2009). The dashed line shows the best linear fit . It can be seen that the linear fit captures a rough trend but still considerably overestimates the outflow flux for models Ideal and Y4. As seen in Table 2, these are the models in which dominates over . This implies that the turbulence in the disk core (which is the source of ) does not contribute to the outflow. In Figure 10(b), we replot the data by replacing with . We find that more tightly correlates with than with . The best linear fit is found to be

 ˙mw=0.016αatm⟨ρ⟩midcs. (26)

This result is consistent with the idea that the outflow is driven at the outer boundaries of the active layers (Suzuki & Inutsuka, 2009), because the dominant contribution to comes from heights very close to .

Although outflow from the simulation box is a general phenomenon in our simulations, it is unclear whether the outflow leaves or returns to the disk. In fact, the outflow velocity observed in our simulations does not exceed the sound speed even at the vertical boundaries. Since the escape velocity is higher than the sound speed, this means that the outflow does not have an outward velocity enough to escape out of the disk. Acceleration of the outflow beyond the escape velocity has not been directly demonstrated by previous simulations as well (Suzuki et al., 2010; Flock et al., 2011). However, Suzuki et al. (2010) point out a possibility that magnetocentrifugal forces and/or stellar winds could accelerate the outflow to the escape velocity. If the escape of the outflow will be confirmed in the future, our scaling formula for will certainly become a useful tool to discuss the dispersal of protoplanetary disks.

### 5.2. Saturation Predictors for the Accretion Stresses

In the previous subsection, we have shown that the amplitudes of various turbulent quantities scale with the vertically integrated stresses and . The next step is to find out how to predict and in the saturated state from the vertical magnetic flux (or equivalently ) and the resistivity profile . As shown in Section 4.2, the turbulent state of a disk depends on the resistivity only through the critical heights of the dead zone, and . Furthermore, the values of and are only weakly affected by the nonlinear evolution of MRI since the fluctuations in and are small inside the dead zone (see Section 4.1). Therefore, we expect that the effect of the resistivity can be well predicted by the values of and in the initial state, i.e., and . With this expectation, we try to derive saturation predictors for and as a function of , , and .

First, we focus on . Figure 11(a) plots versus for all our simulations. We see that scales roughly linearly with . The deviation from the linear scaling is expected to come from the difference in the dead zone size, i.e., and . In Figure 12, we plot the product as a function of . For models except Ideal and Y4, we find that is well predicted by a simple formula

 αcoreβz0=510exp(−0.54hres,0/h). (27)

Figure 11(b) replot the data in Figure 11(a) by replacing with . For models Ideal and Y4, Equation (27) underestimates . As explained in Section 4.2, these models exhibit higher magnetic activity near the midplane than the other models because of no or a thin dead zone . We expect that the higher magnetic activity gives additional contribution to . Taking into account this effect, we arrive at the final predictor function,

 αcore=510exp(−0.54hres,0/h)β−1z0+0.011exp(−3.6hΛ,0/h). (28)

Here, the numerical factors and appearing in the second term have been chosen to reproduce the results of runs Ideal and Y4, respectively. Figure 11(c) compares the final fitting formula with the numerical data. It can be seen that Equation (28) well predicts for all our models. Note that the second term of the predictor function is assumed to have no explicit linear dependence on unlike the first term. In fact, it is possible to reproduce our data by multiplying the second term by a prefactor . However, as we will see below, the absence of the prefactor makes the predictor function consistent with the results of ideal MHD simulations in the literature.

The predictor function for can be obtained in a similar way. Figure 13(a) shows versus for all our runs. We find that a simple linear relation well fits to the data except for models Ideal and Y4. This means that is characterized only by as long as the dead zone is thick . To take into account the cases of thin dead zones, we add a term proportional to as has been done for , and obtain

 αatm=530β−1z0+0.0043exp(−3.6hΛ,0/h), (29)

where the prefactor for the second term has been determined to fit to the result of run Ideal. As seen in Figure 13(b), Equation (29) well predicts the value of for all our runs.

Our predictor functions indicate that the vertically integrated accretion stress is inversely proportional to when a large dead zone is present (