Planetesimal Dynamics in a Turbulent Disk

# Planetesimal and Protoplanet Dynamics in a Turbulent Protoplanetary Disk: Ideal Stratified Disks

## Abstract

Due to the gravitational influence of density fluctuations driven by magneto-rotational instability in the gas disk, planetesimals and protoplanets undergo diffusive radial migration as well as changes in other orbital properties. The magnitude of the effect on particle orbits can have important consequences for planet formation scenarios. We use the local-shearing-box approximation to simulate an ideal, isothermal, magnetized gas disk with vertical density stratification and simultaneously evolve numerous massless particles moving under the gravitational field of the gas and the host star. We measure the evolution of the particle orbital properties, including mean radius, eccentricity, inclination, and velocity dispersion, and its dependence on the disk properties and the particle initial conditions. Although the results converge with resolution for fixed box dimensions, we find the response of the particles to the gravity of the turbulent gas correlates with the horizontal box size, up to 16 disk scale heights. This correlation indicates that caution should be exercised when interpreting local-shearing-box models involving gravitational physics of magneto-rotational turbulence. Based on heuristic arguments, nevertheless, the criterion , where is the horizontal box size and is the distance to the host star, is proposed to possibly circumvent this conundrum. If this criterion holds, we can still conclude that magneto-rotational turbulence seems likely to be ineffective at driving either diffusive migration or collisional erosion under most circumstances.

Instabilities — Magnetohydrodynamics (MHD) — Planets and satellites: formation — Planet-disk interactions — Protoplanetary disks — Turbulence

## 1 Introduction

The gas disk surrounding a young stellar object is usually argued to be a turbulent accretion disk when the star is still in the process of accumulating materials (see, e.g., Frank et al., 2002). While the molecular viscosity of the gas is negligibly small, shear stress resulting from turbulence can effectively drive angular momentum transport in the disk. At the same time, the dynamics of protoplanetary objects embedded in the disk is significantly affected by their interactions with the turbulent gas, at a stage when these objects are still amassing solids to grow in size (e.g., Papaloizou & Terquem, 2006, and references therein). Understanding particle dynamics in a turbulent gas disk, therefore, will provide insight into the paths along which new planets may eventually form.

One of the most promising mechanisms to drive the turbulence is through the magneto-rotational instability (MRI; see, e.g., Balbus & Hawley, 1998, and references therein). A weakly magnetized, differentially rotating gas disk is unstable to linear perturbations, and after nonlinear effects set in, the gas presumably reaches a statistically steady, sustained, turbulent state. The underlying processes that determine the properties of this saturated state are still under active research (e.g., Jamroz et al., 2008; Latter et al., 2009; Vishniac, 2009; Pessah, 2010; Oishi & Mac Low, 2011). Nevertheless, it is known that the strength of this magneto-rotational turbulence depends on the net magnetic flux threading the disk (e.g., Hawley et al., 1995; Johansen et al., 2006; Yang et al., 2009, hereafter YMM09), and thus can be controlled to the desired level and treated as a parameter.

The gravitational torques exerted on an embedded solid object due to density fluctuations in magneto-rotational turbulence have been shown to be stochastic, rendering the motion of the object a random walk (Laughlin et al. 2004; Nelson & Papaloizou 2004; Nelson 2005; Oishi et al. 2007; YMM09; Nelson & Gressel 2010) on top of the ordered migration expected in a laminar disk (Johnson et al., 2006; Adams & Bloch, 2009; Ketchum et al., 2011). Note that the ordered migration may be either inward or outward in different regions of the disk (Paardekooper et al., 2010, 2011; Lyra et al., 2010). In addition to inducing this radial diffusive migration, the same gravitational force drives the diffusion of orbital eccentricity of the solid objects, increasing their velocity dispersion (Ogihara et al. 2007; Ida et al. 2008; YMM09; Nelson & Gressel 2010). The effect on the orbital inclination, however, has not been investigated yet, as that requires the inclusion of the vertical component of the gravity from the host star. Overall, the gravitational influence of the turbulent gas on particle orbital properties may play an important role in the course of planet formation.

To assess the significance of this process in planet formation scenarios, numerical simulations simultaneously evolving magneto-rotational turbulence and particle dynamics can be conducted and the resulting particle orbital evolution can be statistically measured. However, a consistent picture of the actual magnitude of this effect has been elusive. Using the local-shearing-box simulations, we reported in YMM09 that the response of test particles to the gravity of the turbulent gas is systematically lower than what was previously reported in global disk models. We also found in the same study that the effect significantly depends on the horizontal box size of the numerical models. Recently, Nelson & Gressel (2010) have found similar behavior in their local models and argued that large box size might be necessary to see convergence in the stochastic torque and forcing generated by the turbulence.

Extending the local unstratified disk models of YMM09, we now consider disks with vertical density stratification by including the linearized vertical gravity from the central star. We examine whether a stratified disk model gives results consistent with a comparable unstratified model and whether convergence in particle orbital evolution, if any, can be achieved with local-shearing-box simulations. In Section 2, we describe in detail our numerical models. We report the properties of our modeled magneto-rotational turbulence in Section 3 and the statistical evolution of particle orbital properties in Section 4. We analyze the gravitational field generated by the turbulence in Section 5 and discuss the convergence issues of our local models in Section 6. We conclude our discussion in Section 7.

## 2 Numerical Modeling

Directly extending the work we reported in YMM09, we continue to use the Pencil Code2 to model particles moving in magneto-rotational turbulence. We describe in detail the equations and the relevant numerical techniques we implemented in the code.

### 2.1 Magnetohydrodynamics

We use the local shearing box approximation (e.g., Goldreich & Lynden-Bell, 1965; Brandenburg et al., 1995; Hawley et al., 1995). A local shearing box is a small Cartesian box at a large distance from the host star such that the center of the box revolves at the Keplerian angular speed . The box is always oriented with the -axis directed radially and the -axis azimuthally. In contrast to YMM09, we include the linearized vertical gravity from the host star so that the disk is vertically stratified. We again impose a vertical, external magnetic field to maintain a non-zero magnetic flux. The MHD equations then become

 ∂tρ−32ΩKx∂yρ+∇⋅(ρu) = fD, (1) ∂tu−32ΩKx∂yu+u⋅∇u = −1ρ∇p+(2ΩKuy^x−12ΩKux^y−Ω2Kz^z) (2) +1ρJ×(B+Bext)+fV, ∂tA−32ΩKx∂yA = 32ΩKAy^x+u×(B+Bext)+fR, (3)

where is the gas density, is the gas velocity relative to the background shear flow, is the gas pressure, is the electric current density, , and is the vacuum permeability. The terms , , and are numerical dissipation terms, including both hyper-diffusion and shock diffusion, that are needed to stabilize the scheme. The reader is referred to YMM09 for their description.

