© ESO 2021

1 Introduction

The initial stages of the planet formation process are hindered by various barriers. One important hurdle is the radial drift barrier, where solids of millimeter – centimeter – sizes drift as rapidly toward the star as they grow, limiting the overall particle size that can be achieved before reaching the star (Birnstiel & Andrews 2014). A key process in aiding solid concentration is the streaming instability (Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007; Capelo et al. 2019; Schneider et al. 2019). When the back-reaction of the protoplanetary disk solids on the gas is taken into account, local and nearly axisymmetric particle filaments can form. Once these reach the Roche density, their gravitational collapse produces planetesimals of 50 to a few 100 km in size (Johansen et al. 2015; Simon et al. 2016; Schäfer et al. 2017; Klahr & Schreiber 2021).

Protoplanetary disk observations and laboratory experiments that aim to recreate conditions in such settings confirm that solids with a non-monodisperse size distribution coexist in these systems. ALMA observations of the rings in the protoplanetary disk around HL Tau suggest the presence of up to centimeter – sized solids (ALMA Partnership 2015; Zhang et al. 2015). In the case of the TW Hya system, which is one of the closest observed systems with a protoplanetary disk, Menu et al. (2014) study a large set of observational data covering a wide range of near-infrared to millimeter – and centimeter – wavelengths, together with spectral energy distribution data. They show that models that reproduce the observations include grains that range from microns to centimeters in size. Laboratory experiments that study interactions between porous aggregates also show that the outcome of collisions is grains of varying sizes (Güttler et al. 2010; Zsom et al. 2010).

In line with what is mentioned above, many streaming instability models in recent years have included multiple particle species. Krapp et al. (2019) and Paardekooper et al. (2020) find that the linear phase of the instability has longer growth timescales than seen in mono-disperse models. However, McNally et al. (2021), using the methodology described in Paardekooper et al. (2021), expand the parameter space from their nominal set described in Paardekooper et al. (2020) to a larger range of particle sizes and size distributions, which are more in line with those predicted by dust coagulation – fragmentation models (Birnstiel et al. 2011, 2015). They see that a fast growth regime is achieved given that the particle size distribution has a sufficient fraction of large particles, peaks at a friction time of 0.1Ω⁻¹ and has a local dust-to-gas ratio above unity. Similarly, extending the parameter space, Zhu & Yang (2021) also see that given an increased maximum solid size or total dust-to-gas ratio, the multispecies streaming instability is indeed favorable for solid growth. The nonlinear phases of the instability were studied in Bai & Stone (2010a) and Schaffer et al. (2018). Given super-solar metallicity, Bai & Stone (2010a) see the formation of dense filament clumps and in Bai & Stone (2010b), they show that one particle per grid cell is sufficient for numerical convergence.

In this paper, we build on the results presented in Schaffer et al. (2018). We study how increasing the particle number, considering various particle size distributions, and varying the particle size range affects the efficiency of the streaming instability. Our goal is to find the critical metallicities where particle filaments begin to form in order to compare how the previously mentioned parameters influence particle clumping through the streaming instability. We aim to understand if the particle number plays as strong a role in the efficiency of the multispeciesinstability in the stratified case as that what is reported by Krapp et al. (2019) for the unstratified case.

2 Numerical method

Using the Pencil Code¹ we performed two-dimensional shearing box simulations that represent a local segment of a protoplanetary disk. Our system corotates the central star at a distance r with a Keplerian velocity of v_K. Our coordinate axes are oriented with x pointing radially outward, y along direction of the disk rotation and z vertically out of the disk. We neglected both magnetic effects and the self-gravity of the solids. The equations that drive the gas component are the momentum and continuity equations, which are presented in Eqs. (1) and (2) of Schaffer et al. (2018).

Here, the gas density and velocity are defined as ρ_g and u, respectively. The local Keplerian angular frequency is Ω and the dust-to-gas ratio is defined as ε = ρ_p∕ρ_g, with ρ_p denoting the particle density. The dimensionless parameter that sets the large-scale pressure gradient of the disk is defined as $$\eta =-\frac{1}{2}{(\frac{H_{\text{g}}}{r})}^2\frac{\partial \text{ln }P}{\partial \text{ln }r}.$$(1)

As in Schaffer et al. (2018), we non-dimensionalized the pressure gradient as Π = ηv_K∕c_s, where c_s is the sound speed, and chose Π = 0.05, as it represents the pressure support in the inner protoplanetary disk well (Bai & Stone 2010a). In Eq. (1), H_g is the gas scale height and P the pressure. Finally, τ_f is the friction time of the solid component, which in the Epstein regime is defined as $${\tau }_{\text{f}}=\frac{{\rho }_{•}a}{{\rho }_{\text{g}}{c}_{\text{s}}},$$(2)

where ρ_∙ is the internal solid densityand a the solid radius. Using Ω to non-dimensionalize τ_f gives the Stokes number, which becomes $${\tau }_{\text{s}}=\sqrt{2\pi }\frac{{\rho }_{•}a}{{\Sigma }_{\text{g}}},$$(3)

