Turbulence Sets the Length Scale for Planetesimal Formation: Local 2D Simulations of Streaming Instability and Planetesimal Formation

One of the classical ideas in planet formation is to form planetesimals in a direct gravitational collapse of pebble clouds (Safronov 1969; Goldreich & Ward 1973). Most planetesimals were incorporated into planetary bodies, yet the asteroids, Kuiper Belt objects (KBOs), and comets are believed to be leftovers from the initial plethora of planetesimals. A study of these minor bodies in the solar system is therefore a key to understanding planetesimal and ultimately planet formation. One problem is to subtract 4.5 billion years of collisional evolution of planetesimals from the size distribution found today.

Recent observational work (Delbo' et al. 2017) identified a group of asteroids that clearly did not originate as collisional fragments of larger ones. Members of this asteroid group are all larger than 35 km, with a most likely diameter of ∼100 km, confirming the previous assumption that planetesimals are born big (Morbidelli et al. 2009) and that the 100 km bump in the size–frequency distribution (SFD) is primordial. Also the impact size distribution on the surfaces of Pluto and its largest moon Charon, as recently determined in the New Horizons mission (Singer et al. 2019), finds a strong deficiency of KBOs smaller than 1–2 km in size, currently interpreted as an effect of collisional grinding. The bilobed structure of Arrokoth (Stern et al. 2019), with both parts of quite similar material and the general high fraction of binaries among KBOs, is further support for a gravitational collapse scenario for planetesimal formation. As Nesvorný et al. (2019) have shown, binaries formed in the streaming instability scenario have a specific distribution of inclinations, which turns out to be a good match to observed binaries.

The gravitational collapse of pebble clouds, i.e., accumulations of about centimeter-sized solid material, in the solar nebula is indeed a rapid and efficient route to form planetesimals (Johansen et al. 2006, 2007a). The definition of pebbles in planet formation is thus not simply based on their size, but rather given via their aerodynamic properties expressed as a Stokes number St = τ_fΩ, i.e., the product of aerodynamic "friction" or "stopping" time τ_f and the local Keplerian frequency Ω.

The Stokes number is a central parameter in planetesimal formation, because it determines not only the radial drift but also the vertical sedimentation speed v_z  = StΩz at height z above the midplane, as well as the efficiency with which particles couple to the gas turbulence in the disk, which most likely has a correlation time at large scales on the order of 1/Ω. In addition, the strength of instabilities related to the particle feedback onto the gas relates to the Stokes number (Squire & Hopkins 2018), as, for instance, the streaming instability (Youdin & Goodman 2005).

Our numerical simulations (Birnstiel et al. 2012) of particle growth and disk evolution have shown that St = 0.01–0.1 should be the typical value for the largest expected grains, which then dominate the radial influx of dust and ice grains.

The goal of this paper is to put an absolute length or mass scale onto the planetesimal formation process. Because asteroids and KBOs both show a bump in the SFD at 100 km, this scale should not strongly depend on the distance from the central star, but only on the mass of that star, the properties of pebbles in dimension-free Stokes numbers, and the strength of turbulent diffusion also in a dimension-free version like the α description.

If one compares planetesimal formation to star formation, then this would be like the derivation of a Jeans mass for planetesimals. The Jeans mass for stars also reflects local properties of the gas, like the local gas temperature and density, helping us to understand why certain stellar masses are more likely than others. Yet stars come at both higher and lower than 1 Jeans mass. As with planetesimals, there is an initial mass function for stars, generated by physical processes still under debate until today (Offner et al. 2014). Is it competitive accretion after fragmentation of the initial unstable cloud core or the power spectrum of the turbulence that dictates the shape of the initial mass function, or some combination, or something else? Already the determination of the initial stellar mass function is a problem toward the higher masses, as those stars do not live for too long and it may be difficult to identify multiple systems. The same is true for planetesimals, as today one predominantly only finds the leftovers after incorporating most material into planets and having had the planetesimals undergo a 4.5 billion year lasting collisional evolution. Numerical simulations of star and planetesimal formation both find power laws, and the concept of either Jeans mass or the here-derived concept of a diffusion mass puts a scaling on the numerical experiments.

In the following section we summarize our current understanding of planetesimal formation via gravitational instability via turbulent clustering, via trapping in zonal flows, and in the bump-free streaming instability scenario. In Section 3 we derive our Jeans-like stability criterion based on comparing turbulent diffusion and gravitational contraction timescales. We test this stability criterion successfully in Section 4 by a parameter study on planetesimal formation in a streaming instability scenario. We explore the resulting planetesimal sizes for turbulence values in the solar nebula in Section 5 and find indeed 100 km as a realistic value for an equivalent diameter of planetesimal-forming pebble clouds all over the solar nebula. We briefly discuss and interpret our results in Section 6 and shift all technical details to the appendices of this paper.

It is well understood that if local particle traps in a turbulent disk stop or reduce the radial drift sufficiently, then significant local overdensities can occur that will collapse under their own gravitational attraction (Johansen et al. 2006), a process we coin gravito-turbulent planetesimal formation, as it resembles similar processes in star formation. In Johansen et al. (2007a), simulating a fully turbulent disk by means of the magnetorotational instability creating zonal flows, it was shown that trapping leads to a rapid formation of planetesimals, with and without including particle feedback onto the gas. For simulations without the background turbulence and trap formation, the typical solar metallicity was insufficient, as Kelvin–Helmholtz and streaming instability prevents particles from the necessary sedimentation (Johansen et al. 2007a; Gerbig et al. 2020).

Johansen et al. (2014) showed that only a dust enrichment by a factor of 2–3 above the average abundance of solids in the solar nebula will lead to gravitational collapse in the presence of streaming instability without background turbulence and the formation of particle traps in the solar nebula (Klahr et al. 2018). In that case the local dust enrichment could be produced by photoevaporation of the disk gas in the later stages of disk evolution (Carrera et al. 2017).

If the local dust enrichment was the result of different processes, e.g., dust trapping in vortices or zonal flows (Johansen et al. 2006, 2007a; Dittrich et al. 2013; Raettig et al. 2015), then it would start much earlier in the evolution of the solar nebula (Lenz et al. 2019), but it is currently not clear whether this enrichment in pebbles would first have to trigger a streaming instability, or whether the enrichment would directly fragment once it overcomes diffusion by Kelvin–Helmholtz and streaming instability (Gerbig et al. 2020).

It should also be mentioned that if there was no radial pressure gradient in the disk, no streaming instability would occur and one would directly go into gravitational collapse in a laminar disk even for lower dust enrichment at a so-far-undetermined level (Abod et al. 2019). Streaming instability is no precondition for planetesimal formation, but rather can be the controlling agent of formation efficiency and, as we show in this paper, streaming instability sets a threshold mass for planetesimal formation.

Numerical simulations in which a relatively large volume, possibly larger than a trapping region might be, was globally enriched in particles have lead to the derivation of various characteristic power laws in the size distribution of planetesimals (Johansen et al. 2015; Simon et al. 2016, 2017). In the highest-resolution cases a deficiency of small planetesimals seems to appear, which could correspond to the critical mass for collapse as derived in our paper, yet unfortunately the turbulence diffusivity was not measured in those simulations, which would be the controlling parameter as we will show.

The power laws of the size distribution as found in the literature are rather shallow (Johansen et al. 2015; Simon et al. 2016; Schäfer et al. 2017). Most recently Abod et al. (2019) report a mass distribution ${dN}/{{dM}}_{{\rm{P}}}\propto {M}_{{\rm{P}}}^{-p}$ with a value of p ≈ 1.6 for a range of pressure gradients. If one uses an exponentially truncated power law, even p ≈ 1.3 makes a good fit for the low-mass end. Thus, the mass dominating planetesimal size in these simulations is typically the largest object, i.e., several hundred kilometers, reflecting the total mass of pebbles that was initially put into these simulations. As a result, these simulations explore rather the largest possible planetesimals for a favorable scenario of global enrichment in pebbles.

Interestingly also in the cascade model for planetesimal formation by turbulent clustering, in its latest version described in Hartlep & Cuzzi (2020), a locally enhanced amount of pebbles and a reduced headwind as in a zonal flow are needed to form a sufficient number of planetesimals within the lifetime of the solar nebula. In this model it is not trapping of pebbles in pressure maxima or pebble concentration as effect of the streaming instability, but pebbles get stochastically concentrated in the gas turbulence of the disk. Random concentrations that exceed the local Hill density to form a gravitationally bound planetesimal are rare, and it appears to need a fine-tuning of Stokes number, turbulence strength, and, more importantly, a gas density enhancement factor, a solids enhancement factor, and a scale factor for the headwind parameter to get the desired amount and the sizes of planetesimals within 2 million years (Hartlep & Cuzzi 2020).

The three scenarios discussed above, trapping, streaming, and clustering, are not mutually exclusive, as they all involve a local enhancement of material, involving particle traps like zonal flows and vortices. Thus, they can all be parameterized in planetesimal formation rate as a function of local influx of pebbles as we did in Gerbig et al. (2019) and Lenz et al. (2019, 2020), yet at different conversion efficiencies from pebbles to planetesimals.

It is now our claim that in all three scenarios the final collapse of a cloud is regulated by turbulent diffusion on the scales of self-gravitating pebble clouds. In Figure 1 we sketch the general scenario of pebble trapping and converting them locally into planetesimals by gravitational collapse; this general picture should hold for any gravity-assisted planetesimal formation scenario.

Zoom In Zoom Out Reset image size

Building on that paradigm, we ask now for the minimal mass in pebbles that has to be locally accumulated to trigger gravitational collapse. This is different from the ansatz in Drążkowska et al. (2016) and Schoonenberg et al. (2018), as they ask for the critical density or dust-to-gas ratio in the midplane to trigger streaming instability and thus planetesimal formation, whereas one should ask to concentrate pebbles to Hill density and a critical mass. Yet it is a good criterion to ask for the dust-to-gas ratio of one in the midplane, as a critical stage, because one has to overcome this critical limit on the way to Hill density anyway, and once we are past a dust-to-gas ratio of 1, the effect of turbulence gets diminished anyway (Johansen et al. 2007a)

We can formulate a collapse criterion that requires a critical clump size $r={l}_{{\rm{c}}}$ (i.e., the radius of pebble cloud) to be large enough to contract against any sort of underlying turbulent diffusion D. Cuzzi et al. (2008) also looked into the effect of internal turbulent pressure on the diffusion of pebble clouds as a limiting factor for gravitational collapse, but they argued that the headwind a collapsing pebble cloud experiences is stronger than the effect of global turbulence, with a strength of even α = 10⁻³; thus, they neglected the internal diffusion thereafter (Cuzzi et al. 2010). As pointed out in Hartlep & Cuzzi (2020), the untamed headwind in the solar nebula would make the formation of planetesimals as small as 100 km impossible, so they introduced a headwind reduction factor F_β  = 1/30, arguing that planetesimal formation may occur in a zonal flow. This falls pretty much along our point of view, yet as they decrease the effect of the headwind, the assumption that the headwind is stronger than internal diffusion no longer holds. As we will show later, the inclusion of turbulent diffusion should strongly change the results from Hartlep & Cuzzi (2020).

As a side remark, our simulations (see Section 4) contain the full unreduced headwind, which is an essential part of the streaming instability we study. So we see no severe impact of this headwind on the collapse, which is entirely diffusion controlled.

3.1. A Timescale Argument

The particle cloud collapse criterion that we want to derive resembles the Jeans mass (Jeans 1902) criterion, which is the lowest mass threshold for a molecular cloud core to form stars. It is defined by questioning whether internal gas pressure can be overcome by self-gravity, resulting in the cloud collapsing under its own weight. Thus, the Jeans mass is a function of cloud density and temperature. The situation for planetesimal formation, starting from a self-gravitating pebble cloud, is similar, but not the same. The self-gravity of the particle cloud has to overcome two opposing effects: On one hand, the cloud has to withstand the tidal shear forces exerted by the central star. This is a force strong enough to disrupt comets, like Shoemaker-Levy (Asphaug & Benz 1996) in 1993 during its encounter with Jupiter. On the other hand, gravity has to overcome diffusion from inherent turbulent motions of the gas–dust mixture (Shariff & Cuzzi 2015). Hence, objects smaller than the derived critical size cannot collapse but are diffused by gas turbulence, an effect neglected in the original works of planetesimal formation (Safronov 1969; Goldreich & Ward 1973), which only considered gas-free rms velocities among the particles.

We state below two criteria for a gravitational collapse to occur.

3.1.1. First Criterion

The density of a particle cloud has to be larger than the critical density ρ_Hill (Johansen et al. 2014 and Appendix B) that allows the cloud to withstand tidal shear while in orbit around a star of mass M at a distance of R:

$$\begin{eqnarray}&&{\rho }_{\mathrm{Hill}}=\displaystyle \frac{9}{4\pi }\displaystyle \frac{M}{{R}^{3}}.\end{eqnarray} \tag{ 1 }$$

In Figure 2 we plot the Hill density as a function of distance around a solar-mass star. We compare the value to the minimum-mass solar nebula (MMSN; Hayashi 1981) density profile of the gas (1700 g cm⁻³ at 1 au and a radial slope of −1.5) and find that the Hill density exceeds the gas density by a factor of about 200−80 depending on the location in the disk.

Zoom In Zoom Out Reset image size

We also compare the Hill density to a more modern approach of a nebula that would be able to form the solar system at least in terms of planetesimal distribution (Lenz et al. 2020). In that paper the authors apply the aforementioned pebble flux regulated planetesimal formation rate and constrain initial disk parameters, like disk mass, radial slope, and the best-suited α viscosity, to produce an initial planetesimal distribution for the solar nebula, which fits all known constraints to form planets in the current mixed pebble and planetesimal accretion scenario to form terrestrial planets (Walsh et al. 2011) and cores (Raymond & Izidoro 2017) and for the dynamic evolution as in the Nice model (Morbidelli et al. 2007; Levison et al. 2011), among others.

In the most appealing solar nebula (MASN; Lenz et al. 2020) the disk is more massive (0.1M_⊙) and the local gas density is larger than in the MMSN case. The MASN is also subject to viscous evolution; thus, the gas surface density is shallower than in the MMSN case. The MMSN was always too steep to be explained by viscous evolution. In addition, disks in star-forming regions observed in the submillimeter (Andrews et al. 2010) show radial density slopes as shallow as predicted by viscous modeling.

The MASN has an exponential cutoff radius at 20 au, which was the result of fitting the constraint of a low-mass planetesimal disk for the Nice model (Morbidelli et al. 2007). As also shown in Figure 2, the dust-to-gas ratio in the MASN upon reaching Hill density ε_Hill is initially probably between 10 and 100, which is, according to Lenz et al. (2020), when most planetesimals in terms of bulk mass will form.

Streaming instability is still driving turbulence diffusion even at dust-to-gas ratios up to 1000 as shown in Schreiber & Klahr (2018), and we will come back to this issue when discussing planetesimal sizes as a function of local dust-to-gas ratios.

3.1.2. Second Criterion

If the critical density ${\rho }_{\mathrm{Hill}}$ is reached, the gravity at the cloud surface still has to overcome the turbulent diffusion, characterized by the diffusion coefficient D. Now, in order to form a planetesimal, a particle cloud has to contract faster than the diffusion can disperse it. If the particles were large enough to decouple completely from the gas, then the collapse of a cloud at density ρ_Hill would occur on the free-fall time τ_ff:

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}=\sqrt{\displaystyle \frac{3\pi }{32\,G{\rho }_{\mathrm{Hill}}}}=\sqrt{\displaystyle \frac{{\pi }^{2}}{24{{\rm{\Omega }}}^{2}}}=0.64{{\rm{\Omega }}}^{-1}.\end{eqnarray} \tag{ 2 }$$

