© U. Schäfer and A. Johansen 2022

1 Introduction

The streaming instability has emerged in the last decade as the leading candidate among processes to cause the formation of planetesimals (Chiang & Youdin 2010; Johansen et al. 2014), a key step in the growth from sub-micron-sized dust grains to planets. The streaming instability was discovered analytically as a linear instability, that is to say an exponential growth of small linear perturbations, by Youdin & Goodman (2005). Subsequently, it was shown how the instability evolves from a linear into a non-linear phase in numerical simulations (Johansen & Youdin 2007; Bai & Stone 2010a; Kowalik et al. 2013; Mignone et al. 2019). These simulations reveal that the instability in its non-linear regime causes dust to accumulate in largely axisym-metric filaments (Johansen & Youdin 2007; Bai & Stone 2010b; Kowalik et al. 2013; Yang & Johansen 2014; Li et al. 2018). Over-densities in these filaments can be sufficiently large to contract owing to their self-gravity and form planetesimals with typical sizes of tens or hundreds of kilometres (Johansen et al. 2007, 2009, 2015; Simon et al. 2016; Schäfer et al. 2017; Abod et al. 2019).

In all likelihood, the streaming instability is active everywhere dust is present in protoplanetary disks. The instability requires only rotating gas and dust which are coupled via mutual drag, and a gas pressure gradient with respect to the radial distance to the star (Youdin & Goodman 2005). This gradient entails a slight deviation of the orbital speed of the gas and dust components of protoplanetary disks from the Keplerian speed as well as a radial drift of the two components, especially a drift of the dust towards the star (e.g. Adachi et al. 1976; Nakagawa et al. 1986; Brauer et al. 2007). Local dust overdensities drift more slowly than the surrounding dust, resulting in further dust accumulation in the overdensities, a further reduction of their drift speed, and ultimately instability – the streaming instability.

However, the ability of the streaming instability to give rise to dust concentration that is strong enough to trigger planetesimal formation depends on three parameters: the Stokes number of the dust St, the dust-to-gas surface density ratio Σ_d/Σ_g¹, and the strength of the radial gas pressure gradient (Johansen et al. 2009; Bai & Stone 2010b,c; Drazkowska & Dullemond 2014; Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021). The latter is commonly expressed in terms of the dimensionless parameter Π as introduced by Bai & Stone (2010b), which is typically of the order of 0.01 or 0.1. Sekiya & Onishi (2018) propose, and corroborate using simulations, a dependence on the parameter (2π)^1/2 Σ_d/(Σ_gΠ) rather than on Σ_d/Σ_g and Π individually. This parameter expresses the dust-to-gas density ratio in the dust mid-plane layer under the assumption that the thickness of this layer is regulated by the streaming instability. Carrera et al. (2015), Yang et al. (2017), and Li & Youdin (2021) each conducted a suite of simulations to obtain, for a given Stokes number and pressure gradient, a minimum dust-to-gas surface density ratio necessary for planetesimal formation. Ever-smaller minimum values are found in these studies, with Li & Youdin (2021) inferring a value of Σ_d/Σ_g = 2% for St = 0.01 and Π = 0.05 as well as Σ_d/Σ_g = 0.6% for the same pressure gradient but St = 0.1.

These numerical results raise the question of whether the streaming instability on its own can be expected to cause planetesimal formation for the dust sizes and dust-to-gas ratios that are observed in protoplanetary disks. We focus here on Class II disks as we numerically model such disks in this study. Notably, ALMA has enabled measurements of dust sizes and dust-to-gas ratios in these disks in recent years owing to its high spatial resolution and sensitivity.

Maximum dust sizes are most frequently inferred from either the opacity spectral index of thermal dust emission (e.g. Draine 2006; Testi et al. 2014) or from the polarisation of the dust emission owing to self-scattering by other dust grains or grain alignment (Kataoka et al. 2015, Kataoka et al. 2017; Mori & Kataoka 2021). The former method yields millimetre- or centimetre-sizes, with a trend towards smaller sizes at greater distances to the star (Pérez et al. 2012, 2015; Tazzari et al. 2016, 2021a,b; Sierra et al. 2019; Macías et al. 2019, Macías et al. 2021; Carrasco-González et al. 2019; Maucó et al. 2021). In contrast, sizes of not more than a few hundreds of microns are derived applying the latter method (Kataoka et al. 2017; Bacciotti et al. 2018; Ohashi & Kataoka 2019; Ohashi et al. 2020; Lin et al. 2020; Mori & Kataoka 2021).

To explain the inconsistency between these two kinds of measurements, it has been proposed that the sizes obtained from spectral indices might be overestimates that stem from optically thick emission being misinterpreted as optically thin, in particular if self-scattering is neglected (Liu 2019; Zhu et al. 2019; Lin et al. 2020; Ohashi et al. 2020). However, this degeneracy can be broken when multi-wavelength observations are considered, with the sizes inferred from modeling such observations still lying in the millimetre-to-centimetre range (Sierra et al. 2019; Carrasco-González et al. 2019; Macías et al. 2021; Maucó et al. 2021, see also Tapia et al. 2019)². On the other hand, the deviations from spectral index measurements can be alleviated when the shape and chemical composition of dust grains are taken into account in polarisation measurements (Kirchschlager & Bertrang 2020; Yang & Li 2020; Brunngräber & Wolf 2021).

The canonical dust-to-gas ratio in the Milky Way interstellar medium amounts to ~ 1%. The ratio is not constant, however, but varies by factors of a few even within our galaxy (Spitzer 1978; Sodroski et al. 1997) – particularly also in molecular clouds (Liseau et al. 2015) – and by orders of magnitude when considering galaxies with different metallicities (Brinchmann et al. 2013; Rémy-Ruyer et al. 2014). The ratio of dust to gas mass in protoplanetary disks is observed to typically be higher than the canonical interstellar medium value, with ratios of ~10% being common (Ansdell et al. 2016; Miotello et al. 2017; Long et al. 2017; Wu et al. 2018; Macías et al. 2021). Nevertheless, measurements of both dust and gas masses are subject to significant uncertainties. On the one hand, the total dust mass can be reliably inferred from the thermal dust emission only if the emission is optically thin (e.g. Andrews 2020). On the other hand, H₂ – the most abundant gas molecule – is difficult to detect, conversion from the observed mass of CO – the second most abundant molecule – to the total gas mass is complex and prone to error, and observations of the promising tracer molecule HD are still rare (e.g. Miotello et al. 2017; Andrews 2020; Anderson et al. 2022). Based on the composition of the solar photosphere as well as meteorites, Lodders (2003) finds that condensates should have constituted ~1.5% of the mass of the protoplanetary disk which evolved into the Solar System.

Comparing these estimates of the dust sizes and dust-to-gas ratios in protoplanetary disks with the results of the above-mentioned parameter studies (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021) shows that planetesimal formation induced by the streaming instability alone is possible at around or slightly higher than the canonical interstellar medium dust-to-gas ratio. In this context, a variety of mechanisms have been proposed to act in concert with the streaming instability. Among them are two other hydrodynamical instabilities, the subcritical baroclinic instability (Raettig et al. 2015, 2021) and the vertical shear instability (Lehmann & Lin 2022; Schäfer et al. 2020, hereafter SJB20).

The vertical shear instability is a hydrodynamical instability that arises if the gas orbital velocity varies with height – for instance owing to a radial temperature gradient, which is omnipresent in protoplanetary disks (e.g. Andrews 2020) – and the gas cooling timescale is sufficiently short (Nelson et al. 2013; Lin & Youdin 2015). Stoll & Kley (2016) find that the vertical shear instability gives rise to transient pressure fluctuations in which dust accumulates, resulting in overdensities of several times the mean initial dust density. Similarly, dust is concentrated in vortices formed by the instability (Flock et al. 2017, 2020; Lehmann & Lin 2022).

In a previous paper (SJB20), we present evidence that over-densities induced by the vertical shear instability trigger further concentration by the streaming instability. The aim of this paper is to examine to what extent a combination of the vertical shear instability and the streaming instability leads to a reduction of the minimum values of dust-to-gas surface density ratio and dust size or Stokes number required for planetesimal formation compared to if only the streaming instability is considered. To this end, we performed a parameter study adopting the numerical model of SJB20. Our two-dimensional, axisymmetric simulations cover the full scale of protoplanetary disks, while adaptive mesh refinement permits us to locally resolve the formation of dust overdensities in the mid-plane layer owing to the two instabilities.

In Sect. 2, we describe our numerical model. We examine dust concentration and the relation between dust overdensities and pressure bumps in Sect. 3. In Sect. 4, we analyse which metrics can be applied to gauge the potential for planetesimal formation in models such as ours that do not include self-gravity. This is followed by a presentation of the threshold values of dust-to-gas surface density ratio and dust size that are necessary for the streaming instability alone or it and the vertical shear instability in combination to induce dust concentration that is sufficient for planetesimal formation in Sect. 5. We discuss the implications and limitations of our study in Sect. 6 and summarise its results in Sect. 7.

2 Simulations