where Σ_g is the gas surface density. In the rest of paper, we use the Stokes number to represent particle sizes.

The solid component is modeled using super-particles, which represent swarms of physical particles. These follow the equations presented in Eqs. (4) and (5) of Schaffer et al. (2018). The individual particles have velocity v_i and position x_i. The velocitydispersion of the superparticles is captured as the distribution of velocities between the superparticles present within each grid cell. The back-reaction friction force from the particles onto the gas is calculated by mapping the particle density ρ_p onto the grid using a TSC (Triangular Shaped Cloud) scheme (Youdin & Johansen 2007).

The initial conditions for the gas and particle density profile both follow a stratified profile centered around the midplane. The nominal initial particle scale height is H_p,init = 0.015H_g. In order to ensure that the initial particle scale height does not influence the results, we also performed a run with H_p,init = 1H_g. In this case, the initial particle profile becomes unstratified within our simulation domain. The gas component initially has a sub-Keplerian velocity, such that $u(t=0)=-\eta v_{\text{K}}\widehat{y}$ and the initial particle velocity is zero.

We fixed the simulation domain at 0.2H_g × 0.2H_g with a resolution of 128 × 128 grid cells. The boundary conditions in the vertical direction are reflecting. Some small particles hit the vertical boundary and become reflected, however following our findings in Schaffer et al. (2018), we do not expect this to have a significant effect on the results. The boundary conditions in the radial direction are periodic.

To understand how the streaming instability evolves when multiple particle species are considered, instead of performing several runs with different starting metallicities, we removed gas from the system (Carrera et al. 2015). The nominal initial metallicity we implemented is Z_init = 0.01. In addition, we also performed a run with Z_init = 0.005 to ensure that our choice for this value does not influence the outcome. We started removing gas after the first 30 orbits, at which time the particles have saturated to a constant scale height. We decreased the gas density in an exponential manner, such that it is halved every 50 orbits, following $${\rho }_{\text{g}}(t)={\rho }_{\text{g},\text{0}}\exp (-kt).$$(4)

Here, ρ_g,0 is the initial gas density and k = 1∕τ, where τ is the e-folding time and the half-life time is defined as t_rem = τ∕ln(2). We also performed convergence tests, such that the gas density is halved every t_rem = 100 orbits instead (see Appendix A). As the particle density is kept constant overall, this results in continuous metallicity increase with time, where the metallicity is defined as $$Z={\overline{\Sigma }}_{\text{p}}/{\overline{\Sigma }}_{\text{g}},$$(5)

with ${\overline{\Sigma }}_{\text{p}}$ as the particle surface density.

3 Results

3.1 Single-species validation

In order to determine the right criteria for the start of filament formation, we began with single-species simulations. We then compared the results with that of Carrera et al. (2015) and Yang et al. (2017) to ensure that our metallicity limit, Z_crit, holds up against the expectations.

Our criteria for the start of filament formation is ρ_p,max∕ρ_g = 100, where ρ_p,max is the maximum solid density. This criteria is an order of magnitude approximation for the Roche density at 2.5 au in a MMSN disk (Yang et al. 2017). On the top panel of Fig. 1, we show the evolution of the maximum solid density with respect to the gas density of all single species runs, where τ_s = 10⁻³, 10⁻², 10⁻¹, 10⁰, and 10¹. The dashed horizontal line corresponds to ρ_p,max∕ρ_g = 100. As seen on the first panel of Fig. 1, all five particle species reach our prescribed limit. The curve corresponding to τ_s = 10⁻³ has the slowest maximum solid density increase and little small-scale variations, which suggests that the solids sediment to a dense midplane layer without the emergence of significant clumping.

We also studied the importance of the numerical sampling effects by performing runs with the same parameters while increasing the particle number from N_p = 10 to N_p = 20 and to N_p = 100 per cell. The result of this is shown for τ_s,max = 10⁻¹ on the middle and for τ_s,max = 10⁻² on the bottom panel of Fig. 1. Here, we see that both the curves corresponding to N_p = 10 vs. N_p = 20 and N_p = 10 vs. N_p = 100 per cell remain close to unity. The largest variations from unity coincide with the times when the maximum solid density curves shown on the top panel cross the ρ_p,max∕ρ_g = 100 limit. Afterwards, the level of deviation settles down to close to unity again. This implies that numerical parameters such as particle number do not have a significant effect on the evolution of the single-species streaming instability.