In the classical work of gravitational formation of planetesimals (Safronov 1969; Goldreich & Ward 1973), gas drag is neglected and the collapse is controlled by the rms velocity of the colliding particles and the relevant coefficient of restitution. While it is certainly the case that during the collapse of the pebble cloud collisions among the pebbles will eventually dominate (Nesvorný et al. 2010), this is not the case for the onset of collapse (Jansson et al. 2017), and certainly not before the collapse, i.e., as long as diffusion can prevent the collapse. As derived in Appendix E.9, we can safely ignore pebble collisions in our study, as the mean free path for pebbles is larger than the scales of the considered pebble clouds, which is independent from the actual rms velocity. We also derive a critical pebble cloud mass, expressed as an equivalent diameter following the estimates in Nesvorný et al. (2010), ¹ which we plot in Figure 3. We find that only pebble clouds of significantly larger mass (respectively equivalent diameter) than the sizes derived in our paper are subject to collisions during the onset of collapse.

Zoom In Zoom Out Reset image size

Nesvorný et al. (2010) compare the friction time to the collision time, for the onset of collapse of a pebble cloud forming Kuiper Belt binaries of equivalent size of 500 km diameter, and find the collision time to be an order of magnitude shorter than the friction time and therefore neglect gas drag. They assume the pebble cloud to be already virialized and use the virial velocity as rms velocity, yet Jansson et al. (2017) argue that the contraction speed of the pebbled cloud in the presence of gas drag is initially larger than the rms speed; thus, gas friction again wins over collisions for the onset collapse.

Nevertheless, we also added an estimate for the ratio between collision and friction time for our simulations to Appendix E.9 and determined critical pebble cloud masses (respectively, equivalent diameter a_eq) for the MMSN and MASN (see Figure 3), and we find agreement with Nesvorný et al. (2010) that about 500 km is a size to be dominated by collisions in the MMSN, but not in the MASN, where even for 2000 km gas drag may be of relevance. However, we can show that our derived pebble cloud masses are always too small to be dominated by collisions for the onset of collapse.

Translating a given Stokes number into physical grain or pebble diameter results in a range of sizes from submillimeter to several centimeters, depending on their porosity and external conditions like the gas surface density of the disk. In other words, a massive cobweb, a large snowflake, and a small marble can have very different masses and sizes, yet they can still share a common friction time, which is all that matters for both pebble definition and planetesimal formation.

Interestingly, the Stokes numbers of pebbles around 2 au in our models of planetesimal formation in the solar nebula (Lenz et al. 2019, 2020; see Figure 4) are corresponding to solid marbles of several centimeters. Thus, if we assume those big pebbles to be incorporated into planetesimals and material smaller than St = 0.01 to be largely excluded from this process (Lenz et al. 2019), then we would expect St ≈ 0.001 material to be left behind. In the early, more gas-rich stages of the solar nebula, this Stokes number would correspond to roughly 1 cm sized material like calcium-aluminum-rich inclusions (CAIs; Connelly et al. 2012) and at later stages to millimeter-sized pebbles, like the typical chondrules either before or after the flash heating (DeFelice et al. 2019). In that paradigm CAIs and chondrules would be the pebbles that were just a bit too small to be effectively involved in planetesimal formation and the last ones to be subject to pebble accretion. Following Wahlberg Jansson & Johansen (2017), pebbles of St < 1 feel the friction with the gas as the dominant effect controlling the onset of collapse, and the actual contraction time τ_c will become longer than the free-fall time (Shariff & Cuzzi 2015). In the case of the friction time being shorter than the collapse time, the contraction time is inversely proportional to the friction time τ_f, as derived in Appendix C. For St < 0.1 the contraction time is

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=\displaystyle \frac{8}{3{\pi }^{2}}\displaystyle \frac{{\tau }_{\mathrm{ff}}^{2}}{{\tau }_{{\rm{f}}}}.\end{eqnarray} \tag{ 3 }$$

For St > 0.1 the contraction time approaches the free-fall time. As the free-fall time enters this expression as squared, the contraction time is inversely proportional to the actual particle density of the shrinking cloud,

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=\displaystyle \frac{1}{9\mathrm{St}}\displaystyle \frac{{\rho }_{\mathrm{Hill}}}{\rho }\,{{\rm{\Omega }}}^{-1}=\displaystyle \frac{1}{9\mathrm{St}}\displaystyle \frac{r{\left(t\right)}^{3}}{{r}_{0}^{3}}\,{{\rm{\Omega }}}^{-1},\end{eqnarray} \tag{ 4 }$$

which means that once a cloud of radius r₀ is able to contract, the process accelerates with the shrinking cloud radius ∝r(t)³. If the contracting particle cloud of radius r is subject to turbulent motion, then the cloud is diffused on the typical timescale of

$$\begin{eqnarray}&&{\tau }_{{\rm{D}}}=\displaystyle \frac{{r}^{2}}{D}=\displaystyle \frac{1}{\delta }{\left(\displaystyle \frac{r}{H}\right)}^{2}{{\rm{\Omega }}}^{-1},\end{eqnarray} \tag{ 5 }$$

where we use the dimensionless diffusivity δ by scaling D = δHc_s with the vertical disk extent H (aka pressure scale height) and the speed of sound c_s  = HΩ. Regardless of whether the diffusivity measured in δ is due to large-scale gas turbulence in the solar nebula or solely generated by the streaming instability, or some combination of both, this expression holds. We will discuss this issue below when we use real numbers to estimate the resulting mass for planetesimals.

Zoom In Zoom Out Reset image size

The diffusion time scales with r², but the contraction time with r³; thus, the contraction will always win, once started. Once a pebble cloud goes beyond a critical mass, the contraction cannot be halted by diffusion anymore. This is the same effect as the collapse of an isothermal sphere of gas, which cannot be stabilized by the increase of its internal pressure.

Comparing the two timescales of diffusion (Equation (5)) and contraction (Equation (3)) by setting ${\tau }_{{\rm{c}}}={\tau }_{{\rm{D}}}$, we can derive a critical (minimal) cloud radius $r={l}_{{\rm{c}}}$ for a pebble cloud at Hill density to withstand internal diffusion and allow for contraction at

$$\begin{eqnarray}&&{l}_{{\rm{c}}}=\displaystyle \frac{1}{3}\sqrt{\displaystyle \frac{\delta }{\mathrm{St}}}H.\end{eqnarray} \tag{ 6 }$$

This expression can be understood as a critical length for planetesimal formation, similar to the Jeans length in star formation. Clumps of less than the critical density or of radius less than ${l}_{{\rm{c}}}$ cannot collapse, but larger or more massive clumps will (See Figure 5). Moreover, cloud collapse will always set in, once this border of stability is reached; thus, larger or more massive pebble clouds are less likely to form, if accumulation takes longer than the gravitational collapse, thus limiting planetesimal sizes.

Zoom In Zoom Out Reset image size

The classical Jeans length is derived for a homogeneous density distribution, asking when an infinitesimal density perturbation will start to grow and collapse. A derivation of this smallest linear unstable wavelength for pebbles in a Toomre instability fashion gives $\lambda =2\pi {l}_{{\rm{c}}}$ (see Appendix D.1). This means that in the following simulations, in which our simulation size L never covers $L=2\pi {l}_{{\rm{c}}}$, the formation of clumps was triggered by the streaming instability and not by linear growing self-gravity modes, because they do not fit in our simulation domain. The ${l}_{{\rm{c}}}$ criterion is therefore describing the stability of nonlinear density perturbations as created by the streaming instability or any other concentrating effect, when local gas turbulence has to be considered.

3.2. Particle Layer Scale Height

The length scale ${l}_{{\rm{c}}}$ is a more general quantity than we have mentioned so far. If we ask for the vertical scale height of the dust subdisk in the solar nebula upon reaching the Hill density in the midplane, i.e., in equilibrium between vertical diffusion and sedimentation dominated by self-gravity (See Figure 5), we find the functional dependency

$$\begin{eqnarray}&&{\rho }_{d}(z)={\rho }_{H}{\cosh }^{-2}\left(-\displaystyle \frac{z}{\sqrt{2}{l}_{{\rm{c}}}}\right)\approx {\rho }_{H}{e}^{-\tfrac{{z}^{2}}{2{l}_{{\rm{c}}}^{2}}},\end{eqnarray} \tag{ 7 }$$

which for values in z up to one pressure scale height can be approximated with the usual Gaussian distribution of density around the midplane with an error of less than 4% (see Figure D1). This means that our critical length scale is simultaneously the "pressure" scale height of the self-gravitating pebble accumulations. Without self-gravity the scale height of particles h_p (Dubrulle et al. 1995) would be three times larger,

$$\begin{eqnarray}&&{h}_{p}=\sqrt{\displaystyle \frac{\delta }{\mathrm{St}}}H=3{l}_{{\rm{c}}}.\end{eqnarray} \tag{ 8 }$$

With increasing mass of the pebble layer it will also become thinner and thus overproportionally denser.

