© ESO, 2015

1. Introduction

Planetesimals are super-kilometer-sized bodies that are the seeds of terrestrial planets, as well as the cores of ice and gas giants (Safronov 1972; Chiang & Youdin 2010; Johansen et al. 2014). Hence, the production of planetesimals is an important step in planet formation. In the solar system, large asteroids and Kuiper belt objects are the left-over planetesimals that did not become incorporated into a planet. Similar planetesimal belts are known to exist around other stars (Wyatt 2008).

The size distribution of the solar system planetesimals (i.e., large asteroids and Kuiper belt objects) reveals a “knee” at around 100 km in size, with a relative scarcity of bodies smaller than ~100 km. Some authors have argued that this knee cannot be reproduced by bottom-up coagulation scenarios, unless the minimum initial size of a planetesimal was also ~100 km. Therefore, it appears that the size of solids in the protoplanetary disk “jumped” from the submeter scale to ~100 km without passing through the intermediate sizes, and that bodies smaller than 100 km formed later by collisional grinding (Morbidelli et al. 2009). A similar bump is inferred for the trans-Neptunian population (Fraser & Kavelaars 2009; Sheppard & Trujillo 2010; Shankman et al. 2013; Fraser et al. 2014). That said, Weidenschilling (2010) suggests an alternate scenario in which the “knee” is produced in a bottom-up process from an initial population of 0.1 km-sized bodies. In this scenario, the bump at ~100 km is produced by the transition from dispersion-dominated runaway growth to a regime dominated by Keplerian shear.

There are also important theoretical difficulties in building planetesimals in a bottom-up process. The first hurdle is that small particles in the mm-cm size range do not easily stick together to form larger objects. Indeed, the largest grains observed in a protoplanetary disk are mm-cm in size (Testi et al. 2014; van der Marel et al. 2013). To the extent that meter-sized objects do form, they quickly spiral into the central star due to friction with gas orbiting the star at sub-Keplerian speeds (Weidenschilling 1977). Axisymmetric pressure bumps may stop particle drift (Johansen et al. 2009a), but growth rates remain low (Johansen et al. 2008). For example, bodies at 3 au can grow to a maximum of 100 m in 1 Myr (Windmark et al. 2012).

These constraints suggest a different picture of how planetesimals form. It begins with the growth of macroscopic particles by coagulation, followed by the accumulation of these particles by a hydrodynamical process like the streaming instability (see below) that increases the local concentration of solids, $\begin{array}{ccc}\mathit{Z}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathrm{\Sigma }}_{\mathrm{solid}}}{{\mathrm{\Sigma }}_{\mathrm{total}}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathrm{\Sigma }}_{\mathrm{solid}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}\mathrm{+}{\mathrm{\Sigma }}_{\mathrm{solid}}}\mathit{,}& & \end{array}$(1)where Σ_solid and Σ_gas are the surface densities of the solid and gas components of the disk respectively. The solids also sediment toward the midplane of the disk, so that the volume density of solids vs. gas at the midplane can be significantly higher than Z. If the density of solids reaches the Roche density, $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{R}}\mathrm{=}\frac{\mathrm{9}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }\mathrm{G}}\mathit{,}& & \end{array}$(2)where Ω is the Keplerian frequency, then the particle self-gravity overwhelms Keplerian shear, leading to gravitational collapse (Goldreich & Ward 1973). This allows a collection of small objects to coalesce directly into fully formed planetesimals, thereby explaining the large birth sizes.

The main obstacle to achieving high densities is disk turbulence, which is thought to be caused by the magneto-rotational instability (MRI, Balbus & Hawley 1991). But even inside an MRI-free “dead zone” (Gammie 1996), the Kelvin-Helmholtz instability can be very effective at producing turbulence that blocks particle sedimentation to the midplane before they can reach the Roche density (Weidenschilling 1980; Youdin & Shu 2002; Lee et al. 2010). Particle densities high enough for a local gravitational collapse can be obtained by concentration in long-lived gaseous vortices (Barge & Sommeria 1995; Meheut et al. 2012) or pressure bumps caused by the MRI (Johansen et al. 2009a, 2011). These regions can quickly accumulate meter- and decimeter-sized bodies, and for meter-sized boulders, the concentrations can significantly reduce the radial drift caused by friction with the surrounding gas (Johansen et al. 2006, 2007).

The streaming instability is another powerful mechanism for concentrating particles in localized dense clumps and preventing radial drift (Youdin & Goodman 2005; Johansen et al. 2007; Johansen & Youdin 2007; Bai & Stone 2010a). It is driven by the relative drift between the solid and gas components of the disk. The gas component of the disk experiences a radial pressure support that counters the gravitational force of the star. This leads to a speed difference Δv between the gas and solid components of the disk. A useful way to understand the streaming instability is that solids experience a head wind from the gas, and the gas is, in turn, pushed forward by solids. A small overdensity of solids will have a stronger backreaction on the gas, hence a lower radial drift than neighboring regions. The reduced radial drift leads to the further accumulation of solids as neighboring particles drift into the overdensity. In this way, the streaming instability bears some similarity to a traffic jam.

The behavior of solid particles embedded in a gaseous accretion disk is primarily determined by the Stokes number, τ_f = t_f Ω_k, which measures the particle friction time t_f in terms of the Keplerian frequency Ω_k. Some authors have raised the concern that the streaming instability has mainly been studied for relatively large particles, with τ_f ~ 0.1, corresponding to dm-sized particles when applied to the asteroid belt (Shi & Chiang 2013). Particles of this size are inconsistent with coagulation experiments, which show that particles cannot grow beyond mm sizes (Zsom et al. 2010; Güttler et al. 2010). Therefore, the main goal of this paper is to study the streaming instability for particle sizes down to millimeter.

Chondrules make up most of the mass in primitive meteorites, with a typical size of R ~ 0.3mm, and large chondrules reaching R ~ 1mm (e.g., Jacquet 2014). At this size range, chondrules couple with the small-scale turbulent eddies in the disk, and they concentrate in the high-pressure regions between the eddies (Cuzzi et al. 2001). This alone does not lead to particle growth, because a typical collision between two chondrules in the disk will result in bouncing (Güttler et al. 2010). However, a collision between a small μm-sized dust particle and a chondrule does result in sticking. Ormel et al. (2008) have shown that chondrules acquire a porous dust layer around them, which absorbs some of the kinetic energy of collisions and allows the chondrules to stick more easily. Under typical disk conditions, most of the chondrule and dust mass is in the form of millimeter-sized aggregates. For low turbulence regions, the chondrule aggregates have typical sizes around R ~ 4mm, and maximum sizes as high as R ~ 10cm (Ormel et al. 2008). Our work establishes the connection between the streaming instability and chondrules and chondrule aggregates. While we test the limits of the streaming instability at both the small-particle and the large-particle limits, the bulk of our analysis focuses on the size range of chondrules and chondrule aggregates.