We applied the FLASH Code³,⁴ (Fryxell et al. 2000) to conduct simulations of the gas and dust components of protoplanetary disks, including the mutual drag between the two components and the stellar gravity. Gas and dust were modeled on a Eulerian grid and as Lagrangian particles, respectively. We summarise in this section the most important aspects of the simulations, and refer to SJB20 for a more detailed description.

Table 1 lists all simulations and the parameters that distinguish them. The simulations names are composed of, in this order, the instabilities which are simulated, the initial dust-to-gas surface density ratio, as well as the dust size or Stokes number. We discern two kinds of simulations in which both the streaming instability and the vertical shear instability are active: in the scenario SIwhileVSI, both instabilities start to grow at the same time. On the other hand, the vertical shear instability has already saturated before the streaming instability begins to operate in the scenario SIafterVSI. The latter scenario was realised by initialising the dust particles after 50 kyr instead of at the beginning of the simulations. Furthermore, to exclude the streaming instability and simulate only the vertical shear instability we neglected the drag exerted by the dust onto the gas.

2.1 Domain Size and Resolution

The geometry of our two-dimensional simulation domains is cylindrical as this is a natural choice to model protoplanetary disks. The axisymmetric domains range from 10 to 100 au in the radial dimension while encompassing 1 or 2 gas scale heights above and below the disk mid-plane. All simulations involving the vertical shear instability were performed using the larger vertical domain size of 4 scale heights since in SJB20 we show that this size is required to reproduce the turbulent strength which is found in previous numerical studies of this instability. Gas and dust were permitted to leave but not enter the domains through the radial and vertical boundaries, with the pressure being interpolated to ensure vertical hydrostatic equilibrium at the latter boundaries. To avoid that artificial behaviour caused by the inner radial boundary conditions contaminates our results, we excluded the innermost 10 au from all quantitative analysis. (The dust rapidly drifting and settling away from them renders a similar treatment of outer radial and vertical boundaries unnecessary.)

The base resolution of our simulations amounts to ten grid cells per astronomical unit. On top of that, both static and adaptive mesh refinement were employed. We allowed for up to six levels of refinement, with every level corresponding to a doubling of the resolution. Thus, the maximum resolution is equal to 320 cells per astronomical unit. In addition to simulations with these fiducial base and maximum resolutions, we conducted two simulations in which both are doubled to investigate the resolution dependence of our findings. Adaptive mesh refinement was applied to blocks of 10 × 10 cells when the number of dust particles in any cell inside these blocks exceeded ten. Blocks were derefined, on the other hand, if no particles were left in a cell. This allows us to resolve dust concentration owing to the streaming instability and the vertical shear instability in the mid-plane layer while minimising the impact of the layers away from the mid-plane that are void of dust on the computational cost of the simulations.

The usage of static mesh refinement constitutes the only difference between the numerical models we used in SJB20 and in this study. In SJB20, the resolution was increased by one or two refinement levels, respectively, where the gas density exceeds 1% or 10% of the mid-plane density at the inner radial domain boundary. Here, we instead enhanced the resolution by one or two refinement levels, respectively, within 5 au and 1 au above and below the mid-plane. Figure 1 shows that this modification of the static refinement prevents the formation of a significant bump in the gas density at a radius of ~30 au before the dust is introduced in simulations of the vertical shear instability only and of the scenario SIafterVSI. As is evident from Fig. 8 of SJB20, the presence of this bump entails a large gap in the dust density. This gap does not form in the simulations we present in this paper (see Sect. 3.2).

Fig. 1 Gas density ρg in the mid-plane (at a height z = 0) as a function of the radius r. The blue line shows the initial density profile, while the orange and green line depict the profiles after 50 kyr in simulations including the static mesh refinement employed in SJB20 and in this study, respectively. Notably, a considerable bump at a radius of ~30 au is present in the former profile, but not in the latter one.

2.2 Gas

The initial gas density decreases both with the radial distance r to the star,

$$p_g\left(z=0\right)=5.62\times {10}^{-12}{\text{gcm}}^{-3}{\left(\frac{r}{\text{10au}}\right)}^{-9/4},$$(1)

and with the height z above or below the mid-plane,

$$p_g=p_g\left(z=0\right)\exp \left[-\frac{{}_{\gamma }G{M}_s}{c_s^2}\left(\frac{1}{r}-\frac{1}{\sqrt{r^2+z^2}}\right)\right],$$(2)

where G is the gravitational constant, M_S = 1 M_⊙ the stellar mass, c_s = (γRΤ/µ)^1/2 the sound speed, γ the adiabatic index, R the ideal gas constant, Τ the temperature, and μ = 2.33 the mean molecular weight. We use the subscripts g and d to refer to gas and dust, respectively. The radial density gradient is shallower than that in the minimum mass solar nebula model, though the mid-plane density at r = 10 au is about two times higher (Hayashi 1981). The vertical density gradient ensures hydrostatic equilibrium in this dimension. We define the surface density as the integral of the density from −1 gas scale height to 1 gas scale height (rather than from –∞ to ∞). It is thus initially given by

$${\displaystyle {\sum }_g={\displaystyle {\int }_{-1H_g}^{1H_g}{p}_g\text{d}z}={10}^2{\text{gcm}}^{-2}}{\left(\frac{r}{10\text{au}}\right)}^{-1},$$(3)

$$H_g=\sqrt{\frac{c_s^2{r}^3\left(2\gamma G{M}_s-c_s^{2}r\right)}{{\left(c_s^{2}r-\gamma G{M}_s\right)}^2}}=0.846\text{au}{\left(\frac{r}{10\text{au}}\right)}^{5/4}.$$(4)

We adopt the initial gas temperature from the minimum mass solar nebula model (Hayashi 1981),

$$T=88.5\text{K}{\left(\frac{r}{10\text{au}}\right)}^{-1/2}.$$(5)

This radial temperature gradient causes a height-dependence of the orbital velocity and consequently gives rise to the vertical shear instability. In all simulations of this instability (either in combination with the streaming instability or in isolation), an isothermal equation of state was employed,

$$p=\frac{RT}{\mu }{p}_g,$$(6)

because an infinitely short gas cooling timescale provides ideal conditions for the instability (Nelson et al. 2013; Lin & Youdin 2015). Here, P is the pressure. On the other hand, in simulations of the streaming instability only we used an adiabatic equation of state,

$$p=K_{p_g^{\gamma },}$$(7)

where $K=RT{\rho }_g^{1-\gamma }/\mu$ is the polytropic constant – constant in time, but varying in space with the temperature and density distributions – and the adiabatic index γ = 5/3. This is since under these conditions the vertical shear instability is quenched by vertical buoyancy.

The radial gradients in gas density and temperature entail a pressure gradient, which is necessary for the streaming instability to be active. We express the strength of the pressure gradient in terms of the dimensionless parameter (Bai & Stone 2010b)

$$\Pi =-\frac{1}{2c_s{p}_g{\Omega }_K}\frac{\text{d}P}{\text{d}r},$$(8)

where Ω_K is the Keplerian orbital frequency. In the mid-plane, this parameter amounts to

$$\Pi \left(z=0\right)=8.2\times {10}^{-2}{\left(\frac{r}{10\text{au}}\right)}^{1/4}.$$(9)

The initial orbital velocity is chosen such that this pressure gradient and the centrifugal force balance the radial stellar gravity.

2.3 Dust

As is common in numerical models of the streaming instability including Lagrangian particles to represent the dust (Youdin & Johansen 2007; Bai & Stone 2010a), every particle in our simulations possesses the total mass and momentum of a huge number of the dust aggregates that are present in protoplanetary disks, but the drag coupling to the gas of a single such aggregate. Our simulations include 10⁶ dust particles. Initially, these are uniformly distributed in the radial dimension while their vertical positions are randomly sampled from a Gaussian distribution with a scale height equal to 10% of the gas scale height. The noise in the vertical distribution seeds the streaming instability. We find our results to be independent of whether 5 × 10⁵ or 10⁶ particles are simulated as well as of the random seed, and discuss the dependence on the initial scale height in SJB20.

Because of their homogeneous radial distribution, the particle mass is given by

$$m_d=\frac{1}{N_d}{\displaystyle {\int }_{{}_{L_p}}2\pi r{\displaystyle {\sum }_{d}dr},}$$(10)

where N_d = 10⁶ is the number of particles and L_r = 90 au the radial domain extent. Since we choose the dust-to-gas surface density ratio to be constant initially and the gas surface density is inversely proportional to the radius, the mass of every particle thus amounts to

$$m_d=1.27\times {10}^{24}g\left(\frac{{\Sigma }_{\text{d}}/{\Sigma }_{\text{g}}}{0.01}\right).$$(11)

We consider surface density ratios of 0.5%, 1%, and 2%. These ratios are equal to the mean dust-to-gas mass ratio in the Milky Way interstellar medium or smaller or greater by a factor of two. Nonetheless, they are low when compared to the observed mass ratios in protoplanetary disks (see Sect. 1).

All particles in a given simulation are of the same size or Stokes number, respectively, with sizes of 3 mm and 3 cm as well as Stokes numbers of 0.01 and 0.1 being taken into account. As detailed in Sect. 1, these are consistent with the sizes inferred from the opacity spectral index of the dust emission from protoplanetary disks, but larger than the ones derived from the polarisation of the emission. When considering the disk mid-plane, we can convert from dust size to Stokes number and vice versa using the relations