To close the system of equations (1)–(3), we assume an isothermal “equation of state”, , where is the isothermal speed of sound. This is an ideal limit of the energy equation, in which heat diffusion is effectively instantaneous. Gas buoyancy should thus be absent and the vertical mixing in the resulting magneto-rotational turbulence might be artificially enhanced (Bai & Goodman, 2009). Nevertheless, Oishi & Mac Low (2009) showed that vertical motions are not completely suppressed even in the case of adiabatic gas, where no heat diffusion occurs, and the differences in the mean kinetic and the magnetic energies of a vertically stratified gas disk between the isothermal and the adiabatic limits may only amount to about a factor of two.

We set up the gas density so that the gas is in vertical hydrostatic equilibrium initially:

 ρ0(z)=ρmexp(−z2H2), (4)

where is the mid-plane gas density and is the vertical disk scale height.3 We fix the vertical dimension of the computational domain at . The initial magnetic vector potential is set to zero, while the external magnetic field is fixed at such a level that the corresponding plasma is in the mid-plane. Gaussian noise of magnitude in each component of the gas velocity is generated to seed the MRI, where is the orbital period at the center of the shearing box.

Our assumption of a net vertical magnetic field allows us to compare our models directly to the unstratified models presented in YMM09. The strength of this field controls the strength of the MRI-driven turbulence, so imposing it is an effective way to adjust the viscous stress generated by the MRI. A similar procedure has been used in other recent publications including Johansen et al. (2007, 2009), Nelson & Gressel (2010), Suzuki & Inutsuka (2009), and Suzuki et al. (2010). Some justification for this assumption can be found in the argument that at high altitude where the gas is tenuous, the magneto-convective motions tend to turn horizontal magnetic fields into vertical ones (e.g., Proctor & Weiss, 1982; Hurlburt & Toomre, 1988; Brandenburg et al., 1995). Therefore, we take the limit that at the -boundary while allowing the vertical component () to change freely. In terms of the vector potential , this boundary condition translates to . In the horizontal directions, we use the standard sheared periodic boundary conditions for the magnetic field, i.e., the field is periodic after a shift compensating for the Keplerian shear flow (Hawley et al. 1995; see also Section 2.2).

We explored several different vertical boundary conditions for the gas and the field. We found that outflow boundary conditions for both led to numerical instability during the initial growth of the MRI. Although Suzuki & Inutsuka (2009) and Suzuki et al. (2010) did publish models using outflow boundary conditions, their disk lost significant amounts of mass through the boundary during their run. We also found that periodic boundary conditions on both gas and field led to random crashes after tens of orbits of saturated turbulence for reasons that remain unclear. On the other hand, the combination of the vertical field boundary condition with periodic boundary conditions in the gas leads to a stable numerical result for the hundreds of orbits required for us to accurately measure the evolution of particle orbital properties (see also Proctor & Weiss, 1982). The horizontal boundary conditions are again sheared periodic. In this setup, the total mass is conserved while saturated turbulent motion can be maintained near the vertical boundaries at a steady level.

While it is interesting in itself to study the coronal structure and possibly the net flow of the gas (Suzuki & Inutsuka, 2009; Suzuki et al., 2010; Flock et al., 2011), we specifically exclude the modeling of the coronal regions to maintain the steadiness of the magneto-rotational turbulence and thus improve the statistical measurement of the particle orbital properties. The gravitational acceleration due to the host star across the boundary does not present any numerical difficulty. Even though the vertical component of the gravitational acceleration changes sign when a fluid element is wrapped around the vertical boundary, the vertical velocity component of the element retains sign such that its velocity and acceleration always have the same sign after reemerging from the opposite boundary, resembling material ballistically falling back down onto the disk. Numerically, no z-derivative of the gravitational acceleration due to the host star is ever evaluated (see Equation (2)), so the discontinuity of the acceleration across the vertical boundary is not an issue here. Turbulence properties near the disk mid-plane, where the gravitational influence of the turbulent gas on the particles is strongest (Oishi et al., 2007), appear insensitive to the choice of vertical boundary conditions (Johansen et al., 2009; Davis et al., 2010; Guan & Gammie, 2011), as we also demonstrate in Section 4.5.

For large box simulations, concerns have been raised that artificial numerical diffusion varying with may occur due to the radial dependence of the shear velocity and thus corresponding truncation errors on a fixed grid (Johnson et al., 2008). Explicit treatment of the shear advection may be needed to eliminate these truncation errors. In the Pencil Code, such a technique, dubbed shear advection by Fourier interpolation (SAFI), has been implemented (Johansen et al., 2009). In our work, we implement SAFI for all of our simulations to remove this undesirable numerical effect.

### 2.2 Poisson Equation with Isolated Boundary Condition in z-direction

To find the gravitational influence of the turbulent gas on the movement of solid particles, we need to solve the Poisson equation for the fluctuation potential :

 ∇2Φ1=4πGρ1, (5)

where is the gravitational constant and is the density fluctuation with respect to the basic state of the gas stratification , which is given by Equation (4). Note that we have neglected the self-gravity of the gas in Equation (2) on the assumption that the disk is gravitationally stable (see YMM09), so the solution to Equation (5) does not affect the gas dynamics.

Although the vertical boundary conditions for the gas density are assumed to be periodic in the MHD equations (see Section 2.1), we find it inappropriate to assume the same in solving the Poisson equation (5). Consider, for example, density perturbations near the top edge of the computational domain. If periodic boundary conditions in the vertical direction were adopted, particles near the mid-plane would experience the gravity of the same density perturbations at almost equal distance from above and below, resulting in cancellation of the gravitational force exerted by these perturbations on the particles. This would amount to undesirable, potentially large, numerical errors for calculating the gravity of density perturbations at high altitude.

Therefore, we instead adopt isolated or open boundary conditions in the vertical direction to solve Equation (5). With these boundary conditions, we assume that there are no fluctuations in density outside the vertical boundary, i.e., for . The boundary conditions in the horizontal directions remain sheared periodic.

To implement approximate isolated vertical boundary conditions with minimal computational cost, we use a variation of the Hockney & Eastwood (1988) fast algorithm to solve the Poisson equation (5), combining it with the Fourier interpolation technique to achieve sheared periodicity (Brucker et al., 2007; Johansen et al., 2007). The essence of this algorithm is to double the computational domain in the -direction by appending regions with zero density fluctuation4 and apply fast Fourier transforms to this expanded domain. Although the vertical boundary conditions remain periodic with Fourier transforms, the nature of gravity means that influence from the periodic copies, now at least away, of the original computational domain of size onto itself is significantly reduced.

