The fate of planetesimals in turbulent disks. I

# The Fate of Planetesimals in Turbulent Disks with Dead Zones. I. The Turbulent Stirring Recipe

## Abstract

Turbulence in protoplanetary disks affects planet formation in many ways. While small dust particles are mainly affected by the aerodynamical coupling with turbulent gas velocity fields, planetesimals and larger bodies are more affected by gravitational interaction with gas density fluctuations. For the latter process, a number of numerical simulations have been performed in recent years, but a fully parameter-independent understanding has not been yet established. In this study, we present simple scaling relations for the planetesimal stirring rate in turbulence driven by magnetorotational instability (MRI), taking into account the stabilization of MRI due to Ohmic resistivity. We begin with order-of-magnitude estimates of the turbulence-induced gravitational force acting on solid bodies and associated diffusion coefficients for their orbital elements. We then test the predicted scaling relations using the results of recent Ohmic-resistive MHD simulations by Gressel et al. We find that these relations successfully explain the simulation results if we properly fix order-of-unity uncertainties within the estimates. We also update the saturation predictor for the density fluctuation amplitude in MRI-driven turbulence originally proposed by Okuzumi & Hirose. Combination of the scaling relations and saturation predictor allows to know how the turbulent stirring rate of planetesimals depends on disk parameters such as the gas column density, distance from the central star, vertical resistivity distribution, and net vertical magnetic flux. In Paper II, we apply our recipe to planetesimal accretion to discuss its viability in turbulent disks.

dust, extinction – magnetic fields – magnetohydrodynamics – planets and satellites: formation – protoplanetary disks – turbulence
6

## 1. Introduction

Planets are believed to form in circumstellar gas disks called protoplanetary disks. Planet formation begins with coagulation of submicron-sized dust particles through intermolecular forces. This stage is followed by the formation of kilometer-sized planetesimals through the gravitational collapse of microscopic dust aggregates mediated by gravitational (Goldreich & Ward, 1973) or streaming (Youdin & Goodman, 2005; Johansen et al., 2007) instabilities, or through further dust coagulation (e.g., Okuzumi et al., 2012; Windmark et al., 2012). Large planetesimals experience runaway growth mediated by gravitational focusing (Wetherill & Stewart, 1989; Kokubo & Ida, 1996), forming even larger solid bodies called protoplanets (or planetary embryos). In the final stage, protoplanets evolve into gas giants by accreting the disk gas or into larger solid planets through giant impacts during/after the dispersal of the gas disk.

The fate of these formation processes crucially depends on the turbulent state of the gas disk. Turbulence induces a random motion of solid particles smaller than planetesimals through the aerodynamical friction force (e.g., Völk et al., 1980; Ormel & Cuzzi, 2007). The resulting turbulent diffusion acts against accumulation of the solid particles (Carballido et al., 2005; Turner et al., 2006; Fromang & Papaloizou, 2006; Johansen et al., 2006), which limits planetesimal formation via gravitational instability (Youdin, 2011). The enhanced collision velocity may cause catastrophic disruption of the solid bodies, which inhibits direct collisional formation of planetesimals (Johansen et al., 2008; Okuzumi & Hirose, 2012). Turbulence also accumulates solid particles of particular sizes (e.g., Cuzzi et al., 2001; Johansen et al., 2007), but its relevance to planetesimal formation is under debate (Pan et al., 2011).

For planetesimals and larger solid bodies, stochastic gravitational forces induced by gas density fluctuations play a more important role. In turbulent disks, vorticity and nonlinear stress excite gas density fluctuations (Heinemann & Papaloizou, 2009a, b), which give rise to stochastic gravitational forces that act on solid bodies. This particularly affects the motion of large solid bodies that are well decoupled from the gas friction force, causing stochastic orbital migration (Laughlin et al., 2004; Nelson & Papaloizou, 2004; Nelson, 2005; Johnson et al., 2006; Oishi et al., 2007; Rein, 2012) and eccentricity stirring (Nelson, 2005; Ogihara et al., 2007; Ida et al., 2008).

