THE EFFECT OF THE RADIAL PRESSURE GRADIENT IN PROTOPLANETARY DISKS ON PLANETESIMAL FORMATION

Planetesimals are super-kilometer sized bodies that are the building blocks of planets (Safronov 1969; Chiang & Youdin 2010), yet their formation has long been a mystery. Millimeter and centimeter sized particles are routinely observed in protoplanetary disks (PPDs; Wilner et al. 2005; Rodmann et al. 2006; Natta et al. 2007; Ricci et al. 2010; Lommen et al. 2010), but it seems difficult for particles with larger size to form by coagulation (Blum & Wurm 2008; Güttler et al. 2010). Solids close to meter size further suffer from rapid radial drift due to the negative radial pressure gradient (RPG) in PPDs (Weidenschilling 1977). The gravitational instability (GI) scenario of planetesimal formation (Goldreich & Ward 1973) also has difficulties, because even without an external source of turbulence, the Kelvin–Helmholtz instability (KHI) generated from the dusty midplane layer prevents the onset of GI unless the local height-integrated solid to gas mass ratio (Z, hereafter referred to as solid abundance) is about an order of magnitude above solar metallicity (Weidenschilling 1980; Sekiya 1998; Youdin & Shu 2002).

It was found recently that the drag interaction between solids and gas leads to a powerful instability (Goodman & Pindor 2000). This "streaming instability" (SI; Youdin & Goodman 2005) results in spontaneous particle clumping from the equilibrium state between solids and gas (Nakagawa et al. 1986). Numerical simulations demonstrate that the nonlinear saturation of the SI concentrates particles into dense clumps (Johansen & Youdin 2007), promoting planetesimal formation by collective particle self-gravity, bypassing the meter-size barrier. Indeed, Johansen et al. (2007, 2009) showed that planetesimals with sizes of a few hundred kilometers form rapidly by the SI from centimeter to decimeter sized pebbles and rocks, consistent with constraints from the asteroid belt (Morbidelli et al. 2009).

This Letter complements our previous work on the dynamics of particle and gas in the midplane of PPDs (Bai & Stone 2010b, hereafter BS10). We consider a wide size distribution of particles ranging from submillimeter up to meter sizes as an approximation for the outcome of dust coagulation in PPDs (Birnstiel et al. 2010; Zsom et al. 2010). External sources of turbulence such as magnetorotational instability are ignored, as appropriate for the dead zones in PPDs (Gammie 1996; Stone et al. 2000; Bai & Goodman 2009), while SI is the main source of turbulence due to particle settling. Particle self-gravity is ignored in our simulations, as we focus on the precursor of planetesimal formation: particle clumping.

The SI is powered by the RPG in the gaseous disk. The RPG reduces the gas orbital velocity (in the absence of solids) by a fraction η of the Keplerian velocity v_K = Ωr and is characterized by

$$\begin{equation} \Pi \equiv \eta v_K/c_s=\eta r/H_g\, \end{equation} \tag{ 1 }$$

where c_s is the isothermal sound speed and H_g = c_s/Ω is the scale height of the gaseous disk. In BS10, we assumed Π = 0.05 throughout, while in reality Π depends on the parameters of, and location in, the disk. For the widely used minimum mass solar nebular model (MMSN; Hayashi 1981), we have Π ≈ 0.054r^1/4_AU (see Section 2.2 and in particular Equation (7) of BS10 for a full discussion). The strength of the SI turbulence scales with the RPG, which in turn affects the particle–gas dynamics in the disk midplane. In particular, we show in this Letter that particle clumping strongly depends on the RPG, which has important implications for planetesimal formation.

We perform two-dimensional (2D) hybrid simulations of gas and solids using the Athena code (Stone et al. 2008; Bai & Stone 2010a), where gas is treated as a hydrodynamical fluid (without magnetic field) on an Eulerian grid, and the solids are treated as superparticles, each representing a swarm of real particles. We model a local patch of the PPDs using the shearing sheet approximation. The dynamical equations are written in a reference frame corotating at the Keplerian frequency Ω at fiducial radius r. We assume axisymmetry, and the simulations are performed in the radial–vertical (x–z) plane, where Ω is along the z-direction and x points radially outward. The particles are coupled to the gas via aerodynamic drag, characterized by the stopping time t_stop, with momentum feedback included. The gas is assumed to be isothermal, with midplane density ρ_g. Vertical gravity g_z = −Ω²z is included for both particles and the gas.

