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 parameterindependent 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 orderofmagnitude estimates of the turbulenceinduced 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 ohmicresistive MHD simulations by Gressel et al. We find that these relations successfully explain the simulation results if we properly fix orderofunity uncertainties within the estimates. We also update the saturation predictor for the density fluctuation amplitude in MRIdriven turbulence originally proposed by Okuzumi amp Hirose. Combination of the scaling relations and saturation predictor allows us 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: MHD', 'planets and satellites: formation', 'protoplanetary disks', 'turbulence', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Solar and Stellar Astrophysics']
The Fate of Planetesimals in Turbulent Disks with Dead Zones. I. The Turbulent Stirring Recipe
{'THE FATE OF PLANETESIMALS IN TURBULENT DISKS WITH DEAD ZONES. I. THE TURBULENT STIRRING RECIPE': 'S atoshi O kuzumi 1,2,3 and C hris W. O rmel 4,5  \nApJ accepted', 'ABSTRACT': 'Turbulence in protoplanetary disks a ff ects planet formation in many ways. While small dust particles are mainly a ff ected by the aerodynamical coupling with turbulent gas velocity fields, planetesimals and larger bodies are more a ff ected 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 di ff usion coe ffi cients 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.  \nKeywords: dust, extinction - magnetic fields - magnetohydrodynamics - planets and satellites: formation protoplanetary disks - turbulence', '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.  \nThe 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., Volk et al. 1980; Ormel & Cuzzi 2007). The resulting turbulent di ff usion 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  \n- 1 Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro-ku, Tokyo, 152-8551; okuzumi@geo.titech.ac.jp\n- 2 Department of Physics, Nagoya University, Nagoya, Aichi 464-8602, Japan\n- 3 JSPS Superlative Research Fellow\n- 4 Astronomy Department, University of California, Berkeley, CA 94720, USA  \n(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).  \nFor 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 a ff ects 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).  \nThe 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 plan-  \netesimal 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.  \nOne 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 di ff erentially 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 α ∼ 10 -2 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 e ff ects strongly a ff ect 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 e ffi ciently capture ionized gas particles are abundant (e.g., Sano et al. 2000; Ilgner & Nelson 2006). Ambipolar di ff usion has a similar e ff ect 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).  \nRecently, Gressel et al. (2011) first studied the e ff ect 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 e ff ect 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.  \nHowever, 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 0 . 1 µ m-sized dust particles to be 10 % 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 e ff ective 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.  \nThe 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  \nby 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-ofmagnitude 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).  \nThe 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 di ff usion coe ffi cients 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.  \nOur 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 di ff usion coe ffi cients 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 di ff erential (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 F (see Figure 1 (a)). These are assumed to be stochastic variables with vanishing mean values 〈 δρ 〉 = 〈 F 〉 = 0 and nonzero mean square values 〈 δρ 2 〉 and 〈 F 2 〉 ≡ 〈 F 2 〉 . The density perturbation and induced gravitational force are related to each other by Gauss's law for gravity,  \n∇· F = -4 π G δρ, (1)  \nwhere G is the gravitational constant.  \nWe 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 k . With this assumption, we can estimate that |∇ · F | ∼ k 〈 F 2 〉 1 / 2 , where k = | k | is the magnitude of k . Thus, from Equation (1), we have  \n〈 F 2 〉 1 / 2 ∼ 4 π G k 〈 δρ 2 〉 1 / 2 . (2)  \nIf the disk is entirely MRI-turbulent, the characteristic wavenumber (or the inverse of the correlation length) k is of the order ∼ 1 / H , where H is the gas scale height (Guan et al.  \nFigure 1. Schematic illustration showing how density fluctuations create fluctuating gravity fields. The patches indicate density bumps produced by turbulence, while the circle at the center indicates a solid body gravitationally interacting with the density fluctuations. In MRI-active regions, the characteristic wavenumber k 0 of the density fluctuations are of the order ∼ H (panel (a)). In dead zones, density fluctuations have a higher wavenumber because they get sheared out by the background flow when they propagate from active regions (panel (b)).  \n<!-- image -->  \n<!-- image -->  \n(b) w/ shearing-out  \n2009; Heinemann & Papaloizou 2009a,b; Nelson & Gressel 2010). This suggests that  \n〈 F 2 〉 1 / 2 = A 1 GH 〈 δρ 2 〉 1 / 2 , (3)  \nwhere A 1 is an order-of-unity number that represents the overall uncertainty in the above estimate.  \nHowever, 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 k than they had in the active layers because the di ff erentially 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 kr and k φ , respectively. The shearing motion of the gas disks changes the radial wavenumber, and this can be expressed as (Goldreich & Lynden-Bell 1965)  \nkr = kr 0 + 3 2 k φ Ω δ t travel . (4)  \nwhere Ω is the Keplerian frequency, δ t travel is the time passed after the wave is generated, and kr 0 is the initial value of kr . We now assume that a dead zone has a vertical extent | z | ≤ H DZ, where H DZ is the dead-zone half width. Then, for density waves at the midplane, δ t travel should be comparable to the time the waves travel from the active layer to the midplane, i.e., δ t travel ∼ H DZ / cs , where cs is the sound speed. This can be rewritten as Ω δ t travel ∼ H DZ / H since H = cs / Ω . Substituting this into Equation (4) and assuming kr 0 ∼ k φ ∼ 1 / H and k ∼ kr , the characteristic wavenumber k of the density waves at the midplane can be estimated as  \nk ∼ 1 H ( 1 + A 2 H DZ H ) , (5)  \nwhere A 2 is another order-of-unity constant. The assumption kr 0 ∼ k φ may be somewhat inaccurate given the anisotropy of MRI-driven turbulence, but its e ff ect is absorbed in A 2. If we use this in Equation (2) we obtain  \n〈 F 2 〉 1 / 2 = A 1 GH 1 + A 2 H DZ / H 〈 δρ 2 〉 1 / 2 . (6)  \nEquation (6) also applies to ideal MRI-turbulent disks because it reduces to Equation (3) in the limit of H DZ → 0.  \nEquation (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 〈 δρ 2 〉 1 / 2 , and indirectly by enhancing k (shearing-out). The factor (1 + A 2 H DZ / H ) -1 appearing in Equation (6) expresses the second e ff ect.", '2.2. Di ff usion Coe ffi cients': '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 a and eccentricity e of a body during a time interval ∆ t as ∆ a and ∆ e , respectively. A random walk in phase space means that the ensemble averages 〈 ( ∆ a ) 2 〉 and 〈 ( ∆ e ) 2 〉 grow linearly with ∆ t on a timescale much longer than the correlation time τ c of the fluctuating gravity. This process can be characterized by constant di ff usion coe ffi cients  \nDa ≡ 1 2 〈 ( ∆ a ) 2 〉 ∆ t , (7)  \nDe ≡ 1 2 〈 ( ∆ e ) 2 〉 ∆ t . (8)  \nRein & Papaloizou (2009) derived the exact expressions of 〈 ( ∆ a ) 2 〉 and 〈 ( ∆ e ) 2 〉 under stochastic force F . These read  \n〈 ( ∆ a ) 2 〉 = 8 〈 F 2 φ 〉 τ c Ω 2 ∆ t , (9)  \n〈 ( ∆ e ) 2 〉 = ( 〈 F 2 r 〉 + 4 〈 FrF φ 〉 + 4 〈 F 2 φ 〉 ) τ c (1 + Ω 2 τ 2 c ) a 2 Ω 2 ∆ t , (10)  \nwhere 〈 F 2 r 〉 and 〈 F 2 φ 〉 are the mean squared amplitudes of the radial and azimuthal components of F , respectively (see Equations (46) and (47) of Rein & Papaloizou 2009 6 ). The  \n6 Note that Equations (46) and (47) of Rein & Papaloizou (2009) implicitly assume 〈 F 2 r 〉 = 〈 F 2 φ 〉 and 〈 FrF φ 〉 = 0 while our Equation (9) and (10) allow  \ncorresponding di ff usion coe ffi cients are  \nDa = 4 〈 F 2 φ 〉 τ c Ω 2 , (11)  \nDe = ( 〈 F 2 r 〉 + 4 〈 FrF φ 〉 + 4 〈 F 2 φ 〉 ) τ c 2(1 + Ω 2 τ 2 c ) a 2 Ω 2 . (12)  \nFor MRI-driven turbulence, τ c ∼ Ω -1 is suggested by a number of simulations (e.g., Sano et al. 2004; Gressel et al. 2011, 2012). It is also likely that 〈 F 2 r 〉 ∼ 〈 FrF φ 〉 ∼ 〈 F 2 φ 〉 within an accuracy of order unity. Therefore, we anticipate that the di ff usion coe ffi cients are of the forms  \nDa = A a Ω 3 〈 F 2 φ 〉 , (13)  \nDe = A e a 2 Ω 3 〈 F 2 φ 〉 , (14)  \nwhere A a and A e are order-of-unity constants. Note that the uncertainties about the anisotropy of the random force, 〈 F 2 r 〉 / 〈 FrF φ 〉 and 〈 F 2 r 〉 / 〈 F 2 φ 〉 , 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 ( A 1, A 2, A a , and A e ; see Equations (6), (13), and (14)).  \nGNT12 conducted local stratified resistive MHD simulations at 5 AU from the central star with di ff erent sets of the disk mass, ionization strength, and net vertical field strength. The gas surface density was given by Σ = f Σ × 135 g cm -2 with f Σ ( = 1, 2, or 4) being a dimensionless factor. The stellar mass and disk aspect ratio are fixed to M ∗ = M /circledot and H / a = 0 . 05, 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 ( ζ CR), X-rays ( ζ XR), and short-lived radionuclides ( ζ SR), with ζ XR and ζ SR being chosen as 10 and f XR ( = 1 or 20) times the standard values, respectively. The dust-to-gas mass ratio was fixed to 10 -3 , and the size of the dust grains was chosen to be 0 . 1 µ m. 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.  \nOneof the most important parameter is the net vertical magnetic flux 〈 Bz 〉 . 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, 〈 Bz 〉 was chosen in the range 2.68-46 mG. GNT12 also conducted a run with a higher net flux ( 〈 Bz 〉 = 86 mG), but the run resulted in an essentially laminar final state with no sustained MRI turbulence.  \nTable 1 lists the values of f Σ , f XR and 〈 Bz 〉 for all simulations presented by GNT12. Run A1 and B1 are ideal MHD simulations while runs labeled by 'D' include Ohmic di ff usion.", '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.  \nGNT12 calculated the Ohmic resistivity η using the charge reaction network model4 of Ilgner & Nelson (2006). The network consists of free electrons, two species of ions (H + 3 and Mg + ) and charged dust grains. We reproduce the resistivity using the analytic prescription presented by Okuzumi (2009) with the assumption that H + 3 dominates the ions at all heights. 7 To test our calculation, in the upper panel of Figure 2, we plot the vertical profiles of the magnetic Reynolds number H 2 Ω /η 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 ρ = ρ mid exp( -z 2 / 2 H 2 ), where ρ mid = Σ / √ 2 π H is the midplane gas density. Comparing our plot with Figure 1 of GNT12, we confirm that our calculation successfully reproduces their resistivity.  \nWith 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 z the gas motion is unstable if the wavelength λ local( z ) of the most unstable local MRI mode is longer than the gas scale height H . In the presence of Ohmic resistivity, the most unstable wavelength can be approximately given by λ local ≈ max { λ ideal , λ res } , where  \nλ ideal( z ) = 2 π v Az ( z ) Ω (15)  \nλ res( z ) = 2 π η ( z ) v Az ( z ) (16)  \nare the characteristic wavelengths of the unstable modes in the ideal and resistive limits, respectively, with v Az = Bz / √ 4 πρ being the vertical component of the Alfv'en velocity. The most unstable wavelength can be alternatively written as λ local ≈ max { 1 , Λ -1 } λ ideal, where  \nΛ ≡ v 2 Az η Ω = λ ideal λ res (17)  \nis the so-called Elsasser number. The growth rate ν of the local MRI modes is approximately given by ν ≈ min { 1 , Λ -1 } Ω . The instability is strong ( ν ∼ Ω ) when Λ /greatermuch 1, weak ( ν /lessmuch Ω ) when Λ /lessmuch 1, and absent when λ local > H . Thus, the vertical distributions of λ ideal and λ res predict at which height the MRI is unstable.  \n7 In reality, this assumption is not met at low altitudes where the charge transfer from H + 3 to Mg proceeds rapidly. However, this hardly a ff ects the ionization degree (and hence the Ohmic resistivity), since at the low altitudes the recombination mainly occurs on dust grains, for which case the resultant ionization degree is insensitive to the ion composition (see, e.g., Ilgner & Nelson 2006).  \nand  \nTable 1 Model Parameters and Key Observed Quantities of the GNT12 Simulations  \n<!-- image -->  \nFigure 2. Upper panel: magnetic Reynolds number H 2 Ω /η vs. height z for models D1, D1.2, and D1.4 of GNT12 (see also their Figure 1). Lower panel: characteristic MRI wavelength in the laminar state, λ local , 0, for the three models. The solid and dashed segments correspond to λ local , 0 = λ ideal , 0 (or Λ > 1) and λ local , 0 = λ res , 0 (or Λ > 1), respectively. The vertical line marks λ local , 0 = H . The dark and light horizontal bars indicate the ideal and resistive MRI regions defined as H Λ , 0 < | z | < H ideal , 0 and H res , 0 < | z | < H Λ , 0, respectively (see also Figure 1 of OH11).  \n<!-- image -->  \nWhen there is a nonzero net vertical magnetic field 〈 Bz 〉 , it is useful to evaluate λ ideal and λ res assuming that the disk is in the laminar state, i.e., assuming that ρ = ρ mid exp( -z 2 / 2 H 2 ) and Bz = 〈 Bz 〉 at all heights. We denote them by λ ideal , 0 and λ res , 0, respectively. As an example, the lower panel of Figure 2 plots λ local , 0 = max { λ ideal , 0 , λ res , 0 } as a function of z for models D1, D1.2, and D1.4. Note that λ ideal , 0 and λ res , 0 are increasing and decreasing functions of | z | , respectively, because η decreases with | z | while v Az , 0 = 〈 Bz 〉 / √ 4 πρ increases with | z | . For this reason, the region where the MRI is unstable (i.e., λ local < H ) 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 1 . 8 H < | z | < 3 . 2 H , 2 . 1 H < | z | < 3 . 4 H and 2 . 5 H < | z | < 3 . 6 H , respectively.  \nBased on the stability criterion outlined above, OH11 introduced three critical heights H ideal, H Λ , and H res defined by  \nλ ideal( z = H ideal) = H , (18)  \nΛ ( z = H Λ ) = 1 , (19)  \nλ res( z = H res) = H , (20)  \nrespectively. The active layer defined by λ local > H has a vertical extent H res < | z | < H ideal. MRI is stable in the magnetically dominated atmosphere at | z | < H ideal, and in the highη region at | z | < H res. In this paper, we define a dead zone as the region | z | < H res, i.e., H DZ = H res. One can subdivide the active layer into two sublayers H Λ < | z | < H ideal and H res < | z | < H Λ , at which MRI operates strongly and weakly, respectively.  \nAs we did for λ ideal and λ res, we define the critical heights in the laminar state by H ideal , 0, H Λ , 0, and H res , 0. For H ideal , 0, there is an analytic expression (Equation (14) of OH11)  \nH ideal , 0 = [ 2 ln ( β z 0 8 π 2 )] 1 / 2 H , (21)  \nwhere β z 0 ≡ 8 πρ mid / 〈 Bz 〉 2 is the midplane plasma beta measured with the net field strength 〈 Bz 〉 . Table 1 list the values of H ideal , 0, H Λ , 0, and H res , 0 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, H Λ , 0 < | z | < H ideal , 0  \nFigure 3(a) shows 〈 F 2 φ 〉 1 / 2 versus 〈 δρ 2 〉 1 / 2 mid for all the available data. The dashed and dotted lines in Figure 3 show linear scalings (see Equation (3))  \n<!-- image -->  \nDensity Fluctuation  \n/LAngleBracket1  \n∆Ρ  \n/RAngleBracket1  \n1  \n/Slash1  \n2  \nmid  \n/LBracket1  \n/CapOmega  \n/Slash1  \nG  \n/RBracket1  \n<!-- image -->  \nDensity Fluctuation  \n/LAngleBracket1  \n∆Ρ  \n/RAngleBracket1  \n1  \n/Slash1  \n2  \nmid  \n/LBracket1  \n/CapOmega  \n/Slash1  \nG  \n/RBracket1  \nFigure 3. Panel (a): rms azimuthal gravitational force 〈 F 2 φ 〉 1 / 2 vs. rms density fluctuation 〈 δρ 2 〉 1 / 2 mid observed in the GNT12 simulations. The symbols indicate runs B1 (filled circle), D2 (square), D1 and D1. x (open circles), D1WF (triangle), D1-NVFa (plus signs), and D1-NVFb (crosses). The dotted and dashed lines show Equation (22) with A 1 = 1 . 6 and 0 . 16, respectively. Panel (b): Same as panel (a), but here the force amplitude is rescaled by the factor 1 + 4 . 5 H res , 0 / H (see also Figure 4). The solid line indicates Equation (24).  \nand H res , 0 < | z | < H Λ , 0, respectively, for models D1, D1.2, and D1.4.  \nIn principle, the critical heights in a turbulent state can differ from those in the laminar state since λ ideal and λ res depend on ρ and Bz . 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 di ff erence is the largest for H ideal 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 di ff erence is much smaller for H Λ , and is negligible for H res, 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 ( 〈 δρ 2 〉 1 / 2 mid ) and azimuthal gravitational force acting on test particles ( 〈 F 2 φ 〉 1 / 2 ) measured in the GNT12 simulations. These values are taken from Table 3 of GNT12 where they are given in terms of the relative fluctuation 〈 δρ 2 〉 1 / 2 mid /ρ mid and the rms torque 〈 Γ 2 〉 1 / 2 ≡ a 〈 F 2 φ 〉 1 / 2 , respectively.  \n2  \n2  \n2  \n2  \nFigure 4. Ratio 〈 F 2 φ 〉 1 / 2 / 〈 δρ 2 〉 1 / 2 mid vs. dead-zone half width H res , 0. The symbols indicate runs B1 (filled circle), D2 (square), D1 and D1. x (open circles), D1-WF (triangle), D1-NVFa (plus signs), and D1-NVFb (crosses). The dotted, solid, and dashed curves show Equation (23) with A 1 = 1 . 6 and A 2 = 3 . 0, 4 . 5, and 7 . 0, respectively.  \n<!-- image -->  \n〈 F 2 φ 〉 1 / 2 = A 1 GH 〈 δρ 2 〉 1 / 2 mid (22)  \nwith A 1 = 1 . 6 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 shearingout of density waves causes suppression of the random gravity forces.  \nFrom the analysis in Section 2.1, we expect that the e ff ect of the shearing-out can be extracted by taking the ratio between 〈 F 2 φ 〉 1 / 2 and 〈 δρ 2 〉 1 / 2 mid 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 H res defined by Equation (20). Specifically, we here use the value in the laminar state, H res , 0, so that we can calculate it directly from the initial conditions (as we noted in Section 3.1, the value of H res is very insensitive to the presence or absence of turbulence). Figure 4 plots the ratio 〈 F 2 φ 〉 1 / 2 / 〈 δρ 2 〉 1 / 2 mid versus H res , 0. We see a decreasing trend in 〈 F 2 φ 〉 1 / 2 / 〈 δρ 2 〉 1 / 2 mid with increasing H res , 0, 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  \n〈 F 2 φ 〉 1 / 2 〈 δρ 2 〉 1 / 2 mid = A 1 GH 1 + A 2 H res , 0 / H , (23)  \nwhere A 1 and A 2 are the fitting parameters. In Figure 4, the dotted, solid, and dashed curves show Equation (23) with A 1 = 1 . 6 with A 2 = 3 . 0, 4.5, and 7.0, respectively. We find that the set ( A 1 , A 2) = (1 . 6 , 4 . 5) best reproduces the relation between 〈 F 2 φ 〉 1 / 2 and 〈 δρ 2 〉 1 / 2 mid , within an accuracy of factor 2.  \nThus, we have found that the relation between 〈 F 2 φ 〉 1 / 2 and  \n〈 δρ 2 〉 1 / 2 mid can be well represented by  \n〈 F 2 φ 〉 1 / 2 = 1 . 6 GH 1 + 4 . 5 H res , 0 / H 〈 δρ 2 〉 1 / 2 mid , (24)  \nwhich is also shown in the lower panel of Figure 3.', '3.3. Di ff usion Coe ffi cients vs Random Force Amplitude': "The next step is to verify the linear scaling between the diffusion coe ffi cients and 〈 F 2 φ 〉 1 / 2 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 e 0 = 0. They showed that the ensemble averages of ( ∆ a ) 2 and e 2 grow linearly with time ∆ t , indicating a random walk of the particles' motion in the phase space. GNT12 expressed the time evolution of ( ∆ a ) 2 and e 2 in the forms  \n〈 ( ∆ a ) 2 〉 1 / 2 = C σ ( ∆ x ) H ( Ω∆ t 2 π ) 1 / 2 , (25)  \n〈 e 2 〉 1 / 2 = C σ ( e ) H a ( Ω∆ t 2 π ) 1 / 2 , (26)  \nrespectively, where C σ ( ∆ x ) and C σ ( e ) are dimensionless coefficients that depend on the adopted disk model. These values are listed in Table 3 of GNT12.  \nThe di ff usion coe ffi cients Da and De can be read o ff from the values of C σ ( ∆ x ) and C σ ( e ). For Da , we have  \nDa = ( C σ ( ∆ x )) 2 4 π H 2 Ω , (27)  \nwhich directly follows from Equations (7) and (25). For De , we need to convert 〈 e 2 〉 1 / 2 to the eccentricity displacement 〈 ( ∆ e ) 2 〉 1 / 2 for general (nonzero) initial eccentricity. As shown by Yang et al. (2009), this conversion is given by 〈 ( ∆ e ) 2 〉 1 / 2 = √ 2 / (4 -π ) 〈 e 2 〉 1 / 2 . Hence, from Equations (8) and (26), we get  \nDe = 2 ( C σ ( e )) 2 4 π (4 -π ) H 2 Ω a 2 . (28)  \nThe values of Da and De for all the available data are listed in Table 1. The listed values are normalized by H 2 Ω and H 2 Ω / a 2 , respectively, which are the natural units for these di ff usion coe ffi cients in local simulations.  \nNow we calibrate Equations (13) and (14) using the available data. Figure 5 plots Da and De versus the rms azimuthal gravitational force 〈 F 2 φ 〉 1 / 2 measured in the GNT12 simulations. We clearly see the trend Da ∝ De ∝ 〈 F 2 φ 〉 as predicted by Equations (13) and (14). We determine the dimensionless parameters A a and A e 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 A a = 4 . 0 and A e = 1 . 7. The best-fit relations are shown in Figure 5 by the solid lines.  \nTo summarize, we have found that the di ff usion coe ffi cients for the semi-major axis and eccentricity of solid bodies scale with the mean squared amplitude of fluctuating gravity fields as  \nDa = 4 . 0 〈 F 2 φ 〉 Ω 3 , (29)  \nDe = 1 . 7 〈 F 2 φ 〉 a 2 Ω 3 . (30)  \n<!-- image -->  \nFigure 5. Radial di ff usion coe ffi cient Da (upper panel) and eccentricity stirring rate De (lower panel) vs. rms azimuthal gravitational force 〈 F 2 φ 〉 1 / 2 for GNT12 simulations. The symbols indicate runs A1 and B1 (filled circle), D2 (square), D1 and D1. x (open circles), and D1-WF (triangle). The solid lines in the upper and lower panels show Equation (13) with A a = 4 . 0 and Equation (14) with A e = 1 . 7 (or equivalently, Equations (29) and (30)), respectively.  \n<!-- image -->  \nIf we eliminate 〈 F 2 φ 〉 using Equation (23), these scaling relations can be rewritten as a function of 〈 δρ 2 〉 1 / 2 mid and H res , 0. These read  \nDa = 10 (1 + 4 . 5 H res , 0 / H ) 2 ( a 2 H 〈 δρ 2 〉 1 / 2 mid M ∗ ) 2 a 2 Ω , (31)  \nDe = 4 . 4 (1 + 4 . 5 H res , 0 / H ) 2 ( a 2 H 〈 δρ 2 〉 1 / 2 mid M ∗ ) 2 Ω , (32)  \nwhere we have used that Ω = √ GM ∗ / a 3 .", '4. PREDICTING DENSITY FLUCTUATION AMPLITUDES': "Equations (31) and (32) tell us how the di ff usion coe ffi -cients Da and De are related to the amplitude of the gas density fluctuations, 〈 δρ 2 〉 1 / 2 mid . To predict the values of Da and De for given disk parameters, we need to know how 〈 δρ 2 〉 1 / 2 mid depends on these parameters.  \nOH11 provided a simple analytic formula (called the 'saturation predictor') that predicts 〈 δρ 2 〉 1 / 2 mid 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.  \nOH11 Predictor /LBracket1/LAngleBracket1 ∆Ρ 2 /RAngleBracket1 mid 1 /Slash1 2 /RBracket1 OH11 /LBracket1 /CapOmega 2 /Slash1 G /RBracket1  \n<!-- image -->  \nNew Predictor /LBracket1/LAngleBracket1 ∆Ρ 2 /RAngleBracket1 mid 1 /Slash1 2 /RBracket1 new /LBracket1 /CapOmega 2 /Slash1 G /RBracket1  \n<!-- image -->  \nFigure 6. Panel (a): midplane density fluctuation 〈 δρ 2 〉 1 / 2 mid measured in GNT12 simulations vs. the OH11 predictor [ 〈 δρ 2 〉 1 / 2 mid ]OH11 (Equation (33)). The symbols correspond to runs B1 (filled circle), D2 (filled square), D1 and D1. x (open circles), D1-WF (triangle), D1-NVFa (plus signs), and D1-NVFb (crosses). The dotted and dashed lines indicate 〈 δρ 2 〉 1 / 2 mid = [ 〈 δρ 2 〉 1 / 2 mid ]OH11 and 〈 δρ 2 〉 1 / 2 mid = 1 . 2[ 〈 δρ 2 〉 1 / 2 mid ]OH11, respectively. The two open squares mark the runs for which H ideal , ∞ falls below H res , 0. Panel (b): same as panel (a), but here the measured 〈 δρ 2 〉 1 / 2 mid are compared with the updated saturation predictor [ 〈 δρ 2 〉 1 / 2 mid ]new. The solid line shows 〈 δρ 2 〉 1 / 2 mid = [ 〈 δρ 2 〉 1 / 2 mid ]new.", '4.1. The OH11 Predictor versus GNT12 Data': 'The OH11 saturation predictor reads  \n[ 〈 δρ 2 〉 1 / 2 mid ]OH11 = √ 0 . 47 α core ρ mid , (33)  \nwhere α core is a dimensionless coe ffi cient that is proportional to the turbulent accretion stress integrated over height | z | < H ideal (for details, see OH11). OH11 provided an empirical formula that relates α core to key disk parameters such as the net vertical flux 〈 Bz 〉 and the vertical distribution of the resistivity η ( z ). The formula reads  \nα core = 510 β z 0 exp ( -0 . 54 H res , 0 H ) + 0 . 011 exp ( -3 . 6 H Λ , 0 H ) . (34)  \nIn 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.  \nIn Figure 6(a), we compare the measured values of 〈 δρ 2 〉 1 / 2 mid 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,  \nD1. x , D1-WF, and D1-NVFa. However, detailed inspection shows that these observed values are higher than the prediction by ≈ 20% (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).  \nA more important discrepancy can be found for runs D1NVF. In these runs, the net vertical field strength 〈 Bz 〉 was gradually increased from 2.7 mG to 86.0 mG. The results show that 〈 δρ 2 〉 1 / 2 mid initially increases with 〈 Bz 〉 but gets suppressed at 〈 Bz 〉 > 10 . 7 mG. 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 〈 Bz 〉 = 21 . 5 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 〈 δρ 2 〉 1 / 2 mid at high 〈 Bz 〉 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.  \nTo 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 z < H ideal, above which λ local( ≈ λ ideal) exceeds the gas scale height H . The region where local MRI modes exist is also bounded from below by z = H res, below which the Ohmic resistivity stabilizes all MRI modes (i.e., λ local ≈ λ res > H ). Thus, we may measure the vertical width of an MRI-active layer as ∆ z active = H ideal -H res, with H ideal corresponding to the base of the magnetically dominated atmosphere while H res to the the active / dead zone interface.  \nThe 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, H ideal , 0 and H res , 0, as a function of 〈 Bz 〉 . We see that the active layer width in the laminar state ∆ z active , 0 = H ideal , 0 -H res , 0 is already as small as /lessorsimilar H for 〈 Bz 〉 /greaterorsimilar 10 mG. Furthermore, ∆ z active , 0 vanishes at 〈 Bz 〉 ≈ 80 mG, indicating that the magnetically dominated atmosphere completely suppresses the MRI-active layer for this value of 〈 Bz 〉 or larger. This exactly explains what happened in run D1-NVFb, in which the disk returned to a laminar state at 〈 Bz 〉 = 86 mG (see Figure 11 of GNT12).  \nA more quantitative analysis can be made by noting that H ideal further decreases as MRI-driven turbulence develops. Simulations by OH11 show that H ideal measured in the fully turbulent state is smaller than that in the laminar state ( H ideal , 0) by about one scale height (see their Tables 1 and 2). By contrast, H res is nearly independent of the turbulence state  \nNet Vertical Field  \n/LAngleBracket1  \nBz  \n/RAngleBracket1  \n/LBracket1  \nmG  \n/RBracket1  \n<!-- image -->  \nFigure 7. Upper panel: critical heights H ideal , 0 (dashed curve), H res , 0 (dashed curve), and H ideal , ∞ (solid curve) as a function of the net vertical field strength 〈 Bz 〉 for model D1-NVF. Lower panel: midplane density fluctuation δρ mid /ρ mid measured in the D1-NVF runs (plus and cross symbols), compared with the OH11 predictor [ 〈 δρ 2 〉 1 / 2 mid ]OH11 (Equation (33); dashed curve) and the updated saturation predictor [ 〈 δρ 2 〉 1 / 2 mid ]new (Equation (45); solid curve). The two open squares mark the runs for which H ideal , ∞ < H res , 0.  \n<!-- image -->  \nin the active layer because turbulence is always weak on the dead / active zone boundary. The definition of H ideal (Equation (18)) allows us to directly calculate how much H ideal decreases as the turbulence develops at the top of the active layer. If we measure v Az with the rms amplitude of the Bz fields, 〈 B 2 z 〉 1 / 2 , then Equation (18) can be rewritten as  \nρ ideal c 2 s = π 〈 B 2 z 〉 ideal (35)  \nwhere 〈 B 2 z 〉 ideal and ρ ideal are the values of 〈 B 2 z 〉 and ρ at z = H ideal, respectively. We may approximate ρ ideal with the hydrostatic density profile, ρ ideal ≈ ρ mid exp( -H 2 ideal / 2 H 2 ), since the gas pressure dominates over the magnetic pressure at z ≤ H ideal (see also Figure 4(a) of OH11). Thus, the definition of H ideal can be further rewritten as  \nH 2 ideal = 2 H 2 ln ( ρ mid c 2 s π 〈 B 2 z 〉 ideal ) = H 2 ideal , 0 -2 H 2 ln ( 1 + 〈 δ B 2 z 〉 ideal 〈 Bz 〉 2 ) , (36)  \nwhere H ideal , 0 is the value of H ideal in the laminar state (Equation (21)) and 〈 δ B 2 z 〉 ideal ≡ 〈 B 2 z 〉 ideal -〈 Bz 〉 2 is the mean squared amplitude of the fluctuating (turbulent) Bz -fields at z = H ideal. Equation (36) demonstrates that H ideal decreases as turbulence grows.  \nIt is useful to know to what extent H ideal decreases if turbulence is fully developed . The results of the OH11 simulations show that in the fully turbulent state, 〈 δ B 2 z 〉 ideal approximately  \nsatisfies the relation (see Appendix A)  \n〈 δ B 2 z 〉 ideal ∼ 30 〈 Bz 〉 2 . (37)  \nInserting this relation into Equation (36), we find that H ideal should decrease to H ideal , ∞ , where  \nH 2 ideal , ∞ ≈ H 2 ideal , 0 -7 H 2 . (38)  \nIn the upper panel of Figure 7, the solid curve shows H ideal , ∞ versus 〈 Bz 〉 for model D1-NVF. The values of H ideal , ∞ for all the GNT12 simulations are listed in Table 1. Since H res is hardly a ff ected by turbulence, we can estimate the width of the active layer in the fully turbulent state as ∆ z active , ∞ = H ideal , ∞ -H res , 0. As seen in the figure, ∆ z active , ∞ becomes negative for 〈 Bz 〉 /greaterorsimilar 10 mG. This suggests that turbulence cannot be fully developed in run D1-NVFb with 〈 Bz 〉 = 21 . 5 and 43.0 mG;if turbulence were fully developed, then the magnetically dominated atmosphere would completely suppress the active layer.  \nThe 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 e ff ect in Figure 8. This e ff ect is not taken into account in the OH11 saturation predictor (Equation (33)) as H ideal , ∞ was larger than H res , 0 for all the OH11 simulations. This explains why the OH11 predictor overestimates the amplitude of the density fluctuations at high 〈 Bz 〉 .', "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 〈 Bz 〉 .  \nThe model is based on two assumptions. Firstly, we assume that if H ideal , ∞ < H res , 0 then turbulence grows until H ideal reaches H res , 0. This means that the saturated value of 〈 δ B 2 z 〉 ideal is given by Equation (36) with H ideal = H res , 0. Solving the equation with respect to 〈 δ B 2 z 〉 ideal, we get  \n〈 δ B 2 z 〉 ideal =        exp ( H 2 ideal , 0 -H 2 res , 0 2 H 2 ) -1        〈 Bz 〉 2 . (39)  \nCombining this with Equation (37), we obtain the saturation predictor for 〈 δ B 2 z 〉 ideal for general cases,  \n〈 δ B 2 z 〉 ideal = 30 L〈 Bz 〉 2 , (40)  \nwhere L is defined by  \nL = min { 1 , 1 30        exp        H 2 ideal , 0 -H 2 res , 0 2 H 2        -1        } . (41)  \nSecondly, we assume that 〈 δρ 2 〉 mid is proportional to 〈 δ B 2 z 〉 ideal; namely, if 〈 δ B 2 z 〉 ideal is suppressed by factor L , 〈 δρ 2 〉 mid is suppressed by the same factor. This assumption can be expressed as  \nThis expresses that turbulence in the active layer is limited at a low level ( L < 1) when H ideal , ∞ < H res , 0 (see also Figure 8). We will call L the 'saturation limiter.'  \n〈 δρ 2 〉 mid = L〈 δρ 2 〉 mid , ∞ , (42)  \nwhere 〈 δρ 2 〉 mid , ∞ is the value of 〈 δρ 2 〉 mid for fully developed MRI turbulence ( L = 1). We take 〈 δρ 2 〉 1 / 2 mid , ∞ =  \nFigure 8. Schematic description of the physics behind the saturation limiter L (Equation (41)). The gray regions indicate the MRI-active layer defined by H res < z < H ideal. The upward arrows show the magnetic fields, while the downward wavy arrows represent propagation of density waves from the active layer to the midplane. In the laminar state, the upper boundary of the active layer is located at z = H ideal , 0 (Equation (21)). As turbulence develops, H ideal decreases because of the buildup of fluctuating Bz fields, and reaches H ideal , ∞ (Equation (38)) when the turbulence is fully developed. If H ideal , ∞ > H res (case (a)), the active layer can have a finite thickness even with H ideal = H ideal , ∞ , and hence a fully turbulent state is realized ( L = 1). If H ideal , ∞ < H res (case (b)), the turbulence stops developing at the point where H ideal reaches H res, and hence gets saturated at a low level ( L < 1).  \n<!-- image -->  \n1 . 2[ 〈 δρ 2 〉 1 / 2 mid ]OH11 based on the results of the GNT12 simulations for low 〈 Bz 〉 (see Section 4.1). This second assumption can be validated with the following argument. OH11 showed that that the internal energy density of fluctuation, c 2 s 〈 δρ 2 〉 / 2 ρ , is nearly constant along the vertical direction (see also Figure 13 of GNT12). This implies that  \n〈 δρ 2 〉 mid ρ mid ∼ 〈 δρ 2 〉 ideal ρ ideal . (43)  \nIt can also be shown that the internal and magnetic energy densities of fluctuation are nearly equipartitioned at | z | ≈ H ideal, i.e.,  \nc 2 s 〈 δρ 2 〉 ideal 2 ρ ideal ∼ 〈 δ B 2 z 〉 ideal 8 π (44)  \n(see Appendix A for the supporting data). Equations (43) and (44) imply that c 2 s 〈 δρ 2 〉 mid / 2 ρ mid ∼ 〈 δ B 2 z 〉 ideal / 8 π , and hence 〈 δρ 2 〉 mid ∝ 〈 δ B 2 z 〉 ideal.  \nSubstituting 〈 δρ 2 〉 1 / 2 mid , ∞ = 1 . 2[ 〈 δρ 2 〉 1 / 2 mid ]OH11 and Equation (33) into Equation (42), we arrive at the new saturation predictor  \n[ 〈 δρ 2 〉 1 / 2 mid ]new = √ 0 . 68 L α core ρ mid , (45)  \nwhere the factor √ 0 . 68 comes from 1 . 2 × √ 0 . 47. We stress again that the saturation limiter is relevant only when the net flux is so high as to significantly change H ideal. If the net field is not so strong and hence L = 1 holds, then the new predictor di ff ers from the previous OH11 predictor only by factor 1 . 2.  \nTo see how the saturation limiter L operates at high 〈 Bz 〉 , we compare in the bottom panel of Figure 7 the midplane density fluctuation measured in runs D1-NVF with the OH11 predictor [ 〈 δρ 2 〉 1 / 2 mid ]OH11 and updated predictor [ 〈 δρ 2 〉 1 / 2 mid ]new. Comparison with the whose data are shown in Figure 6(b). As we see, [ 〈 δρ 2 〉 1 / 2 mid ]OH11 predicts a monotonic increase in the fluctuation amplitude toward high 〈 Bz 〉 , while [ 〈 δρ 2 〉 1 / 2 mid ]new predicts a flat amplitude at 〈 Bz 〉 /greaterorsimilar 10 mG because of the saturation limiter. The prediction by [ 〈 δρ 2 〉 1 / 2 mid ]new is remarkably consistent with the measurements by GNT12. The new predictor also naturally explains the fact that the disk returns to a laminar state (i.e., the density fluctuations vanish) at 〈 Bz 〉 = 86 mG. We stress that the saturation limiter  \nhas been constructed without further calibration with the numerical data of GNT12 . It is remarkable that it nevertheless explains the observed tendency of the density fluctuation amplitude at high net vertical fields. However, given the limited number of the data points, further support from simulations is desirable to warrant its general validity.  \nThe same refinement should be applicable to the saturation predictor for the gas velocity fluctuation amplitude 〈 δv 2 〉 1 / 2 since δv ∝ δρ for sound waves. The OH11 predictor for the midplane velocity fluctuation amplitude is given by  \n[ 〈 δv 2 〉 1 / 2 mid ]OH11 = √ 0 . 78 α core cs . (46)  \nAs we did for the density fluctuation amplitude, we multiply [ 〈 δv 2 〉 1 / 2 mid ]OH11 by 1 . 2 √ L to obtain a new predictor  \n[ 〈 δv 2 〉 1 / 2 mid ]new = √ 1 . 1 L α core cs . (47)  \nIt should be noted, however, that the accretion stress (or e ff ective viscosity) is not necessarily limited since nonfluctuating, large-scale magnetic fields can contribute to it (Turner & Sano 2008; Gressel et al. 2011; Bai & Stone 2013a,b). Indeed, GNT12 observed no suppression in the accretion stress at high 〈 Bz 〉 where the density fluctuation amplitude already reaches the ceiling (see their Figure 9).", '5.1. The Predictor Functions for the Di ff usion Coe ffi cients': 'In Section 3, we have found how the di ff usion coe ffi cients for the orbital elements of planetesimals, Da and De , are related to the amplitude of density fluctuations at the midplane, 〈 δρ 2 〉 1 / 2 mid . We have also obtained in Section 4 the saturation predictor for 〈 δρ 2 〉 1 / 2 mid as a function of the net vertical flux 〈 Bz 〉 and vertical distribution of the Ohmic resistivity η ( z ). Here, we combine these relations to provide the predictor functions for Da and De . Substitution of Equation (45) into Equations (31) and (32) gives  \nDa = 1 . 1 L α core (1 + 4 . 5 H res , 0 / H ) 2 ( Σ a 2 M ∗ ) 2 a 2 Ω , (48)  \nDe = 0 . 47 L α core (1 + 4 . 5 H res , 0 / H ) 2 ( Σ a 2 M ∗ ) 2 Ω , (49)  \nrespectively, where we have used that ρ mid = Σ / √ 2 π H . The factor L α core come from the fact that the di ff usion coe ffi cients are proportional to 〈 δρ 2 〉 mid. The factor (1 + 4 . 5 H res , 0 / H ) -2 accounts for the suppression of the random gravity due to the shearing-out of the density fluctuations in the presence of a dead zone. These formulae together with the predictor for α core (Equation (34)) allow us to compute the turbulent diffusion coe ffi cients as a function of the net vertical flux and vertical distribution of Ohmic resistivity.', '5.2. Comparison with Previous Recipes Based on Ideal MHD Simulations': 'We check the consistency between our stirring recipe and previous ones proposed by Ida et al. (2008) and Yang et al. (2009, 2012). The previous recipes assume ideal MHD, so we will take H res , 0 = H Λ , 0 = 0 in the following comparison. Then, Equations (48) and (49) reduce to  \nDa ≈ 5 . 5 × 10 -3 ( α 10 -2 )( Σ a 2 M ∗ ) 2 a 2 Ω , (50)  \nDe ≈ 2 . 4 × 10 -3 ( α 10 -2 )( Σ a 2 M ∗ ) 2 Ω , (51)  \nrespectively. Here, we have used that α core ≈ α/ 2 in the absence of a dead zone, where α is the Shakura-Sunyaev viscosity parameter. 8 The above forms are useful when comparing the above equations with previous recipes.', '5.2.1. Ida et al. (2008)': "Ogihara et al. (2007) and Ida et al. (2008) simulated turbulent stirring of planetesimals using a random gravity field model that mimics MRI-driven turbulence originally proposed by Laughlin et al. (2004). Based on the results of these simulations, Ida et al. (2008) proposed a simple formula for the eccentricity growth,  \n〈 e 2 〉 1 / 2 ≈ 0 . 1 γ ( M /circledot M ∗ )( Σ 2400 g cm -2 )( a 1 AU ) 2 ( Ω∆ t 2 π ) 1 / 2 , (52)  \nwhere γ is a dimensionless parameter that characterizes the amplitude of the modeled random gravity fields (see Equation (6) of Ogihara et al. 2007). Here, we augmented the factor M /circledot / M ∗ to the original formula, since Ida et al. (2008) fixed M ∗ to be M /circledot while the magnitude of the modeled random gravity actually scales as M -1 ∗ (see Equation (5) of Ogihara et al. 2007). Ida et al. (2008) expected γ ∼ 10 -2 -10 -3 for ideal MRI-turbulent disks from the results of MHD simulations by Laughlin et al. (2004).  \nTo enable comparison, we rewrite Equation (52) in terms of the eccentricity di ff usion rate defined by Equation (8). If we use Equation (14) together with the relationship 〈 e 2 〉 ∼ 〈 ( ∆ e ) 2 〉 (see Section 3.3), the above formula can be rewritten as  \n[ De ]IGM08 ≈ 0 . 01 ( γ 10 -3 ) 2 ( Σ a 2 M ∗ ) 2 Ω . (53)  \nComparison between [ De ]IGM08 and our De (Equation (51)) allows us to know how the dimensionless parameter γ should  \nbe related to the strength of turbulence, α . We find  \nγ ≈ 5 × 10 -4 ( α 10 -2 ) 1 / 2 . (54)  \nAs we will see in Section 5.2.2, Yang et al. (2012) obtained a consistent result for the case of α ∼ 10 -2 . Baruteau & Lin (2010) also obtained a similar conclusion based on hydrodynamical simulations with Laughlin et al. (2004)'s random gravity model. In the absence of a dead zone, α takes a value of ∼ 10 -2 or larger depending on the strength of the net vertical magnetic flux (Davis et al. 2010; Suzuki et al. 2010). Therefore, our recipe suggests that γ /greaterorsimilar 10 -3 for ideal MRI-turbulent disks, supporting the expectation by Ida et al. (2008).  \nHowever, we stress again that a simple relationship between γ and α like Equation (54) does not apply in the presence of a dead zone. Comparison between Equations (49) and (53) shows that γ must be interpreted as  \nγ ≈ 7 × 10 -4 L 1 / 2 1 + 4 . 5 H res , 0 / H ( α core 10 -2 ) 1 / 2 . (55)  \nThe factor α core express the level of the accretion stress at low altitudes and therefore crudely corresponds to α in Equation (54). A dead zone reduces the value of α core as expressed by Equation (34). However, this factor does not capture all the roles of a dead zone. A dead zone induces the shearing-out of density waves, and thereby further reduces the planetesimal stirring rate as expressed by the prefactor (1 + 4 . 5 H res , 0 / H ) -1 in Equation (55). Suppression of MRI activity at very high 〈 Bz 〉 also reduces γ through the saturation limiter L .  \n5.2.2. Yang et al. (2012)  \nYang et al. (2009, 2012) studied planetesimal stirring in local ideal MHD simulations, with an emphasis on the dependence of the results on the horizontal box size adopted in the simulations. For stratified disks, Yang et al. (2012) proposed analytic expressions for 〈 ( ∆ a ) 2 〉 1 / 2 and 〈 ( ∆ e ) 2 〉 1 / 2 ,  \n〈 ( ∆ a ) 2 〉 1 / 2 = 6 . 6 × 10 -5 ( Lh √ 2 H ) 1 . 35 ξ √ 2 H ( Ω∆ t 2 π ) 1 / 2 , (56)  \n〈 ( ∆ e ) 2 〉 1 / 2 = 7 . 2 × 10 -5 ( Lh √ 2 H ) 1 . 08 ξ √ 2 H a ( Ω∆ t 2 π ) 1 / 2 , (57)  \nwhere ξ ≡ 4 π G ρ mid(2 π/ Ω ) 2 and Lh is the horizontal box size. In terms of the di ff usion coe ffi cients (Equation (7) and (8)), these expressions can be written as  \n[ Da ]YMM12 ≈ 0 . 03 ( Lh 20 H ) 2 . 16 ( Σ a 2 M ∗ ) 2 a 2 Ω , (58)  \n[ De ]YMM12 ≈ 0 . 01 ( Lh 20 H ) 2 . 7 ( Σ a 2 M ∗ ) 2 Ω . (59)  \nComparing the above expressions with Equations (50) and (51), we find that our stirring recipe is reasonably consistent with those of Yang et al. (2012) given that GNT12 adopted the local box of azimuthal extent 12-16 H and that α ∼ 10 -2 in the simulations of Yang et al. (2012, see their Figure 2, bottom panel).  \nWe here note that the horizontal box size dependence appearing in Equations (58) and (59) does not a ff ect the validity  \nof our recipe at a practical level. Yang et al. (2012) suggest that the box size of local simulations needs to be as large as the orbital radius to reproduce the results of global simulations. However, the azimuthal box size of the GNT12 simulations was as large as the orbital radius, since a ∼ 20 H in a typical protoplanetary disk. In fact, Yang et al. (2012) suggested γ ≈ 6 × 10 -4 for α ∼ 10 -2 under the criterion Lh ≈ a (see their Section 6.2), which is precisely consistent with our conclusion (Equation (54)).  \nInvoking the results of Ida et al. (2008) for γ ∼ 10 -3 , Yang et al. (2012) concluded that planetesimals are able to survive mutual collisional destruction even in fully developed MRI turbulence. However, this conclusion should be interpreted with care since this is only true in inner regions of protoplanetary disks where a ∼ 1 AU (see Figure 4 of Ida et al. 2008). It is also important to note that a high net vertical magnetic flux can give a higher α and hence a higher γ . Moreover, growth of planetesimals to protoplanets is not obvious even in the inner disk regions, because the condition for gravitational runaway growth is generally severer than that for surviving collisional destruction. In fact, according to the results of Ida et al. (2008), the growth condition can break down in the entire part of protoplanetary disks for γ ∼ 10 -3 unless the planetesimal size is larger than 100 km. We will discuss this issue in more detail in Paper II.", '5.3. Implications for the Excitation of Planetesimals in Protoplanetary Disks': 'In disks harboring MRI-inactive dead zones, the planetesimal stirring rate De will be reduced considerably with respect to the ideal MHD limit. Two mechanisms contribute to reduce De :  \n- 1. The e ff ective α , or rather, the amplitude of the density fluctuations, is smaller when the MRI is not ideal.\n- 2. The shearing-out will distort the geometry of the density fluctuations, rendering them less e ff ective to stir planetesimals at the midplane.  \nThe first e ff ect is primarily a function of the net vertical magnetic flux 〈 Bz 〉 (see OH11). At low 〈 Bz 〉 , the magnitude of the density fluctuations becomes low; that is, the disks become less turbulent. However, the density fluctuation amplitude cannot be too large even for very high 〈 Bz 〉 as the magnetically dominated atmosphere will then suppress the active layer (this is expressed by the saturation limiter L we introduced in Section 4.3). The second e ff ect is also very important. This will reduce the random torque by up to factor ∼ 10 (see Figure 4) and consequently the stirring rate by up to ∼ 100 when disk possess dead zones.  \nIn this work we have presented a general machinery to capture these e ff ects and to quantify the excitation rate De . Unfortunately, in contrast to Equation (51), in the presence of a dead zone the e ff ective excitation rate is no longer a simple function of the viscosity α . Apart from 〈 Bz 〉 , it also becomes a function of the resistivity profile, η ( z ), which reflects the ionization fraction of the gas. And to calculate the ionization fraction one requires to know the ionization sources, their rates, and the properties of small dust grains.  \nPrevious studies investigating planetesimal accretion (e.g., Morbidelli et al. 2009; Ormel et al. 2010; Weidenschilling 2011; Meschiari 2012) have included a prescription for excitation of planetesimals by density fluctuations. They generally  \nfind (consistent with Ida et al. 2008, Nelson & Gressel 2010, Gressel et al. 2011, GNT12, and our findings) that under ideal MRI conditions (high α and γ ) planetesimals do not accrete, but fragment. To investigate positive outcomes these studies have artificially reduced γ , motivated by a lower α value that may be applicable in dead zones. However, equations like (51) no longer apply in this limit.  \nA physically motivated way to obtain the e ff ective γ would follow the prescriptions outlined in this work, which are fully consistent with the detailed MHD calculations performed by GNT12, but expands it to more general conditions. An example of such an application will be presented in Paper II. In closing, we also want to emphasize that the stirring conditions during the planet formation epoch do not need to be time-independent. For example, planetesimal erosion or fragmentation may inject a large number of small grains, which potentially enlarges the dead zone. Similarly, the net vertical magnetic flux 〈 Bz 〉 may evolve over time due to disk accretion (Rothstein & Lovelace 2008) and / or turbulent magnetic di ff usion (Lubow et al. 1994).', '6. SUMMARY': 'In this study, we have presented a recipe for turbulent stirring of planetesimals in MRI-driven turbulence. From orderof-magnitude estimates, we have derived scaling relations that link the turbulent stirring rates to the amplitude of the density fluctuations and other relevant disk parameters (Section 2). The scalings do not rely on a specific choice of disk parameters, and hence allow to generalize the results of numerical simulations to wider parameter spaces. Our model also accounts for the e ff ects of the shearing-out of density waves on the resulting stirring rates in the presence of a dead zone. We have tested the predicted scalings using the published data of MHDsimulations by GNT12 (Section 3). We have shown that our scaling relations successfully explain the observed data if we fix order-of-unity uncertainties within the relations (Equations (24), (29), and (30)).  \nWe also have updated the saturation predictor for the density fluctuation amplitude in MRI-driven turbulence proposed by OH11 (Section 4). We find that the OH11 predictor overestimates the amplitude when the MRI-active upper layer is significantly suppressed from above by the strongly magnetized, MRI-stable atmosphere. To account for this e ff ect, we have constructed a correction function, which we call the saturation limiter (Equation (41)), on the basis of a layered MRI-turbulent disk model by OH11. The updated predictor function (Equation (45)) successfully reproduces the saturated amplitude of the density fluctuations observed in the GNT12 simulations.  \nCombinations of the scaling relations and saturation predictor (Equations (48) and (49)) enable us to know how the turbulent stirring rate of planetesimals generally depends on disk parameters such as the gas column density, distance from the central star, vertical resistivity distribution, and net vertical magnetic flux. An example of such application will be presented in Paper II, where we examine if runaway growth of planetesimals is possible in turbulent disks. However, given that the simulations our recipe is based on only cover a limited range of the parameter space, it is yet to be warranted if our recipe applies to arbitrary disk conditions. We encourage further testing of its general validity.  \nFinally, we comment that our recipe only takes into account Ohmic di ff usion and neglects other non-ideal MHD e ff ects. Recently, Bai & Stone (2013b) have performed local stratified  \nFigure 9. Mean squared amplitude of the Bz fluctuations measured at z = H ideal, 〈 δ B 2 z 〉 ideal, vs. the squared net flux strength 〈 Bz 〉 2 for all OH11 simulations. The open circles show the data for runs with H res , 0 > 0 while the cross symbols are for runs with H res , 0 = 0. The line indicates 〈 δ B 2 z 〉 = 30 〈 Bz 〉 2 ideal (see Equation (37)).  \n<!-- image -->  \nMHD simulations with both Ohmic resistivity and ambipolar di ff usion and concluded that ambipolar di ff usion can dramatically suppress the turbulent motion of the disk gas. The e ff ect of Hall di ff usion is yet to be examined with stratified simulations, but Wardle & Salmeron (2012) suggest that Hall di ff usion can increase or decrease the vertical extent of the MRIactive layers by an order of magnitude in mass. Inclusion of these non-Ohmic e ff ects will be an interesting extension of this work.  \nWe thank Shigenobu Hirose for sharing with us the results of his numerical simulations. We are also grateful to Hiroshi Kobayashi, Takayuki Muto, Oliver Gressel, Neal Turner, Je ff Cuzzi, Xuening Bai, Hidekazu Tanaka, Takeru Suzuki, Shuichiro Inutsuka, and the anonymous referee for useful comments. S.O. was supported by the Grant-in-Aid for JSPS Fellows (22 · 7006) from MEXT of Japan. C.W.O. acknowledges support for this work by NASA through Hubble Fellowship grant #HST-HF-51294.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.', 'VERIFYING AUXILIARY SCALING LAWS': 'We have derived the new predictor function [ 〈 δρ 2 〉 mid]new assuming two auxiliary scaling relations, Equations (37) and (44). In this appendix, we verify these relations with the simulation data provided by OH11.  \nFigure 9 shows 〈 δ B 2 z 〉 ideal versus 〈 Bz 〉 2 observed in the OH11 simulations. Each data point corresponds to a run with a different resistivity profile. The open circles show the data for runs with H res , 0 > 0 while the cross symbols are for runs with H res , 0 = 0. We see that the data for H res , 0 > 0 are well explained by Equation (37) with an accuracy of factor 2. This is not the case for H res , 0 = 0, for which the high activity of MRI near the midplane influences the activity at high altitudes. However, this fact does not invalidate the use of Equation (37) in deriving the saturation limiter L because the limiter operates only when a large dead zone is present.  \nShown in Figure 10 is c 2 s 〈 δρ 2 〉 ideal / 2 ρ ideal versus 〈 δ B 2 z 〉 ideal / 8 π observed in the OH11 simulations. The  \nFigure 10. Internal energy density of the ρ fluctuations at z = H ideal, c 2 s 〈 δρ 2 〉 ideal / 2 ρ ideal, vs. the energy density of the Bz fluctuations at the same height, 〈 δ B 2 z 〉 ideal / 8 π , for OH11 simulations. The open circles show the data for runs with H res , 0 > 0 while the cross symbols for runs with H res , 0 = 0. The line indicates c 2 s 〈 δρ 2 〉 ideal / 2 ρ ideal = 2 〈 δ B 2 z 〉 ideal / 8 π .  \n<!-- image -->  \ndata are well fitted by c 2 s 〈 δρ 2 〉 ideal / 2 ρ ideal = 2 〈 δ B 2 z 〉 ideal / 8 π , which supports Equation (44). A similar scaling is also observed in local unstratified simulations (Sano et al. 2004).', 'REFERENCES': '- Bai, X.-N. 2011, ApJ, 739, 50\n- Bai, X.-N., & Stone, J. M. 2013a, ApJ, 767, 30\n- Bai, X.-N., & Stone, J. M. 2013b, ApJ, in press (arXiv:1301.0318)\n- Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214\n- Baruteau C., & Lin D. N. C. 2010, ApJ, 709, 759\n- Carballido, A., Stone, J. M., & Pringle, J. E. 2005, MNRAS, 358, 1055\n- Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496\n- Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52 Dzyurkevich, N., Turner, N. J., Henning, Th., & Kley, W. 2013, ApJ, 765, 114\n- Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751\n- Gammie, C. F. 1996, ApJ, 457, 355\n- Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125\n- Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051\n- Gressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291\n- Gressel, O., Nelson, R. P., & Turner, N. J. 2012, MNRAS, 422, 1240  \n(GNT12)  \n- Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, ApJ, 694, 1010\n- Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742\n- Heinemann, T., & Papaloizou, J. C. B. 2009a, MNRAS, 397, 52\n- Heinemann, T., & Papaloizou, J. C. B. 2009b, MNRAS, 397, 64\n- Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292\n- Ilgner, M., & Nelson, R. P. 2006, A&A, 445, 205\n- Jin, L. 1996, ApJ, 457, 798\n- Johansen, A., Brauer, F., Dullemond, C., Klahr, H., & Henning, T. 2008, A&A, 486, 597\n- Johansen, A., Klahr, H., & Mee, A. J. 2006, MNRAS, 370, L71\n- Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022\n- Johnson, E. T., Goodman, J., & Menou, K. 2006, ApJ, 647, 1413\n- Kokubo, E., & Ida, S. 1996, Icarus, 123, 180\n- Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489\n- Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 267, 235\n- Meschiari, S. 2012, ApJ, 752, 71\n- Mohanty, S., Ercolano, B, & Turner, N. J. 2013, ApJ, 764, 65\n- Morbidelli, A., Bottke, W. F., Nesvorny, D., & Levison, H. F. 2009, Icarus, 204, 558\n- Nelson, R. P. 2005, A&A, 443, 1067\n- Nelson, R. P., & Gressel, O. 2010, MNRAS, 409, 639\n- Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849\n- Ogihara, M., Ida, S., & Morbidelli, A. 2007, Icarus, 188, 522\n- Oishi, J. S., Mac Low, M.-M., & Menou, K. 2007, ApJ, 670, 805\n- Okuzumi, S. 2009, ApJ, 698, 1122\n- Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 (OH11)\n- Okuzumi, S., & Hirose, S. 2012, ApJ, 743, L8\n- Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106\n- Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413\n- Ormel, C. W., & Okuzumi, S. 2013, ApJ, in press (Paper II)\n- Ormel, C. W., Dullemond, C. P., & Spaans, M. 2010, Icarus, 210, 507  \n- Pan, L., Padoan, P., Scalo, J., Kritsuk, A. G., & Norman, M. L. 2011, ApJ, 740, 6\n- Perez-Becker, D., & Chiang, E. 2011a, ApJ, 727, 2\n- Perez-Becker, D., & Chiang, E. 2011b, ApJ, 735, 8\n- Rein, H. 2012, MNRAS, L526\n- Rein, H., & Papaloizou, J. C. B. 2009, A&A, 497, 595\n- Rothstein, D. M., & Lovelace, R. V. E. 2008, ApJ, 677, 1221\n- Sano, T., Inutsuka, S., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321\n- Sano, T., & Miyama, S. M. 1999, ApJ, 515, 776\n- Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486\n- Suzuki, T. K., Muto, T., & Inutsuka, S. 2010, ApJ, 718, 1289\n- Turner, N. J., & Sano, T. 2008, ApJ, 679, L131  \nTurner, N. J., Willacy, K., Bryden, G., & Yorke, H. W. 2006, ApJ, 639, 1218 Volk, H. J., Jones, F. C., Morfill, G. E., & Roser S. 1980, A&A, 85, 316 Wardle, M.. & Salmeron, R. 2012, MNRAS, 422, 2737 Weidenschilling, S. J. 2011, Icarus, 214, 671 Wetherill, G. W., & Stewart, G. R. 1989, Icarus, 77, 330 Windmark, F., Birnstiel, T., Guttler, C., et al. 2012, A&A,540, A73 Yang, C.-C., Mac Low, M.-M., & Menou, K. 2009, ApJ, 707, 1233 Yang, C.-C., Mac Low, M.-M., & Menou, K. 2012, ApJ, 748, 79 Youdin, A. N. 2011, ApJ, 731, 99  \nYoudin, A. N., & Goodman, J. 2005, ApJ, 620, 459'}