We now describe the algorithm to solve Equation (5) in detail. After expanding the vertical dimension of the computational domain and initializing the expanded region with zero density fluctuations, is Fourier transformed in the -direction into , where is the -wavenumber. The result is phase shifted to recover periodicity in -direction: , where is the time step. Then is Fourier transformed in the -direction into , where is the -wavenumber. For the modes , is Fourier transformed in -direction to find , where is the -wavenumber. The fluctuation potential in Fourier space is calculated using the usual convolution theorem:

 ˘Φ1(kx,ky,kz)=−4πG(kx+3ΩKkyδt/2)2+k2y+k2z˘ρ1(kx,ky,kz),for k2x+k2y>0. (6)

It is then inverse Fourier transformed in to find .

In principle, the mode can similarly be treated by Fourier transforms in the vertical direction and convolved by Equation (6) as above. However, this mode has an analytical solution that we can employ with marginal computational cost to improve the solution to Equation (5).5 It represents a horizontally infinite, vertically thin mass layer of constant density at a given altitude , where the density is equal to the horizontal average of the density fluctuation at . (The layer can be considered as either a slab of finite size and constant volume density or as an infinitely thin sheet of surface density ; the solutions outside the layer are the same by symmetry and Gauss’ law.) The gravitational acceleration due to this mode is constant above and below . Therefore, the resultant fluctuation potential can be calculated by

 ~Φ1(kx=0,ky=0,z)=2πG∫~ρ1(kx=0,ky=0,z′)|z−z′|dz′, (7)

where we have arbitrarily defined the reference potential to be zero at each altitude. The discretized version of Equation (7) reads

 ~Φ1(kx=0,ky=0,zj)=2πG∑k~ρ1(kx=0,ky=0,zk)|j−k|Δz2. (8)

Equation (8) remains exact for this mode given that we do not consider reconstruction of the density field within each cell. Finally, we reverse the process to inverse Fourier transform in and to derive the fluctuation potential in real space , incorporating the corresponding phase shifts with opposite sign.

In practice, it is not necessary to allocate storage space for the whole extended domain in the -direction. The appended zero density is only involved in the calculation between the forward and inverse Fourier transforms in and . In addition, the convolutions in (Equations (6) and (8)) for different modes are independent of each other. Therefore, these convolutions can be distributed among processors and performed sequentially by allocating only a single one-dimensional working array along the -direction.

### 2.3 Particle Dynamics

We continue to use the zero-mass approximation to model our solid particles as in YMM09. In this approximation, particles behaves as test particles and only respond to the gravity of the host star and the gas. We ignore the viscous drag forces between particles and gas. This remains a good approximation for kilometer-sized planetesimals (e.g., Oishi et al., 2007). We also ignore mutual gravitational interactions between particles, which helps us isolate the net effect induced by hydromagnetic turbulence.

The equations of motion for each particle, therefore, become

 dxpdt = up−32ΩKxp^y, (9a) dupdt = (2ΩKup,y^x−12ΩKup,x^y−Ω2Kzp^z)+g0−∇Φ1. (9b)

The vector is the position of the particle in the shearing box, while is the velocity of the particle relative to the background shear flow. In Equation (9b), the terms inside the parentheses stem from the linearized gravity of the host star and the Coriolis and centrifugal forces in the co-rotating frame of the shearing box. The remaining terms on the right-hand side represent the acceleration attributed to the gas: is the gravitational acceleration due to the basic state of the gas stratification and has the analytical expression

 g0(z)=−2π3/2GρmHerf(zH)^z, (10)

while is the gravitational potential due to density fluctuations with respect to the basic state and is the solution of the Poisson Equation (5) (Section 2.2). Note that by separating the density perturbation from the equilibrium density stratification and treating the gravity of the equilibrium state exactly, we improve the accuracy of the gravitational potential resulting from the turbulence.

The position and velocity of each particle is evolved in time by integrating the equations of motion (2.3) simultaneously with the third-order Runge-Kutta steps advancing the MHD equations. In addition to the Courant conditions set by the MHD equations, the time-step is limited by the absolute maximum of velocity so that no particles can cross more than half a grid zone each timestep. We calculate the gradient of the fluctuation potential on the grid and then quadratically interpolate it onto each particle.

We uniformly distribute particles in a horizontal plane. We do not allow these particles to move until after 20 orbital periods, when the turbulence has reached a statistically steady state. The particles can have an initial eccentricity by setting an initial velocity so that they are at the apogee of their orbits (see YMM09). They can also have an initial inclination by placing the initial particle plane at a distance to the mid-plane. Particles are wrapped around when crossing any of the six boundary planes of the shearing box.

## 3 Properties of Magneto-Rotational Turbulence

In this section, we present several properties of the magneto-rotational turbulence in our numerical models and discuss their convergence with the box size and the resolution.

Figure 1 shows the horizontally averaged density fluctuation, inverse plasma , and parameter in the mid-plane () as a function of time, for resolutions up to 64 points per and horizontal box dimensions up to . The magnitude of the density fluctuation is indicated by the rms value relative to the equilibrium density , the plasma is the ratio of the equilibrium thermal pressure to the magnetic pressure, and is calculated by (e.g., Brandenburg, 1998)

 α=√23⟨ρuxuy−BxBy/μ0⟩ρ0c2s, (11)

which includes both Reynolds and magnetic stresses.6

The magneto-rotational turbulence in the mid-plane of our models reaches a statistically steady state at about . In this saturated state, we see in general that the larger the horizontal size of the box, the smaller the amplitude of oscillation in the turbulence properties. Little difference exists between an 884 box and a 16164 box at the same resolution, indicating convergence with horizontal box dimension. On the other hand, the amplitude of oscillation at saturation level increases with resolution for fixed box dimensions. There exists some evidence that the magnetic activity and thus the value increase with resolution, especially for the case of 224 box, but this effect seems to become smaller with larger box sizes. We also note that there exists a curious trend of increase in density perturbation with time for the case of an 884 box at the resolution of 64 points/. In spite of these trends, the time average of each turbulence properties remains roughly the same for different resolutions and thus shows convergence. We note that the turbulence properties at the saturated state in the mid-plane of our vertically stratified disks agrees with those in our unstratified disks of YMM09 (see also Guan & Gammie, 2011). The means and standard deviations of these properties for each set of resolution and box dimensions are reported in Table 1, which includes the magnetic stress and the Reynolds stress as well.

Figure 2 shows the time-averaged vertical profiles of density perturbation, inverse plasma , and parameter at saturation level over a period of . Each of the three properties increases with , indicating stronger turbulence at higher altitude. The increasing activity with altitude is related to the increasing importance of the external magnetic field, which can be quantified by (Section 2.1). This agrees with the trend found in models of unstratified disks (Hawley et al. 1995; Johansen et al. 2006; YMM09). Notice that at our highest altitude, and thus the corona is not modeled in our simulations. Given that our vertical height only amounts to 2, the gas flow in our computational domain remains marginally stable against magnetic buoyancy during all stages of the development of the MRI (e.g., Shi et al., 2010; Guan & Gammie, 2011; Simon et al., 2011). As shown in Figure 2, the time-averaged inverse plasma beta near the vertical boundary amounts to only about 0.4, and thus on average, thermal pressure still dominates magnetic pressure there.