The dynamical equations and simulation setup are identical to those in BS10. Specifically, we consider a particle size distribution that is discretized into a number of particle size bins, with each bin covering half a dex in τ_s ≡ Ωt_stop, bounded by the minimum and maximum stopping time τ_min and τ_max. In what follows, we label our simulations with names of the form Rmn, where m, n are integers obtained by τ_min = 10^−m and τ_max = 10⁻ⁿ, thus run Rmn uses 2(m − n) + 1 particle size bins (or particle species). We use a variety of grid resolution and box sizes to capture the fastest growing modes of the SI, and we use 10⁵ particles per species in all our simulations. We assume uniform mass distribution across all particle size bins, with total solid abundance Z. As in BS10, we consider four groups of runs, R41, R21, R30, and R10.

For each particle size distribution, we perform simulations with three different values of the RPG parameter: Π = 0.025, 0.05, and 0.10. According to Equations (4) and (7) in BS10, the dependence of Π on disk parameters such as temperature and mass is relatively weak, thus our range of Π covers a large parameter space of disk models. For each set of runs, we perform a series of simulations with different Z values, starting at Z = 0.01 and increasing the value by 0.01 for each new run, until strong particle clumping occurs or Z = 0.07. Our simulation run parameters are summarized in Table 1. They are identical to the 2D run parameters in BS10, but use different values for Π and a larger range in Z. Since the natural length scale of the SI is ηr, we use smaller (bigger) simulation box sizes for smaller (larger) Π values so that ηr is resolved by an equal number of grid cells in each series of runs. All simulations are run for at least 900Ω⁻¹.

Three-dimensional (3D) simulation with the inclusion of the azimuthal dimension is necessary to capture the KHI (Chiang 2008; Barranco 2009), which is mainly caused by the vertical shear in the gas azimuthal velocity. Our simulations are 2D rather than 3D for two reasons. First, for our adopted particle size distribution, the turbulence generated from the SI stops particles from settling before the onset KHI (BS10), at least for the Π = 0.05 case. Second, in order to properly resolve the SI, relatively high resolution is required (Bai & Stone 2010a), and particle clumping does depend on resolution.¹ Simulations in 3D with such high resolution are too costly due to the small box size (thus time step), especially for Π = 0.025.

The simulations saturate in about 50–100 orbits. Particle settling triggers the SI, and particles with different stopping times are maintained at different heights determined by the balance between settling and turbulent diffusion. Most interestingly, SI efficiently concentrates particles into dense clumps when the solid abundance is sufficiently high (Johansen et al. 2009; BS10). In Figure 1, we plot the evolution of maximum particle density ρ_p,max for the R30 runs with different values of Π and Z. There is a clear dichotomy in the evolution: either the maximum particle density stays at a relatively small value (ρ_p,max ≲ 50ρ_g) or it reaches as high as 10³ρ_g, indicative of strong clumping. A particle clump becomes gravitationally bound when its density exceeds the Roche density (see Equation (18) of BS10), which also is of the order 10³ρ_g for typical PPDs. Therefore, particle clumping is a prelude to planetesimal formation via gravitational collapse. For each value of Π, the transition from non-clumping to clumping is sharp as Z gradually increases, and we can define a critical solid abundance Z_crit, where strong particle clumping occurs for Z>Z_crit. Comparing different panels indicates that Z_crit monotonically increases with Π.

Zoom In Zoom Out Reset image size

To better demonstrate the effect of the disk RPG on planetesimal formation, we plot the maximum particle density as a function of solid abundance Z for each group of simulations with all values of Π in Figure 2. The maximum density is taken from the largest value of ρ_p,max over the last 20 orbits of the simulations, when all the runs are fully saturated. We take ρ_p = 500ρ_g (which is the same order as the Roche density) as a rough indicator of planetesimal formation, as shown by the dotted lines. For each value of Π, one can determine Z_crit at the intersection of the dashed and corresponding solid lines, which is given in the last column of Table 1. We see that for all four groups of runs, Z_crit monotonically increases with Π. Moreover, Z_crit depends on Π more sensitively when the particle size is on average smaller. When Π = 0.025, Z_crit is about 0.015 for all our four groups of runs. For Π = 0.1, run R41 does not show particle clumping below Z = 0.07, for runs R21 and R30 Z_crit is about 0.06, while for runs R10 Z_crit drops to below 0.05.

Zoom In Zoom Out Reset image size

Particle clumping is a highly nonlinear effect due to the SI, and is more likely to develop when the average disk midplane solid to gas mass ratio is large. If one assumes D to be the midplane vertical diffusion coefficient of the SI turbulence, particles with stopping time τ_s would settle to a layer with thickness of the order $H_p\approx \sqrt{D/\Omega \tau _s}$. The diffusion coefficient D depends on the particle size distribution, ≈ ZH_g/H_p and Π. Without vertical gravity, the only length scale in the problem is ηr = ΠH_g, thus D ∝ Π². If one assumes D ∝ ^α, we obtain by generalizing the toy model in BS10 (see their Equations (16) and (17))²