The turbulence-induced eccentricity stirring severely constrains the formation of protoplanets at a fundamental level. In order for gravitational runaway growth to set in, the velocity dispersion of planetesimals must be smaller than their escape velocity (Wetherill & Stewart, 1989). However, in a fully turbulent disk, this requirement is unlikely to be satisfied for planetesimals smaller than 100 km in size (Ida et al., 2008; Nelson & Gressel, 2010). This indicates that runaway growth could be significantly delayed depending on the turbulent state of the disks (Nelson, 2005; Ormel et al., 2010). Moreover, the high turbulence-driven relative velocity can make a collision between planetesimals disruptive rather than accumulative, especially in outer regions of the disks (Nelson, 2005; Ida et al., 2008; Nelson & Gressel, 2010; Yang et al., 2009, 2012). Thus, to understand the fate of planetesimal growth and succeeding planet formation, it is essential to know how gas turbulence is driven in protoplanetary disks, and how its strength depends on the disk environment.

One mechanism that can drive strong turbulence is the magnetorotational instability (MRI; Balbus & Hawley, 1991). This is an MHD instability resulting from the coupling between a differentially rotating gas disk and magnetic fields. In an ideal case where the coupling is strong enough, the MRI drives strong gas turbulence with the Shakura–Sunyaev parameter or even larger depending on the strength of the net vertical magnetic fields (e.g., Davis et al., 2010; Suzuki et al., 2010). However, because the ionization degree of protoplanetary disks is generally very low, non-ideal MHD effects strongly affect the actual level of the turbulence. For example, a high Ohmic resistivity near the disk midplane prevents the coupling between the gas and magnetic fields and thereby creates a “dead zone” where MRI is inactive (Gammie, 1996). The size of the dead zone depends on the ionization degree of the disk gas, and is generally large when tiny dust particles that efficiently capture ionized gas particles are abundant (e.g., Sano et al., 2000; Ilgner & Nelson, 2006). Ambipolar diffusion has a similar effect on MRI, but at higher altitudes where the gas density is low (Bai, 2011; Perez-Becker & Chiang, 2011a, b; Mohanty et al., 2013; Dzyurkevich et al., 2013).

Recently, Gressel et al. (2011) first studied the effect of an Ohmic dead zone on turbulent planetesimal stirring. They performed local stratified MHD simulations at 5 AU taking into account a high Ohmic resistivity provided by abundant small dust particles. They showed that the resulting large dead zone considerably suppresses the planetesimal stirring rate. The effect has been more extensively studied in their latest paper (Gressel et al. 2012; henceforth GNT12) for various values of the net vertical magnetic flux. They concluded that planetesimal growth beyond the disruption barrier is possible in a dead zone if the net flux is so weak that upper MRI-active layers do not generate strong density waves. This indicates that a dead zone can provide a safe haven for planetesimals.

However, there still remain two open issues. First, how much dust is needed to maintain a large enough dead zone? Gressel et al. (2011, 2012) fixed the amount of -sized dust particles to be in mass of the total solids in the disk. However, it is unclear whether this amount is reasonable in late stages of planet formation where a significant fraction of solids in the disk should have been incorporated into planetesimals. In principle, tiny particles can be resupplied when planetesimals undergo collisional fragmentation or erosion. However, such tiny particles are usually removed immediately through their mutual sticking and/or sweep up by larger dust particles. Thus, the amount of residual dust is determined by the balance between these competing processes, and therefore cannot be determined a priori. Second, can a dead zone act as a safe haven at every location in protoplanetary disks? The results of Gressel et al. (2011, 2012) only apply to 5 AU from the central star, but turbulent planetesimal stirring is generally more effective further out in disks (Ida et al., 2008). In order to study whether the dead zone is beneficial for planetesimal growth in general circumstances, a model that does not rely on a specific choice of disk parameters is desirable.