In Appendix B, we present the space-time diagrams of the single species runs, where τ_s = 10⁻³, 10⁻², 10⁻¹, 10⁰, and 10¹ and N_p = 10 per cell. As shown in Fig. B.1 by the orange horizontal lines, the threshold metallicity for clumping, Z_crit, is reached in all five cases. The figure also shows that indeed, the filaments above the Z_crit limit are clearly visible and permanent, except in the case of the first panel, where τ_s = 10⁻³. Here, the generated turbulence is not strong enough and the particles settle into a thin, dense layer thus increasing the solid density without the formation oflarge filaments. As Fig. B.2 shows, here we see evidence of some vertically-elongated structures in the vertical component of the gas velocity, which could be the result of the vertical shear instability (Lin 2021). However, overall we see that our Z_crit limit catches the particle filaments well, so we apply this criteria moving forward. Also, we prescribed this criteria not necessarily to give absolute values for solid clumping, but rather to allow us to make comparisons between the importance that various parameters play in the streaming instability. As high resolution runs are not computationally feasible for the large parameter study presented in this paper, we focused on determining the relative differences in Z_crit and thus the parameter space in question.

In Fig. 2, we summarize the Z_crit values for all five particle species. We also show the critical metallicity values presented in Carrera et al. (2015) and Yang et al. (2017) with red and see that our values are somewhat higher. This is especially true in the case of τ_s = 10⁻³ and 10⁻². One reason for the discrepancy is likely that we set the initial metallicity to Z_init = 0.01, whereas in both Carrera et al. (2015) and Yang et al. (2017) Z_init = 0.005. In addition, we also have a different way of defining the beginning of filament formation. Carrera et al. (2015) define particle clumps as temporally stable overdensities in the solid surface density, while Yang et al. (2017) perform high resolution studies with a fixed metallicity to find the critical metallicity by tracking the maximum local solid density. In Fig. 2 we also show the effect of doubling the gas removal time, t_rem. In the case of both τ_s = 10⁻¹ and τ_s = 10⁻² the change in Z_crit is about a factor of one half and in the direction toward the results of Carrera et al. (2015) and Yang et al. (2017).

Fig. 1 Top panel: maximum solid density evolution as a function of time for single species runs with τs = 10−3, 10−2, 10−1, 100 and 101. The metallicity is shown on the top axis. The horizontal black dashed line corresponds to ρp,max ∕ρg = 100, which represents our adopted metallicity threshold for filament formation. Middle and bottom panels: Effect of changes in Np for τs = 10−1 and 10−2, respectively.Here, the solid lines show the relative maximum density for 10 vs. 20 particles per cell while the dashed lines show the relative maximum density for 10 vs. 100 particles per cell. In both cases, the level of deviation from unity is low, implying that the particle number does not have a significant effect on the maximum density of the single-species streaming instability. The colors correspond to the legend in the top panel.

Fig. 2 Critical metallicity threshold for filament formation, Zcrit, with respect to Stokes number for gas removal times of 50 and 100 orbits. With red we show the critical metallicity values measured by Carrera et al. (2015) and Yang et al. (2017) for comparison. These values are somewhat lower due to differences in initial conditions, resolution and in the definitions of the start of filament formation.

3.2 Multispecies results

In this section, we studied how the streaming instability behaves once multiple particle species are considered. We investigated the significance of the number of simulation particles, N_p, maximum Stokes number, τ_s,max, solid size distribution exponent, q, initial metallicity, Z_init, and initial particle scale height, H_p,init, on the resulting metallicity threshold, Z_crit.

We summarized the run parameters in Table 1. The first two columns show the minimum and maximum values of the particle size range. Between these values, the particles are distributed in a continuous way, such that each particle in a given run has a unique size chosen randomly to sample the assumed size distribution. The third column shows q, the exponent of the particle size distribution described as dN∕da ∝ a^−q. Here, N is the particle number density and a the particle size.The value of q can take up different values, based on the physical processes considered. In this paper, we considered q = 3.5, which applies in the case of self-similar fragmentation cascade (Dohnanyi 1969; Williams & Wetherill 1994). In addition, in order to stay consistent with various possible outcomes of solid-gas dynamics, we also implemented solid size distribution exponents of q = 2.5 and q = 3 (Birnstiel et al. 2011). The fourth column in Table 1 shows the number of simulation particles, N_p, per cell. The fifth column shows the initial metallicity, Z_init and the final column the initial particle scale height, H_p,init, of the system.

3.2.1 Particle settling

In order to understand how the multispecies streaming instability behaves, we began by tracking the evolution of the midplane dust-to-gas ratio and the average particle scale height as a function of gas scale height. We present the evolution of ϵ = ρ_p∕ρ_g and H_p ∕H_g in the top and bottom panels of Fig. 3, respectively. The top panel shows that changes in either q or τ_s,max do not affect the evolution of the dust-to-gas ratio significantly, even though all these systems saturate to different scale heights. On Fig. 4 we show the average particle number per cell in the midplane for the system corresponding to the pink curves in Fig. 3. We see that as the particles sediment toward the midplane, the particle number increases here significantly compared to the initial value of N_p = 70 per cell. This shows that the particle number per cell in the mid-plane where the streaming instability evolves is much higher than the mean particle number per cell over the box. Other than our nominal initial particle scale height of H_p,init = 0.015H_g, we also show here the model, where we set H_p,init = 1H_g. The top andbottom panel of Fig. 3 show that both the average dust-to-gas ratio in the midplane and the particle scale height evolve similarly, given both initial particle scale height values. Thus, we find that both an initially stratified and unstratified particle profile produce the same outcome.