$$\begin{equation} H_p\propto Z^{\alpha /(\alpha +2)}\Pi ^{2/(\alpha +2)}. \end{equation} \tag{ 2 }$$

Since one expects α < 0 for relatively large ≳ 1, we see that H_p not only depends sensitively on Z, as shown in BS10, it depends even more sensitively on Π. This critical dependence makes the average particle to gas mass ratio at midplane quickly increase as Π decreases, promoting strong particle clumping, which explains the trend in Figure 2. The more sensitive dependence of Z_crit on Π for smaller particles can be interpreted as a result of the power-law index α tending to be more negative for smaller τ_s.

Our 2D simulations do not contain the azimuthal dimension, which suppresses KHI. Nevertheless, we can check whether KHI would occur before particles settle into a sufficiently thin layer to trigger strong SI by plotting the vertical profile of the Richardson number for simulations at critical solid abundance for all Π values, see Figure 3. The Richardson number is calculated by the method described in Section 3.1 of BS10. Although the Richardson number alone does not determine the stability when Coriolis force (Gómez & Ostriker 2005) and radial shear (Barranco 2009) are present, one may still take Ri_crit = 0.1 as an approximation for the critical value for instability (Chiang 2008; Lee et al. 2010a). We see in Figure 3 that for most of these runs, the Richardson number across the disk is above the critical value, meaning that the system is KH stable at Z_crit. Moreover, as Π increases, the Richardson number at Z_crit decreases. In particular, for R30 and R10 runs with Π = 0.1, the Richardson number at the disk midplane drops below 0.1, which may be subject to KHI. Therefore, the condition for strong clumping for these runs may be more stringent than shown in Figure 2.

Zoom In Zoom Out Reset image size

In this Letter, we have demonstrated that a small radial pressure gradient (RPG) strongly favors particle clumping and planetesimal formation by the streaming instability (SI) via three effects. First, for a fixed particle size distribution, the critical solid abundance Z_crit above which strong clumping occurs is reduced. Second, at fixed solid abundance, strong clumping occurs for smaller particles. Third, the Kelvin–Helmholtz Instability is less likely to be triggered (and therefore will not suppress particle clumping) at Z = Z_crit. Moreover, a smaller RPG also implies smaller radial drift and collision velocities, promoting grain growth, strengthening our second point above. One caveat from our 2D simulations is that the condition for particle clumping is slightly more stringent in 3D than in 2D (BS10). Therefore, the values of Z_crit obtained in this Letter may be considered as a lower bound. Nevertheless, the conclusions in this Letter are robust. In fact, Johansen et al. (2007, see their supplemental material) also found that planetesimal formation is eased with smaller pressure gradient using 3D simulations.

RPG provides the free energy for the SI, therefore, a small RPG weakens the SI turbulence and promotes particle settling to higher densities, which partially explains our results. The nonlinear dependence of particle clumping on local particle density due the SI also plays a crucial role. Similarly, a small RPG also favors the gravitational instability scenario of planetesimal formation (relevant when all particles are strongly coupled to the gas), because there is less free energy available in the vertical shear to prevent dust settling. However, the threshold abundance Z_GI for gravitational instability to operate is generally larger than Z_crit obtained here. For example, under the condition of the minimum mass solar nebular, Z_GI is around 0.1 for Π = 0.05 (Youdin & Shu 2002). More recent simulations by Lee et al. (2010b) found that Z_GI ≈ 0.06 for Π = 0.025.

Combined with the results in BS10, we have shown that strong particle clumping is favored when there is: (1) a small RPG, (2) large solid abundance, (3) large solid size, and (4) a dead zone. These results indicate that planetesimals preferentially form in specific locations in PPDs. The global structure and evolution of PPDs is of crucial importance: it determines the RPG profile, the radial transport of particles (Youdin & Shu 2002; Haghighipour & Boss 2003a, 2003b), and grain coagulation that determines the particle size distribution. Recent models on the structure and evolution of PPDs (e.g., Jin & Sui 2010; Zhu et al. 2010) generally show much more complicated structures than the simple minimum mass solar nebular or α-disk models due to non-steady state accretion as well as the presence of dead zones. In addition, RPG may reach zero in local pressure bumps at the snow line (Kretke & Lin 2007) or the inner edge of the dead zone (Dzyurkevich et al. 2010) which can also be the preferred site for planetesimal formation. These results all-together provide useful constraints on the initial conditions for the formation of planet embryos and planets (Kokubo & Ida 2000, 2002).

We are grateful to Anders Johansen and D.N.C. Lin for helpful discussions. We thank Andrew Youdin and an anonymous referee for useful comments. This work is supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program Grant NNX09AQ90H awarded to X.N.B.