$$\begin{array}{lll}St\left(z=0\right)\hfill & =\hfill & t_{d,stop}\left(z=0\right){\Omega }_K\left(z=0\right)=\frac{a{p}_s}{c_s{p}_g\left(z=0\right)}{\Omega }_K\left(z=0\right)\hfill \\ \hfill & =\hfill & 6\times {10}^{-3}\frac{1}{\sqrt{\gamma }}\left(\frac{St}{0.01}\right){\left(\frac{r}{10au}\right)}^{-1}\hfill \end{array}$$(12)

$$a\left(z=0\right)=5.02mm\sqrt{\gamma }\left(\frac{St}{0.01}\right){\left(\frac{r}{10au}\right)}^{-1},$$(13)

where t_d,stop is the dust stopping time, a the dust size, and ρ_s = 1 g cm⁻³ is the dust material density. Here, we utilise that the dust size is smaller than the gas mean free path length at all gas densities in our model and therefore calculate the dust stopping time in the Epstein regime. Figure 2 illustrates these relations exemplarily for our model with an isothermal equation of state. The particles initially orbit with the Keplerian speed.

Fig. 2 Dust size a (left ordinate) and Stokes number St (right ordinate) as functions of the radius r. Conversion relations between dust size and Stokes number in the mid-plane of our protoplanetary disk model (Eqs. (12)–(13)) are shown for the two dust sizes (orange lines) and the two Stokes numbers (blue lines) that we simulated, assuming an isothermal equation of state.

3 Dust Concentration

3.1 Synopsis of Relevant Results from Schäfer et al. (2020)

As in SJB20, we investigate three scenarios: one in which the streaming instability operates in isolation and two in which both this instability and the vertical shear instability are active. In SIwhileVSI, both instabilities develop simultaneously. This scenario is similar to that including only the streaming instability in that the turbulence in the mid-plane dust layer is predominantly driven by the streaming instability, with the dust scale height amounting to ~1% of the gas scale height.

In SIafterVSI, on the other hand, the vertical shear instability has saturated before the streaming instability begins to grow. The vertical shear instability remains the primary source of turbulence in the dust layer in this scenario, though the streaming instability causes turbulence locally in dust overdensities. Since the vertical shear instability causes stronger turbulent motions in the vertical dimension than the streaming instability, the dust scale height is equal to ~10% of the gas scale height in this scenario.

Nonetheless, in SJB20 we show that dust overdensities are more prominent and the maximum dust-to-gas volume density ratio is higher in the scenario SIafterVSI than in the scenario SIwhileVSI. We conjecture that dust overdensities caused by the vertical shear instability trigger the streaming instability and are reinforced by it. In this section, we investigate dust concentration in the three scenarios in more detail.

3.2 Spatial Dust Distribution

To begin with, we explore the vertical distribution of the dust in our models. Figure 3 depicts the dust-to-gas volume density ratio in simulations of the streaming instability only, of the scenario SIwhileVSI, and of the scenario SIafterVSI with the same initial dust-to-gas surface density ratio of 1% and Stokes number of 0.1. In all three simulations, local maxima of the volume density ratio are not concentrated in the mid-plane, but present at all heights inside the dust layer. Nevertheless, in the simulations of the streaming instability alone and of the scenario SIwhileVSI the dust layer resembles a Gaussian distribution, with the radial average of the volume density ratio reaching its maximum near or in the mid-plane.

In contrast, in the SIafterVSI simulation the shape of the layer is wave-like and markedly deviates from a Gaussian. While the dust scale height is considerably greater in this scenario than in the other two, the vertical thickness of the dust layer is not. This illustrates that stronger large-scale turbulence regulates the scale height while weaker small-scale turbulence drives the internal diffusion in the layer. The millimetre-sized dust in the three-dimensional models of the vertical shear instability that are presented by Flock et al. (2017, 2020) possesses a similar wave-shaped distribution, both in the radial-vertical plane as in our two-dimensional model and in the azimuthal-vertical plane (see Figs. 8–9 of Flock et al. 2017 and Figs. 7 and D1 of Flock et al. 2020.)

As shown also in SJB20, we find that dust overdensities are greater in the scenario SIafterVSI than in the scenario SIwhileVSI. This is illustrated in the upper panels of Fig. 4, which depict the dust-to-gas surface density ratio in the same simulations as shown in Fig. 3. From the figure, it can further be seen that after 5 kyr and at radii less than 40 au dust over-densities are larger in SIwhileVSI than in the simulation of the streaming instability only. This can be explained by the vertical shear instability growing over several thousand years – it takes about 30 orbital periods to saturate (Stoll & Kley 2014; Flock et al. 2017), longer than the streaming instability (SJB20) – until it starts to influence the dust mid-plane layer and enhance dust overdensities in it.

The question arises whether comparable simulations of the scenarios SIwhileVSI and SIafterVSI would eventually evolve into the same state if they were continued sufficiently long. To address this question, we show the dust-to-gas volume density ratio in simulations of these two scenarios with an initial dust-to-gas surface density ratio of 2% and a dust size of 3 cm in Fig. 5. In the figure, the dust distribution is depicted 10 kyr or ~ 100 local orbits after the dust is introduced into the simulations (at the beginning in the case of the SIwhileVSI simulation and after 50 kyr in the case of the SIafterVSI simulation). As can be seen, the dust layer in SIwhileVSI continues to possess a Gaussian shape, while the dust layer in SIafterVSI still does not. We cannot exclude that the two simulations would develop into a similar state on even longer timescales, though.

In summary, the vertical shear instability and the streaming instability together induce stronger dust concentration than the streaming instability alone, independent of whether the vertical shear instability (scenario SIafterVSI) or the streaming instability (scenario SIwhileVSI) is the dominant source of turbulence in the dust layer. The concentration is strongest in the scenario SIafterVSI, despite the vertical shear instability inducing a larger dust scale height in this scenario than the streaming instability in the other two scenarios as this larger scale height is not synonymous with a greater vertical thickness of the dust layer.

Fig. 3 Ratio of dust ρd to gas volume density ρg as a function of radius r and height z, the latter given in units of gas scale heights Hg (left panels). The average over the radius range shown in the left panels is depicted as a function of height in the right panels. The top, middle, and lower panels show the dust-to-gas volume density ratio at the end of simulations of the streaming instability in isolation, the scenario SIwhileVSI, and the scenario SIafterVSI, respectively, with an initial dust-to-gas surface density ratio of 1% and a Stokes number of 0.1. The horizontal white lines in the lower panel illustrate the range of heights depicted in the upper and middle panels. In all simulations, local maxima of the volume density ratio are distributed over the entire vertical dimension of the dust layer. While the vertical dust distribution in the simulations of the streaming instability only and of the scenario SIwhileVSI is largely Gaussian, the dust layer possesses the form of a wave rather than a Gaussian shape in the SIafterVSI simulation.

Fig. 4 Ratio of dust Σd to gas surface density Σg (upper panels) and difference between azimuthal gas velocity υg,ϕ and Keplerian velocity υK (lower panels). Both quantities are depicted as a function of radius r and time t after the initialisation of the dust particles at td,init. The velocity difference is averaged over the vertical domain size, weighted by the dust mass, and expressed as a Mach number. In the absence of turbulence, this quantity is equal to the negative of the parameter Π as defined in Eq. (8). Black dots indicate that the azimuthal velocity is equal to or greater than the Keplerian velocity. Enhancements in the dust-to-gas surface density ratio and in the difference between azimuthal gas velocity and Keplerian velocity are larger at late times (t – td,init > 5 kyr) and small radii (r < 40 au) in SIwhileVSI than in the simulation of the streaming instability in isolation, and greatest at all times and radii in SIafterVSI. Surface density ratio and velocity difference can be seen to generally be correlated.

Fig. 5 Dust-to-gas volume density ratio ρd/ρg as a function of radius r and height z in simulations of the scenarios SIwhileVSI (top panel) and SIafterVSI (bottom panel) with an initial dust-to-gas surface density ratio of 2% and a dust size of 3 cm. The volume density ratio is depicted at the end of the simulations 10 kyr after the initialisation of the dust particles. The vertical dust distribution in the SIwhileVSI simulation resembles a Gaussian, but is wave-like in the SIafterVSI simulation.

Fig. 6 Dust-mass-weighted mean deviation of azimuthal gas velocity from Keplerian velocity υg,ϕ – υK, normalised by the sound speed cs, as a function of radius r and time after the dust particle initialisation t – td,init. The left and right panels depict simulations of the vertical shear instability only and of the scenario SIafterVSI, respectively. Local increases in the azimuthal gas velocity are stronger, with the velocity more frequently exceeding the Keplerian velocity, and especially at r > 30 au also more persistent in SIafterVSI.

3.3 Connection Between Dust Concentration and Pressure Bumps

The lower panels of Fig. 4 depict the difference between the azimuthal gas velocity and the Keplerian velocity. A correlation of this velocity difference and the dust-to-gas surface density ratio, which is shown in the upper panels, is apparent in our model of the streaming instability in isolation (left panels). This can be explained by the drag exerted by the dust on the gas causing the azimuthal gas velocity to be closer to the Keplerian velocity where the dust-to-gas volume density ratio is enhanced.

