1 Methods

We simulate the motion of gas and solids in a local, corotating patch of a protoplanetary disc in the shearing box approximation. The axes are oriented such that x points outwards along the cylindrical radius, y points along the main rotation direction, and z points vertically out of the disk parallel to the Keplerian rotation vector \( {\varvec \Upomega} \) . The gas is modelled on a fixed grid, while the solid particles are treated as superparticles, each having a space and velocity vector (Johansen et al. 2006). We apply two-way drag forces between the gas and the particles through a particle-mesh scheme that ensures momentum conservation and allows to follow both the linear and the non-linear evolution of the streaming instability (Youdin et al. 2007, Johansen et al. 2007).

Since we do not consider magnetic fields in this work, the most important dynamics comes through the streaming instability that develops in the relative motion of gas and solids (Youdin et al. 2005) and through Kelvin-Helmholtz instabilities arising in the velocity difference of the gas between the dust-rich mid-plane layer and the dust-free regions above and below the mid-plane (Weidenschilling 1980; M. Sekiya 1998; Youdin et al. 2002; Gómez et al. 2005; Johansen et al. 2006). The streaming instability has its highest growth rate for particle sizes that couple by friction to the gas on an orbital time-scale, and at r = 5 AU in the Minimum Mass Solar Nebula this corresponds to approximately particles of 30 cm in radius.

In the presented simulations we consider four different particle sizes simultaneously, with \(\Upomega \tau_{\rm f}=0.25,0.50,0.75,1.00\). Here \(\tau_{\rm f}\) is the friction time, the time-scale upon which the particle velocity would approach the gas velocity is absence of other forces. The growth from \(\mu\hbox{m}\) dust grains to cm-sized rocks likely occurs by coagulation (Weidenschilling 1997; Dullemond et al. 2005; Brauer et al. 2008). This important growth step is not considered in this work, but experiments on sticking properties and growth of dust grains show that it is possible to form macroscopic bodies by dust coagulation (Blum et al. 2000; Wurm et al. 2001).

2 Results

We first consider the non-linear evolution of the streaming instability without the self-gravity of the particles.

The effect of varying the global solids-to-gas ratio ε is shown in Fig. 1. The top panels show space-time plots of particle density versus radial coordinate x and time t, averaging over the azimuthal and vertical directions. For the canonical solids-to-gas ratio ε = 0.01 there is not a trace of persistent overdensities in the solids, However, as the solids-to-gas ratio is increased to ε = 0.02 and ε = 0.03 strong overdensities occur in the particle layer. One sees clearly from the angle between the structures and the x-axis how the radial drift is locally reduced in these overdense regions, leading to almost completely vertical structures for ε = 0.03. The maximum density in the box actually increases by a dramatic two orders of magnitude as the solids-to-gas ratio is increased from 0.01 to 0.02, far more than expected from the modest increase in the global amount of solids. Rather the increased back-reaction drag force on the gas instigates traffic jam like pile-ups, due to a non-linear manifestation of the streaming instability already found in Johansen et al. (2006, 2007).

Fig. 1
figure 1

The particle density versus radial coordinate x and time t, for 3-D simulations of multiple particle sizes \({\varvec \Upomega}\tau_{\rm f}=0.25,0.50,0.75,1.00\), shown for three different values of the global solids-to-gas ratio \(\epsilon={\varvec \Upsigma}_{\rm p}/{\varvec \Upsigma}_{\rm g}\) (top panels). Traffic jams in the radial drift flow set in from ε = 0.02, but are absent in the ε = 0.01 case where the back-reaction drag force is too weak. The bottom panels show the maximum particle density in the box (black line) and the scale height of the particle mid-plane layer (red/gray line) as a function of time. Increasing the solids-to-gas ratio above 0.01 leads to a dramatic increase in maximum density by up to two orders of magnitude. The mid-plane layer also gets slightly thinner (right axis)

