DYNAMICS OF SOLIDS IN THE MIDPLANE OF PROTOPLANETARY DISKS: IMPLICATIONS FOR PLANETESIMAL FORMATION

We consider local PPD models and formulate the equations of gas and solids using the shearing sheet approximation (Goldreich & Lynden-Bell 1965). We choose a local reference frame located at a fiducial radius, corotating at the Keplerian angular velocity Ω. The dynamical equations are written using Cartesian coordinates, with $\hat{\bm {x}},\hat{\bm {y}},\hat{\bm {z}}$ denoting unit vectors pointing to the radial, azimuthal, and vertical direction, where ${\bOmega }$ is along the $\hat{\bm {z}}$-direction. The gas density and velocities are denoted by ρ_g, u in this non-inertial frame. We include a distribution of particles coupled with gas via aerodynamic drag, where the velocity of particle i is denoted by v_i. The drag force is characterized by stopping time t_stop, and equals (u − v)/t_stop per unit particle mass. Particles with different sizes have different stopping times, labeled by subscript "k." Backreaction from the particles to gas is included, which is necessary for the study of KHI and SI. In this non-inertial frame, the equations for the gas read

$$\begin{equation} \frac{\partial \rho _g}{\partial t}+\nabla \cdot (\rho _g\bm {u})=0\,, \end{equation} \tag{ 1 }$$

$$\begin{eqnarray} &&\frac{\partial \rho _g\bm {u}}{\partial t} + \nabla \cdot (\rho _g\bm {uu}+P\bm {I}) \nonumber\\ &&\quad = \rho _g\bigg [2{\bm u}\times {\bOmega }+3\Omega ^2x\hat{\bm {x}}-\Omega ^2z\hat{\bm {z}}+ \sum _{k}\epsilon _k\frac{\overline{\bm {v}}_k-\bm {u}}{t_{{\rm stop},k}}\bigg ], \ \end{eqnarray} \tag{ 2 }$$

where the source terms include Coriolis force, radial tidal potential as well as disk vertical gravity. The last term in the momentum equation represents the backreaction (or momentum feedback) from particles to gas: _k and $\overline{{\bm v}_k}$ denote the local mass density and velocity of particles of type k. In this paper, we neglect the effect of magnetic fields and focus on the interaction between gas and solids in the dead zone of PPDs (Gammie 1996). An isothermal equation of state for the gas is used throughout this paper, where P = ρ_gc²_s and c_s is the isothermal sound speed.

Similarly, the equation of motion for particle i of type k can be written as

$$\begin{equation} \frac{d\bm {v}_i}{d t}=-2\eta v_K\Omega \hat{\bm {x}}+2{\bm v}_i\times {\bOmega }+3\Omega ^2x_i\hat{\bm {x}} -\Omega ^2z_i\hat{\bm {z}}-\frac{\bm {v}_i-\bm {u}}{t_{{\rm stop},k}}.\ \ \end{equation} \tag{ 3 }$$

In the above equation, we have added an inward force term $-2\eta v_K\Omega \hat{\bm {x}}$ to mimic the effect of an outward radial pressure gradient in the gas (Bai & Stone 2010a), where ηv_K is the difference between gas velocity and the Keplerian velocity in the absence of particles. This term will shift both gas and particle azimuthal velocities by ηv_K relative to those in the real system. To avoid confusion, we always use u and v to denote velocities that correspond to the real system (i.e., subtracting the azimuthal velocity component from the simulation by ηv_K). Particle self-gravity is ignored as we focus on the dynamics in the midplane of the PPD dead zone and precursor of the planetesimal formation.

In our simulations, we have applied an orbital advection algorithm for both gas and particles (Stone & Gardiner 2010; Bai & Stone 2010a), and the actual velocities used in the simulation are measured relative to the linearized Keplerian shear flow: ${\bm u}^{\prime }={\bm u}+(3/2)\Omega x\hat{\bm {y}}$ for gas flow, and ${\bm v}^{\prime }_i={\bm v}_i+(3/2)\Omega x\hat{\bm {y}}$ for individual particles.

Measuring velocities in units of the sound speed, time in units of Ω⁻¹, and length in units of the gas scale height H_g ≡ c_s/Ω, the parameters in the problem are reduced to the following.

• 1.  
  The dimensionless particle stopping time τ_k ≡ Ωt_stop,k for particle species k.
• 2.  
  The solid abundance parameter Z_k for each particle species, which measures the height-integrated particle-to-gas mass ratio.
• 3.  
  The parameter characterizing the strength of the radial pressure gradient Π ≡ ηr/H_g = ηv_K/c_s.

Below, we apply a disk model and provide the scaling relation between the disk model parameters and these dimensionless parameters used in our simulation.

We adopt a generalized solar nebular model where the disk is vertically isothermal and all the disk quantities have a power-law dependence on the radius (Youdin & Shu 2002)

$$\begin{eqnarray} \Sigma _g &=& 1700f_gr_{\rm AU}^{-b}{\rm \ g\ cm}^{-2}\,,\nonumber\\ T &=& 280f_Tr_{\rm AU}^{-c}{\rm \ K}\,,\\ M_* &=& f_MM_{{\odot }}\,,\nonumber \end{eqnarray} \tag{ 4 }$$

where Σ_g is the gas surface mass density, T is the disk temperature, M_* is the mass of the central star, and r_AU ≡ r/1 AU. These parameters fix the disk model. Although the global disk profile may not follow the simple power-law form, we can always approximate a local patch of the disk in the above form, which is very general. In the standard minimum-mass solar nebular (MMSN) model (Hayashi 1981), we have b = 3/2,  c = 1/2,  f_T = f_g = f_M = 1. The radial profiles of other physical quantities are

$$\begin{eqnarray} \Omega &=& 2\pi f_M^{1/2}r_{\rm AU}^{-3/2}{\rm \ yr}^{-1}\,,\nonumber\\ v_K &=& 30f_M^{1/2}r_{\rm AU}^{-1/2}{\rm \ km\ s}^{-1}\,,\\ c_s &=& f_T^{1/2}r_{\rm AU}^{-c/2}{\rm \ km\ s}^{-1}\,,\nonumber\\ H_g &=& 3.4\times 10^{-2}f_T^{1/2}f_M^{-1/2}r_{\rm AU}^{(3-c)/2}{\rm AU}\,, \nonumber \end{eqnarray} \tag{ 5 }$$

where in the calculation of the sound speed, we assume the mean molecular weight μ = 2.33.

The background gas density profile is

$$\begin{equation} \rho _{g,b}(r,z)=\frac{\Sigma _g}{\sqrt{2\pi }H_g}\exp \big({-}z^2/2H_g^2\big)\,, \end{equation} \tag{ 6 }$$

where subscript "b" denotes "background." Using this gas density and sound speed, one can derive the radial pressure gradient in the gaseous disk, thus obtain the amount of reduction ηv_K in the gas rotation velocity. After some algebra, we can derive the pressure length scale parameter

