Issue 
A&A
Volume 579, July 2015



Article Number  A43  
Number of page(s)  20  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201425120  
Published online  24 June 2015 
How to form planetesimals from mmsized chondrules and chondrule aggregates^{⋆}
Lund ObservatoryDepartment of Astronomy and Theoretical Physics, Lund
University,
Box 43,
22100
Lund,
Sweden
email: danielc@astro.lu.se; anders@astro.lu.se; mbd@astro.lu.se
Received: 7 October 2014
Accepted: 23 April 2015
The size distribution of asteroids and Kuiper belt objects in the solar system is difficult to reconcile with a bottomup formation scenario due to the observed scarcity of objects smaller than ~100 km in size. Instead, planetesimals appear to form topdown, with large 100−1000 km bodies forming from the rapid gravitational collapse of dense clumps of small solid particles. In this paper we investigate the conditions under which solid particles can form dense clumps in a protoplanetary disk. We used a hydrodynamic code to model the interaction between solid particles and the gas inside a shearing box inside the disk, considering particle sizes from submillimetersized chondrules to metersized rocks. We found that particles down to millimeter sizes can form dense particle clouds through the runaway convergence of radial drift known as the streaming instability. We made a map of the range of conditions (strength of turbulence, particle massloading, disk mass, and distance to the star) that are prone to producing dense particle clumps. Finally, we estimate the distribution of collision speeds between mmsized particles. We calculated the rate of sticking collisions and obtain a robust upper limit on the particle growth timescale of ~10^{5} years. This means that mmsized chondrule aggregates can grow on a timescale much smaller than the disk accretion timescale (~10^{6}−10^{7} years). Our results suggest a pathway from the mmsized grains found in primitive meteorites to fully formed asteroids. We speculate that asteroids may form from a positive feedback loop in which coagualation leads to particle clumping driven by the streaming instability. This clumping, in turn, reduces collision speeds and enhances coagulation. Future simulations should model coagulation and the streaming instability together to explore this feedback loop further.
Key words: minor planets, asteroids: general / planets and satellites: formation / planets and satellites: terrestrial planets / protoplanetary disks / hydrodynamics / methods: numerical
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2015
1. Introduction
Planetesimals are superkilometersized 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 leftover 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 bottomup 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 transNeptunian 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 bottomup process from an initial population of 0.1 kmsized bodies. In this scenario, the bump at ~100 km is produced by the transition from dispersiondominated runaway growth to a regime dominated by Keplerian shear.
There are also important theoretical difficulties in building planetesimals in a bottomup process. The first hurdle is that small particles in the mmcm size range do not easily stick together to form larger objects. Indeed, the largest grains observed in a protoplanetary disk are mmcm in size (Testi et al. 2014; van der Marel et al. 2013). To the extent that metersized objects do form, they quickly spiral into the central star due to friction with gas orbiting the star at subKeplerian 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, (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, (2)where Ω is the Keplerian frequency, then the particle selfgravity 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 magnetorotational instability (MRI, Balbus & Hawley 1991). But even inside an MRIfree “dead zone” (Gammie 1996), the KelvinHelmholtz 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 longlived 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 decimetersized bodies, and for metersized 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 dmsized 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 smallscale turbulent eddies in the disk, and they concentrate in the highpressure 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 μmsized 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 millimetersized 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 smallparticle and the largeparticle 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 icedust 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 latestage 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 (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 nonzero 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 (4)where ρ is the gas density, R is the particle size, and v_{th} is the mean thermal speed of gas molecules, 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}, The friction time 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 (5)where ρ is the gas density. We commonly express t_{f} in terms of the Keplerian frequency Ω_{k} to obtain the Stokes number The Stokes number is a scalefree measure of the stopping time. Since the disk scale height is H = c_{s}/ Ω_{k}, the Stokes number can be written as (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 (7)With Eq. (6), this gives the particle Stokes number in terms of the local surface density (8)Combined with Eq. (3), we obtain the relation between R and τ_{f} for the MMSN (9)Note that ρ_{•} = 3.5gcm^{3} is a typical density for a mmsized 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 subKeplerian speed. The difference between Keplerian speed v_{k} and and gas speed u_{φ} can be written as (10)where the pressure gradient parameter η is given by (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 scalefree numerical simulations, it is helpful to normalize the headwind speed by the sound speed and write (12)This Δ is a free parameter that measures the strength of the pressure support and the headwind. The pressure is . Equation (7)gives the gas density at the midplane, which combined with Ω_{k} = c_{s}/H allows us to write The isothermal sound speed c_{s} is determined primarily by the local temperature (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 Therefore, the pressure P ∝ r^{3}c_{s} follows a power law with exponent Equation (12)can now be written as (14)Using the parameters from the MMSN and the ChiangYoudin nebula (CY) we get 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 largescale 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 where σ = 2 × 10^{15}cm^{2} is the collision cross section of an H_{2} molecule, and n is the particle number density, which is 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^{1} to simulate the dynamics of gas and solid particles in a 2dimensional vertical slice of the protoplanetary disk (Youdin & Johansen 2007). We model a box of dimensions 0.2H × 0.2H in the radialvertical 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 “superparticles” representing the solids. Each superparticle 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 overestimate 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 (17)where 0.2H is the width of the box and v_{r} is the particle radial velocity, which is given by With this one can show that t_{cross}< 50 orbits , 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.
Fig. 1 Particle scale height H_{p} (top; see Eq. (20)), and rootmeansquare of the vertical component of the particle speed v_{rms} (bottom; see Eq. (21)) at the end of the sedimentation phase. The sound speed is set to c_{s} = 800ms^{1}. Half of the simulations (red squares) have a standard pressure support (Δ = 0.05) and the other half (blue triangles) have a reduced pressure support (Δ = 0.025). For particles smaller than τ_{f} ~ 0.3, lower Δ results in stronger sedimentation (lower H_{p}). 

Open with DEXTER 
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 rootmeansquare of the particle positions (z coordinate) and vertical speed, 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.
Fig. 2 Particle scale height H_{p} for small particles (size τ_{f} ≤ 0.002, R ≤ 1.5 mm) over the course of the simulation. The top figure shows the entire simulation, including the sedimentation phase. After 250 orbits, H_{p} has reached a steady state and the numerical experiment begins. The bottom figure shows only the gas removal phase. Gas is removed exponentially. In the bottom plot, the simulation time is replaced by the solid concentration Z = Σ_{solid}/ Σ_{total}. For particles in this size scale, the scale height is determined primarily by Δ, and secondly by Z. 

Open with DEXTER 
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 overdensities 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.
Fig. 3 Particle density ρ_{p} as a function of solid concentration (left) and height (right). The solid concentration is Z = Σ_{solid}/ Σ_{total}. The left figure shows ρ_{p} at the midplane (height z = 0) as a function of Z. The right figure shows ρ_{p} at Z = 0.08 as a function of the z coordinate. The particle density is measured in terms of the initial gas density ρ_{0}. All the simulations have particle size τ_{f} = 0.002. Two runs have a standard pressure support (Δ = 0.05) and two have a reduced pressure support (Δ = 0.025). Two runs have 128 × 128 grid resolution and two are 256 × 256. There is no evidence of the density cusp predicted by Sekiya (1998) and Youdin & Shu (2002). The Δ = 0.05 and Δ = 0.025 runs should have clear cusps after Z ~ 0.10 and Z ~ 0.02 respectively (Youdin & Shu 2002). These simulations may have an additional source of mixing that inhibits the formation of the SekiyaYoudinShu cusp. The result is the same for smaller, τ_{f} = 0.001 particles (not shown). 

Open with DEXTER 
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 wellcoupled particles such as these, there is a critical Z beyond which turbulent mixing becomes ineffective, resulting in a highdensity “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 (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 SekiyaYoudinShu 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, longlived 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 pileup process is somewhat analogous to a traffic jam: An initial overdensity of solids has a stronger backreaction on the gas. This reduced Δ, lowers the radial drift, leading to a pileup 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.
Fig. 4 Spacetime diagrams showing the solid surface density Σ_{solid} as a function of the radial coordinate x and simulation time. The right hand axis shows the mean solid concentration Z = ⟨ Σ_{solid} ⟩ / ⟨ Σ_{total} ⟩. The figure shows selected runs with τ_{f} ≤ 0.03. The top row has runs with standard pressure support (Δ = 0.05) and the bottom row has runs with a reduced pressure support (Δ = 0.025). The surface density is shown by a color scale. Some of the runs form visible filaments. For Δ = 0.05 there is a clear trend where lower τ_{f} requires higher Z to form filaments. For Δ = 0.025 the behavior is less predictable – while the particle concentration is higher, the low radial drift speed interferes with the streaming instability. 

Open with DEXTER 
Fig. 5 Same as Fig. 4, but for larger particle sizes (τ_{f} = 0.1 to 3). For τ_{f} = 0.1 and Δ = 0.025 the simulation cannot resolve the minimum solid concentration (Z = Σ_{solid}/ Σ_{total}) that leads to clumping. All simulations begin at Z = 0.005 and that is sufficient to produce visible clumps for that τ_{f} and Δ. Also, for τ_{f} = 1 and Δ = 0.05, the filaments that form are short lived. 

Open with DEXTER 
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 5m 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 nonlinear 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 shortlived and suffer strong radial drift. Therefore, our strategy is to integrate the particle density over a 25orbit interval and use the KolmogorovSmirnov (KS) test to distinguish the resulting particle distribution from a uniform distribution. We test several 25orbit periods, spaced in steps of 10 orbits so that they overlap. The result of the KS test is a pvalue, 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.
Fig. 6 Radial drift speed v_{r} of particles as a function of the particle size τ_{f} and solid concentration Z = Σ_{solid}/ Σ_{total}. Speed is measured in terms of the sound speed c_{s}. All solids experience gas drag. Particles of size τ_{f} ~ 1 experience very rapid drift, which decreases for larger particles. In addition, higher Z forces the gas component to orbit closer to the Keplerian speed, reducing  v_{r} . However, in all cases, the drift speed remains high. For reference, for  v_{r}  = 0.005c_{s} the particles cross the disk after a few thousand years. 

Open with DEXTER 
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 colorcoded: 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.
Fig. 7 Surface density of solids Σ_{solid} as a function of the radial coordinate x, evaluated at two points in time 25 orbits apart (first blue, then green). The top plot is in the suspension regime (particle size τ_{f} = 0.001) and has solid concentration Z = ⟨ Σ_{solid} ⟩ / ⟨ Σ_{total} ⟩ = 0.095. In this regime, there is very little clumping even at high Z. The middle plot is the beginning of the streaming regime (τ_{f} = 0.003) and Z = 0.095. In this regime, clumps form and are stable. The bottom plot is in the radial drift regime (τ_{f} = 3) with Z = 0.05. In this regime particles form large overdensities that experience rapid radial drift. 

Open with DEXTER 
Fig. 8 Region of the particle size vs. concentration phase space where the streaming instability is active. Particle size is measured in the stopping time τ_{f} and the particle concentration is Z = Σ_{solid}/ Σ_{total}. The colors mark the probability that particle clumps can form, where red is “unlikely”, magenta is “somewhat likely”, blue is “likely”, and green is “very likely”. See Appendix A for the precise definitions and methods. When different simulations give different results, the two extreme values are shown. For example, if two simulations give “likely” and one gives “unlikely”, this is expressed as a square that is half red, half green. The symbols indicate the number of simulations available for the result. A triangle means three simulation, a cross means two, and a circle means one. Regions without symbols are extrapolations. The black lines indicate our (subjective) assessment of main regions where the behavior of particles is qualitatively different. 

Open with DEXTER 
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 μmsized 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, (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, mmsized 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 breakup or icedust aggregates (Sirono 2011), as well as largescale 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 (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, (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}.
Fig. 9 Location of the asteroid forming region in terms of disk mass M_{disk} and semimajor axis r. This can also be thought of as a “disk age vs. r” phase space. The disk mass is measured in terms of the minimummass solar nebula (M_{MMSN}). The top plot shows the region where R = 4mm chondrule aggregates can form clumps. The solid and dashed lines are lines of constant particle Stokes number τ_{f} = 0.0032, 0.0056, 0.01, and 0.0178. To form clumps, these particles require a solid concentration of Z = 0.069, 0.038, 0.026, and 0.021 respectively (see Fig. 8). The bottom plot shows the same regions for R = 1mm chondrules. As the disk evolves, and the disk mass drops, the planetesimal forming region moves closer to the star. 

Open with DEXTER 
Fig. 10 Maximum likelihood estimate â of the characteristic speed between solid particles. The particle size is τ_{f} = 0.003 and the solid concentration (Z = Σ_{solid}/ Σ_{total}) is 0.074 for all the simulations. The procedure is to select particle pairs with separation less than D and compute â from Eq. (27). Extrapolating to D → 0 gives a rough estimate of the characteristic collision speed â_{0}. The region in green corresponds to 0 ≤ â_{0} ≤ 2σ_{a} and the 1σ confidence interval for the slope. 

Open with DEXTER 
Fig. 11 Probability that the relative speed between two solid particles lies below the sticking threshold v_{crit} ~ 0.03 cm/s (Güttler et al. 2010). The particle size is τ_{f} = 0.003 and Z = Σ_{solid}/ Σ_{total} = 0.074 in all the simulations. Given a set of particle pairs with separations less than D, one can fit a Maxwellian distribution (Figs. 10 and 12) determine the fraction of speeds below v_{crit}. Extrapolating to D → 0 gives an estimate of the fraction of particle collisions that result in sticking. The region in green corresponds to the 1σ confidence region marked in Fig. 10. 

Open with DEXTER 
5. Coagulation for mmsized 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 mmsized, τ_{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 MaxwellBoltzmann distribution, (26)
Fig. 12 Maxwellian distribution of relative speeds for mmsized particles (τ_{f} = 0.003, pressure support Δ = 0.05 and solid concentration Z = Σ_{solid}/ Σ_{total} = is 0.074). For particles separated by distance d< 10^{5}H (H is the disk scale height) the Maxwellian distribution has characteristic speed â ~ 0.3 cm/s (red). Extrapolating as d → 0 we obtain a conservative estimate of â_{0} = 0.05 cm/s (blue). For reference, we include the sticking, bouncing and fragmenting regions from Güttler et al. (2010), and the rootmeansquared speed v_{rms} for all the particles in the simulation. Note that the shape of the curve is distorted by the log scale; in particular, the great majority of collisions are in the bouncing regime. Also note that v_{rms} is a poor estimate for the the collision speeds. The value â is obtained by the maximum likelihood method. 

Open with DEXTER 
The MaxwellBoltzmann distribution is parametrized by the characteristic speed a. Given a series of relative speeds , the Maximum Likelihood Estimate for a is (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, â = â_{0} + mD, the extrapolated value of â is â_{0} = −0.0003 ± 0.05 cm/s. Since the true â_{0} cannot be negative, we take the range (28)Figure 10 shows the extrapolation of â for â_{0} in this interval, and for the 1σ confidence region for the slope m. We are interested in the probability that a collision of two mmsized 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 (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 (midrange value for â_{0}). The bulk of the collisions result in bouncing, and a small fraction result in sticking. For comparison, the figure also shows the rootmeansquared vertical speed v_{rms} for all the particles in the grid. The figure shows how v_{rms} greatly overestimates 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 (30)where n_{p} is the particle number density, and is the mean relative particle speed. In these simulations, n_{p} ~ 4 × 10^{3}m^{3} at the midplane, and for a Maxwellian, . Figure 13 shows the results for t_{stick}. As D → 0, we get t_{stick} ~ 10^{3}−10^{5} years, which is a small fraction of the disk lifetime of ~10^{7} 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^{5} 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.
Fig. 13 Time between sticking collisions. The figure shows the three 128^{2} runs from Fig. 11 plus two higherresolution runs. In all runs the particle size is τ_{f} = 0.003 and Z = Σ_{solid}/ Σ_{total} = 0.074. We fit a Maxwellian distribution to to all particle pairs with separation less than D. Extrapolating the fit toward small D we estimate the time between sticking collisions. The points on the figure appear to follow a power law t_{stick} ∝ D^{2}. The region in green corresponds to the 1σ confidence region marked in Fig. 10. The figure shows that t_{stick}< 10^{5} years is a robust upper bound on t_{stick}. 

Open with DEXTER 
Fig. 14 Characteristic collision speed â_{0} between solid particles at different points in our simulations, with 1σ error bars, for particle size τ_{f} = 0.003. The solid concentration (Z = Σ_{solid}/ Σ_{total}) has little effect on â_{0}. The value â_{0}< 0.1 cm/s (red dashed line) is a robust upper bound on â_{0}. This corresponds to coagulation timescale of t_{stick}< 10^{5} years, which is significantly less than the disk lifetime (10^{6}−10^{7} years). Note that the fitting method can give negative values because it is a simple linear extrapolation (see main text). 

Open with DEXTER 
6. Conclusion
The formation of planetesimals remains a difficult problem. The traditional bottomup coagulation model faces important challenges from theory and observation. There is growing evidence that planetesimals must form topdown, 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 submillimetersized chondrules to metersized 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} = Ω_{k}t_{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 breakup of icedust 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 asteroidforming 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 mmsized particles (τ_{f} = 0.003). We find a robust upper limit on the timescale for these particles to stick and grow to larger sizes (~10^{5} years) that is significantly smaller than the lifetime of the disk (~10^{7} 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 mmsized chondrules, particularly if weak turbulence facilitates the growth of chondrule aggregates with sizes of a few millimeters.
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 20103710 and 20113991) and the European Research Council Starting Grant 278675PEBBLE2PLANET 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
 Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, MNRAS, 369, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010a, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010b, ApJS, 190, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010c, ApJ, 722, L220 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: Clarendon) [Google Scholar]
 Chiang, E., & Youdin, A. N. 2010, Ann. Rev. Earth Planet. Sci., 38, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Fraser, W. C., & Kavelaars, J. J. 2009, AJ, 137, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, ApJ, 782, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayashi, C. 1981, Progr. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hughes, D. W. 1980, Earth Planet. Sci. Lett., 51, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Jacquet, E. 2014, Icarus, 232, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Klahr, H., & Henning, T. 2006, ApJ, 636, 1121 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Johansen, A., Brauer, F., Dullemond, C., Klahr, H., & Henning, T. 2008, A&A, 486, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009a, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Mac Low, M.M. 2009b, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 547 [Google Scholar]
 Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, A. T., Chiang, E., AsayDavis, X., & Barranco, J. 2010, ApJ, 718, 1367 [NASA ADS] [CrossRef] [Google Scholar]
 Meheut, H., Meliani, Z., Varniere, P., & Benz, W. 2012, A&A, 545, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., Cuzzi, J. N., & Tielens, A. G. G. M. 2008, ApJ, 679, 1588 [NASA ADS] [CrossRef] [Google Scholar]
 Safronov, V. S. 1972, Evolution of the protoplanetary cloud and formation of the earth and planets (Keter Publishing House), 212 [Google Scholar]
 Sekiya, M. 1998, Icarus, 133, 298 [CrossRef] [Google Scholar]
 Shankman, C., Gladman, B. J., Kaib, N., Kavelaars, J. J., & Petit, J. M. 2013, ApJ, 764, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Sheppard, S. S., & Trujillo, C. A. 2010, ApJ, 723, L233 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, J.M., & Chiang, E. 2013, ApJ, 764, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Sirono, S.I. 2011, ApJ, 733, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Smirnov, N. 1948, Ann. Math. Statist., 19, 279 [CrossRef] [Google Scholar]
 Suzuki, T. K., & Inutsuka, S.i. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 339 [Google Scholar]
 van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wall, J. V., & Jenkins, C. R. 2003, Practical Statistics for Astronomers, eds. R. Ellis, J. Huchra, S. Kahn, G. Rieke, & P. B. Stetson (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 2010, ApJ, 722, 1716 [NASA ADS] [CrossRef] [Google Scholar]
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wyatt, M. C. 2008, ARA&A, 46, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
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, shortlived, or have too much radial drift. We chose 25 orbits because it is half of the 50orbit 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 shortlived 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 largescale measure of clumping. The KolmogorovSmirnov (Smirnov 1948; Wall & Jenkins 2003) test gives a standard way to compare F and G. It computes a pvalue that measures probability that the data set that produced F follows the distribution Gwhere D is the maximum distance between F and G, and n is the number of independent observations. We caution the reader against interpreting these pvalues as strict probabilities; we only use p as a convenient metric to measure clumping. With that in mind, we divide the pvalues 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 chisquared test on ⟨ ρ_{p} ⟩ _{zt}, without sorting the grid cells. We use the reduced chisquared statistic, (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 σ^{2} is the expected variance between different measurements. If , that indicates that the difference between observations and the model is greater than can be explained by the variance σ^{2}. This can be because the model is incorrect, or because the true variance between measurements is greater than σ^{2}. We divide into four intervals:

Clumping is unlikely:

Clumping is somewhat likely:


Clumping is very likely:
That leaves σ^{2} as the only unspecified parameter. When the simulations start, σ^{2} is almost zero, giving unrealistically high values. As sedimentation proceeds, the value of σ^{2} increases. We opted for a σ^{2} that is typical for the last 25 orbits of the sedimentation phase. Because is sensitive to σ^{2}, there is necessarily some amount of subjectivity in σ^{2}. In the end, we chose σ^{2} so that the τ_{f} = 0.003 simulations wold give the same result as in Fig. 8. Therefore, the true test is whether the method remains consistent with the KS method across all the other particle sizes.
Figure A.2 shows the final results for the 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.
Fig. A.1 Particle density distribution for three particle sizes. Particle size is measured in the stopping time τ_{f}, where τ_{f} = 0.001, 0.003 and 3 correspond to the suspension regime, streaming regime and radial drift regime respectively. Row a) is the Surface density of solids Σ_{solid} as a function of the radial coordinate x, evaluated at two points in time 25 orbits apart (first blue, then green). Row b) averages Σ_{solid} over the 25orbit interval ⟨ Σ_{solid} ⟩ _{t}. Row c) shows ⟨ Σ_{solid} ⟩ _{t} with the grid cells sorted from highest to lowest density. A steep curve indicates stable particle clumps. Row d) shows the cumulative distribution of ⟨ Σ_{solid} ⟩ _{t} in row c) (solid red curve) along with the cumulative distribution of a perfectly uniform particle distribution (dashed blue line). The distance between these two curves serves as a measure of particle clumping. In these plots the particle concentrations Z = ⟨ Σ_{solid} ⟩ / ⟨ Σ_{total} ⟩ are (left to right) 0.095, 0.095 and 0.05. 

Open with DEXTER 
Fig. A.2 As in Fig. 8, but using the method to estimate the likelihood of clumping instead of the KSderived method. The figure marks the region of the particle size vs. concentration phase space where the streaming instability is active. Particle size is measured in the stopping time τ_{f} and the particle concentration is Z = Σ_{solid}/ Σ_{total}. The colors mark the probability that particle clumps can form, where red is “unlikely”, magenta is “somewhat likely”, blue is “likely”, and green is “very likely”. When different simulations give different results, the two extreme values are shown. The symbols indicate the number of simulations available. The circle, cross and triangle indicate, respectively, one, two and three simulations. Regions without symbols are extrapolations. To facilitate the comparison with Fig. 8, the boundary marked by the black lines has been copied from Fig. 8 without modification. 

Open with DEXTER 
Appendix B: Convergence tests
Appendix B.1: Threedimensional box
Fig. B.1 Spacetime diagram of two 3D runs with particle size τ_{f} = 0.03 and Δ = 0.05. The color indicates the column density Σ_{solid}/ ⟨ Σ_{solid} ⟩ using the same color bar as Figs. 4 and 5. Both runs begin with Z = 0.01 and have Z increased only for a short interval. In one run (left), Z grows to 0.02 and in the other (right) Z grows to 0.03. Clumps are visible for Z ≥ 0.02, consistent with the 2D runs. 

Open with DEXTER 
Fig. B.2 A snapshot of run 3D.Z2 taken at t = 150 orbits (Z = 0.02). Image a) is the view from the side (x − z plane). Image b) is the view from above (x − y plane). The color scale marks the column solid density normalized to the initial gas density (Σ_{solid}/ Σ_{gas,0}). The two images have a different color scale, since the column density in the azimuthal direction is much higher. Note that the particle structure is nearly axisymmetric. 

Open with DEXTER 
Fig. B.3 A sequence of vertical snapshots (xz plane) showing the formation of distinct particle clumps in one of our 2D simulations (left) and our two 3D simulations (middle, right). For each run, the images are taken 10 orbits apart. The time and Z values are shown in each plot. The times were chosen to show the formation of the first distinct clumps, and to give a similar peak density on the final image (ρ_{p,max} = 1.8,1.5,2.1left to right). The color indicates the column density using the same color scale as in Fig. B.2a. The other 2D runs give similar results to the one shown. The 2D run is consistent with the 3D runs. In particular, the 3D runs do not show any evidence of additional turbulence stirring compared to the 2D runs. 

Open with DEXTER 
Fig. B.4 Maximum particle density (ρ_{p,max}) for one of the 2D simulations and the two 3D simulations. The top plot compares run 3D.Z2 (blue) and a 2D run (red). The bottom plot compares run 3D.Z3 (blue) and the same 2D run (red). In each plot the purple line shows ρ_{p,max} for the 3D run after averaging along the azimuthal direction. The vertical dashed lines mark the places where the 3D runs transition from Z = 0.01 to Z = 0.019 (for 3D.Z2) or Z = 0.029 (for 3D.Z3). Similarly, the vertical solid lines mark the progressive increase of Z during the course of the 2D runs. 

Open with DEXTER 
Our simulations are axisymmetric, they neglect the effect of the KelvinHelmholtz Instability (KHI). Whether the KHI is present is dictated by the Richardson number (Chandrasekhar 1961), where u is the gas speed, g = Ω^{2}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 longterm 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.
Fig. B.5 Spacetime diagrams showing the solid surface density Σ_{solid} as a function of the radial coordinate x and simulation time τ_{f} = 0.003 and τ_{f} = 0.010 particles. The surface density is shown by color, using the same color scale as in Figs. 4 and 5. The figure shows both 128 × 128 simulations and the corresponding 256 × 256 simulations. In all simulations Δ = 0.05. The higherresolution runs produce results consistent with the 128 × 128 runs. 

Open with DEXTER 
All Figures
Fig. 1 Particle scale height H_{p} (top; see Eq. (20)), and rootmeansquare of the vertical component of the particle speed v_{rms} (bottom; see Eq. (21)) at the end of the sedimentation phase. The sound speed is set to c_{s} = 800ms^{1}. Half of the simulations (red squares) have a standard pressure support (Δ = 0.05) and the other half (blue triangles) have a reduced pressure support (Δ = 0.025). For particles smaller than τ_{f} ~ 0.3, lower Δ results in stronger sedimentation (lower H_{p}). 

Open with DEXTER  
In the text 
Fig. 2 Particle scale height H_{p} for small particles (size τ_{f} ≤ 0.002, R ≤ 1.5 mm) over the course of the simulation. The top figure shows the entire simulation, including the sedimentation phase. After 250 orbits, H_{p} has reached a steady state and the numerical experiment begins. The bottom figure shows only the gas removal phase. Gas is removed exponentially. In the bottom plot, the simulation time is replaced by the solid concentration Z = Σ_{solid}/ Σ_{total}. For particles in this size scale, the scale height is determined primarily by Δ, and secondly by Z. 

Open with DEXTER  
In the text 
Fig. 3 Particle density ρ_{p} as a function of solid concentration (left) and height (right). The solid concentration is Z = Σ_{solid}/ Σ_{total}. The left figure shows ρ_{p} at the midplane (height z = 0) as a function of Z. The right figure shows ρ_{p} at Z = 0.08 as a function of the z coordinate. The particle density is measured in terms of the initial gas density ρ_{0}. All the simulations have particle size τ_{f} = 0.002. Two runs have a standard pressure support (Δ = 0.05) and two have a reduced pressure support (Δ = 0.025). Two runs have 128 × 128 grid resolution and two are 256 × 256. There is no evidence of the density cusp predicted by Sekiya (1998) and Youdin & Shu (2002). The Δ = 0.05 and Δ = 0.025 runs should have clear cusps after Z ~ 0.10 and Z ~ 0.02 respectively (Youdin & Shu 2002). These simulations may have an additional source of mixing that inhibits the formation of the SekiyaYoudinShu cusp. The result is the same for smaller, τ_{f} = 0.001 particles (not shown). 

Open with DEXTER  
In the text 
Fig. 4 Spacetime diagrams showing the solid surface density Σ_{solid} as a function of the radial coordinate x and simulation time. The right hand axis shows the mean solid concentration Z = ⟨ Σ_{solid} ⟩ / ⟨ Σ_{total} ⟩. The figure shows selected runs with τ_{f} ≤ 0.03. The top row has runs with standard pressure support (Δ = 0.05) and the bottom row has runs with a reduced pressure support (Δ = 0.025). The surface density is shown by a color scale. Some of the runs form visible filaments. For Δ = 0.05 there is a clear trend where lower τ_{f} requires higher Z to form filaments. For Δ = 0.025 the behavior is less predictable – while the particle concentration is higher, the low radial drift speed interferes with the streaming instability. 

Open with DEXTER  
In the text 
Fig. 5 Same as Fig. 4, but for larger particle sizes (τ_{f} = 0.1 to 3). For τ_{f} = 0.1 and Δ = 0.025 the simulation cannot resolve the minimum solid concentration (Z = Σ_{solid}/ Σ_{total}) that leads to clumping. All simulations begin at Z = 0.005 and that is sufficient to produce visible clumps for that τ_{f} and Δ. Also, for τ_{f} = 1 and Δ = 0.05, the filaments that form are short lived. 

Open with DEXTER  
In the text 
Fig. 6 Radial drift speed v_{r} of particles as a function of the particle size τ_{f} and solid concentration Z = Σ_{solid}/ Σ_{total}. Speed is measured in terms of the sound speed c_{s}. All solids experience gas drag. Particles of size τ_{f} ~ 1 experience very rapid drift, which decreases for larger particles. In addition, higher Z forces the gas component to orbit closer to the Keplerian speed, reducing  v_{r} . However, in all cases, the drift speed remains high. For reference, for  v_{r}  = 0.005c_{s} the particles cross the disk after a few thousand years. 

Open with DEXTER  
In the text 
Fig. 7 Surface density of solids Σ_{solid} as a function of the radial coordinate x, evaluated at two points in time 25 orbits apart (first blue, then green). The top plot is in the suspension regime (particle size τ_{f} = 0.001) and has solid concentration Z = ⟨ Σ_{solid} ⟩ / ⟨ Σ_{total} ⟩ = 0.095. In this regime, there is very little clumping even at high Z. The middle plot is the beginning of the streaming regime (τ_{f} = 0.003) and Z = 0.095. In this regime, clumps form and are stable. The bottom plot is in the radial drift regime (τ_{f} = 3) with Z = 0.05. In this regime particles form large overdensities that experience rapid radial drift. 

Open with DEXTER  
In the text 
Fig. 8 Region of the particle size vs. concentration phase space where the streaming instability is active. Particle size is measured in the stopping time τ_{f} and the particle concentration is Z = Σ_{solid}/ Σ_{total}. The colors mark the probability that particle clumps can form, where red is “unlikely”, magenta is “somewhat likely”, blue is “likely”, and green is “very likely”. See Appendix A for the precise definitions and methods. When different simulations give different results, the two extreme values are shown. For example, if two simulations give “likely” and one gives “unlikely”, this is expressed as a square that is half red, half green. The symbols indicate the number of simulations available for the result. A triangle means three simulation, a cross means two, and a circle means one. Regions without symbols are extrapolations. The black lines indicate our (subjective) assessment of main regions where the behavior of particles is qualitatively different. 

Open with DEXTER  
In the text 
Fig. 9 Location of the asteroid forming region in terms of disk mass M_{disk} and semimajor axis r. This can also be thought of as a “disk age vs. r” phase space. The disk mass is measured in terms of the minimummass solar nebula (M_{MMSN}). The top plot shows the region where R = 4mm chondrule aggregates can form clumps. The solid and dashed lines are lines of constant particle Stokes number τ_{f} = 0.0032, 0.0056, 0.01, and 0.0178. To form clumps, these particles require a solid concentration of Z = 0.069, 0.038, 0.026, and 0.021 respectively (see Fig. 8). The bottom plot shows the same regions for R = 1mm chondrules. As the disk evolves, and the disk mass drops, the planetesimal forming region moves closer to the star. 

Open with DEXTER  
In the text 
Fig. 10 Maximum likelihood estimate â of the characteristic speed between solid particles. The particle size is τ_{f} = 0.003 and the solid concentration (Z = Σ_{solid}/ Σ_{total}) is 0.074 for all the simulations. The procedure is to select particle pairs with separation less than D and compute â from Eq. (27). Extrapolating to D → 0 gives a rough estimate of the characteristic collision speed â_{0}. The region in green corresponds to 0 ≤ â_{0} ≤ 2σ_{a} and the 1σ confidence interval for the slope. 

Open with DEXTER  
In the text 
Fig. 11 Probability that the relative speed between two solid particles lies below the sticking threshold v_{crit} ~ 0.03 cm/s (Güttler et al. 2010). The particle size is τ_{f} = 0.003 and Z = Σ_{solid}/ Σ_{total} = 0.074 in all the simulations. Given a set of particle pairs with separations less than D, one can fit a Maxwellian distribution (Figs. 10 and 12) determine the fraction of speeds below v_{crit}. Extrapolating to D → 0 gives an estimate of the fraction of particle collisions that result in sticking. The region in green corresponds to the 1σ confidence region marked in Fig. 10. 

Open with DEXTER  
In the text 
Fig. 12 Maxwellian distribution of relative speeds for mmsized particles (τ_{f} = 0.003, pressure support Δ = 0.05 and solid concentration Z = Σ_{solid}/ Σ_{total} = is 0.074). For particles separated by distance d< 10^{5}H (H is the disk scale height) the Maxwellian distribution has characteristic speed â ~ 0.3 cm/s (red). Extrapolating as d → 0 we obtain a conservative estimate of â_{0} = 0.05 cm/s (blue). For reference, we include the sticking, bouncing and fragmenting regions from Güttler et al. (2010), and the rootmeansquared speed v_{rms} for all the particles in the simulation. Note that the shape of the curve is distorted by the log scale; in particular, the great majority of collisions are in the bouncing regime. Also note that v_{rms} is a poor estimate for the the collision speeds. The value â is obtained by the maximum likelihood method. 

Open with DEXTER  
In the text 
Fig. 13 Time between sticking collisions. The figure shows the three 128^{2} runs from Fig. 11 plus two higherresolution runs. In all runs the particle size is τ_{f} = 0.003 and Z = Σ_{solid}/ Σ_{total} = 0.074. We fit a Maxwellian distribution to to all particle pairs with separation less than D. Extrapolating the fit toward small D we estimate the time between sticking collisions. The points on the figure appear to follow a power law t_{stick} ∝ D^{2}. The region in green corresponds to the 1σ confidence region marked in Fig. 10. The figure shows that t_{stick}< 10^{5} years is a robust upper bound on t_{stick}. 

Open with DEXTER  
In the text 
Fig. 14 Characteristic collision speed â_{0} between solid particles at different points in our simulations, with 1σ error bars, for particle size τ_{f} = 0.003. The solid concentration (Z = Σ_{solid}/ Σ_{total}) has little effect on â_{0}. The value â_{0}< 0.1 cm/s (red dashed line) is a robust upper bound on â_{0}. This corresponds to coagulation timescale of t_{stick}< 10^{5} years, which is significantly less than the disk lifetime (10^{6}−10^{7} years). Note that the fitting method can give negative values because it is a simple linear extrapolation (see main text). 

Open with DEXTER  
In the text 
Fig. A.1 Particle density distribution for three particle sizes. Particle size is measured in the stopping time τ_{f}, where τ_{f} = 0.001, 0.003 and 3 correspond to the suspension regime, streaming regime and radial drift regime respectively. Row a) is the Surface density of solids Σ_{solid} as a function of the radial coordinate x, evaluated at two points in time 25 orbits apart (first blue, then green). Row b) averages Σ_{solid} over the 25orbit interval ⟨ Σ_{solid} ⟩ _{t}. Row c) shows ⟨ Σ_{solid} ⟩ _{t} with the grid cells sorted from highest to lowest density. A steep curve indicates stable particle clumps. Row d) shows the cumulative distribution of ⟨ Σ_{solid} ⟩ _{t} in row c) (solid red curve) along with the cumulative distribution of a perfectly uniform particle distribution (dashed blue line). The distance between these two curves serves as a measure of particle clumping. In these plots the particle concentrations Z = ⟨ Σ_{solid} ⟩ / ⟨ Σ_{total} ⟩ are (left to right) 0.095, 0.095 and 0.05. 