In addition to particle size, the streaming instability also depends on the magnitude of the radial pressure gradient (Bai & Stone 2010c), as well as the solid concentration (Johansen et al. 2009b; Bai & Stone 2010a). In particular, Bai & Stone (2010a) have shown that the streaming instability is effective in producing particle clumps for particle size τ_f ≳ 0.1, and solid concentration Z ≳ 0.02.

A number of global processes increase the local solid concentration, including radial drift of particles from the outer disk into the inner disk (Youdin & Shu 2002), large scale pressure bumps or vortices (Johansen et al. 2009a), and lastly, there is a region near the ice line where ice-dust aggregates (“dirty snowballs”) break up on a timescale comparable to the infall timescale, leading to a significant accumulation of solids by as much as a factor of 6.7 (Sirono 2011). The second way to increase Z is to remove gas from the disk, either through late-stage photoevaporation (e.g., Alexander et al. 2006a,b) or through disk winds (Bai & Stone 2013; Suzuki & Inutsuka 2009). Gas removal has the double effect of increasing τ_f and Z simultaneously. Hence, gas removal may be a powerful way to trigger the streaming instability in the late stages of the disk.

This paper is structured as follows. In Sect. 2 we present our numerical model, including the simulation setup and initial conditions. In Sect. 3 we present our core result – a map of the region in the Z vs. τ_f phase space that is consistent with particle clumping (Fig. 8). In Sect. 4 we review the implications and offer new constraints on planetesimal formation models (Fig. 9). In Sect. 5 we measure the collision speed for τ_f = 0.003 particles and find that particle growth occurs faster than the disk accretion timescale. In Sect. 6 we present our conclusions.

2. Numerical model

We use the Pencil Code (Youdin & Johansen 2007) to model a small vertical slice of the protoplanetary disk in which the solar system planets formed using a stratified model that includes vertical gravity. The canonical model for this disk is the minimum mass solar nebula (MMSN, Hayashi 1981). In it, the surface density of the gas component of the disk follows the power law $\mathrm{\Sigma }\mathrm{=}\mathrm{1700}\mathrm{g}{\mathrm{cm}}^{-2}\hspace{0.17em}{\left(\frac{\mathit{r}}{\mathrm{au}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{·}$(3)In addition to the gas, the disk also contains solid particles that are initially μm in size and represent about 1% of the disk mass. Over time, particle sizes can grow by coagulation, and the mass ratio between gas and solids may change if particles migrate through the disk, or if the gas becomes depleted.

2.1. Aerodynamics of solid particles

The gas component of the disk produces a drag force on the particles whenever the particles have a non-zero relative speed v_rel with respect to the gas. The form of the gas drag law depends on the size of the particle relative to the mean free path λ of gas particles. Small particles, with size R< 9λ/ 4 experience Epstein drag, while larger particles experience Stokes drag. We will show in Sect. 2.4 that the the smaller particles in our simulations lie in the Epstein regime. Particles in the Epstein regime experience the drag force ${\mathit{F}}_{\mathrm{drag}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}{\mathit{R}}^{\mathrm{2}}{\mathit{v}}_{\mathrm{th}}\hspace{0.17em}{\mathit{v}}_{\mathrm{rel}}$(4)where ρ is the gas density, R is the particle size, and v_th is the mean thermal speed of gas molecules, $\begin{array}{ccc}{\mathit{v}}_{\mathrm{th}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\sqrt{\frac{\mathrm{8}\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathit{T}}{\mathit{\pi }\hspace{0.17em}\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}}}\mathrm{·}& & \end{array}$Here k is the Boltzmann constant, T is the local temperature, μ is the mean molecular weight of the gas, and m_p is the proton mass (Armitage 2010). It is convenient to write v_th in terms of the isothermal sound speed c_s, $\begin{array}{ccc}{\mathit{c}}_{\mathrm{s}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\sqrt{\frac{\mathit{k}\hspace{0.17em}\mathit{T}}{\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}}}\mathrm{\Rightarrow }{\mathit{v}}_{\mathrm{th}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}\sqrt{\frac{\mathrm{8}}{\mathit{\pi }}}\mathrm{·}& & \end{array}$The friction time $\begin{array}{ccc}{\mathit{t}}_{\mathrm{f}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathit{m}\hspace{0.17em}{\mathit{v}}_{\mathrm{rel}}}{{\mathit{F}}_{\mathrm{drag}}}& & \end{array}$measures how long it takes for a particle with mass m to have its relative speed changed by order unity as a result of gas drag. For a particle in the Epstein regime (Eq. (4)) with material density ρ_•, we can write the friction time as ${\mathit{t}}_{\mathrm{f}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{•}}\hspace{0.17em}\mathit{R}}{\mathit{\rho }\hspace{0.17em}{\mathit{v}}_{\mathrm{th}}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{•}}\hspace{0.17em}\mathit{R}}{\mathit{\rho }\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}}\sqrt{\frac{\mathit{\pi }}{\mathrm{8}}}\mathit{,}$(5)where ρ is the gas density. We commonly express t_f in terms of the Keplerian frequency Ω_k to obtain the Stokes number $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{f}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{\Omega }}_{\mathrm{k}}\hspace{0.17em}{\mathit{t}}_{\mathrm{f}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{\Omega }}_{\mathrm{k}}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{•}}\hspace{0.17em}\mathit{R}}{\mathit{\rho }\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}}\sqrt{\frac{\mathit{\pi }}{\mathrm{8}}}\mathrm{·}& & \end{array}$The Stokes number is a scale-free measure of the stopping time. Since the disk scale height is H = c_s/ Ω_k, the Stokes number can be written as ${\mathit{\tau }}_{\mathrm{f}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{•}}\hspace{0.17em}\mathit{R}}{\mathit{\rho }\hspace{0.17em}\mathit{H}}\sqrt{\frac{\mathit{\pi }}{\mathrm{8}}}\mathrm{·}$(6)

2.2. Size of particles in the Epstein regime