To test the derived ${l}_{{\rm{c}}}$ criterion, we perform shearing box simulations of self-gravity-induced collapse in a gas and dust mixture starting from fully developed streaming instability turbulence. The idea is that only when a cloud of diameter 2l_c would fit into the simulation box with size L, then collapse should occur. Thus, for the numerical simulations the collapse criterion is

$$\begin{eqnarray}&&L\gt 2{l}_{{\rm{c}}}.\end{eqnarray} \tag{ 9 }$$

We perform simulations for two particle sizes, that is, for St = 0.1 and St = 0.01 particles and choosing an average dust-to-gas ratio of ε₀ = 3, suggesting that trapping and sedimentation have already achieved this level of local dust concentration. We used the Pencil Code in setups based on Schreiber & Klahr (2018) in a radial-azimuthal setup to save computation time. See Table 1 for the simulation parameters. All simulations have the same numerical resolution of 256² grid cells and the same initial dust-to-gas ratio ${\varepsilon }_{0}=3$. The physical domain size L is altered around the predicted critical length scale $2{l}_{{\rm{c}}}$. Two simulations around $2{l}_{{\rm{c}}}$ with St = 0.1 are additionally altered in gas pressure gradient η. Maximum dust density time series can be found in Figure 6. The z-dimension has only one grid cell. Table 2 gives the measured turbulence and diffusivity, including the predicted length scale ${l}_{{\rm{c}}}$.

Zoom In Zoom Out Reset image size

Table 1. Overview of All Simulations

Name	Lx, Ly	dx,y	η	${T}_{\max }\left[{T}_{\mathrm{orb}}\right]$
Ae3L002	0.02 H	7.81 × 10−5	0.05	2.82
Ae3L001	0.01 H	3.91 × 10−5	0.05	3.67
Ae3L0005	0.005 H	1.95 × 10−5	0.05	4.24
Ae3L0005lp	0.005 H	1.95 × 10−5	0.025	13.06
Ae3L0005hp	0.005 H	1.95 × 10−5	0.1	32.78
Ae3L0003	0.003 H	1.17 × 10−5	0.05	10.03
Ae3L0003lp	0.003 H	1.17 × 10−5	0.025	14.47
Ae3L0002	0.002 H	7.81 × 10−6	0.05	3.83
Ae3L0001	0.001 H	3.91 × 10−6	0.05	3.50
Be3L005	0.05 H	1.95 × 10−4	0.05	12.57
Be3L003	0.03 H	1.17 × 10−4	0.05	26.22
Be3L002	0.02 H	7.81 × 10−5	0.05	50.93
Be3L001	0.01 H	3.91 × 10−5	0.05	31.83
Be3L0005	0.005 H	1.95 × 10−5	0.05	16.84
Be3L0003	0.003 H	1.17 × 10−5	0.05	11.58

Note. The name indicates the Stokes number: A—St = 0.1; B—St = 0.01. The rest of the name refers to the initial dust-to-gas ratio (always the same) and to domain size L in the following column, grid spacing d_x,y, gas sub-Keplerianicity η, and maximum simulation run time in orbits. Self-gravity is turned on at T = 1.59T_orb in the A runs and at T = 4.77 ${T}_{\mathrm{orb}}$ for the B runs. Additional simulations were performed with variation in the pressure gradient η by a factor of 2 (hp = high pressure) or by a factor of $\tfrac{1}{2}$ (lp = low pressure).

Download table as:  ASCIITypeset image

Table 2. Simulation Results: A Runs with St = 0.1 Particles and B Runs with St = 0.01 Particles

Name	Np.	δx	Δδx	${l}_{{\rm{c}}}$	urms	urms,x	vrms	vrms,x
Ae3L002	2	1.23 × 10−5	4.80 × 10−8	7.40 × 10−3	5.77 × 10−3	6.93 × 10−3	4.02 × 10−3	5.02 × 10−3
Ae3L001	8	8.34 × 10−6	6.94 × 10−8	6.09 × 10−3	4.24 × 10−3	6.78 × 10−3	3.05 × 10−3	4.88 × 10−3
Ae3L0005	1	5.86 × 10−6	2.95 × 10−8	5.10 × 10−3	3.55 × 10−3	4.39 × 10−3	2.74 × 10−3	3.22 × 10−3
Ae3L0005lp	3	2.25 × 10−6	1.19 × 10−6	3.16 × 10−3	2.27 × 10−3	2.60 × 10−3	1.79 × 10−3	2.00 × 10−3
Ae3L0005hp	0	1.48 × 10−5	9.09 × 10−6	8.12 × 10−3	6.20 × 10−3	8.08 × 10−3	4.76 × 10−3	6.44 × 10−3
Ae3L0003	0	2.26 × 10−6	1.36 × 10−8	3.17 × 10−3	2.55 × 10−3	4.20 × 10−3	1.71 × 10−3	2.96 × 10−3
Ae3L0003lp	1	1.04 × 10−6	5.61 × 10−7	2.15 × 10−3	1.59 × 10−3	2.52 × 10−3	1.23 × 10−3	2.07 × 10−3
Ae3L0002	0	2.00 × 10−6	1.58 × 10−8	2.98 × 10−3	1.62 × 10−3	2.61 × 10−3	1.34 × 10−3	1.82 × 10−3
Ae3L0001	0	1.31 × 10−6	8.04 × 10−9	2.41 × 10−3	1.79 × 10−3	4.56 × 10−3	0.88 × 10−3	2.84 × 10−3
Be3L005	2	2.36 × 10−5	8.25 × 10−6	3.24 × 10−2	5.68 × 10−3	6.81 × 10−3	5.53 × 10−3	6.61 × 10−3
Be3L003	1	1.81 × 10−5	9.23 × 10−6	2.84 × 10−2	4.27 × 10−3	5.05 × 10−3	4.07 × 10−3	4.83 × 10−3
Be3L002	0	1.28 × 10−5	7.40 × 10−6	2.39 × 10−2	3.97 × 10−3	4.91 × 10−3	3.75 × 10−3	4.67 × 10−3
Be3L001	0	5.09 × 10−6	1.44 × 10−6	1.50 × 10−2	3.44 × 10−3	3.54 × 10−3	3.22 × 10−3	3.32 × 10−3
Be3L0005	0	2.85 × 10−6	9.35 × 10−7	1.13 × 10−2	3.54 × 10−3	3.09 × 10−3	2.40 × 10−3	2.88 × 10−3
Be3L0003	0	1.49 × 10−6	6.58 × 10−7	8.13 × 10−3	3.36 × 10−3	2.77 × 10−3	2.56 × 10−3	2.46 × 10−3

Note. Diffusivities and velocities are measured by tracking the radial position of 10⁴ particles for several orbits and treating it similarly to a turbulence-driven random walk. Diffusivity and rms velocities are measured in the nongravitating fully turbulent situation. The number of planetesimals N_p is the number of objects we find in our final snapshots.

Download table as:  ASCIITypeset image

Different resolution and different sizes of the box change the strength of the streaming instability and thus the turbulent diffusion (Schreiber & Klahr 2018). Therefore, each setup has a different critical l_c , even for the same pressure gradient, dust-to-gas ratio, and Stokes number. As a consequence, we do not determine here the ultimate value for δ for streaming instability in general, for which one would need global 3D high-resolution studies with a wide range of Stokes numbers, but we focus on testing the validity of the 2l_c  < L criterion by varying the box size L. An alternative method is to keep L fixed but to alter the pressure gradient, which we have shown by means of additional simulations (Ae3L0005lp, Ae3L0005hp, Ae3L0003lp).

In comparison to other work (Johansen et al. 2015; Simon et al. 2016), we were able to ensure to resolve the critical length scale ${l}_{{\rm{c}}}$ by 64 grid cells or more. We varied the size of the simulation domain, L, in a set of models at Hill density, and following our prediction, only boxes larger than $2{l}_{{\rm{c}}}$ collapsed (Figure 7), i.e., when a cloud of radius ${l}_{{\rm{c}}}$ would have fit into the box. The radial diffusivity δ_x is measured in the situation of saturated streaming instability, but before gravity is switched on, see additional plots in Appendix E. In this measurement, the diffusivity increases with simulation domain size L, since larger modes of the streaming instability are stronger diffusing particles, which are suppressed in smaller simulation domain sizes. As shown in this paper, the simulations Ae3L0003, Ae3L0002, and Ae3L0001 are the ones not collapsing from our A parameter set. Even Ae3L0003 is not collapsing after more than eight orbits (see Figure 8). On this scale, diffusion acts faster than collapse, whereas on scales larger than L ≥ 0.005H planetesimals did form.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As predicted, the run Ae3L0005 is the smallest simulation still capable of producing a planetesimal, and in fact there is only one forming (see Figure 8). For the next largest simulation Ae3L001 we could count eight bound objects of different appearing size, varying by a only a factor of 2 in size. Two of them are in a bound binary system; see Movie 2 online. The Ae3L002 simulation produces also the formation of several planetesimals, which start colliding and merging; thus, only two planetesimals survive. This collision and merging has to be taken with a grain of salt, as we do not allow our pebble clouds to contract to solid density. We refer to dedicated simulations of cloud collapse as performed by Nesvorný et al. (2010).

All rms velocities are measured in a nongravitating fully SI turbulent snapshot. Gas rms velocities are measured by using the grid data, particle rms velocities by using the complete particle data set. Diffusivities are calculated by tracking a set of 10⁴ particles (see Appendix E.6). The measured values are summarized in Table 2.