The aim of this study is to provide a general recipe for planetesimal stirring in MRI-driven turbulence. We construct scaling relations that clarify how the turbulent quantities relevant to planetesimal stirring depend on each other and on basic disk parameters. This is an extension of recent work by Okuzumi & Hirose (2011, henceforth OH11). They performed a systematic set of local stratified MHD simulations with a dead zone, and provided an analytic prescription for the amplitude of the gas density fluctuations as a function of the net vertical flux, vertical resistivity profile, and other disk parameters. In this paper, we begin with an order-of-magnitude estimate to derive scaling relations that link the density fluctuation amplitude to the turbulent stirring rate of solid bodies. We then calibrate them using the published data by GNT12. We also update the density fluctuation recipe of OH11 using the same published data. An application of our recipe to runaway planetesimal growth will be presented in Paper II (Ormel & Okuzumi, 2013).

The plan of this paper is as follows. In Section 2, we present order-of-magnitude estimates that predict relationships among the density fluctuation, random gravity, and orbital diffusion coefficients for planetesimals. In Section 3, we compare our predictions with the simulation results presented by GNT12 to present calibrated prescriptions for planetesimal stirring. The predictor function for the density fluctuations is presented in Section 4. Comparison with previous results relying on ideal MHD and implication for planetesimal stirring in protoplanetary disks is given in Section 5. A summary of this study is given in Section 6.

## 2. Order-of-Magnitude Estimates

In order to clarify how the turbulent stirring rate of planetesimals generally depends on disk parameters, we begin with deriving scaling relations between relevant turbulent quantities from order-of-magnitude arguments. Verification and calibration of the derived relations will be done in Section 3.

Our estimation follows two steps. First, we relate gas density fluctuations to random gravitational forces on planetesimals using Gauss’s law for gravity. We then relate the random gravity to the diffusion coefficients for planetesimals. The second step is based on recent work by Rein & Papaloizou (2009) that regards the equation of motion for planetesimals as a stochastic differential (or Langevin) equation.

### 2.1. Random Gravity

We denote the gas density perturbation by and the induced gravitational force on planetesimals (per unit mass) by (see Figure 1 (a)). These are assumed to be stochastic variables with vanishing mean values and nonzero mean square values and . The density perturbation and induced gravitational force are related to each other by Gauss’s law for gravity,

 ∇⋅F=−4πGδρ, (1)

where is the gravitational constant.

We want to estimate the amplitude of the random gravity for given density fluctuation amplitude. This can be done by assuming that the density fluctuations have a characteristic wavenumber . With this assumption, we can estimate that , where is the magnitude of . Thus, from Equation (1), we have

 ⟨F2⟩1/2∼4πGk⟨δρ2⟩1/2. (2)

If the disk is entirely MRI-turbulent, the characteristic wavenumber (or the inverse of the correlation length) is of the order , where is the gas scale height (Guan et al., 2009; Heinemann & Papaloizou, 2009a, b; Nelson & Gressel, 2010). This suggests that

 ⟨F2⟩1/2=A1GH⟨δρ2⟩1/2, (3)

where is an order-of-unity number that represents the overall uncertainty in the above estimate.

However, if there is a dead zone at the midplane, we need to take into account the shearing-out of the density fluctuations. In the presence of a dead zone, the sources of the density fluctuation at the midplane are density waves that have propagated from the upper MRI-active layers. At the midplane, these waves have a higher than they had in the active layers because the differentially rotating background flow shears them out during their propagation (see Figure 1 (b)). The importance of the shearing-out has first been pointed out by GNT12 (see their Section 5.1), and we will quantify this with the following argument. Let us denote the radial and azimuthal wavenumbers of a density wave by and , respectively. The shearing motion of the gas disks changes the radial wavenumber, and this can be expressed as (Goldreich & Lynden-Bell, 1965)

 kr=kr0+32kϕΩδttravel. (4)

where is the Keplerian frequency, is the time passed after the wave is generated, and is the initial value of . We now assume that a dead zone has a vertical extent , where is the dead-zone half width. Then, for density waves at the midplane, should be comparable to the time the waves travel from the active layer to the midplane, i.e., , where is the sound speed. This can be rewritten as since . Substituting this into Equation (4) and assuming and , the characteristic wavenumber of the density waves at the midplane can be estimated as

 k∼1H(1+A2HDZH), (5)

