Requirements for Gravitational Collapse in Planetesimal Formation—The Impact of Scales Set by Kelvin–Helmholtz and Nonlinear Streaming Instability

If the cloud is embedded in turbulent gas, and its particles are coupled to the gas, i.e., St < 1, particles will be subject to turbulent diffusion. We quantify the strength of this diffusion with the dimensionless parameter δ, such that the diffusion coefficient is given by

$$\begin{eqnarray}&&D=\delta {c}_{{\rm{s}}}H.\end{eqnarray} \tag{ 4 }$$

Here, we introduced the sound speed c_s and gas pressure scale height

$$\begin{eqnarray}&&H=\displaystyle \frac{{c}_{{\rm{s}}}}{{\rm{\Omega }}}.\end{eqnarray} \tag{ 5 }$$

Diffusion is assumed to be isotropic. We discuss and test this assumption in Appendix A. Following Fick's second law for diffusion, particles are diffused over the length scale r with the diffusion timescale (Youdin & Lithwick 2007)

$$\begin{eqnarray}&&{t}_{\mathrm{diff}}=\displaystyle \frac{{r}^{2}}{\delta {c}_{{\rm{s}}}H}.\end{eqnarray} \tag{ 6 }$$

Thus, cloud collapse will be prevented by turbulent diffusion if t_diff < t_con or

$$\begin{eqnarray}&&\displaystyle \frac{{r}^{2}}{\delta {c}_{{\rm{s}}}H}\lt \displaystyle \frac{{\rm{\Omega }}}{4\pi G{\rho }_{{\rm{c}}}\mathrm{St}}.\end{eqnarray} \tag{ 7 }$$

Assuming that the cloud extends to the vertical particle scale height H_p, i.e., if H_p ≈ r, vertical integration of ρ_c yields the cloud particle column density Σ_p,c, i.e.,

$$\begin{eqnarray}&&{\rho }_{{\rm{c}}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}{2r}.\end{eqnarray} \tag{ 8 }$$

Thus, we can transform Equation (7) to a diffusion-limited stability criterion