$$\begin{eqnarray} \Pi \approx -\frac{1}{2}\frac{d\ln {P}}{d\ln {r}}\frac{c_s}{v_K}&=&\bigg (\frac{3+2b+c}{4}- \frac{3-c}{4}\frac{z^2}{H_g^2}\bigg)\frac{c_s}{v_K}\nonumber\\ &\approx & 0.054f_T^{1/2}f_M^{-1/2}r_{\rm AU}^{1/4}. \end{eqnarray} \tag{ 7 }$$

Note that Π = Π(r, z) depends on both radius and height. Nevertheless, in this paper, our simulation box is concentrated in the disk midplane where z ≪ H_g, therefore we can neglect the dependence of Π on z. In the last equation of the above formula, we have applied the power-law indices of the MMSN model. The dependence on disk temperature f_T, stellar mass f_M as well as disk radius r is relatively weak. It is worth mentioning that the dependence of Π on disk mass is only through the surface density profile parameter b, free from the scaling parameter f_g. Therefore, the value Π ≈ 0.05 should apply to a wide range of disk models.

Next we consider the scaling relations for the dimensionless stopping time. Because the gas motion in PPDs is expected to be subsonic, the relevant drag laws from the gas to the solids in PPDs are the Epstein drag law (Epstein 1924), which applies when particle size is smaller than the gas mean free path, and the Stokes drag law, which applies for larger bodies. We assume all solid bodies have spherical shapes, then the stopping time in these two regimes can be expressed as (Weidenschilling 1977)

$$\begin{equation} t_{\rm stop}=\left\lbrace \begin{array}{@{}l@{\quad }l@{}} {\displaystyle {\rho _sa \over \rho _gc_s}}\, & \mbox{$a<9\lambda _m/4$ (Epstein regime)}, \\ \ms {\displaystyle {4\rho _sa^2 \over 9\rho _gc_s\lambda _m}}\, & \mbox{$a>9\lambda _m/4$ (Stokes regime)}, \\ \end{array}\right. \end{equation} \tag{ 8 }$$

where ρ_s ≈ 3 g cm⁻³ and a are the density and radius of the solid body, λ_m = (n_gσ)⁻¹ = μm_H/ρ_gσ is the mean free path of the gas, and σ ≈ 2 × 10⁻¹⁵ cm² is the molecular collision cross section (Chapman & Cowling 1970). From the above equations, we see that the particle stopping time depends linearly on gas density in the Epstein regime. Nevertheless, the gas density can be regarded as constant near the disk midplane where we study. Therefore, in our local simulations, we can safely take t_stop as depending on particle size a only.

To better handle the relation between particle size and its corresponding stopping time, we express the relation between τ_s ≡ Ωt_stop and a by applying our disk model. The result is

$$\begin{eqnarray} \tau _s &=&\max \big [ 4.4\times 10^{-3}a_{\rm cm}f_g^{-1}r_{\rm AU}^b,\nonumber\\ && 1.4\times 10^{-3}a_{\rm cm}^2f_T^{-1/2}f_Mr_{\rm AU}^{(c-3)/2}\big ], \end{eqnarray} \tag{ 9 }$$

where a_cm is the particle radius measure in centimeter. In the MMSN model, at 1 AU, particles smaller than 3 cm are in the Epstein regime. At larger radii, the Epstein regime applies to much larger particles.

Fiducially, we consider the MMSN model at 1 AU, and set the pressure length scale parameter Π = 0.05. This parameter is kept fixed in all our simulations. Instead of considering a particle size distribution in radius, we consider the distribution in τ_s. Then one can easily translate it into particle radius given the parameters of the disk model. We discretize a continuous particle size distribution into a number of bins. Each bin covers half a dex in τ_s in the logarithmic scale. For simplicity, we assume a uniform particle mass distribution across the bins, that is, all the particle bins (or particle species) have equal amount of mass. The parameters for the size distribution are therefore the minimum and maximum dimensionless stopping time τ_min and τ_max (translated to a_min and a_max, respectively). Physically, our assumption means that most of the mass of the solids resides in the size range between a_min and a_max and roughly follows a flat distribution in logarithmic scale. To control the total particle mass, we use the total solid abundance parameter

$$\begin{equation} Z=\sum _{k=1}^{N_{\rm type}}Z_k\,,\quad {\rm with}\quad Z_k=Z/N_{\rm type}\,, \end{equation} \tag{ 10 }$$

where N_type is the number of particle types (bins). Currently, the best estimate of the solar metallicity is about 0.015 (Lodders 2003). A substantial fraction of the metal elements may reside in dust grains and grow into larger bodies. In our simulations, we consider three abundance values Z = 0.01, 0.02, and 0.03. This choice covers a relatively wide range of disk metallicities. Moreover, because our simulation focuses on a local patch in a PPD, the local abundance may not necessarily be equal to the averaged value in the PPD.

As we explained in Section 1, we perform simulations in both 2D and 3D. Our 2D simulations are in the $\hat{\bm x}$–$\hat{\bm z}$ plane (i.e., axisymmetric). Details of the implementation and code tests of the particle–gas hybrid scheme are given in Bai & Stone (2010a). Our simulations use the standard shearing box approach (Hawley et al. 1995), where the radial boundary condition is periodic with azimuthal shear. Azimuthal boundary conditions are periodic. Vertical gravity is included in our simulations, and we choose reflection boundary condition in the $\hat{\bm z}$-direction, which is the same as that in Johansen et al. (2009). In general, we use 256 cells in the radial (and azimuthal, if applicable) direction. Guided by Bai & Stone (2010a), properly resolving the SI with τ_s = 0.1 requires about 128 cells per pressure length scale ηr. With this required resolution, our simulation box size is typically small, spanning only about 2–4ηr. Such small box size is also necessary to capture the typical wavelength of the KHI, if present (Johansen et al. 2006a). In our simulations, we generally use N_p = 65536 particles per type for 2D simulations and N_p = 3145728 particles per type in 3D runs (in which cases N_type = 7). Larger N_p are used when N_type is smaller to keep the total number of particles similar in all our simulations. Our choice of particle number guarantees at least one particle per cell per particle type around the disk midplane, as required for numerical convergence (Bai & Stone 2010a).

In our simulations, we set the initial particle density profile to be a Gaussian centered on disk midplane with scale height H_d = 0.015 H_g for all particle types. The particle and gas velocities are computed from a multi-species Nakagawa–Sekiya–Hayashi (NSH) equilibrium, where the classical single-species NSH equilibrium (Nakagawa et al. 1986) solution is generalized to include multiple species of particles (see Appendix A). Note that different particles have different velocities, and the velocities of particles and gas depend on z.

The choice of our simulation box size and boundary conditions in the vertical direction merit further discussion. In the simulations, gas–particle interaction in the disk midplane generates turbulence and excites vertical motions in the gas. Ideally, the vertical box size should extend to a few H_g, similar to what is used for MRI simulations (e.g., Stone et al. 1996), however, this would make 3D simulations too expensive. We have conducted a series of tests in 2D with a single particle species τ_s = 1 using different vertical box sizes and either reflecting or periodic boundary conditions. In both cases, particles settle to the disk midplane with a spatial distribution reminiscent of sinusoidal waves that slowly drift in the radial direction. We find that the particle scale height is more intermittent when using periodic boundary conditions. Moreover, periodic boundary conditions appear to suppress asymmetric modes in the gas azimuthal velocity around the disk midplane. Using reflecting boundary conditions, we find essentially no difference between the particle scale heights and clumping properties obtained from different vertical box sizes once the box height is much larger than the particle scale height, although it takes longer for the system to reach a quasi-steady state when a larger vertical box size is used. The drift velocities of the wave-like pattern of particles do differ when different vertical box sizes are used, but they are unlikely to affect the properties discussed in Sections 3–6.⁴ Guided by these results, as long as the vertical boundary of our simulation box is well above the scale height of all particle species, one should get converged results from the simulations.

