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/0004-6361/201425120 | |
Published online | 24 June 2015 |
How to form planetesimals from mm-sized chondrules and chondrule aggregates⋆
Lund ObservatoryDepartment of Astronomy and Theoretical Physics, Lund
University,
Box 43,
22100
Lund,
Sweden
e-mail: 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 bottom-up formation scenario due to the observed scarcity of objects smaller than ~100 km in size. Instead, planetesimals appear to form top-down, 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 submillimeter-sized chondrules to meter-sized rocks. We found that particles down to millimeter sizes can form dense particle clouds through the run-away convergence of radial drift known as the streaming instability. We made a map of the range of conditions (strength of turbulence, particle mass-loading, disk mass, and distance to the star) that are prone to producing dense particle clumps. Finally, we estimate the distribution of collision speeds between mm-sized particles. We calculated the rate of sticking collisions and obtain a robust upper limit on the particle growth timescale of ~105 years. This means that mm-sized chondrule aggregates can grow on a timescale much smaller than the disk accretion timescale (~106−107 years). Our results suggest a pathway from the mm-sized 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 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, (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 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 = tf Ωk, which measures the particle friction time tf 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 (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 vrel 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 vth 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 mp is the proton mass (Armitage 2010). It is convenient to write vth in terms of the isothermal sound speed cs,
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 tf in terms of the Keplerian frequency Ωk to obtain the Stokes number
The Stokes number is a scale-free measure of the stopping time. Since the disk scale height is H = cs/ Ω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 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 vk 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 vk. 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
(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 = cs/H allows us to write
The isothermal sound speed cs is determined primarily by the local temperature
(13)where k is the Boltzmann constant, T is the temperature, mp 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-3cs follows a power law with exponent
Equation (12)can now be written as
(14)Using the parameters from the MMSN and the Chiang-Youdin 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 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 where σ = 2 × 10-15cm2 is the collision cross section of an H2 molecule, and n is the particle number density, which is
Using Eq. (7), along with H = cs/ Ω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 cs ≡ Ω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 Code1 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 tcross, which is (17)where 0.2H is the width of the box and vr is the particle radial velocity, which is given by
With this one can show that tcross< 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 Hp (top; see Eq. (20)), and root-mean-square of the vertical component of the particle speed vrms (bottom; see Eq. (21)) at the end of the sedimentation phase. The sound speed is set to cs = 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 Hp). |
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 (Hp) and the degree of turbulence using the root-mean-square of the particle positions (z coordinate) and vertical speed, where zi is the z coordinate of the ith particle, and vz,i is the vertical component of the particle’s velocity. Figure 1 shows Hp and vrms 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 Hp 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, Hp 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. |
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 Hp, including the sedimentation phase. The long sedimentation phase ensures that Hp 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 Sekiya-Youdin-Shu cusp. The result is the same for smaller, τf = 0.001 particles (not shown). |
Figure 2 shows that, for τf< 0.003, the particle size has very little effect on Hp. Instead, Hp 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 (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.
![]() |
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. |
![]() |
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. |
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 | vr | ~ 0.01cs. 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.
![]() |
Fig. 6 Radial drift speed vr 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 cs. 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 | vr |. However, in all cases, the drift speed remains high. For reference, for | vr | = 0.005cs the particles cross the disk after a few thousand years. |
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.
![]() |
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 over-densities that experience rapid radial drift. |
![]() |
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. |
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, (23)where cs 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 (Mdisk) 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 Mdisk has decreased.
These scenarios are not exclusive. The formation of planetesimals may well rely on a combination of, say, disk evolution to decrease Mdisk and pressure bumps to increase Z. We wish to investigate which combination of (Z, α, Mdisk) 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 MMMSN 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 rpl 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 Mdisk.
![]() |
Fig. 9 Location of the asteroid forming region in terms of disk mass Mdisk 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 minimum-mass solar nebula (MMMSN). 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. |
![]() |
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. |
![]() |
Fig. 11 Probability that the relative speed between two solid particles lies below the sticking threshold vcrit ~ 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 vcrit. 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. |
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, (26)
![]() |
Fig. 12 Maxwellian distribution of relative speeds for mm-sized particles (τf = 0.003, pressure support Δ = 0.05 and solid concentration Z = Σsolid/ Σtotal = is 0.074). For particles separated by distance d< 10-5H (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 root-mean-squared speed vrms 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 vrms is a poor estimate for the the collision speeds. The value â is obtained by the maximum likelihood method. |
The Maxwell-Boltzmann 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 mm-sized particles will result in sticking. For a material density of ρ• ~ 3.5gcm-3, the critical speed is around vcrit ~ 0.03 cm/s (Güttler et al. 2010). The probability that v<vcrit is given by the cumulative distribution function
(29)Figure 11 shows the probability of a sticking collision F(vcrit) 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 â0). 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 vrms for all the particles in the grid. The figure shows how vrms greatly over-estimates the particle collision speeds.
Finally, we use F(vcrit) 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 np is the particle number density, and
is the mean relative particle speed. In these simulations, np ~ 4 × 10-3m-3 at the midplane, and for a Maxwellian,
. Figure 13 shows the results for tstick. As D → 0, we get tstick ~ 103−105 years, which is a small fraction of the disk lifetime of ~107 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 tstick< 105 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 1282 runs from Fig. 11 plus two higher-resolution 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 tstick ∝ D2. The region in green corresponds to the 1σ confidence region marked in Fig. 10. The figure shows that tstick< 105 years is a robust upper bound on tstick. |
![]() |
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 tstick< 105 years, which is significantly less than the disk lifetime (106−107 years). Note that the fitting method can give negative values because it is a simple linear extrapolation (see main text). |
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 = Ωktf, where tf 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 Mdisk 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 (~105 years) that is significantly smaller than the lifetime of the disk (~107 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 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 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, (A.3)where k = n − 1 is the degrees of freedom, Oj is the jth observation (value of ⟨ ρp ⟩ zt on grid cell j), Ej = ⟨ ρ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 25-orbit 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. |
![]() |
Fig. A.2 As in Fig. 8, but using the |
Appendix B: Convergence tests
Appendix B.1: Three-dimensional 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. |
![]() |
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. |
![]() |
Fig. B.3 A sequence of vertical snapshots (x-z 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. |
![]() |
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. |
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), where u is the gas speed, g = Ω2z 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.
![]() |
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 higher-resolution runs produce results consistent with the 128 × 128 runs. |
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
- 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 [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 [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 [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., Asay-Davis, 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 [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]
All Figures
![]() |
Fig. 1 Particle scale height Hp (top; see Eq. (20)), and root-mean-square of the vertical component of the particle speed vrms (bottom; see Eq. (21)) at the end of the sedimentation phase. The sound speed is set to cs = 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 Hp). |
In the text |
![]() |
Fig. 2 Particle scale height Hp 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, Hp 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. |
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 Sekiya-Youdin-Shu cusp. The result is the same for smaller, τf = 0.001 particles (not shown). |
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. |
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. |
In the text |
![]() |
Fig. 6 Radial drift speed vr 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 cs. 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 | vr |. However, in all cases, the drift speed remains high. For reference, for | vr | = 0.005cs the particles cross the disk after a few thousand years. |
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 over-densities that experience rapid radial drift. |
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. |
In the text |
![]() |
Fig. 9 Location of the asteroid forming region in terms of disk mass Mdisk 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 minimum-mass solar nebula (MMMSN). 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. |
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. |
In the text |
![]() |
Fig. 11 Probability that the relative speed between two solid particles lies below the sticking threshold vcrit ~ 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 vcrit. 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. |
In the text |
![]() |
Fig. 12 Maxwellian distribution of relative speeds for mm-sized particles (τf = 0.003, pressure support Δ = 0.05 and solid concentration Z = Σsolid/ Σtotal = is 0.074). For particles separated by distance d< 10-5H (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 root-mean-squared speed vrms 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 vrms is a poor estimate for the the collision speeds. The value â is obtained by the maximum likelihood method. |
In the text |
![]() |
Fig. 13 Time between sticking collisions. The figure shows the three 1282 runs from Fig. 11 plus two higher-resolution 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 tstick ∝ D2. The region in green corresponds to the 1σ confidence region marked in Fig. 10. The figure shows that tstick< 105 years is a robust upper bound on tstick. |
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 tstick< 105 years, which is significantly less than the disk lifetime (106−107 years). Note that the fitting method can give negative values because it is a simple linear extrapolation (see main text). |
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 25-orbit 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. |
In the text |
![]() |
Fig. A.2 As in Fig. 8, but using the |
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. |
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. |
In the text |
![]() |
Fig. B.3 A sequence of vertical snapshots (x-z 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. |
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. |
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 higher-resolution runs produce results consistent with the 128 × 128 runs. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.