We are interested in the formation of planetesimals at the midplane of the disk, where the gas density is $\mathit{\rho }\mathrm{=}\frac{\mathrm{\Sigma }}{\mathit{H}\hspace{0.17em}\sqrt{\mathrm{2}\mathit{\pi }}}\mathrm{·}$(7)With Eq. (6), this gives the particle Stokes number in terms of the local surface density ${\mathit{\tau }}_{\mathrm{f}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{•}}\hspace{0.17em}\mathit{R}}{\mathrm{\Sigma }}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{·}$(8)Combined with Eq. (3), we obtain the relation between R and τ_f for the MMSN $\mathit{R}\mathrm{~}\mathrm{78}\mathrm{cm}{\mathit{\tau }}_{\mathrm{f}}{\left(\frac{\mathit{r}}{\mathrm{2.5}\hspace{0.17em}\mathrm{au}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}{\left(\frac{{\mathit{\rho }}_{\mathrm{•}}}{\mathrm{3.5}\mathrm{g}{\mathrm{cm}}^{-3}}\right)}^{-1}\mathrm{·}$(9)Note that ρ_• = 3.5gcm^-3 is a typical density for a mm-sized chondrule (Hughes 1980), and r = 2.5au lies in the inner part of the present day asteroid belt. Note also that this equation is only valid within the Epstein regime. In Sect. 2.4 we show that this corresponds to τ_f ≤ 1. Our two largest particle sizes, τ_f = 3 and 10, lie in the Stokes regime. These particles are a few meters in size.

2.3. Pressure support

The gas component in a protoplanetary disk experiences a radial pressure support that partially counters the gravitational force from the central star. As a result, the gas in the disk orbits at sub-Keplerian speed. The difference between Keplerian speed v_k and and gas speed u_φ can be written as $\mathrm{\Delta }\mathit{v}\mathrm{\equiv }{\mathit{v}}_{\mathrm{k}}\mathrm{-}{\mathit{u}}_{\mathit{\phi }}\mathrm{=}\mathit{\eta }{\mathit{v}}_{\mathrm{k}}\mathit{,}$(10)where the pressure gradient parameter η is given by $\mathit{\eta }\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathit{v}}_{\mathrm{k}}}\right)}^{\mathrm{2}}\frac{\mathit{\partial }\ln \mathit{P}}{\mathit{\partial }\ln \mathit{r}}\mathit{,}$(11)where P is the gas pressure (Nakagawa et al. 1986). In the absence of gas, the solid component of the disk would orbit at the Keplerian speed v_k. In the presence of gas, the solids experience a “headwind” from the gas component, with a wind speed of Δv. This headwind is a source of drag that can lead to radial drift. For scale-free numerical simulations, it is helpful to normalize the headwind speed by the sound speed and write $\mathrm{\Delta }\mathrm{\equiv }\frac{\mathrm{\Delta }\mathit{v}}{{\mathit{c}}_{\mathrm{s}}}\mathrm{=}\mathit{\eta }\frac{{\mathit{v}}_{\mathrm{k}}}{{\mathit{c}}_{\mathrm{s}}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathit{v}}_{\mathrm{k}}}\right)\frac{\mathit{\partial }\ln \mathit{P}}{\mathit{\partial }\ln \mathit{r}}\mathrm{·}$(12)This Δ is a free parameter that measures the strength of the pressure support and the headwind. The pressure is $\mathit{P}\mathrm{=}\mathit{\rho }\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$. Equation (7)gives the gas density at the midplane, which combined with Ω_k = c_s/H allows us to write $\begin{array}{ccc}\mathit{P}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathrm{\Sigma }}{\mathit{H}\hspace{0.17em}\sqrt{\mathrm{2}\mathit{\pi }}}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathrm{\Sigma }\hspace{0.17em}{\mathrm{\Omega }}_{\mathrm{k}}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}}{\sqrt{\mathrm{2}\mathit{\pi }}}\hspace{0.17em}\mathrm{\propto }\hspace{0.17em}{\mathit{r}}^{-3}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}\mathit{.}& & \end{array}$The isothermal sound speed c_s is determined primarily by the local temperature ${\mathit{c}}_{\mathrm{s}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\sqrt{\frac{\mathit{k}\hspace{0.17em}\mathit{T}}{\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{60}\mathrm{m}{\mathrm{s}}^{-1}{\left(\frac{\mathit{\mu }}{\mathrm{2.3}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{\mathit{T}}{\mathrm{K}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(13)where k is the Boltzmann constant, T is the temperature, m_p is the proton mass, and μ is the mean molecular weight. In the MMSN the disk is assumed to be optically thin, so that its temperature profile is T = 280K (r/ au)^{− 1 / 2}. However, Chiang & Youdin (2010) proposed a colder disk, with T = 120K (r/ au)^{− 3 / 7}. For the moment, then, we will simply write the temperature profile as $\begin{array}{ccc}\mathit{T}\mathrm{=}{\mathit{T}}_{\mathrm{0}}\hspace{0.17em}{\left(\frac{\mathit{r}}{\mathrm{au}}\right)}^{\mathrm{-}\mathit{\beta }}\mathrm{\Rightarrow }{\mathit{c}}_{\mathrm{s}}\mathrm{=}\mathrm{60}\mathrm{m}{\mathrm{s}}^{-1}\sqrt{{\mathit{T}}_{\mathrm{0}}}{\left(\frac{\mathit{r}}{\mathrm{au}}\right)}^{\mathrm{-}\mathit{\beta }\mathit{/}\mathrm{2}}\mathrm{·}& & \end{array}$Therefore, the pressure P ∝ r^-3c_s follows a power law with exponent $\begin{array}{ccc}\frac{\mathit{\partial }\ln \mathit{P}}{\mathit{\partial }\ln \mathit{r}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\left(\mathrm{3}\mathrm{+}\frac{\mathit{\beta }}{\mathrm{2}}\right)\mathrm{·}& & \end{array}$Equation (12)can now be written as $\mathrm{\Delta }\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{10}}^{-3}\sqrt{{\mathit{T}}_{\mathrm{0}}}\left(\mathrm{3}\mathrm{+}\frac{\mathit{\beta }}{\mathrm{2}}\right){\left(\frac{\mathit{r}}{\mathrm{au}}\right)}^{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}\mathit{/}\mathrm{2}}\mathrm{·}$(14)Using the parameters from the MMSN and the Chiang-Youdin nebula (CY) we get $\begin{array}{ccc}& & {\mathrm{\Delta }}_{\mathrm{MMSN}}\mathrm{=}\mathrm{0.068}\hspace{0.17em}{\left(\frac{\mathit{r}}{\mathrm{2.5}\mathrm{au}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}\\ & & {\mathrm{\Delta }}_{\mathrm{CY}}\mathrm{=}\mathrm{0.046}\hspace{0.17em}{\left(\frac{\mathit{r}}{\mathrm{2.5}\mathrm{au}}\right)}^{\mathrm{2}\mathit{/}\mathrm{7}}\mathrm{·}\end{array}$Given the uncertainty in the disk model, we choose Δ = 0.05 as a typical value for r ~ 2.5 au (present day asteroid belt). This is also the value chosen by Bai & Stone (2010a), so that using Δ = 0.05 facilitates comparison with their work. We also ran simulations with Δ = 0.025 to study the effect of a reduced pressure support. This can occur as part of a large-scale pressure bump, which could form around the snow line (Kretke & Lin 2007).

2.4. Epstein vs. Stokes regimes

Particles with R< 9λ/ 4 experience Epstein drag, where the mean free path λ is given by $\begin{array}{ccc}\mathit{\lambda }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{n}\hspace{0.17em}\mathit{\sigma }}\mathit{,}& & \end{array}$where σ = 2 × 10^-15cm² is the collision cross section of an H₂ molecule, and n is the particle number density, which is $\begin{array}{ccc}\mathit{n}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathit{\rho }}{\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}}\mathrm{·}& & \end{array}$Using Eq. (7), along with H = c_s/ Ω_k and Eq. (13), one can compute the mean free path. For the MMSN, the mean free path is λ ~ 43.6 cm at 2.5 AU. That means that particles smaller than R ~ 98cm lie in the Epstein regime. From Eq. (9)we find that all our particles with τ_f ≤ 1 lie in the Epstein regime.

2.5. Simulation setup

In our simulations we use c_s ≡ Ω_k ≡ 1 as a natural scaling unit for the problem. The free parameters for the problem are Δ and τ_f given by Eqs. (8)and (12). We use the Pencil Code¹ to simulate the dynamics of gas and solid particles in a 2-dimensional vertical slice of the protoplanetary disk (Youdin & Johansen 2007). We model a box of dimensions 0.2H × 0.2H in the radial-vertical plane centered on the midplane of the disk. The box has periodic boundary conditions. Our baseline simulations have a 128 × 128 grid resolution with 16 384 “super-particles” representing the solids. Each super-particle represents a large number of physical particles. We run a total of 102 baseline simulations as we vary the key simulation parameters:

• We test seventeen particle sizes spaced equally in log scale fromτ_f = 0.001 to τ_f = 10.

• We test both the canonical pressure support Δ = 0.05, and a reduced pressure support Δ = 0.025.

• We run each (τ_f,Δ) simulation three times, with different initial particle positions. The particle positions are chosen from a random distribution, uniform within the box.

Every simulation begins with a solid fraction of Z = 0.005 (Eq. (1)). First, we allow the solids to sediment to the midplane. In most cases, the sedimentation is complete within 50 orbits. For the τ_f< 0.003 simulations we allowed 250 orbits for the sedimentation phase. Once the particles reach an equilibirum scale height, we begin the numerical experiment: we reduce the gas density ρ exponentially, so that it halves every 50 orbits. We keep τ_f constant, which is equivalent to reducing the particle size R at the same rate as ρ (see Eq. (6)). This procedure is also equivalent to increasing the particle density ρ_p while keeping ρ and R fixed. The objective of the experiment is to determine the value Z needed to produce particle clumps for a given τ_f.

2.6. System response time

It is desirable that the rate of gas removal be slower than the response time of the system. Otherwise our procedure will tend to over-estimate the Z value needed to produce particle clumps. The response time of the system is in the order of the particle crossing time t_cross, which is ${\mathit{t}}_{\mathrm{cross}}\mathrm{=}\frac{\mathrm{0.2}\mathit{H}}{\mathrm{\left|}{\mathit{v}}_{\mathrm{r}}\mathrm{\right|}}\mathit{,}$(17)where 0.2H is the width of the box and v_r is the particle radial velocity, which is given by $\begin{array}{ccc}& & {\mathit{v}}_{\mathrm{r}}\mathrm{=}\frac{\mathrm{-}\mathrm{2}\mathit{\eta }{\mathit{v}}_{\mathrm{k}}}{{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}{\mathit{\tau }}_{\mathrm{f}}^{-1}}\mathrm{=}\frac{\mathrm{-}\mathrm{2}{\mathit{c}}_{\mathrm{s}}\mathrm{\Delta }}{{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}{\mathit{\tau }}_{\mathrm{f}}^{-1}}\mathit{,}\\ & & \mathrm{\Rightarrow }{\mathit{t}}_{\mathrm{cross}}\mathrm{=}\frac{\mathrm{0.1}\mathit{H}}{{\mathit{c}}_{\mathrm{s}}\mathrm{\Delta }}\mathrm{(}{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}{{\mathit{\tau }}_{\mathrm{f}}^{-1}}^{\mathrm{)}}\mathit{.}\end{array}$With this one can show that t_cross< 50 orbits $\mathrm{\iff }{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}{\mathit{\tau }}_{\mathrm{f}}^{-1}\mathit{<}\mathrm{50}\mathit{\pi }$, for Δ = 0.05. This roughly corresponds to the interval 0.0064 <τ_f< 157. This means that for most of our simulations the system response time is less than 50 orbits, but for τ_f< 0.0064 we will overestimate the Z needed for particle clumps. For τ_f< 0.0064, our results should be considered robust upper limits. This is a necessary limitation because these simulations also suffer from small time steps which force long computation times.

3. Results

A reduced pressure support (Δ = 0.025) is associated with slower radial drift, lower gas turbulence, and a correspondingly higher concentration of solids at the midplane. We can quantitatively measure the scale height of particles (H_p) and the degree of turbulence using the root-mean-square of the particle positions (z coordinate) and vertical speed, $\begin{array}{ccc}& & {\mathit{H}}_{\mathrm{p}}\mathrm{=}\sqrt{\mathrm{⟨}{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{n}}\sum _{\mathit{i}}{\mathit{z}}_{\mathit{i}}^{\mathrm{2}}}\mathit{,}\\ & & {\mathit{v}}_{\mathrm{rms}}\mathrm{=}\sqrt{\mathrm{⟨}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{n}}\sum _{\mathit{i}}{\mathit{v}}_{\mathit{z,i}}^{\mathrm{2}}}\mathit{,}\end{array}$where z_i is the z coordinate of the ith particle, and v_z,i is the vertical component of the particle’s velocity. Figure 1 shows H_p and v_rms across the entire range of simulations. For particles smaller than τ_f ~ 0.3, the lower Δ is associated with both a lower scale height and lower vertical motion. Particles larger than τ_f ~ 0.3 are poorly coupled with the gas, and are largely unaffected by Δ. The effect of particle size and concentration is more complex.

We found it helpful to divide the particles into three different size bins, based on their qualitative behavior: (1) particles smaller than τ_f ~ 0.003 form a suspension, where particles cannot decouple from the gas even at high concentration and particle clumps do not occur. (2) Between τ_f ~ 0.003 and τ_f ~ 0.3 is the optimal range for particles to engage in the streaming instability. (3) Particles larger than τ_f ~ 0.3 suffer strong radial drift, and the over-densities that do form are quickly destroyed by gas erosion.

3.1. Suspension regime (τ_f < 0.003)

Our smallest particles have τ_f = 10^-3 and τ_f = 10^-2.75 ≈ 0.002. These particles do not sediment easily because they are strongly coupled to the gas. Therefore, we wait 250 orbits for the solids to sediment before we begin the experiment. Figure 2 shows the evolution of H_p, including the sedimentation phase. The long sedimentation phase ensures that H_p has converged before we start removing gas.

Figure 2 shows that, for τ_f< 0.003, the particle size has very little effect on H_p. Instead, H_p is determined by Δ. For Δ = 0.05, the particle concentration is also important for Z> 0.04. Some authors have proposed that, for well-coupled particles such as these, there is a critical Z beyond which turbulent mixing becomes ineffective, resulting in a high-density “cusp” at the midplane (Sekiya 1998; Youdin & Shu 2002). The critical value of Z depends on the disk conditions, and range from Z ~ 0.02 for a cold disk, and Z ~ 0.10 for a MMSN (Youdin & Shu 2002). To test this idea, we ran 256 × 256 simulations with τ_f = 10^-2.75. The particle density can be written as ${\mathit{\rho }}_{\mathrm{p}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathrm{\Sigma }}_{\mathrm{solid}}}{{\mathit{H}}_{\mathrm{p}}\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(\frac{\mathrm{-}{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{H}}_{\mathrm{p}}^{\mathrm{2}}}\right)\mathit{,}$(22)where z is the vertical coordinate. Figure 3 shows the midplane particle density ρ_p(z = 0) as a function of Z, and a snapshot of ρ_p(z) at Z = 0.08. Although we find an increase in ρ_p(z = 0) with Z, it is not nearly of the magnitude proposed by Youdin & Shu (2002). We suspect that our simulations have an additional source of mixing that interferes with the formation of the Sekiya-Youdin-Shu cusp. For r = 1 au, our 128 × 128 grid cells have a physical size around 5000−8000 km. Both the 128 × 128 and 256 × 256 runs should be able to resolve the cusp seen in Fig. 1 of Youdin & Shu (2002). On the other hand, the discrepancy between the 256 × 256 runs and the 128 × 128 runs indicates that our results have not converged, and higher resolutions may reveal higher densities.

3.2. Streaming regime (0.003 ≤τ_f≤ 0.3)

After τ_f ~ 0.003, solid particles begin to decouple from the gas, leading to a regime where the radial drift is sufficient to make the streaming instability effective. Figures 4 and 5 show spacetime diagrams for a selection of runs across all regimes, from τ_f = 10^-3 to τ_f = 3. For a standard pressure gradient (Δ = 0.05), there are visible, long-lived filaments for 0.003 ≤ τ_f ≤ 0.3. For τ_f = 1 the filaments are short lived. Somewhere between τ_f = 0.3 and τ_f = 1 particle clumps start to become difficult again.

When Δ = 0.05, the smallest particles that form visible filaments are have τ_f = 0.003. The pile-up process is somewhat analogous to a traffic jam: An initial overdensity of solids has a stronger back-reaction on the gas. This reduced Δ, lowers the radial drift, leading to a pile-up as other solids (with faster radial drift) enter the filament and increase the density even further. Somewhere between τ_f = 0.03 and 0.1, clumping starts to become difficult again, as the smaller solid cross section (fewer particles, higher mass per particle) cannot lower Δ as effectively. After around τ_f ≳ 1, traffic jams are no longer easily observable.

3.3. Radial drift regime (τ_f > 0.3)

The regime with τ_f> 0.3 is not a good environment for asteroid formation. Although particles in this regime do clump, they continue to rapidly drift toward the star, and the clumps that form are very transient as they are easily eroded by wind in the disk. Figure 6 shows the radial speed for various particle sizes and concentrations. Higher concentrations of solids are associated with a decrease in radial speed, but even at high Z the radial speed remains high. In this regime, larger particles are more able to resist gas drag and exhibit lower radial drift. However, even at τ_f ~ 10, the radial drift is as high as | v_r | ~ 0.01c_s. At that speed, these particles would cross the disk in around a thousand years.

Als note that the overdensities that do form in this regime are very short lived, as they are easily destroyed by erosion from the gas component. But in any case, it seems unlikely that this regime is at all important for realistic disks. As noted earlier, the “particles” in this regime are large. A τ_f = 10 particle is actually a 5-m boulder. The inefficiency of sticking (Güttler et al. 2010), and the streaming regime, probably present formidable barriers for any solids to reach the radial drift regime.

3.4. Conditions for particle streaming

To locate particle clumps we look for stable peaks in the solid surface density Σ_solid. The linear stability analysis of Youdin & Goodman (2005) shows that the streaming instability grows only in modes which have both a radial and a vertical variation. The non-linear simulations of Johansen & Youdin (2007) also show significant vertical variation in the particle density. However, simulations which include the vertical gravity on the particles are stratified; they have a dense midplane layer and little additional variation in the vertical direction. For this reason, the vertical average captures very well the magnitude of particle overdensities in the stratified simulations that we present here.

Figure 7 illustrates our strategy. We are interested in identifying localized peaks in particle density that are stable over multiple orbits. That is to say, an overdensity needs to have low radial drift and high longevity. For example, for τ_f = 0.001 there are no strong overdensities, and for τ_f = 3, the overdensities are short-lived and suffer strong radial drift. Therefore, our strategy is to integrate the particle density over a 25-orbit interval and use the Kolmogorov-Smirnov (KS) test to distinguish the resulting particle distribution from a uniform distribution. We test several 25-orbit periods, spaced in steps of 10 orbits so that they overlap. The result of the KS test is a p-value, which measures the probability that two observed data sets arose from the same underlying distribution. The method is discussed in Appendix A. We divide the results of the KS test into four categories:

• Clumping is unlikely:

  p ≥ 0.45

• Clumping is somewhat likely:

  0.25 ≤ p< 0.45

• Clumping is likely:

  0.10 ≤ p< 0.25

• Clumping is very likely:

  p< 0.10.

Figure 8 summarizes the result of the KS test for all the simulations. The symbols in the figure mark the points where the test was performed. The figure is color-coded: for each point, the surrounding square has two colors (possibly the same color twice) that mark the extreme values of p. In this way, the figure indicates the range of values observed across the simulations. The figure shows that, within the streaming regime, there is a distinct region in the Z − τ_f phase space where particle clumps are consistently likely to form. Of special interest is the fact that particles as small as τ_f = 0.003 can form distinct particle clumps, at least some of the time, for solid concentrations around Z ~ 0.065. This value is high, but as we noted earlier, there is evidence that the ice line is associated with a particle concentration at least that high (Sirono 2011), even before we consider turbulent eddies, disk winds or photoevaporation.

Appendix B has convergence tests. Figure B.1 shows the spacetime diagram of a 3D simulation with τ_f = 0.03; Fig. B.5 shows the spacetime diagrams of higher resolution simulations with τ_f = 0.003 and τ_f = 0.010, and grid size 256 × 256. In all cases, the results are consistent with the spacetime diagrams in Figs. 4 and 5.

4. Constraints on planetesimal formation models

We consider the formation of dense particle clumps, as the precursors of planetesimals. Ormel et al. (2008) have found that μm-sized dust sticks to chondrules to form a porous rim, which absorbs some of the kinetic energy from collisions. This allows the chondrules to form small aggregates, whose size depends primarily on the speed of the turbulent eddies in the disk, ${\mathit{v}}_{\mathrm{eddy}}\mathrm{~}{\mathit{c}}_{\mathrm{s}}\sqrt{\mathit{\alpha }}$(23)where c_s is the local sound speed, and α is a proportionality constant that measures the strength of the disk turbulence. Ormel et al. (2008) found that, for typical disk turbulence (α = 10^-4), the chondrule aggregates have R ~ 1mm in size. For low turbulence (α = 10^-6), the chondrule aggregates have R ~ 4mm, with many aggregates as large as R ~ 7mm.

For r> 10au, mm-sized aggregates have τ_f> 0.01 and the streaming instability occurs easily. For r< 5au, Fig. 8 imposes some constraints on the plausible scenarios that can produce dense particle clumps. These constraints are illustrated in Fig. 9. In the figure, the “disk mass vs. r” phase space is divided into regions where different disk models can form particle clumps. Broadly speaking, there are three avenues to form planetesimals in the inner solar system:

• 1.

  Figure 8 shows that the streaming instability can be pushed to particles as small as τ_f = 0.003 if the solids can reach sufficiently high mass fractions. There are various disk processes that may enhance Z sufficiently to induce clumping. They include the break-up or ice-dust aggregates (Sirono 2011), as well as large-scale pressure bumps (Johansen et al. 2009a).

• 2.

  An alternative to high particle concentration is low disk turbulence. Ormel et al. (2008) have shown that, for α ~ 10^-6, chondrule aggregates can reach R ~ 4−5 mm. Figure 9 shows how low α can be much more effective than high Z in triggering particle clumps.

• 3.

  Finally, since τ_f is fundamentally a measure of the particle stopping time, as the disk clears and the gas density decreases, τ_f increases for all particles. Figure 9 shows how the disk mass (M_disk) affects the location of the planetesimal forming region. The effect is most dramatic in the inner ~1−2 AU, where particle clumping is not possible until M_disk has decreased.

These scenarios are not exclusive. The formation of planetesimals may well rely on a combination of, say, disk evolution to decrease M_disk and pressure bumps to increase Z. We wish to investigate which combination of (Z, α, M_disk) are compatible with the formation of particle clumps in the disk. To do this, we generalize the Hayashi nebula to allow for different disk masses $\mathrm{\Sigma }\mathrm{=}\mathrm{1700}\mathrm{g}{\mathrm{cm}}^{-2}\hspace{0.17em}{\left(\frac{\mathit{r}}{\mathrm{au}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}\left(\frac{{\mathit{M}}_{\mathrm{disk}}}{{\mathit{M}}_{\mathrm{MMSN}}}\right)\mathit{,}$(24)where M_MMSN is the disk mass for the Hayashi nebula. Equation (8)relates τ_f to Σ, the particle size R, and material density ρ_•. Take a typical material density of ρ_• = 3.6gcm^-3, and we obtain the minimum distance r_pl where asteroids can form by the streaming instability, ${\mathit{r}}_{\mathrm{pl}}\mathrm{=}\mathrm{4.33}\mathrm{au}{\left(\frac{{\mathit{M}}_{\mathrm{disk}}}{{\mathit{M}}_{\mathrm{MMSN}}}\hspace{0.17em}\frac{{\mathit{\tau }}_{\mathrm{f}\mathit{,}\min }}{\mathrm{0.003}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}{\left(\frac{\mathit{R}}{\mathrm{mm}}\right)}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\mathit{,}$(25)where τ_f,min is the minimum Stokes number needed to trigger the streaming instability, which is determined by the particle concentration. In a disk where solids can concentrate to Z = 0.065, Fig. 8 gives a minimum Stokes number of τ_f,min = 0.003. Otherwise, take τ_f,min = 0.006 and Z = 0.04. Equation (25) can be thought of as the location of the “planetesimal line”, or the minimum distance where planetesimals can form. This line moves inward as the disk evolves. Figure 9 shows the location of the planetesimal line for different choices of Z, α, and M_disk.

5. Coagulation for mm-sized particles

Particles with R = 1mm (τ_f = 0.003) may be a critical stage in the formation of planetesimals. These particles lie at the upper limits of coagulation, and the lower limits of the streaming instability. The formation of planets may hinge on the ability of these particles to grow from τ_f = 0.003, and move deeper into the streaming regime, where gravitational collapse can occur (see Eq. (2)and Bai & Stone 2010b). In this section we explore the ability of mm-sized, τ_f = 0.003 to grow faster than the disk accretion timescale.

Our simulations are too coarse to directly measure collision speeds between particles, but with some extrapolation, it is possible to produce robust upper limits. To do this, we note that in a small region of the disk the particle velocities should be randomized and their relative speeds should follow a Maxwell-Boltzmann distribution, $\mathit{f}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\left(\frac{{\mathit{v}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{3}}}\right)\exp \left(\frac{\mathrm{-}{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}{\mathit{a}}^{\mathrm{2}}}\right)\mathrm{·}$(26)

The Maxwell-Boltzmann distribution is parametrized by the characteristic speed a. Given a series of relative speeds $\mathrm{\{}{\mathit{v}}_{\mathit{i}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}$, the Maximum Likelihood Estimate for a is $\mathit{â}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathrm{3}\mathit{n}}\sum _{\mathit{i}\mathrm{=}\mathrm{0}}^{\mathit{n}}{\mathit{v}}_{\mathit{i}}^{\mathrm{2}}}\mathit{.}$(27)The strategy, then, is to extract a collection of particle pairs that are closer than some distance D and compute â from the relative speeds. This is shown in Fig. 10 for several choices of D, and for three different simulations with τ_f = 0.003. If we extrapolate â for D → 0, we obtain a rough estimate of the Maxwellian distribution of collision speeds. If we fit a linear relation, â = â₀ + mD, the extrapolated value of â is â₀ = −0.0003 ± 0.05 cm/s. Since the true â₀ cannot be negative, we take the range $\mathrm{0}\mathrm{cm}{\mathrm{s}}^{-1}\mathrm{\le }\mathrm{â}\mathrm{0}\mathrm{\le }\mathrm{0.1}\mathrm{cm}{\mathrm{s}}^{-1}\mathit{.}$(28)Figure 10 shows the extrapolation of â for â₀ in this interval, and for the 1σ confidence region for the slope m. We are interested in the probability that a collision of two mm-sized particles will result in sticking. For a material density of ρ_• ~ 3.5gcm^-3, the critical speed is around v_crit ~ 0.03 cm/s (Güttler et al. 2010). The probability that v<v_crit is given by the cumulative distribution function $\mathit{F}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{v}}\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{x}\mathrm{=}\mathrm{erf}\left(\frac{\mathit{v}}{\mathit{a}\sqrt{\mathrm{2}}}\right)\mathrm{-}\frac{\mathit{v}}{\mathit{a}}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\exp \left(\frac{\mathrm{-}{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}{\mathit{a}}^{\mathrm{2}}}\right)\mathrm{·}$(29)Figure 11 shows the probability of a sticking collision F(v_crit) for the â values from Fig. 10. As D → 0, the steady decrease in particle speeds leads to a significant increase in the fraction of collisions that result in sticking.

Figure 12 shows the Maxwellian distributions for â = 0.3 cm/s (typical value for D = 10^-5 H) and for â = 0.05 cm/s (mid-range value for â₀). The bulk of the collisions result in bouncing, and a small fraction result in sticking. For comparison, the figure also shows the root-mean-squared vertical speed v_rms for all the particles in the grid. The figure shows how v_rms greatly over-estimates the particle collision speeds.

Finally, we use F(v_crit) to determine the typical time needed before a 1 mm grain has a sticking collision. The key question is whether the particles can grow faster than the disk evolution timescale. The time needed for a sticking collision is given by ${\mathit{t}}_{\mathrm{stick}}\mathrm{=}{\mathrm{\left(}\mathit{\pi }\mathrm{\right(}\mathrm{2}\mathit{r}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathit{v̅}\hspace{0.17em}{\mathit{n}}_{\mathrm{p}}\hspace{0.17em}\mathit{F}\mathrm{(}{\mathit{v}}_{\mathrm{crit}}{\mathrm{\left)}}^{\mathrm{\right)}}}^{-1}\mathit{,}$(30)where n_p is the particle number density, and is the mean relative particle speed. In these simulations, n_p ~ 4 × 10^-3m^-3 at the midplane, and for a Maxwellian, $\mathit{v̅}\mathrm{=}\mathrm{2}\hspace{0.17em}\mathit{a}\hspace{0.17em}\sqrt{\mathrm{2}\mathit{/}\mathit{\pi }}$. Figure 13 shows the results for t_stick. As D → 0, we get t_stick ~ 10³−10⁵ years, which is a small fraction of the disk lifetime of ~10⁷ years. We repeated this experiment for several particle concentrations from Z = 0.01 to 0.1. Figure 14 shows that â< 0.1cms^-1 is a robust upper bound, largely independent of Z. We conclude that the particle growth timescale t_stick< 10⁵ years is also robust. Therefore, we conclude that τ_f = 0.003 particles can double in mass within a small fraction of the disk lifetime. As the particles grow, the streaming instability becomes more effective and the density of the particle clumps increases.

6. Conclusion

The formation of planetesimals remains a difficult problem. The traditional bottom-up coagulation model faces important challenges from theory and observation. There is growing evidence that planetesimals must form top-down, from the gravitational collapse of dense clumps of solids (e.g., Morbidelli et al. 2009). In this paper we present the results from over 100 hydrodynamic simulations where we modeled the dynamics of solid particles inside a protoplanetary disk. We model particle sizes from submillimeter-sized chondrules to meter-sized rocks, and evaluate the effect of particle concentration and gas pressure gradient. Our key results can be summarized as:

• 1.

  We present the particle properties as a size vs. concentration phase space, and we map the region where particle clumps can be expected to form (Fig. 8). We measure the particle size as τ_f = Ω_kt_f, where t_f is the friction time and Ω_k is the Keplerian frequency. We measure particle concentration as Z = Σ_solids/ Σ_total, where Σ denotes surface density. We find that particles as small as τ_f = 0.003 can participate in the streaming instability and form dense clumps at Z ~ 0.065. Although this solid concentration is quite high, it is smaller than what can be produced by the break-up of ice-dust aggregates near the snow line (Sirono 2011). More generally, high Z values may be achieved by the radial drift of solids from the outer disk into the inner disk (Youdin & Shu 2002), and by large scale pressure bumps or vortices (Johansen et al. 2009a).

• 2.

  We find constraints on asteroid formation models and present them in Fig. 9. We map the region in the M_disk vs. r phase space that is consistent with the formation of dense particle clumps. The location of the asteroid-forming region is primarily a function of disk mass and turbulent viscosity, as low α leads to larger chondrule aggregates Ormel et al. (2008). After that, the most important factor is the disk’s ability to increase the solid concentration locally.

• 3.

  We estimate the probability of sticking collisions between mm-sized particles (τ_f = 0.003). We find a robust upper limit on the timescale for these particles to stick and grow to larger sizes (~10⁵ years) that is significantly smaller than the lifetime of the disk (~10⁷ years). A high concentration of particles is not required for this result. As the particles grow, the streaming instability becomes more effective, leading to the formation of planetesimals.

Altogether we find that particle concentration by the streaming instability provides a viable path to forming asteroids directly from mm-sized chondrules, particularly if weak turbulence facilitates the growth of chondrule aggregates with sizes of a few millimeters.

Online material

Appendix A: A measure of particle clumping

Particle clumping is difficult to measure objectively. Previous authors have measured the peak particle density, but this is a very local measure that contains limited information about the particle behavior. In particular when small particles (τ_f ≤ 0.003) form filaments that are clearly visible in a spacetime diagram (Fig. 4), no single grid cell achieves a very high density. At the same time, large particles (τ_f> 3) have may have a few cells with very high density, but the density peaks are unstable and short lived.

Figure A.1 illustrates our approach. For each time step, we take the particle surface density Σ_solid (top row). Particle overdensities are clearly visible on this plot. We then average Σ_solid over 25 orbits (⟨ Σ_solid ⟩ _t, second row). This serves to to smooth out overdensities that are unstable, short-lived, or have too much radial drift. We chose 25 orbits because it is half of the 50-orbit timescale for Z to increase. We considered using the radial speed to separate clumps that are drifting into the star. The appeal of radial speed is that it provides an instantaneous measure. However, radial speed does not capture the longevity of a particle clump. For example, as in a traffic jam, the pattern speed of the clump may be lower than the particle drift speed. At the same time, particles with low radial speed may produce short-lived clumps. Therefore, it is necessary to compare the location of a clump across time. Averaging Σ_solid over time is a simple way to achieve this goal.

Finally, we sort the grid cells of ⟨ Σ_solid ⟩ _t in order from highest to lowest density, and compute the cumulative distribution F of the sorted grid cells. If the particle distribution is uniform, then F forms a straight diagonal line G from (0, 0) to (1, 1). We use the difference F−G as a large-scale measure of clumping. The Kolmogorov-Smirnov (Smirnov 1948; Wall & Jenkins 2003) test gives a standard way to compare F and G. It computes a p-value that measures probability that the data set that produced F follows the distribution G$\begin{array}{ccc}& & \mathit{Q}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{2}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{j}\mathrm{-}\mathrm{1}}\exp \mathrm{(}\mathrm{-}\mathrm{2}{\mathit{z}}^{\mathrm{2}}{\mathit{j}}^{\mathrm{2}}\mathrm{)}\\ & & \mathit{p}\mathrm{=}\mathit{Q}\mathrm{\left(}\mathit{D}\sqrt{\mathit{n}}\mathrm{\right)}\mathit{,}\end{array}$where D is the maximum distance between F and G, and n is the number of independent observations. We caution the reader against interpreting these p-values as strict probabilities; we only use p as a convenient metric to measure clumping. With that in mind, we divide the p-values into four broad categories:

• Clumping is unlikely:

  p ≥ 0.45

• Clumping is somewhat likely:

  0.25 ≤ p< 0.45

• Clumping is likely:

  0.10 ≤ p< 0.25

• Clumping is very likely:

  p< 0.10.

Figure 8 shows the result of this analysis across our entire range of simulations. The figure is broadly consistent with the spacetime diagrams in Figs. 4 and 5.

Appendix A.1: Alternate measures of particle clumping

To confirm our results, we devised an alternate method to measure particle clumping, based as much as possible on very different principles. Our second method is to perform a chi-squared test on ⟨ ρ_p ⟩ _zt, without sorting the grid cells. We use the reduced chi-squared statistic, ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathit{k}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{k}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{(}{\mathit{O}}_{\mathit{j}}\mathrm{-}{\mathit{E}}_{\mathit{j}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\mathit{,}$(A.3)where k = n − 1 is the degrees of freedom, O_j is the jth observation (value of ⟨ ρ_p ⟩ _zt on grid cell j), E_j = ⟨ ρ_p ⟩ is the jth expected density for a uniform distribution, and σ² is the expected variance between different measurements. If ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{>}\mathrm{1}$, that indicates that the difference between observations and the model is greater than can be explained by the variance σ². This can be because the model is incorrect, or because the true variance between measurements is greater than σ². We divide ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ into four intervals:

• Clumping is unlikely:

  ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{<}\mathrm{0.7}$

• Clumping is somewhat likely:

  $\mathrm{0.7}\mathrm{\le }{\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{<}\mathrm{1.0}$

• $\mathrm{1.0}\mathrm{\le }{\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{<}\mathrm{1.3}$

• Clumping is very likely:

  ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\ge }\mathrm{1.3.}$

That leaves σ² as the only unspecified parameter. When the simulations start, σ² is almost zero, giving unrealistically high ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ values. As sedimentation proceeds, the value of σ² increases. We opted for a σ² that is typical for the last 25 orbits of the sedimentation phase. Because ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ is sensitive to σ², there is necessarily some amount of subjectivity in σ². In the end, we chose σ² so that the τ_f = 0.003 simulations wold give the same result as in Fig. 8. Therefore, the true test is whether the ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ method remains consistent with the KS method across all the other particle sizes.

Figure A.2 shows the final results for the ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ method. The similarity between this figure and Fig. 8 is striking. Considering how different the two techniques are, the agreement between Figs. 8 and A.2 indicates that our core results are robust.

Appendix B: Convergence tests

Appendix B.1: Three-dimensional box

Our simulations are axisymmetric, they neglect the effect of the Kelvin-Helmholtz Instability (KHI). Whether the KHI is present is dictated by the Richardson number (Chandrasekhar 1961), $\begin{array}{ccc}\mathrm{Ri}\mathrm{=}\frac{\mathit{g}}{\mathit{\rho }}\frac{\mathrm{(}\mathit{\partial \rho }\mathit{/}\mathit{\partial z}\mathrm{)}}{\mathrm{(}\mathit{\partial u}\mathit{/}\mathit{\partial z}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}& & \end{array}$where u is the gas speed, g = Ω²z is the vertical gravitational acceleration, and ρ is the effective fluid density (treating gas and solids as a single fluid). Bai & Stone (2010a) have shown that the streaming instability is able to keep the Richardson number above the critical value needed for the KHI, so that the KHI is absent in this problem. To confirm this, we performed two 3D simulations with τ_f = 0.03 and Δ = 0.05. The experiment performed in our 3D simulation is slightly different from that of our 2D simulations. Instead of increasing Z gradually, we only test the behavior of the system at three discrete values: Z = 0.01, Z = 0.02 and Z = 0.03. The reason for this is that 3D simulations are very computationally expensive, especially for small particles, so that a long-term simulation is prohibitive. Our 3D runs have three phases:

• Phase I: each run starts with Z = 0.01, and Z is held constant for 30 orbits.

• Phase II: the value of Z is increased quickly. In run “3D.Z2” we increase Z to Z = 0.02, and in run “3D.Z3” we increase Z to Z = 0.03.

• Phase III: the value of Z is once again held constant for the remainder of the run.

Figure B.1 shows the spacetime diagrams for both runs. The figure shows that both runs produce visible particle clumps for Z ≥ 0.02 and no clumps for Z = 0.01. It is also clear that clumps form more readily with Z = 0.03. These results are consistent with our 2D simulations (Fig. 4). For particles smaller than τ_f = 0.03, the small integration steps make 3D simulations prohibitive. However, as explained in Sect. 2.6, the long response time of very small particles (τ_f ≤ 0.0064) means that our simulations will overestimate Z; so our results can be considered robust upper bounds for this size range.

Figure B.2 gives a side view (x − z axis) and a top view (x − y axis) of the 3D.Z2 run at t = 150 orbits. The clumps and their filamentary structure are very prominent. An important feature of the 3D runs is that they form visible clumps more easily than the 2D runs. Figure B.3 shows the formation of clumps in the two 3D runs and one of the 2D runs. This figure shows that the 3D runs are consistent with the 2D run. In particular, the clumps are readily visible at Z = 0.02 and the 3D runs show no evidence of additional stirring from the KHI.

Figure B.4 shows the maximum particle density (ρ_p,max) for the two 3D simulations and one of the 2D simulations. One salient feature of the figure is that in the early phase of the simulations, before clumping occurs, the 3D runs consistently show a peak density 2−3 times greater than the 2D runs. This occurs because particles can also accumulate along the azimuthal direction. The figure also shows ρ_p,max for the 3D runs after averaging along the azimuthal direction (purple lines). This averaging removes the apparent discrepancy between the 2D and 3D runs. In the later stages the 2D runs reach higher peak densities because they have very high Z values.

Appendix B.2: Resolution

We performed two simulations at 256 × 256 grid resolution with τ_f = 0.003 particles. Figure B.5 shows the spacetime diagrams for both sets of simulations. Since the behavior of the 256 × 256 runs is almost identical to the corresponding 128 runs, we conclude that, in terms of resolution, our simulations are fully converged.

---

Acknowledgments

The authors would like to thank Axel Brandenburg, Andrew Youdin, Lennart Lindegren, Mike Alexandersen, the editor and the anonymous referee for their comments that helped improve the manuscript. We also acknowledge the support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grants 2010-3710 and 2011-3991) and the European Research Council Starting Grant 278675-PEBBLE2PLANET that made this work possible. Computer simulations were performed using the Alarik cluster at Lunarc Center for Scientific and Technical Computing at Lund University. Some simulation hardware was purchased with grants from the Royal Physiographic Society of Lund.

References

All Figures