$$\begin{eqnarray}&&r\lt \displaystyle \frac{\delta }{\mathrm{St}}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{2\pi G{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}=: {L}_{\mathrm{diff}}.\end{eqnarray} \tag{ 9 }$$

If the radius r of a cloud with column density (mass per area) ${{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}$ falls below L_diff, collapse will be prevented by turbulent diffusion.

Whether a particle is gravitationally bound to the cloud or disrupted by tidal shear originating from the central star may be assessed by comparing the cloud's self-gravitational force to the force of tidal gravity. Self-gravity dominates when

$$\begin{eqnarray}&&\displaystyle \frac{{{GM}}_{c}}{{r}^{2}}\gt 3\displaystyle \frac{{{GM}}_{* }}{{a}^{2}}\displaystyle \frac{r}{a},\end{eqnarray} \tag{ 10 }$$

where a is the separation between the cloud and the star. This criterion may be conveniently rewritten using the Hill-radius

$$\begin{eqnarray}&&{r}_{{\rm{H}}}=a{\left(\displaystyle \frac{{M}_{{\rm{c}}}}{{M}_{* }}\right)}^{1/3}.\end{eqnarray} \tag{ 11 }$$

The cloud will be stable against tidal shear if its density, ρ_c, is larger then its local Hill density, i.e., (Klahr & Schreiber 2016)

$$\begin{eqnarray}&&{\rho }_{{\rm{c}}}\gt {\rho }_{{\rm{H}}}:= {M}_{{\rm{c}}}{\left(\displaystyle \frac{4}{3}\pi {r}_{{\rm{H}}}^{3}\right)}^{-1}=\displaystyle \frac{9}{4\pi }\displaystyle \frac{{M}_{* }}{{a}^{3}}.\end{eqnarray} \tag{ 12 }$$

Using Equation (8), the cloud radius r for a given cloud column density ${{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}$ must not exceed

$$\begin{eqnarray}&&r\lt \displaystyle \frac{2\pi }{9}\displaystyle \frac{G{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}{{{\rm{\Omega }}}^{2}}=: {L}_{\mathrm{hill}}\end{eqnarray} \tag{ 13 }$$

for the cloud to be stable against tidal shear.

We combine the diffusion criterion in Equation (9) and the tidal shear criterion in Equation (13) to conclude that for local instability there must exist an r, such that L_diff < r < L_hill, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{\mathrm{St}}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{2\pi G{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}\lt r\lt \displaystyle \frac{2\pi }{9}\displaystyle \frac{G{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}{{{\rm{\Omega }}}^{2}}.\end{eqnarray} \tag{ 14 }$$

This is equivalent to

$$\begin{eqnarray}&&\displaystyle \frac{3}{2}\sqrt{\displaystyle \frac{\delta }{\mathrm{St}}}\lt \displaystyle \frac{\pi G{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}{{c}_{{\rm{s}}}{\rm{\Omega }}}.\end{eqnarray} \tag{ 15 }$$

The right-hand side of Equation (15) is very reminisent of the inverse Toomre Q parameter (Toomre 1964)

$$\begin{eqnarray}&&Q:= \displaystyle \frac{{c}_{{\rm{s}}}{\rm{\Omega }}}{\pi G{{\rm{\Sigma }}}_{{\rm{g}}}},\end{eqnarray} \tag{ 16 }$$

where Σ_g is the gas surface density. For Q < 1, the gas disk is subject to gravitational instability and fragmentation, since neither gas pressure nor tidal forces may prevent the gas disk's collapse (see, e.g., Baehr et al. 2017).

Next, we can define the total metallicity

$$\begin{eqnarray}&&Z:= \displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{p}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}},\end{eqnarray} \tag{ 17 }$$

where Σ_p and Σ_g are spatially averaged particle and gas column densities, respectively, i.e.,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{p}}}=\left\langle {\int }_{-\infty }^{\infty }{\rho }_{{\rm{p}}}{dz}\right\rangle \,\mathrm{and}\,{{\rm{\Sigma }}}_{{\rm{g}}}=\left\langle {\int }_{-\infty }^{\infty }{\rho }_{{\rm{g}}}{dz}\right\rangle ,\end{eqnarray} \tag{ 18 }$$

where ρ_p and ρ_g are particle and gas volume densities, respectively. We express the local metallicity of the particle clump as

$$\begin{eqnarray}&&{Z}_{{\rm{c}}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}{{{\rm{\Sigma }}}_{g}}=\displaystyle \frac{\epsilon {{\rm{\Sigma }}}_{{\rm{p}}}}{{{\rm{\Sigma }}}_{g}}=\epsilon Z.\end{eqnarray} \tag{ 19 }$$

The dimensionless quantity characterizes potential radial and azimuthal enhancement in particle column density due to, for example, streaming instability. We note that for the purposes of our paper, Z is defined by Equation (17) at each location in the disk and, thus, does not necessarily correspond to the stellar [Fe/H] and could in fact vary temporarily and spatially due to dust drift and growth processes. For example, the water ice line has been prominently featured as a location where drift or sublimation processes can produce enhanced local metallicities (e.g., Saito & Sirono 2011; Hyodo et al. 2019).

We can transform Equation (15) to

$$\begin{eqnarray}&&1\gt \displaystyle \frac{3}{2}\displaystyle \frac{Q}{\epsilon Z}\sqrt{\displaystyle \frac{\delta }{\mathrm{St}}}=: {Q}_{{\rm{p}}}.\end{eqnarray} \tag{ 20 }$$

Equation (20) is a criterion for instability. It neatly highlights the ambivalent effect of streaming instability on planetesimal formation, of both enhancing local particle densities (Johansen et al. 2007; Chiang & Youdin 2010; Simon et al. 2016), thus increasing while also providing the cloud with turbulent diffusion δ (Schreiber & Klahr 2018). Furthermore, high metallicities favor collapse, which is in line with, e.g., Youdin & Shu (2002) and Johansen et al. (2009). The critical value of Q_p = 1, which is equivalent to ${\rho }_{{\rm{c}}}={\rho }_{{\rm{H}}}$ and L_hill = L_diff, corresponds to the diffusion-limited collapse criterion of Klahr & Schreiber (2016, 2020) and the attached marginally critical length scale, given by

$$\begin{eqnarray}&&{r}_{\mathrm{crit}}=\displaystyle \frac{1}{3}\sqrt{\displaystyle \frac{\delta }{\mathrm{St}}}H.\end{eqnarray} \tag{ 21 }$$

This can be derived by setting the three terms in Equation (14) equal to each other. However, for Q_p < 1, more modes become unstable leading to a range of unstable scales that are both larger and smaller than the critical length scale. This may lead to the initial planetesimal size distributions seen in, e.g., Simon et al. (2016, 2017), Schäfer et al. (2017), Abod et al. (2019), and Nesvorny et al. (2019).

Metallicity Z and dominating Stokes number St can be informed from dust-evolution models that contain growth, drift, and fragmentation (e.g., Birnstiel et al. 2010, 2012) and are therefore often treated as input parameters in planetesimal formation models, as is the Toomre Q, which parameterizes disk mass. On the contrary, the diffusion coefficient δ and metallicity enhancement are local properties that cannot be treated as input parameters. However, both can be related directly to the scales regulating local particle densities in planetesimal-forming disks. These scales are set by drag instabilities, in particular, the particle KHI (Sekiya 1998; Johansen et al. 2006; Chiang 2008) and the nonlinear saturation of streaming instability (Youdin & Goodman 2005; Yang & Johansen 2014; Squire & Hopkins 2018a). Thus, investigating how these two instabilities interplay and regulate local enhancement and diffusion is crucial for understanding planetesimal formation in the scenario of local gravitational collapse.

Equations (28) and (29) can only describe particle and gas orbital velocities correctly as long as the two are not significantly subject to frictional coupling. For strongly coupled particles with St < 1, the equilibrium azimuthal particle velocity v_{0, y} of the combined fluid is (e.g., Nakagawa et al. 1986; Sekiya 1998; Youdin & Shu 2002)

$$\begin{eqnarray}&&{v}_{0,\phi }(a,z)={v}_{{\rm{K}}}\left(1-\displaystyle \frac{\eta }{1+\mu (z)}\right).\end{eqnarray} \tag{ 35 }$$

We note that Equation (35) depends on the local dust-to-gas ratio, μ, such that the azimuthal velocity limits to Equation (29) when gas dominates and to Equation (28) when dust dominates. The restriction St < 1 is valid in the context of the herein presented simulations and is also justified physically by dust-evolution models like those of Birnstiel et al. (2010, 2012).

Since particles are not vertically supported by a pressure gradient, they settle toward the mid-plane, resulting in a strongly z-dependent dust-to-gas ratio μ(z). In the mid-plane, where μ(z) is maximal, Equation (35) implies that both the particles and the gas will orbit at velocities close to Keplerian, whereas in layers above the mid-plane, their orbital velocities are reduced (by ηv_K in dust-free layers). The resulting vertical gradient in orbital velocity leads to KHIs, sometimes known as particle vertical shearing instabilities, that vertically mix and stir up particles (Sekiya 1998; Sekiya & Ishitsu 2000, 2001; Youdin & Shu 2002; Chiang 2008). We will follow Chiang (2008) in their analytical consideration of the KHI.

Whether or not the dust layer is subject to KHI is commonly assessed with the Richardson number Ri, which compares buoyancy oscillations (for incompressible perturbations given by the Brunt–Väisälä frequency) to the rate of vertical shearing (Chandrasekhar 1961). It is given by

$$\begin{eqnarray}&&\mathrm{Ri}:= \displaystyle \frac{\left({g}_{z}/\rho \right)\left(\partial \rho /\partial z\right)}{{\left(\partial {v}_{\phi }/\partial z\right)}^{2}},\end{eqnarray} \tag{ 36 }$$

where ${g}_{z}=-{{\rm{\Omega }}}^{2}z$ is the vertical component of stellar gravity, and $\rho ={\rho }_{{\rm{g}}}+{\rho }_{{\rm{p}}}$. If the Richardson number falls below a critical value, i.e., if

$$\begin{eqnarray}&&\mathrm{Ri}\lt {\mathrm{Ri}}_{\mathrm{crit}}=\displaystyle \frac{1}{4}\end{eqnarray} \tag{ 37 }$$

for Cartesian flows (Chandrasekhar 1961; Howard & Maslowe 1973; Li et al. 2003; Chiang 2008), KHIs are triggered and dust parcels are overturned and mixed. The exact value of the critical Richardson number is problem dependent. Johansen et al. (2006) found a critical value of Ri_crit ∼ 1 when numerically studying self-sustained KHI in protoplanetary disks.

As dust settles much closer to the mid-plane than gas, ∂ρ/∂z ≈ ∂ρ_p/∂z. Hence, the Richardson number in protoplanetary disks becomes (Chiang 2008)

$$\begin{eqnarray}&&\mathrm{Ri}=-\displaystyle \frac{{{\rm{\Omega }}}^{2}z}{{\rho }_{{\rm{p}}}+{\rho }_{{\rm{g}}}}\displaystyle \frac{\left(\partial {\rho }_{{\rm{p}}}/\partial z\right)}{{\left(\partial {v}_{\phi }/\partial z\right)}^{2}}=-\displaystyle \frac{1}{{\left(\eta a\right)}^{2}}\displaystyle \frac{{\left(1+\mu \right)}^{3}}{\partial \mu /\partial z}z,\end{eqnarray} \tag{ 38 }$$

where we plugged in Equation (35). For $\mathrm{Ri}=\mathrm{const}$, which is expected if settling times are much longer than KHI growth rates (Johansen et al. 2006), integration of Equation (38) yields a vertical profile for the dust-to-gas ratio μ(z) dependent on the mid-plane dust-to-gas ratio μ₀, i.e., (Chiang 2008)

$$\begin{eqnarray}&&\mu (z)={\left[\displaystyle \frac{1}{{\left(1+{\mu }_{0}\right)}^{2}}+\displaystyle \frac{1}{\mathrm{Ri}}\displaystyle \frac{{z}^{2}}{{\left(\eta a\right)}^{2}}\right]}^{-1/2}-1.\end{eqnarray} \tag{ 39 }$$

Solving μ(z) = 0 for z yields the maximum extent of the particle layer (Chiang 2008)

$$\begin{eqnarray}&&{z}_{\max }:= \sqrt{\mathrm{Ri}}\displaystyle \frac{\sqrt{{\mu }_{0}^{2}+2{\mu }_{0}}}{1+{\mu }_{0}}\eta a=\displaystyle \frac{\sqrt{\mathrm{Ri}}}{2}\displaystyle \frac{\sqrt{{\mu }_{0}^{2}+2{\mu }_{0}}}{1+{\mu }_{0}}\beta H.\end{eqnarray} \tag{ 40 }$$

Equation (40) predicts the vertical height of the particle layer for a given mid-plane dust-to-gas ratio and constant Richardson number. For μ₀ ≳ 1, z_max limits to $0.5\sqrt{\mathrm{Ri}}\beta H$.

Particles are expected to settle until the mid-plane is dense enough and the critical Richardson number is reached. Mixing and vertical excitation are now energetically favorable and counteract further settling, and thus, particles are lifted to z_max. When particles are stirred up, the Richardson number of the flow rises to subcritical values again such that the KHI cannot grow further. This leads to a self-regulating steady state where the flow Richardson number returns to the critical value, and the maximum extent is again given by Equation (40). As Equation (40) assumes a constant Ri-flow, we treat the steady-state flow Richardson number as a parameter and investigate which value best fits our numerical setup.

Note that this consideration does not account for the high-density mid-plane cusp found by Sekiya (1998), Youdin & Shu (2002), and Gómez & Ostriker (2005) for super-solar metallicities. Indeed, for sufficiently high metallicities, Youdin & Shu (2002) showed that this over-dense particle mid-plane can be massive enough to gravitationally fragment and form planetesimals (also see, e.g., Chiang & Youdin 2010). These enhancements in metallicity may be found in the inner disk due to radial drift of particles (Youdin & Shu 2002), in particle traps such as pressure bumps (e.g., Rice et al. 2006; Pinilla et al. 2012; Taki et al. 2016) or as the disk dissipates due to photoevaporation (Gorti et al. 2015).

However, none of these processes are thought to generically operate and sufficiently increase metallicities across a broad range of disk environments.

The streaming instability (Youdin & Goodman 2005), therefore, may offer a promising path of ubiquitously forming planetesimals (Simon et al. 2016; Nesvorny et al. 2019). This section briefly summarizes the nature of streaming instability, recent developments, and its importance in the context of the collapse criterion presented in Equation (20).

The streaming instability is a linear instability that occurs when the particles' back-reaction onto gas is taken into account (Youdin & Goodman 2005; Jacquet et al. 2011; Zhuravlev 2019). For marginally coupled particles and local dust-to-gas ratios exceeding unity, the SI displays rapid growth rates, thus, strongly enhancing local particle densities (Carrera et al. 2015; Yang et al. 2017) and triggering planetesimal formation (Bai & Stone 2010a; Johansen et al. 2015; Simon et al. 2016; Schäfer et al. 2017; Nesvorny et al. 2019). Whether or not this critical dust-to-gas ratio can be achieved depends on the metallicity and pressure gradient (Sekiya & Onishi 2018), which determine the vertical particle height set by KHI. The functionality of the SI can be understood in the context of resonant drag instabilities, a class of instabilities where the relative dust-gas motion is in resonance with a pressure perturbation in the gas (Squire & Hopkins 2018b). As showed by Squire & Hopkins (2018a), for the SI, the equilibrium solution for the relative dust-gas streaming velocity provided by Nakagawa et al. (1986) is in resonance with an epicyclic perturbation in the gas.

Recently, the role of the linear SI in the context of planetesimal formation was cast into doubt by two additional considerations that were omitted in the stability analysis by Youdin & Goodman (2005) and, likewise, in most numerical work. First, the inclusion of external turbulence strongly suppressed SI growth rates, in particular, on small scales (Umurhan et al. 2020). Second, Krapp et al. (2019) showed that the inclusion of multiple particle species also damps growth rates significantly. Still, in dead zones, where turbulence is entirely self-induced by SI and KHI, and if the particle size distribution is sufficiently top-heavy and narrow as suggested by recent dust-evolution models (e.g., Birnstiel et al. 2011, 2016; Okuzumi et al. 2012), SI remains an important processes and continues to be a key aspect of global models for planetesimal formation, such as that of Drążkowska et al. (2016).

Therefore, understanding the ambivalent role of SI for gravitational collapse is crucial. On the one hand, its linear phase and the subsequent onset of turbophoresis⁴ can concentrate particles and increase local densities to Hill density and beyond. On the other hand, the self-induced turbulence intrinsic to the nonlinear saturation of SI also diffuses particles apart and, thus, limits collapse on small scales (Klahr & Schreiber 2016, 2020). In fact, the fastest linearly growing modes are always on the smallest scales (Youdin & Goodman 2005; Squire & Hopkins 2018a; Umurhan et al. 2020) and, therefore, typically diffusion-limited. As the complex interplay of these oppositional processes during 3D nonlinear SI remains to be studied in detail, we propose a straightforward estimate for a characteristic scale of SI-induced over-densities.

Three-dimensional simulations of nonlinear SI universally feature large-scale azimuthally elongated filaments as linear modes are sheared apart (see, e.g., Johansen et al. 2012; Yang & Johansen 2014; Simon et al. 2016; Abod et al. 2019; Li et al. 2018, or Section 6). Previously, their typical separation (feeding zone) has been investigated numerically by Yang & Johansen (2014), and their radial extent has been connected to the operating scale of SI (Umurhan et al. 2020). In the following, we apply the Richardson criterion used to derive z_max for the KHI in Equation (40) to SI filaments to derive an analogous estimate for a characteristic radial scale.

To an order of magnitude, we may neglect Keplerian shear, and the radial gradient in azimuthal velocity follows from a radially dependent dust-to-gas ratio μ = μ(x) that, by construction, is present in SI filaments. Motivated by the findings of Squire & Hopkins (2018a), we take the restoring force per unit mass as g_x = −Ω²x from radial epicyclic oscillations. Under the assumption of a constant Richardson number, we can exactly match the argument presented in Section 3.1 and arrive at a maximum radial extent x_max for SI filaments given by

$$\begin{eqnarray}&&{x}_{\max }={z}_{\max }({\mu }_{0}\geqslant 1)=\displaystyle \frac{\sqrt{\mathrm{Ri}}}{2}\beta H=: {L}_{\max }.\end{eqnarray} \tag{ 41 }$$

We point out that while the simplification ${\mu }_{0}\geqslant 1$ is per definition true, as otherwise, the SI would not have produced an over-dense filament in the first place, the assumption of a constant Richardson number is likely not correct due to the asymmetry that arises when Keplerian shear is taken into account. On the scales of SI filaments, the linearized Keplerian shear velocity is ∼ηv_K, comparable to the velocity difference due to the particle density gradient, and adds to the relative dust-gas streaming velocity on one side of the filament and subtracts from the other. Also, as peak densities during the nonlinear SI often significantly exceed unity, it seems likely that an analogous structure to the high-density mid-plane cusp found by Sekiya (1998) forms, which is neither considered in the vertical nor the radial Richardson criterion. For these reasons, we stress that this radial scale has to be taken as a rough estimate.

Still, we show in Section 7 that having an analytic prediction for the radial scale, even if approximate, is very useful for retracing local diffusivity and enhancement to metallicity and pressure gradient and, thus, for evaluating our collapse criterion in Equation (20) without a detailed knowledge of the turbophoretic state of the nonlinear SI.

For all simulations, we consider a non-magnetized gas obeying an isothermal equation of state

$$\begin{eqnarray}&&P={c}_{{\rm{s}}}^{2}{\rho }_{{\rm{g}}}.\end{eqnarray} \tag{ 44 }$$

This is a good assumption for disk regions that are dominated by irradiation instead of accretion heating (Kratter & Murray-Clay 2011). Initially, the gas is evenly distributed in the radial and azimuthal directions, while vertically stratified due to the vertical component of stellar gravity causing mid-plane sedimentation. We achieve sub-Keplerian gas velocities by artificially imposing a radial pressure gradient to support the gas. This is quantified by β defined in Equation (32), which is added to the gas Euler equation (see Equation (46)). Boundary conditions are shear-periodic, periodic, and periodic in $x,y$, and z, respectively. We defer to Li et al. (2018) for a discussion of the effect of vertical boundary conditions on streaming instability.

We seed the box with Lagrangian super-particles following a Gaussian with a scale height of ${H}_{{\rm{p}},0}$. We test different values for ${H}_{{\rm{p}},0}$ and discuss our choices for subsequent simulation runs in Appendix B.1. Each super-particle is characterized by its position ${{\boldsymbol{x}}}_{{\rm{p}}}=({x}_{{\rm{p}}},{y}_{{\rm{p}}},{z}_{{\rm{p}}})$ and velocity ${\boldsymbol{v}}=({v}_{x},{v}_{y},{v}_{z})$ (measured with respect to ${u}_{0,y}$) and represents a swarm of identical physical solids interacting with the gas as a group. The triangular-shaped cloud scheme (Hockney & Eastwood 1988; Youdin & Johansen 2007) is applied to smooth out super-particle properties to the neighboring grid cells. Particle properties are characterized by the input parameters of particle Stokes number St from Equation (2) and disk metallicity Z from Equation (17). Note that the Pencil Code up-scales particle masses such that the total dust-to-gas ratio in the boxes matches Z even if the vertical domain size does not include the entire gas disk.

We consider particles with a fixed St. For streaming instability with a mixture of different particle species, we defer to Schaffer et al. (2018).

We also note that while the particles may display significant clumping in our simulations, the gas is comparatively static throughout all simulations.

In our Pencil Code setup, the gas obeys the continuity equation

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

where ${f}_{{\rm{D}}}\left({\rho }_{{\rm{g}}}\right)$ represents an artificial hyper-diffusivity. The gas velocity ${\boldsymbol{u}}$ relative to ${u}_{0,y}$ is evolved via the equation of motion

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}t}+\left({\boldsymbol{u}}\cdot {\rm{\nabla }}\right){\boldsymbol{u}}+{u}_{0,y}\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}y}=-{c}_{{\rm{s}}}^{2}{\rm{\nabla }}{\rm{l}}{\rm{n}}{\rho }_{{\rm{g}}}+{\rm{\Omega }}h\beta \hat{{\boldsymbol{x}}}\\ \,\,+\,\left(2{\rm{\Omega }}{u}_{y}\hat{{\boldsymbol{x}}}-\displaystyle \frac{1}{2}{\rm{\Omega }}{u}_{x}\hat{{\boldsymbol{y}}}-{{\rm{\Omega }}}^{2}z\hat{{\boldsymbol{z}}}\right)+\displaystyle \frac{\mu {\boldsymbol{w}}}{{t}_{{\rm{s}}}}+{f}_{\nu }\left({\boldsymbol{u}},{\rho }_{{\rm{g}}}\right).\end{array}\end{array}\end{eqnarray} \tag{ 46 }$$

Here, the terms from left to right are: the velocity time derivative, advection due to velocity perturbation, advection due to shear flow, pressure gradient, artificially imposed centrifugal support due to global radial pressure gradient, the combination terms from linearized stellar gravity, centrifugal force, Coriolis force, the back-reaction of the drag force exerted on gas by particles, and an artificial hyper-viscosity.

Hyper-viscosity and diffusivity are used to stabilize the code and smooth out steep gradients while preserving the power at the larger scales (see Appendix B of Yang & Krumholz 2012).

Similarly, the evolution equations for the particles read

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{x}}}{\partial t}=-q\omega {x}_{p}\hat{{\boldsymbol{y}}}+{\boldsymbol{v}}\end{eqnarray} \tag{ 47 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{\boldsymbol{v}}}{{\rm{\partial }}t}=\left(2{\rm{\Omega }}{v}_{y}\hat{{\boldsymbol{x}}}-\displaystyle \frac{1}{2}{\rm{\Omega }}{v}_{x}\hat{{\boldsymbol{y}}}-{{\rm{\Omega }}}^{2}z\hat{{\boldsymbol{z}}}\right)-\displaystyle \frac{{\boldsymbol{w}}}{{t}_{{\rm{s}}}}\end{eqnarray} \tag{ 48 }$$

with the key difference being that particles are not subject to pressure gradient.

Self-gravity is implemented by solving the Poisson equation with the fast Fourier transform algorithm (Gammie 2001) containing gravitational softening. We consider the gravitational potential ${\rm{\Phi }}({\boldsymbol{r}})$ at arbitrary position ${\boldsymbol{r}}$ of both gas and particles. The relative strength of tidal shearing compared to self-gravity is parameterized with the Toomre Q in Equation (16) (Toomre 1964).

This implies a value for the gravitational constant G in code units, which is then plugged into the right-hand side of the Poisson equation, i.e.,

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Phi }}({\boldsymbol{r}})=4\pi G\rho ({\boldsymbol{r}})=\sqrt{\displaystyle \frac{8}{\pi }}\displaystyle \frac{{{\rm{\Omega }}}^{2}}{Q{\rho }_{{\rm{g}},0}}\rho ({\boldsymbol{r}}).\end{eqnarray} \tag{ 49 }$$

Note that Q relates to the self-gravity parameter $\hat{G}$ used by Johansen et al. (2012), Simon et al. (2016), and Schäfer et al. (2017) to parameterize self-gravity via

$$\begin{eqnarray}&&\hat{G}=\sqrt{\displaystyle \frac{8}{\pi }}\displaystyle \frac{1}{Q}\approx \displaystyle \frac{2}{Q}.\end{eqnarray} \tag{ 50 }$$

Lastly, we note that in our simulations, even for Q < 1, the gas does not collapse as the chosen domain size does not include the critical unstable wavelength of 2πH. In Section 8.1, we briefly discuss the applicability of our work to disks subject to the Toomre instability.

In our simulations, we set the code units for angular frequency Ω = 1, isothermal sound speed c_s = 1, and initial mid-plane gas density ${\rho }_{{\rm{g}},0}=1$. While the Pencil Code is agnostic to the choice of unit, we deem it useful to express our results in units that have physical meaning. Thus, we choose the orbital period, which is related to the dynamical time via 2πΩ⁻¹, the gas scale height H, and the initial mid-plane gas density ρ_g,0 as time, length, and density units, respectively. As the orbital frequency Ω is given by

$$\begin{eqnarray}&&{\rm{\Omega }}=\sqrt{\displaystyle \frac{{{GM}}_{* }}{{a}^{3}}},\end{eqnarray} \tag{ 51 }$$

we can calculate scaling relations of our code units, i.e.,

$$\begin{eqnarray}&&2\pi {{\rm{\Omega }}}^{-1}=\cdot {\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{\tfrac{3}{2}}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\mathrm{sun}}}\right)}^{-\tfrac{1}{2}}\mathrm{yr}\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&H=0.027\cdot {\left(\displaystyle \frac{T}{180{\rm{K}}}\right)}^{\tfrac{1}{2}}{\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{\tfrac{3}{2}}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\mathrm{sun}}}\right)}^{-\tfrac{1}{2}}\mathrm{au},\end{eqnarray} \tag{ 53 }$$

where we used the isothermal sound speed given by

$$\begin{eqnarray}&&{c}_{{\rm{s}}}=\sqrt{\displaystyle \frac{{k}_{{\rm{B}}}T}{\bar{m}}},\end{eqnarray} \tag{ 54 }$$

with Boltzmann constant k_B and a mean molecular weight of $\bar{m}=2.33{m}_{{\rm{p}}}$, where m_p is the proton mass.

While the density scales freely with ${\rho }_{{\rm{g}},0}$ for simulations where self-gravity is disabled, Equation (16) introduces a scaling relation once self-gravity is introduced. Thus, the density unit can be expressed via

$$\begin{eqnarray}&&{\rho }_{{\rm{g}},0}=3\cdot {10}^{-7}\cdot {\left(\displaystyle \frac{Q}{1}\right)}^{-1}{\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{-3}\left(\displaystyle \frac{{M}_{* }}{{M}_{\mathrm{sun}}}\right){\rm{g}}\,{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 55 }$$

We conduct multiple simulations with various different parameter choices, displayed in Table 1. We follow Simon et al. (2016) and choose a domain size of L_x = L_y = L_z = 0.4H, which allows us to sufficiently capture multiples of the characteristic scale βH. The majority of our simulations are based on a relatively small numerical grid with a size of 64 × 64 × 64, thus, sacrificing resolution in favor of computational time. Additionally, we explore a higher resolution of 128 × 128 × 128 for selected simulation setups. We denote the number of grid cells in the radial, azimuthal, and vertical directions with N_x, N_y, and N_z, respectively.

For these choices, the grid cell sizes are $\sim 0.006H$ and ∼0.003H, respectively, which, for β = 0.1 and typical viscosity values, may not be enough to sufficiently resolve the fastest growing streaming instability wavelength (Umurhan et al. 2020). However, the marginally unstable scale for local gravitational collapse r_crit in Equation (21) (Klahr & Schreiber 2020) is for typical values just resolved (see Section 5.2). Nevertheless, the fact that our simulations, in particular, those with low resolution, are likely not fully converged is to be kept in mind when evaluating our numerical results. We discuss the implications in Appendix. B.2. For better-resolved studies of pure streaming instability as well as gravitational collapse, we refer to, e.g., Yang et al. (2017), Schreiber & Klahr (2018), Sekiya & Onishi (2018), Nesvorny et al. (2019), and Klahr & Schreiber (2020), respectively.

For low-resolution simulations, the grid is seeded with one super-particle per grid cell, on average. High-resolution simulations are seeded with 0.1 super-particles per grid cell to save on computation time. As particles are expected to settle to a layer of thickness ∼0.04H after Equation (40), corresponding to about a tenth of the domain size, the mid-plane layer will nevertheless have about one super-particle per grid cell.

To appropriately relate particle scale heights to z_max, it is worthwhile to first check whether the assumption of a constant Richardson number is appropriate for our system. Due to azimuthal and radial symmetry of the settling simulations, particle positions are projected onto the vertical axis. We calculate the Richardson number at height z after Equation (38) via

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Ri}(z) & = & -\displaystyle \frac{4z}{{\beta }^{2}{H}^{2}}\displaystyle \frac{{\left[1+\mu \left(z\right)\right]}^{3}}{\partial \mu /\partial z}\\ & \approx & -\displaystyle \frac{4z}{{\beta }^{2}{H}^{2}}{\left[1+\mu \left(z\right)\right]}^{3}\displaystyle \frac{{\rm{\Delta }}z}{{\rm{\Delta }}\mu }.\end{array}\end{eqnarray} \tag{ 58 }$$

This calculation is limited to our resolution. As such, Δμ is the difference in the dust-to-gas ratio of two vertically adjacent slices, Δz = L_z/N_z, and z-positions of interest fulfill $z\in n\cdot {L}_{z}/{N}_{z}$ with ${N}_{z}\gt n\in {\mathbb{N}}$. To avoid statistical effects, we only calculate the Richardson number at a given z if the two adjacent slices contain more particles than some threshold value, which was chosen to be 500.

The height-dependent Richardson number after Equation (58) for the final snapshot of high-resolution settling simulations (N_x,y,z = 128, St = 0.005) for three different metallicities Z = 0.0002, Z = 0.002, and Z = 0.02 is shown in the middle panel of Figure 2. The top panel depicts the number of particles per slice in z. Note that although the mid-plane contains about three times as many particles for Z = 0.0002 compared to Z = 0.02, the mid-plane dust-to-gas ratio is smaller, as individual super-particle masses are lower by two orders of magnitude.

Zoom In Zoom Out Reset image size

Although the Z = 0.0002 simulation is not resolved well enough to convincingly determine a trend, the vertical profile of the flow Richardson number for Z = 0.002 and Z = 0.02 qualitatively matches results by, e.g., Johansen et al. (2006) and Bai & Stone (2010a). Our estimate of the Richardson number also supports findings by Johansen et al. (2006) that the appropriate critical Richardson number is indeed greater than the literature value of one-fourth. Thus, while we continue to treat the Richardson number as a parameter, we deem the assumption of it being constant required to derive z_max in Equation (40) reasonable.

Next, we relate the vertical particle extent z_max to vertical diffusivity and show than rms velocities of particles behave accordingly. For a discussion of the radial diffusivity and our assumption of spherically symmetric diffusion, see Appendix A.

The dimensionless diffusion coefficient δ is commonly measured using the rms velocity, the vertical component of which is given by (Carrera et al. 2015; Schreiber & Klahr 2018)

$$\begin{eqnarray}&&{v}_{\mathrm{rms},{\rm{z}}}=\sqrt{\langle {v}_{z}^{2}\rangle }=\sqrt{\displaystyle \frac{1}{{N}_{\mathrm{par}}}\sum _{i}{v}_{z,i}^{2}},\end{eqnarray} \tag{ 59 }$$

where v_{z, i} is the vertical component of the velocity of particle i. It relates to the dimensionless diffusion coefficient via (Johansen et al. 2006b; Schreiber & Klahr 2018)

$$\begin{eqnarray}&&{\delta }_{{\rm{z}}}={\tau }_{\mathrm{corr}}{\left(\displaystyle \frac{{v}_{\mathrm{rms},{\rm{z}}}}{{c}_{{\rm{s}}}}\right)}^{2},\end{eqnarray} \tag{ 60 }$$

where τ_corr is the dimensionless correlation time. For St > 1 particles, this correlation time is just given by the Stokes number St (Youdin & Lithwick 2007). However, in our case, for well-coupled particles with St < 1, particle dispersion does not originate from random movement but from gas turbulence. Thus, eddy turnover times will dictate correlation times (Youdin & Lithwick 2007). For Kolmogorov turbulence, where particle turbulent kinetic energy is dominated by the largest eddies, whose turnover time is just the orbital time, then τ_corr = 1. In numerical simulations, however, this is often not a good assumption. For example, Schreiber & Klahr (2018) found values of 0.1 < τ_corr < 1 depending on their simulation setup. Hence, we deem it worthwhile to measure appropriate correlation times for our simulation setup.

The equating settling time in Equation (56) and diffusion time in Equation (6) yield a well-known expression for the particle scale height (Youdin & Lithwick 2007)

$$\begin{eqnarray}&&{H}_{{\rm{p}}}=\sqrt{\displaystyle \frac{{\delta }_{z}}{\mathrm{St}}}H.\end{eqnarray} \tag{ 61 }$$

In fact, Equation (61) is the definition of ${\delta }_{z}$ (and via our isotropy assumption, also δ_x) in the context of our work. Plugging in Equation (60) yields an expression for the dimensionless correlation time

$$\begin{eqnarray}&&{\tau }_{\mathrm{corr}}=\mathrm{St}\cdot {\left(\displaystyle \frac{{H}_{{\rm{p}}}/H}{{v}_{\mathrm{rms},{\rm{z}}}/{c}_{{\rm{s}}}}\right)}^{2}.\end{eqnarray} \tag{ 62 }$$

We determine v_rms according to Equation (59) and measure the particle scale height H_p using (Carrera et al. 2015)

$$\begin{eqnarray}&&{H}_{{\rm{p}}}=\sqrt{\langle {z}^{2}\rangle }=\sqrt{\displaystyle \frac{1}{{N}_{\mathrm{par}}}\sum _{i}{\left|{z}_{i}-\langle z\rangle \right|}^{2}},\end{eqnarray} \tag{ 63 }$$

where z_i is the vertical position of particle i, and $\langle z\rangle$ is its mean over all i (here, $\langle z\rangle =0$).

The dimensionless correlation times required to match rms velocities of the three high-resolution settling simulations at the final snapshot to the respective particle scale heights are ${\tau }_{\mathrm{corr}}(Z=0.0002)\approx 2$, τ_corr(Z = 0.002) ≈ 0.3, and τ_corr(Z = 0.02) ≈ 0.1. This is in line with typical correlation times measured by, e.g., Schreiber & Klahr (2018).

The top panel of Figure 3 plots the final snapshot correlation times versus metallicity for both low- and high-resolution simulations. The panel also distinguishes between simulations where the correlation time was able to reach a constant, non-evolving value within the simulation run time (steady state) and simulations where the correlation time was still evolving (no steady state), i.e., the particles are still in the settling phase, and the flow has likely not yet reached the KHI-regulated equilibrium. The two simulations with super-solar metallicities of Z = 0.1 and Z = 0.2, which develop the high-density mid-plane cusp, are excluded from this classification.

Zoom In Zoom Out Reset image size

Having verified that (vertical) rms velocities correlate with the (vertical) diffusivity and, thus, the particle scale height as anticipated, we proceed with Equation (61) and use H_p to measure diffusivity.

The bottom panel of Figure 3 plots the final snapshot particle scale height versus metallicity Z. In order to overlay the analytic prediction from the self-regulated KHI steady state, we require a measurement for the mid-plane dust-to-gas ratio μ₀. We find that the approximation μ₀ ≈ (H/H_p)Z well matches our numerical result for μ₀ (also compare to Equation (3) of Sekiya & Onishi 2018).

Moreover, to account for the fact that z_max, by construction, is the maximum particle extent, we must divide by an order unity factor to yield a quantity that relates to the scale height. The bottom panel of Figure 3 suggests that z_max/e is a good match for low and intermediate metallicities that do not possess a high-density mid-plane cusp.

Together with Equation (61), H_p ≈ z_max/e implies

$$\begin{eqnarray}&&\sqrt{\displaystyle \frac{{\delta }_{{\rm{z}}}}{\mathrm{St}}}=\displaystyle \frac{{H}_{{\rm{p}}}}{H}\approx \displaystyle \frac{1}{e}\displaystyle \frac{{z}_{\max }}{H}=\displaystyle \frac{\sqrt{\mathrm{Ri}}}{2e}\displaystyle \frac{\sqrt{{\mu }_{0}^{2}+2{\mu }_{0}}}{1+{\mu }_{0}}\beta .\end{eqnarray} \tag{ 64 }$$

Note, again, that the first equality in Equation (64) is just the definition of ${\delta }_{z}$ in the context of our work. Relating δ_z to rms velocities requires knowledge of correlation times, which we will defer to future work. The second equality relates the measured particle scale height to its analytic prediction from Chiang (2008), which was reviewed in Section 3.1.

We can combine Equation (64) with Equation (21) to find that, under the assumption of spherically symmetric diffusion (see Appendix A), the critical radius scale is just

$$\begin{eqnarray}&&{r}_{\mathrm{crit}}\approx \displaystyle \frac{1}{3}{H}_{{\rm{p}}}=\displaystyle \frac{\sqrt{\mathrm{Ri}}}{6e}\displaystyle \frac{\sqrt{{\mu }_{0}^{2}+2{\mu }_{0}}}{1+{\mu }_{0}}\beta H,\end{eqnarray} \tag{ 65 }$$

which, for β = 0.1 and typical values for μ₀ ≳ 1 and Ri ≈ 1 (Johansen et al. 2006, or Figure 2), is r_crit ≈ 0. 06β H = 0.006H and, thus, just resolved for the N_x,y,z = 64 simulations (also see Appendix B.2).

The top panels of Figures 4 and 5 display azimuthally integrated particle densities. The wave-like structure characteristic to diagonal streaming instability wave modes emerges (compare to, e.g., Li et al. 2018). As this feature needs to be accounted for when calculating the particle scale height, ${H}_{{\rm{p}},{\rm{j}}}$ is calculated via Equation (63) at each slice j in x (N_x slices with radial extent L_x/N_x) and then averaged, i.e.,

$$\begin{eqnarray}&&{H}_{{\rm{p}}}=\displaystyle \frac{1}{{N}_{x}}\sum _{j}{H}_{{\rm{p}},{\rm{j}}}.\end{eqnarray} \tag{ 66 }$$

The resulting final snapshot particle scale heights are depicted in Figure 6 against pressure gradient β and Z. The high-resolution settling simulations are plotted for comparison. The top panel of Figure 6 shows an approximately linear scaling of H_p in β, as expected from Equation (40). For the trivial case of β = 0, the relative dust-gas velocity is zero, and both KHI and SI lose their energy source such that razor-thin settling occurs. The bottom panel of Figure (6) reproduces the bottom panel of Figure 3. For low and intermediate metallicities, z_max describes H_p well, even in the presence of SI. For Z > 0.01, the high-density mid-plane cusp develops, and the particle scale height decreases below the expected height. Thus, we expect the approximation in Equation (64) to hold both in the presence and absence of streaming instability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The bottom panels of Figures 4 and 5 display vertically integrated particle densities. Over-dense streaming instability filaments develop for Z ≥ 0.02 and 0 < β ≤ 0.13, which is in line with, e.g., Bai & Stone (2010b).

To quantify the width of a filament, we define the slice density as the vertically and azimuthally integrated density, i.e.,

$$\begin{eqnarray}&&{\widetilde{{\rm{\Sigma }}}}_{{\rm{p}}}(x):= \int {{\rm{\Sigma }}}_{{\rm{p}}}(x,y){dy}.\end{eqnarray} \tag{ 67 }$$

The temporal evolution of the radial distribution of the slice density for the high-resolution fiducial run is shown in Figure 7. The slope of the filaments correlates with the local dust-to-gas ratio according to Nakagawa et al. (1986). We also identify radial epicyclic oscillations, especially at early times.

Zoom In Zoom Out Reset image size

The radial extent of the a filament Δx can be defined as

$$\begin{eqnarray}&&{\rm{\Delta }}x={x}_{+}-{x}_{-},\end{eqnarray} \tag{ 68 }$$

with x₋ < x < x₊ such that ${\widetilde{{\rm{\Sigma }}}}_{{\rm{p}}}({x}_{-}),{\widetilde{{\rm{\Sigma }}}}_{{\rm{p}}}({x}_{+})$<${\widetilde{{\rm{\Sigma }}}}_{{\rm{p}},0}\lt {\widetilde{{\rm{\Sigma }}}}_{{\rm{p}}}(x)$ holds true for all x₋ < x < x₊. Here, ${\widetilde{{\rm{\Sigma }}}}_{{\rm{p}},0}$ is the mean particle slice density per grid cell extent x_cell = L_x/N_x. Note that we must limit the radial resolution of this procedure to the gas grid cell size to remain self-consistent.

The top panel of Figure 8 shows the radial extent of the most dense filament at the final snapshot versus pressure gradient β for those simulations that developed distinct streaming instability filaments with Z = 0.02, i.e., FidRun_128, Beta_128 with $\beta =0.05$, FidRun_64, and Beta_64 and β ∈ {0.05, 0.07, 0.1, 0.13} (see also Table 1, and Figures 4 and 5). The numerical data well match L_max in Equation (41). To the extent that the vertical and radial scale heights of over-dense particle filaments are set by the same scale, Equation (61) implies that vertical and radial diffusivity must also be approximately equal and the approximation in Equation (64) is valid for radial diffusivity as well.

Zoom In Zoom Out Reset image size

The resulting dimensionless enhancement parameter defined in Equation (19) is estimated as

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{p}},{\rm{c}}}}{{{\rm{\Sigma }}}_{{\rm{p}}}}\approx {\left(\displaystyle \frac{1}{{\rm{\Delta }}x{\widetilde{{\rm{\Sigma }}}}_{{\rm{p}},0}}{\int }_{x-}^{{x}_{+}}{\widetilde{{\rm{\Sigma }}}}_{{\rm{p}}}(x){dx}\right)}^{p},\end{eqnarray} \tag{ 69 }$$

which, by construction, is larger than unity. Here, we introduced a correction coefficient p ≥ 1.

While the term in parenthesis assumes an azimuhthally isotropic filament, with radially uniformly distributed material, the correction exponent p accounts for potential additional azimuthal enhancement, as well as the fact that a radially uniform distribution tends to underestimate the local enhancement. Note that while the latter effect could be corrected via a multiplicative factor, the former requires an exponential coefficient. As it is challenging to account for these corrections rigorously, we treat them with a single correction exponent and set it to best match our simulation results. This procedure should be regarded as a first step and could be refined in future assessments of the collapse criterion.

The bottom panel of Figure 8 shows the evolution of for the three simulations, for which self-gravity is turned on in Section 7 at 40 orbits (FidRun_64, Z = 0.02, and LowZ_128, Z = 0.01) or 60 orbits (FidRun_128, Z = 0.02). We set p = 2.5 for the two fiducial runs with Z = 0.02. Our methodology of measuring enhancement in Equation (69) is specifically designed for enhancements within azimuthally elongated filaments. As such, in setups where those are not distinctly visible, as is the case for Z = 0.01, our procedure does not work, since the projection onto the x-axis removes all information about density enhancements that are rotated in the x–y-plane. Thus, we use p = 4.5 for the low-metallicity run with Z = 0.01, which best reproduces locally measured dust-to-gas ratios, to still get comparable values and assess the collapse criterion (see Section 7).

As expected for Z = 0.01, the local enhancement never exceeds an order unity factor of  ∼ 5. This is in line with Carrera et al. (2015) and Yang et al. (2017), who concluded that streaming instability cannot effectively concentrate particles for this metallicity. For Z = 0.02, the enhancement parameter fluctuates around 10 ≲ ≲ 100. There is no qualitative difference between the two tested resolutions to be found. We attribute the peak in at 55 orbits for FidRun_128 to statistical fluctuations typical to the nonlinear phase of streaming instability.

We numerically showed, albeit only for a limited range of parameters, that local gravitational collapse to form a planetesimal requires the diffusion- and tidal-shear-limited collapse criterion presented in Equation (20) to hold true. Applying the criterion to global disk models is not obvious, as both and δ are local properties that are potentially dependent on resolution and included physics. However, in Section 5, we showed that the ratio δ/St can be related to L_max via Equation (64). Thus, the collapse criterion in Equation (20) can be expressed as

$$\begin{eqnarray}&&1\gtrsim \displaystyle \frac{3}{2}\displaystyle \frac{Q}{\epsilon Z}\displaystyle \frac{\sqrt{\mathrm{Ri}}}{2e}\displaystyle \frac{\sqrt{{\mu }_{0}^{2}+2{\mu }_{0}}}{1+{\mu }_{0}}\beta \approx \displaystyle \frac{1}{10}\displaystyle \frac{Q}{\epsilon }\displaystyle \frac{\beta }{Z}={Q}_{{\rm{p}}}.\end{eqnarray} \tag{ 71 }$$

If Q_p ≲ 1, i.e., for small pressure gradients and Q values, and for large metallicities and local enhancements, we expect collapse and planetesimal formation. This Toomre-like formulation of the collapse criterion highlights that cloud collapse and planetesimal formation do not fundamentally depend on the presence of streaming instability, and interpreting streaming instability as a required catalyst for planetesimal formation is not appropriate. Instead, collapse is assessed by comparing disk mass (quantified by Q), local dust content Z, and pressure gradient β. In fact, Sekiya & Onishi (2018) recently showed that combining pressure gradient and metallicity to a single parameter appropriately describes streaming instability clumping, so it is no surprise that this ratio also appears in our collapse criterion.

The criterion in Equation (71) can also be utilized to assess if it is easier to collapse the gas disk or local particle clumps. If

$$\begin{eqnarray}&&\displaystyle \frac{{Q}_{{\rm{p}}}}{Q}\approx \displaystyle \frac{1}{10}\displaystyle \frac{\beta }{\epsilon Z}\lesssim 1,\end{eqnarray} \tag{ 72 }$$

particle collapse is easier. For our fiducial pressure gradient value of β = 0.1, and for no local enhancements ( = 1), this would be the case for Z ≳ 0.01.

We further deem it worth commenting on the applicability of our criterion to gas disks that form spiral arms due to gravitational instability. As seen in, e.g., Baehr & Klahr (2019), simulations able to resolve the gravitational instability of the gas and with Q < 1 can produce particle-rich regions in gas over-densities (such as spiral arms) on the scale of the domain size of our simulations. Whether such particle concentrations are present in gravitoturbulent Q > 1 disks and if they are subject to the criterion of this paper is the focus of ongoing research. As such, even if Q_p/Q > 1 initially, active gravitational instability may easily provide the right conditions for planetesimal formation to proceed, as Q is small by construction, and gas collapse/concentration can locally enhance Z.

By replacing diffusivity with particle scale height, one also eliminates the explicit dependency of the collapse criterion on Stokes number. We point out that the Stokes number still implicitly effects the collapse criterion by influencing the local enhancement due to streaming instability. Our work is most appropriate for scenarios of narrow, top-heavy particle size distributions. The validity of the collapse criterion for different Stokes numbers as well as for a particle size distribution (compared to, e.g., Schaffer et al. 2018), which may strongly damp SI growth rates (Krapp et al. 2019), remains to be studied. Similarly, we defer a high-resolution test to future work (also see Appendix B.2).

Although the physical origin of KHI investigated in Chiang (2008) and in our work is vertical shear, we explicitly abstain from referring to it as vertical shear instability. This term is commonly used for a gas-only instability active in protoplanetary disks that are able to cool sufficiently fast to allow for the development of vertical gas perturbations (see, e.g., Urpin & Brandenburg 1998; Nelson et al. 2013; Lin & Youdin 2015; Pfeil & Klahr 2019). Here, the vertical gradient in orbital frequency is not caused by dust-gas coupling but instead by radial variations in temperature and entropy (Barker & Latter 2015), leading to vertically elongated modes (Arlt & Urpin 2004; Klahr et al. 2018). The vertical shear instability is explicitly not studied in our work, as we are exclusively considering locally isothermal scenarios. We refer to Schäfer et al. (2020), who study the effect of vertical shear instability on a particle layer subject to SI.

Further, thermally driven (subcritical) baroclinic instabilities (Klahr & Bodenheimer 2003; Petersen et al. 2007a, 2007b; Klahr et al. 2018) and their linear phase, the convective overstability (Klahr & Hubbard 2014; Lyra 2014; Latter 2016), are based on radially perturbed gas parcels carrying entropy and have previously been shown to affect the particle mid-plane layer. There, inward moving gas radiatively thermalizes with its environment, which leads to buoyancy oscillations that in turn transport entropy outward. This mechanism was shown to form vortices that can trap particles and, thus, significantly alter the particle distribution in the mid-plane layer (e.g., Klahr et al. 2018).

The gas vertical shear instability, the subcritical baroclinic instability, and the convective overstability, as well as others not yet mentioned, like the magneto-rotational instability (MRI; see, e.g., Balbus & Hawley 1991, 1998; Davis et al. 2010; Bai 2011; Béthune et al. 2016), or zombie vortex instabilities (Marcus et al. 2015, 2016; Lesur & Latter 2016; Umurhan et al. 2016), can introduce additional turbulence (Klahr et al. 2018), which, as shown by Umurhan et al. (2020), reduces growth rates of SI modes and, thus, may render particle clumping and planetesimal formation via SI-induced over-densities more challenging (also see Chen & Lin 2020).

Thus, while the diffusion-limited collapse criterion from Klahr & Schreiber (2016, 2020) and Equation (20) is generally applicable for the gravitational collapse of over-dense clumps, the herein presented numerical verification thereof and, specifically, the expression in Equation (71) are most appropriate for protoplanetary disk MRI dead zones. In particular, our criteria apply when δ ≳ α, where the level of turbulence generated by other processes is given by the commonly used parameter α (Shakura & Sunyaev 1973) such that the turbulent velocity dispersion of the gas is ${v}_{\mathrm{turb}}={\alpha }^{1/2}{c}_{{\rm{s}}}$ (Cuzzi et al. 2001). For the simulations in this paper, δ ∼ 10⁻⁷–10⁻⁵ depending on metallicity (see Figures 2, 6) and Stokes number after Equation (61). Gas and particle diffusivity are related via the dimensionless Schmidt number Sc = α/δ (e.g., Youdin & Johansen 2007), which, for typical Schmidt numbers of order unity (see, e.g., Schreiber & Klahr 2018), leads to weaker gas turbulence than typically produced in, e.g., MRI simulations, where α can reach values up to ∼10⁻³ (see, e.g., Bai 2015).

Recently, the DSHARP survey (see, e.g., Andrews et al. 2018) and, likewise, the Ophiuchus disk survey employing ALMA (Cieza et al. 2019) showed that rings are an abundant substructure in protoplanetary disks, which may hint at the existence of planets carving gaps in the gas disk, which leads to a pressure bump collecting the radially drifting dust particles (Pinilla et al. 2012). Recently, Stammler et al. (2019) were also able to explain the observed rings as a byproduct of ongoing planetesimal formation. Regardless, it seems likely that pressure bumps constitute a common phenomena in protoplanetary disks (Dullemond et al. 2018), and one must therefore ask about the applicability of our collapse criterion therein. Pressure bumps lead to a local decrease in pressure gradient β. As a result, the relative dust-gas velocity decreases, and KHI weakens, leading to a thinner particle layer. As this effect is reflected in the proportionality of Q_p in β in Equation (71), we expect pressure bumps to favor planetesimal formation, which agrees with findings by Dittrich et al. (2013) who investigated planetesimal formation in MRI-induced zonal flows and pressure bumps.

Lastly, we acknowledge that the evolution of formed clumps has not been studied in our work. This includes planetesimal growth, mergers, and accretion (e.g., Kokubo & Ida 2012; Johansen & Bitsch 2019; Liu et al. 2019; San Sebastián et al. 2019), possible fragmentation (Wakita et al. 2017; Gerbig et al. 2019), migration (e.g., Goldreich & Tremaine 1979; Murray-Clay & Chiang 2006; Kley & Nelson 2012), pebble accretion (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012; Bitsch et al. 2015; Rosenthal & Murray-Clay 2019), and fragmenting and gravitoturbulent disks (Gibbons et al. 2012; Booth & Clarke 2016; Baehr & Klahr 2019).

Our results show that particle concentration by streaming instability and planetesimal formation through gravitational collapse are not equivalent. As such, our simulations were able to produce planetesimals for a metallicity of Z = 0.01 for Q ≲ 2, even though no prior streaming instability clumping occurred and, likewise, did not fragment for high Q values regardless of significant prior SI-induced particle concentration.

Instead, the diffusion-limited collapse criterion in Equation (71) implies that whether or not planetesimals form depends only on disk mass via the Toomre parameter Q, the ratio of pressure gradient and metallicity β/Z (compare to Sekiya & Onishi 2018), and the local enhancement in dust content .

The herein presented Toomre-like formulation of the diffusion-limited collapse criterion also may be able to shed new light on the initial planetesimal mass function. Equation (14) predicts a range of length scales subject to nonlinear gravitational instability, determined by the value of the right-hand side of Equation (20) (or Equation (71)). By studying the dispersion relation in more detail, one can attach growth rates to each unstable scale, which can then be used to analytically inform an initial planetesimal mass function.

In past numerical studies of the initial mass function, e.g., Simon et al. (2016), Schäfer et al. (2017), and Nesvorny et al. (2019), the enhancement is the only one of the quantities determining cloud collapse that is evolving in time. Therefore, the range of unstable scales and also initial planetesimal masses in past numerical studies may have been predominantly affected by parameter choices instead of characteristic properties of streaming and KHI. This interpretation could explain why Simon et al. (2016) formed planetesimals when self-gravity was active from the start—the collapse criterion in Equation (20) was already fulfilled by the initial condition—and why they found that higher self-gravity parameters lead to more massive planetesimals: more scales are unstable for lower Q values. Real disks however, are not expected to start off unstable to gravitational collapse, so one must ask how numerical simulations can inform initial planetesimal masses more physically.

Our collapse criterion in Equation (71) offers a possible solution to this conundrum.

In principle, there are four ways a system can evolve into fulfilling the collapse criterion. However, Q and global pressure gradient β in the absence of pressure bumps (see Section 8.2) tend to be non-evolving. As such, softly turning on self-gravity as performed by Schäfer et al. (2017) does not necessarily reflect real conditions in protoplanetary disks, in particular, since the strength of self-gravity does not explicitly affect SI growth rates. This is in contrast to the metallicity and Stokes number, which both highly influence streaming instability clumping (e.g., Carrera et al. 2015) but also inform the collapse criterion directly (for the metallicity) and implicitly via the enhancement parameter .

Thus, to prevent the initial condition from collapsing, a simulation could be set up with a temporally increasing metallicity and having self-gravity turned on from the start. Then, planetesimals will only be formed once the collapse criterion is fulfilled. For high-mass disks, we expect the mid-plane to start collapsing before the metallicity is high enough for streaming instability to start (compare to our Z = 0.01 simulation). This is of particular interest in the light of recent results by Powell et al. (2019) suggesting that disks may tend to be more massive than previously appreciated.

On the other hand, for low-mass disks, we expect the onset of streaming instability to occur before the collapse criterion is fulfilled. In this case, streaming instability can locally increase particle densities, such that collapse occurs as soon as is high enough. A similar effect can be achieved by gradually increasing the Stokes number of super-particles. Both approaches mimic drift and growth processes in the sense that the abundance of particles with Stokes numbers appropriate for streaming instability increases over time.

Regardless of if the entire mid-plane or only locally enhanced regions collapse, we expect that, unless the growth timescale of the perturbation is shorter than its collapse timescale, marginally unstable modes will collapse first (i.e., for the right-hand side of Equation (20) equal to unity). This would prevent the collapse of larger (or smaller) scales and dictate a single characteristic and dominant planetesimal size determined by the critical length scale in Equation (21) (see Klahr & Schreiber 2016, 2020) and, therefore, potentially lead to a initial planetesimal mass function that is qualitatively different to the power law seen in, e.g., Johansen et al. (2015) and Simon et al. (2016). Instead, the resulting initial mass function may be similar to those derived via cascade models of turbulent concentration (Cuzzi et al. 2010; Hartlep & Cuzzi 2020). We defer the detailed investigation of the planetesimal initial mass function in the context of the herein presented criterion to future work.

To conclude, our analytical considerations and numerical simulations suggest that the gravitational collapse of over-dense clouds equally depends on metallicity, pressure gradient, Toomre Q, and potential local enhancements. These quantities should therefore be evaluated together, when assessing planetesimal formation. Streaming and KHI are crucial in regulating local over-densities. Our collapse criterion remains to be tested for higher resolutions, different Stokes numbers, and pressure gradients and for a larger range of metallicities.

To verify that the vertical extent of the particle layer is independent of initial dust distribution, we performed particle settling tests with varying widths of the initial Gaussian distribution. For all simulations, we chose N_x,y,z = 64, Z = 0.02, St = 0.01, and β = 0.1. We therefore expect a settling time of τ_set ≈ 15 orbits (see Equation (56)). The result of this test is shown in Figure B1. If particles are initialized below the expected extent, which, for our choice of parameters, is z_max ≈ 0.25H after Equation (40), they are excited to the expected scale, whereas they settle to z_max if they are initialized above this scale. We note that the bottom left two panels in Figure B1 show more turbulent features than the bottom right panel, because here, the flow was initialized with a subcritical Richardson number, and the KHI had to become active sooner to stir particles up.

Zoom In Zoom Out Reset image size

We conclude that the initial particle scale height does not affect the final particle scale height, and ${H}_{{\rm{p}},0}$ can, in principle, be chosen arbitrarily. In practice, for simulations with St = 0.005, we initialize particles with a scale height close to the expected value in Equation (40) in order to minimize computational time. For St = 0.2, the settling time in Equation (56) is very short, and we can choose a higher initial height (also see Table 1).

In this paper, we present simulations with a resolution of 64 × 64 × 64 and 128 × 128 × 128 (see Table 1). In particular, the lower resolution is not sufficient to resolve the fastest growing streaming instability modes, which are always on the smallest scales. Still, due to being numerically inexpensive, it allowed us to conduct a comparatively extensive parameter study. Moreover, the diffusion-limited collapse criterion presented in Section 2 and by Klahr & Schreiber (2020) suggests that not the fastest growing SI wave-mode but the critical radius in Equation (21) needs to be resolved to accurately present gravitational collapse and planetesimal formation. The fact that the low-resolution fiducial run including self-gravity required a lower Q value for fragmentation to occur than the corresponding high-resolution simulation (see Figures 9 and 10) suggests that the low-resolution simulations were not in fact able to sufficiently resolve r_crit. Since we do not test our criterion for a resolution exceeding N_x,y,z = 128, we cannot be sure that our high-resolution simulations are converged. Still, our estimate for the critical radius of r_crit ≈ 0.006H covers two grid cells for N_x,y,z = 128. It therefore stands that those simulations are, if not fully converged, close to numerical convergence.

We note that due to the strong spatial and temporal fluctuations in particle concentration inherit to the nonlinear saturation of streaming instability, numerical convergence is challenging to test by solely evaluating fragmentation. For example, the evolution of the maximum particle density shown in Figure 9 before self-gravity is turned on is qualitatively similar for the two fiducial runs. Due these challenges and the large computational cost of performing a comprehensive convergence test, we leave this exercise for future work.