## 4 Orbital Properties of Massless Particles

We now report the response of zero-mass particles to the gravity of the density fluctuations of the statistically steady, numerically convergent magneto-rotational turbulence described in Section 3. The reference time in the following discussion is the time at which the turbulent gas reaches its saturated state and the particles are allowed to move.7

The evolution of the mean orbital radius of one particle can be found by averaging the radial position of the particle over each epicycle motion. For the case of ideal unstratified disks, the distribution of particles in terms of the orbital radius change can be described by a time-dependent normal distribution centered at zero:

 f(Δx,t)=1√2πσx(t)exp[−Δx22σ2x(t)], (12)

where is the orbital radius change from the initial radius at and is the time-dependent standard deviation. We reported in YMM09 that depends on the properties of the gas disk and can be concisely expressed by

 σx(t)=CxξH(tP)1/2, (13)

where is a dimensionless proportionality constant8 and is a dimensionless quantity indicating the strength of the gas gravity. Equations (12) and (13) demonstrate the diffusive nature of particle radial migration driven by magneto-rotational turbulence.

For the case of ideal stratified disks studied here, we again find the evolution of particle orbital radius can be described by Equations (12) and (13). The dimensionless quantity is now defined by

 ξ≡4πGρmP2=4(2π)3/2/Qg (14)

in terms of the mid-plane gas density , where is the Toomre stability parameter for the gas. In Table 2, we report the measured values of the constant from our simulations at different resolutions and horizontal box sizes for the case of particles with zero initial inclination. Comparing the value from the 224 stratified model at a resolution of 64 points/ with that from the unstratified model at the same resolution and horizontal box size (see Equation (15) of YMM09), we find excellent agreement in the two values.

As can be seen from Table 2 and Figure 3, the value of remains roughly the same with different resolutions for fixed box dimensions (except the anomaly shown by the 224 box at a resolution of 64 points/). Conversely, it is significantly dependent on the horizontal box size . The larger the size, the stronger the influence of the turbulence on the particle orbital radius. This relationship can be represented by the following power-law fit:

 Cx≃6.6×10−5(Lh/H)1.35. (15)

We find no evidence of convergence with horizontal box size up to , the largest size we have investigated.

Equation (13) can be transformed into the diffusion coefficient introduced by Johnson et al. (2006) for describing the radial random walks of orbiting particles induced by turbulent torques, in terms of the Keplerian orbital angular momentum . The reader is referred to YMM09 for a detailed description of this procedure. We emphasize here that this transformation involves no assumption about the correlation time of the stochastic torques and is thus a direct measurement of . Using a heuristic choice of dimension for the diffusion coefficient, Johnson et al. (2006) defined a dimensionless parameter to represent the magnitude of . We report our measured values of in Table 2. The dependence of on the horizontal box size in our models for is shown in Figure 4 and can be written as

 ϵ≃6.5×10−6(Lh/H)2.69. (16)

### 4.2 Eccentricity

The amplitude of each epicyclic oscillation of a particle gives the instantaneous eccentricity of the particle orbit. We reported in YMM09 that in ideal unstratified disks, the distribution of particles in terms of the eccentricity deviation should be a time-dependent normal distribution centered at zero, as long as the particles have non-negligible initial eccentricity (c.f. Equation (12)):

 f(Δe,t)=1√2πσe(t)exp[−Δe22σ2e(t)], (17)

where is the eccentricity deviation from the initial eccentricity at and the time-dependent standard deviation can be written as

 σe(t)=Ceξ(HR)(tP)1/2, (18)

where is a dimensionless proportionality constant. We emphasize that the time-dependent Rayleigh distribution9 found for particles with zero (or negligible) initial eccentricity

 f(e,t)=eσ2e(t)exp[−e22σ2e(t)] (19)

is a manifestation of Equation (17) since the eccentricity is a positive definite quantity. Equations (17) and (19) share the same time-dependent parameter , and there exists no evidence that depends on the initial eccentricity .

We find that the same evolution of particle distribution in eccentricity deviation also holds for the case of ideal stratified disks. The measured values of the constant for different resolutions and horizontal box sizes when particles have zero initial inclination are listed in Table 2. The same comparison performed in Section 4.1 for the orbital radius evolution indicates that the eccentricity evolution in a stratified disk again agrees very well with that in an unstratified disk (see Equation (17) of YMM09). These comparisons demonstrate that a local unstratified disk model gives results consistent with the mid-plane of a local stratified disk model that has the same physical conditions except vertical stratification.

As in the case of orbital radius discussed in Section 4.1, the eccentricity evolution of the particles does not noticeably depend on the resolution for given box dimensions, while it is strongly affected by the horizontal box size (Table 2 and Figure 3). A power-law regression gives

 Ce≃7.2×10−5(Lh/H)1.08, (20)

which is close to a linear relation to the horizontal box size. We discuss the box-size effect on both the orbital radius and the eccentricity in Section 6.

The eccentricity driven by magneto-rotational turbulence enhances orbital crossing among planetesimals and thus increases the chance of collisions between them. Ida et al. (2008) defined a dimensionless parameter to represent the strength of this effect. Equation (18) can be used to estimate the value of , and the reader is referred to YMM09 for a detailed description of this procedure. We report in Table 2 our measured values of . As shown in Figure 4, the dependence of on the horizontal box size in our models for can be described by

 γ≃3.6×10−5(Lh/H)1.08. (21)

### 4.3 Inclination

The only orbital property of a particle that cannot be measured in an unstratified disk model is the inclination . In a stratified disk, vertical linear gravity from the host star (the third term in the parentheses on the right-hand side of Equation (9b)) provides a restoring force such that particles oscillate in the -direction about the mid-plane, as demonstrated in Figure 5. Note that because the particles also experience the gravity of the gas, the period of the oscillation is in the linear limit (derived from Equations (9b) and (10)), which is slightly shorter than the orbital period . We can calculate the induced inclination (in radians) for a single particle by , in which and are the maximum and the minimum vertical positions in one oscillation, respectively, provided that .

Figure 6(a) shows the histograms of particles with zero initial inclination sorted in bins of instantaneous inclination at three different times. Similar to the eccentricity distribution for particles with described by Equation (19), the inclination distribution resembles a time-dependent Rayleigh distribution

 f(i,t)=iσ2i(t)exp[−i22σ2i(t)], (22)

where is a time-dependent parameter denoting the width of the distribution.