There is, however, some deviation in how both ϵ and H_p∕H_g evolve as a function of the maximum Stokes number and the slope of the particle size distribution. In terms of the particle size range, the systems with τ_s = 10⁻³−10⁻² show the least amount of short-term variations in the development of their dust-to-gas ratios. The particle scale height in both the q = 3 and q = 3.5 systems implies that the particles slowly sink toward the midplane and are not stirred significantly.

Given our nominal particle size range of τ_s = 10⁻³−10⁻¹, there is no significant deviation in how the dust-to-gas ratio evolves based on whether q = 2.5, 3 or 3.5. However, the particle scale height corresponding to the system where q = 3.5 shows closer settling to the midplane, while the system with q = 2.5 tends to saturate farther away from the midplane. The smaller the value of the size distribution exponent, the more mass there is in the larger solids. This implies that the system with q = 2.5 results in more turbulence than the system with q = 3 and even more than the one with q = 3.5, meaning that the particles are more efficiently diffused to larger scale heights.

The evolution of ϵ does not change significantly, whether the maximum particle size is τ_s,max = 10⁰ or τ_s,max = 10¹ for either q = 3 or q = 3.5. In all four systems, there is quite a lot of short-term variations in ϵ during the first approximately 200 orbits, which then levels out as the particles settle closer to the midplane. Overall, the level of short-term variations in the dust-to-gas ratio seems to increase with the maximum Stokes number.

Fig. 3 Average dust-to-gas ratio in the midplane (top panel) and particle scale height as a function of gas scale height (bottom panel) with respect to time for all Stokes number ranges and size distribution exponents with Np = 70 per cell. The initial particle scale height is Hp,init = 0.015Hg, except in thecase of the gray curves, where Hp,init = 1Hg. The metallicity is shown on the top axis. The top panel shows that ϵ is not significantly affected by the slope of the size distribution, but the short-time variations increase with τs,max. The bottom panel shows that the particle scale heights saturate closer to the midplane for small size ranges, but there is no clear trend for how Hp ∕Hg changes as a function of q and τs,max.

Fig. 4 Average particle number per cell in the midplane as a function of time. Initially, the particle number is Np = 70 per cell, the particle size range is τs = 10−3–10−1 and the size distribution exponent is q = 3.5. As the particles sediment toward the midplane, the particle number and thus the dust-to-gas ratio increases here significantly. The particle number is around 300 per cell, even in the early stages of the simulation. This is much higher than the mean particle number per cell, since the scale height of the dust layer is much smaller than the vertical extent of the box.

Fig. 5 Space-time plots showing the solid surface density for runs with q = 3.5, τs = 10−3−10−1 and Np = 7, 70, 140, 280 vs. 700 per cell. The last panel, indicated with a dagger, corresponds to the run where Np = 700 per cell and the particles are initialized with a different randomization as in our nominal 700 per cell model, shown on the left-hand side of the last row. The orange line marks the limit of threshold metallicity, Zcrit, where ρp,max∕ρg = 100. The y-axis on the right-hand side shows the metallicity, Z. In all four systems, the critical metallicity is reached after about 100–120 orbits, signifying that increasing the particle number does not make the multispecies streaming instability less efficient in forming filaments.

Fig. 6 Maximum solid density evolution as a function of time comparing the effect of particle number for runs with τs = 10−3−10−1 and q = 3, q = 3.5. The dagger symbol notes the model where we initialized the particle position with a different randomization than that of the system corresponding to the orange curve, with the same particle number. The metallicity is shown on the top axis. The dashed line corresponds to ρp,max ∕ρg = 100, which represents the limit of solid clumping. All curves show similar trends, independent of the particle number and the slope of the solid size distribution.

3.2.2 Influence of particle number

Here, we looked at how the particle number affects the efficiency of filament formation. In the space-time plots of Fig. 5 we compared how N_p = 7, 70, 140, 280 versus N_p = 700 per cell changes the solid surface density evolution of systems with a size distribution exponent of q = 3.5 and particle size range of τ_s = 10⁻³–10⁻¹. We see that the critical metallicity is reached after approximately 100–120 orbits in all four systems, independent of the particle number. Figure 6 shows that the maximum particle density of these systems follows very similar trends as well. In both the case of the N_p = 140 (pink curve) and the nominal N_p = 700 per cell (orange curve) runs, there is a small dip in ρ_p,max∕ρ_g at around 100 orbits, but this is followed by increase again. We note, that the saturation of the maximum particle density seems to be somewhat delayed as the particle number is increased to N_p = 140 per cell (pink curve) and N_p = 700 per cell (orange curve). In Fig. 6 we also show the maximum solid density evolution of the systems where q = 3. Here, we compare N_p = 7 with N_p = 70 per cell and find that there is no significant difference between the corresponding curves.