Table 1 lists the parameters of all of our simulations. Our runs are labeled using names with the form RxyZz–nD, where x, y are integers corresponding to τ_min = 10^−x, and τ_max = 10^−y, z ≡ 100Z represents the solid abundance, and n = 2 (n = 3) denotes 2D (3D) simulations. When referring to simulations with fixed x and y but all possible values of z and/or n, we omit the Zz, and/or the nD, parts of the names. We focus on two groups of runs. In the first group, the maximum particle stopping time is τ_max = 0.1. We use seven particle species to span three orders of magnitude in stopping time (down to τ_min = 10⁻⁴) for the series of runs labeled R41, while in the series labeled R21, we use three particle species to span one order of magnitude in stopping time (down to τ_min = 10⁻²). In the second group of runs, the maximum particle stopping time is τ_max = 1.0, and the minimum stopping time is chosen to be τ_min = 10⁻³ (R30) or 0.1 (R10). In each series of runs (R41, R21, R30, R10), we perform three 2D simulations with Z = 0.01, 0.02, and 0.03, and two 3D simulations with Z = 0.01 and Z = 0.03. Because of a smaller τ_max in the first group, higher resolution is needed to resolve the SI.

To determine when a saturated state is reached in each simulation, we monitor the particle vertical scale height H_p,k for each particle species k, defined as the rms value of the z coordinate of all particles. Saturation occurs when particle settling and turbulent diffusion are in balance, so that the scale height of all particle species is steady. In Figure 1, we show the time evolution of the vertical scale height for each particle species (marked by different colors) in all our runs. Solid and dashed curves represent 2D and 3D simulations, respectively. We see that most of the 2D runs saturate within about 50 orbits.⁵ The 3D simulations are very time consuming, so we run them for shorter periods. From Figure 1, all 3D runs saturate before we terminate the simulations, although some just barely so.

Zoom In Zoom Out Reset image size

In the last column of Table 1, we provide the time of saturation T_s (in parentheses) for each simulation. Unless otherwise stated, we will perform data analysis in the time interval between the saturation time T_s and the end time of the simulation T_e. In the R41Z3-2D and R21Z3-2D runs, there are sudden jumps in particle heights followed by settling, and this process repeats over time quasi-periodically. Averaging over many cycles is required to reduce the influence of these intermittent "bursts." The vertical distribution of the smallest particles in the R41-3D and R21-3D runs with Z = 0.03 are not fully saturated at the end of our simulations. Nonetheless, the scale heights of the largest particles (which dominate the dynamics) in these runs have reached steady state, therefore we consider them to be saturated.

Before presenting a detailed data analysis, we show the distribution of particles at the end of our simulations in Figure 2. Results from 3D runs are shown by projecting particle positions in three orthogonal directions. The number of particles plotted is much less than the actual number of particles used in the simulation. The trends in particle scale height evident in Figure 1 can be clearly seen: particles with small τ_s are diffused to larger heights. Note that we overplot larger particles on top of small particles, so that small particles near the midplane are less visible. The SI is present in all the simulations, and we will discuss various aspects of Figure 2 in the following sections.

Zoom In Zoom Out Reset image size

The source of turbulence responsible for stirring up the particles can in principle be due to both KHI and SI. It is important to decipher which instability is the dominant process. Generally speaking, the onset of SI requires the averaged particle-to-gas mass ratio ≳ 1. The strength of the instability decreases as the averaged particle size becomes smaller, and vanishes as τ_s → 0, for which the dust and gas behave as a single fluid. The onset of KHI requires a steep vertical profile of the gas azimuthal velocity, which generally corresponds to larger dust-to-gas mass ratio at disk midplane. In our simulations, a substantial fraction of the particles have a relatively large stopping time with τ_s > 10⁻², and SI clearly plays an important role in the generation of disk midplane turbulence. It remains to study whether KHI is present and whether KHI is dynamically important.

The classical result on the onset of KHI in a vertically stratified disk is based on the Richardson number criteria (Chandrasekhar 1961)

$$\begin{equation} {{\rm Ri}}_{x,y}\equiv \frac{g}{\rho }\ \frac{(\partial \rho /\partial z)}{(\partial u_{x,y}/\partial z)^2}\,, \end{equation} \tag{ 11 }$$

where we define the Richardson number from radial and azimuthal velocity shear, as indicated by subscripts x, y. In the above equation, g = Ω²z is the vertical gravitational acceleration, and ρ is the effective fluid density (see discussion below). The Richardson number measures the amount of work required to overturn the fluid (numerator) in comparison to the amount of free energy available in the vertical shear (denominator). For Cartesian flow with no rotation, the necessary condition for instability is given by Ri < Ri^crit = 1/4. This criteria no longer hold when rotation (Coriolis force) and radial shear (differential rotation) are included, especially when the rotation frequency Ω is comparable to the Brunt–Väisälä frequency of buoyant oscillations. Generally speaking, the Coriolis force distablizes the fluid (Gómez & Ostriker 2005), while radial shear acts to stablize the fluid. Lee et al. (2010a) found that Ri^crit is typically smaller than 1/4 and is roughly proportional to dust-to-gas mass ratio at disk midplane. In this paper, we adopt the critical Richardson number to be Ri^crit = 0.1 as suggested by Chiang (2008).

The Richardson number criterion is based on a single fluid, in which case ρ simply represents fluid density. With the addition of perfectly coupled dust, the dust–gas system behaves as a single fluid, where the dust contributes to the mass but not the pressure of the fluid, thus ρ = ρ_g + ρ_p. When particles are not perfectly coupled, the definition of ρ becomes somewhat ambiguous, but we expect ρ_g < ρ < ρ_g + ρ_p. Below we provide a simple formula for ρ in this regime that reduces the above two limiting cases when ρ_p = 0 and when τ_s → 0.

In the absence of any turbulence and vertical gravity, the equilibrium state between gas and dust (with fixed stopping time) is described by the NSH solution (Nakagawa et al. 1986). In particular, the azimuthal gas velocity relative to Keplerian velocity is given by

$$\begin{equation} u^{\prime }_y=-\bigg [1-\frac{\epsilon (1+\epsilon)}{(1+\epsilon)^2+\tau _s^2}\bigg ] \eta v_K\,, \end{equation} \tag{ 12 }$$