One might expect the inclination distribution of particles with nonzero initial inclination would evolve as a time-dependent normal distribution with fixed center and increasing width, similar to the eccentricity distribution of particles with nonzero initial eccentricity. However, Figure 6(b) shows different behavior: the distribution increasingly deviates from a normal distribution with time. The peak leans toward the mid-plane and asymmetry develops with more particles on the left side of the peak (i.e., smaller inclination). This indicates the diffusion is stronger near the mid-plane, which in fact is consistent with vertical stratification of the gas density. It would be interesting to measure the dependence of the diffusion coefficient on vertical height, but given our limited sample of initial inclinations, this remains to be investigated. We therefore focus our attention on the case of particles with zero initial inclination.

Figure 7 shows for disks with varying gravity parameter and particles with varying initial eccentricity but zero initial inclination . Note that when is normalized by , all the curves roughly coincide. This indicates that is linearly dependent on both , but independent of , similar to what we reported for the evolution of the mean radius and the eccentricity. The results can thus be summarized by the following expression:

 σi(t)=Ciξ(HR)(tP)1/2, (23)

where is a dimensionless proportionality constant. Table 2 lists our measured values of the constant for different resolutions and box dimensions when particles have zero initial inclination.

We note that in contrast to the orbital radius and the eccentricity, the evolution of the orbital inclination is not significantly affected by the horizontal box size. The inclination of the particle orbits is mainly affected by the vertical structures of the density field, which in turn is related to the vertical scale height of the density stratification. The horizontal structures in the magneto-rotational turbulence may have little correlation with the vertical structures so that the horizontal box size has little effect on the vertical motions of the particles.

### 4.4 Velocity Dispersion

Finally, all three components of the velocity dispersion among the particles as a function of time can be measured in our stratified disk models. Similarly to what we reported in YMM09, each component assumes the same form:

 σu,i(t)=Siξcs(tP)1/2, (24)

where the index is either , , or and is the corresponding dimensionless proportionality constant. We emphasize that always holds because of the fixed ratio of amplitudes in and directions in the epicycle motions of the particles, which has been verified in our simulations. The values of and at different resolutions and box dimensions for the case of particles with zero initial inclination are listed in Table 2. Note that the value of for the 224 box at a resolution of 64 points/ is consistent with that from the unstratified disk with the same horizontal size and resolution we reported previously (see Equation (18) of YMM09), a further confirmation of the consistency between stratified and unstratified models discussed in Sections 4.1 and 4.2.

As shown by Figure 8, we do not see evident dependence of velocity dispersion on resolution of our models. On the other hand, the horizontal component (and thus ) significantly depends on the horizontal box size while the vertical component does not. This is in accordance with the dependence of the three orbital properties found in previous subsections. From our measured values for , we quantify and with the following expressions:

 Sx ≃ (1.0×10−4)(Lh/H)1.08, (25) Sz ≃ 2.8×10−4. (26)

### 4.5 Boundary Condition Effects

In order to confirm that our numerically convenient but perhaps physically inconsistent assumption of vertical field and periodic gas boundary conditions in the -direction does not change our results significantly, we compare turbulence properties and particle dynamics from models with different vertical boundary conditions for the magnetic fields run at a lower resolution of 16 points per disk scale height in large boxes with size . The gas density and velocity remain periodic. Figure 9 shows that within about one disk scale height, the vertical structure of the rms density perturbation, the plasma , and the parameter are in quantitatively close agreement between the models. Note that the model with periodic boundary conditions creates a coronal region as low as , where . More importantly, Figure 10 shows from top to bottom the time evolution in the standard deviations of particle orbital radius, eccentricity and inclination for the two models. The orbital radius and eccentricity agree very well irrespective of the vertical boundary conditions, and the effect on inclination differs by only 30%, likely due to increased activity at high altitudes for the model with periodic magnetic fields. Thus, we believe that our results on particle dynamics are robust.

## 5 Properties of the Stochastic Force Exerted by the Gas

We have demonstrated in Section 3 that various properties of the saturated magneto-rotational turbulence in our simulations converge with both resolution and box dimensions. However, while converging with resolution, the response of massless particles to the gravity of the turbulent gas does not similarly converge with the horizontal box size, as presented in Section 4. This result raises serious questions about the validity of using the local-shearing-box approximation to simulate the dynamics of any particles under gravitational influence of magneto-rotational turbulence. In order to understand the lack of convergence with box dimension in the orbital evolution of massless particles, we further study the properties of the gravitational force exerted by the turbulent gas.

### 5.1 Magnitude and Fourier Amplitudes of the Force

Table 3 lists the radial and the azimuthal components of the gravitational force exerted by the turbulent gas, and , respectively. Since the force is stochastic, we ensemble average the root mean square of each force component over all particles, and then report it as time average and variation scaled by , where is the mass of a particle. Note that in our models, particles only follow the gravitational potential generated by the gas and the host star, so the acceleration of each particle is independent of its mass. This is an advantage in that the measurements reported in this paper can be used to gauge the gravitational influence of magneto-rotational turbulence on solid particles or planetary objects of any mass and thus its significance compared with other relevant interactions.

As can be seen in Table 3, both the radial and the azimuthal components of the force generally increase with the horizontal box size. It is not clear if there exists any convergence in the radial component towards larger box size, since the 16164 box at a resolution of 32 points per is the only model that contradicts the general increasing trend. Certainly, the azimuthal component (i.e., the torque) shows no sign of convergence up to the largest box size we have investigated, with a regression power-law index of about 0.6. Nelson & Gressel (2010) also reported in their recent local models a slight trend of increasing torque with increasing horizontal box size (see their Figure 8 and the corresponding discussion).

Since the density structure of the turbulent gas determines its gravitational force on the particles, it is informative to consider how the structure changes with resolution and horizontal box size. The density structure of the gas is best analyzed in Fourier space, using the Fourier amplitudes defined in Section 2.2.10 Figure 11 plots the time-averaged Fourier amplitudes of the gas density in the mid-plane along either the radial or the azimuthal direction. The top-left panel shows the amplitude as a function of for . In general, the largest amplitude resides at the longest wavelength while monotonically decreasing with increasing wavenumber (Johansen et al., 2009). For any given box dimensions, increasing resolution has little effect on each amplitude, apart from extending the profile toward higher wavenumber. Flattening of the profile for wavelengths 8 is hinted for 16164 boxes. We find no evidence of a downward turn toward even longer wavelength. On the other hand, increasing the horizontal box size lowers the overall profile of the Fourier amplitudes. Interestingly, the amplitude of the longest wavelength mode remains roughly constant.

The bottom-left panel of Figure 11 plots the summation of Fourier amplitudes over all at any given . In contrast to the amplitudes for , the horizontal box size has little effect on the amplitudes for most of the wavelengths. Nevertheless, increasing resolution indeed increases the inertial range and resolves more power toward shorter wavelength.