Overall, Fig. 6 also shows that given the fixed maximum Stokes number of τ_s,max = 10⁻¹, neither the particle number, nor the size distribution exponent makes a significant difference on when the limit for solid clumping is reached. We also note that the slope of the curves presented here is set not only by the efficiency of filament formation but also by the continuous removal of gas, through which we increased the metallicity.

To show the role of the particle number in the efficiency of filament formation, we present Fig. 7. The first panel shows the critical metallicity given q = 2.5 for our nominal particle number of 70 per cell. The second panel corresponds to the case of q = 3, where we found that a factor of ten increase from 7 to 70 particles per cell results in only about a 0.005 increase in Z_crit. In the third panel we show the systems with q = 3.5, where we see that increasing the particle number does not automatically mean an increase in Z_crit. The threshold metallicity for filament formation is increased for 7 compared to 70 particles per cell. As we increased the particle number to 140 per cell, Z_crit is increased slightly but once N_p = 280 per cell, the critical metallicity is decreased again. Finally, when the particle number is increased to 700 per cell, the critical metallicity remains in the same range as in the case of fewer particles per cell, regardless of the method we chose to initialize the particle positions. Most importantly, Z_crit changed by only about a factor of unity as N_p changes by a factor of one hundred. All this implies that numerical parameters, such as particle number, do not have a significant effect on the time of filament formation, and thus the evolution of the streaming instability. We note here, that these metallicity values are calculated given our criteria that solid clumping begins once ρ_p,max∕ρ_g = 100.

We summarized the importance of the particle number in Fig. 8. We clearly see that given q = 3 and 3.5, N_p does not have a significant effect on the critical metallicity of particle clumping and that an increase in particle number does not directly translate to higher Z_crit.

3.2.3 Influence of particle size range

We also examined the effect of the size distribution on the clumping efficiency for different particle size ranges. Figure 9 compares the evolution of runs with τ_s,max = 10⁻², 10⁻¹, 10¹, 10⁰ and q = 2.5, 3, 3.5. As before, the black dashed line shows the critical metallicity of filament formation, namely where ρ_p,max∕ρ_g = 100.

Given a maximum Stokes number of τ_s,max = 10⁻² (first panel),both the system with q = 3 and q = 3.5 evolve at about the same rate for about the first 170 orbits. Then, the q = 3 system overtakes the density growth and the maximum particle density crosses the dashed line at approximately 175 orbits. In the case of the system with less mass in larger solids (q = 3.5), the filament formation begins about 25 orbits later.

Given our nominal Stokes number range of τ_s = 10⁻³–10⁻¹ however, the maximum solid density evolution follows the same trend independent of the size exponent, as shown in the second panel of Fig. 9. The critical metallicity is reached much earlier compared to the case of τ_s,max = 10⁻², after about 100 orbits.

Moving on to the larger particle size range of τ_s = 10⁻³–10⁰, the metallicity threshold for filament formation is reached after approximately 120 orbits given q = 3 and about 50 orbits later given q = 3.5, as shown in the third panel of Fig. 9. Based on the comparison between all size distribution exponents, and maximum Stokes numbers, it is the run with q = 3 and τ_s,max = 10⁰ that is the most successful in forming permanent filaments early on.

Finally, given the largest size range of τ_s = 10⁻³–10¹, filament formation starts almost simultaneously independent of q. The two curves on the bottom panel of Fig. 9 show that the maximum solid densities evolve differently once Z_crit is reached.

We compared the critical metallicities reached in our systems as a function of particle number in Fig. 7. We see on the second and third panels that given the smallest maximum Stokes number of 10⁻², filaments form once the metallicity has increased significantly. When comparing q = 3 and q = 3.5 for the same particle range of τ_s = 10⁻³–10⁻², we see that the latter case corresponds to higher Z_crit values. This is not surprising, as the system with an exponent of q = 3 contains more mass in larger species, which are more successful in driving the streaming instability. As found by Bai & Stone (2010a), it is particles of τ_s ≈ 10⁻²–10⁻¹ that actively participate in the streaming instability and thus aid the formation of filaments. In the case of our nominal size range of τ_s = 10⁻³–10⁻¹, similar Z_crit values are reached independent of whether q = 2.5, 3 or 3.5. We also find that changing the initial particle scale height to H_p,init = 1, such that the particle distribution is unstratified within our simulation domain, results in a similar critical metallicity value as in the case of the stratified system with the same parameters. We vary the initial metallicity as well and see that Z_init = 0.005 corresponds to a Z_crit that agrees well with that of the same system with Z_init = 0.01. Once the maximum Stokes number is increased to 10⁰, there is quite a large change in the resulting critical metallicity when comparing the systems with q = 3 and q = 3.5, as shown in the second and third panels.However, when τ_s,max is increased to 10¹ the Z_crit values are nearly identical for both size exponents that we examined. Overall, there is no clear trend for how the combinationof the particle size distribution exponent and the Stokes number range affect the outcome of the evolution.