where = ρ_p/ρ_g, and prime means Keplerian velocity is subtracted. For convenience, we define Δu_y ≡ −u'_y. For perfectly coupled particles, τ_s = 0 and we find ρ_g + ρ_p = ρ_gηv_K/Δu_y. For particles with finite stopping time, Δu_y becomes closer to ηv_K, which reflects the fact that the particle–gas coupling is weaker so that gas velocity shifts toward the dust-free value. Therefore, Δu_y can be regarded as an indicator of particle–gas coupling. In this spirit, we define the effective gas density as

$$\begin{equation} \rho ^{\rm eff}\equiv \rho _g\frac{\eta v_K}{\Delta u_y}. \end{equation} \tag{ 13 }$$

It is trivial to check that in the limit → 0, ρ^eff → ρ_g, and ρ^eff → ρ_g + ρ_p when τ_s → 0. In the calculation of the Richardson number, we substitute ρ by ρ^eff. Since ρ_g is nearly constant over the height of our simulation box, Equation (11) becomes

$$\begin{equation} {{\rm Ri}}_{x,y}=-\frac{\Omega ^2z}{\Delta \overline{u}_y}\frac{(\partial \overline{u}_y/\partial z)}{(\partial \overline{u}_{x,y}/\partial z)^2}\,, \end{equation} \tag{ 14 }$$

where the overbar means averaging over the horizontal plane. Note that Ri depends on z.

Before calculating the Richardson number profile from our simulations, we first return to the spatial distribution of particles in Figure 2. In 2D simulations, we see that the distribution of particles around the disk midplane is highly non-uniform, and exhibits wave patterns in the x–z plane that are almost stationary over time. Results from 3D simulations show very similar features in the x–z plane. In particular, in runs R30Z1-3D and R10Z1-3D, there is a clear segregation of particles with different stopping times, and their wave patterns have a phase shift relative to each other. However, in the y–z plane, there is no coherent structure in the projected distribution of particles in any of our 3D simulations. This contrasts with the expectations from the KHI, where the particle layer kinks and breaks into clumps (Johansen et al. 2006a; Barranco 2009). Based on this observation, we infer that in our 3D simulations, KHI is not present in the azimuthal direction. Moreover, in the x–y plane, we see azimuthally elongated stripes of the large particles (in black). This feature, together with the standing wave structure in the x–z plane, is most likely to be due to SI. KHI resulting from the vertical shear in the gas radial velocity is another possibility, however, we have found that Ri_x is always larger than Ri_y from our simulations, therefore the KHI is unlikely to play a role in the simulations presented here.

In Figure 3, we show the Richardson number profile associated from u_y calculated from the saturated states of all our simulations. The Richardson number is generally smallest in the disk midplane, and increases with height. In almost all our 3D simulations (dashed curves), Ri_y is greater than the critical value (0.1), therefore, the dusty midplane layer is expected to be stable against vertical shear, consistent with the spatial distribution of particles discussed above. Given the fact that Ri does not solely determine stability, this observation does not entirely exclude the possibility that Ri_y could be maintained by KHI. However, it is important to note that KHI is suppressed in 2D. We see that Ri_y from all our 2D simulations (solid curves) are generally close to their 3D counterpart. This means that the SI itself is able to maintain Ri above the critical value, and suggests that the KHI is indeed absent in all our simulations.

Zoom In Zoom Out Reset image size

The main reason that we do not observe KHI is that the turbulence generated from the SI is strong enough to prevent particles from settling sufficiently to trigger KHI. We note that the strength of the SI turbulence decreases as the particle stopping time τ_s decreases (as expected from the linear analysis of Youdin & Goodman 2005, and as confirmed by our numerical experiments). The turbulence in our simulations is mainly generated from relatively large particles with τ_s ≳ 0.01 (see also the next subsection). We have not explored the regime where all particles are strongly coupled to the gas. However, in this regime, we expect the SI to be generated on much smaller spatial scales with much lower amplitude, so that the particles settle until the KHI is triggered. In this regime, the dust–gas system behaves as a single fluid, where the dust contributes to the mass density but not the pressure of the fluid. This is the approach adopted by Chiang (2008), Barranco (2009), and Lee et al. (2010a, 2010b) to study the KHI.

Figure 4 shows the vertical density profiles for particles of different types from all our 3D simulations, calculated by binning the particles into vertical grid cells and averaging over time after saturation. Results from 2D simulations are generally similar.

Zoom In Zoom Out Reset image size

The vertical density profile of particles is determined by the balance between particle settling and turbulent diffusion. Unlike studies of passive particles under the influence of homogeneous external turbulence (Youdin & Lithwick 2007), the turbulence from our simulations is self-generated, and is non-homogeneous (strongest at the disk midplane). To study the properties of turbulent diffusion, one approach would be to assume some functional form for the vertical profile of the diffusion coefficient D_g,z(z), and fit the particle density profiles. However, after several experiments we found it difficult to fit the density profile of all particle species simultaneously with any simple functional form of D_g,z(z).⁶ In fact, the wave patterns in the x–z plane shown in Figure 2 suggest that the classical turbulent diffusion scenario may be too simple.

Instead of fitting the vertical profile of the turbulent diffusion coefficient in the gas, we pose the question in another way: what is the effective vertical diffusion coefficient at the disk midplane for the particles that are driving the turbulence? Since we have identified the SI as the source of the midplane turbulence, one expects particles with relatively large stopping times to drive the turbulence both from a theoretical point of view (Youdin & Goodman 2005) and from non-stratified simulations of SI (Johansen & Youdin 2007, Bai & Stone, unpublished). To address these questions more quantitatively, we find the following approach particularly useful.

We fit the horizontally averaged vertical density profile of the largest particles τ_s = τ_max in each simulation using the classical picture of turbulent diffusion. Since these particles (as well as particles with slightly smaller τ_s) actively drive the disk turbulence, the gas turbulent diffusion coefficient across this particle layer can be regarded as constant. Therefore, the vertical density profile of these particles is expected to be Gaussian, with scale height (Youdin & Lithwick 2007)

$$\begin{equation} H_p(\tau _s)=\sqrt{\frac{D_{g,z}(0)}{\Omega \tau _s}} \sqrt{\frac{\tau _s+\tau _e}{\tau _s+\tau _e+\tau _s\tau _e^2}}\,, \end{equation} \tag{ 15 }$$

where τ_e = Ωt_eddy is the turnover time of the largest eddies. The basic assumption behind this formula is stochastic turbulent forcing on passive particles with the autocorrelation function of the turbulence P(t) = exp(−t/t_eddy)/2π, corresponding to a Kolmogorov spectrum. We do not have much knowledge of t_eddy for SI turbulence, but expect it to be comparable to the orbital time (the only timescale of the problem), and take τ_e = 1. The exact value of τ_e does not matter much, since it only gives an order unity correction to H_p.

By fitting the vertical density profile of the largest particles with a Gaussian, we obtain D_g,z(0) for all simulations, and the results are summarized in Table 2. For 3D runs, the results are also plotted in Figure 4 as dashed lines. We see that the vertical profiles of the largest particles are well fitted with a Gaussian. In addition, we predict the vertical density profile for other particle species, using Equation (15) and assuming a diffusion coefficient which is constant with height. Obviously, this will overpredict the scale heights for small particles, since they respond to the turbulence passively. However, for particles that actively participate in the instability, we expect their density profile to be comparable to the predicted profile, since they are driving turbulence to maintain D_g,z close to D_g,z(0) across their scale heights. In this way, we are able to identify the particle species that are responsible for the disk turbulence (hereafter termed as "active" particles).