The right-hand column of Figure 11 plots the azimuthal counterpart of the left-hand column. We find similar features in the azimuthal amplitudes as those found in the radial ones. The only difference is that the flat part of the spectrum in the azimuthal direction is short compared with that in the radial direction, as can be seen in the bottom-right panel. Power that is not captured in the long-wavelength range of the spectrum in the 224 boxes is probably the cause for artificially higher amplitude at (the top-right panel), which may in turn be responsible for the anomaly found in the particle dynamics shown in Section 4.

The Fourier amplitudes shown in Figure 11 help us understand the relationship between the gravitational force exerted by the gas and the horizontal box size. The radial and the azimuthal force components for any given (horizontal) Fourier mode are related to the corresponding Fourier amplitude by (see Equation (6))

 Fx(k)∼4πG(kxk2)|˘ρ(k)|≲4πG|˘ρ(k0,0)|k0∝Lh, (27) Fy(k)∼4πG(kyk2)|˘ρ(k)|≲4πG|˘ρ(0,k0)|k0∝Lh, (28)

where we have used the facts that the dominant modes in the radial and the azimuthal directions are the longest wavelength ones and , and both are about constant irrespective of resolution and box size. Therefore, both components should be at most linearly proportional to the horizontal box size in our simulations. Although involving some simplifications, Equations (27) and (28) at least qualitatively explain the general trend of increasing stochastic force with horizontal box size.

### 5.2 Correlation Time

We next evaluate the correlation time of the stochastic torque. This quantity is in fact not well defined in the literature, and different authors have different preferences for its evaluation. For the purpose of consistency with the framework of radial diffusive migration of particles driven by random torque, we adopt the definition of Johnson et al. (2006):

 τc≡D(J)/¯¯¯¯¯¯¯¯δΓ2, (29)

where is the radial diffusion coefficient as a function of the Keplerian orbital angular momentum , is the fluctuating part of the torque, and the overline denotes an ensemble average. We repeat the relationship between and the dimensionless parameter here (see Johnson et al., 2006, and YMM09):

 D(J)=2.1×10−316πϵξ2(HR)2(J2P). (30)

To calculate the correlation time with our definition, both the diffusion coefficient and the magnitude of the stochastic torque are needed. As emphasized in YMM09, we have all the information about the particle distribution against each orbital property as a function of time, and thus we can directly evaluate the diffusion coefficient in the context of a diffusion equation, which is reported in Section 4.1 and Table 2. Furthermore, we also have the ensemble-averaged rms value of the stochastic torque, which is reported in Section 5.1 and Table 3. Equations (29) and (30) can then be used to calculate the correlation time , and the result is listed in Table 3.

As indicated in Table 3, the correlation time of the stochastic torque does not appreciably depend on the resolution but it is significantly affected by the horizontal box size. The larger the horizontal box size, the longer the correlation time of the torque; this trend was also reported by Nelson & Gressel (2010). Up to the largest box size we have investigated (), it is not clear if will converge with . This behavior is similar to that of the torque magnitude, as discussed in Section 5.1.

As a side note, one may want to evaluate the correlation time of the stochastic torque without referring to stochastic particle orbital evolution and invert Equation (29) to estimate the radial diffusion coefficient indirectly. We have shown in YMM09, for the case of ideal unstratified disks, that the following formula, motivated by the definition of under the framework of the Fokker-Planck formalism, serves this purpose:

 τc≈∫∞0¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ACF(τ)dτ2¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ACF(0), (31)

where is the autocorrelation function of the stochastic torque as a function of the time lag . Here we repeat the exercise and use Equation (31) to estimate for the case of ideal stratified disks. The result is listed in Table 3 and shows remarkable consistency with the values obtained directly from the evolution of particle orbital distribution. Therefore, Equation (31) remains a useful estimator of the correlation time .

To further understand the trend of increasing correlation time with increasing horizontal box size, we plot the autocorrelation function of the stochastic torque in Figure 12. As also demonstrated by Nelson & Gressel (2010), the larger the box, the less oscillatory the autocorrelation function is and the less pronounced is the negative contribution to the radial diffusion of the particles. This behavior in turn is consistent with increasing values of the correlation time and the radial diffusion coefficient.

In addition to the stochastic torque, a similar analysis may be conducted for the radial component of the gravitational force exerted by the gas. We find that the correlation time of the radial force also increases with the horizontal box size. However, our models indicate a much longer correlation timescale, which may amounts to several tens of orbital periods. Given that our simulations only lasted for about 100, the magnitude of the correlation time is not well calculated and much longer simulations may be required to cover the full spectrum of the autocorrelation function for the radial force.

Johansen et al. (2009) have discussed the dominant role of the longest wavelength, purely radial Fourier mode in great detail and suggested that these kind of persistent structures are similar to the zonal flows found in many other astrophysical systems. As noted by the authors, these axisymmetric structures can not generate torques that affect the orbital radius of the embedded particles. While the radial forces due to the zonal flows can still affect the eccentricity of the particles, their long correlation time indicates that they may not be responsible for the diffusive evolution of the particle eccentricity, as shown in YMM09 and Section 4.2.11 On the other hand, the longest wavelength, purely azimuthal Fourier mode in the gas structure12 is significant enough that both the radial diffusive migration and the eccentricity evolution of the particles are affected by the large-scale azimuthal density perturbations and thus the box dimensions.

## 6 Relation of Box Size to Particle Orbital Evolution

As described in Section 5, the magnitude and the correlation time of the gravitational force exerted by the turbulent gas correlates with the horizontal size of our local shearing box. This behavior essentially explains the correlation of the stochastic evolution of the particle orbital radius and eccentricity with the horizontal box size reported in Section 4. These results challenge the validity of the local shearing box to describe large-scale density structures in magneto-rotational turbulence.

By measuring the two-point correlation functions of density, velocity, and magnetic field in magneto-rotational turbulence, Guan et al. (2009) found extended density features in contrast to well localized magnetic structures. The authors suggested that the density features are propagating acoustic waves excited by the turbulence (see also Heinemann & Papaloizou, 2009a, b). It is not clear yet what the dissipation length scale for these waves is, which is crucial to understanding whether a local shearing box can capture the largest-scale structures in the density fluctuations excited by magneto-rotational turbulence.

On the other hand, a global disk model may require high resolution in order to self-consistently produce large-scale structures. For example, Johansen et al. (2009) argued that the inverse cascade of magnetic energy from small scales to large scales might be responsible for ultimately launching zonal flows. To confirm this mechanism in a global context, a model that is capable of resolving at least the correlation lengths in the magnetic structures might be necessary.

### 6.1 Bridge to a Global Disk Model?

Given the correlation with horizontal box size reported in this work, it is not clear how a local shearing box can be connected to a global disk model so that both models give consistent results on the particle dynamics under the gravitational influence of density fluctuations in magneto-rotational turbulence.