where is another order-of-unity constant. The assumption may be somewhat inaccurate given the anisotropy of MRI-driven turbulence, but its effect is absorbed in . If we use this in Equation (2) we obtain

 ⟨F2⟩1/2=A1GH1+A2HDZ/H⟨δρ2⟩1/2. (6)

Equation (6) also applies to ideal MRI-turbulent disks because it reduces to Equation (3) in the limit of .

Equation (6) indicates that in the presence of a dead zone no simple linear scaling applies to the relation between the magnitudes of the random gravity fields and density fluctuations, as observed by GNT12. A dead zone reduces the random forcing in two ways: directly by suppressing , and indirectly by enhancing (shearing-out). The factor appearing in Equation (6) expresses the second effect.

### 2.2. Diffusion Coefficients

In fluctuating gravity fields, the motion of solid bodies can be described as a random walk in phase space. Rein & Papaloizou (2009) formulated this process by treating the equation of motion for the bodies as a Langevin equation with stochastic forcing. To describe the results of Rein & Papaloizou (2009), we denote the changes in the semi-major axis and eccentricity of a body during a time interval as and , respectively. A random walk in phase space means that the ensemble averages and grow linearly with on a timescale much longer than the correlation time of the fluctuating gravity. This process can be characterized by constant diffusion coefficients

 Da≡12⟨(Δa)2⟩Δt, (7)
 De≡12⟨(Δe)2⟩Δt. (8)

Rein & Papaloizou (2009) derived the exact expressions of and under stochastic force . These read

 ⟨(Δa)2⟩=8⟨F2ϕ⟩τcΩ2Δt, (9)
 ⟨(Δe)2⟩=(⟨F2r⟩+4⟨FrFϕ⟩+4⟨F2ϕ⟩)τc(1+Ω2τ2c)a2Ω2Δt, (10)

where and are the mean squared amplitudes of the radial and azimuthal components of , respectively (see Equations (46) and (47) of Rein & Papaloizou 20097). The corresponding diffusion coefficients are

 Da=4⟨F2ϕ⟩τcΩ2, (11)
 De=(⟨F2r⟩+4⟨FrFϕ⟩+4⟨F2ϕ⟩)τc2(1+Ω2τ2c)a2Ω2. (12)

For MRI-driven turbulence, is suggested by a number of simulations (e.g., Sano et al., 2004; Gressel et al., 2011, 2012). It is also likely that within an accuracy of order unity. Therefore, we anticipate that the diffusion coefficients are of the forms

 Da=AaΩ3⟨F2ϕ⟩, (13)
 De=Aea2Ω3⟨F2ϕ⟩, (14)

where and are order-of-unity constants. Note that the uncertainties about the anisotropy of the random force, and , are absorbed in these constants.

## 3. Calibration with GNT12 Data

In the previous section, we have predicted how turbulent quantities relevant to the orbital evolution of planetesimals should be related to each other. Here, we test these predictions using the published data of MHD simulations by GNT12. Our goal is to determine the order-of-unity constants involved in the predicted relationships (, , , and ; see Equations (6), (13), and (14)).

GNT12 conducted local stratified resistive MHD simulations at 5 AU from the central star with different sets of the disk mass, ionization strength, and net vertical field strength. The gas surface density was given by with (= 1, 2, or 4) being a dimensionless factor. The stellar mass and disk aspect ratio are fixed to and , respectively. The Ohmic resistivity was calculated from the balance between external ionization and recombination in the gas phase and on dust grains. The ionization rate was the sum of the contributions from cosmic rays (), X-rays (), and short-lived radionuclides (), with and being chosen as 10 and (= 1 or 20) times the standard values, respectively. The dust-to-gas mass ratio was fixed to , and the size of the dust grains was chosen to be . The results of all these simulations are summarized in Table 3 of GNT12. We will use these data to test and calibrate our order-of-magnitude relations.

One of the most important parameter is the net vertical magnetic flux . This is a conserved quantity in a local-box simulation, and determines the strength of MRI turbulence in the saturated state (e.g., Hawley et al., 1995; Sano et al., 2004; Suzuki et al., 2010; Okuzumi & Hirose, 2011). In the GNT12 simulations, was chosen in the range 2.68–46 mG. GNT12 also conducted a run with a higher net flux (), but the run resulted in an essentially laminar final state with no sustained MRI turbulence.