From the R41 runs, we see that active particles range from τ_s = 0.1 (for R41Z1) to τ_s ≳ 0.01 (R41Z3). Active particles for R21 runs have τ_s ≳ 0.03. For R30 runs, particles with τ_s ≳ 0.03 are active, while for R10 runs, all particles are active. We see that although there is a diversity in the size range of active particles, which depends on both solid abundance and particle size distribution, the minimum size of active particles for most of our runs is about τ_s = 0.01–0.03. For run R41Z1, although we have identified somewhat larger τ_s values for active particles, particles with τ_s = 0.01–0.1 must actively participate in the instability because the abundance of τ_s = 0.1 particles alone is too small to trigger SI.

Next we study the midplane diffusion coefficient from our simulations. We emphasize that the strength of the turbulence (hence D_g,z) is self-regulated: the settling of particles continues until the turbulence they generate is sufficient to stop the settling. To better interpret our results, we construct a toy model describing the self-regulated turbulence. In this model, we assume all particles are active, and that the particles are single sized, with fixed stopping time τ_s. Since all particles are active, their vertical density profile can be approximated by a Gaussian, so that the particle-to-gas mass ratio at the disk midplane is given by = ZH_g/H_p. The midplane diffusion coefficient D_g,z depends on both τ_s and . For simplicity, we parameterize the dependence as

$$\begin{equation} D_{g,z}=\epsilon ^\alpha f(\tau _s)H_g^2\Omega \,, \end{equation} \tag{ 16 }$$

where D is normalized to H²_gΩ, f(τ_s) is a coefficient that incorporates the dependence of D on τ_s, and α is a power-law index that reflects the sensitivity of the dependence of D on . We note that D_g,z → 0 at both → 0 and → ∞, therefore, we expect α>0 when is small and α < 0 for large . Using Equation (15) and neglecting the second square root (which is order unity) on the right-hand side, we obtain

$$\begin{eqnarray} H_p &=&\bigg (\frac{f}{\tau _s}Z^\alpha \bigg)^{\frac{1}{2+\alpha }}\cdot H_g\nonumber\\ D_{g,z} &=&(fZ^\alpha)^{\frac{2}{2+\alpha }} \tau _s^{\frac{\alpha }{2+\alpha }}\cdot H_g^2\Omega. \end{eqnarray} \tag{ 17 }$$

In the above equations, the dependence of particle scale height and diffusion coefficient on Z is reflected in the index α. When α is positive, increasing Z leads to larger H_p and larger D_g,z. When α is negative, the situation reverses. Below, we apply this simple model to our results. Since our simulations contain multiple particle species, we may take to represent the contribution from all particle species participating in the SI (i.e., with τ_s ≳ 0.01).

Our R30 and R10 runs show similar behavior between 2D and 3D simulations with respect to vertical diffusion properties. Increasing Z from 0.01 to 0.02 produces stronger turbulence, while further increasing Z to 0.03 dramatically reduces H_p. This corresponds to the transition from α>0 to α < 0 at a threshold (hence threshold Z = Z_th). Beyond Z_th, H_p sensitively depends on Z because the corresponding power-law index α/(2 + α) quickly drops to large negative values once α turns negative. Consequently, a small increase in Z results in strong particle settling and greatly enhances midplane particle density. This result has important implications for particle clumping discussed in the next section.

In our 2D R41 and R21 runs, we see that D_g,z monotonically decreases with Z, suggesting α < 0 for Z ⩾ 0.01. Based on this result, we infer that the strength of the SI for a particle size range τ_s = 0.01–0.1 is a decreasing function of for ≳ 0.5. The 3D simulations give somewhat different results. For both 3D R41 and R21 runs, D_g,z slightly increases with Z at least in the range Z ⩽ 0.03, indicating α ⩾ 0. It is very likely that the threshold abundance Z_th is above 0.03, which is substantially larger than their 2D counterparts. We note that the behavior of the SI turbulence for τ_s ≲ 0.1 particles in 3D is different from that in 2D in non-stratified simulations (Johansen & Youdin 2007). Our results indicate that the difference remains when vertical gravity is included, and 3D simulations are needed to better catch the dynamics of small particles.

In our toy model, all of our ignorance on the dependence of D_g,z on τ_s is encapsulated in the unknown function f(τ_s). From Table 2, we see that the R30 and R10 runs generally have larger D_g,z than R41 and R21 runs. This result implies that turbulence generated from larger particles τ_s ∼ 1 is stronger than that from smaller particles, i.e., f(τ_s) is an increasing function of τ_s in this range, consistent with results from non-stratified simulations (Johansen & Youdin 2007).

In sum, we have identified that particles actively participating in SI generally have stopping time τ_s ≳ 0.01. The strength of the turbulence largely depends on the density of these active particles at disk midplane. We find that the particle scale height (thus the turbulent diffusion coefficient) strongly depends on solid abundance. Such strong dependence is caused by a sharp drop in the strength of the turbulence with increasing particle-to-gas mass ratio when is larger than a certain threshold value.

Probably the most interesting property of the SI is the concentration of particles. The degree of particle concentration strongly depends on the mass distribution of solids in PPDs. In our simulations, we normalize particle density to the background gas density at the disk midplane ρ_g,b(r, z = 0). A useful scale to measure particle concentration is the Roche density, above which the particle clump can be considered as gravitationally bound (Binney & Tremaine 2008)

$$\begin{equation} \rho _{\rm roche}(r)\approx \frac{3M_*}{r^3}=1.34\times 10^3\frac{(f_Mf_T)^{1/2}}{f_g}r_{\rm AU}^{\frac{2b-c-3}{2}}\rho _{g,b}(r,0). \end{equation} \tag{ 18 }$$

The normalized Roche density (relative to the background gas density at midplane) scales as the square root of stellar mass and disk temperature, and is inversely proportional to disk mass, meaning that the Roche density is easier to reach for massive disks (with large f_g). In the MMSN model, the Roche density is of the order ρ_roche = 10³ρ_g,b, and only weakly depends on r as r^−1/4.

In Figure 5, we show the time evolution of maximum particle density ρ_p,max from all our simulations. We first look at results from 2D simulations. For all the four run series, ρ_p,max increases with solid abundance Z. However, the dependence of ρ_p,max on Z is highly nonlinear. For run series R21, R30, and R10, there is no significant clumping of particles for Z = 0.01 and 0.02. However, significant clumping occurs at Z = 0.03, with maximum particle density reaching 10³ times the background gas density, comparable to the Roche density (Equation (18)). This trend is consistent with the results by Johansen et al. (2009) (see also the supplemental information in Johansen et al. 2007), who considered particles with stopping time in the range of τ_s = 0.1–0.4. As emphasized in the previous section, there is a sharp enhancement of averaged midplane particle density with increasing Z once Z exceeds some threshold value. This density enhancement further favors strong concentration of particles by SI, which explains the trend we have observed in Figure 5.