From the maximum dust density time series for the models using St = 0.1 (Figure 6) one finds that with decreasing box size it takes longer to form a planetesimal (blue lines); even so, diffusivity was getting weaker in those runs. When crossing the border of stability $2{l}_{{\rm{c}}}$ planetesimal formation stalls and the particle cloud remains in a turbulent state. When we decreased (or increased) the radial pressure gradient (see Appendix F) a weaker (stronger) streaming instability led to a smaller (larger) critical length l_c. Thus, we find that our criterion reliably predicts the outcome of our numerical experiments for different Stokes numbers and pressure gradients. Please see also our online content: Movie 1 (https://youtu.be/gkHiluqH8HY) compares simulations Ae3L0005 and Ae3L0005. Both use St = 0.1 particles, but only the larger box shows collapse and planetesimal formation. In Movie 2 (https://youtu.be/nA87-9_trUc) we show the evolution of St = 0.1 pebbles for all six different box sizes in Table 1, and in Movie 3 (https://youtu.be/CCywDPKVU8w) we show the same for the 10 times smaller particles with St = 0.01. One clearly sees that for the smaller pebbles the contraction into planetesimals takes longer, as expected. But most importantly, as can be read from Table B1, our collapse criterion always gave the right prediction on whether a simulation would lead to collapse or not.

Zoom In Zoom Out Reset image size

As also can be seen from Table 2 and Figure 7, each simulation found a different critical length scale l_c range for the cases that lead to collapse. For St = 0.1 particles it was l_c = (7.4−5.1) × 10⁻³, whereas smaller Stokes numbers St = 0.01 had l_c = (3.2−2.8) × 10⁻². The first trend is that smaller Stokes numbers lead to larger planetesimals, as expected, yet also the diffusivity apparently changed; thus, the St = 0.01 unstable cloud was not 10 but only about 5 times larger than the average St = 0.1 pebble cloud. This has to be taken with a grain of salt; note that we only did 2D simulations of a pretty restricted local size to test the l_c criterion. In more global simulations one finds stronger diffusion and interestingly in some cases even a scaling of δ ∼ St, which leads to an l_c independent of an explicit δ and St (Schreiber & Klahr 2018). More investigations on δ as a function of St in more global setups are desperately needed to constrain l_c directly on (a range of) particle sizes, radial pressure gradient, and local dust load in the disk. Still our result holds: if you run a turbulent particle and gas simulation for a certain Stokes number and determine the diffusivity before turning on self-gravity, then our L > 2l_c criterion can tell you whether collapse and planetesimal formation will occur.

It makes sense to ask how much mass is in the gravitational pebble cloud of radius l_c . This characteristic mass is the initial condition for the further collapse into one or multiple planetesimals. As for a given diffusivity δ, Stokes number St, and pressure scale height H/R, l_c scales linearly with distance R to the star, and the volume of the pebble cloud scales as R³. And at the same time the Hill density drops as R⁻³, which indicates that the mass m_c included in our pebble cloud of size l_c is first-order independent from the distance to the star,

$$\begin{eqnarray}&&{m}_{c}=\displaystyle \frac{4\pi }{3}{l}_{{\rm{c}}}^{3}{\rho }_{\mathrm{Hill}}=\displaystyle \frac{1}{9}{\left(\displaystyle \frac{\delta }{\mathrm{St}}\right)}^{\tfrac{3}{2}}{\left(\displaystyle \frac{H}{R}\right)}^{3}{M}_{\odot },\end{eqnarray} \tag{ 10 }$$

but only a function of pebble size, turbulence, aspect ratio of the disk, and stellar mass. Note that there are no additional dependencies on the actual metallicity of the solar nebula (respectively protoplanetary disk), nor on its mass or density profile.

When we now calculate a diameter a_eq for a solid body of equivalent mass as this cloud of mass m_c , this is not to claim that the cloud will collapse at 100% efficiency into precisely one planetesimal; it rather shall express what size ranges are possible in the collapse and potential subsequent fragmentation. To accommodate for this uncertainty, we can incorporate the collapse efficiency q, which describes what fraction of the pebble cloud ends up in one individual planetesimal.

To convert mass into a size, one needs a density, which makes a difference whether you form some fluffy comet with ${\rho }_{\bullet }\,=0.5\,{\rm{g}}/{\mathrm{cm}}^{3}$ or an asteroid with a mean density of ${\rho }_{\bullet }=2\,{\rm{g}}/{\mathrm{cm}}^{3}$. For convenience we therefore use the density of our Sun as density ${\rho }_{\bullet }=1\,{\rm{g}}/{\mathrm{cm}}^{3}$, which makes our estimate much simpler. The error in size we introduce is thus ±20%, which is currently beyond the precision of our theory anyway. This is also consistent with the conversion of mass into equivalent size as done in the main part of Nesvorný et al. (2010). ² Thus, our equivalent body will have a diameter of

$$\begin{eqnarray}&&{a}_{\mathrm{eq}}=2{l}_{{\rm{c}}}{\left(\displaystyle \frac{{\rho }_{\mathrm{Hill}}}{{\rho }_{\bullet }}\right)}^{1/3},\end{eqnarray} \tag{ 11 }$$

and the possible forming planetesimals will be slightly smaller, as q enters the size only weakly:

$$\begin{eqnarray}&&{a}_{{\rm{c}}}={q}^{1/3}{a}_{\mathrm{eq}}.\end{eqnarray} \tag{ 12 }$$

Combining Equations (6) and (7) and expressing a_eq in terms of solar radii, one arrives at the relation

$$\begin{eqnarray}&&{a}_{\mathrm{eq}}=\sqrt[3]{\displaystyle \frac{8{\rho }_{\odot }\,}{9{\rho }_{\bullet }}}\sqrt{\displaystyle \frac{\delta }{\mathrm{St}}}\displaystyle \frac{H}{R}{R}_{\odot },\end{eqnarray} \tag{ 13 }$$

where the first term is of order unity (mean density of the Sun being similar to the mean density of planetesimals) and could be neglected for order-of-magnitude estimates.

In Lenz et al. (2020) we constrain the parameter space for the solar nebula via the influence of disk properties on planetesimal formation in our paradigm of pebble flux regulated planetesimal formation. The idea is similar to the MMSN (Hayashi 1981), but in contrast to that model, dust does not locally grow into planetary cores, but pebbles drift large distances, before being converted into planetesimals. Therefore, the nebulae in this paradigm can be shallower in surface density profile of gas and still form centrally concentrated distributions of planetesimals. At the same time, this shallower profile is consistent with viscous accretion disk theory for constant α values and also fits better observations of disks around young stars (Andrews et al. 2010).

The best parameter set to produce a planetesimal population that could explain the formation of the solar system (see Lenz et al. 2020 for details) is a relatively massive yet gravitational stable disk mass of 0.1M, a viscosity of α = 3 × 10⁻⁴, an exponential cutoff radius at 20 au, and a fragmentation speed of v_frag = 200 cm s⁻¹ that determines the Stokes numbers (see Figure 4) of the largest pebbles for the assumed background turbulence α. This leads to the gas density distribution as seen in Figure 2. The u-shape in the ε_Hill is a direct result from the truncation radius of 20 au.

From Schreiber & Klahr (2018) we know that the diffusivity δ₀ = 2.7 × 10⁻⁶ for St = 0.1 and ε = 10 roughly scales inversely with ε in the range of interest 10–100 and about linear with the Stokes number; thus, we use the prescription

$$\begin{eqnarray}&&\delta ={\delta }_{0}\displaystyle \frac{10}{{\varepsilon }_{\mathrm{Hill}}}\displaystyle \frac{\mathrm{St}}{0.1}.\end{eqnarray} \tag{ 14 }$$

The locally dominating Stokes number we can estimate from the fragmentation limit (Birnstiel et al. 2012). The particle size is determined by global turbulence (α), where pebbles spend most of their time, before being locally concentrated to Hill density; therefore, ε does not affect St:

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{frag}}=\displaystyle \frac{1}{3}\displaystyle \frac{{v}_{\mathrm{frag}}^{2}}{\alpha {c}_{s}^{2}}.\end{eqnarray} \tag{ 15 }$$

We plot the Stokes numbers in Figure 4. Thus, as we are in the range of validity for Equation (14), we see that δ/St simplifies to δ₀ 10/ε_Hill. Now we receive the length scale prediction:

$$\begin{eqnarray}&&{l}_{c}\propto \sqrt{{\delta }_{0}}\sqrt{\displaystyle \frac{10}{{\varepsilon }_{\mathrm{Hill}}}}H.\end{eqnarray} \tag{ 16 }$$

From this we can calculate the mass of the unstable cloud with Equation (10) and from that again the equivalent diameter a_eq if that mass is compressed into a solid body of roughly density ${\rho }_{\bullet }=1\,{\rm{g}}/{\mathrm{cm}}^{3}$ (Equation (13)). For instance, for the MMSN with Σ ∝ R^−1.5 and H/R = 0.025 · (R/au)^1/4 this leads to a_eq ∝ R^−3/8, explaining the worst-case size difference between, for instance, 3 and 30 au by a factor of 2.3.

In Figure 10 we plot the resulting equivalent diameters. The predicted sizes are not constant with radius, yet vary surprisingly little from 80 to 140 km for a region from 3 to 50 au. We also show predicted sizes when we adopt the MMSN and a 3 times MMSN, and the results are also in the 30−100 km range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The size range of around 80 km for the asteroid belt fits nicely with the measurements by Delbo' et al. (2017). We also find sizes on the order of 100 km in the Kuiper Belt region, which opens the question of how to form smaller objects. Thus, we also calculate the critical sizes for a nebula that has dropped with viscous evolution to 10% and 1% of its initial mass, and we see that the equivalent size will also shrink over time. Thus, at the current location of, for instance, Arrokoth of 40 au, the equivalent size will shrink from 50 km down to 13 km, which would argue for a late formation of Arrokoth (Stern et al. 2019) with its equivalent radius of about 20 km. Of course, depending on the collapse efficiency and number of multiples that formed, the birth cloud of Arrokoth might also have had a larger mass.

The presented simulations in our paper were 2D. Meanwhile, we tested our criterion of stability $2{l}_{{\rm{c}}}\lt L$ also in a limited set of 3D simulations (Klahr & Schreiber, AAS25612) and find a conformation of our findings from the present paper.

Yet eventually one has to test our paradigm for a range of dust-to-gas ratios, pebble sizes, and total mass of the disk (Gerbig et al. 2020), which all have the effect of the critical mass triggering collapse. Such a detailed study does not exist yet. Studies like Schäfer et al. (2017), Simon et al. (2017), Abod et al. (2019), and several more use, in fact, very large boxes and form hundreds of planetesimals to study their size distribution as a function of the nebula conditions. But unfortunately, diffusivity was not measured in these simulations to check for the applicability of our criterion. Vertical diffusion could be measured in post processing for those existing simulations by measuring the dust scale height in the turbulent state with and without self-gravity (Johansen et al. 2007b). This method is unfortunately not possible for radial diffusion, as there is no equilibrium state with gravity balanced by diffusion. If radial and vertical diffusion were equal, one could rely on the vertical diffusion to estimate l_c , yet turbulence from streaming and Kelvin–Helmholtz instability is known to be rather anisotropic (Johansen & Youdin 2007; Schreiber & Klahr 2018; Gerbig et al. 2020). So eventually one has to repeat those simulations on the size distribution of planetesimals formed via self-gravity and streaming instability (Schäfer et al. 2017; Simon et al. 2017; Abod et al. 2019) and then apply a particle tracker as we describe in Appendix E.6 to determine radial diffusion. In Klahr & Schreiber (AAS25612) we present such a study on 3D streaming instability and the measurement of radial and vertical diffusion, yet on much smaller scales than in the aforementioned studies.

If the pebble cloud were to collapse into a single object, we could directly use the equivalent size as planetesimal size. This is, of course, not likely, as angular momentum conservation in the spinning and collapsing cloud can lead to fragmentation into multiple planetesimals, just like it does in star formation, and also the headwind may drain some material from the collapsing cloud. Yet If only q = 1/8 of the collapsing cloud is converted into an individual planetesimal, its size will be smaller by just a factor of 2 in comparison to our q = 1 estimate. This effect alone may be sufficient to explain a size spread among asteroids.

5.1. The Effect of Global Turbulence

If we use global turbulence in the solar nebula instead of the locally generated turbulence by the streaming instability, we can also estimate a critical size and mass scale for gravitational collapse. First, we have to consider that the global turbulence will have to cascade down from the integral or driving length scale of turbulence to the scales of l_c . Assuming a Kolmogorov cascade, as the most optimistic model, we derived the diffusivity $\delta ({l}_{{\rm{c}}})$ for the case in which the global turbulence would mix proportional to the traditional α value, acting on the integral scale of turbulence $L=\sqrt{\alpha }H$ (Schreiber & Klahr 2018):

$$\begin{eqnarray}&&\delta ({l}_{{\rm{c}}})=\displaystyle \frac{{\alpha }^{\tfrac{1}{3}}}{1+{\varepsilon }_{\mathrm{Hill}}}{\left(\displaystyle \frac{{l}_{{\rm{c}}}}{H}\right)}^{\displaystyle \frac{4}{3}},\end{eqnarray} \tag{ 17 }$$

considering the mass load of material (ε_Hill > 1) in the dust clump reaching Hill density (Johansen et al. 2007a). As can be seen in Figure 2, the Hill density has a different slope than the gas density in the midplane, and the latter will also change over time. So depending on the mass and profile of the solar nebula, ε_Hill may vary between 200 at 1 au and 20 at 30 au. The same may be true for α, which may be quite different at 1 au versus 100 au.

If we use this local diffusivity stemming from global turbulence in the estimate for the critical length scale Equation (6), we find

$$\begin{eqnarray}&&{l}_{{\rm{c}}}=\displaystyle \frac{1}{27}\displaystyle \frac{{\alpha }^{1/2}}{{\left(1+\varepsilon \right)}^{\tfrac{3}{2}}}{\mathrm{St}}^{-\tfrac{3}{2}}H.\end{eqnarray} \tag{ 18 }$$

Note the strong dependency on the Stokes number in this case, which comes from the effect that the strength of diffusion got a length scale dependence. In the case for particle-induced turbulence, when diffusivity depends on the Stokes number, St canceled out from the equations. Yet here we have to use an explicit St, as it follows from the fragmentation limit (see Equation (15)):

$$\begin{eqnarray}&&{l}_{{\rm{c}}}=\frac{1}{\sqrt{27}}\frac{{\alpha }^{2}}{{\left(1+\varepsilon \right)}^{\frac{3}{2}}}\frac{{c}_{s}^{3}}{{v}_{\mathrm{frag}}^{3}}H.\end{eqnarray} \tag{ 19 }$$

For our α = 3 × 10⁻⁴, which we also used in the previous part, we receive equivalent sizes that are generally smaller than the ones for the pure streaming case (see Figure 10), and thus global turbulence does not play a role here to set the smallest scales. Yet note that the influence of global turbulence depends much more strongly on α; thus, already a global alpha of 1 × 10⁻³ as used in Hartlep & Cuzzi (2020) would lead to 10 times larger sizes. Also the v_frag has a dramatic effect in this case. That alone should rule out pure external turbulence as setting the length scales for planetesimal formation in our solar system. Detailed 3D simulations of this scenario are, of course, still missing, especially as the source and related the overall strength of turbulence in the solar nebula and the shape and the extent of the turbulent spectrum are heavily under debate (Klahr et al. 2018; Pfeil & Klahr 2019).