Table 1 lists the values of , and for all simulations presented by GNT12. Run A1 and B1 are ideal MHD simulations while runs labeled by ‘D’ include Ohmic diffusion.

### 3.1. Characterization of the Vertical Structure

We have predicted in Section 2.1 that the relation between the density fluctuation amplitude and associated random force depends on the vertical extent of a dead zone. In order to verify this, we need to define the dead zone in advance and in a way that does not depend on any specific choice of the model parameters. In this study, we follow OH11 and define dead and active zones in terms of the linear perturbation theory of MRI with Ohmic resistivity.

GNT12 calculated the Ohmic resistivity using the charge reaction network model4 of Ilgner & Nelson (2006). The network consists of free electrons, two species of ions ( and ) and charged dust grains. We reproduce the resistivity using the analytic prescription presented by Okuzumi (2009) with the assumption that dominates the ions at all heights.10 To test our calculation, in the upper panel of Figure 2, we plot the vertical profiles of the magnetic Reynolds number for models D1, D1.2, and D1.4 of GNT12. Here, the disk is assumed to be vertically in hydrostatic equilibrium, and the the vertical profile of the disk gas is given by , where is the midplane gas density. Comparing our plot with Figure 1 of GNT12, we confirm that our calculation successfully reproduces their resistivity.

With the information of , we characterize the vertical turbulent structure of the disk using the four-layer description proposed by Okuzumi & Hirose (2011). The characterization is based on the linear stability analysis of MRI in vertically stratified disks in the presence of Ohmic resistivity (Jin, 1996; Sano & Miyama, 1999). The linear analysis shows that at each height the gas motion is unstable if the wavelength of the most unstable local MRI mode is longer than the gas scale height . In the presence of Ohmic resistivity, the most unstable wavelength can be approximately given by , where

 λideal(z)=2πvAz(z)Ω (15)

and

 λres(z)=2πη(z)vAz(z) (16)

are the characteristic wavelengths of the unstable modes in the ideal and resistive limits, respectively, with being the vertical component of the Alfvén velocity. The most unstable wavelength can be alternatively written as , where

 Λ≡v2AzηΩ=λidealλres (17)

is the so-called Elsasser number. The growth rate of the local MRI modes is approximately given by . The instability is strong () when , weak () when , and absent when . Thus, the vertical distributions of and predict at which height the MRI is unstable.

When there is a nonzero net vertical magnetic field , it is useful to evaluate and assuming that the disk is in the laminar state, i.e., assuming that and at all heights. We denote them by and , respectively. As an example, the lower panel of Figure 2 plots as a function of for models D1, D1.2, and D1.4. Note that and are increasing and decreasing functions of , respectively, because decreases with while increases with . For this reason, the region where the MRI is unstable (i.e., ) is bounded from both below and above. In this paper, we will refer to such regions as the active layers. For models D1, D1.2, and D1.4, the active layers are located at , and , respectively.

Based on the stability criterion outlined above, OH11 introduced three critical heights , , and defined by

 λideal(z=Hideal)=H, (18)
 Λ(z=HΛ)=1, (19)
 λres(z=Hres)=H, (20)

respectively. The active layer defined by has a vertical extent . MRI is stable in the magnetically dominated atmosphere at , and in the high- region at . In this paper, we define a dead zone as the region , i.e., . One can subdivide the active layer into two sublayers and , at which MRI operates strongly and weakly, respectively.

As we did for and , we define the critical heights in the laminar state by , , and . For , there is an analytic expression (Equation (14) of OH11)

 Hideal,0=[2ln(βz08π2)]1/2H, (21)

where is the midplane plasma beta measured with the net field strength . Table 1 list the values of , , and for all the GNT12 simulations. In the lower panel of Figure 2, the dark and light horizontal bars indicate the ideal and resistive MRI regions, and , respectively, for models D1, D1.2, and D1.4.

In principle, the critical heights in a turbulent state can differ from those in the laminar state since and depend on and . One can see this by comparing the critical heights in the initial (laminar) and time-averaged (turbulent) states measured in the OH11 simulations (see Tables 1 and 2 of OH11). The difference is the largest for due to strong fluctuating magnetic fields at the top of active layers, and we will discuss this in more detail in Section 4.2. In contrast, the difference is much smaller for , and is negligible for , because magnetic activity is weak at the boarder of dead/active regions.

### 3.2. Random Gravitational Force vs Density Fluctuation

To test Equation (6), we compare the rms values of the random gravitational force and density fluctuations measured by GNT12. Table 1 lists the rms amplitudes of the density fluctuations at the midplane () and azimuthal gravitational force acting on test particles () measured in the GNT12 simulations. These values are taken from Table 3 of GNT12 where they are given in terms of the relative fluctuation and the rms torque , respectively.

Figure 3(a) shows versus for all the available data. The dashed and dotted lines in Figure 3 show linear scalings (see Equation (3))

 ⟨F2ϕ⟩1/2=A1GH⟨δρ2⟩1/2mid (22)

with and 0.16, respectively. It is clearly seen that no linear scaling can explain the whole data. GNT12 pointed out this using the data for D1 runs (open circles and crosses in Figure 3). We find that this can be seen more clearly by adding the data for runs B1 (filled circle) and D2 (filled square), for which the dead zone is absent and smaller than that in the D1 runs, respectively. This fact strengthens the idea that shearing-out of density waves causes suppression of the random gravity forces.

From the analysis in Section 2.1, we expect that the effect of the shearing-out can be extracted by taking the ratio between and and comparing it with the half width of the dead zone. As we stated in Section 3.1, we measure the dead zone half width with the critical height defined by Equation (20). Specifically, we here use the value in the laminar state, , so that we can calculate it directly from the initial conditions (as we noted in Section 3.1, the value of is very insensitive to the presence or absence of turbulence). Figure 4 plots the ratio versus . We see a decreasing trend in with increasing , which is consistent with the idea that the random force is weakened as the density waves travel from the active layer to the midplane (see Section 2.1). Following Equation (6), we fit the data shown in Figure 4 with a function of the form

 ⟨F2ϕ⟩1/2⟨δρ2⟩1/2mid=A1GH1+A2Hres,0/H, (23)

where and are the fitting parameters. In Figure 4, the dotted, solid, and dashed curves show Equation (23) with with , 4.5, and 7.0, respectively. We find that the set best reproduces the relation between and , within an accuracy of factor 2.

Thus, we have found that the relation between and can be well represented by

 ⟨F2ϕ⟩1/2=1.6GH1+4.5Hres,0/H⟨δρ2⟩1/2mid, (24)

which is also shown in the lower panel of Figure 3.

### 3.3. Diffusion Coefficients vs Random Force Amplitude

The next step is to verify the linear scaling between the diffusion coefficients and as predicted by Equations (13) and (14). GNT12 measured the change in the semimajor axis and eccentricity of particles without gas friction and with initial eccentricity . They showed that the ensemble averages of and grow linearly with time , indicating a random walk of the particles’ motion in the phase space. GNT12 expressed the time evolution of and in the forms

 ⟨(Δa)2⟩1/2=Cσ(Δx)H(ΩΔt2π)1/2, (25)
 ⟨e2⟩1/2=Cσ(e)Ha(ΩΔt2π)1/2, (26)

respectively, where and are dimensionless coefficients that depend on the adopted disk model. These values are listed in Table 3 of GNT12.

The diffusion coefficients and can be read off from the values of and . For , we have

 Da=(Cσ(Δx))24πH2Ω, (27)

which directly follows from Equations (7) and (25). For , we need to convert to the eccentricity displacement for general (nonzero) initial eccentricity. As shown by Yang et al. (2009), this conversion is given by . Hence, from Equations (8) and (26), we get

 De=2(Cσ(e))24π(4−π)H2Ωa2. (28)

The values of and for all the available data are listed in Table 1. The listed values are normalized by and , respectively, which are the natural units for these diffusion coefficients in local simulations.

Now we calibrate Equations (13) and (14) using the available data. Figure 5 plots and versus the rms azimuthal gravitational force measured in the GNT12 simulations. We clearly see the trend as predicted by Equations (13) and (14). We determine the dimensionless parameters and so that the maximum logarithmic error between the data and predictions from each of the equations is minimized. We find that the best-fit parameters are and . The best-fit relations are shown in Figure 5 by the solid lines.

To summarize, we have found that the diffusion coefficients for the semi-major axis and eccentricity of solid bodies scale with the mean squared amplitude of fluctuating gravity fields as

 Da=4.0⟨F2ϕ⟩Ω3, (29)
 De=1.7⟨F2ϕ⟩a2Ω3. (30)

If we eliminate using Equation (23), these scaling relations can be rewritten as a function of and . These read

 Da=10(1+4.5Hres,0/H)2(a2H⟨δρ2⟩1/2midM∗)2a2Ω, (31)
 De=4.4(1+4.5Hres,0/H)2(a2H⟨δρ2⟩1/2midM∗)2Ω, (32)

where we have used that .

## 4. Predicting Density Fluctuation Amplitudes

Equations (31) and (32) tell us how the diffusion coefficients and are related to the amplitude of the gas density fluctuations, . To predict the values of and for given disk parameters, we need to know how depends on these parameters.

OH11 provided a simple analytic formula (called the “saturation predictor”) that predicts as a function of key disk parameters. In this section, we test whether the formula accurately predicts the density fluctuation amplitude observed in the GNT12 simulations.

### 4.1. The OH11 Predictor versus GNT12 Data

 [⟨δρ2⟩1/2mid]OH11=√0.47αcoreρmid, (33)

where is a dimensionless coefficient that is proportional to the turbulent accretion stress integrated over height (for details, see OH11). OH11 provided an empirical formula that relates to key disk parameters such as the net vertical flux and the vertical distribution of the resistivity . The formula reads

 αcore=510βz0exp(−0.54Hres,0H)+0.011exp(−3.6HΛ,0H). (34)

In fact, the numerical prefactor appearing in Equation (33) weakly depends on the numerical resolution adopted in simulations; below we will refine the prefactor in accordance with the data of GNT12 simulations that adopted a higher resolution.

In Figure 6(a), we compare the measured values of for the GNT12 simulations with the prediction from Equation (33). We see that the OH11 predictor reasonably reproduces most of the observed values, especially for runs D1, D1., D1-WF, and D1-NVFa. However, detailed inspection shows that these observed values are higher than the prediction by (see the dotted line in Figure 6(a)). This discrepancy can be attributed to the fact that GNT12 adopted a higher numerical resolution than OH11. As shown by both OH11 and GNT12, the density fluctuation amplitude increases slowly with improving the numerical resolution (see Figure 18 of OH11 and Figure A1 of GNT12). We expect that the values measured by GNT12 are well converged to the true values since an even higher resolution does not give any significant change in the density fluctuation amplitude (see Appendix A of GNT12).

A more important discrepancy can be found for runs D1-NVF. In these runs, the net vertical field strength was gradually increased from 2.7 mG to 86.0 mG. The results show that initially increases with but gets suppressed at . This can be seen in Figure 7 of GNT12, and we also show this in the lower panel of our Figure 7. The OH11 predictor does not reproduce this suppression (as shown by the gray dashed line in Figure 7) and consequently overestimates the density fluctuation amplitude for and 43.0 mG (marked by the open squares in Figure 6).

### 4.2. Why are the Density Fluctuations Suppressed at High ⟨Bz⟩?

GNT12 explained the suppressed at high as a consequence of narrowed MRI-active layers. In a stratified disk, an MRI-active layer is bounded from above by a magnetically dominated atmosphere. As the magnetic fields become stronger, the base of the atmosphere moves down to the midplane, narrowing the active layer beneath. In the limit of high fields, the atmosphere will erode most of the active layer, and will in turn suppress the generation of density fluctuations.