Figure 10 serves as comparison between the effect of τ_s,max and q (given N_p = 70 per cell). We also show here how our single-species results presented in Fig. 2, compared to these values. It is clear that in the case of τ_s,max = 10⁻² the critical metallicity threshold approximately doubles as we go from the single-species case to q = 3.5. This transition can be explained by decreasing the amount of large particles in the system and having less active species available. In the case of τ_s,max = 10⁻¹, changing the size distribution has little effect on Z_crit. This could be due to the fact that the size dependency for the streaming instability is generally weak in this parameter space. For τ_s,max > 10⁻¹, changing q does not contribute to significant changes in Z_crit, as in these systems moving mass to smaller species is still keeping a significant fraction of the mass in active particles of τ_s ≈ 10⁻²–10⁻¹ (Bai & Stone 2010b). So overall, there is no clear trend for how changing the Stokes number range changes the efficiency of filament formation.

Fig. 7 Critical metallicity for filament formation, Zcrit, as a function of Stokes number. The initial metallicity and particle scale height are Zinit = 0.01 and Hp,init = 0.015, unless otherwise specified. The colorbar corresponds to the number of particles used in the given simulation and the color of the box edges to the exponent of the size distribution, as shown in the legend. Some of the bar are shifted horizontally for better visibility. The legend on the left-most panel marks the runs, where the particle positions were initialized using different randomization. The horizontal lines cutting through the boxes represent the actual Zcrit values.

Fig. 8 Critical metallicity for filament formation, Zcrit, as a functionof particle number. Here, τs,max is fixed as 10−1, while q = 3 and q = 3.5. There is no significant change in Zcrit with change in particle number for either size distribution exponent.

Fig. 9 Maximum solid density evolution as a function of time comparing the effect of the size distribution exponents and the particle size ranges given Np = 70 per cell. The initial particle scale height is Hp,init = 0.015Hg, except in thecase of the yellow curve, where Hp,init = 1Hg. The metallicity is shown on the top axis. The dashed line corresponds to ρp,max ∕ρg = 100, which represents our assumed threshold for solid clumping. Independent of q, the maximum solid density evolves similarly in the first, second and last panels until the critical metallicity limit is reached. Given that τs,max = 100 however, there is a deviation of about 50 orbits between the time of filament formation, depending on whether q = 3 or q = 3.5. There is also no clear trend on how ρp,max∕ρg changes as both a function of q and τs,max.

Fig. 10 Critical metallicity for filament formation, Zcrit, as a functionof particle number and maximum Stokes number with Hp,init = 0.015 and Zinit = 0.01. We show the effect of changing q, while keeping the particle number fixed as Np = 70 per cell. We also show the single-species values from Fig. 2 for comparison and see that the systems with more mass in large particles (q = 3) correspond to lower Zcrit values than those with q = 3.5. This is however not the case when τs,max = 10−1 or 101. Given these maximum Stokes numbers (or Stokes numbers given the single-species case), the critical metallicity values are similar, independent of the size distribution exponent.

4 Conclusion

In this paper, we studied how the multispecies streaming instability behaves given change in the following parameters: particle number, particle size range, and number density distribution exponent. In order to compare how each of these parameters influence the outcome, we compared the critical metallicity where filament formation begins. We point out that these Z_crit values are not necessarily absolute, given the relatively low resolution of our simulations, but rather serve as a way to compare the effects of the previously mentioned parameters on particle clumping. In Figs. 7 and 8, we summarized these results and showed that varying the particle number does not lead to a significant change in the efficiency of filament formation, given our chosen criteria. We note however, that as Fig. 5 shows, whether the formed filaments merge or not, can vary as a function of particle number. On the other hand, if self-gravity was included, these filaments would likely form planetesimals around the time the critical metallicity is reached. We also found that there is no correlation between how changing the size distribution exponent affects the outcome.

To understand the dynamics of small particles well, 3D simulations are needed (Bai & Stone 2010a) and in the future the systems discussed in this paper should be tested in 3D as well. In this work, we were limited by fixed resolution box size and somewhat short simulation times. However, keeping such numerical parameters fixed makes comparisons between others more straightforward.

Our results imply that given the realistic case of multiple particle species, the particle size distribution, initial metallicity, and initial particle scale height do not have a significant effect on the outcome of the efficiency and timing of filament formation, given our chosen criteria. Our systems reached similar dust-to-gas ratios (see Fig. 3) – which is set naturally as a consequence of the particle sedimentation – and thus reached the critical midplane density for clumping around the same time. Another result of particle sedimentation is universally high-density conditions in the midplane. These are conditions under which linear multispecies streaming instability models also find short growth timescales (Krapp et al. 2019; Paardekooper et al. 2020, 2021; McNally et al. 2021). As a consequence, the multispecies streaming instability could indeed be successful in forming solid clumps.