Nevertheless, we conjecture at this point that the criterion might provide a physical length scale for a local shearing box. At this scale, the local-shearing-box approximation formally breaks down since it assumes . The curvature terms become important and this might trigger turbulent eddies to damp propagating waves. If this conjecture could be verified, Equations (15), (20) and (25) would prove useful in evaluating the gravitational influence of the turbulent gas on particle dynamics. Substituting for in these equations, our models predict that

 ⎧⎪ ⎪⎨⎪ ⎪⎩Cx≃6.6×10−5(H/R)−1.35Ce≃7.2×10−5(H/R)−1.08Sx≃1.0×10−4(H/R)−1.08% for α∼10−2. (32)

Nelson & Gressel (2010) presented the most recent work on stochastic planetesimal orbital evolution in a global disk model. Their measurement of the orbital radius and radial velocity dispersion can be directly compared with ours. These authors reported in their global model G3, whose was the closest to ours, that and . In this specific model, a disk aspect ratio13 of and a strength of gas gravity of at 5 AU were adopted. In comparison, Equations (13), (24), and (32) give and . Our results are only about a factor of two lower than those of Nelson & Gressel (2010). Therefore, the criterion for a local shearing box seems to offer some consistency between global and local models.

These comparisons are not meant to be accurate and conclusive, though, since the proposed criterion invites significant uncertainty. Further investigation into the large-scale structures of magneto-rotational turbulence and the discrepancy between global and local disk models is still needed.

### 6.2 Implications for Planet Formation and Migration