Nevertheless, an increased azimuthal gas velocity is not only an indicator of dust overdensities, but also of pressure bumps, that is local deviations from the global gas pressure gradient. This is since the velocity depends linearly on the pressure gradient and is equal to the Keplerian velocity if the pressure attains a local maximum and the gradient vanishes (assuming that turbulence is negligible). We note that Li et al. (2018) show that the streaming instability causes weak pressure bumps, but independently of dust concentration. Moreover, Yang & Johansen (2014) find a weak anti-correlation between dust and gas density in their streaming instability simulations, while dust overdensities and local increases in the azimuthal gas velocity are correlated in our simulations.

On the other hand, previous studies show that dust accumulates in vortices (Flock et al. 2017, 2020; Lehmann & Lin 2022) and short-lived pressure fluctuations (Stoll & Kley 2016) induced by the vertical shear instability. While vortices cannot form in our two-dimensional model, it can already be seen from Fig. 1 that the vertical shear instability indeed gives rise to pressure bumps. A comparison of the left and middle panels of Fig. 4 reveals that thus both local enhancements in the azimuthal gas velocity and dust overdensities are more pronounced at late times and small radii – when the vertical shear instability has grown sufficiently to influence the dynamics in the dust mid-plane layer – in the scenario SIwhileVSI than if only the streaming instability is simulated.

Like the dust-to-gas surface density, the azimuthal gas velocity is most strongly locally enhanced in the scenario SIafterVSI (right panels) as compared with the other two scenarios. While in these scenarios turbulence in the dust layer is predominantly driven by the streaming instability, it is mainly caused by the vertical shear instability in SIafterVSI. Nonetheless, augmentations of the azimuthal gas velocity are greater and more long-lived in this scenario than the augmentations that result from the transient pressure bumps caused by the vertical shear instabiltity. This can be seen from Fig. 6, which shows the deviation of the azimuthal gas velocity from the Keplerian velocity in our simulation of only the vertical shear instability and the simulation of the scenario SIafterVSI with the same initial dust-to-gas surface density ratio and dust size.

That is, local enhancements in the azimuthal gas velocity are stronger in the scenario SIafterVSI than if either only the streaming instability or only the vertical shear instability are considered. In line with what we speculate in SJB20, we conclude that in SIafterVSI dust accumulates in pressure bumps caused by the vertical shear instability, with these accumulations acting as seeds for and being reinforced by the streaming instability. Indeed, in SIafterVSI we show that while turbulence in the dust layer is primarily driven by the vertical shear instability in this scenario, it is induced by the streaming instability locally in dust overdensities.

We do not further discuss our model of only the vertical shear instability, and in particular dust concentration in it, in this paper. This is because simulating only this instability requires neglecting the drag exerted by the dust on the gas since otherwise the streaming instability, too, would be active. Excluding this drag is only justified if the dust density is much less than the gas density. However, this condition cannot be reconciled with our finding that dust concentration is high when the drag of the dust onto the gas is taken into account.

Fig. 7 Maximum (orange lines) as well as dust-mass-weighted average (blue lines) and standard deviation (blue-shaded areas) of dust-to-gas volume density ratio ρd/ρg as a function of time after the dust particle initialisation t – td,init The maximum, mean, and standard deviation are well-correlated, with the Pearson correlation coefficients6 of maximum and mean being given in the top right of the panels. Three phases in the evolution of the volume density ratio in the SIwhileVSI simulation can be distinguished: an initial phase of sedimentation, followed by a quasi-steady phase, and a phase of strong clumping characterised by strong variations. In comparison, the simulation of the streaming instability never evolves into the strong-clumping phase, while the SIafterVSI simulation skips the quasi-steady phase.

4 Metrics to Establish Planetesimal Formation

The formation of planetesimals can be observed in three-dimensional simulations of the streaming instability including dust self-gravity (e.g. Johansen et al. 2007, 2009, 2011). However, because the computational cost of such three-dimensional simulations is prohibitive, two-dimensional models without self-gravity are employed to explore which combinations of dust-to-gas surface density ratio and dust size or Stokes number provide the necessary conditions for the streaming instability to induce planetesimal formation (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021).

These parameter studies – including the one presented in this paper – aim to establish whether dust concentration owing to the streaming instability is sufficiently strong that it would lead to planetesimal formation in equivalent three-dimensional models with self-gravity. Carrera et al. (2015) assume that this is the case if the time-averaged dust surface density distribution deviates sufficiently from a uniform distribution. Yang et al. (2017), on the other hand, investigate whether strong dust clumping can be seen from the spatial distribution of the dust-to-gas volume density ratio. Finally, Li & Youdin (2021) examine whether the maximum dust density exceeds the Roche density.

However, an important caveat to any approach involving the maximum dust-to-gas volume density ratio is that this quantity is inherently stochastic. A measured high density ratio might be a numerical artifact, attained in very few or only a single grid cell. And even if it is physical, it can only be used to establish if planetesimals would potentially form, but does not reveal any information about the properties of the planetesimals. Moreover, when comparing to the Roche density one needs to bear in mind that overdensities which barely exceed the Roche density might only be transient if they are not sufficiently gravitational unstable to overcome turbulent diffusion⁵. We address these caveats below.

4.1 Correlation Between Maximum and Mean Dust-to-gas Density Ratio and Strong Clumping

Figure 7 depicts the evolution of the maximum, mean, and standard deviation of the dust-to-gas volume density ratio, where the mean is weighted by the dust mass, in the same simulations as shown in Figs. 3–4. Similar to in the simulations of the streaming instability presented by Yang et al. (2017) and Li & Youdin (2021), the evolution of the maximum volume density ratio can be divided into three phases: dust settling and formation of an equilibrium mid-plane layer, a subsequent quasi-steady state, and finally a phase of strong clumping. The latter is characterised by greater maxima and large variations that reflect the formation and dissolution of prominent overdensities, as is evident from comparing this figure with the upper panels of Fig. 4.

It is this strong-clumping phase that is associated with the formation of planetesimals (Johansen et al. 2015; Yang et al. 2017; Li & Youdin 2021). Only the simulations of SIwhileVSI and SIafterVSI evolve into such a phase, however, while the simulation of the streaming instability alone remains in an equilibrium state. Indeed, as we show in Sect. 3.2, dust overdensities caused by the vertical shear instability and the streaming instability in combination are larger than the ones induced by the streaming instability in isolation. While in the SIafterVSI simulation sedimentation is directly followed by strong clumping, the SIwhileVSI simulation undergoes a quasi-steady phase until strong clumping sets in when the vertical shear instability starts to affect the dust layer.

There is a strong correlation between the maximum of the dust-to-gas volume density ratio and its mean and standard deviation in all models. This shows that the maximum is in fact representative of the dust clumping behaviour not only in a single or a few cells, but globally in our models that cover the mid-plane layer on the scale of entire protoplanetary disks. As can be seen in Fig. 7, after sedimentation the mean and maximum are nearly constant in the simulation of the streaming instability only, with the Pearson correlation coefficient⁶ of the two amounting to R = 0.819. In the SIafterVSI simulation, on the other hand, the correlation coefficient is even greater with a value of R = 0.888 despite greater variations in mean and maximum. This is because the building and breaking up of over-densities is reflected not only in the maximum, but also in the average. We find this to be generally true in our models since the median correlation coefficient is equal to R = 0.823 for simulations which develop into a strong-clumping phase, compared with R = 0.857 for simulations which do not advance from a quasi-steady state. Simon et al. (2016) find a similar correlation between the maximum dust density and a measure of the dust density dispersion, the ratio of root mean square to mean density.

Fig. 8 Maximum dust-to-gas volume density ratio (ρd/ρg)max as a function of time after the dust particle initialisation t – td,init. In the left, middle, and right panels, respectively, simulations of the streaming instability only, of the scenario SIwhileVSI, and of the scenario SIafterVSI with different initial dust-to-gas surface density ratios (Σd/Σg)init, fixed dust sizes a or Stokes numbers St, and resolutions (fiducial or double) can be seen. Increasing any of the three parameters causes the maximum volume density ratio to increase as well. In addition, the maximum is larger if a dust size of a = 3 cm rather than a Stokes number of 0.1 is considered, as is evident from the left and middle panels. Neither the simulation of the streaming instability with the doubled resolution nor the corresponding one with the fiducial resolution develop into a strong-clumping phase, with the maximum volume density ratio being greater by a factor of a few in the former simulation. In comparison, the maximum is enhanced by more than an order of magnitude during the strong-clumping phase in the double-resolution simulation of SIafterVSI than in respective fiducial-resolution simulation.

4.2 Parameter Study of Maximum Dust-to-gas Density Ratio

We now discuss the dependence of the maximum dust-to-gas volume density ratio on two physical parameters, the dust size and the initial dust-to-gas surface density ratio, as well as one numerical parameter, the simulation resolution. Figure 8 shows the evolution of the maximum volume density ratio in simulations of the streaming instability alone and of the scenarios SIwhileVSI and SIafterVSI with different combinations of these three parameters. Additionally, the average maximum during either the strong-clumping phase for all simulations that evolve into such a phase or during the quasi-steady phase for all other simulations is listed in Table 2.

To begin with, the values of all three parameters being equal, the maximum volume density ratio is higher in the scenario SIafterVSI than in the scenario SIwhileVSI, and lowest in the model of only the streaming instability. This can be gathered from Table 2, and from comparing the maximum volume density ratios in the simulations of the three scenarios with an initial surface density ratio of 1% and a dust size of 3 cm that are depicted as orange lines in Fig. 8. It underlines that the vertical shear instability and the streaming instability in concert cause stronger dust concentration than the streaming instability in isolation.

Furthermore, in agreement with previous studies of the streaming instability (Bai & Stone 2010b; Johansen et al. 2015; Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021), we find that an increase in either the initial surface density ratio or the dust size generally results in a larger maximum volume density ratio. In addition, it is higher if the dust size is fixed at 3 cm or 3 mm, corresponding to Stokes numbers ranging from ~0.05 or ~0.005 at r = 10 au to ~0.5 or ~0.05 at r = 100 au in the mid-plane of our model (see Fig. 2 and Eq. (12)), than if the Stokes number is fixed at 0.1 or 0.01. We attribute this to dust piling up at small radii in the former case but not in the latter one. The speed of the radial dust drift can be expressed as

$$v_{d,r}=\frac{2St{f}_g^2\Pi c_s}{1+S{t}^2{f}_g^2},$$(14)

where f_g is the local ratio of gas density to total density of dust and gas (Nakagawa et al. 1986). Since Π(z = 0) ∝ r^1/4 (see Eq. (8)) and c_s ∝ r^−1/4 in our model, the drift speed is independent of the radius in our simulations with a fixed Stokes number (if f_g is radially constant). On the other hand, in simulations with a fixed dust size the Stokes number increases with the radius, and thus also the drift speed.

An enhanced resolution also entails an increase in the maximum volume density ratio, both during the quasi-steady phase and while strong clumping occurs. The former is demonstrated by our model of the streaming instability with an initial surface density ratio of 1% and a Stokes number of 0.1 which is shown in the left panel of Fig. 8, with the maximum volume density ratio being greater by a factor of a few if the resolution is twice as high. Similar increases are found by Johansen et al. (2015, see their Fig. 1) and Yang et al. (2017). On the other hand, doubling the resolution in the model of the scenario SIafterVSI with the same initial surface density ratio but a dust size of 3 mm that is depicted in the right panel of the figure results in an enhancement of the maximum volume density ratio by more than an order of magnitude during the strong-clumping phase. This is again in agreement with the enhancement in previously presented models of the streaming instability (Yang & Johansen 2014; Johansen et al. 2015; Yang et al. 2017).

Neither our simulation of the streaming instability with the fiducial resolution nor the one with doubled resolution develops into a strong-clumping phase, so the increased resolution does not lead to strong clumping by itself. While this is consistent with the findings by Li & Youdin (2021), Yang & Johansen (2014) and Yang et al. (2017) show that enhancing the resolution can indeed cause such a transition to strong clumping⁷. It is thus important to keep in mind that whether or not strong clumping and potentially planetesimal formation arise in models of the streaming instability (and the vertical shear instability) can be dependent on the numerical resolution.

4.3 Mid-plane Dust-to-gas Density Ratio

As noted in Sect. 3.2, dust overdensities are spread out over the entire vertical extent of the dust layer. Therefore, even though local maxima in the dust-to-gas volume density ratio are of the order of ten or a hundred (see Fig. 7), the mid-plane density ratio in general remains less than one. Figure 9 shows the evolution of the mid-plane density ratio, averaged over the range of radii depicted in Fig. 3, in the same simulations as can be seen in that figure as well as in the corresponding simulations with an initial dust-to-gas surface density ratio of 2%.

We find the mid-plane density ratio not to be indicative of dust concentration in the scenario SIafterVSI. It is lowest in this scenario, whereas dust overdensities are greatest. In addition, the fluctuations of the mid-plane density ratio are likely the result of oscillations of the wave-shaped dust layer rather than the building up and breaking up of dust overdensities. The mid-plane density ratio increases with the dust-to-gas surface density ratio because dust-induced buoyancy increasingly suppresses the vertical shear instability (Lin 2019; SJB20), which is the main driver of turbulence in the dust layer in this scenario.

Compared with the scenario SIafterVSI, the dust scale height is smaller (see Sect. 3.2) and the mid-plane density ratio therefore higher in the model of the streaming instability in isolation and in the scenario SIwhileVSI. It is comparable in these two models since in both of them the turbulence in the dust layer is predominantly caused by the streaming instability.

Nonetheless, as in the case of the scenario SIafterVSI, the mid-plane density ratio does not reflect dust concentration in the model of the streaming instability alone and in the scenario SIwhileVSI. This can be gathered from the fact that the mid-plane density ratio reaches its maximum early owing to sedimentation and remains close to constant at a lower value afterwards, while maxima in the dust-to-gas surface density ratio (see Fig. 4) as well as in the maximum dust-to-gas volume density ratio (see Fig. 8) are attained later. The simulations of the streaming instability presented by Flock & Mignone (2021, see their Figs. 2–3) exhibit a similar discrepancy between the evolution of the mid-plane dust-to-gas mass ratio and the maximum dust-to-gas mass ratio as well as the dust surface density.

Furthermore, the mid-plane density ratio is almost exactly twice as high if the initial surface density ratio is doubled from 1% to 2%, while the maximum volume density ratio increases non-linearly for these surface density ratios (see Fig. 8; Johansen et al. 2009, 2015). In the streaming instability simulations by Li & Youdin (2021), the mid-plane density ratio as well increases only approximately linearly with the surface density ratio, while the maximum volume density ratio can be enhanced by orders of magnitude.

Fig. 9 Dust-to-gas density ratio ρd/ρg in the mid-plane as a function of time after the dust particle initialisation t – td,init. The mid-plane density ratio is computed as the mean over r = 55 au to 65 au and the rolling average over 100 yr. Simulations of the streaming instability only as well as the scenarios SIwhileVSI and SIafterVSI with initial dust-to-gas surface density ratios of 1% (dashed lines) or 2% (solid lines) and a Stokes number of 0.1 are shown. The mid-plane density ratio is similar in the simulations of the streaming instability and the scenario SIwhileVSI, and lower in the scenario SIafterVSI. While the ratio fluctuates throughout the latter simulation, it initially attains a maximum before evolving into a quasi-steady state in the former two. In all three cases, doubling the initial surface density ratio entails an increase by about a factor of two also in the mid-plane density ratio.

4.4 Fraction of Dust Mass in Roche-unstable Overdensities

In Sect. 4.2, we show that the maximum dust-to-gas volume density ratio is higher in models which undergo a strong-clumping phase than in models which do not. However, it remains to be investigated whether dust overdensities could undergo gravitational collapse and form planetesimals if self-gravity were included in either or both kinds of models. To this end, in Fig. 10 we show the fraction of the total dust mass that is associated with overdensities of at least once or twice the Roche density. Our simulations of only the streaming instability and of the scenarios SIwhileVSI and SIafterVSI that can be seen also in Figs. 3, 4 and 7 as well as a simulation of the streaming instability with the same initial dust-to-gas surface density ratio and Stokes number but with doubled resolution are depicted.

The figure demonstrates that overdensities are robustly gravitationally unstable in models that experience strong clumping, but not in ones that remain in a quasi-steady state. During the strong-clumping phase in the simulations of the scenarios SIafterVSI and SIwhileVSI, Roche-unstable overdensities comprise up to 1% of the dust mass, and overdensities which are twice as dense a fraction that is only marginally less. The formation and dissolution of these overdensities induces variations of the fraction that are as large as an order of magnitude or more especially in the SIafterVSI simulation, though.

In comparison, a considerably smaller fraction of about 0.01% and 0.1%, respectively, of the dust mass is part of overdensities that barely exceed the Roche density in the streaming-instability-only simulations with the fiducial resolution and with the doubled resolution, which both do not undergo strong clumping (see Fig. 8). On top of that, the fraction that is contained in overdensities of at least twice the Roche density is less by an order of magnitude or more in these two simulations.

4.5 Planetesimal Formation Rate

Assuming that every overdensity that exceeds the Roche density collapses in its entirety to form one or multiple planetesimals, the fraction of the total dust mass that is part of such overden-sities at a given time is equivalent to a momentary planetesimal formation efficiency. However, since gravitational collapse and planetesimal formation are not actually included in our model, each dust particle can be part of multiple such overdensities during the course of the simulations.

We therefore randomly select 10000 dust particles at the start of the strong-clumping phase in each simulation which develops into such a phase, and for each of these particles track how much time passes until the dust density at its momentary location for the first time is greater than the Roche density. Figure 11 shows the cumulative fraction of the total dust mass that has been part of a Roche-unstable overdensity in our simulations of the streaming instability only as well as the scenarios SIwhileVSI and SIafterVSI with an initial dust-to-gas surface density ratio of 2% and a Stokes number of 0.1. Since we generally find a roughly linear increase of this cumulative fraction with time, we compute average fractions that become associated with a Roche-unstable overdensity per unit time. We refer to these as planetesimal formation rate, with a rate of 1 kyr⁻¹ indicating that all particle mass would be converted to planetesimal mass within 1 kyr. For all simulations in which strong clumping occurs, we list these planetesimal formation rates in Table 2.

Analogous to the trends we describe for the maximum dust-to-gas volume density ratio in Sect. 4.2, we find these rates to be higher if any of the initial surface density ratio, the initial dust size or Stokes number, and the resolution are greater. There is one notable difference, though, which can be gathered both from the figure and the table: While the maximum volume density ratio is greatest in the scenario SIafterVSI, the rates are by tendency largest in the scenario SIwhileVSI – the overall highest rate is attained in a simulation of the streaming instability alone, though. The reason for this discrepancy probably lies in dust diffusion being stronger in the scenario SIafterVSI, where the vertical shear instability is the main source of turbulence in the dust layer, than in the scenario SIwhileVSI or in our model of the streaming instability in isolation, where it is (predominantly) caused by this instability. Overdensities are thus more prone to dispersing and forming anew. This explanation is supported by a comparison of the upper panels of Fig. 4 as well as both by the fluctuations in the fraction of the total dust mass that is comprised in Roche-unstable overdensities being greatest in the scenario SIafterVSI, as can be seen from Fig. 10, and by the standard deviations of the planetesimal formation rates being comparatively large in this scenario.

Fig. 10 Fraction of total dust mass fd that is located where the dust density ρd is greater than a threshold density ρthres as a function of time after the dust particle initialisation t – td,init. The fraction is calculated as the rolling average over 100 yr. In addition to the simulations shown in Figs. 4 and 7, an analogous simulation of the streaming instability with doubled resolution is depicted. We choose threshold densities equal to once (solid lines) and twice the Roche density (dashed lines). In both the SIwhileVSI and the SlafterVSI simulation, Roche-unstable overdensities temporarily comprise as much as 1% of the dust mass, and overdensities that exceed twice the Roche density a slightly smaller fraction. Nevertheless, the fraction fluctuates by an order of magnitude or more over time. In contrast, in the simulations of the streaming instability with the doubled and the fiducial resolution only 0.1% and 0.01%, respectively, of the dust mass is part of overdensities greater than the Roche density, and a fraction that is less by at least an order of magnitude is contained in overdensities of twice the Roche density or more.

Fig. 11 Cumulative fraction of the total dust mass fd,tot that has encountered a location where the dust density ρd exceeds the Roche density ρR as a function of the time t after the start of the strong-clumping phase td,clump. Simulations of the streaming instability in isolation (blue lines) as well as of the scenarios SIwhileVSI (orange lines) and SlafterVSI (green lines) with an initial dust-to-gas surface density ratio of 2% and a Stokes number of 0.1 are depicted. The cumulative fractions (solid lines) increase close to linearly (dashed lines) with time, with the increase being strongest in the scenario SIwhileVSI and weakest if only the streaming instability is taken into account.

5 Thresholds for Planetesimal Formation

Based on what we discussed in the previous sections, we conclude that if strong clumping occurs in a simulation of ours, then planetesimals would form in the simulation if self-gravity were included. This is since overdensities in simulations that undergo strong clumping robustly exceed the Roche density. Whether a simulation develops into a strong-clumping phase is evident, for instance, from the evolution of the maximum dust-to-gas volume density ratio. The maximum is well-correlated with the mean and standard deviation of the volume density ratio and thus representative of dust concentration at large. From Fig. 12 and Table 2, it can be gathered in which of our simulations strong clumping and thus potentially planetesimal formation occur.

Overall, we find that the vertical shear instability and the streaming instability together cause stronger dust concentration than the streaming instability in isolation. This is reflected in these two instabilities inducing strong clumping for combinations of (initial dust-to-gas surface density ratio, dust size or Stokes number) as small as (0.5%, St = 0.1) or (1%, a = 3 mm) in the scenario SIafierVSI, while in the scenario SIwhileVSI at least (1%, St = 0.1) is required. In comparison, the thresholds are higher when only the streaming instability is considered, with (1%, a = 3 cm) or (2%, St = 0.1) being necessary.

In agreement with our findings, planetesimal formation arises for (1%, St ≈ 0.1), but not for (1%, St ≈ 0.01), in the models of the streaming instability including a prescribed pressure bump that are presented by Carrera et al. (2021) and Carrera & Simon (2022). On the other hand, Flock & Mignone (2021) and Li & Youdin (2021) find the streaming instability to lead to planetesimal formation for (1%, St = 0.1) also in the absence of a pressure bump.

In addition to a differentiation between simulations which do or do not experience strong clumping, Fig. 12 depicts trends in the planetesimal formation rate. Half-filled circles represent simulations in which the rate does not exceed 0.1 kyr⁻¹, that is to say after 1 kyr less than 10% of the total dust mass would have been converted to planetesimal mass, while filled circles represent simulations with a rate equal to or greater than this threshold value. As we discuss in the previous section, the rates are generally highest in the scenario SIwhileVSI and exceed the threshold in all cases in this scenario. In the scenario SIafterVSI, on the other hand, the rates are greater than the threshold only in simulations with dust size of 3 cm, and in our model of the streaming instability in isolation only for this dust size and an initial dust-to-gas surface density ratio of 2%.

Fig. 12 Thresholds for strong clumping as well as planetesimal formation rates. Initial dust-to-gas surface density ratios (Σd/Σg)init and Stokes numbers St or dust sizes a are plotted on the ordinate and the abscissa, respectively. Blue, orange, and green colors represent our models of the streaming instability in isolation, the scenario SIwhileVSI and the scenario SIafterVSI. Simulations that do not experience strong clumping are depicted as empty circles; ones in which the planetesimal formation rate during the strong-clumping phase is less than 0.1 kyr−1 are represented by half-filled circles and ones in which the rate is greater than this value as filled circles. Lines separate combinations of (surface density ratio, dust size or Stokes number) for which strong clumping does or does not occur. The minimum values of these two parameters required for strong clumping are smallest in the scenario SIafterVSI, with the combinations (0.5%, St = 0.1) and (1%, a = 3 mm) being sufficient. The planetesimal formation rates, on the other hand, are by tendency highest in the scenario SIwhileVSI.

6 Discussion

6.1 General Implications

We study dust concentration and the potential for planetesimal formation in two-dimensional global simulations of the streaming instability and the vertical shear instability. These simulations represent three different scenarios, two in which the two instabilities coexist and one involving only the streaming instability, and include different initial dust-to-gas surface density ratios, dust sizes or Stokes numbers, respectively, and resolutions.

In the scenario SIafterVSI, the vertical shear instability attains a saturated state before the streaming instability begins to grow, while in the scenario SIwhileVSI both instabilities start to develop at the same time. Both scenarios are plausible since the growth from micron-sized dust grains to the millimetre- or centimetre-sized dust aggregates that we simulated takes thousands of orbital periods (Zsom et al. 2010; Lorek et al. 2018), but it is unclear at which point during this growth protoplanetary disks have evolved into a state similar to the one in our model which is favourable to the development of the vertical shear instability.

Our simulations show that the vertical shear instability causes the formation of gas pressure bumps that promote dust concentration by the streaming instability. This is consistent with both the findings that the vertical shear instability induces dust concentration in pressure bumps (Stoll & Kley 2016) and that pressure bumps - be they prescribed (Lenz et al. 2019; Carrera et al. 2021, 2022; Carrera & Simon 2022; Lehmann & Lin 2022; Xu & Bai 2022b) or caused by the magnetorotational instability under the assumption of ideal magnetohydrodynamics (Johansen et al. 2007, 2011) – facilitate dust accumulation and planetesimal formation owing to the streaming instability. The same has been found for vortices induced by the vertical shear instability (Lehmann & Lin 2022) and by the subcritical baroclinic instability (Raettig et al. 2015, 2021).

We note that Xu & Bai (2022a) present simulations of the magnetorotational instability including ambipolar diffusion in which the instability gives rise to pressure bumps where dust accumulates. In contrast to in our model, though, they find that this dust accumulation is a consequence of the mutual drag between gas and dust, but not of the streaming instability, because it occurs in pressure maxima where there is no pressure gradient to drive the streaming instability (but see Auffinger & Laibe 2018; Lin & Hsu 2022).

A key result of our study is that dust overdensities are sufficient to form planetesimals for lower dust-to-gas surface density ratios and smaller dust sizes if both the vertical shear instability and the streaming instability are considered than when only the streaming instability is taken into account. Both our and previous work (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021) has shown that either dust-to-gas ratios that are greater than the canonical interstellar medium value of 1% or dust sizes which are larger than what is observed in protoplanetary disks are required for the streaming instability in isolation to lead to planetesimal formation.

On the other hand, in our scenario SIafterVSI – in which the vertical shear instability has saturated before the streaming instability begins to grow – we find planetesimal formation to be possible for a surface density ratio of 1% and a dust size of 3 mm, in agreement with the sizes derived from observed opacity spectral indices (e.g. Sierra et al. 2019; Macías et al. 2019, Macías et al. 2021; Carrasco-González et al. 2019), though larger than the ones inferred from polarisation measurements (e.g. Ohashi et al. 2020; Mori & Kataoka 2021). That is, if this scenario applies to some or all protoplanetary disks, planetesimal formation via the vertical shear instability and the streaming instability should be omnipresent in the parts of these disks where the vertical shear instability is active, which roughly correspond to the region between 10 au and 100 au covered by our simulation domains (Lin & Youdin 2015; Malygin et al. 2017; Pfeil & Klahr 2019).

At first glance, our results seem to contradict previous studies finding that turbulence reduces the growth rate of the linear streaming instability (Umurhan et al. 2020; Chen & Lin 2020) and that driven Kolmogorov-like turbulence inhibits planetesimal formation owing to the non-linear instability (Gole et al. 2020). We note, though, that in these studies turbulence is purely a source of isotropic diffusion and viscosity, while our and the above-mentioned work evinces that instabilities which drive turbulence also give rise to pressure bumps and vortices. In addition, turbulence is not generally isotropic, and purely vertical diffusion is not detrimental to radial dust concentration (Yang et al. 2018; SJB20).

6.2 Implications for One-dimensional Models Including Planetesimal Formation

In one-dimensional models of protoplanetary disks that include a prescription of the formation of planetesimals via the streaming instability, it is often assumed that the mid-plane dust-to-gas density ratio needs to exceed unity for planetesimal formation to occur (e.g. Drążkowska & Dullemond 2014; Drążkowska et al. 2016; Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017; Schoonenberg et al. 2018; Stammler et al. 2019). This condition is based on the simulations by Johansen & Youdin (2007), which notably do not include the vertical stellar gravity.

However, while this condition might be sufficient, our model shows that it is not necessary. We find that dust concentration is sufficiently strong for planetesimal formation in all of our simulations with an initial dust-to-gas surface density ratio of 2% and a Stokes number of 0.1, particularly also in the one of the streaming instability alone (see Fig. 12). Nevertheless, it can be seen from Fig. 9 that the mid-plane density ratio remains less than one in these simulations.

We therefore recommend to base prescriptions of planetesimal formation owing to the streaming instability on whether the dust-to-gas surface density ratio and the dust size exceed the threshold values given by our Fig. 12 if both the streaming instability and the vertical shear instability are modeled; and by this figure, by Eqs. (8) and 9 of Yang et al. (2017), or by Eq. (10) of Li & Youdin (2021) if only the streaming instability is taken into consideration. If planetesimal formation is established to occur, our Table 2 lists planetesimal formation rates, rates at which dust mass is converted to planetesimal mass, for a given combination of surface density ratio and dust size in each of our three scenarios.

6.3 Limitations

While we investigate which conditions, in terms of dust-to-gas surface density ratio and dust size, are necessary for the vertical shear instability and streaming instability in combination or the streaming instability in isolation to induce planetesimal formation, we cannot actually model this process since our two-dimensional simulations do not include self-gravity. As with previous similar parameter studies (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021), this is because it is not computationally feasible to cover a significant fraction of the parameter space with three-dimensional simulations.

Nonetheless, we discuss in detail in Sect. 4 which metrics are applicable to two-dimensional models to gauge the potential for planetesimal formation in equivalent three-dimensional ones. Furthermore, it has been shown that the surface density ratios and dust sizes that are necessary for dust concentration to be sufficient for planetesimal formation are comparable in two-and three-dimensional simulations of the streaming instability (Yang et al. 2017; Li & Youdin 2021) and vertical shear instability (Lehmann & Lin 2022). Nevertheless, the non-linear regimes of the streaming instability (Kowalik et al. 2013) and the vertical shear instability (Nelson et al. 2013; Stoll & Kley 2014), in contrast to their linear regimes (Youdin & Goodman 2005; Nelson et al. 2013; Barker & Latter 2015), are not axisymmetric. In particular, the vertical shear instability gives rise to vortices in which dust accumulates (Flock et al. 2020; Lehmann & Lin 2022).

More generally, our model could be improved upon by deviating from the assumption of a constant initial surface density ratio and dust size or Stokes number, and instead considering the structures of rings and gaps with varying dust and gas surface densities and maximum dust sizes that are observed in protoplanetary disks (e.g. Macías et al. 2019, 2021; Carrasco-González et al. 2019; Andrews 2020). This could entail a need for other parameters than the surface density ratio and dust size, for instance the dust flux (Lenz et al. 2019; Flock & Mignone 2021), to describe the conditions that are required for planetesimal formation to occur. We further do not consider dust size distributions (Bai & Stone 2010b; Schaffer et al. 2018, Schaffer et al. 2021; Krapp et al. 2019; Zhu & Yang 2021; Yang & Zhu 2021). In addition, the gas disk is in our model is only affected by the stellar gravity as well as the vertical shear instability and the streaming instability, processes like disk winds and other instabilities are not taken into account; we simulated a simple isothermal or adiabatic equation of state rather than the various disk heating and cooling processes; and we do not consider dust and gas chemistry.

7 Summary

We employ two-dimensional, axisymmetric adaptive mesh refinement simulations of the vertical shear instability and the streaming instability, which cover the outer regions of protoplanetary disks on a global scale, to identify the threshold values of dust-to-gas surface density ratio and dust size or Stokes number that are required for the two instabilities in conjunction or the streaming instability in isolation to induce planetesimal formation. Similar parameter studies are presented by Carrera et al. (2015), Yang et al. (2017), and Li & Youdin (2021), though all of these are based on two-dimensional local shearing box simulations of the streaming instability only.

Since self-gravity is not included in our two-dimensional simulations, planetesimal formation is not actually modeled. We therefore dedicate a detailed discussion to metrics that can be applied to differentiate between simulations in which it would occur and ones in which it would not:

• The maximum dust-to-gas volume density ratio is well-correlated with the mean and standard deviation of the volume density ratio and thus representative of the strength of dust concentration at large.

• The maximum volume density ratio increases if either the initial dust-to-gas surface density ratio, the dust size, or the resolution is enhanced. This is in agreement with previous work on the streaming instability (Bai & Stone 2010b; Yang & Johansen 2014; Johansen et al. 2015; Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021). Moreover, it is higher in our models with a fixed dust size than in ones with a comparable fixed Stokes number. This is because the drift speed increases with the radial distance to the star in the former case (but not in the latter), resulting in dust piling up in the radial dimension.

• While some simulations remain in a quasi-steady state after the dust has settled to a mid-plane layer, others evolve into a phase of strong dust clumping (Johansen et al. 2015; Yang et al. 2017; Li & Youdin 2021). This phase is characterised by comparatively high maxima and large fluctuations in the maximum volume density ratio, which are a consequence of the formation and dissolution of strong overdensities. While whether or not strong clumping occurs does not depend on the resolution in our study, a transition to strong clumping is found in previous studies of the streaming instability if the resolution is increased (Yang & Johansen 2014; Yang et al. 2017).

• Only in simulations that undergo a strong-clumping phase are overdensities robustly Roche-unstable, that is to say only in these simulations is a significant fraction of the total dust mass comprised in overdensities that exceed once and even twice the Roche density.

• The mid-plane dust-to-gas density ratio is not suitable as an indicator of local dust concentration because local maxima of the volume density ratio can be found at all heights within the dust layer.

This leads us to the conclusion that models in which planetesimal formation would be possible if self-gravity were taken into account distinguish themselves from ones in which it would not in that they develop into a strong-clumping phase. From Fig. 12 and Table 2, it can be gathered which of our simulations do or do not experience strong clumping. Consistent with the trends described above for the maximum volume density ratio, these are simulations with high initial surface density ratios and large dust sizes or Stokes numbers.

Importantly, we find that the minimum surface density ratios and dust sizes or Stokes numbers which are necessary for strong clumping and thus potentially planetesimal formation are lower in our models of both the vertical shear instability and the streaming instability than in our model of only the latter. They are smallest if the vertical shear instability has saturated before the streaming instability begins its growth, and higher if both instabilities start their growth at the same time. The reason for this lies in the vertical shear instability giving rise to pressure bumps in which dust accumulates, with these accumulations in turn seeding further dust concentration by the streaming instability.

Since the cumulative fraction of the dust mass that has been part of a Roche-unstable overdensity increases largely linearly with time in our simulations, we calculate planetesimal formation rates. These rates are listed in Table 2 and indicate which fraction of the dust mass becomes associated with such an over-density and would be converted to planetesimal mass per unit time. They exceed 0.1 kyr⁻¹ in several cases and 0.25 kyr⁻¹ in the best case, and are again by tendency greater in our models of both instabilities in concert than in our model of the streaming instability in isolation, though largest if both instabilities grow simultaneously.

Acknowledgements