Acknowledgements

We thank the anonymous referee for their comments that helped improve the manuscript. N.S. was funded by the “Bottlenecks for particle growth in turbulent aerosols” grant from the Knut and Alice Wallenberg Foundation (2014.0048). N.S. is thankful to Daniel Carrera for useful discussions. A.J. thanks the Swedish Research Council (grant 2018-0486), the Knut and Alice Wallenberg Foundation (grant 2017.0287) and the European Research Council (ERC Consolidator Grant 724687-PLANETESYS) for research support. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at LUNARC in Lund University.

Appendix A Gas removal time

Fig. A.1 Maximum solid density evolution as a function of time comparing the effect of the time of gas removal, such that the gas density is halved every 50 or 100 orbits, respectively. Here, q = 3.5, τs = 10−3–10−1 and both systems have Np = 7 particles per cell. The dashed line corresponds to ρp,max∕ρg = 100, which represents the critical maximum solid density needed for filament formation. There is a difference of about 0.02 in the critical metallicity, depending on the gas removal time.

Here, we tested whether the speed at which the gas is removed from the system matters for when the critical metallicity of filament formation is reached. Other than the nominal gas removal timescale of 50 orbits, we also studied systems where the gas density is halved every 100 orbits instead. In Table A.1, we compare the critical metallicity for both cases. We expect that the resulting Z_crit values do not significantly depend on the removal time. We see that ΔZ_crit ≈ 0.02 given a 50 orbit delay between t_rem = 50 orbits versus t_rem = 100 orbits. Overall, as also seen in Fig. 2 our nominal removal timescale of t_rem = 50 orbits serves as a higher limit for the Z_crit of filament formation.

Appendix B Single species solid density

In Fig. B.1, we show the space-time figures of the runs for particle sizes of τ_s = 10⁻³, 10⁻², 10⁻¹, 10⁰, and 10¹. The orange lines mark the critical metallicity where the ρ_p,max∕ρ_g limit is met.All five runs meet this criteria and for a more detailed discussion see Sect. 3.1.

In Fig. B.2, we show a snapshot of the vertical component of the gas velocity as a function of sound speed at t ~ 260 2πΩ⁻¹. Here, the particles are in the process of settling and form a thick particle layer around the midplane. We found the presence of vertically-elongated structures in the vertical component of the gas velocity above the diffused particle layer (Lin 2021).

Fig. B.1 Space-time plot of runs with τs = 10−3, 10−2, 10−1, 100 and 101 with respect to radial coordinate, x and orbital time. The y-axis on the right-hand side shows the metallicity, Z. The colors correspond to the solid surface density and the orange line to ρp,max ∕ρg = 100. This criteria fits well with the start of filament formation in all systems, except in the first panel, where τs = 10−3. Here, the limit is reached without the formation of any visible filaments.

Fig. B.2 Vertical component of the gas velocity as a function of sound speed of the system with τs = 10−3 at t ~ 260 2πΩ−1. Above the sedimenting particle layer, we find the presence of vertically-elongated structures in uz.

References

All Tables

All Figures

Fig. 1 Top panel: maximum solid density evolution as a function of time for single species runs with τs = 10−3, 10−2, 10−1, 100 and 101. The metallicity is shown on the top axis. The horizontal black dashed line corresponds to ρp,max ∕ρg = 100, which represents our adopted metallicity threshold for filament formation. Middle and bottom panels: Effect of changes in Np for τs = 10−1 and 10−2, respectively.Here, the solid lines show the relative maximum density for 10 vs. 20 particles per cell while the dashed lines show the relative maximum density for 10 vs. 100 particles per cell. In both cases, the level of deviation from unity is low, implying that the particle number does not have a significant effect on the maximum density of the single-species streaming instability. The colors correspond to the legend in the top panel.

In the text

Fig. 2 Critical metallicity threshold for filament formation, Zcrit, with respect to Stokes number for gas removal times of 50 and 100 orbits. With red we show the critical metallicity values measured by Carrera et al. (2015) and Yang et al. (2017) for comparison. These values are somewhat lower due to differences in initial conditions, resolution and in the definitions of the start of filament formation.

In the text

Fig. 3 Average dust-to-gas ratio in the midplane (top panel) and particle scale height as a function of gas scale height (bottom panel) with respect to time for all Stokes number ranges and size distribution exponents with Np = 70 per cell. The initial particle scale height is Hp,init = 0.015Hg, except in thecase of the gray curves, where Hp,init = 1Hg. The metallicity is shown on the top axis. The top panel shows that ϵ is not significantly affected by the slope of the size distribution, but the short-time variations increase with τs,max. The bottom panel shows that the particle scale heights saturate closer to the midplane for small size ranges, but there is no clear trend for how Hp ∕Hg changes as a function of q and τs,max.

In the text