Zoom In Zoom Out Reset image size

The particle clumping also depends on the particle size distribution. In the R41 run series, where the majority of the particle mass resides in strongly coupled particles τ_s < 10⁻², we see that there is no significant clumping of particles up to Z = 0.03. As noted in the previous section, particles that effectively participate in SI are those with relatively large stopping times τ_s ≳ 10⁻². These particles are also the ones that actively participate in the clumping (see the next subsection). For R41 runs, the abundance of these "active" particles is much smaller than our R21, R30, and R10 runs, which makes the critical (total) abundance for strong particle clumping larger. In fact, we do observe strong clumping when we increase the total abundance to Z = 0.05. Based on the discussion above, we conclude that in order for the SI to efficiently concentrate particles, the mass of the solids with stopping time τ_s ≳ 10⁻² should exceed a critical value Z_crit. The results from 2D simulations suggest that $\sum _{\tau _k\geqslant 10^{-2}}Z_k \gtrsim Z_{\rm crit}\approx 0.02$ is necessary for significant particle clumping.⁷

The 3D simulations show similar trends as in 2D, but the condition for strong particle clumping is more stringent. Among the eight 3D runs, strong clumping occurs only in run R10Z3-3D. The maximum density for all other runs remains small in the saturated state (ρ_p,max ≲ 50ρ_g,b). In particular, the 3D R21Z3 and R30Z3 runs do not show clumping as in their 2D counterparts, and both of them have larger D_g,z. Since KHI is unlikely to be present in these simulations, the different results between our 2D and 3D simulations should be attributed to the different behavior of the SI in 2D and 3D. It appears that the formation of dense particle clumps favors the mass distribution of particles to be dominated by larger particles than in 2D, or larger values of Z_crit is needed.

Interestingly, in run R30Z3-3D, a very dense clump (actually a nearly axisymmetric stripe) forms at about t = 150 Ω⁻¹. The composition of this (transient) clump is similar to its counterpart R30Z3-2D (see next subsection). It lasts for about 10 orbital times and then is gradually dissolved. Both the Richardson number profile and particle distribution disfavor the presence of KHI during the process. Nor is there any significant vorticity generation in the vicinity of the clump which might indicate KHI. By comparing with Figure 1, we see that the period during which the clump is dissolved is accompanied by an increase of the height of relatively small particles with τ_s ≲ 0.1. It is likely that the formation of the transient clump is due to our unrealistic initial condition.⁸

The results we have obtained show a clear dichotomy on the particle concentration properties. Specifically, the maximum density is either very small with ρ_p,max ≲ 50ρ_g,b, or very large with ρ_p,max ≳ 1000ρ_g,b. Self-gravity becomes important when the particle density approaches the Roche density (Equation (18)). This means that for our simulations that do not show signature of strong clumping, adding self-gravity will not change the picture qualitatively.⁹ For simulations with strong clumping, the maximum particle density is already comparable with the Roche density, and in this case we expect the formation of a few planetesimals from the simulations as in Johansen et al. (2009).

Particle concentration properties are known to depend on numerical resolution. To assess the validity of our results, we have also performed the same set of simulations with half our standard resolution. We find the same dichotomy between strong clumping and no clumping. The only exception is the R30Z3-3D run: it shows strong particle clumping in the low-resolution run which does not dissolve as in our standard resolution run. The reason is that the turbulence generated from the lower resolution run is weaker, thus particles settle more which favors clumping. This test justifies the necessity of conducting high-resolution simulations. In the mean time, it suggests that the critical abundance for particle clumping in this run may be only slightly larger than 0.03. Therefore, the particle clumping properties from 2D and 3D simulations is not dramatically different when τ_max = 1.

In this subsection, we discuss more details of the three simulations that exhibit strong particle clumping. First, we examine the composition of these dense clumps by plotting the cumulative probability distribution function (CPDF) of particle densities for different particle species P(ρ_p > ρ). The CPDF measures the probability of a particle residing in a region with total particle density larger than ρ. In Figure 6, we plot the CPDFs of the three runs: R21Z3, R30Z3, and R10Z3. At relatively high densities with ρ_p ≳ 10²ρ_g, we see that in all three cases, the dense regions are composed of particles with the largest stopping times. In run R21Z3-2D, the mass fraction of different particle species in the dense clumps is increasing with particle stopping time τ_s, and is completely dominated by the largest particles τ_s = 0.1. In the case of R30Z3-2D and R10Z3-2D, where the largest particles have τ_s = 1.0, the composition of the clumps is dominated by the two largest particle species. Contribution from other particle species to the clumps is almost negligible by mass.

Zoom In Zoom Out Reset image size

For R21Z3 and R30Z3 runs, 3D simulations do not show particle clumping, therefore, the resulting CPDFs differ substantially from those in 2D runs. Nevertheless, these CPDFs provide typical examples for simulations without clumping. The shapes of the CPDFs from different particles are very similar, and curves for larger particles are located to the right of those for smaller particles, consistent with the vertical stratification of particles. For run R10Z3-3D, the particle clumping is stronger than the 2D case, and the densest clumps are almost equally made of particles with τ_s = 1 and τ_s = 10^−1/2.

Next, we consider the motion of the dense clumps. In Figure 2, we mark the location of the densest point with a red dot in runs with strong particle clumping. By monitoring the location of the densest point with time, we find that it wanders slowly. Another useful way of studying the dynamics of the clumps is by tracking the radial trajectories x_i(t) of a sample of particles. We relocate the particle positions when they cross the radial boundaries of our simulation box so that their trajectories are continuous. By tracing a large number of particles in the saturated state of our runs, we obtain the distribution of x(t +  Δt) −  x(t) for each particle species at time interval Δt. In Figure 7, we show the probability distribution of x(t +  Δt) −  x(t) for a number particle species from our run R10-3D. When Z = 0.01, no particle clumping occurs. The distribution of x(t +  Δt) −  x(t) is close to a Gaussian (or a parabola in logarithmic scale) and the width increases with Δt, consistent with undergoing a random walk. Meanwhile, the center of the distribution drifts inward with time (see Section 5 for more discussion). However, when particle clumps are present, as in the Z = 0.03 case, the shape of the distribution deviates substantially from a Gaussian, especially for particles that make up the clumps (the largest particles, shown in the blue and green curves). For these clump-making particles, the width of particle distribution still increases with Δt, as expected from turbulent diffusion, but a substantial fraction of these particles stay nearly stationary without drifting (near x = 0), making the resulting distribution more and more elongated with time. The leftmost location of the particle distribution moves inward with time, and is set by the radial drift velocity. More interestingly, we see almost evenly separated multiple peaks in the distribution function. In fact, the separation between these peaks equals the radial size of our simulation box. The physical picture becomes clear that the clumps stop some of the particles from drifting radially, and particles are kept in the clump for a few orbits or more before leaving for the next clump. Similar behavior is observed for other runs with particle clumping.

Zoom In Zoom Out Reset image size