We are thankful to the anonymous referee for their comments that helped to improve in particular the clarity of this paper. To analyse and visualise the simulations presented in this paper, the Python packages yt (http://yt-project.org; Turk et al. 2011), Matplotlib (https://matplotlib.org; Hunter 2007), and NumPy (https://numpy.org; Oliphant 2006) were used. The FLASH Code has in part been developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. Computational resources employed to conduct the simulations were provided by the Regionales Rechenzentrum at the University of Hamburg, by the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN), and by the SCIENCE HPC Center at the University of Copenhagen. U.S. and A.J. are thankful for funding from the European Research Foundation (ERC Consolidator Grant 724687-PLANETESYS). A.J. further gratefully acknowledges funding from the Knut and Alice Wallenberg Foundation (Wallenberg Scholar Grant 2019.0442), the Swedish Research Council (Project Grant 2018-04867), the Danish National Research Foundation (DNRF Chair Grant DNRF159), and the Göran Gustafsson Foundation.

References

All Tables

All Figures

Fig. 1 Gas density ρg in the mid-plane (at a height z = 0) as a function of the radius r. The blue line shows the initial density profile, while the orange and green line depict the profiles after 50 kyr in simulations including the static mesh refinement employed in SJB20 and in this study, respectively. Notably, a considerable bump at a radius of ~30 au is present in the former profile, but not in the latter one.

In the text

	Fig. 2 Dust size a (left ordinate) and Stokes number St (right ordinate) as functions of the radius r. Conversion relations between dust size and Stokes number in the mid-plane of our protoplanetary disk model (Eqs. (12)–(13)) are shown for the two dust sizes (orange lines) and the two Stokes numbers (blue lines) that we simulated, assuming an isothermal equation of state.
In the text

Fig. 3 Ratio of dust ρd to gas volume density ρg as a function of radius r and height z, the latter given in units of gas scale heights Hg (left panels). The average over the radius range shown in the left panels is depicted as a function of height in the right panels. The top, middle, and lower panels show the dust-to-gas volume density ratio at the end of simulations of the streaming instability in isolation, the scenario SIwhileVSI, and the scenario SIafterVSI, respectively, with an initial dust-to-gas surface density ratio of 1% and a Stokes number of 0.1. The horizontal white lines in the lower panel illustrate the range of heights depicted in the upper and middle panels. In all simulations, local maxima of the volume density ratio are distributed over the entire vertical dimension of the dust layer. While the vertical dust distribution in the simulations of the streaming instability only and of the scenario SIwhileVSI is largely Gaussian, the dust layer possesses the form of a wave rather than a Gaussian shape in the SIafterVSI simulation.

In the text

Fig. 4 Ratio of dust Σd to gas surface density Σg (upper panels) and difference between azimuthal gas velocity υg,ϕ and Keplerian velocity υK (lower panels). Both quantities are depicted as a function of radius r and time t after the initialisation of the dust particles at td,init. The velocity difference is averaged over the vertical domain size, weighted by the dust mass, and expressed as a Mach number. In the absence of turbulence, this quantity is equal to the negative of the parameter Π as defined in Eq. (8). Black dots indicate that the azimuthal velocity is equal to or greater than the Keplerian velocity. Enhancements in the dust-to-gas surface density ratio and in the difference between azimuthal gas velocity and Keplerian velocity are larger at late times (t – td,init > 5 kyr) and small radii (r < 40 au) in SIwhileVSI than in the simulation of the streaming instability in isolation, and greatest at all times and radii in SIafterVSI. Surface density ratio and velocity difference can be seen to generally be correlated.

In the text

Fig. 5 Dust-to-gas volume density ratio ρd/ρg as a function of radius r and height z in simulations of the scenarios SIwhileVSI (top panel) and SIafterVSI (bottom panel) with an initial dust-to-gas surface density ratio of 2% and a dust size of 3 cm. The volume density ratio is depicted at the end of the simulations 10 kyr after the initialisation of the dust particles. The vertical dust distribution in the SIwhileVSI simulation resembles a Gaussian, but is wave-like in the SIafterVSI simulation.

In the text

Fig. 6 Dust-mass-weighted mean deviation of azimuthal gas velocity from Keplerian velocity υg,ϕ – υK, normalised by the sound speed cs, as a function of radius r and time after the dust particle initialisation t – td,init. The left and right panels depict simulations of the vertical shear instability only and of the scenario SIafterVSI, respectively. Local increases in the azimuthal gas velocity are stronger, with the velocity more frequently exceeding the Keplerian velocity, and especially at r > 30 au also more persistent in SIafterVSI.

In the text

Fig. 7 Maximum (orange lines) as well as dust-mass-weighted average (blue lines) and standard deviation (blue-shaded areas) of dust-to-gas volume density ratio ρd/ρg as a function of time after the dust particle initialisation t – td,init The maximum, mean, and standard deviation are well-correlated, with the Pearson correlation coefficients6 of maximum and mean being given in the top right of the panels. Three phases in the evolution of the volume density ratio in the SIwhileVSI simulation can be distinguished: an initial phase of sedimentation, followed by a quasi-steady phase, and a phase of strong clumping characterised by strong variations. In comparison, the simulation of the streaming instability never evolves into the strong-clumping phase, while the SIafterVSI simulation skips the quasi-steady phase.

In the text

Fig. 8 Maximum dust-to-gas volume density ratio (ρd/ρg)max as a function of time after the dust particle initialisation t – td,init. In the left, middle, and right panels, respectively, simulations of the streaming instability only, of the scenario SIwhileVSI, and of the scenario SIafterVSI with different initial dust-to-gas surface density ratios (Σd/Σg)init, fixed dust sizes a or Stokes numbers St, and resolutions (fiducial or double) can be seen. Increasing any of the three parameters causes the maximum volume density ratio to increase as well. In addition, the maximum is larger if a dust size of a = 3 cm rather than a Stokes number of 0.1 is considered, as is evident from the left and middle panels. Neither the simulation of the streaming instability with the doubled resolution nor the corresponding one with the fiducial resolution develop into a strong-clumping phase, with the maximum volume density ratio being greater by a factor of a few in the former simulation. In comparison, the maximum is enhanced by more than an order of magnitude during the strong-clumping phase in the double-resolution simulation of SIafterVSI than in respective fiducial-resolution simulation.

In the text

Fig. 9 Dust-to-gas density ratio ρd/ρg in the mid-plane as a function of time after the dust particle initialisation t – td,init. The mid-plane density ratio is computed as the mean over r = 55 au to 65 au and the rolling average over 100 yr. Simulations of the streaming instability only as well as the scenarios SIwhileVSI and SIafterVSI with initial dust-to-gas surface density ratios of 1% (dashed lines) or 2% (solid lines) and a Stokes number of 0.1 are shown. The mid-plane density ratio is similar in the simulations of the streaming instability and the scenario SIwhileVSI, and lower in the scenario SIafterVSI. While the ratio fluctuates throughout the latter simulation, it initially attains a maximum before evolving into a quasi-steady state in the former two. In all three cases, doubling the initial surface density ratio entails an increase by about a factor of two also in the mid-plane density ratio.

In the text

Fig. 10 Fraction of total dust mass fd that is located where the dust density ρd is greater than a threshold density ρthres as a function of time after the dust particle initialisation t – td,init. The fraction is calculated as the rolling average over 100 yr. In addition to the simulations shown in Figs. 4 and 7, an analogous simulation of the streaming instability with doubled resolution is depicted. We choose threshold densities equal to once (solid lines) and twice the Roche density (dashed lines). In both the SIwhileVSI and the SlafterVSI simulation, Roche-unstable overdensities temporarily comprise as much as 1% of the dust mass, and overdensities that exceed twice the Roche density a slightly smaller fraction. Nevertheless, the fraction fluctuates by an order of magnitude or more over time. In contrast, in the simulations of the streaming instability with the doubled and the fiducial resolution only 0.1% and 0.01%, respectively, of the dust mass is part of overdensities greater than the Roche density, and a fraction that is less by at least an order of magnitude is contained in overdensities of twice the Roche density or more.

In the text

Fig. 11 Cumulative fraction of the total dust mass fd,tot that has encountered a location where the dust density ρd exceeds the Roche density ρR as a function of the time t after the start of the strong-clumping phase td,clump. Simulations of the streaming instability in isolation (blue lines) as well as of the scenarios SIwhileVSI (orange lines) and SlafterVSI (green lines) with an initial dust-to-gas surface density ratio of 2% and a Stokes number of 0.1 are depicted. The cumulative fractions (solid lines) increase close to linearly (dashed lines) with time, with the increase being strongest in the scenario SIwhileVSI and weakest if only the streaming instability is taken into account.

In the text

Fig. 12 Thresholds for strong clumping as well as planetesimal formation rates. Initial dust-to-gas surface density ratios (Σd/Σg)init and Stokes numbers St or dust sizes a are plotted on the ordinate and the abscissa, respectively. Blue, orange, and green colors represent our models of the streaming instability in isolation, the scenario SIwhileVSI and the scenario SIafterVSI. Simulations that do not experience strong clumping are depicted as empty circles; ones in which the planetesimal formation rate during the strong-clumping phase is less than 0.1 kyr−1 are represented by half-filled circles and ones in which the rate is greater than this value as filled circles. Lines separate combinations of (surface density ratio, dust size or Stokes number) for which strong clumping does or does not occur. The minimum values of these two parameters required for strong clumping are smallest in the scenario SIafterVSI, with the combinations (0.5%, St = 0.1) and (1%, a = 3 mm) being sufficient. The planetesimal formation rates, on the other hand, are by tendency highest in the scenario SIwhileVSI.

In the text