The width of the mid-plane layer is shown in the right axis of the bottom panels of Fig. 1. As the solids-to-gas ratio is increased, the mid-plane layer contracts accordingly, as the turbulent motions get weaker. However, there is only around a factor of three decrease in the scale height of the solids when going from ε = 0.01 to ε = 0.03. Together with the increased amount of solids this would lead to a factor ten increase in the density of solids, if the turbulence had otherwise remained unchanged. However the occurrence of traffic jams gives a maximum density at least an order of magnitude higher than this.

For the simulation with ε = 0.02 we turned on the self-gravity of the particles just before a clumping event occurred. The overdense filament formed as usual, due to the streaming instability, but the self-gravity caused the condensation of several gravitationally bound objects out of the particle layer. In Fig. 2 we show a snapshot of the particle positions around two orbits after the concentration event started. Five newly born planetesimals are clearly seen in the sea of dust particles. Each of the planetesimals has a mass corresponding to a solid body of several hundred kilometers in radius. The planetesimals show clear sign of accretion of rocks and boulders from the surroundings and indeed grow rapidly in mass even after the initial contraction.

Fig. 2
figure 2

Five newly born planetesimals embedded in the sea of dust whence they formed. This plot was created with Partiview and looks down on the mid-plane from a slightly inclined angle. The planetesimals, all with radii of a few 100 km, show clear signs of active accretion of the remaining rocks and boulders

2.1 Conclusions

A curious feature of the presented simulations is a strong, non-linear dependence of the turbulent state on the global value of the solids-to-gas ratio. At a canonical solids-to-gas ratio of 1%, spread over four discrete bins covering a bit less than an order of magnitude in particle size, a strong streaming instability develops and puffs up the mid-plane layer to such widths that clumping is strongly suppressed. This environment is about the least beneficial imaginable for particle growth, with low densities and high collision speeds. However, already a doubling the solids-to-gas ratio to 0.02 leads to strong clumping with local particle densities reaching a thousand times the gas density. Such a distinctly non-linear behaviour likely comes about because the turbulent state of the streaming instability has a counterintuitive feature: the strength of the turbulence increases with decreasing solids-to-gas ratio. As turbulence develops, the mid-plane layer is puffed up, leading to stronger turbulence and to an ever increasing width of the mid-plane layer.

The solids-to-gas ratio dependence indicates that the overall metallicity of the protoplanetary disc matters. There is a strong correlation between high metallicity and the probability of finding close-in extrasolar giant planets (Gonzalez 1997; Santos et al. 2001). Our findings certainly support that planet formation is facilitated in discs with a high metallicity, but other effects, such as core accretion and planet migration, may benefit equally from higher metallicity, so it is not per se clear how important a metallicity dependent planetesimal formation is for this picture.

An important conclusion of the current study is that magnetised turbulence is not crucial for planetesimal formation by self-gravity. Long-lived high pressure regions that form spontaneously in magnetorotational turbulence can concentrate particles by up to two orders of magnitude (Johansen et al. 2006). By a combination of high pressure trapping and streaming instability it was demonstrated in Johansen et al. (2007) that gravitational collapse could occur at r = 5 AU in a disc with four times the Minimum Mass Solar Nebula surface density and a solids-to-gas ratio of 0.01. For the pure streaming instability, on the other hand, a slightly higher solids-to-gas ratio is needed, ε = 0.02, but planetesimals can form for much lower surface densities, around the Minimum Mass Solar Nebula value or lower. This has important prospects for the formation in planetesimals in magnetically dead regions in protoplantary discs, occurring e.g. around the mid-plane in the inner 10–20 AU (Sano et al. 2000; Fromang et al. 2002; Semenov et al. 2004).

The question of collisional fragmentation has not been addressed in this work. Preliminary studies indicate collision speeds between particles in the range 1–5 m per second during the concentration events and during the collapse and accretion phases. This is approximately a factor two times lower than the collision speeds found in magnetorotational turbulence (Johansen et al. 2007) and supports the idea that weakly turbulent discs provide a good environment for planetesimal formation. However, the full treatment of the sticking and collisional fragmentation of various particle sizes during the gravitational collapse will require improved numerical models as well as better experimental knowledge of the outcome of collisions between rocks and even boulders (Blum et al. 2008).