We calculate the averaged radial drift velocities for each particle species from all our runs, and the results are shown in Figure 8. The measured mean drift velocities are shown in squares (2D) and circles (3D). We have also plotted the 1σ limits for particle drift velocity based on the rms fluctuations, which are indicated in blue and red vertical bars. In the figure, the velocities are normalized to ηv_K. Clearly, the radial drift velocity monotonically decreases with particle stopping time, and the drift is fastest for marginally coupled particles.

Zoom In Zoom Out Reset image size

The classical result on the radial drift of particles is the NSH equilibrium solution (Nakagawa et al. 1986). It describes the equilibrium state between solids and gas in unstratified (neglecting vertical gravity) Keplerian disks, where gas is partially supported by radial pressure gradient. In the NSH equilibrium, the drift speed is given by

$$\begin{equation} v_x=-\frac{2\tau _s}{(1+\epsilon)^2+\tau _s^2}\eta v_K. \end{equation} \tag{ 19 }$$

We emphasize that the conventional NSH solution is obtained by considering a single species of solids. Equation (19) does not simply generalize to the case with multiple species of particles by replacing to _k for each particle species k. In Appendix A, we provide the generalized formula for multi-species NSH equilibrium, and the solution involves evaluation of an inverse matrix of order 2N_type. It reflects the fact that although different particle species do not interact directly with each other, they are indirectly coupled via their interactions with gas.

In Figure 8, the bold solid lines show the expected radial drift velocities from single-species NSH equilibrium. We see that there are large deviations from the measured mean drift velocities, with two notable features. First, for relatively large particles, the drift velocities are reduced from single-species NSH values. The reduction is strongest for runs with the largest Z. Second, the smallest particles drift outward, rather than inward as expected from the single-species NSH solution.

To calculate the expected radial drift velocity from a multi-species equilibrium, we first use the particle density profiles extracted from Section 3.2 and calculate the drift velocity in each vertical bin. The drift velocity is then weighted by particle density in each bin to yield the mean drift velocity. The results are plotted in dashed and dash-dotted lines (for 2D and 3D runs, respectively) in Figure 8. We see that these curves provide an excellent fit to the measured mean radial drift velocities in all simulations. In fact, the two features mentioned above are natural consequences of the multi-species solution. Due to the sub-Keplerian motion of the gas, particle drag increases gas angular momentum, leading to outward drift of gas. In the presence of both weakly coupled and strongly coupled particles, the strongly coupled particles are tied to the gas and therefore drift outward with the gas. Marginally coupled particles still drift inward, but due to the influence of the smaller particles, these particles feel a weaker headwind (i.e., the gas azimuthal velocity is closer to the Keplerian value), resulting in a smaller drift velocity compared with the single-species solution. With increasing Z, thus higher midplane particle density, the gas becomes more entrained by the solids, leading to stronger reduction of the drift velocity for large particles.

The residuals from the multi-species NSH solution fit to the measured mean drift velocities are largest for particles with largest τ_s, likely due to their participation in SI, and/or clumping. In the non-stratified simulation of Johansen & Youdin (2007), it was shown that in the saturated state of SI, the radial drift velocity is either increased or decreased depending on run parameters. In our simulations, these effects are secondary compared with the multi-species effect. The measured drift velocities from 2D (squares) and 3D (circles) simulations generally agree with each other. The (small) differences can be attributed to the differences in the particle vertical density profiles.

So far we have focused on the mean radial drift velocities. In the saturated state of our simulations, the particle radial drift velocities follow a distribution, due to the SI. We see in Figure 8 that in most of the runs, the fluctuation level is about (0.05–0.15)ηv_K. This fact is closely related to the radial diffusion of particles discussed in the next subsection. Based this observation, we can estimate the particle radial diffusion coefficient to be D_x ∼ (0.1ηv_K)²/Ω ∼ 2.5 × 10⁻⁵ c_sH_g.

The radial diffusion of particles is generally characterized by the radial diffusion coefficient D_x. From our simulations, we can measure D_x for different particle species based on the random walk model of particle diffusion. We calculate the distribution of shift in the particle radial position at various time intervals Δt as in Figure 7, and measure the width (rms) of the distribution σ as a function of Δt. The spreading due to a random walk results in an Gaussian distribution, and σ is related to the diffusion coefficient by

$$\begin{equation} D_x=\frac{1}{2}\frac{d\sigma _x^2}{dt}. \end{equation} \tag{ 20 }$$

For each particle species, we measure σ²_x for different Δt, and fit the slope in of the σ² − Δt curve by linear regression. The results are summarized in Figure 9. The range of the radial diffusion coefficient is consistent to within an order of magnitude of the estimate in the last subsection based on the spread of radial drift velocities. It is also comparable with the vertical diffusion coefficient at disk midplane estimated in Section 3.2 (see Table 2). Below we discuss these results further.

Zoom In Zoom Out Reset image size

First, the above procedure for measuring the diffusion coefficient does not apply to runs that show strong particle clumping. As we see in Figure 7, the distribution of x(t + Δt) − x(t) deviates strongly from a Gaussian due to the influence of the clumps. The measured width of the distribution is about half the distance traveled by the fastest drifting particles (those that are not confined in the clumps), and we observe that σ²_x scales as Δt² rather than Δt from our measurement. Therefore, the measured D_x from R21Z3-2D, R30Z3-2D, and R10Z3 (both 2D and 3D) runs for those clump-making particles (or the largest two particle species in the run) is not valid. In Figure 9, we see the measured D_x for these particles have anomalously large values. Such particles can reside in the disk for much longer than if there were no clumping.

Next, we discuss diffusion of non-clumping particles. In each simulation, the measured D_x generally approaches an asymptotic value for particles with τ_s ≲ 10⁻², but is different between different particle species for particles with τ_s > 10⁻². This can be due to multiple reasons. First, similar to the vertical diffusion of particles, the radial diffusion coefficient also depends on the vertical position in the disk, and the radial diffusion in the disk midplane is expected to be the strongest. Our measured D_x can be considered as a vertically averaged quantity. Therefore, D_x is expected to be larger for particles with larger τ_s, since they stay closer to the midplane. This trend is observed in runs R41 and R21. Second, different particles react differently to the turbulence. In the case of Kolmogorov turbulence, the particle diffusivity scales as (1 + τ²_s)⁻¹ (Youdin & Lithwick 2007). This may be responsible for the decrease of D_x toward τ_s = 1 in R30 and R10 runs with Z = 0.01 and 0.02. Third, different particles participate in the SI in different ways (i.e., actively or passively). The SI may strongly affect the transport properties of the active particles, with the extreme example being the clump-making particles discussed above. Despite the different values of D_x for different particle species, one may take the asymptotic value of D_x as measured from the smallest particles as characteristic of the radial diffusion coefficient in the gas. These asymptotic values correlates with the vertical diffusion coefficient well (see Table 2).

To address the effectiveness of radial diffusion compared with radial drift, we denote the mean radial drift velocity to be v_r = κηv_K, and the diffusion coefficient to be D_x = (βηv_K)²/Ω. After time t, the ratio

$$\begin{equation} \zeta \equiv \frac{v_rt}{\sigma (t)}=\frac{\kappa }{\beta }\sqrt{\Omega t/2}\,, \end{equation} \tag{ 21 }$$