Open with DEXTER  
In the text 
Fig. A.2 As in Fig. 8, but using the method to estimate the likelihood of clumping instead of the KSderived method. The figure marks the region of the particle size vs. concentration phase space where the streaming instability is active. Particle size is measured in the stopping time τ_{f} and the particle concentration is Z = Σ_{solid}/ Σ_{total}. The colors mark the probability that particle clumps can form, where red is “unlikely”, magenta is “somewhat likely”, blue is “likely”, and green is “very likely”. When different simulations give different results, the two extreme values are shown. The symbols indicate the number of simulations available. The circle, cross and triangle indicate, respectively, one, two and three simulations. Regions without symbols are extrapolations. To facilitate the comparison with Fig. 8, the boundary marked by the black lines has been copied from Fig. 8 without modification. 

Open with DEXTER  
In the text 
Fig. B.1 Spacetime diagram of two 3D runs with particle size τ_{f} = 0.03 and Δ = 0.05. The color indicates the column density Σ_{solid}/ ⟨ Σ_{solid} ⟩ using the same color bar as Figs. 4 and 5. Both runs begin with Z = 0.01 and have Z increased only for a short interval. In one run (left), Z grows to 0.02 and in the other (right) Z grows to 0.03. Clumps are visible for Z ≥ 0.02, consistent with the 2D runs. 

Open with DEXTER  
In the text 
Fig. B.2 A snapshot of run 3D.Z2 taken at t = 150 orbits (Z = 0.02). Image a) is the view from the side (x − z plane). Image b) is the view from above (x − y plane). The color scale marks the column solid density normalized to the initial gas density (Σ_{solid}/ Σ_{gas,0}). The two images have a different color scale, since the column density in the azimuthal direction is much higher. Note that the particle structure is nearly axisymmetric. 

Open with DEXTER  
In the text 
Fig. B.3 A sequence of vertical snapshots (xz plane) showing the formation of distinct particle clumps in one of our 2D simulations (left) and our two 3D simulations (middle, right). For each run, the images are taken 10 orbits apart. The time and Z values are shown in each plot. The times were chosen to show the formation of the first distinct clumps, and to give a similar peak density on the final image (ρ_{p,max} = 1.8,1.5,2.1left to right). The color indicates the column density using the same color scale as in Fig. B.2a. The other 2D runs give similar results to the one shown. The 2D run is consistent with the 3D runs. In particular, the 3D runs do not show any evidence of additional turbulence stirring compared to the 2D runs. 

Open with DEXTER  
In the text 
Fig. B.4 Maximum particle density (ρ_{p,max}) for one of the 2D simulations and the two 3D simulations. The top plot compares run 3D.Z2 (blue) and a 2D run (red). The bottom plot compares run 3D.Z3 (blue) and the same 2D run (red). In each plot the purple line shows ρ_{p,max} for the 3D run after averaging along the azimuthal direction. The vertical dashed lines mark the places where the 3D runs transition from Z = 0.01 to Z = 0.019 (for 3D.Z2) or Z = 0.029 (for 3D.Z3). Similarly, the vertical solid lines mark the progressive increase of Z during the course of the 2D runs. 

Open with DEXTER  
In the text 
Fig. B.5 Spacetime diagrams showing the solid surface density Σ_{solid} as a function of the radial coordinate x and simulation time τ_{f} = 0.003 and τ_{f} = 0.010 particles. The surface density is shown by color, using the same color scale as in Figs. 4 and 5. The figure shows both 128 × 128 simulations and the corresponding 256 × 256 simulations. In all simulations Δ = 0.05. The higherresolution runs produce results consistent with the 128 × 128 runs. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.