In the clustering model (Hartlep & Cuzzi 2020) the authors use α = 0.001 and a Stokes number of St = 0.04 to form planetesimals at 3 au in the solar nebula, indicating a large fragmentation velocity of about 600 cm s⁻¹. In one of the models they assume a 10-fold MMSN, leading to a lower necessary concentration for collapse of ε_Hill = 25 in agreement with our estimate (see Figure 2). In that case we get a threshold size for the pebble cloud equivalent to a planetesimal diameter of 2300 km, which is much larger than their derived lower threshold based on ram pressure as argued for in that paper of about 10 km. But in order to have such a small ram pressure threshold, they had to assume to be in a zonal flow where the headwind was reduced by a factor of 30. If the headwind was not decreased, then they would also have received a threshold of more than 1000 km, because, as can be seen in their Figure 9, the lower limit due to ram pressure (blue curve) would move upward by a factor of 30, and then the intersection with the upper limit for mass loading (red curve) would fall at a size of more than 2000 km.

So whereas it is justified to reduce the headwind in a zonal flow, the effect of turbulent diffusion will not be reduced, as it is an integral part of turbulent clustering. The argument that turbulent diffusion can be neglected in comparison to ram pressure (Cuzzi et al. 2008) does not hold, if one is reducing the ram pressure.

Streaming instability in our simulations has a dual role in the process of planetesimal formation, both being contrary to each other. On larger scales the streaming instability helps to form planetesimals by concentrating dust into dense clouds and to reach Hill density, yet on small scales it prevents the formation of arbitrarily small planetesimals by diffusing collapsing clumps faster than they can collapse.

We derived a critical pebble cloud mass to undergo gravitational collapse in the presence of turbulent diffusion, which may be driven either by streaming instability or by global gas turbulence. We showed the validity of our criterion in 2D streaming instability simulations of planetesimal formation. The resulting critical pebble cloud masses for turbulence values typical for streaming instability correspond to equivalent diameters of a_eq ≈ 100 km and are thus compatible to those from individual planetesimals of up to 100 km in diameter or several smaller ones, depending on the efficiency of the final contraction and subsequent fragmentation into multiple planetesimal systems.

Global turbulence on low levels, as needed for streaming instability in the first place, seems not to have a strong impact on small scales, beyond setting the Stokes number St. But as long as St is in a range to trigger the streaming instability, the ratio of diffusivity over Stokes number is approximately constant as far as we know, and thus neither α nor the value for the fragmentation speed v_frag has an influence on the equivalent diameter a_eq.

Under the assumptions that streaming instability leads to diffusion inversely proportional to the Stokes number and proportional to the dust-to-gas ratio upon reaching Hill density, as found in Schreiber & Klahr (2018), the equivalent diameter depends only on the local scale height ratio ∝H/R and the inverse square root of the local dust-to-gas ratio of pebbles at Hill density $\propto {\varepsilon }_{\mathrm{Hill}}^{-1/2}$. As H/R and ε_Hill in a gas model of the solar nebula increase slowly with distance to the Sun, the value for the equivalent diameter a_eq may vary by only a factor of 2 between 3 and 30 au. Considering a steeper gas profile in the outer part of the nebula can become even smaller than at 3 au.

Global turbulence alone as setting the size of planetesimals will introduce a huge error bar, as the equivalent diameter would strongly depend on fragmentation speed v_frag and α; thus, we can ignore this effect as long as the equivalent sizes are smaller than those set by streaming instability.