If the criterion holds, Equations (16) and (21) imply

 {ϵ≃6.5×10−6(H/R)−2.69γ≃3.6×10−5(H/R)−1.08% for α∼10−2. (33)

Consider a region at  AU in a minimum-mass solar nebula (MMSN; Hayashi, 1981), where . Equations (33) give and . This value of the parameter is about one order of magnitude higher than was reported in YMM09 for the case of a 222H unstratified box, while the value of the parameter is only about a factor of three larger. Note that according to Equations (33), both and decrease with increasing disk aspect ratio , which often increases with increasing radial distance .

We suggested in YMM09 that in a typical protoplanetary disk, the radial diffusive migration of protoplanets induced by magneto-rotational turbulence may be unimportant compared to secular migration. According to Johnson et al. (2006), for the diffusive migration to be able to dominate over type I migration, the parameter should be greater than or about 0.1–1 for the MMSN. Our new measurement remains orders of magnitude smaller than this transition value. Therefore, our previous conclusions on the insignificance of planetary diffusive migration may still hold even though is increased by an order of magnitude. However, note that recent calculations for type I migration in non-isothermal disks (e.g., Lyra et al., 2010; Paardekooper et al., 2010, 2011) may negate the dominance of secular migration over stochastic migration driven by magneto-rotational turbulence.

We also argued in YMM09 that kilometer-sized planetesimals moving in magneto-rotational turbulence survive mutual collisional destruction, except in the inner region of a young protoplanetary disk. This was based on the results of Ida et al. (2008) for the cases of and . Since our new measurement does not fall outside this range, the same conclusion still applies.

In contrast, Nelson & Gressel (2010) argued that the velocity dispersion of kilometer-sized planetesimals excited by magneto-rotational turbulence might be so large that these planetesimals should be erosive to each other. Before a consistent scenario for the survivability of planetesimals could be assembled, however, several factors remain to be accounted for. First, eccentricity damping due to tidal interaction between a planetesimal and its surrounding gas acts to circularize the orbits of the planetesimals (e.g., Goldreich & Tremaine, 1980). This mechanism was absent in the models of Nelson & Gressel (2010) and thus their measurement of the velocity dispersion of planetesimals should be considered as an upper limit.14 On the other hand, using analytical arguments to estimate the equilibrated eccentricity of planetesimals, Ida et al. (2008) took three possible eccentricity damping effects into account — tidal interaction, gas drag, and inelastic collisions — and thus their estimate of the velocity dispersion of planetesimals might be more realistic. Even though significant uncertainty was involved in the assumed strength of turbulent excitation of eccentricity in the models of Ida et al. (2008), it could be improved by more accurate measurement of the parameter, as provided by the present work.

Other interesting progress in the study of mutual collisions of planetary objects includes different outcomes from different impact angles (head-on vs. hit-and-run scenarios; see, e.g., Asphaug, 2009; Marcus et al., 2009) and improved calculations on material properties (e.g., Leinhardt & Stewart, 2009; Stewart & Leinhardt, 2009; Wettlaufer, 2010). Given all these new developments, a reconsideration of the survivability of kilometer-sized planetesimals moving in magneto-rotational turbulence seems warranted.

## 7 Summary

Directly extending our previous publication (YMM09), we continue to study massless particles moving under the gravitational influence of density fluctuations due to saturated magneto-rotational turbulence in a local, isothermal, Keplerian gas disk. We include linearized vertical gravity from the host star and thus vertical stratification of the gas disk. For comparison, the conditions in the mid-plane of the vertically stratified disks are exactly the same as those in the unstratified disks of YMM09.

In order to accurately measure the gravitational effect of the turbulent gas, we separate the gas density into two components: the basic state for the vertical hydrostatic equilibrium (Equation (4)) and the density deviation from this basic state. We use the exact gravitational acceleration due to the basic state (Equation (10)) and only solve the Poisson equation for the gravitational potential due to the density deviation (Equation (5)). We emphasize that since the Poisson equation is linear in density, this approach does not assume small density fluctuations. Furthermore, we implement isolated boundary conditions in the vertical direction and thus any density fluctuation outside the vertical computational domain is neglected.

By imposing a weak, uniform external magnetic field, we maintain a constant level of saturated magneto-rotational turbulence in the disk mid-plane. Several turbulence properties demonstrate convergence with both resolution up to 64 points per disk scale height and horizontal box size up to 16. The Shakura & Sunyaev (1973) parameter in the mid-plane of our models is controlled at the level of 10.

However, even though the properties of the turbulent gas appear numerically convergent, the dynamics of massless particles moving under the gravity of this turbulent gas does not converge with the horizontal box size . The larger the horizontal box size, the stronger the gravitational effect of the gas on the particles. Specifically, the evolution of the orbital radius, the eccentricity, and the horizontal velocity dispersion of the particles is roughly linearly dependent on up to 16. This trend was also found in our unstratified models (YMM09), and we find consistency between the unstratified models and the mid-plane of the stratified models. In contrast to the horizontal components of the particle movement, we find that the evolution of the inclination and the vertical velocity dispersion is not significantly affected by .

The dependence of particle dynamics on the horizontal box size can be traced back to the density structure of the gas. Consistent with the large-box models studied by Johansen et al. (2009), the longest wavelength Fourier mode dominates the density spectrum along the radial direction in our models. Furthermore, we find that the longest wavelength Fourier mode in the azimuthal direction also strongly influences the particle dynamics, leading to the diffusive evolution of both the orbital radius and the eccentricity of the particles. The spectral amplitudes of these longest wavelength modes are roughly constant against the horizontal box size. Using a simple single-mode analysis, we show that the linear dependence of the particle response is a natural outcome of these findings for the density spectrum.

Correlation of particle dynamics with box size poses a major difficulty for the interpretation of local-shearing-box models involving gravitational physics of magneto-rotational turbulence. We can nevertheless conjecture that , where is the distance of the box center to the host star, might be a natural scale of choice for a local model to approach reality. If this conjecture holds, we find that our previous conclusions in YMM09 on the unimportance of radial diffusive migration for protoplanets as well as the survivability of kilometer-sized planetesimals under collisional destruction may still be valid. Ultimately, high-resolution global disk models and detailed comparisons with large-box local models might be necessary to settle this issue.

We thank the anonymous referee for critical comments that significantly improved the rigor and clarity of this paper, and thank Charles Gammie, Anders Johansen, Roman Rafikov, Clément Baruteau, Martin Pessah, Richard Nelson, and Oliver Gressel for their insightful discussion on this research. This research was supported in part by the Perimeter Institute for Theoretical Physics. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. Partial support of this work was provided by the NASA Origins of Solar Systems Program under grant NNX07AI74G.

### Footnotes

1. affiliation: Present address: Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA.
2. The Pencil Code is publicly available at http://code.google.com/p/pencil-code/.
3. Note that there exists a factor of difference between our definition of disk scale height (Equation (4)) and some other authors’ in the literature.
4. In our case, we append for so the vertical domain now covers .
5. We note that use of Equation (6) would yield a 50% relative error in the vertical acceleration due to a horizontal uniform slab at a distance away (after domain doubling) when compared to the exact acceleration derived from Equation (7).
6. The bracket denotes the horizontal average of the enclosed variable at any given altitude and thus is a function of .
7. We choose this reference time to be as reported in Section 3. See Figure 1.
8. as well as other dimensionless constants introduced in the following discussions may depend on the parameter. This dependency is not investigated in this work.
9. Note that the standard deviation of a Rayleigh distribution is equal to .
10. Note that is strictly periodic in sheared coordinates and the wavenumber is slightly different from that of a direct Fourier decomposition of ; see Section 2.2.
11. It remains possible that the stochastic radial forces due to small-scale axisymmetric density perturbations drive the diffusive evolution of the particle eccentricities.
12. Note that this mode does not represent a single shearing wave, and the wavenumber is independent of time; see footnote 5.1 and Section 2.2.
13. See footnote 4.
14. This is not true if mutual gravitational scattering between planetesimals dominates over turbulence-driven excitation. However, this might not be the case in a typical protoplanetary environment; see YMM09.

### References

1. Adams, F. C., & Bloch, A. M. 2009, ApJ, 701, 1381
2. Asphaug, E. 2009, Ann. Rev. Earth Plan. Sci., 37, 413
3. Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737
4. Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1
5. Brandenburg, A. 1998, in Theory of Black Hole Accretion Disks, ed. M. A. Abramowicz et al. (Cambridge: Cambridge Univ. Press), 61
6. Brandenburg, A., Nordlund, Å., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741
7. Brucker, K. A., Isaza, J. C., Vaithianathan, T., & Collins, L. R. 2007, J. Comput. Phys., 225, 20
8. Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52
9. Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, in press (arXiv:1104.4565)
10. Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (3rd ed.; Cambridge: Cambridge Univ. Press)
11. Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125
12. Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425
13. Guan, X., & Gammie, C. F. 2011, ApJ, 728, 130
14. Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, ApJ, 694, 1010
15. Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
16. Hayashi, C. 1981, Prog. Theor. Phys. Supp., 70, 35
17. Heinemann, T., & Papaloizou, J. C. B. 2009a, MNRAS, 397, 52
18. Heinemann, T., & Papaloizou, J. C. B. 2009b, MNRAS, 397, 64
19. Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles (London: IOP)
20. Hurlburt, N. E., & Toomre, J. 1988, ApJ, 327, 920
21. Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292
22. Jamroz, B., Julien, K., & Knobloch, E. 2008, Astron. Nach., 329, 675
23. Johansen, A., Klahr, H., & Mee, A. J. 2006, MNRAS, 370, L71
24. Johansen, A., Oishi, J. S., Mac Low, M.-M., Klahr, H., Henning, T., & Youdin, A. 2007, Nature, 448, 1022
25. Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269
26. Johnson, B. M., Guan, X., & Gammie, C. F. 2008, ApJS, 177, 373
27. Johnson, E. T., Goodman, J., & Menou, K. 2006, ApJ, 647, 1413
28. Ketchum, J. A., Adams, F. C., & Bloch, A. M. 2011, ApJ, 726, 53
29. Latter, H. N., Lesaffre, P., & Balbus, S. A. 2009, MNRAS, 394, 715
30. Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489
31. Leinhardt, Z. M., & Stewart, S. T. 2009, Icarus, 199, 542
32. Lyra, W., Paardekooper, S.-J., & Mac Low, M.-M. 2010, ApJ, 715, L68
33. Marcus, R. A., Stewart, S. T., Sasselov, D., & Hernquist, L. 2009, ApJ, 700, L118
34. Nelson, R. P. 2005, A&A, 443, 1067
35. Nelson, R. P., & Gressel, O. 2010, MNRAS, 1392
36. Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849
37. Ogihara, M., Ida, S., & Morbidelli, A. 2007, Icarus, 188, 522
38. Oishi, J. S., & Mac Low, M.-M. 2009, ApJ, 704, 1239
39. Oishi, J. S., & Mac Low, M.-M. 2011, ApJ, in press (arXiv:1102.5093)
40. Oishi, J. S., Mac Low, M.-M., & Menou, K. 2007, ApJ, 670, 805
41. Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950
42. Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293
43. Papaloizou, J. C. B., & Terquem, C. 2006, Rep. Prog. Phys., 69, 119
44. Pessah, M. E. 2010, ApJ, 716, 1012
45. Proctor, M. R. E., & Weiss, N. O. 1982, Rep. Prog. Phys., 45, 1317
46. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
47. Shi, J., Krolik, J. H., & Hirose, S. 2010, ApJ, 708, 1716
48. Simon, J. B., Hawley, J. F., & Beckwith, K. 2011, ApJ, 730, 94
49. Stewart, S. T., & Leinhardt, Z. M. 2009, ApJ, 691, L133
50. Suzuki, T. K., & Inutsuka, S.-i. 2009, ApJ, 691, L49
51. Suzuki, T. K., Muto, T., & Inutsuka, S.-i. 2010, ApJ, 718, 1289
52. Vishniac, E. T. 2009, ApJ, 696, 1021
53. Wettlaufer, J. S. 2010, ApJ, 719, 540
54. Yang, C.-C., Mac Low, M.-M., & Menou, K. 2009, ApJ, 707, 1233 (YMM09)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters