Streaming Instability for Particle-Size Distributions

TL;DR

This study investigates how a distribution of particle sizes affects the growth and efficiency of the streaming instability, a key process in planetesimal formation, revealing that multi-species systems grow more slowly than single-species ones.

Contribution

First analysis of the linear streaming instability considering a particle-size distribution, showing slower growth rates and implications for planet formation models.

Findings

Multi-species streaming instability grows on longer timescales than single-species.

High dust-to-gas ratios can lead to growth times comparable to secular instabilities.

Particle segregation may enhance instability growth.

Abstract

The streaming instability is thought to play a central role in the early stages of planet formation by enabling the efficient bypass of a number of barriers hindering the formation of planetesimals. We present the first study exploring the efficiency of the linear streaming instability when a particle-size distribution is considered. We find that, for a given dust-to-gas mass ratio, the multi-species streaming instability grows on timescales much longer than those expected when only one dust species is involved. In particular, distributions that contain close-to-order-unity dust-to-gas mass ratios lead to unstable modes that can grow on timescales comparable, or larger, with those of secular instabilities. We anticipate that processes leading to particle segregation and/or concentration can create favourable conditions for the instability to grow fast. Our findings may have important implications for a large number of processes in protoplanetary disks that rely on the streaming instability as usually envisioned for a unique dust species. Our results suggest that the growth rates of other resonant-drag-instabilities may also decrease considerably when multiple species are considered.

Full text

protoplanetary disks – instabilities – hydrodynamics

††software: NumPy (Walt et al., 2011), Matplotlib (Hunter, 2007)

1 Introduction

The building blocks of planetary bodies are kilometer-sized planetesimals, which are believed to form and grow from small dust-particles present in the protoplanetary disk. When growing from micro- to meter-sized objects, the particles need to overcome the so-called “growth barriers” (see Testi et al., 2014, for a review). In particular, under typical conditions expected in protoplanetary disks, the radial-drift barrier prevents particles from growing beyond centimeter scales. This corresponds roughly to the particle size for which the timescale involved in particle growth is comparable to the timescale associated with their radial drift (Whipple, 1972; Weidenschilling, 1977; Nakagawa et al., 1986). One mechanism envisioned to overcome this barrier is the streaming instability (Youdin & Goodman, 2005; Youdin & Johansen, 2007; Jacquet et al., 2011; Auffinger & Laibe, 2018).

Even though there have already been a handful of papers reporting numerical simulations exploring the effects of dust-size distributions in the non-linear outcome of the streaming instability (see e.g. Bai & Stone, 2010a; Schaffer et al., 2018), a systematic study addressing its linear regime is yet to be presented. Such a study is necessary for a number of reasons. Conceptually, it provides a more sensible framework to assess the efficiency of the streaming instability in more realistic models of protoplanetary disks. From a computational perspective, it provides valuable benchmarks against which to test numerical codes (e.g. Benítez-Llambay et al., 2019), as well as an idea of the numerical requirements to recover the proper evolution of the instability during its early (linear) phase.

In this Letter, we present the first study of the linear phase of the streaming instability involving a distribution of dust-particles of different sizes. Our systematic exploration of parameter space allows us to provide the growth rate of the most unstable mode as a function of the dust-to-gas mass ratio, particle-size range, and number of dust species considered for describing a particle-size distribution.

2 Multi-species Streaming Instability

2.1 Dust-size Distribution

We consider an underlying (continuous) number-density distribution of particles, which is a power-law of index $q$ in the particle size, $a$ (Dohnanyi, 1969; Mathis et al., 1977). In what follows, we take $q=-3.5$ for definitiveness111We have obtained similar results to those presented here using $q=-4$, which corresponds to a distribution with equal mass per logarithmic bin in Stokes number. Our analysis can be easily generalized to accommodate for more sophisticated distributions. (see e.g., Garaud et al., 2004). We consider the Epstein regime in which the Stokes number is proportional to the particle size, i.e., $T_{\rm s}\propto a$.

A discrete approximation of this dust-size distribution is characterized by the total gas-to-dust mass ratio, $\epsilon$, a range of Stokes numbers properly bound by a minimum and a maximum, $\Delta T_{\rm s}=[T_{\rm s,min},T_{\rm s,max}]$, and the total number of species $N$ associated with a distinct Stokes number, $T_{{\rm s},i}$, with $i=1,\ldots,N$. The total mass of the distribution remains constant when varying $N$, i.e., $\sum_{i=1}^{N}\epsilon_{i}=\epsilon$, where $\epsilon_{i}$ is the dust-to-gas mass ratio associated with a given dust species. We define the Stokes numbers characterizing the distribution evenly in logarithmic scale222The results of this work do not change qualitatively if we distribute the Stokes numbers so that they all contain the same dust-to-gas mass ratio per species..

A robust discrete approximation of the underlying particle-size distribution should in principle lead to a dynamical model that converges as the number of dust species considered increases, i.e., as the continuous limit is approached. Therefore, it is of particular interest to understand the sensitivity of the results obtained with respect to the number of dust species, $N$, used to describe the underlying dust-size distribution.

2.2 Linear Modes in Fourier Space

The equations describing the dynamics of a gas coupled to $N$ dust species via drag forces in the framework of the shearing box, together with the analytical steady-state background solution, have been recently derived in Benítez-Llambay et al. (2019), see Section 3.5. Linear axisymmetric perturbations with respect to the steady-state background lead to $4(N+1)$ (continuity and momentum) equations describing the multiple-species streaming instability, presented in their Appendix E. A given eigenmode of this linear system evolves in space and time according to $\textrm{Re}[\delta\hat{f}(k_{x},k_{z})\,e^{i(k_{x}x+k_{z}z)-\omega t}]$. Here, $\delta\hat{f}(k_{x},k_{z})$ are the $4(N+1)$-dimensional (complex) eigenvectors in Fourier space, spanned by the wavenumbers $(k_{x},k_{z})$, associated with the (complex) eigenvalue $\omega(k_{x},k_{z})$. In the context of the streaming instability, it is customary to work with dimensionless wavenumbers $K=H_{0}^{2}k/R_{0}$, where $H_{0}$ is the disk sale-height at the fiducial radius $R_{0}$, where the shearing box is centered, and to use the Keplerian angular frequency, $\Omega_{0}\equiv\Omega_{\rm K}(R_{0})$, to scale the eigenvalues.

The early evolution of the instability is governed by the unstable modes – i.e., those with $\textrm{Re}(\omega)<0$ – with maximum growth rate333i.e., $\sigma=\textrm{max}\{|\textrm{Re}(\omega)|:\textrm{Re}(\omega)<0\}$. $\sigma$. In order to identify these modes, we consider the space spanned by the set $(K_{x},K_{z})$ when each normalized wavenumber takes values in the range $[10^{-1},10^{3}]$. Our fiducial grid is evenly spaced in logarithmic scale and contains $260$ cells in each direction. Given a dust-size distribution, the dynamical evolution of a specific mode is completely determined by the spectrum of $4(N+1)$ complex eigenvalues $\omega$. This spectrum is found using the function eig of NumPy (Walt et al., 2011), which uses LAPACK routines for complex nonsymmetric matrices (Anderson et al., 1999).

2.3 Fastest Growing Modes – Two Test Cases

We consider two discrete particle-size distributions both with $\epsilon=1$, but spanning two overlapping ranges of Stokes numbers: $\Delta T_{\rm s}^{\rm I}=\left[10^{-4}\,,10^{-1}\right]$ and $\Delta T_{\rm s}^{\rm II}=\left[10^{-4}\,,1\right]$. We compute the growth rate using $N\in\{16,128,512\}$ dust species. These considerations lead to six different eigenvalue problems that are solved to find the fastest growing modes as a function of $(K_{x},K_{z})$. The results corresponding to $\Delta T_{\rm s}^{\rm I}$ and $\Delta T_{\rm s}^{\rm II}$, for each of the adopted $N$-values, are shown in the upper and lower panels of Fig. 1, respectively. For $\Delta T_{\rm s}^{\rm I}$, the upper panels show a maximum growth rate that converges with increasing dust species to $\sigma\simeq 1.6\times 10^{-2}\Omega_{0}$ (see Fig. 2). For this distribution, the set of modes that grow fastest converge to a confined region close to the center of the explored domain in $(K_{x},K_{z})$. In contrast, for $\Delta T_{\rm s}^{\rm II}$ the maximum growth rate decreases monotonically from $\sigma\simeq 6.7\times 10^{-2}\Omega_{0}$ for 16 species to $\sigma\simeq 0.33\times 10^{-2}\Omega_{0}$ for 512 species (see also Fig. 2).

The sensitivity of the results obtained for the fastest growth rate with respect to the number of species $N$ can be better appreciated in the leftmost panel in Fig. 2, which shows the growth rates of the most unstable modes for the two distributions with $\Delta T_{\rm s}^{\rm I}$ (orange line) and $\Delta T_{\rm s}^{\rm II}$ (blue line), when the number of dust species doubles from $N=2$ to $N=2048$. For the case $\Delta T_{\rm s}^{\rm I}$ the maximum growth rate converges when using $64$ dust species. This is not the case for $\Delta T_{\rm s}^{\rm II}$, for which the maximum growth rate decreases below $\sigma<10^{-3}\Omega_{0}$. However, in a region around $K_{x}=10^{-1}$, the growth rate converges to $\sim 5\times 10^{-4}\Omega_{0}$ for $N\geq 64$. This value sets the time scale of the linear instability for $N>2048$ species.

In order to shed some light on the strikingly different behaviour exhibited by $\Delta T^{\rm I}$ and $\Delta T^{\rm II}$, we show in Fig. 3 high resolution maps in $(K_{x},K_{z})$ zooming-in the neighborhood of the fastest growing modes. The left and right panels show the maximum growth rate for 64 and 128 species, respectively. The upper and lower panels correspond to $\Delta T_{\rm s}^{\rm I}$ and $\Delta T_{\rm s}^{\rm II}$ (orange and blue curves in Fig. 2), respectively. These maps reveal a pattern with fringes whose number increases linearly with $N$, as they split unstable regions whose growth rates decay also linearly with $N$. When these fringes merge, i.e., their separation is smaller than their width, the growth rates converge. It is worth stressing that the imaginary parts corresponding to the most unstable eigenvalues (shown in the middle panel of Fig. 2) do converge as $N$ increases in both cases.

The decay observed in the growth rate of the most unstable modes when $N$ increases in the case $\Delta T_{\rm s}^{\rm II}$ is in stark contrast with the behavior observed for $\Delta T_{\rm s}^{\rm I}$. This warrants a systematic exploration of parameter space including also the dust-to-gas mass ratio. Before embarking on this, and given the complexity of the equations involved in the linear mode analysis describing the multi-species streaming instability, we provide an independent check of our solutions below.

2.4 Verification of the Linear Mode Analysis

We test the time evolution of the most unstable modes using the publicly available multifluid code FARGO3D (Benítez-Llambay & Masset, 2016; Benítez-Llambay et al., 2019), following the procedure described in section 3.5.4 of Benítez-Llambay et al. (2019).

We consider four representative cases from Fig. 2 (red and green filled circles) and use the corresponding eigenvectors to initialize four numerical simulations. The rightmost panels of Fig. 2 show the time evolution for one of the components (the dust density for species 1, $\delta\rho_{1}$) for each of these four modes. The solutions obtained with FARGO3D are shown with red and green unfilled circles for 16 and 128 dust-species, respectively. The black solid lines are the solutions obtained from our linear mode analysis described in Section 2.2. The excellent agreement between the time evolution of the selected eigenmodes provides additional support to our linear calculations. This critically reduces the possibility of potential issues in several steps of our analysis including the derivation of the background equilibrium, the linearization of the perturbed system, and the method used to find the eigenvalues and eigenvectors.

2.5 Systematic Parameter Space Exploration

We seek the growth rate of the most unstable mode given a particle-size distribution characterized by the dust-to-gas mass ratio, $\epsilon$, a range of Stokes numbers, $\Delta T_{\rm s}=[T_{\rm s,min},T_{\rm s,max}]$, and the total number of species $N$. We consider four different mass ratios, $\epsilon=\{0.01,0.1,0.5,1\}$ and two sets of intervals in Stokes numbers for which either the minimum is fixed and the maximum varies, i.e., $\Delta T_{{\rm s,min}}=[10^{-4},T_{{\rm s,max}}]$ with $T_{{\rm s,max}}\in\{10^{-3},10^{-2},10^{-1},1\}$, or the maximum is fixed and the minimum varies, i.e., $\Delta T_{\rm s,max}=[T_{\rm s,min},1]$ with $T_{\rm s,min}=\{10^{-1},10^{-2},10^{-3},10^{-4}\}$. For each of these intervals in Stokes numbers, we consider an increasing number of dust-species by doubling $N$ from $2$ to $2048$ while keeping the dust-to-gas mass ratio characterizing the distribution constant. This procedure leads to $4\times(4\times 2-1)\times 11=308$ independent discrete distributions.

2.6 Results

The maximum growth rates for each of the distributions defined above are shown in Fig. 4. Each panel corresponds to a different dust-to-gas mass ratio $\epsilon$. The rows and columns correspond to a given $N$ and $\Delta T_{\rm s}$, respectively. Each cell is color-coded according to the logarithm of the maximum growth rate obtained in $(K_{x},K_{z})$-space, following the method described in Section 2.

The most relevant outcomes are (i) the growth rate of the most unstable modes corresponding to the majority of the distributions with low dust-to-gas mass ratios $\epsilon\lesssim 0.1$ have not converged and decreases below $10^{-3}\Omega_{0}$, independently of $T_{\rm s,max}$. In particular, when $T_{\rm s,min}=10^{-4}$ is fixed, (Fig. 4, upper panels), the upper bound for the grow rate decreases from $10^{-3}\Omega_{0}$ to $10^{-5}\Omega_{0}$ as $T_{\rm s,max}$ decreases. (ii) The range of Stokes numbers for which convergence of the growth rate with the number of species is reached increases with $\epsilon$ when $T_{\rm s,min}=10^{-4}$. (iii) When fixing $T_{\rm s,max}=1$ (Fig. 4, lower panels) convergence of the growth rate, with $N=2048$, is achieved for none of the cases considered but one. The only exception is the case that corresponds to $\epsilon=1$ and $T_{\rm s,min}=10^{-1}$, for which the most unstable mode has a growth rate $\sigma\simeq 6\times 10^{-3}\Omega_{0}$. In all other cases, the growth rate decreases below $10^{-3}\Omega_{0}$, independently of the dust-to-gas mass ratio, $\epsilon$.

3 Outcome and Implications

We have provided the first systematic study of the linear growth of the multi-species streaming instability. We found two different types of behaviors. On the one hand, there are distributions for which convergence of the growth rates is reached by considering between a handful and a couple of hundred dust species. In the majority of cases we considered, however, we were only able to find upper limits to the growth rates, which are, in many cases, well below the values obtained when only one dust species is involved. This result is better appreciated in Fig. 5, where we show the maximum growth rates for the classical (gas and one dust-species) streaming instability (leftmost panel) together with those obtained for the distributions with $2048$ species studied in Section 2.5 (center and rightmost panels).

Below, we briefly discuss some of the most relevant aspects and consequences of our findings.

Growth Rate Decay and Connection with Resonant Drag Instabilities

For a large fraction of the particle distributions we studied, the maximum growth rate $\sigma$ decreases as the number of species $N$ increases (and the mass per species decreases). A concrete example is illustrated in Fig. 2, where the growth for the distribution with $\epsilon=1$ and $\Delta T_{\rm s}^{\rm II}=\left[10^{-4}\,,1\right]$ decreases as $1/N$ for sufficiently large $N$. It is natural to compare this to the growth obtained when considering a single dust species with a dust-to-gas mass ratio $\epsilon_{1}\sim\epsilon/N$ corresponding to the species with Stokes number that in isolation leads to fastest growth.

We found that for $\epsilon<10^{-2}$ the maximum growth rate of the distribution is in good agreement with the value obtained when considering in isolation the species with largest Stokes number, $T_{\rm s}=1$ in this case (see Fig. 6). However, as the total mass of the distribution considered is increased, a significant difference between these growth rates exists. This is due to two effects that are difficult to disentangle; i) the background drift-velocity for each individual dust-species is modified because it is sensitive to the total mass of the particle-size distribution and not just the mass per bin (see Benítez-Llambay et al. 2019) and ii) as the ensuing gas perturbation increases, the coupling between species increases and interference among them may not be negligible.

In the limit of very small dust-to-gas mass ratio, $\epsilon\ll 1$, we have checked that the $N$ dust-species streaming instability can be well described as the superposition of $N$ different two-fluid instabilities occurring in a seemingly independent way. This suggests that, in this regime, the resonant drag instability (RDI) framework (Squire & Hopkins, 2018) can help provide insight into the behavior of the multi-species streaming instability. At sufficiently low dust-to-gas mass ratios the maximum growth rate is expected to decay as $\sqrt{\epsilon_{1}}\sim\sqrt{\epsilon/N}$. This is indeed obtained for distributions with dust-to-gas mass ratios smaller than $\epsilon=0.01$ (see Fig. 6).

Single-dust Species Models of Streaming Instability

We report the emergence of an unstable region for $K_{x}\gg 1$ and $K_{z}\geq 0$, which has not been observed before (see Youdin & Goodman, 2005; Youdin & Johansen, 2007). A good example is provided by the two distributions studied in Fig. 1. While a comprehensive study of this new unstable region is beyond the scope of this work, we found that the necessary condition for it to appear is the presence of at least two dust-species with opposite background drift-velocities. This observation suggests that it may not be possible to capture the full dynamics of multi-species streaming instabilities using single dust-species models (e.g., Laibe & Price, 2014; Lin & Youdin, 2017).

Simulations of Multi-Species Streaming Instability

The successful recovery of known solutions is a key benchmark for any numerical code. Previous numerical studies of the non-linear evolution of the streaming instability with one dust species have been reported to recover, for example, its linear phase (see e.g. Johansen & Youdin, 2007; Balsara et al., 2009; Kowalik et al., 2013; Chen & Lin, 2018; Riols & Lesur, 2018; Benítez-Llambay et al., 2019). Even though the linear results have been derived using an Eulerian formalism, it has been shown that numerical codes evolving Lagrangian particles agree very well during the early phases of the streaming instability (Youdin & Johansen, 2007). This suggests that our findings will also hold in the Lagrangian framework. It is possible that a significant decay of the growth rate has not yet been observed in multiple dust-species simulations because of the relatively low number of species that have been used so far (see e.g. Bai & Stone, 2010b; Schaffer et al., 2018).

Critical Dust-to-gas Mass Ratio

Our study suggests that particle-size distributions with $\epsilon\geq 0.5$ are required to allow the multi-species streaming-instability to grow on timescales shorter than $10^{3}\Omega_{0}^{-1}$ (see Fig. 5). This, however, depends on the range of Stokes numbers defining the distribution. For example, if $T_{\rm s,max}=1$ the growth rates are smaller than $10^{-3}\Omega_{0}$, even for large dust-to-gas mass ratios (i.e., $\epsilon=1$). When $\epsilon=0.5$ and $T_{\rm s,min}=10^{-4}$ the maximum growth rate converges to $\sigma\simeq 10^{-3}\Omega_{0}$ for $T_{\rm s,max}=10^{-1}$, and the instability can grow faster for particle-size distributions with $T_{\rm s,max}\leq 10^{-2}$. It can also grow faster if the dust-to-gas mass ratio increases to $\epsilon\simeq 1$ for those ranges of Stokes numbers. We additionally found that, if the total mass of the distribution decreases below $\epsilon=0.5$, the instability develops on timescales of the order of $10^{5}\Omega_{0}^{-1}$, or even longer, depending on the range of Stokes numbers spanned by the particle distribution.

Planetesimal Formation

We anticipate that the multi-species streaming instability could still be an efficient mechanism to enable planetesimal formation if dust-particles are filtered/segregated according to their size and accumulated somewhere in the disk. This will naturally produce regions with large concentrations of dust with distributions characterized by specific particle-sizes. For instance, vertical sedimentation affected by the presence of winds (e.g Riols & Lesur, 2018) or turbulence sustained by the vertical shear instability (e.g., Lin, 2019) can favor the accumulation of larger grains at the mid-plane of protoplanetary disks. Some other potential mechanisms for such filtering/segregation are vortices (e.g., Barge & Sommeria, 1995; Raettig et al., 2015; Ragusa et al., 2017), zonal flows (e.g., Johansen et al., 2009; Dittrich et al., 2013; Béthune et al., 2016; Krapp et al., 2018), planet-induced pressure bumps (e.g. Zhu et al., 2012; Pinilla et al., 2012; Weber et al., 2018), and planetary torques (Benítez-Llambay & Pessah, 2018; Chen & Lin, 2018). Another candidate in this regard is the growth (via dust coagulation) to sizes that are limited by particle-drift, fragmentation or bouncing (Testi et al., 2014). Such mechanisms naturally appear to favour accumulating significant amount of mass into (near-)monodisperse populations of particles.

We conclude that a properly resolved particle-size distribution can significantly affect the linear phase of the streaming instability. Depending on the dust-size distribution and dust-to-gas mass ratio, the multi-species instability may only grow on timescales much larger than those expected from the classical (gas and one dust-species) case when approaching the continuum limit. In particular, distributions that contain moderate to high dust-to-gas mass ratios (i.e., $\epsilon\lesssim 1$) will only grow on timescales comparable with those of secular instabilities. Taken at face value, our results imply that the scope of the streaming instability may be narrowed down profoundly. Nevertheless, processes leading to particle segregation and/or concentration may create favourable conditions for the instability to develop. Because the growth rate of the multi-species streaming instability depends sensitively on the number of dust species used to represent a distribution, our results may also have important implications for the wider class of resonant-drag-instabilities discussed in Squire & Hopkins (2018).

We thank the referee, Jeremy Goodman, for encouraging us to connect our findings, and their implications, with the more general class of resonant-drag-instabilities. We thank Troels Haugbølle for useful discussions and helpful suggestions. The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 638596) (LK, OG). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 748544 (PBLL). MEP gratefully acknowledges support from the Independent Research Fund Denmark (DFF) via grant no. DFF 8021-00400B.