reflects the relative importance between radial drift and turbulent diffusion, where $\sigma (t)=\sqrt{2D_xt}$. Diffusion is important when ζ ≲ 1. From Figure 8, we see that for the largest particles, κ ≳ 0.1. From Figure 9, we have β ≲ 1. Therefore, the effect of radial diffusion of particles becomes negligible compared with radial drift beyond 100 orbital periods. Again, this discussion does not apply to the situation when particle clumping is present, where large particles can be retained in the clumps and some of them may survive the radial drift.

The main purpose of this paper is to study the dynamics of solids and gas in the midplane of PPDs using hybrid simulations. The solids and gas are coupled aerodynamically, characterized by the dimensionless stopping time τ_s = Ωt_stop. We consider a wide size distribution of solids as an approximation to the outcome of grain growth in PPDs, ranging from submillimeter to meter size. The key ingredient of our simulations is the inclusion of feedback from particles to gas. Feedback is important when the local particle-to-gas mass ratio exceeds order unity. Moreover, it is essential for the generation of SI and KHI. In our simulations, we assume no external source of turbulence, as an approximation for the dead zone of PPDs. Turbulence in the disk midplane is generated self-consistently from the SI (driven by the radial pressure gradient in the gas) and/or KHI (driven by vertical shear). Our simulations are local, since very high numerical resolution is essential to resolve the SI and KHI. Self-gravity is ignored, as we focus on the particle–gas dynamics before the formation of planetesimals.

Our simulations are characterized by three sets of dimensionless parameters, namely, the particle size distribution τ_k, solid abundance Z, and a parameter Π characterizing the radial pressure gradient. In this paper, we fix Π = 0.05, as appropriate for a wide range of disk model parameters (see Section 2.2). The dependence of the particle clumping properties on Π is presented in a separate paper (Bai & Stone 2010b). We consider a flat mass distribution in logarithmic bins in τ_s, and vary Z from 0.01 to 0.03 (see Table 1). We conduct both 2D and 3D simulations, where 2D simulations are performed in the radial–vertical plane in order for the SI to be actively generated. We run the simulations for 40–200 orbits and study the properties of the particles and gas in the saturated state. The main results are summarized below.

• 1.  
  SI plays the dominant role in the dynamics of PPD midplane when the largest solids have stopping times τ_s ≳ 10⁻². Particles with τ_s ≳ 10⁻² actively participate in SI, while smaller particles behave passively. KHI is not observed in all our simulations, which suggests that it may be important only when all particles have τ_s ≲ 10⁻².
• 2.  
  The strength of the turbulence generated by the SI and the scale height of the particle layer are self-regulated. There exists some threshold solid abundance, above which increasing Z will result in weaker turbulence, which promotes particle settling, leading to rapid drop of the thickness of the particle layer and strong particle clumping.
• 3.  
  SI can concentrate particles into dense clumps with solid density exceeding the Roche density, which acts as the prelude of planetesimal formation. The particle clumping generally requires the presence of relatively large particles with τ_s ≳ 10⁻². It also sensitively depends on solid abundance, in favor of super-solar metallicity.
• 4.  
  The dense particle clumps are mostly made of the largest particles with size range spanning less than one order of magnitude. These particles are trapped in the clumps for several orbital times before leaving the clumps, providing a way for large particles to survive radial drift.
• 5.  
  The mean radial drift velocity for each particle species agrees well with a multi-species NSH equilibrium solution (see Appendix A). Strongly coupled particles drift outward, and the radial drift velocity for particles with larger τ_s is strongly reduced relative to the conventional single-species NSH value, especially at large Z. This can increase the lifetime of the largest particles by a factor of a few.
• 6.  
  Turbulence generated by the SI leads to radial diffusion of particles, but the diffusion is slow and its effect is negligible compared with radial drift after about 100 orbital periods (for the largest particles). Particle clumping effectively enhances radial diffusion by retaining a fraction of large particles in the clumps.
• 7.  
  Mutual collision velocity between τ_s ≳ 10⁻² particles is dominated by the difference in their radial drift velocities, and agrees well with calculations using the multi-species NSH equilibrium. The collision velocity is strongly reduced toward large disk metallicity relative to predictions from single-species NSH solution. Collision velocity induced by SI turbulence is only secondary.

In this subsection, we combine the results summarized in the previous subsection and discuss various implications for planetesimal formation. In particular, we emphasize the importance of the local enrichment of solid materials in PPDs on planetesimal formation (two feedback loops, see Sections 7.2.2 and 7.2.3). Our logical chain is summarized in Figure 11, and we elaborate various aspects of this diagram in the following.

Zoom In Zoom Out Reset image size

7.2.1. Conditions for Strong Particle Clumping

7.2.2. Enrichment of Local Solid Abundance

7.2.3. Implications for Grain Growth

7.2.4. Influence from Magnetic Activity

Our simulations take a local patch from a simple global disk model in which all physical quantities follow a power-law dependence on disk radii. The global structure of PPDs may be more complicated. In particular, the presence of a dead zone in PPDs changes the steady-state disk surface density profile, and may lead to local pressure maxima at the snow line (Kretke & Lin 2007; Brauer et al. 2008b), and the inner edge of the dead zone (Dzyurkevich et al. 2010). These pressure bumps are able to trap particles very efficiently. Moreover, solids can also be trapped in long-lived vortices (Barge & Sommeria 1995; Klahr & Henning 1997; Klahr & Bodenheimer 2003; Johansen et al. 2004; Lyra et al. 2009), although their existence is in debate. Finally, the structure of PPDs may be very non-steady, and undergo periods of outburst (Zhu et al. 2009). Global models take into account large-scale variations of disk structure and can follow the disk evolution. Currently, the numerical resolution required for resolving the SI is prohibitively high for running 3D global simulations. Nevertheless, one can perform local simulations in different disk environments, and piece them together to form a global picture, as in Tilley et al. (2010).

Our simulations focus on the dynamics in the vicinity of the disk midplane with very limited radial and vertical box sizes and no magnetic field. This is mainly constrained by the fine grid resolution required for this study. In reality, MRI turbulence in the active layers may excite vertical oscillations in the disk midplane (Fleming & Stone 2003; Oishi & Mac Low 2009). Moreover, turbulent diffusion of ions may produce magnetic activity at the disk midplane, making it "undead" (Turner & Sano 2008). Finally, turbulent mixing of particles may become important when the active layer is relatively thick (Turner et al. 2010). Including these effects into numerical simulations of planetesimal formation requires enlarging our box size in all dimensions by a factor of at least 10. Moreover, 3D rather than 2D simulations are necessary to maintain sustained MRI turbulence. Such simulations are computationally expensive, however recently the static mesh refinement (SMR) algorithm in the Athena MHD code has been parallelized. The cost of hybrid (particles and gas) simulations of layered disks can be substantially reduced by using fine resolution at the disk midplane (to capture the SI) and in the active layers (to capture the MRI), while using coarse resolution everywhere else. With SMR, it becomes feasible to study the effect of (non-ideal) MRI turbulence in the active layer on particle dynamics and planetesimal formation in the disk midplane. This is planned as future work.