To confirm the hypothesis raised by GNT12, we estimate the vertical extent of the active layer using the critical heights we introduced in Section 3.1. As explained there, local MRI modes can exist only at , above which exceeds the gas scale height . The region where local MRI modes exist is also bounded from below by , below which the Ohmic resistivity stabilizes all MRI modes (i.e., ). Thus, we may measure the vertical width of an MRI-active layer as , with corresponding to the base of the magnetically dominated atmosphere while to the the active/dead zone interface.

The active layer width defined above naturally explains the suppression of the density fluctuation amplitude at high net fields. In the upper panel of Figure 7, we plot the critical heights in the laminar state, and , as a function of . We see that the active layer width in the laminar state is already as small as for . Furthermore, vanishes at , indicating that the magnetically dominated atmosphere completely suppresses the MRI-active layer for this value of or larger. This exactly explains what happened in run D1-NVFb, in which the disk returned to a laminar state at (see Figure 11 of GNT12).

A more quantitative analysis can be made by noting that further decreases as MRI-driven turbulence develops. Simulations by OH11 show that measured in the fully turbulent state is smaller than that in the laminar state () by about one scale height (see their Tables 1 and 2). By contrast, is nearly independent of the turbulence state in the active layer because turbulence is always weak on the dead/active zone boundary. The definition of (Equation (18)) allows us to directly calculate how much decreases as the turbulence develops at the top of the active layer. If we measure with the rms amplitude of the fields, , then Equation (18) can be rewritten as

 ρidealc2s=π⟨B2z⟩ideal (35)

where and are the values of and at , respectively. We may approximate with the hydrostatic density profile, , since the gas pressure dominates over the magnetic pressure at (see also Figure 4(a) of OH11). Thus, the definition of can be further rewritten as

 H2ideal = 2H2ln(ρmidc2sπ⟨B2z⟩ideal) (36) = H2ideal,0−2H2ln(1+⟨δB2z⟩ideal⟨Bz⟩2),

where is the value of in the laminar state (Equation (21)) and is the mean squared amplitude of the fluctuating (turbulent) -fields at . Equation (36) demonstrates that decreases as turbulence grows.

It is useful to know to what extent decreases if turbulence is fully developed. The results of the OH11 simulations show that in the fully turbulent state, approximately satisfies the relation (see Appendix A)

 ⟨δB2z⟩ideal∼30⟨Bz⟩2. (37)

Inserting this relation into Equation (36), we find that should decrease to , where

 H2ideal,∞ ≈ H2ideal,0−7H2. (38)

In the upper panel of Figure 7, the solid curve shows versus for model D1-NVF. The values of for all the GNT12 simulations are listed in Table 1. Since is hardly affected by turbulence, we can estimate the width of the active layer in the fully turbulent state as . As seen in the figure, becomes negative for . This suggests that turbulence cannot be fully developed in run D1-NVFb with and 43.0 mG; if turbulence were fully developed, then the magnetically dominated atmosphere would completely suppress the active layer.

The above analysis confirms the hypothesis by GNT12 that the magnetically dominated atmosphere limits the saturation level of turbulence in the active layer. We summarize the mechanisms of this effect in Figure 8. This effect is not taken into account in the OH11 saturation predictor (Equation (33)) as was larger than for all the OH11 simulations. This explains why the OH11 predictor overestimates the amplitude of the density fluctuations at high .

### 4.3. Refining the OH11 Predictor with the “Saturation Limiter”

Based on the above consideration, we construct a toy model that accounts for the suppression of the density fluctuation amplitude at high .

The model is based on two assumptions. Firstly, we assume that if then turbulence grows until reaches . This means that the saturated value of is given by Equation (36) with . Solving the equation with respect to , we get

 ⟨δB2z⟩ideal=[exp(H2ideal,0−H2res,02H2)−1]⟨Bz⟩2. (39)

Combining this with Equation (37), we obtain the saturation predictor for for general cases,

 ⟨δB2z⟩ideal=30L⟨Bz⟩2, (40)

where is defined by

 L=min{1,130[exp(H2ideal,0−H2res,02H2)−1]}. (41)

This expresses that turbulence in the active layer is limited at a low level () when (see also Figure 8). We will call the “saturation limiter.”

Secondly, we assume that is proportional to ; namely, if is suppressed by factor ,