The derived criterion supports the idea that there is a preferred initial planetesimal birth size possibly with a Gaussian distribution and not yet the power-law distributions observed today. The power-law size distribution is then the outcome of planetesimal collisions and pebble accretion (Johansen et al. 2015). In fact, there is recent observational evidence that the initial size distribution of asteroids was much shallower than presumed (Tsirvoulis et al. 2018), and this could be reproduced indeed by a Gaussian initial distribution with a width of 45 km centered around a diameter of 80–85 km (Delbo' et al. 2017). The herein-explained diffusion-regulated gravitational collapse of a pebble cloud is so far the only prediction for a narrow initial size distribution of planetesimals instead of a wider power-law distribution.

Planetesimals of significantly smaller size (<10 km) should form with a lower likelihood in this process, because turbulence can destroy initial pebble clouds of an equivalently low mass before they collapse, and thus they can only be by-products of bigger planetesimals forming. And much larger planetesimals are less frequent, because during the slowly increasing local accumulation of pebbles the lowest possible mass will already lead to collapse.

We have shown that for a range of assumptions for the initial solar nebula this leads to equivalent radii of 80–140 km, for the regions of interest, explaining why planetesimals and thus asteroids and classical KBOs have a kink in their distribution at about the same size. In this paradigm, larger asteroids, as well as giant planet cores, are the result of secondary growth processes like pebble accretion (Klahr & Bodenheimer 2006; Ormel & Klahr 2010; Johansen et al. 2015) and collisions (Kobayashi et al. 2016), whereas smaller objects are either the outcome of a collisional fragmentation cascade (Morbidelli et al. 2009) or products of the pebble cloud fragmenting in a size range of objects. These evolution processes explain the currently observed power laws above and below the initial size, and thus the characteristic shape of the mass distribution of minor objects in the solar system of today is an imprint of the initial size that we explain in this work.

As the gas mass in the solar nebula decreases over time, the dust-to-gas ratio upon reaching Hill density for the pebbles will increase, which will lead to smaller planetesimal sizes. So in general the trend will be to first form large planetesimals and then later allow for smaller ones. Possibly at very late times with little gas left in the nebula the streaming instability was weak enough to allow for the formation of smaller planetesimals, which may be the origin of 1–10 km sized comets. Their mass contribution should then be lower when compared to the bigger planetesimals, as at late times also the pebble reservoir of the disk runs empty. Further research will have to clarify this scenario.

But keep in mind that we derived in this paper the equivalent to the Jeans mass in a cloud core before star formation. The gravitational collapse and the formation of multiple systems (Nesvorný et al. 2010) with an initial mass function below and above the critical pebble cloud mass are a whole different story and deserve more attention.

Online content includes three animations of the simulations, which are listed below, provided on YouTube, and preserved in Zenodo under a Creative Commons (Attribution) license: doi:10.5281/zenodo.3971104.

Movie 1: https://youtu.be/gkHiluqH8HY. Simulations Ae3L0005 and Ae3L0005. Both use St = 0.1 particles, but only the larger box shows collapse and planetesimal formation.

Movie 2: https://youtu.be/nA87-9_trUc. The evolution of St = 0.1 pebbles for all six different box sizes in Table 1. See Figure 11 for the end state of these simulations.

Movie 3: https://youtu.be/CCywDPKVU8w. The evolution of St = 0.01 pebbles.

We have to thank Bill Bottke for initiating this project by asking us how to explain 100 km sized objects all over the solar system. We are indebted to Anders Johansen, Marco Delbo, Allesandro Morbidelli, Jeff Cuzzi, Wladimir Lyra, Hans Baehr, Christian Lenz, and Karsten Dittrich for many fruitful discussions and technical advise. Many thanks also to David Nesvorny` and Andrew Youdin for discussing their estimates of collision time versus friction time with us. And last but not least, thanks to our anonymous referee, who suggested to include a detailed size estimate of objects in the solar nebula. A.S. has been supported by the Studienstiftung des deutschen Volkes. This research was funded by the Deutsche Forschungsgemeinschaft Schwerpunktprogramm (DFG SPP) 1385 "The First Ten Million Years of the Solar System" under contract KL 1469/4-(1-3), by DFG SPP 1833 "Building a Habitable Earth" under contract KL 1469/13-(1-2), by DFG SPP 1992: "Exoplanet Diversity" under contract KL 1469/17-1, by DFG Research Unit FOR 2544 "Blue Planets around Red Stars" under contract KL 1469/15-1, and by experimental work with our colleagues in Duisburg under contract KL 1469/14-1. We also acknowledge support from the DFG via the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany's Excellence Strategy (grant EXC-2181/1-390900948). We received additional support by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe." This research was performed in part at KITP Santa Barbara by the National Science Foundation under grant No. NSF PHY11-25915. The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time for a GCS Large-Scale Project (additional time through the John von Neumann Institute for Computing (NIC)) on the GCS share of the supercomputer JUQUEEN (Stephan & Docter 2015) at Jülich Supercomputing Centre (JSC). GCS is the alliance of the three national supercomputing centers HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich), and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK), and Nordrhein-Westfalen (MIWF). Additional simulations were performed on the THEO and ISAAC clusters of the MPIA and the COBRA, HYDRA, and DRACO clusters of the Max-Planck-Society, both hosted at the Max-Planck Computing and Data Facility in Garching (Germany).

Shi & Chiang (2013) discuss the criteria for a gravitational instability of small particles embedded in gas. They consider two cases: one in which the sound crossing time across a self-gravitating particle cloud is longer than the stopping time (τ_sound > ${\tau }_{{\rm{s}}}$) between particles and gas, and one case in which it is shorter. In the first case, the mixture behaves like a suspension and the clump gets stabilized by the pressure gradient of the gas as it is getting compressed. In this case one has to consider the stability of the gas and dust mixture (Cuzzi et al. 2008; Shi & Chiang 2013). In the other extreme, dubbed as the "small clump" regime, when the stopping time is longer than the sound crossing time (τ_sound < ${\tau }_{{\rm{s}}}$), one can neglect the effect of the gas being compressed. If we consider particles with a Stokes number of St = 0.1 and typical dust-to-gas ratios of ε = 3–100 for collapse, then the sound crossing distance is $\lambda =H\ \tfrac{\mathrm{St}}{\sqrt{\varepsilon }}\gt {10}^{-2}H$; see Equation (39) in Shi & Chiang (2013). This distance is larger than the clumps we consider in the main part of the paper for gravitational collapse l_c  ≈ 4 × 10⁻³ H, and the gas can be treated as incompressible in our considerations. In the case of the smaller particles (St = 0.01) we are already entering the "large clump" regime, which we discuss in Section E.8, as the gas is getting slightly compressed during the collapse of the particle cloud. Nevertheless, the collapse criteria that we derive for the "small clump" regime still hold.

Tidal forces and shear forces exerted by the host star are able to disrupt clumps if their density is less than a critical density ${\rho }_{{\rm{c}}}$. For an overview of all used symbols and quantities in this paper we refer to Table B1. For instance, a comet gets disrupted in close vicinity to a star or a planet, under the condition of differential gravity (tidal force) from the central object being stronger than the internal gravitational binding. This Roche criterion can be expressed via a two-sphere problem with radii a/2, mass m, and separation of a, both spheres being located at a mean distance R from the central object of mass M. Note here that the Roche criterion assumes neither orbital motion nor rotation, as well as no underlying gas flow. Thus, the maximal separation between the two spheres to overcome the tidal forces is

$$\begin{eqnarray}&&{a}_{{\rm{c}}}=R{\left(\displaystyle \frac{m}{16M}\right)}^{\displaystyle \frac{1}{3}}.\end{eqnarray} \tag{ B1 }$$

This derivation directly leads to a critical breakup density for a spherical particle cloud

$$\begin{eqnarray}&&{\rho }_{{\rm{c}}}=2.5\displaystyle \frac{M}{{R}^{3}},\end{eqnarray} \tag{ B2 }$$

whereas the more detailed work of Chandrasekhar (1967) gives a value for the Roche density of

$$\begin{eqnarray}&&{\rho }_{\mathrm{Roche}}=3.5\displaystyle \frac{M}{{R}^{3}}.\end{eqnarray} \tag{ B3 }$$

The Roche criterion is useful to study the breakup of bodies in a close encounter but less suited for the stability analysis of a self-gravitating particle cloud that is in an orbital motion. For this situation it is necessary to include a centrifugal potential around the primary object. This is the Hill criterion, in which a test particle stays bound to a secondary orbiting object; in our case it is the center of mass of the particle cloud, with distance a between test particle and center of mass. Assuming a spherical cloud of homogeneous density distribution with a total mass of m rotating at distance R around a central star with mass M, this leads to an expression equivalent to the Hill sphere with radius

$$\begin{eqnarray}&&{a}_{\mathrm{Hill}}=R\sqrt[3]{\displaystyle \frac{m}{3M}}.\end{eqnarray} \tag{ B4 }$$

A test particle can only stay bound if its distance a from the center of mass m is a < a_Hill. Based on that, one can derive a critical density of that particle cloud, the Hill density:

$$\begin{eqnarray}&&{\rho }_{\mathrm{Hill}}=\displaystyle \frac{9}{4\pi }\displaystyle \frac{M}{{R}^{3}}\approx 0.72\displaystyle \frac{M}{{R}^{3}}.\end{eqnarray} \tag{ B5 }$$

This value is smaller than the Roche density by a factor of 5, because the bound rotation of the whole particle cloud around its host star gives an additional stabilizing effect to it.

Table B1. Used Symbols and Quantities

Symbol	Definition	Description
a	2a = d	pebble radius and diameter
${\tau }_{{\rm{s}}}$	${\tau }_{{\rm{s}}}=\tfrac{a{\rho }_{\bullet }\pi }{2{\rm{\Omega }}{{\rm{\Sigma }}}_{\mathrm{gas}}}$	friction time/stopping time
$\mathrm{St}$	${\tau }_{{\rm{s}}}{\rm{\Omega }}$	Stokes number
${l}_{{\rm{c}}}$	${l}_{{\rm{c}}}={r}_{\mathrm{crit}}$	critical length scale (radius) of a particle cloud
mc	${m}_{c}=\tfrac{4\pi }{3}{l}_{{\rm{c}}}^{3}{\rho }_{\mathrm{Hill}}$	critical mass of a particle cloud
aeq	${a}_{\mathrm{eq}}=2{l}_{{\rm{c}}}{\left(\tfrac{{\rho }_{\mathrm{Hill}}}{{\rho }_{\bullet }}\right)}^{1/3}$	equivalent contracted diameter of a particle cloud
u, v		gas and dust velocity
vfrag		dust fragmentation velocity
Ω		orbital frequency
R,z		semimajor axis, vertical distance to midplane
${c}_{{\rm{s}}}$		sound speed
H,hp		gas and particle scale height
Torb	Torb = 2π/Ω	orbital period
H	$H={c}_{{\rm{s}}}/{\rm{\Omega }}$	gas disk scale height
R		heliocentric distance
${R}_{\odot }$, ρSun		stellar radius and solar density
${\rho }_{\mathrm{Hill}}$	${\rho }_{\mathrm{Hill}}=9M/4\pi {R}^{3}$	Hill density
${\rho }_{{\rm{c}}}$		shear stable could density, expressed in Hill density
${\rho }_{\bullet }$		solid-body density
${\rho }_{{\rm{d}},0}$		initial mean dust density (simulation)
${\rho }_{{\rm{g}},0}$		initial mean gas density (simulation)
ε	$\varepsilon ={\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}$	dust-to-gas density ratio
${\varepsilon }_{\max }$, ${\varepsilon }_{0}$		maximum and initial dust-to-gas ratio (simulation)
εHill	${\varepsilon }_{\mathrm{Hill}}={\rho }_{\mathrm{Hill}}/{\rho }_{{\rm{g}}}$	dust-to-gas ratio to reach Hill density
D, δ	$\delta =D/{{Hc}}_{{\rm{s}}}$	(dimensionless) local diffusion coefficient
α		(dimensionless) large-scale viscosity and diffusion coefficient
${\tau }_{\mathrm{ff}}$		free-fall time
${\tau }_{{\rm{c}}}$		contraction time (incl. friction)
${\tau }_{{\rm{D}}}$	${\tau }_{{\rm{D}}}={r}^{2}/D$	diffusion time
${\lambda }_{c}$	${\lambda }_{c}=2\pi {l}_{{\rm{c}}}$	Jeans length scale for our simulations
r		particle cloud radius
${\lambda }_{\mathrm{free}}$		mean free path
L	[1/H]	simulation domain size
dx,y	[1/H]	simulation grid resolution
η	$\eta =\tfrac{1}{2}{\left(\tfrac{H}{R}\right)}^{2}\tfrac{\mathrm{dln}\rho }{\mathrm{dln}R}$	pressure gradient parameter
$\hat{G}$		self-gravity parameter (simulation)

Download table as:  ASCIITypeset image

The Roche density is derived for gravity only; thus, the gravitational acceleration difference by the Sun δg across a body of diameter a at location R scales as $\delta g=-{{\rm{\Omega }}}^{2}{R}^{2}\left(\tfrac{1}{{R}_{o}^{2}}-\tfrac{1}{{R}_{i}^{2}}\right)$ with ${R}_{o,i}=R\pm \tfrac{1}{2}a$, whereas the Hill density is calculated for an object in circular orbit. Then, the effective potential due to rotation reduces the difference in radial acceleration to δg* =δg + Ω² (R_o  − R_i ) = δg + Ω² a; thus, a lower density is sufficient to prevent the tidal disruption of a body.

Sekiya (1983) defines his critical density for an axisymmetric 3D annulus of particles and finds a value of

$$\begin{eqnarray}&&{\rho }_{\mathrm{Sekiya}}=0.62\displaystyle \frac{M}{{R}^{3}},\end{eqnarray} \tag{ B6 }$$

in fact very close to our simple derivation.

As a conclusion of our derivation for a diffusion-limited collapse, we find that only the sixth root of the particle cloud density enters the resulting planetesimal size. Hence, all above-mentioned estimates will give nearly the same result, which is sufficient for an order-of-magnitude estimate. In the following we will use the ${\rho }_{\mathrm{Hill}}$ in our calculations as critical density. In several works on planetesimal formation in fact the Hill density is used, yet it is often labeled as Roche density (e.g., Johansen et al. 2014). So we wanted to elude a little the subtle differences in the definition.

As our gas is effectively incompressible during the collapse of the pebble cloud of the St = 0.1–0.01 particles, the dust–gas friction indeed alters the contraction time to longer times than the free-fall time. Here we derive a contraction timescale ${\tau }_{{\rm{c}}}$ with dependency on Stokes number St = ${\tau }_{{\rm{s}}}$Ω, ignoring pressure effects and assuming collapse at terminal velocity.

C.1. Contraction Time for Frictional Particles

A pressure-free sphere of density ρ collapses under its own gravity within a free-fall timescale

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}=\sqrt{\displaystyle \frac{3\pi }{32G\rho }},\end{eqnarray} \tag{ C1 }$$

which for the case of Hill density $\rho ={\rho }_{\mathrm{Hill}}$ (see Equation (B5)) is 1/10 of an orbital period T_orb = 2πΩ⁻¹:

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}=0.64{{\rm{\Omega }}}^{-1}\approx 0.1{T}_{\mathrm{orb}}.\end{eqnarray} \tag{ C2 }$$

Yet, in the case of stopping time ${\tau }_{{\rm{s}}}$ by particle–gas friction being shorter than the free-fall time, particles can maximally fall at their terminal velocity (Cuzzi et al. 2008):

$$\begin{eqnarray}&&{v}_{{\rm{t}}}(r)=-{\tau }_{{\rm{s}}}\displaystyle \frac{{mG}}{{r}^{2}}.\end{eqnarray} \tag{ C3 }$$

We calculate this new frictional contraction time from Equation (C3) via integration:

$$\begin{eqnarray}&&r(t)=\sqrt[3]{{r}_{0}^{3}-3{\tau }_{{\rm{s}}}{mGt}}\quad r\left({\tau }_{{\rm{c}}}\right)=0\ \Rightarrow \ {\tau }_{{\rm{c}}}=\displaystyle \frac{{r}_{0}^{3}}{3{\tau }_{{\rm{s}}}{mG}}.\end{eqnarray} \tag{ C4 }$$

Hence, a clump of size r₀ and mass density ρ is expected to collapse within a collapse time

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=\displaystyle \frac{1}{4\pi {\tau }_{{\rm{s}}}\rho G}.\end{eqnarray} \tag{ C5 }$$

We can now express this collapse time in terms of free-fall time (Cuzzi et al. 2008) and combine both terms for long (Equation (C1)) and short stopping times (Equation (C5)) to

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}={\tau }_{\mathrm{ff}}\left(1+\displaystyle \frac{8{\tau }_{\mathrm{ff}}}{3{\pi }^{2}{\tau }_{{\rm{s}}}}\right).\end{eqnarray} \tag{ C6 }$$

In the case in which the Stokes number of the particles is smaller than the critical value of

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{crit}}=\displaystyle \frac{8}{3{\pi }^{2}}{\tau }_{\mathrm{ff}}=0.172,\end{eqnarray} \tag{ C7 }$$

the simple expression

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=\displaystyle \frac{1}{9\mathrm{St}}\,{{\rm{\Omega }}}^{-1}\end{eqnarray} \tag{ C8 }$$

is a sufficient approximation (see Figure C1).

Zoom In Zoom Out Reset image size

In the following section we compare our analytic estimates to a numerical integration of the settling process.

C.2. Numerical Test of Contraction Time and Analytic Fit

The differential equation governing the settling process is

$$\begin{eqnarray}&&{\partial }_{t}v=-\displaystyle \frac{{mG}}{{r}^{2}}-\displaystyle \frac{v}{{\tau }_{{\rm{s}}}}.\end{eqnarray} \tag{ C9 }$$

We time integrate this equation with a Leap Frog algorithm with the initial condition r(t = 0) = r₀ and v(t = 0) = 0. The parameters G and m are chosen in a way to initialize the cloud at Hill density, i.e., spreading the mass m evenly over the volume $V=\tfrac{4}{3}\pi {r}_{0}^{3}$. We performed a set of simulations for different single particle sizes ranging from St = 10⁻³ to St = 10; see Figure 1. We find that the simple fit from Equation (C6) is perfectly suited for all particle sizes, and even Equation (C8) gives good results up to St = 0.1.

Planetesimal formation happens in a shearing environment; hence, one needs to ensure tidal disruption of a particle cloud as the primary condition for planetesimal formation via gravitational collapse. As in the main paper, we start from setting diffusion time and collapse time equal. Here, lets consider small particles with St ≪ 1, falling at terminal velocity:

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{{\rm{c}}} & = & {\tau }_{{\rm{D}}}\\ \displaystyle \frac{{r}_{{\rm{c}}}^{3}}{3{\tau }_{{\rm{s}}}{mG}} & = & \displaystyle \frac{1}{4\pi {\tau }_{{\rm{s}}}{\rho }_{\mathrm{int}}G}=\displaystyle \frac{{r}_{{\rm{c}}}^{2}}{D}=\displaystyle \frac{{r}_{{\rm{c}}}^{2}}{\delta {{Hc}}_{{\rm{s}}}},\end{array}\end{eqnarray} \tag{ D1 }$$

where $m=4/3\pi {r}_{{\rm{c}}}^{3}{\rho }_{\mathrm{int}}$ is the bulk mass of the cloud with radius ${r}_{{\rm{c}}}$. The diffusion timescale stems from Fick's second law of diffusion, $\dot{\rho }=D{{\rm{\nabla }}}^{2}\rho$. We simplify by expressing the internal cloud density in terms of Hill density via a scaling parameter f:

$$\begin{eqnarray}&&{\rho }_{\mathrm{int}}=f\cdot {\rho }_{\mathrm{Hill}}=\displaystyle \frac{9f}{4\pi }\displaystyle \frac{M}{{R}^{3}}.\end{eqnarray} \tag{ D2 }$$

With this simplification, the critical cloud diameter is

$$\begin{eqnarray}&&{r}_{{\rm{c}}}={l}_{{\rm{c}}}=\displaystyle \frac{1}{3}\sqrt{\displaystyle \frac{\delta }{f\cdot \mathrm{St}}H}.\end{eqnarray} \tag{ D3 }$$

This expression is valid for all f as long as the condition for shear and tidal stability is given.

D.1. Jeans Length for Planetesimal Formation

A full analysis of the stability of dust under self-gravity embedded in gas would lead to a Toomre analysis. There one performs a linear analysis of the problem and derives a stability criterion from a dispersion relation, following a mixed case of Goldreich & Ward (1973) and Safronov (1969), similar yet not identical to the secular gravitational instability (Ward 2000). The resulting Toomre criterion would tell us whether there was a fastest growing mode, which, as we would see, is larger than our simulation domain. Toomre combines two obstacles for gravitational collapse: (A) tidal forces, i.e., angular momentum conservation on large scales, and (B) thermal pressure on small scales.

Thus, we will focus on the small scales, which is just the Jeans length, and we will see that even the Jeans length is larger than our box size.

The difference in the derivation here is that instead of the thermal pressure, we use the diffusion flux of particles j = −D∇ρ to be in equilibrium with sedimentation. Starting with the continuity equation for the dust particles,

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}v\rho =0.\end{eqnarray} \tag{ D4 }$$

The flux ρv is given by diffusion and sedimentation under self-gravity with potential Φ,

$$\begin{eqnarray}&&{\partial }_{t}\rho -{\rm{\nabla }}\left(D{\rm{\nabla }}\rho +\tau \rho {\rm{\nabla }}{\rm{\Phi }}\right)=0,\end{eqnarray} \tag{ D5 }$$

where we use the terminal velocity ansatz $v=\tau g=-\tau g{\rm{\nabla }}{\rm{\Phi }}$, with g the gravitational acceleration. With a linearization in density $\rho ={\rho }_{0}+\rho ^{\prime}$, which will also lead to a linearization in Φ, we can simplify this using ∇ρ₀ = 0 and ∇Φ₀ = 0 to

$$\begin{eqnarray}&&{\partial }_{t}\rho ^{\prime} -D{{\rm{\nabla }}}^{2}\rho ^{\prime} -\tau {\rho }_{0}{{\rm{\nabla }}}^{2}{\rm{\Phi }}^{\prime} =0.\end{eqnarray} \tag{ D6 }$$

We replace ${\rm{\Phi }}^{\prime}$ via the Poisson equation

$$\begin{eqnarray}&&{\partial }_{t}\rho ^{\prime} -D{{\rm{\nabla }}}^{2}\rho ^{\prime} -\tau {\rho }_{0}4\pi G\rho ^{\prime} =0.\end{eqnarray} \tag{ D7 }$$

With the usual plane wave ansatz $\rho ^{\prime} ={\rho }_{a}{e}^{-i\left(\omega t-{kx}\right)}$ we get

$$\begin{eqnarray}&&-i\omega +{k}^{2}D-\tau {\rho }_{0}4\pi G=0,\end{eqnarray} \tag{ D8 }$$

and it is obvious that all waves with $k\lt {k}_{c}=\sqrt{4\pi G{\rho }_{0}\tfrac{\tau }{D}}$ will be unstable and collapse. If we put in the Hill density, we receive

$$\begin{eqnarray}&&{k}_{c}=3\sqrt{\displaystyle \frac{\tau }{D}}{\rm{\Omega }},\end{eqnarray} \tag{ D9 }$$

which again we express in wavelength:

$$\begin{eqnarray}&&{\lambda }_{c}=2\pi \frac{1}{3}\sqrt{\frac{\delta }{\mathrm{St}}}H=2\pi {l}_{c};\end{eqnarray} \tag{ D10 }$$

thus, our numerical setup is linear stable to self-gravity, because L < 2πl_c , but once streaming instability has created nonlinear perturbations, those can collapse to planetesimals, if diffusion is weak enough as stated by L > 2l_c .

D.2. Particle Scale Height

We can also ask for the scale height that a particle layer would have if being in equilibrium between sedimentation and turbulent diffusion. For no self-gravity and particles with Stokes numbers larger than the dimensionless diffusivity δ this would be the well-known result

$$\begin{eqnarray}&&h=\sqrt{\displaystyle \frac{\delta }{\mathrm{St}}}H,\end{eqnarray} \tag{ D11 }$$

that is, if vertical gravity stems purely from the star. But in the case of reaching Hill density in the midplane, self-gravity is an order of magnitude stronger than the stellar gravity. We reuse the above condition for equilibrium from Equation (D5):

$$\begin{eqnarray}&&{\partial }_{t}\rho =-{\partial }_{z}\left(D{\partial }_{z}\rho +\tau \rho {\partial }_{z}{\rm{\Phi }}\right)=0,\end{eqnarray} \tag{ D12 }$$

which leads to the differential equation:

$$\begin{eqnarray}&&{\partial }_{z}^{2}\mathrm{ln}\rho =-\displaystyle \frac{\tau }{D}4\pi G\rho .\end{eqnarray} \tag{ D13 }$$

If we express ρ in Hill density as above and combine the remaining terms in our above-defined critical length ${l}_{{\rm{c}}}$, this is simply

$$\begin{eqnarray}&&{\partial }_{z}^{2}\mathrm{ln}\rho =-\displaystyle \frac{1}{{l}_{{\rm{c}}}^{2}}\displaystyle \frac{\rho }{{\rho }_{\mathrm{Hill}}},\end{eqnarray} \tag{ D14 }$$

which has the analytic solution

$$\begin{eqnarray}&&\rho (z)={\rho }_{\mathrm{Hill}}\left[1-{\tanh }^{2}\left(-\displaystyle \frac{z}{\sqrt{2}{l}_{c}}\right)\right]=\displaystyle \frac{{\rho }_{\mathrm{Hill}}}{{\cosh }^{2}\left(-\tfrac{z}{\sqrt{2}{l}_{c}}\right)}.\end{eqnarray} \tag{ D15 }$$

Thus, l_c is a truly versatile value for dust layers and clumps likewise (see Figure D1). Note that collapse here can only occur in 2D or 3D, because then the collapse time shrinks faster with length than diffusion time can increase. But in 1D a flat layer would reexpand if compressed below its vertical equilibrium height l_c because the gravitational potential at the surface of a flat sheet does not depend on the thickness of that sheet.

Zoom In Zoom Out Reset image size

E.1. Used Method: Pencil Code

For our numerical investigations we use the Pencil Code ³ (see Brandenburg 2001; Brandenburg & Dobler 2002, 2005; Youdin & Johansen 2007 for details). The Pencil Code is a general numerical solver, here used on a finite-difference hydrodynamical code using sixth-order symmetric spatial derivatives and a third-order Runge–Kutta time integration. The simulations are done in the shearing-sheet approximation (see Goldreich & Lynden-Bell 1965; Hawley & Balbus 1992; Brandenburg et al. 1995), a Cartesian coordinate system corotating with Keplerian frequency Ω at distance R from the star. Thus, all quantities have to be interpreted as being local, e.g., the shear is linearized via

$$\begin{eqnarray}&&{u}_{y}^{(0)}=-(3/2){\rm{\Omega }}x,\end{eqnarray} \tag{ E1 }$$

with x the radial coordinate in the simulation frame. All simulations are dimension free; hence, time and scaling can be chosen arbitrarily, e.g., by defining the distance to the star. The coordinate system $\left({{\boldsymbol{e}}}_{x},{{\boldsymbol{e}}}_{y},{{\boldsymbol{e}}}_{z}\right)$ can be identified as $\left({{\boldsymbol{e}}}_{r},{{\boldsymbol{e}}}_{\varphi },{{\boldsymbol{e}}}_{z}\right)$. The boundary conditions are periodic in y-/z-direction and shear-periodic in x-direction. We perform our simulations in 3D (Johansen et al. 2007a), but with only one grid cell in z-direction; see Section E.5. This means that when one suppresses modes in the vertical direction, diffusion by SI is weaker in this direction anyway and also avoids the necessity to use a thin-disk approximation or even a 2D gravity approach. To ensure dissipation on grid scale, sixth-order hyperdissipation terms are used (Lyra et al. 2008, 2009), since the Pencil Code high-order scheme has only marginally numerical dissipation.

All particles used in the simulations are Lagrangian superparticles, each representing a swarm of identical particles interacting with the gas as a bulk. Their properties, e.g., density, are smoothed out to the neighboring grid cells via the Triangular Shaped Cloud (TSC) scheme (Youdin & Johansen 2007).

E.2. Physical Model

The simulations solve the Navier–Stokes equation for the gas and the particle motion in a shearing box approximation. The gas velocity ${\boldsymbol{u}}$ relative to the Keplerian shear is evolved via

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+({\boldsymbol{u}}\cdot {\rm{\nabla }}){\boldsymbol{u}}+{u}_{y}^{(0)}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial y}=2{\rm{\Omega }}{u}_{y}\hat{x}-\displaystyle \frac{1}{2}{\rm{\Omega }}{u}_{x}\hat{y}+{{\rm{\Omega }}}^{2}z\hat{z}\end{eqnarray} \tag{ E2 }$$

$$\begin{eqnarray}&&\ -\ \displaystyle \frac{1}{{\rho }_{{\rm{g}}}}{c}_{s}^{2}{\boldsymbol{\nabla }}{\rho }_{{\rm{g}}}-\displaystyle \frac{{\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}}{{\tau }_{{\rm{s}}}}\left[{\boldsymbol{u}}-{\boldsymbol{v}}\right]+{f}_{\nu }({\boldsymbol{u}},{\rho }_{{\rm{g}}}),\end{eqnarray} \tag{ E3 }$$

where the second and third terms on the left-hand side are the advection terms by the perturbed velocity and by the shear flow, respectively. On the right are the terms for Coriolis force, the pressure gradient (with ${\boldsymbol{P}}={c}_{s}^{2}{\boldsymbol{\nabla }}\rho$), the particle–gas drag interface, and the viscosity term. The pressure gradient is split up into a global enforced pressure gradient via η (see Table 1), which is acting on the gas rather than the particles (compare with Athena code) and in the local contribution from actual evaluated gas density in the simulation domain.

The gas density is evolved with the continuity equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial t}+\left({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}\right){\rho }_{{\rm{g}}}+{u}_{y}^{\left(0\right)}\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial y}=-{\rho }_{{\rm{g}}}{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}+{f}_{{\rm{D}}}\left({\rho }_{{\rm{g}}}\right).\end{eqnarray} \tag{ E4 }$$

The particles are evolved via

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{D}}{x}^{\left(i\right)}}{{\rm{D}}t}={{\boldsymbol{v}}}^{\left(i\right)}+{v}_{y}^{\left(0\right)}\hat{y},\end{eqnarray} \tag{ E5 }$$

with Keplerian orbital velocity ${v}_{y}^{\left(0\right)}$ and particle velocity ${{\boldsymbol{v}}}^{\left(i\right)}$, which is evolved similarly to the gas

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{D}}{{\boldsymbol{v}}}^{(i)}}{{\rm{D}}t}=2{\rm{\Omega }}{v}_{y}^{(i)}\hat{x}-\displaystyle \frac{1}{2}{\rm{\Omega }}{v}_{x}^{(i)}\hat{y}-{{\rm{\Omega }}}^{2}z\hat{z}+\displaystyle \frac{1}{{\tau }_{{\rm{s}}}}[{{\boldsymbol{v}}}^{(i)}-{\boldsymbol{u}}({x}^{(i)})]\end{eqnarray} \tag{ E6 }$$

but without the gas pressure gradient acting on it. The interface between gas and particles is determined by the gas and dust densities ${\rho }_{{\rm{g}}}$ and ${\rho }_{{\rm{d}}}$ and friction time ${\tau }_{{\rm{s}}}$. As it is typically, in this paper the friction time is expressed in orbital periods, called the Stokes number $\mathrm{St}={\tau }_{{\rm{s}}}{\rm{\Omega }}$.

The gravitational potential is calculated by solving the nondimensional form of the Poisson equation

$$\begin{eqnarray}&&{\left(H{\rm{\nabla }}\right)}^{2}{\rm{\Phi }}{/{c}_{{\rm{s}}}}^{2}=\hat{G}\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\end{eqnarray} \tag{ E7 }$$

via the Fourier method (Johansen et al. 2007a); hence, ${\rm{\Phi }}\left({\boldsymbol{x}}\right)=\sum _{k}{{\rm{\Phi }}}_{k}\exp \left({\rm{i}}{\boldsymbol{k}}\cdot {\boldsymbol{x}}\right)$ with spatial wavenumber ${\boldsymbol{k}}$ and ${{\rm{\Phi }}}_{k}=-4\pi \hat{G}{\tilde{\rho }}_{k}/{\left|{\boldsymbol{k}}\right|}^{2}$. Here ${\tilde{\rho }}_{k}$ is the Fourier amplitude and $\hat{G}$ is the self-gravity parameter that one gets by adopting the Pencil Code unit system for shearing box simulations of ${c}_{{\rm{s}}},\gamma ,{\rho }_{{\rm{g}},0},{\rm{\Omega }},H=1$.

E.3. Numerical Model

First and foremost, we are interested in the particle diffusivity δ since it will allow us to predict whether collapse can occur or not. Therefore, we start with gravity switched off in order to get the simulation in a saturated streaming instability state. In this gravity-free state a particle tracking scheme can be used to measure the pure diffusivity of the streaming instability, as explained in Section E.6. Once this is achieved, gravity is switched on in a fashion that sets the initial dust density to be the Hill density, which is sufficient to ensure collapse if streaming instability does not prevent it.

Since we demand a certain dust-to-gas ratio for our study, the only way to set the initial simulation dust density to Hill density is by altering the gravitational constant $\hat{G}$ such that

$$\begin{eqnarray}&&f\cdot {\rho }_{\mathrm{Hill}}=\displaystyle \frac{9f}{4\pi }\displaystyle \frac{M}{{D}^{3}}={\rho }_{\mathrm{tot}}={\rho }_{{\rm{d}}},\end{eqnarray} \tag{ E8 }$$

with stellar mass M and distance of the particle cloud from the central star D. We set the total density to the particle density, since the gas density stays constant throughout the collapse and thus does not contribute to the gravity acting in the simulation. Here the parameter f is introduced to alter the internal density in terms of Hill density. This equation can be simplified by using the dust-to-gas ratio $\varepsilon =\tfrac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}$ and by using $M=\tfrac{{D}^{3}{{\rm{\Omega }}}^{2}}{\hat{G}}$, resulting in

$$\begin{eqnarray}&&\displaystyle \frac{9f}{4\pi }\displaystyle \frac{{{\rm{\Omega }}}^{2}}{\hat{G}}={\rho }_{{\rm{g}}}\cdot \varepsilon .\end{eqnarray} \tag{ E9 }$$

Solving for $\hat{G}$ results in

$$\begin{eqnarray}&&\hat{G}=\displaystyle \frac{9f}{4\pi }\displaystyle \frac{{{\rm{\Omega }}}^{2}}{{\rho }_{{\rm{g}}}\cdot \varepsilon }.\end{eqnarray} \tag{ E10 }$$

In our case with ${\varepsilon }_{0}=3$, Ω = 1, ${\rho }_{{\rm{g}},0}=1$, and f = 1 we have to set the gravitational constant to $\hat{G}=0.2387$.

E.4. Model Setup

E.5. Collapse Criteria Validity in Our 2D Simulations

Since full 3D simulations are highly expensive compared to 2D, we here use a setup were the z-dimension has a single grid cell. Consequently, we use the same gravitational force as in a full 3D setup, but instead of evaluating the collapse of a 3D sphere we evaluate the collapse of an infinitely extended 3D cylinder. Nevertheless, here we show that the free-fall and collapse times are in fact identical in both cases.

The gravitational force on the surface of a 3D sphere with radius R and mass $m=4/3\pi {\rho }_{\mathrm{int}}{R}^{3}$ is

$$\begin{eqnarray}&&{F}_{{\rm{g}},\mathrm{sph}}=-\displaystyle \frac{{GmM}}{{R}^{2}}\hat{r},\end{eqnarray} \tag{ E11 }$$

with unit vector $\hat{r}$, since we can collapse the whole cylinder to a single line of mass M. The gravitational force on the surface of a 3D cylinder, around this line of mass with linear density λ = M/L, with cylinder length L, one gets by calculating the gravitational potential ΔΦ = 4πG ${\rho }_{\mathrm{int}}$:

$$\begin{eqnarray}&&{F}_{{\rm{g}},\mathrm{cyl}}=-\displaystyle \frac{2{Gm}\lambda }{R}\hat{r}.\end{eqnarray} \tag{ E12 }$$

Following Appendix C, we get the equation of motion for a particle on the cylinder surface as it collapses as

$$\begin{eqnarray}&&r\left(t\right)=\sqrt{{r}_{0}^{2}-4\pi {\tau }_{{\rm{s}}}{\rho }_{\mathrm{int}}{Gt}}.\end{eqnarray} \tag{ E13 }$$

Already at this point one can see a clear parallel to Equation (C4), since they only differ in the exponent of the root function. This dependence is then eliminated when solving for collapse time τ_c, and for both spheres and cylinders the collapse time is

$$\begin{eqnarray}&&{\tau }_{{\rm{c}},\mathrm{cyl}}=\displaystyle \frac{1}{4\pi {\tau }_{{\rm{s}}}{\rho }_{\mathrm{int}}G}={\tau }_{{\rm{c}},\mathrm{sph}}.\end{eqnarray} \tag{ E14 }$$

We want to stress that the point of our 2D model is rather to show that our analytic criterion of balancing the particle cloud contraction with diffusion is properly predicting the outcome of these nonlinear simulations. The fact that contraction time is identical for 2D and 3D configurations explains why the criterion is also suited for our 2D simulations.

E.6. Measuring Particle Diffusivity in a Shear Flow

The critical quantity preventing collapse is diffusivity D of the streaming instability, which can be expressed in disk units of orbits Ω and sound speed c_s

$$\begin{eqnarray}&&\delta =\displaystyle \frac{D}{{c}_{{\rm{s}}}^{2}/{\rm{\Omega }}}.\end{eqnarray} \tag{ E15 }$$

The diffusion is measured by tracking the position of a sample of at least 10⁴ superparticles and recording their travel distance with time. The time derivative of the variance of the resulting travel distance histogram gives directly the diffusion D by using

$$\begin{eqnarray}&&D=\displaystyle \frac{1}{2}\displaystyle \frac{\partial {\sigma }_{\mathrm{Gauss}}^{2}}{\partial t},\end{eqnarray} \tag{ E16 }$$

with Gaussian variance ${\sigma }_{\mathrm{Gauss}}^{2}$ of the distribution, as introduced in Johansen & Youdin (2007). This leads to a mean travel distance from the initial particle positions of $\left\langle {r}^{2}\left(t\right)\right\rangle ={Dt}$ after a time t. The diffusivity is measured in the saturated phase of the streaming instability for each simulation before gravity is switched on (see Figure E1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

E.7. Error Bar Estimation

The error in diffusivity Δδ is estimated by calculating the standard deviation of diffusivity time series $D\left(t\right)$; see Equation (E16). From this one gets the error in the critical length scale via

$$\begin{eqnarray}&&{\rm{\Delta }}{l}_{{\rm{c}}}=\displaystyle \frac{2}{6}{\rm{\Delta }}\delta \cdot {\delta }^{-\tfrac{1}{2}}.\end{eqnarray} \tag{ E17 }$$

E.8. Increasing Gas Pressure during the Collapse

Gas pressure might increase within the collapse phase owing to friction of particles acting on the gas, dragging it along while collapsing. The reason is that the collapse phase is a situation of high dust concentration, meaning that the momentum of the dust is large, and the Stokes number is low, so its motion is well coupled onto the gas. Shariff & Cuzzi (2015) describe this effect in numerical 1D models. They claim that it can lead to oscillations in internal dust density and particle cloud core size, hence delaying the collapse for a certain parameter range, i.e., initial dust-to-gas ratio of ε = 10–100.

In our simulations we also check for changes in gas pressure. Since we perform our simulation in the ideal gas limit, and since $P={\rho }_{{\rm{g}}}{c}_{{\rm{s}}}^{2}$, we have to check our simulations for an increase in gas density that correlates with particle cloud collapse. Figure E2 shows the time series of gas and dust density for the critical collapsing cases for both investigated Stokes numbers: Ae3L0005 and Be3L003. We find for St = 0.1 no change in gas density. The strongest change in gas density is happening far after the planetesimal has formed. For St = 0.01 we indeed find a correlated increasing gas pressure, being slowly built up while the dust cloud is collapsing. But the change in pressure is with Δp ≈ 0.01 rather small, and consequently for this setup we can assume to not be in the suspension regime; though gas pressure might have an influence at the unresolved scales, it will not prevent the collapse.

Zoom In Zoom Out Reset image size

E.9. Effects of Particle Collisions during the Collapse

To justify that we are allowed to neglect particle–particle collisions in our numerical experiments, we have to estimate the collision timescale and compare it to the collapse timescale. The collision time per particle is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{coll}}={\lambda }_{\mathrm{free}}/{v}_{\mathrm{rms}},\end{eqnarray} \tag{ E18 }$$

i.e., the ratio of mean free path of a particle and the particle bulk rms velocity. The mean free path is a function of particle number density n with a certain size a and their combined cross section 4πa²:

$$\begin{eqnarray}&&{\lambda }_{\mathrm{free}}=\displaystyle \frac{1}{4\pi {{na}}^{2}}.\end{eqnarray} \tag{ E19 }$$

The number density has to be calculated from the solid density ρ_• of the particle and the mass density of the pebble cloud, which we express as multiples f of the Hill density ${\rho }_{\mathrm{dust}}=f{\rho }_{\mathrm{Hill}}$

$$\begin{eqnarray}&&n=\displaystyle \frac{f{\rho }_{\mathrm{Hill}}}{\tfrac{4}{3}\pi {a}^{3}{\rho }_{\bullet }},\end{eqnarray} \tag{ E20 }$$

and thus

$$\begin{eqnarray}&&{\lambda }_{\mathrm{free}}=\displaystyle \frac{a}{3f}\displaystyle \frac{{\rho }_{\bullet }}{{\rho }_{\mathrm{Hill}}},\end{eqnarray} \tag{ E21 }$$

which could be explicitly calculated if we knew the actual particle size. This is only possible if one could define all physical parameters entering the relation between Stokes number and particle size, i.e., stellar mass, distance to the star, density, gas temperature, and the porosity of the dust. But we know that the mass density of the pebble cloud equals the Hill density or multiples of it, plus the Stokes number is St = 0.1. With

$$\begin{eqnarray}&&\mathrm{St}={\tau }_{{\rm{s}}}{\rm{\Omega }}=\displaystyle \frac{a{\rho }_{\bullet }{\rm{\Omega }}}{\rho {c}_{{\rm{s}}}}\end{eqnarray} \tag{ E22 }$$

this gives a size of

$$\begin{eqnarray}&&a=\displaystyle \frac{H\mathrm{St}{\rho }_{\mathrm{Hill}}}{{\rho }_{0}\varepsilon },\end{eqnarray} \tag{ E23 }$$

where we express the gas density as Hill density per dust-to-gas ratio ε. Combining both expressions results in

$$\begin{eqnarray}&&{\lambda }_{\mathrm{free}}=\displaystyle \frac{H\ \mathrm{St}}{3f\varepsilon }.\end{eqnarray} \tag{ E24 }$$

With our run parameters ε = 3 and St = 0.1 this relates to

$$\begin{eqnarray}&&{\lambda }_{\mathrm{free}}=\displaystyle \frac{0.01}{f}H.\end{eqnarray} \tag{ E25 }$$

This means that for all simulations, inside their initial homogeneous particle distribution, the mean free path is larger than the smallest expected critical length of ${l}_{c}\simeq 0.004H$. Within late-stage particle overdensities with f ≥ 10 this now changes to ${\lambda }_{\mathrm{free}}\geqslant 0.001H$, but the length scale of the overdensities is still smaller around l ≈ 0.0001H.

We conclude that in all clumps found in our simulations the mean free path is equal to or larger than the clump size itself. When the mean free path indeed becomes comparable to the clump size but the particle rms speed is less than the collapse velocity of the clump, then the collision timescale will still be longer than the collapse timescale. Our derivation here is equivalent to the discussion by Youdin & Lithwick (2007).

E.10. Comparison to Estimates in the Literature

Nesvorný et al. (2010) find that for a KBO with a radius of 250 km at 30 au in an MMSN (Hayashi 1981) with 10 g cm⁻² local surface density the ratio between collision time and friction time should be

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{\mathrm{coll}}}{{\tau }_{{\rm{f}}}}\approx 0.05\sqrt{\displaystyle \frac{{R}^{2}}{30\,\mathrm{au}}}\displaystyle \frac{250\,\mathrm{km}}{{R}_{\mathrm{eq}}}{f}_{H}^{7/2},\end{eqnarray} \tag{ E26 }$$

which would define the radius at which friction and collisions are equal to 12.5 km or 25 km in diameter. This is smaller than we would have estimated above, so we recapitulated their estimate. They used solid density of 2 g cm⁻³ and some additional order-of-magnitude shortcuts.

In communication with the authors of Nesvorný et al. (2010) we found that for the nebula models in this paper the critical size to have collisions dominate over friction is larger than 100 km (See Figure 10) and smaller than the 500 km considered in Nesvorný et al. (2010).

It is therefore safe to neglect collisions in the present work (with ε = 3). In follow-up 3D studies we will treat them correctly in order to get a better understanding of the final outcome of planetesimals (e.g., multiplicity and spin rate) from the described process of self-gravity.

F.1. Varying the Pressure Gradient

F.2. Varying the Initial Dust-to-gas Ratio