Fig. 4 Average particle number per cell in the midplane as a function of time. Initially, the particle number is Np = 70 per cell, the particle size range is τs = 10−3–10−1 and the size distribution exponent is q = 3.5. As the particles sediment toward the midplane, the particle number and thus the dust-to-gas ratio increases here significantly. The particle number is around 300 per cell, even in the early stages of the simulation. This is much higher than the mean particle number per cell, since the scale height of the dust layer is much smaller than the vertical extent of the box.

In the text

Fig. 5 Space-time plots showing the solid surface density for runs with q = 3.5, τs = 10−3−10−1 and Np = 7, 70, 140, 280 vs. 700 per cell. The last panel, indicated with a dagger, corresponds to the run where Np = 700 per cell and the particles are initialized with a different randomization as in our nominal 700 per cell model, shown on the left-hand side of the last row. The orange line marks the limit of threshold metallicity, Zcrit, where ρp,max∕ρg = 100. The y-axis on the right-hand side shows the metallicity, Z. In all four systems, the critical metallicity is reached after about 100–120 orbits, signifying that increasing the particle number does not make the multispecies streaming instability less efficient in forming filaments.

In the text

Fig. 6 Maximum solid density evolution as a function of time comparing the effect of particle number for runs with τs = 10−3−10−1 and q = 3, q = 3.5. The dagger symbol notes the model where we initialized the particle position with a different randomization than that of the system corresponding to the orange curve, with the same particle number. The metallicity is shown on the top axis. The dashed line corresponds to ρp,max ∕ρg = 100, which represents the limit of solid clumping. All curves show similar trends, independent of the particle number and the slope of the solid size distribution.

In the text

Fig. 7 Critical metallicity for filament formation, Zcrit, as a function of Stokes number. The initial metallicity and particle scale height are Zinit = 0.01 and Hp,init = 0.015, unless otherwise specified. The colorbar corresponds to the number of particles used in the given simulation and the color of the box edges to the exponent of the size distribution, as shown in the legend. Some of the bar are shifted horizontally for better visibility. The legend on the left-most panel marks the runs, where the particle positions were initialized using different randomization. The horizontal lines cutting through the boxes represent the actual Zcrit values.

In the text

	Fig. 8 Critical metallicity for filament formation, Zcrit, as a functionof particle number. Here, τs,max is fixed as 10−1, while q = 3 and q = 3.5. There is no significant change in Zcrit with change in particle number for either size distribution exponent.
In the text

Fig. 9 Maximum solid density evolution as a function of time comparing the effect of the size distribution exponents and the particle size ranges given Np = 70 per cell. The initial particle scale height is Hp,init = 0.015Hg, except in thecase of the yellow curve, where Hp,init = 1Hg. The metallicity is shown on the top axis. The dashed line corresponds to ρp,max ∕ρg = 100, which represents our assumed threshold for solid clumping. Independent of q, the maximum solid density evolves similarly in the first, second and last panels until the critical metallicity limit is reached. Given that τs,max = 100 however, there is a deviation of about 50 orbits between the time of filament formation, depending on whether q = 3 or q = 3.5. There is also no clear trend on how ρp,max∕ρg changes as both a function of q and τs,max.

In the text

Fig. 10 Critical metallicity for filament formation, Zcrit, as a functionof particle number and maximum Stokes number with Hp,init = 0.015 and Zinit = 0.01. We show the effect of changing q, while keeping the particle number fixed as Np = 70 per cell. We also show the single-species values from Fig. 2 for comparison and see that the systems with more mass in large particles (q = 3) correspond to lower Zcrit values than those with q = 3.5. This is however not the case when τs,max = 10−1 or 101. Given these maximum Stokes numbers (or Stokes numbers given the single-species case), the critical metallicity values are similar, independent of the size distribution exponent.

In the text

Fig. A.1 Maximum solid density evolution as a function of time comparing the effect of the time of gas removal, such that the gas density is halved every 50 or 100 orbits, respectively. Here, q = 3.5, τs = 10−3–10−1 and both systems have Np = 7 particles per cell. The dashed line corresponds to ρp,max∕ρg = 100, which represents the critical maximum solid density needed for filament formation. There is a difference of about 0.02 in the critical metallicity, depending on the gas removal time.

In the text

Fig. B.1 Space-time plot of runs with τs = 10−3, 10−2, 10−1, 100 and 101 with respect to radial coordinate, x and orbital time. The y-axis on the right-hand side shows the metallicity, Z. The colors correspond to the solid surface density and the orange line to ρp,max ∕ρg = 100. This criteria fits well with the start of filament formation in all systems, except in the first panel, where τs = 10−3. Here, the limit is reached without the formation of any visible filaments.

In the text

	Fig. B.2 Vertical component of the gas velocity as a function of sound speed of the system with τs = 10−3 at t ~ 260 2πΩ−1. Above the sedimenting particle layer, we find the presence of vertically-elongated structures in uz.
In the text