Streaming Instability in Turbulent Protoplanetary Disks

TL;DR

This paper re-examines the streaming instability in turbulent protoplanetary disks using an alpha turbulence model, revealing how turbulence intensity affects instability scales and viability within specific parameter ranges.

Contribution

It provides a detailed analysis of how turbulence impacts the growth and scale of streaming instabilities, identifying conditions under which they are viable in protoplanetary disks.

Findings

Turbulence reduces streaming instability growth rates compared to laminar models.

Moderate turbulence shifts instability scales toward larger, vertically oriented sheets.

Streaming instability is viable in a narrow alpha–Stokes number parameter space with significant growth rates.

Abstract

The streaming instability for solid particles in protoplanetary disks is re-examined assuming the familiar alpha ($\alpha$) model for isotropic turbulence. Turbulence always reduces the growth rates of the streaming instability relative to values calculated for globally laminar disks. While for small values of the turbulence parameter, $\alpha < 10^{-5}$, the wavelengths of the fastest-growing disturbances are small fractions of the local gas vertical scale height $H$, we find that for moderate values of the turbulence parameter, i.e., $\alpha \sim 10^{-5}-10^{-3}$, the lengthscales of maximally growing disturbances shift toward larger scales, approaching $H$. At these moderate turbulent intensities and for local particle to gas mass density ratios $\epsilon < 0.5$, the vertical scales of the most unstable modes begin to exceed the corresponding radial scales so that the instability appears in the form of vertically oriented sheets extending well beyond the particle scale height. We find that for hydrodynamical turbulent disk models reported in the literature, with $\alpha = 4\times 10^{-5} - 5\times 10^{-4}$, together with state of the art global evolution models of particle growth, the streaming instability is predicted to be viable within a narrow triangular patch of $\alpha$--$\tau_s$ parameter space centered on Stokes numbers, $\tau_s \sim 0.01$ and $\alpha \sim 4\times 10^{-5}$ and, further, exhibits growth rates on the order of several hundred to thousands of orbit times for disks with 1 percent ($Z= 0.01$) cosmic solids abundance or metallicity. Our results are consistent with, and place in context, published numerical studies of streaming instabilities.

Figures

b). During this initial bounce phase ($t<50$P*orb*) the number of filamentary structures is impossible to assess. As such, we examined the behavior and character of this simulation after the initial response phase passed (after $t=50$P*orb*), when the number of filaments was first countable. For the corresponding value of $H_{p}/H$ ($\approx 0.01$), the predicted $\lambda_{xm}$ better matches the corresponding measured $\lambda_{{}_{x,{{\rm obs}}}}$.

Tables

Table 1. Table 1 Critical values of α 𝛼 \alpha and τ s subscript 𝜏 𝑠 \tau_{s} as a function of Z 𝑍 Z for δ = 0.05 𝛿 0.05 \delta=0.05 .

$Z$	${\tau }_c$	${\alpha }_c$	${ϵ}_c$
0.001†	0.044	$4.4\times {10}^{-8}$	$\approx 1$
0.010	0.45	$3.7\times {10}^{-5}$	1.10
0.020	1.05	$2.0\times {10}^{-4}$	1.44
0.030	1.41	$4.3\times {10}^{-4}$	1.72
0.040	1.71	$7.1\times {10}^{-4}$	1.96
0.080	3.50	$2.5\times {10}^{-3}$	3.74

Table 2. Table 2 Simulations of Yang et al. ( 2017 ) compared to theoretical predictions

$\begin{array}{c}\mathrm{run}\\ \mathrm{type}\end{array}$	$\delta$	${\tau }_s$	$Z$	${\left(\frac{H_p}{H}\right)}^{†}$	$\frac{{\alpha }_{\mathrm{eff}}}{{10}^{-7}}$	${ϵ}_{\mathrm{eff}}$	$\frac{t_{{}_{gm}}}{P_{\mathrm{orb}}}$	${\frac{t_{{}_o}}{P_{\mathrm{orb}}}}^{‡}$	$\frac{{\lambda }_{{}_{xm}}}{H}$	$\frac{{\lambda }_{{}_{x,\mathrm{obs}}}}{H}$	Identifiers♯
2D	0.05	0.01	0.01	—	13.8a	$\sim {0.90}^b$	11000	—	1.2	—	${\left(0.2H\right)}^2,2560H^{-1}$
	0.05	0.01	0.02	0.012	14.5	1.66	166		0.13	$\sim$ 0.10-0.20	${\left(0.2H\right)}^2,2560H^{-1}$
	0.05	0.01	0.02	0.011	12.1	1.81	106		0.10	0.08	${\left(0.4H\right)}^2,2560H^{-1}$
	0.05	0.01	0.04	0.010	10.0	4.0	65		0.06	0.05	${\left(0.2H\right)}^2,1280H^{-1}$
	0.05	0.01	0.04	0.0035	1.27	11.2	16		0.01	n/a c	${\left(0.4H\right)}^2,1280H^{-1}$
				0.010d	10.0	4.0	65		0.06	0.04
	0.05	0.001	0.03	—	5.4a	$\sim {1.3}^b$	39000	—	10	—	${\left(0.2H\right)}^2,1280H^{-1}$
	0.05	0.001	0.04	0.016	2.25	2.66	1850		0.14	0.10-0.20	${\left(0.2H\right)}^2,1280H^{-1}$
	0.05	0.001	0.04	0.015	2.56	2.5	2300		0.17	0.10-0.20	${\left(0.4H\right)}^2,1280H^{-1}$
3D	0.05	0.01	0.02	0.010	10.0	2.0	71		0.08	0.067e	${\left(0.2H\right)}^3,640H^{-1}$
	0.05	0.001	0.04	0.011	2.0	3.6	800		0.07	$\sim {0.067}^f$	${\left(0.2H\right)}^3,640H^{-1}$

Table 3. Table 3 Selected simulations of Li et al. ( 2018 ) compared to theoretical predictions. Periodic boundary condition runs shown. Table heading are defined as in Table LABEL:Yang_etal_2017_simulations .

$\begin{array}{c}\mathrm{run}\\ \mathrm{type}\end{array}$	$\delta$	${\tau }_s$	$Z$	$\left(\frac{H_p}{H}\right)$	$\frac{{\alpha }_{\mathrm{eff}}}{{10}^{-7}}$	${ϵ}_{\mathrm{eff}}$	$\frac{t_{{}_{gm}}}{P_{\mathrm{orb}}}$	${\frac{t_{{}_o}}{P_{\mathrm{orb}}}}^{\oplus }$	$\frac{{\lambda }_{{}_{xm}}}{H}$	${\frac{{\lambda }_{{}_{x,\mathrm{obs}}}}{H}}^{‡}$	Identifiers
3D	0.05	0.314	0.02	0.001	3.14	20.0	0.5	$\sim 5$	0.007	$\sim 0.015$	${\left(0.2H\right)}^3,320H^{-1}$
	0.05	0.314	0.02	0.001	3.14	20.0	0.5	$\sim 5$	0.007	$\sim 0.016$	${\left(0.4H\right)}^3,320H^{-1}$
	0.05	0.314	0.02	0.001	3.14	20.0	0.5	$\sim 5$	0.007	$\sim 0.016$	${\left(0.8H\right)}^3,320H^{-1}$

Table 4. Table 4 Selected simulations of Gerbig et al. ( 2020 ) compared to theoretical predictions. Table heading are defined as in Table LABEL:Yang_etal_2017_simulations .

$\begin{array}{c}\delta \end{array}$	${\tau }_s$	$Z$	${\left(\frac{H_p}{H}\right)}^{▷}$	$\frac{{\alpha }_{\mathrm{eff}}}{{10}^{-7}}$	${ϵ}_{\mathrm{eff}}$	$\frac{t_{{}_{gm}}}{P_{\mathrm{orb}}}$	${\frac{t_{{}_o}}{P_{\mathrm{orb}}}}^{\oplus }$	$\frac{{\lambda }_{{}_{xm}}}{H}$	${\frac{{\lambda }_{{}_{x,\mathrm{obs}}}}{H}}^{‡}$	Identifiers
0.05	0.2	0.0002	0.0020	8	0.10	7.3		0.13	0.13	${\left(0.4H\right)}^3,160H^{-1}$
			$\approx 0^{♭}$	${\infty }^{\phi }$	—	—		—
			0.0055♯	60.5	0.036	17.4		0.45
0.05	0.2	0.01	0.0155	481	0.65	68.0	—	1.1	—	${\left(0.4H\right)}^3,320H^{-1}$
			0.0140♭	392	0.71	90.0		1.1
			0.0170♯	578	0.59	57.6		1.1
0.05	0.2	0.02	0.0125	313	1.60	13.9		0.18	0.13ℵ, 0.10γ	${\left(0.4H\right)}^3,320H^{-1}$
			0.0100♭	200	2.0	4.72		0.09
			0.0150♯	450	1.34	64.6		0.44
0.05	0.2	0.03	0.0085	144	3.55	2.77		0.06	$0.10$	${\left(0.4H\right)}^3,320H^{-1}$
			0.0070♭	98	4.28	2.07		0.05
			0.0100♯	200	3.0	3.84		0.08
0.05	0.2	0.20	0.0055	60	36.4	29.1		0.04	$0.065$	${\left(0.4H\right)}^3,160H^{-1}$
			0.0030♭	18	66.6	16.5		0.01
			0.0080♯	128	25.0	38.2		0.09
0.025	0.2	0.02	0.0030ϖ	18	6.66	2.16		0.018	$0.021$	${\left(0.4H\right)}^3,320H^{-1}$
0.10	0.2	0.02	0.0280	1568	0.71	90.9	—	2.15	—	${\left(0.4H\right)}^3,160H^{-1}$
			0.0245♭	1201	0.82	165.0		2.31
			0.0315♯	1984	0.64	60.45		66.5

Table 5. Table 5 Selected simulations of Yang et al. ( 2018 ) compared to theoretical predictions. Only ideal MHD simulations analyzed (see text). Table heading are defined as in Table LABEL:Yang_etal_2017_simulations .

$\delta$	${\tau }_s$	$Z$	${\left(\frac{H_p}{H}\right)}^{♯}$	$\frac{{\alpha }_{\mathrm{eff}}}{{10}^{-3}}$	${ϵ}_{\mathrm{eff}}$	$\frac{t_{{}_{gm}}}{P_{\mathrm{orb}}}$	$\frac{t_{{}_o}}{P_{\mathrm{orb}}}$	$\frac{{\lambda }_{{}_{xm}}}{H}$	${\frac{{\lambda }_{{}_{x,\mathrm{obs}}}}{H}}^{♭}$	Identifiers†
0.05	0.10	0.01	0.30a	9.0	0.07	787	—	38.0	—	$4H\times 8H\times 8H,16H^{-1}$
0.05	0.10	0.02	0.26b	6.8	0.08	642	—	30.2	—	$4H\times 8H\times 8H,16H^{-1}$
0.05	0.10	0.04	0.26b	6.8	0.16	870	—	31.7	—	$4H\times 8H\times 8H,16H^{-1}$
0.05	0.10	0.08	0.20a	4.1	0.41	9640	—	32.6	—	$4H\times 8H\times 8H,16H^{-1}$

Table 6. Table 6 Maximum ∗ achievable Stokes numbers ( τ s subscript 𝜏 𝑠 \tau_{s} ) in global disk evolution models for Z = 0.01 𝑍 0.01 Z=0.01

	$\begin{array}{c}\mathrm{Inner}\mathrm{Disk}\\ {\left(R=3AU\right)}^{△}\end{array}$		$\begin{array}{c}\mathrm{SnowLine}\\ {\left(R=5AU\right)}^{△}\end{array}$		$\begin{array}{c}\mathrm{Outer}\mathrm{Disk}\\ {\left(R=30AU\right)}^{△}\end{array}$
$\frac{\alpha }{{10}^{-4}}$	——		——		——
	${\tau }_s$	$a$(cm)∗	${\tau }_s$	$a$(cm)∗	${\tau }_s$	$a$(cm)∗
100†	$7.0\times {10}^{-4}$	0.03	$9.3\times {10}^{-3}$	0.43	$6.9\times {10}^{-3}$	0.07
40††	$1.8\times {10}^{-3}$	0.38	$1.1\times {10}^{-2}$	1.6	$6.2\times {10}^{-3}$	0.09
10†	$2.5\times {10}^{-3}$	0.76	$1.9\times {10}^{-2}$	4.6	$4.8\times {10}^{-3}$	0.11
4††	$3.4\times {10}^{-3}$	1.9	$1.3\times {10}^{-2}$	5.3	$7.7\times {10}^{-3}$	0.12
1†	$5.7\times {10}^{-3}$	5.8⋄	$1.2\times {10}^{-2}$	13	$1.1\times {10}^{-2}$	0.30
0.1†	$2.2\times {10}^{-2}$	12⋄	$4.9\times {10}^{-2}$	39⋄	$1.6\times {10}^{-2}$	0.37
10‡	$3.1\times {10}^{-3}$	0.79	$4.5\times {10}^{-3}$	0.89	$1.9\times {10}^{-3}$	0.05

Equations

$$\frac{L_{{}_{\mathrm{SI}}}}{r}=\frac{\partial P/\partial r}{2{\rho }_{{}_g}{\Omega }^{2}r}\sim {\delta }^2.$$

$$\displaystyle t_{s}=\left\{\begin{array}[]{lr}\displaystyle\frac{\rho_{{}_{\bullet}}a}{\rho_{g}c_{s}};&a\leq(4/9)\lambda_{mfp};\\ \displaystyle\frac{8\rho_{{}_{\bullet}}a}{3\rho_{g}C_{d}\Delta V_{pg}};&a>(4/9)\lambda_{mfp}\end{array}\right.,$$

$$\frac{\alpha c_{s}H}{1+{\tau }_s^2}{\nabla }^2{\rho }_p,$$

$$-c_d^2\nabla {\rho }_p;c_d^2=\frac{\alpha c_s^2}{1+{\tau }_s^2},$$

$${\rho }_g\left({\partial }_t+u_g{\partial }_x+w_g{\partial }_z\right)u_g-2\Omega {\rho }_g{v}_g$$

$$\displaystyle\rho_{g}\left(\partial_{t}+u_{g}\partial_{x}+w_{g}\partial_{z}\right)v_{g}+\small\mbox{$\frac{1}{2}$}\Omega\rho_{g}u_{g}$$

$${\rho }_g\left({\partial }_t+u_g{\partial }_x+w_g{\partial }_z\right)w_g$$

$${\partial }_x{u}_g+{\partial }_z{w}_g$$

$$-\partial \tilde{P}/\partial r=2\eta {\Omega }^{2}r=2\delta c_s\Omega .$$

$${\rho }_p\left({\partial }_t+u_p{\partial }_x+w_p{\partial }_z\right)u_p-2{\rho }_p\Omega v_p=-c_d^2{\partial }_x{\rho }_p+{\rho }_p{\rho }_g\mu (u_g-u_p),$$

$$\displaystyle\rho_{p}\left(\partial_{t}+u_{p}\partial_{x}+w_{p}\partial_{z}\right)v_{p}+\small\mbox{$\frac{1}{2}$}\rho_{p}\Omega u_{p}=\rho_{p}\rho_{g}\mu(v_{g}-v_{p}),$$

$${\rho }_p\left({\partial }_t+u_p{\partial }_x+w_p{\partial }_z\right)w_p=-c_d^2{\partial }_z{\rho }_p+{\rho }_p{\rho }_g\mu (w_g-w_p),$$

$${\partial }_t{\rho }_p+{\partial }_x{\rho }_p{u}_p+{\partial }_z{\rho }_p{w}_p=\frac{\alpha c_{s}H}{1+{\tau }_s^2}{\nabla }^2{\rho }_p,$$

$$U_{{}_{g0}}$$

$$V_{{}_{g0}}$$

$$U_{{}_{p0}}$$

$$V_{{}_{p0}}$$

$$\Delta U\equiv U_{{}_{g0}}-U_{{}_{p0}}={\tau }_s\left(\frac{2+2ϵ-3{\tau }_s\alpha /2}{(1+ϵ{)}^2+{\tau }_s^2}\right)\delta c_{{}_s}.$$

$$u_g\to U_{{}_{g0}}+u_g^{\prime },v_g\to V_{{}_{g0}}+v_g^{\prime },w_g\to w_g^{\prime },p\to p^{\prime },$$

$$u_p\to U_{{}_{p0}}+u_p^{\prime },v_p\to V_{{}_{p0}}+v_p^{\prime },w_p\to w_p^{\prime },{\rho }_p\to {\overline{\rho }}_p\left(1+{\Delta }_d^{\prime }\right),$$

$$u_g^{\prime }={\partial }_z{\psi }^{\prime },w_g^{\prime }=-{\partial }_x{\psi }^{\prime }.$$

$${\omega }_g^{\prime }\equiv {\partial }_z{u}_g^{\prime }-{\partial }_x{w}_g^{\prime },{\omega }_p^{\prime }\equiv {\partial }_z{u}_p^{\prime }-{\partial }_x{w}_p^{\prime }.$$

$$\psi^{\prime}\rightarrow\hat{\psi}\exp\big{(}ik_{{}_{x}}x+ik_{{}_{z}}z-i\omega t\big{)}+{\rm c.c.},$$

$$F(\omega ;k_x,k_z,{\tau }_s,ϵ,\delta ,\alpha )=0,$$

$$ϵ=ϵ(\alpha ,{\tau }_s,Z)=Z\sqrt{\frac{\alpha +{\tau }_s}{\alpha }},$$

$${\omega }^2=\frac{{\Omega }^2{k}_z^2}{k_x^2+k_z^2},$$

$$\omega =k_x\left(U_{{}_{p0}}-U_{{}_{g0}}\right)$$

$$k_z^2=\frac{k_x^4}{k_a^2-k_x^2},k_a\equiv \frac{(1+ϵ{)}^2+{\tau }_s^2}{2\delta (1+ϵ){\tau }_{s}H}.$$

$$t_{{}_d}^{-1}\equiv \frac{\left|U_{{}_{p0}}\right|}{r}=\frac{2{\tau }_s{\delta }^2\Omega }{(1+ϵ{)}^2+{\tau }_s^2}.$$

$$\alpha\tau_{s}>{U_{f}^{2}}\Big{/}{2c_{s}^{2}}.$$

$$-i\omega\left(\begin{array}[]{c}\psi^{\prime}\\ v_{g}^{\prime}\\ u_{p}^{\prime}\\ v_{p}^{\prime}\\ w_{p}^{\prime}\\ \Delta_{p}^{\prime}\end{array}\right)+\mathtt{M}\left(\begin{array}[]{c}\psi^{\prime}\\ v_{g}^{\prime}\\ u_{p}^{\prime}\\ v_{p}^{\prime}\\ w_{p}^{\prime}\\ \Delta_{p}^{\prime}\end{array}\right)=0,$$

$$\mathtt{M}\equiv\left(\begin{array}[]{cccccc}a&\frac{2ik_{z}}{K^{2}}&\frac{ik_{z}\epsilon}{K^{2}\tau_{s}}&0&-\frac{ik_{x}\epsilon}{K^{2}\tau_{s}}&-\frac{2ik_{z}\epsilon(\epsilon+1)}{K^{2}\left((\epsilon+1)^{2}+\tau_{s}^{2}\right)}\\ \frac{ik_{z}}{2}&a&0&-\frac{\epsilon}{\tau_{s}}&0&-\frac{\epsilon\tau_{s}}{(\epsilon+1)^{2}+\tau_{s}^{2}}\\ -\frac{ik_{z}}{\tau_{s}}&0&\frac{1}{\tau_{s}}-\frac{2ik_{x}\tau_{s}}{(\epsilon+1)^{2}+\tau_{s}^{2}}&-2&0&ik_{z}c_{d}^{2}\\ 0&-\frac{1}{\tau_{s}}&\frac{1}{2}&\frac{1}{\tau_{s}}-\frac{2ik_{x}\tau_{s}}{(\epsilon+1)^{2}+\tau_{s}^{2}}&0&0\\ \frac{ik_{x}}{\tau_{s}}&0&0&0&\frac{1}{\tau_{s}}-\frac{2ik_{x}\tau_{s}}{(\epsilon+1)^{2}+\tau_{s}^{2}}&ik_{z}c_{d}^{2}\\ 0&0&ik_{x}&0&ik_{z}&\frac{K^{2}\alpha/\delta^{2}}{\tau_{s}^{2}+1}-\frac{2ik_{x}\tau_{s}}{(\epsilon+1)^{2}+\tau_{s}^{2}}\\ \end{array}\right),$$

$$a\equiv \alpha K^2/{\delta }^2+\frac{ϵ}{{\tau }_s}+\frac{2i{k}_xϵ{\tau }_s}{(ϵ+1{)}^2+{\tau }_s^2}.$$

$$\det \left(-i\omega \mathtt{I}+\mathtt{M}\right)=0,$$

Full text

1. Introduction

The formation of the first 100 km size planetesimals remains a poorly understood chapter in the standard story of planetary formation. A given radial zone in a protoplanetary disk will initially contain $\mu$m size particles of material which is solid under local conditions. These grains grow by sticking until they become sub-mm to mm-sized aggregates (Birnstiel et al., 2012; Estrada et al., 2016). In the inner nebula, some (still mysterious) heating process melts these (ice-poor) aggregates and forms chondrules. These chondrules apparently can form few-cm-scale aggregates depending on local disk properties (Simon et al., 2018). The picture is less clear in the outer nebula, where chondrule formation may never occur and the material remains ice-rich.

However, what remains elusive is understanding how these aggregate mm-to-cm-size clusters eventually coalesce into 100 km sized planetesimals, especially if the nebula is even weakly turbulent. This is because further particle evolution by “incremental growth” (by sticking) encounters several barriers - including (but not limited to) the radial-drift barrier, the bouncing barrier, and the fragmentation barrier. For a more comprehensive discussion of these barrier mechanisms see the discussion found in Brauer et al. (2008); Zsom et al. (2010); Birnstiel et al. (2012), and Estrada et al. (2016). Additional barriers to incremental growth in turbulence reappear at 1-10km size, due to gravitational stirring of such objects by fluctuations in the gas density, much like giant molecular clouds scatter stars in our galaxy (Ida et al., 2008; Yang et al., 2009; Nelson, & Gressel, 2010; Gressel et al., 2011; Yang et al., 2012). This latter realization has led to a growing suspicion that planetesimals may have been “born big”, close to the typical 100km sizes we see today (Johansen et al., 2007; Cuzzi et al., 2008; Morbidelli et al., 2009). Moreover, recent meteroritical work suggests that a substantial amount of the accretion that formed the solar nebula’s first planetesimals, some of them leaving behind only their molten Fe-Ni cores, and even the initial 20-50 M⊕ core of Jupiter, occurred in less than 0.5 Ma after the formation of the first remaining solids (the so-called Calcium-Aluminum refractory oxide inclusions, Kruijer et al., 2017). Thus, formation of sizeable planetesimals seemed to have started early, well inside the snowline, and after Jupiter’s initial core formed (which probably required snowline planetesimals as precursors). It is in this context that current planetesimal formation theories must be assessed. How can planetesimals be born big, starting very early (and continuing for several Myr), directly from mm-cm size objects?

One popular hypothesis is that clumps of small particles are collected into 100km or larger planetesimals by the streaming instability (SI) and ultimately gravitational collapse. SI is a linear instability (grows without limit from small perturbations, under the right conditions) that can enhance the concentration of particles in protoplanetary disks (Youdin & Goodman, 2005; Youdin & Johansen, 2007) (YG2005 and YJ2007 hereafter). The dynamics involves the resonant momentum exchange between the disk gas and its embedded particles treated as a pressure-free second fluid (Squire & Hopkins, 2018a, b; Hopkins & Squire, 2018a, b) – see also Section 2.1. In protoplanetary disks, the SI is strongest for axisymmetric disturbances and the growth rates are most rapid when the local volume mass densities in the gas ($\rho_{g}$) and particle components ($\rho_{p}$) are comparable, i.e., when $\epsilon\equiv\rho_{p}/\rho_{g}\geq 1$. Linear stability analyses indicate that for laminar Keplerian flows, the SI grows fastest for small wavelengths and for near-unity Stokes numbers $\tau_{s}\equiv\Omega t_{s}$, where $t_{s}$ is the particle gas drag stopping time and $\Omega$ is the local disk rotation time (also see section 2.2)111In fact, the inviscid calculation indicates that the growth rates asymptote to finite values as the wavelength of the vertical disturbance approaches zero provided the radial wavelength is larger than some minimum value. This suggests the problem may be ill-posed in the inviscid limit. However this short wavelength catastrophe is averted when viscosity is included (see our general results below). Note, in our discussion throughout this work we sometimes refer to $\tau_{s}$ by its more commonly used symbolic designator, “St”.

These features of the SI suggest that this process may play an important role either (a) in the late stage of a protoplanetary disk’s evolution when the disk has lost most of its gas because of photodissociation-or-MHD-driven winds from the star or disk (e.g., Ercolano et al., 2017; Carrera et al., 2017) , and/or accretion onto the central object or (b) if the disk is laminar (nonturbulent), in which case the particle component can settle to high $\epsilon$ near the disk midplane. Several detailed numerical simulations of laminar disks have shown that the SI rapidly enhances local particle concentrations (e.g., Bai & Stone, 2010; Lyra & Kuchner, 2013; Carrera et al., 2015; Yang & Johansen, 2014; Yang et al., 2017; Schreiber & Klahr, 2018; Li et al., 2018, among many others). Particle concentration is further helped along if particle self-gravity is included in models – e.g., as done in Johansen et al. (2007); Simon et al. (2016, 2017) and most recently Gerbig et al. (2020) – and can drive particle enhancements to the precipice of gravitational collapse and onwards (Johansen et al., 2015; Schäfer et al., 2017; Li et al., 2019).

However, regions of protoplanetary disks in which particle growth is of greatest interest (i.e., 1-100 AU) are possibly weakly-to-moderately turbulent (Turner et al., 2014; Lyra, & Umurhan, 2019). That is, recent theoretical advances suggest that non-ionized regions of protoplanetary disks support several instability processes that can lead to sustained turbulence: the Vertical Shear Instability (VSI), Convective Overstability (COV) and Zombie Vortex Instability (ZVI).

Turbulence in disks is often thought of in terms of a zero-order closure “alpha-disk” model, wherein gas turbulence is represented by an enhanced viscosity quantified by a turbulent (kinematic) viscosity coefficient, $\nu_{t}\equiv\alpha c_{s}H$,where $c_{s}$ and $H$ are the local isothermal soundspeed and the vertical pressure scale heights, respectively (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974).222A turbulent viscosity assumes that downscale momentum transfer occurs in the inertial range of a fully developed statistically steady turbulent fluid. The “$\alpha$-model” scales this in terms of the typical speed- and length-scales encountered in locally rotating sections of protoplanetary disks ($c_{s}$ and $H$). As such, the quantity $\alpha$ is typically interpreted to be the inverse of the local turbulent Reynolds number, Ret and might be thought of as a measure of the amplitude of the turbulent velocity field compared to the local sound speed. For a more pedagogical review appropriate to astrophysical fluids see Regev et al. (2016). The alpha-disk model, notwithstanding its crudeness, does a good job in characterizing disk structure and consequent flow in a turbulent protoplanetary disk medium. Most recently Stoll et al. (2017) examined the turbulent state of the VSI and found the emergence of large-scale meridional flow to be well predicted by an $\alpha$-model albeit with effectively non-isotropic diffusion stresses owing to its characteristic large vertical motions. However, MHD turbulence might not be as well represented by an $\alpha$-model, and we keep this in mind throughout. Numerical analyses of the three above-mentioned turbulence generating mechanisms report values of $\alpha$ in the range of $10^{-5}-10^{-3}$. This turbulence may also, in principle, diffuse particles away from the disk midplane, reducing the values of the density ratio ($\epsilon$) near the midplane, while also radially diffusing and dispersing radial concentrations of particles before they can grow appreciably (Fromang & Papaloizou, 2006; Okuzumi & Hirose, 2011; Zhu et al., 2015; Riols & Lesur, 2018; Yang et al., 2018) . On the other hand, some numerical simulations seem to show SI manifesting even in the presence of turbulence (see below).

So, what really happens to the SI in the presence of turbulence? YG2005 presented a brief and mostly qualitative discussion of the possibilities, but declined to pursue them on the grounds that protoplanetary disks were probably nonturbulent. There are a few numerical simulations examining the fate of the SI in either an axisymmetric or fully 3D model of a protoplanetary disk experiencing some level of turbulent motions, whether self-excited or driven by some other mechanism (including but not limited to, Johansen et al., 2007; Balsara et al., 2009; Bai & Stone, 2010; Tilley et al., 2010; Carrera et al., 2015; Yang et al., 2017, 2018; Li et al., 2018).333 Note that Fromang & Papaloizou (2006)’s two-fluid turbulent set-up is a possible “precursor” turbulent SI analysis, although this is not yet been demonstrated. Some of these studies examined the fate of particle clumping in the presence of magnetorotational turbulence, with turbulent intensity $\alpha\sim 10^{-3}-10^{-2}$, and for a range of the two particle parameters $\tau_{s}$ and $\epsilon$ (Fromang & Papaloizou, 2006; Johansen et al., 2007; Balsara et al., 2009; Tilley et al., 2010), while others modeled only the self-generated “midplane” turbulence surrounding a layer of particles that had settled toward the midplane of an otherwise laminar disk flow (e.g., Bai & Stone, 2010; Carrera et al., 2015; Yang et al., 2017; Li et al., 2018), and most recently in a forced driven model of turbulence (Gole et al., 2020)444This study appeared during the revision phase of this work..

While the final state of the SI subject to these different kinds of turbulence indeed varies, a common point of agreement between all of these investigations is that interesting solids clumping and instability emerges for parameter values in which $\epsilon\sim 1$ and, most importantly, when $0.1<\tau_{s}<3$. The results reported in Johansen et al. (2007) are particularly noteworthy as they show that the SI emerges, in turbulence, with a preferred radial lengthscale of about one pressure scale height ($H$) and has a growth time scale of about 10 local orbital periods. A question facing these and other previous numerical studies of SI in turbulence is whether such combinations of initial conditions - large particles that have somehow grown without being disrupted in such moderate levels of turbulence - are self-consistent (see section 7.7). Meanwhile, several other studies have shown that quite small particles can undergo SI in disks that are not turbulent at all, globally; the particles experience only a tiny amount of local turbulence, called “midplane turbulence”, generated by the densely settled particle layer (Barranco , 2009; Carrera et al., 2015; Yang et al., 2017). Our results explain these different outcomes in a unified and consistent way.

Well-resolved numerical experiments of two-fluid processes are expensive. A theoretical prediction for the fate of the SI under arbitrary turbulent protoplanetary disk conditions would be a useful tool both in developing some quantitative estimate for the expected length and timescales of growth in such a nebula, and in planning future detailed numerical experiments. We present such a theory, extending the SI analysis done in YG2005 and YJ2007 with the addition of a simple $\alpha$-model of disk turbulence, that provides an effective turbulent viscosity acting on the gas as well as a prescription for the effective turbulent particle diffusion resulting from the statistically steady stirring of particles by inertial scale gaseous eddies. The basic assumption is that the fundamental processes driving turbulence are unaffected by the particles, and lead to a statistically steady isotropic turbulence in the gas555This tactic was hinted at in Section 3.2.2. of YG2005. Also, an assumption similar to this was employed in studying Tolmien-Schlichting waves of cold (non-MHD) turbulent protoplanetary disks (Umurhan, & Shaviv, 2009). (see however Lin, 2019, and Section 7.6).

We present a linear stability analysis of the SI in such an $\alpha$-disk model of protoplanetary disk turbulence. We determine the growth rate and wavelength of the fastest growing modes as functions of various properties including gas disk opening angle, particle Stokes numbers, local particle-to-gas mass density ratio, and the intensity of turbulence. r During the revision phase of this manuscript, Chen & Lin (2020) released a study with similar aims and we find that our results are in mutual agreement. In Section 2 we review the basic properties of the SI, including recent theoretical developments regarding resonant drag instabilities. In Section 3 we motivate our model equations, describe the steady state, introduce infinitesimal perturbations, and discuss various caveats with respect to our representation of turbulence. Section 4 is concerned with verification. In Section 5 we survey the results of the stability analysis, focussing on the most rapidly growing modes. Section 6 applies our theory to four recently published numerical studies of the SI. In Section 7 we identify regions in the parameter space of particle Stokes numbers and disk turbulent intensity for which the SI remains a feasible path to planetesimal formation at various locations of a model disk with 1 percent metallicity. We consider these in light of various known barriers to particle growth. In Section 8 we summarize our findings, we discuss various issues spurred by our analysis, and point to future directions.

2. SI Mode review

2.1. Broad physical picture

Because a pressure-supported gaseous disk orbits the central object at sub-Keplerian speeds, momentum exchange between gas and particles via drag forces induces a relative radial drift. In steady state, for example, a radially diminishing steady pressure profile typically causes the gas to spiral out while causing a single-size particle to spiral inward (e.g., see Eqns. [13-14]). If multiple particle size species are considered, one or several (but not all) of their smaller components can spiral outwards with the gas (e.g., Estrada et al., 2016; Benítez-Llambay et al., 2019).

The SI arises from perturbations in this relative drifting steady state and how it modifies oscillatory motions in the disk: Momentum exchange is generally modeled as a function of the relative velocities between the two fluids multiplied by the product of the two fluid densities times a drag coefficient representing the type of physically appropriate drag mechanism. The momentum channeled from the mean state and into perturbations through modifications of the drag exchange term due to particle density fluctuations is the root of the linear instability. YG2005 and YJ2007 show that these density fluctuations draw momentum from the mean drift state and destabilize oscillating disk inertial waves (e.g., see YJ2007). The insightful precursor toy model of Goodman & Pindor (2000) argues that this sort of mean-momentum wave-phase sensitive “tapping” via gas drag can generically lead to instability in otherwise damped oscillating systems.

A recent comprehensive theoretical study (Squire & Hopkins, 2018a, b) demonstrates that the SI is a member of a particular class of resonant drag instabilities (RDI). A two-fluid system, in which one component is pressure free and streaming with velocity ${\bf w}_{s}$ with respect to the second (non-zero pressure) fluid, is potentially resonantly unstable to any wave phenomenon with wavevector ${\bf k}$ supported by the fluid if ${\bf w}_{s}\cdot{\bf k}$ equals the oscillation frequency of the wave phenomenon. In this broad framework the SI is an RDI arising from the particle stream’s resonance with the inertial waves supported by the gas. We apply this prescription in rationalizing the trends contained in the inviscid models discussed in the verification section 4.

Generically speaking, however, the potential for instability holds for any class of waves that the fluid system might support including sound waves, gravity waves, magnetosonic waves as well as Rossby waves and potentially many others (Hopkins & Squire, 2018a, b). For example, Schreiber & Klahr (2018) recently have shown that the SI occurs for non-axisymmetric vertically restricted disturbances in simulations of disks which means that the waves with which the particle stream becomes resonant are not axisymmetric inertial oscillations but, instead, either non-axisymmetric inertial oscillations or Rossby waves. Similar effects seem to be characterizing the particle-vortex numerical experiments recently reported in Surville & Mayer (2018). Three key ingredients for resonance are that (a) there exists some means of momentum/energy exchange between the two fluid systems, for example, whether it be by means of classical fluid drag (as it is for the SI), or via dynamical drag if the two fluids are self-gravitating (e.g., see Chapter 13 of Chandrasekhar, 1961), (b) a relative drift velocity between the particle and gas components manifests, and (c) the fluid component supports some kind of wave phenomenon.

2.2. Some physical properties

We review some of the basic physical properties of the SI based on the analysis of YJ2007. The analysis here and throughout this paper is based on a (nearly) point-analysis performed at some disk cylindrical radial position $r$. The disk position is assumed to be locally isothermal666 “Local” in the sense common in protoplanetary disk literature, namely, that it is only a function of the radial coordinate. characterized by a temperature $T$ and soundspeed $c_{s}\equiv\sqrt{{\cal R}T/\mu}$, where the gas constant ${\cal R}\approx 8314$ J/kg/K and $\mu$ is the gas mean molecular weight. The local rotation rate of the disk is $\Omega$, the Keplerian rotation speed is $v_{{}_{K}}=r\Omega$, and the effective thickness of the disk is measured by the disk-opening angle parameter $\delta\equiv c_{{}_{s}}/v_{{}_{K}}=H/r$.

We assume there exists a global pressure field ($P$) whose gradient varies on the radial disk scale, i.e., \partial P/\partial r={\cal O}\left({P\big{/}r}\right). The SI operates on length scales $L_{{}_{{\rm SI}}}$ given dimensionally (not quantitatively) by

[TABLE]

This means that $L_{{}_{{\rm SI}}}\sim\delta H$; In this work, we adopt the definition $L_{{}_{{\rm SI}}}\equiv\delta H$. We refer to $\delta$ as measuring the relative flow of the particles past the gas at the midplane typically $\sim\delta c_{s}$ (see below).777The lengthscale $\delta H$ is typical also of the most unstable VSI modes (Umurhan et al., 2016), and the quantity $\delta^{2}$ is the same as the pressure gradient parameter $\eta$ of Nakagawa et al (1986). The disk opening angle is also equal to $\delta$.

The analysis of YJ2007 (and YG2005) assumes that vertical variation of gas density plays no role and there is no momentum or mass diffusion (ie., no turbulent viscosity or diffusivity). They show that the SI is the primary instability of axisymmetric inertial modes. The stability of a given mode with radial and vertical wavenumbers, $k_{{}_{x}}$ and $k_{{}_{z}}$ respectively, is a function of the local value of $\epsilon$ and the stopping time $t_{s}$ defined according to whether the particles are in the Epstein or Stokes regimes respectively:

[TABLE]

where $\rho_{{}_{\bullet}}$ is the density of a given particle, $a$ is the particle radius, and $\lambda_{mfp}$ is the molecular mean-free-path. The Epstein regime is appropriate for particles whose radii $a\leq(4/9)\lambda_{mfp}$ while the stopping times for $a>(4/9)\lambda_{mfp}$ are for the Stokes regime, in which $C_{d}$ is a particle drag coefficient and $\Delta V_{pg}$ is the relative speed between a particle and the gas (Weidenschilling 1977). Depending on the size of the particle and $\Delta V_{pg}$, $C_{d}$ may itself be a function of $\Delta V_{pg}$ (e.g., see Estrada et al., 2016). As noted previously, the stopping times are scaled by the local orbital time, giving the “Stokes number” $\tau_{s}\equiv\Omega t_{s}$.

For a given pair of input parameters ($\tau_{s},\epsilon$), instabilities are typically expected for values of $k_{x}<k_{z}$ and growth rates are found to be maximal for values of $L_{{}_{{\rm SI}}}k_{x}\approx L_{{}_{{\rm SI}}}k_{z}\sim 100$ or more (see Figure 1 of YJ2007 and Figure 2 of YG2005). That is, the wavelengths of fastest growth are usually much smaller than $L_{\rm SI}=\delta H$. Instability is most favorable for values of $\epsilon\sim 1$. Disturbingly, the problem as set up appears somewhat ill-posed, in that instability appears to persist for certain finite values of $k_{x}$ as $k_{z}\rightarrow\infty$. In cases of most physical interest, the instability growth timescales must be much faster than the radial drift rates (see Figure 8 of YG2005). The analysis of YJ2007 assumes the gaseous component is compressible, yet they demonstrate that the contribution of gas compressibility is negligible to the instability mechanism (see also Section 4). As such, they show that the particle fluctuations are the key ingredient of the instability.

3. Model equations

3.1. Assumptions

We build upon the inviscid model setup of YG2005. We review the assumptions made and indicate what is new to this paper:

The gas component is incompressible. As examined in the original studies, the growth rate of the strongest inviscid SI is on the orbit timescale, while its typical lengthscale $\sim\delta H$. This means that acoustic disturbances propagate across those relevant lengthscales on correspondingly shorter timescales $\sim\delta P_{orb}$, where $P_{orb}\equiv 2\pi/\Omega$ is the local disk orbit time. 2. 2.

Spatial variations in all disk quantities are negligible except for the mean Keplerian shear, which is assumed constant. 3. 3.

There are no disk vertical density variations (no vertical density stratification). 4. 4.

The Stokes numbers are constant. 5. 5.

We consider only particles consisting of a single mass-dominant size and, thus, ours is a “2-fluid” model. Particle growth evolutionary models suggest that the assumption of a mass-dominant size is not bad (Zsom et al., 2010; Estrada et al., 2016). Because Krapp et al. (2019) have recently shown that the streaming instability is diminished for disk models containing most reasonable particle distributions in even weakly turbulent disks 888Krapp et al. (2019) show that only for particle distributions with a very restricted range in particle sizes and mass loading, i.e., $0.0001<{\rm St}<0.01$ with $\epsilon>0.4$, is the growth rates of the SI enhanced compared to the two-fluid approach typically taken. See their Figure 5. , our assumption is conservatively favorable to SI. 6. 6.

Scales of interest are small enough so that it is appropriate to use the shearing box assumption, which neglects large scale effects including curvature. 7. 7.

Gas and particle perturbations are axisymmetric.. 8. 8.

(New to this paper) Turbulence is assumed isotropic and is represented in the gas momentum equation by the standard $\alpha$-disk model in which turbulence is modeled as an enhanced kinematic viscosity, $\nu_{t}=\alpha c_{s}H$, controlled by the non-dimensional parameter $\alpha$. In other words, internal stresses arising from momentum exchange due to turbulence is modeled with the term, $\nu_{T}\rho_{g}\nabla^{2}{\bf u}_{g}$, on the righthand side of the gas momentum equation (e.g. Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974). We also include the effect of radial accretion of the gas component due to the underlying turbulence-driven viscous evolution. 9. 9.

(New to this paper) Turbulence causes stirring of the particle component, represented by a turbulent diffusion term as a source term in the particle mass conservation equation (Cuzzi et al., 1993; Dobrovolskis et al., 1999; Youdin & Lithwick, 2007; Carballido et al., 2011; Estrada et al., 2016):

[TABLE]

which captures the effect of diminished stirring of particles with large inertias ($\tau_{s}\gg 1$). 10. 10.

(New to this paper) Since collisions between particles are not important the particle phase is typically assumed to be pressure free. However, turbulent stirring introduces an effective pressure gradient upon the particle momentum conservation in the form of an effective particle pressure term,

[TABLE]

(Cuzzi et al., 1993; Dobrovolskis et al., 1999; Jacquet et al., 2011).The effects of this “particle pressure” term are small, except for the wavelengths of the fastest growing mode.

We note a possible shortcoming of adopting the $\alpha$-disk model, namely the assumption that the turbulence is isotropic. It is known that numerical studies of at least one of the proposed hydrodynamical mechanisms that may drive turbulence in protoplanetary disk Ohmic (“dead”) zones (i.e., the VSI) have seen turbulent stresses that are clearly non-isotropic (Stoll et al., 2017), as well as in Dead Zones whose activity is driven by sandwiching MHD turbulent layers like in Yang et al. (2018). Incorporating such higher-order effects would require formulating a Reynolds averaging type of mixing scale model that takes into account the anisotropy of the shear stresses following the approaches found in Cuzzi et al. (1993); Dobrovolskis et al. (1999). This should be revisited in future analyses.

3.2. Equations of Motion and Steady State

We write the fundamental equations of motion in the local frame rotating at $\Omega$. The radial coordinate is $x$, the azimuthal coordinate $y$ and $z$ is the vertical coordinate. We represent the mean azimuthal shear as a departure from a mean Keplerian state $V_{0}=-(3/2)\Omega(r-r_{0}){\bf{\hat{y}}}$ (Umurhan & Regev, 2004). The disk gas velocity relative to the mean state is ${\bf v}_{g}=u_{g}{\bf{\hat{x}}}+v_{g}{\bf{\hat{y}}}+w_{g}{\bf{\hat{z}}}$ while the corresponding relative particle velocity is ${\bf v}_{p}=u_{p}{\bf{\hat{x}}}+v_{p}{\bf{\hat{y}}}+w_{p}{\bf{\hat{z}}}$. By writing $P=\tilde{P}+p$, we split the pressure field into a sum of the background field ($\tilde{P}$) – i.e., that drives the relative mean motion between the gas and particles – and perturbation field ($p$). The axisymmetric equations of motion for the gas are

[TABLE]

Momentum exchange between gas and particle phases arises in terms above of the type $\rho_{p}(u_{g}-u_{p})/t_{s}$, where we have also defined the parameter $\mu\equiv 1/\rho_{g}t_{s}$. It is assumed that the conditions in the disk (i.e., gas density and temperature as well as particle abundance) fall into the Epstein regime (though see Sec. 7.5). The background disk pressure gradient, the primary driver of the SI, is given by

[TABLE]

In order to account for the viscous torque induced by the background alpha-disk model, we include on the LHS of equation (5) the appropriate background viscous forcing in the form of an acceleration $-\alpha c_{s}H\rho_{g}\nabla^{2}{\bf v}_{g}$. The equations of motion of the the particle phase are,

[TABLE]

where the typical kinetic energy per unit mass of the particles induced by the turbulent stirring of the gas is given by c_{d}^{2}=\alpha c_{s}^{2}\big{/}\left(1+\tau_{s}^{2}\right). Steady uniform solutions of Eqns. (4-12) are sought assuming no vertical velocities and constant steady gas and particle densities, $\overline{\rho}_{g}$ and $\overline{\rho}_{p}$. Following the procedures of Nakagawa et al. (1986) and YJ2007 (their Eqns 7-8) we have that uniform gas velocities $U_{{}_{g0}},V_{{}_{g0}}$ are

[TABLE]

and the uniform particle velocities are

[TABLE]

(cf., Dipierro et al., 2018). Given the importance of the relative radial velocity between the gas and particle phases in rationalizing the SI in terms of resonant conditions (a la Squire & Hopkins, 2018a, b; Hopkins & Squire, 2018a, b), we find

[TABLE]

3.3. Linearized perturbations

We linearly perturb equations (4-12) around this steady state according to

[TABLE]

for the gas quantities and

[TABLE]

for the particles. Since the gas is incompressible, we further write the quantities $u_{g}^{\prime}$ and $w_{g}^{\prime}$ as derived from a perturbation streamfunction $\psi^{\prime}$:

[TABLE]

One can formally define the azimuthal gas and particle “fluid” vorticity fields as:

[TABLE]

The gas vorticity is related to the stream function via $\omega_{g}^{\prime}=\nabla^{2}\psi_{g}^{\prime}$. The perturbed quantities are then Fourier decomposed. For example, the streamfunction is written as

[TABLE]

where $k_{x}$ and $k_{z}$ are the radial and vertical wavenumbers (respectively), the frequency is $\omega$, and $\hat{\psi}$ is the normalmode amplitude.

$\omega_{i}\equiv{\rm Im}(\omega)>0$ indicates growth with a corresponding e-folding timescale, $t_{g}\equiv 2\pi/\omega_{i}$. We restrict our consideration to positive values of $k_{x}$ and $k_{z}$.999 In performing double Fourier transforms of real quantities, one is free to restrict attention to positive wavenumber values for one chosen spatial direction, here we choose to be the radial one ($x$). The linearized PDE’s possess z-reflection symmetry since imposing $z\rightarrow-z$ together with $w_{g}\rightarrow-w_{g}$ and $w_{p}\rightarrow-w_{p}$ leaves the equations invariant. The eigenvalues $\omega$ are insensitive to the sign of $k_{{}_{z}}$, confirming the same reflections made in YG2005. Thus, we may safely restrict attention to positive values of $k_{z}$ as well. Since $k_{x}>0$, values where ${\rm Re}(\omega)>0$ indicate outwardly propagating patterns. We are reminded that no such symmetry characteristic exists in the radial direction due to the imposed symmetry breaking provided by both the presence of turbulence and an externally imposed radially dependent pressure field $\tilde{P}$.

The combined system reduces (using Mathematica) to a generalized expression for the dispersion relation of the form

[TABLE]

in which $F$ is a sixth order algebraic equation in $\omega$. In the inviscid limit, the six temporal modes correspond to two inertial waves in the gas phase, two inertial waves in the particle phase, and two “zero”-temperature acoustic modes in the particle phase (e.g., see discussion in Chapter 10 of Chandrasekhar, 1961)101010 In the absence of turbulent stirring a pressure-less fluid could be viewed as a zero-temperature gas. With $c_{d}^{2}\neq 0\ll c_{s}^{2}$ the particle phase effectively behaves as a low temperature compressible fluid. . The algebraic equation for $\omega$ is solved using standard root-finding methods found in Matlab 2017a. The eigenvalue $\omega$ depends upon five parameter expressions: the two wavenumbers $k_{{}_{x}},k_{{}_{z}}$, the Stokes number $\tau_{{}_{s}}$, the density ratio $\epsilon=\overline{\rho}_{p}/\overline{\rho}_{g}$, and the ratio of the turbulent intensity parameter to the disk opening ratio squared: $\alpha/\delta^{2}$. In all of our parameter scans, we assume $\delta=0.05$, a typical value for nominal disk temperatures and orbit velocities. Dynamic lengthscales are normalized by $\delta c_{s}/\Omega=\delta H=L_{{}_{{\rm{SI}}}}$ and growth rates normalized by $\Omega$.111111In YG2005 and YJ2007, the lengthscales are quoted as “$\eta r$”, where $\eta$ is the radial pressure parameter of Nakagawa et al (1986), which expression is equivalent to $\delta H=c_{s}H/v_{k}$.

3.4. Turbulent dilution model

Finally, we restrict our choice of $\epsilon$ to be physically consistent with the idea that a turbulent disk will loft particles away from the midplane and, consequently, result in dilution of $\epsilon$ near the midplane. Following previous authors (Dubrulle et al., 1995; Carballido et al., 2006, 2011; Estrada et al., 2016), we estimate an effective particle scale height $H_{p}$ from the balance between particle settling toward the midplane and upward lofting by turbulent motions. Thus, given an initial local ratio of the particle surface mass density to gas surface mass density, $Z$, then in a region of vertical thickness $H_{p}$ we broadly assign a local particle to gas mass density ratio via the simple relationship:

[TABLE]

which, henceforth, is referred to as the turbulent dilution model (TDM). In this form, $H_{p}=\sqrt{\alpha/(\alpha+\tau_{s})}H$. A disk with a cosmic abundance of about 1 percent would correspond to $Z=0.01$. Even if a protoplanetary disk starts out with a cosmic abundance of solids, any given radial location within that disk may have values of $Z$ that vary with the disk’s evolution (Birnstiel et al., 2012; Estrada et al., 2016; Sengupta et al., 2019). As such, we allow for values of $Z$ departing from the fiducial cosmic abundance value.

4. Verification

As a robustness test, we set $\alpha\rightarrow 0$ and were able to faithfully reproduce all individual eigenvalues and eigenvectors quoted in Table 1 of YJ2007 as well as the growth rate diagrams shown in their Fig. 1. We indeed confirm that the SI is an instability of inertial modes of mixed character (i.e., particle-gas modes). We note that the calculation in YJ2007 assumed a compressible gas component while our model assumes the gas to be incompressible. The fact that the growth rates are essentially identical in both calculations strongly suggests that the SI (as considered in their study) is practically insensitive to gas compressibility and, further, it would be sound to examine its evolution with the assumption of an incompressible gas.

Figure 1 shows the maximum growth rates as a function of $k_{x}$ and $k_{z}$ which reproduces the quality reported in YG2005. There exists a combination of $k_{x}$ and $k_{z}$ for which the growth rates are locally maximal. According to Squire & Hopkins (2018a), there exists a wave-drift resonance relationship identifying this combination as those values of $k_{x}$-$k_{z}$ for which some collective fluid mode has a projected phase speed that resonates with the relative drift velocity between the gas and the particles. In the case of the classic SI, one such wavemode is an inertial wave. In cases for which both the mass-loading is weak (small $\epsilon$) and the coupling between gas and particles is strongish ($\tau_{s}<1$) inertial modes in a collective gas-particle medium can be approximated by the wave response in the gas assuming no coupling to the particles. Figure 1 shows the actual growth rates determined for such a strongly coupled weakly mass-loaded model. In this extreme case, it is elementary to show that

[TABLE]

(Lyra, & Umurhan, 2019). Since there are no $y$ disturbances, identifying the radial component of the inertial wave to the relative drift velocities means equating

[TABLE]

Inserting (22) and the appropriate steady radial drift expressions from (15) with $\alpha=0$ into the above expression, we find that the desired $k_{x}$-$k_{z}$ relationship is

[TABLE]

The wave-drift resonance relationship expressed in (24) is shown as a dashed line over the growth rates showcased in Figure 1. We see clearly that the resonance relationship follows the maximum growth rates as one scans along $k_{z}$. This lends confidence that the resonance condition is a very good predictor for identifying conditions corresponding to maximal growth of this instance of the RDI in the $\epsilon\ll 1$ regime.

5. Results

Our closed-form solutions (see Appendix A) permit a finely-resolved sweep in parameter space varying both the Stokes number $\tau_{s}$ and turbulence parameter $\alpha$. The particle-to-gas volume mass density ratio $\epsilon$ is automatically determined as a function of $\alpha$ and $\tau_{s}$ based on the global (unsettled) solids mass fraction $Z$ and the turbulent dilution model (TDM), eq. (21). For most parameter sweeps, we usually set $Z$ at nominal cosmic abundance, $Z=0.01$, but we also consider other values of $Z=$0.02,0.03, 0.04, and 0.08. For these values of $Z$ but with a fixed disk opening angle $\delta=0.05$, we study the fastest growing mode as a function of $\alpha$ and $\tau_{s}$.

We caution against cavalierly linking the results of this work to the original theoretical studies (e.g.,YG2005 and YJ2007) which surveyed the inviscid linear stability calculation for ${\cal O}\left({1}\right)$ values of $\epsilon$. We observe that by adopting the TDM the inviscid limit is singular as it predicts $\epsilon\rightarrow\infty$. Indeed the TDM assumes that some type of quasi-steady state has emerged between the particles and the surrounding turbulent state. The only permissible pathway linking results of this model to those of the inviscid limit is to take the double limit $\alpha\rightarrow 0$, $\tau_{s}\rightarrow 0$, while setting $\tau_{s}/\alpha$ to a finite constant. The latter allows for choosing an arbitrary value of $\epsilon$, but the limit corresponds to dynamics involving particles instantaneously responding to the gas motions with no relative drift. In other words, this limit represents regular inviscid incompressible single-fluid gas dynamics with a slightly enhanced mean density.

Our most basic result is this: isotropic turbulence, as measured by $\alpha$, causes the growth rates of the SI to diminish, while also increasing the wavelengths corresponding to fastest growth. This is not surprising, because the shortest wavelength modes are eaten away by turbulent diffusion of momentum and particle concentration.

5.1. Individual model results

In this section we show the properties of individual models. In particular, Figure 2 displays contour plots of growth rates as a function of $k_{x}$ and $k_{z}$ for three values of $\alpha$ together with $\tau_{s}=1$ fixed.

For weak turbulence there exist wavelengths corresponding to maximum growth that are both very short, and have growth rates on the order of the disk rotation frequency. For example, the top panel of Fig. 2 shows results for the very low value of $\alpha\approx 4\times 10^{-6}$ with $\tau_{s}=1$, perhaps corresponding to “midplane turbulence” around a settled particle layer in a globally laminar nebula. The wavelengths of maximum growth are $\lambda_{x}({\rm max})=2\pi/k_{x}({\rm max})\approx 0.037H$ and $\lambda_{z}({\rm max})=2\pi/k_{z}({\rm max})\approx 0.027H$; thus, only a little smaller than $\delta H$. This fastest growing mode has an e-folding timescale $t_{gm}\approx 0.5{\rm P}_{orb}$, where $t_{gm}\equiv t_{g}({\rm max})\equiv\left[\Omega/\left(2\pi\omega_{i}({\rm max})\right)\right]{\rm P}_{orb}$. As the intensity of turbulence increases (middle and lower panels), the wavelengths corresponding to maximal growth get larger and the corresponding growth rates diminish.

The middle panel of Fig. 2 shows a wavenumber survey for $\alpha\approx 8\times 10^{-5}$ and $\tau_{s}=1$ – a so-called weak-to-moderately turbulent model, even if the particle layer is rather densely settled. The fastest growing wavelengths in both directions are nearly equal with $\lambda_{x}({\rm max})\approx 0.25H$ and $\lambda_{z}({\rm max})\approx 0.34H$, and start becoming of the same order of magnitude as the local disk scale height. The corresponding growth timescale is now considerably longer with $t_{gm}\approx 6.9{\rm P}_{orb}$.

The bottom panel of Fig. 2 similarly shows a wavenumber survey for $\alpha\approx 10^{-3}$ and $\tau_{s}=1$, a model we term strongly turbulent. Even here, according to the TDM, the particle layer has settled to a thickness of only $H/30$ because of the large $\tau_{s}$. The wavelengths of fastest growth become even larger, and the relative ($x,z$) length scales become more disparate with $\lambda_{x}({\rm max})\approx 1.55H$ and $\lambda_{z}({\rm max})\approx 7.85H$ (implying vertical sheet-like disturbances). The corresponding e-folding growth timescale is $t_{gm}\approx 10.6{\rm P}_{orb}$.

The pattern propagation of the turbulent SI also depends upon the degree of turbulence. We measure the radial pattern speed to be given by $c_{r}\equiv\omega_{r}/k_{x}$, where $\omega_{r}={\rm Re}(\omega)$. The pattern propagation of the fastest growing mode is denoted by $c_{r}({\rm max})$ and equal to $\omega_{r}({\rm max})/k_{x}({\rm max})$. With reference to Fig. 2 we see that the pattern propagation is inward for the two largest values of $\alpha$ shown while it is outward for the nominally weakly turbulent model.

The most strongly turbulent model shown in Fig. 2 closely corresponds to the conditions modeled in Johansen et al. (2007). Figure 1(b-c) of Johansen et al. (2007) depicts the growth of the SI in a MRI driven turbulent setting in which $\alpha={\cal O}\left({10^{-3}}\right)$ and where the Stokes number $\tau_{s}\sim 1$. A cursory examination of the growing modes in those simulations during the early linear growth phase ($t<30$Porb) shows that radial wavelength of the most prominently growing structure is $\sim$1.2$H$, with a growth timescale of 10-15 Porb. Furthermore, our model predicts that the pattern speed $c_{r}({\rm max})=0.12H/$Porb and inward. The apparent propagation of the growing mode in Johansen et al (2007) is also inward with a pattern speed $~{}0.1H/$Porb (we discuss pattern speeds further in sec. 5.3). While this certainly does not prove that our simple model is sufficiently predictive, it is encouraging that it predicts features that are in both qualitative and approximate quantitative agreement with previously published numerical simulations.

5.2. Growth: general survey

Figure 3 shows the growth timescale of the fastest growing mode, as a function of $\tau_{s}$ and $\alpha$ for $Z=0.01$. There are several notable results. There exists a critical branch line, defined nearly by $\epsilon=1$, in which the growth timescales are infinite. This critical curve extends from $\tau_{s}=0$ and terminates at a critical value of the Stokes parameter which we denote by $\tau_{c}$. For the parameter combination considered here ($Z=0.01$, $\delta=0.05$), $\tau_{c}\approx 0.45$ at $\alpha=\alpha_{c}\approx 3.3\times 10^{-5}$. Although not perceptible for the $Z=0.01$ case shown in Fig. 3, the actual location of $\tau_{c}$ corresponds to a value of $\epsilon\approx 1.1$ based on the turbulent dilution model. Similar growth timescale plots, calculated for several different values of $Z$, are shown in Figure 4. The corresponding value of $\tau_{c}$ more clearly corresponds to increasingly larger values of $\epsilon$ as $Z$ increases. The critical line appears to hug the $\epsilon=1$ line until $\tau_{s}$ begins to approach $\tau_{c}$ from below, whereupon the curve bends downward in $\alpha$ forming a beak-like shape (this is most starkly apparent for $Z=0.08$ in Fig. 4). This critical point always corresponds to values of $\epsilon$ larger than 1. We have summarized the observed trends in Table LABEL:Critical_Values_alpha. Generally speaking, $\tau_{c}$ is some function of $Z$ and $\delta$, but a theory clarifying the meaning of this point and its mapping as a function of these and other parameters is not undertaken here.

The character of the turbulent SI is sensitive to whether or not $\tau_{s}<\tau_{c}$ or $\tau_{s}>\tau_{c}$. In the case where $\tau_{s}<\tau_{c}$, the growth rate dramatically depends on which side of the branch line one is on. For the region below the branch line (i.e., for $\epsilon>1$) – the so-called laminar zone (Zone I) – the growth timescales are generally fairly short (less than orbit times) in broad accordance with the low turbulence results depicted in the top panel of Fig. 2 as well as in line with expectations based on published numerical and theoretical studies examining the SI in the inviscid limit (also see introductory discussion of section 6).

On the other hand, for $\tau_{s}<\tau_{c}$ and above the branch line (i.e., $\epsilon<1$), in the so-called turbulent regime (Zone II), the growth timescales are extremely long - anywhere from tens to thousands of local orbit times. For optically thick disks in which $\tau_{s}=0.01$ and where the VSI is operating $\alpha\sim 2\times 10^{-4}$ (Estrada et al., 2016; Malygin et al., 2017), the growth timescales are just under $10^{4}$ orbit times (at Jupiter this corresponds to about $10^{5}$ years). Approaching the $\epsilon=1$ line from either side of the branch line results in growth timescales approaching infinity. In other words, at $\epsilon=1$ the mode is marginal, neither growing nor decaying.

For $\tau_{s}>\tau_{c}$ the fate of the linear SI is different. The branch line $\epsilon=1$ ceases to have any consequence for the growth rates. In this region the growth rates are relatively fast, anywhere between tenths and tens of orbit times. There appears to be a boundary that separates the turbulent zone from the laminar zone (in this range around $\alpha\sim 10^{-4}$) but this is apparent only when looking at the neutral pattern propagation speed line (discussed further below). For $\alpha>10^{-4}$ (Zone III) the growth timescales are about tens of orbit times while for $\alpha<10^{-4}$ (Zone IV) the growth timescales are even 10-100 times shorter. Zone III is further distinguished from Zone IV in the character of the pattern speeds and propagation directions (next section). Finally, in a broad sense, we identify Zones II-III as comprising the “turbulent regime” since they embody regions of relatively large values of $\alpha$, while we refer to Zones I,IV as comprising the “laminar regime”.

We surveyed our results to test whether or not the incompressibility assumption remains valid. For a given maximally growing mode characterized by wavelength $\lambda$ we find that both the growth timescale and pattern propagation timescale (across radial distance $\lambda$) are always much longer than the sound propagation time across the same lengthscale, consistent with the physical basis for neglecting compressibility in the gaseous phase.

5.3. Pattern Speeds – general survey

Figure 5 depicts the pattern propagation speed of the fastest growing modes for the same parameter sweep discussed above. We restrict our attention to $Z=0.01$ noting that the qualitative character we report here carries over to the other values of $Z$. We immediately observe that the pattern speed is 0 along the branch line ($\epsilon=1$) separating Zone I from II. This is consistent with predictions made for the inviscid limit (YG2005) wherein the pattern speed is zero for $\epsilon=1$. However, beyond the critical point $\tau_{c}$, the 0 pattern speed curve does not lie on the $\epsilon=1$ line but, instead, follows the curve designating 0 pattern speed which roughly separates Zone III and IV – see dashed line of Fig. 5.

In analyzing the inviscid SI, YG2005 predict that its pattern speed depends on whether or not $\epsilon>1$. We find that those trends carry over to this turbulent model: the pattern speeds are outward ($c_{r}>0$) in laminar Zone I ($\epsilon>1$) while they are inward ($c_{r}<0$) for turbulent Zone II. Pattern speeds in this regime are typically less than $0.01H/P_{{\rm orb}}$. A stark qualitative distinction appears in examining $c_{r}$ in Zones III and IV. Within the relatively turbulent Zone III $c_{r}<0$; that is, inward, just as in Zone II. However the pattern speeds are very high, on the order 0.1-0.2 $H/P_{{\rm orb}}$. This high drift rate was already observed in our earlier analysis (sec. 5.1) of Johansen et al. (2007)’s simulations. Meanwhile, in Zone IV the outwardly propagating patterns also drift with relatively high speeds ( $<0.05H/P_{{\rm orb}}$) that are still, however, slower than those of Zone III. In either case, such high pattern speeds means that such structures might drift out of the region of interest on relatively short timescales, before they are able to produce planetesimals– see further discussion in sec. 7.2.

5.4. Mode structure – general survey

Similar to the previous sub-section, we restrict our discussion to $Z=0.01$ noting that the qualitative character we report here carries over to the other values of $Z$. In Figures 6 and 7 we plot the structure of the fastest growing modes, where Fig. 6 shows its radial wavelength $\lambda_{x}\equiv 2\pi/k_{x}({{\rm max}})$ while Fig. 7 displays its aspect ratio defined to be $\lambda_{z}/\lambda_{x}\equiv k_{x}({{\rm max}})/k_{z}({{\rm max}})$.

In accordance with the trends predicted in the inviscid limit, within the weakly turbulent regions (Zones I,IV) the radial wavelengths of the fastest growing modes are on the order of $H$ and get increasingly shorter as $\alpha$ decreases. In Fig. 7 we indicate the location where $\lambda_{z}/\lambda_{x}\approx 1$; this lies mostly in Zone IV with some spillover into Zone I. Scanning around this region we see that $\lambda_{z}/\lambda_{x}$ indeed remains ${\cal O}\left({1}\right)$ implying that the mode structure is that of axially symmetric ring structures in this rough patch of parameter space.

On the other hand, in Zone II, above the branch line $\epsilon=1$, $\lambda_{x}$ steadily increases with $\alpha$. This trend continues into Zone III where the turbulence is still relatively strong but $\tau_{s}>\tau_{c}$. Moreover, the aspect ratio rapidly goes from flattened or tubular rings into vertically oriented sheets as $\lambda_{z}/\lambda_{x}$ steadily grows with increasing turbulent intensity. In both Figures 6 and 7, we have indicated regions for which $\lambda_{z}/\lambda_{x}>5$; this property represents the entirety of Zones II and III. Naturally, caution should be exercised in interpreting the results in the limit that $\lambda_{z}$ greatly exceeds several scale heights. That is, when the predicted vertical lengthscales greatly exceed $H$ throughout the majority of Zone II, the growth timescales indicated throughout Zone II of Figures 3-4 may be lower bounds because there may be no expression of the SI for this parameter regime.

It is, however, encouraging that the trends predicted here appear to be consistent with the simulation results reported in Balsara et al. (2009). In their Figure 7, for particle sizes in which the SI appears to manifest, the vertical structure in the particle density appears uniform up to the scale height of the particles themselves.

6. Comparison to some recent published studies

It is revealing to compare the predictions of the theory with published numerical simulations where SI is (filled symbols), or is not (open symbols), seen to operate. The values of $(\tau_{s},\alpha)$ for the turbulent runs are given explicitly in the papers cited (Johansen et al., 2007; Balsara et al., 2009) and we have indicated their approximate locations on the growth timescale plot Fig. 3. In particular, for these simulations (with $\delta=0.05,Z=0.01$) SI is only manifest squarely inside of Zone III, for relatively large $\tau_{s}$. We note that these simulations considered several particle species with differing Stokes numbers while the growth rates shown in the Figure correspond to the linear theory inferred from assuming a single size species into which all mass is assigned. We consider our theoretical predictions of growth rates to be upper bounds for simulations conducted with several species simultaneously present, based on the recent findings of Krapp et al. (2019). In the following we examine several recent studies that provide enough information for us to compare their results with our theoretical predictions.

6.1. Yang et al. (2017)

The numerical experiments run by Yang et al. (2017) are more nuanced. They considered a single particle species of a given Stokes number initiated in a globally nonturbulent nebula - effectively laminar flow.121212The simulations of Bai & Stone (2010) are similarly set up, however we do not analyze them here because no spacetime plots of vertically integrated density were provided. A follow-up study examining these in more detail is warranted. They initialize the particles with a Gaussian vertical distribution symmetric about the midplane with a corresponding effective scale height $H_{p}=0.2H$ and allow the particles to settle toward the midplane. In these simulations and others like it (Bai & Stone, 2010) the layer collapses until the particle scale height reaches a minimum “bounce” value, out of which the SI begins to become manifest. We conjecture that this first bounce corresponds to the rapid transition of the fluid into an unsteady (possibly turbulent) state. Whether or not the dynamically active fluid state during this first bounce is driven by the SI or some other mechanism like the Kelvin-Helmholtz instability remains to be established. Bai & Stone (2010) state that the long-time development of the unsteady particle layer under their conditions is not maintained by Kelvin-Helmholtz overturn but, rather, by the SI itself. However, in a footnote, they allow that this would not be the case for smaller St. Indeed, the recent study of Gerbig et al. (2020) shows that whether SI or the Kelvin-Helmholtz instability maintain the vertical diffusion of particles depends upon the value of St in concert with $\delta$, with the latter process being important for values of $Z$ in the vicinity of $0.01$ for $\delta=0.05$, or toward higher values of $Z$ as $\delta$ increases .131313 This study was released during the revision phase of this article. Also, the streaming parameter $\delta$ is $2\beta$ in their study. Deeper understanding of this so-called early bounce phase requires further analysis. In either case, we conjecture that it is out of an unsteady-possibly-turbulent state that the subsequent SI develops in these simulations, and it may in turn drive yet another component of turbulence. We do not think these candidate mechanisms can be, nor need to be, distinguished at this stage. To keep our terminology as generic and descriptive as possible, we refer to the fluid state during this early bounce phase as “midplane turbulence”.

In light of the preceding discussion, for all the simulations presented Yang et al. (2017) record the characteristic value of $H_{p}/H$ at every time step. In practice, the particles settle toward a nominal minimum base state value of $\left(H_{p}/H\right)_{{\rm s}}$ (the aforementioned “bounce”) representing the first instance of balance between turbulent stirring and midplane directed settling, and it is from this base state that the SI begins to emerge. For values of $\tau_{s}\approx 0.01$ this state appears to be reached by 20Porb while for $\tau_{s}\approx 0.001$ the base state is reached at more like 80Porb.

Thus, to place the “detected SI” results of Yang et al. (2017) onto the theoretical growth timescale plots of Figs. 3-4, we had to infer an effective value of the $\alpha$-parameter, henceforth $\alpha_{{\rm eff}}$, characterizing the near-midplane turbulence. Following the relationship between $H_{p}$ and $H$ found after Eq. (21) defining the TDM, under the assumption $\tau_{s}\gg\alpha_{{\rm eff}}$, we approximate $\alpha_{{\rm eff}}\approx\tau_{s}\left({H_{p}}/{H}\right)_{{\rm s}}^{2}$, yielding a corresponding $\epsilon_{{\rm eff}}\approx Z\sqrt{\tau_{s}/\alpha_{{\rm eff}}}$. This procedure works well for all of the simulations presented in Yang et al. (2017) that present a time-series for ${H_{p}}/{H}$. In the two cases where this information is not provided, we estimate an average midplane value of $\epsilon_{{\rm est}}=\rho_{p}/\rho_{g}$ by matching the color of the quantity against the provided color bar. Equating $\epsilon_{{\rm eff}}$ to $\epsilon_{{\rm est}}$, and subsequently inserting it into the rewritten TDM, leads to \alpha_{{\rm eff}}=\tau_{s}Z^{2}/\big{(}\epsilon_{{\rm eff}})^{2}. We have indicated the parameter locations of most of the simulations reported in Yang et al. (2017) throughout Figs. 3-4.

With $\alpha_{{\rm eff}}$ we computed the e-folding timescale of the fastest growing mode, $t_{gm}$, and its associated radial wavelength $\lambda_{{}_{xm}}\equiv\lambda_{x}({{\rm max}})$. The observed e-folding growth timescales, $t_{o}$, in these simulations are difficult to assess based on the graphs provided. Thus we estimate $t_{o}$ to be less than the observed saturation timescales (hereafter, $t_{{\rm sat}}$) found tabulated in Table 2 of Yang et al. (2017). We visually determined an effective lengthscale of structures observed to emerge during the linear phase (hereafter, $\lambda_{{}_{x,{{\rm obs}}}}$) according to the following procedure: Figures 3 and 7 of Yang et al. (2017) show spacetime plots for the azimuthally averaged, vertically integrated, scaled particle density which they denote by \big{<}\Sigma_{p}\big{>}/\Sigma_{g,0}. Approximately after the time particles have settled near the midplane ($t=20$Porb or $t=100$Porb depending upon $\tau_{s}$, see above) but well before $t_{{\rm sat}}$, we count the number ($N_{p}$) of peaks in \big{<}\Sigma_{p}\big{>}/\Sigma_{g,0}, across the radial size of the simulation domain, $L_{x}$, and then we estimate $\lambda_{{}_{x,{{\rm obs}}}}\approx L_{x}/N_{p}$. In practice, we focused on counting resolvable peaks within the first two predicted e-folding growth timescales because once coherent filamentary structures emerge they appear to nonlinearly interact, resulting in a series of mergers before the saturated state becomes manifest. In some cases, it was difficult to resolve an unambiguous number of peaks and this uncertainty was noted. In all cases we counted peaks as soon as the structures appeared coherent to the eye.

The results of this activity are summarized in Table LABEL:Yang_etal_2017_simulations. Despite the crudeness of this approximate analysis we find that our theory predicts the general properties reported in Yang et al. (2017). In particular, it is clear for those Yang et al (2017) initial conditions which did manifest SI, that the effective $\epsilon_{{\rm eff}}$’s were always greater than 1, consistent with our predictions that growth is relatively strong for such conditions. In the two cases discussed by Yang et al. (2017) in which the SI was not observed, we predict that the growth times for those input parameter conditions are much longer than the time for which those simulations were conducted. Of interest is the case $\tau_{s}=0.001,Z=0.03$ where no structures were observed to emerge up to simulation time t=5000Porb. According to our predictions, if the simulation were run past $t=5\times 10^{4}$Porb then the SI should become apparent. Similar reasoning also applies to their $\tau_{s}=0.01,Z=0.01$ simulation. We also note that their 2D simulation for the case $\tau_{s}=0.01,Z=0.04$ for the box $\sim 0.4H\times 0.4H$ experiences a brief dip in $H_{p}/H$ to $\approx 0.004$ before leveling out to $H_{p}/H\approx 0.01$ at $t\approx 80$Porb (see red dotted time series shown in middle panel row of their Fig. 2 as well as right panel of their Fig. 3b). During this initial bounce phase ($t<50$Porb) the number of filamentary structures is impossible to assess. As such, we examined the behavior and character of this simulation after the initial response phase passed (after $t=50$Porb), when the number of filaments was first countable. For the corresponding value of $H_{p}/H$ ($\approx 0.01$), the predicted $\lambda_{xm}$ better matches the corresponding measured $\lambda_{{}_{x,{{\rm obs}}}}$.

6.2. Li et al. (2018)

Li et al. (2018) examine the fate of the SI under similar circumstances and initial setups considered in Yang & Johansen (2014) but using three different boundary conditions – periodic, reflecting and outflow – and they show that the nonlinear SI state is largely insensitive to the boundary conditions employed. Using a fixed grid resolution ($\sim H/320\approx 0.003H$) Li et al. (2018) checked for the emergence of filaments for three different box sizes: $(0.2H)^{3},(0.4H)^{3},(0.8H)^{3}$. Their simulations were run with $\tau_{s}=0.314$ and they quote time series quantities similar to $H_{p}/H$. During the early linear “first-bounce” phase, all simulations collapsed into a layer with $H_{p}/H\sim 0.001$ within $t=1$Porb. The predicted maximum growth timescale for these input parameters are short, i.e., $t_{gm}\approx 0.5$Porb. This layer then begins to develop the SI just as in the simulations reported in both Yang & Johansen (2014) and Yang et al. (2017). Coherent filaments undergoing epicyclic oscillations are clearly discernible by $t=3$Porb and it is also at this point they undergo nonlinear merging.

We perform the same crude analysis as above based on the periodic boundary condition runs (see Fig. 10 of Li et al., 2018) and the results are summarized in Table LABEL:Li_etal_2018_simulations. We estimated $\lambda_{{}_{x,{{\rm obs}}}}$ by counting the number of filaments around $t=2$Porb. We choose this time because filaments are not easily identifiable at any time before this point, while for times after this point the filaments begin their process of merging. For all simulations run, we find $\lambda_{{}_{x,{{\rm obs}}}}\approx 0.015H$, while our theory predicts a value approximately half of that, $\lambda_{xm}\approx 0.007H$. This means that the observed fastest growing mode encompasses about 5 grid points while our predicted value would not be resolvable (at between 2 and 3 of their gridpoints). Indeed, we note that their measured value for $H_{p}/H$ at $t\approx 1$Porb is much smaller than their grid resolution. We conjecture that their simulations are not sufficiently resolved and should be run for at least 2-3 times the resolution originally used. Convergence would perhaps be indicated if during this early organization phase ($t<1$Porb) the minimum value of $H_{p}/H$ remains fixed with increasing resolution. What is not yet clear is how the response of the particles, and the turbulence they churn up as they approach the midplane, depend upon resolution as such a systematic study remains to be done.

6.3. Gerbig et al. (2020)

Gerbig et al. (2020) consider the fate of SI under setups and conditions similar to those of the previous two studies. One of their main aims is to disentangle to what degree the midplane turbulence churned up by the settling dust particles is driven either by the Kelvin-Helmoltz instability or the SI, and they seek to shed light on this outcome as a function of the dust streaming parameter, $\delta$ (which is equal to $\beta/2$ in their analysis), $Z$, and for a fixed value of $\tau_{s}=0.2$. The results they uncover are subtle but it is evident from their reported simulations that when SI is not present, the dust layer thickens consistently with what one might predict based on the Richardson criterion for Kelvin-Helmholtz overturn under the influence of gravity. They provide the results for a suite of runs in which they quote a measured dust layer thickness plus its standard deviation, together with snapshots of the particle accumulation as viewed from several perspectives. For one run, $Z=0.02,\tau_{s}=0.2,\delta=0.05$, they also provide a spacetime (radius-time) plot of the vertically integrated , azimuthally averaged particle density. Focusing on only a subset of simulations reported in their study – all before self-gravity is turned on – we compare their measured and tabulated quantities against predictions made using our theory. As per the procedure defined earlier, we estimate an $\alpha_{{\rm eff}}$ based on their reported value of $H_{p}$. We also apply our theory to corresponding “high” and “low” values of $\alpha_{{\rm eff}}$ based on their quoted one standard deviation values of $H_{p}$. The results of this exercise are summarized in Table LABEL:Gerbig_etal_2020_simulations.

We find reasonable agreement between our theoretical predictions and their reported study. The predicted versus observed wavelengths are mutually consistent with one another within one standard deviation of their reported values of $H_{p}$. Most of the filament wavelength structure that emerges in their simulations when SI is present is shown at late times ($t=40$P), likely well after a significant amount of filament merging has taken place. For this reason we suppose that the observed wavelength/average-separation ($\lambda_{{}_{x,{{\rm obs}}}}$) should be larger than the predicted wavelength of our theory. Given the error bounds on their reported $H_{p}$, this appears to be a plausible interpretation for the simulation with $Z=0.03,\tau_{s}=0.2,\delta=0.05$, in which our predicted wavelength ($0.05H<\lambda_{{}_{xm}}<0.08H$) is smaller than the observed value ($\lambda_{{}_{x,{{\rm obs}}}}=0.10H$). We also observe that for the simulation with $Z=0.01,\tau_{s}=0.2,\delta=0.05$, in which SI is not observed, our theory predicts that the fastest growing wavelength should be $\lambda_{{}_{xm}}=1.1H$ with a growth timescale of $\sim 70$P. Given that their simulation was run in a box with radial length equal to $0.4H$, we predict that if the same simulation was rerun with a radial extent of at least 1.1$H$ or longer, then SI should appear. A similar prediction can be made for their $Z=0.02,\tau_{s}=0.2,\delta=0.10$ simulation, although the radial scale of their box should be at least twice as much, if not more (see Table LABEL:Gerbig_etal_2020_simulations).

6.4. Yang et al. (2018)

Yang et al. (2018) examine the response of the SI in a setup that supports two scenarios: one in which the full layer is turbulent due to the MRI (i.e., Ideal MHD, and hereafter “iMHD”), and the other as a Dead Zone (“DZ”, hereafter) in which the midplane is Ohmic and effectively MHD inactive but where the disk gas about 1 gas scale height away from the midplane is MRI active. They consider a particle component with $\tau_{s}=0.1$ and metallicities of $Z=0.01,0.02,0.04$ and $0.08$. In the following we restrict our attention to the iMHD models. We do not consider their DZ models because the SI physics contained in them are probably strongly influenced by coherent structures, a feature which is outside purview of our theory – see Appendix B for further discussion.

For the iMHD models we repeat the calculations done in the previous subsections and the results are summarized in Table LABEL:Yang_etal_2018_simulations. We estimate $H_{p}$ for each simulation based on what we were best able to surmise to be the particle collective’s first bounce based on the time series found in the first column of their Figure 9. The value of $\alpha_{eff}$ is different than the value quoted in their Table 1 which is based on the state of the flow at very late stages of their simulation. For all of the metallicities considered we find that the fastest growing modes have wavelengths ($\sim 30$) much larger than the radial computational domain size ($\sim 4$H). The predicted growth timescales are several hundred orbit times (at the very least) which is far longer than the 100 P timescales on which their simulations were conducted.

However, inspection of their results shows that the iMHD model achieves very modest increases in particle enhancements over and above the nominal base value of either the vertically integrated particle densities, i.e., $\Sigma_{p}/\Sigma_{p,0}$ where enhancements range from a factors of 2 (on average) to no more than factors of 5 in extremes (see Table 3 of their study). The particles are prevented from settling to the midplane ($H_{p}/H\sim 0.3$) owing to the large vertical fluid stresses arising from the MRI turbulence with effective values of $\alpha\sim 0.008$. The top row of Figure 10 of Yang et al. (2018) depicts a spacetime diagram showing $\Sigma_{p}$ as a function of radial position and time for the four iMHD models of differing $Z$ compared against a sample run in which there is no back-reaction of the particles back onto the fluid. Except for the $Z=0.08$ model, inspection of the remaining figures show that there is no obvious distinction between those simulations with and without back-reaction, as there is no clear nor sustained emergent radial structure beyond the periodic box scale. In the $Z=0.08$ model (far right, top-row panel of their Figure 10), there is a brief period (between $t/P_{orb}=20$ and $40$) in which strong accumulation appears to take place but then it eventually dissipates.

We interpret the possibility of ephemeral bursts of accumulation to be the result of concentrations driven into place by large scale coherent structures of the flow field of the turbulent state, or perhaps even by turbulent concentration (Hartlep & Cuzzi, 2020). For example, the driving motions of the MRI (channel modes) and its secondary turbulent transition (parasitic/Kelvin-Helmholtz overturn) are primarily axisymmetric which means that the corresponding fluid eddies and waves of the driving scales are azimuthally aligned – meaning that there are coherent structures in the flow and pressure field being imposed upon the flow by the unstable MRI activity. Such motions may be responsible for any temporary accumulations seen in these iMHD models, instead of SI per se. Moreover, our theory is unable to treat nor predict what happens to the SI in the presence of unsteady yet coherent structures: the simple turbulence model like we have employed here assumes the motions are uncorrelated (implicitly placing it within the inertial range of a turbulent flow). See more discussion in Appendix B. .

7. Turbulent SI validity regime constrained by realistic protoplanetary disk models

We now put these results into the context of realistic global disk evolution and particle growth models, and other observed constraints. Recall that our theory assumes a single dominant, mass-carrying particle size. Benítez-Llambay et al. (2019) and Krapp et al. (2019) as well as others note that when multiple particle species are present in the mix, the overall growth rate of the SI may be significantly modified (probably diminished) as compared to single-sized particles (Krapp et al., 2019).

Thus, we ask the questions: for what turbulent disk conditions and particle properties does the SI provide a direct path to planetesimal formation, and, are these initial conditions realistic and self-consistent? We consider this question for three general locations in the solar nebula, each of which nominally representing: (i) the inner disk at $R\approx 3$AU, $T=265$K with $c_{s}\approx 1050$ m/s (assuming a H2 gas), (ii)the snowline around $R\approx 5$AU, $T=140$K with $c_{s}\approx 760$ m/s, and (iii) the outer disk at $R=30$AU, $T=73$K with $c_{s}\approx 550$ m/s.

7.1. Disk Lifetime Constraints

For the inner disk, we rule out SI models that predict growth timescales that are significantly longer than $\sim 1$ Ma, based on the evidence that planetesimals were abundantly forming before that time (Kruijer et al., 2017; Scott et al., 2018, see also Section 1). This translates to some number of equivalent Porb depending upon where the turbulent SI theory is being modeled. For example, we nominally place the snowline at $5AU$ which means that model parameters that predict growth time scales in excess of $(0.4-2.0)\times 10^{5}$ Porb are ruled out.

7.2. Particle and Mode Radial Drift

Growing particles will drift radially at different rates due to variable coupling with the nebula gas. The drift rate increases with Stokes number $\tau_{s}$, reaching a peak when $\tau_{s}=t_{s}\Omega\sim 1$, and then decreases for larger sizes. A Stokes number of unity corresponds to different size particles at different places in the protoplanetary disk, depending on the local gas surface density. For the standard Minimum Mass Solar Nebula (MMSN), meter-sized particles drift the fastest in the terrestrial planet region, whereas further out in the disk, where gas densities are low and the pressure gradients can be strong (especially if there is a strong gradient in the gas density at the disk outer edge), much smaller (mm-size) particles have $\tau_{s}=1$. In fact, inward drifting particles may drift faster than they can grow in size by sticking; this is the so-called “radial drift barrier” (e.g., Brauer et al., 2008; Birnstiel et al., 2012; Estrada et al., 2016; Sengupta et al., 2019).

A similar criterion can be applied to the SI. Using the mean radial velocity of the particle component (eq. 14) to estimate the time $t_{d}$ a particle takes to traverse the scale of the disk, we find (for $\tau_{s}\gg\alpha$)

[TABLE]

For a given set of parameters, the SI is considered “viable” if the derived growth rate is faster than the drain rate, i.e., when $\omega_{i}>t_{{}_{d}}^{-1}$. Growth rates depend very much on $\epsilon$, but sufficiently high solids-to-gas ratios that allow the solids to drive the gas motions are hard to achieve in turbulence (e.g., see Estrada et al., 2016), unless one imposes arbitrary trapping mechanisms such as pressure bumps to thwart radial drift (e.g., Drazkowska et al., 2013, 2014). Fractal particle growth leads to highly porous particles which drift radially much more slowly (Ormel et al., 2007; Estrada & Cuzzi, 2008; Okuzumi et al., 2012; Estrada et al., 2020); this may provide a means to weaken the radial drift barrier and generate the necessary solids enhancements, but decreasing $\tau_{s}$ also shifts the case to the left in Fig. 3, generally weakening SI for any $\alpha$.

As in the simulations reported in Johansen et al. (2007), when $\tau_{s}\sim{\cal O}\left({1}\right)$ the pattern drift of growing modes is relatively fast and there emerges the possibility that the growing mode drifts in toward the star faster than the overdensity can grow sufficiently for gravitational instability to take root. This concerns the notion of “convective instability” familiar in hydrodynamics (Drazin & Reid, 2004; Regev et al., 2016) and plasma physics (see Chapter 18 of Bers , 2016). As such, we assess the conditions in which the pattern drift timescale (to reach the star) is much shorter than the unstable growth timescale. The time it takes for inwardly propagating structures to drift into the star is $t_{w}=r/\left|c_{r}\right|$ (provided $c_{r}<0$); so if $\omega_{i}>t_{{}_{w}}^{-1}$, we consider the SI as viable.

7.3. Fragmentation Barrier

everal studies (Brauer et al., 2008; Birnstiel et al., 2012; Estrada et al., 2016; Sengupta et al., 2019) show that the fragmentation barrier in a turbulent medium may be estimated by assessing the inequality

[TABLE]

The above expression represents the condition in which the kinetic energy of colliding particles per unit mass as driven by turbulent eddies, approximately $2\tau_{s}\alpha c_{s}^{2}$ (Voelk et al., 1980; Ormel & Cuzzi, 2007) exceeds their binding energy per unit mass, which is quantified by $U_{f}$, the fragmentation speed for the particle (in reality a loose aggregate) in question. For silicate aggregates $U_{f}\approx 1.4$ m/s (Zsom et al., 2010; Güttler et al., 2010).141414We note, however, that it has been suggested that the fragmentation speeds of relatively high temperature ($T\sim 200-400$ K) aggregates composed of organic covered silicate sub-micron grains may approach 100’s of m/s (Homma et al., 2019). For the above quoted temperature range, a 0.03 $\mu$m monomer Si aggregate with about equal amount (by radius) of organic covering has a fragmentation velocity of $\sim 100$m/s , while a 1$\mu$m sized Si aggregate with about 3% organic coating (by radius) has a smaller fragmentation velocity diminish of about 10 $m/s$ (see their Figure 4), with the assertion that this reduction in critical speeds continues for larger masses. While it is unlikely that such sub-micron sized aggregates collide with one another with such high speeds in the turbulent nebula, this structural strengthening feature of organic covered Si grains should nonetheless be incorporated into future analyses.

The value of $U_{f}$ for H2O ice aggregates is more nuanced. Laboratory studies of material properties suggest that sticking is significantly more effective for H2O ice than for silicates (Bridges et al., 1996; Wada et al., 2009), allowing ice particles to grow larger, faster and more porous, with effective fragmentation velocity thresholds an order of magnitude or more higher than for silicates (Okuzumi et al., 2012; Wada et al., 2009, 2013). On the other hand, Musiolik & Wurm (2019) report on experimental results testing the surface energy of $\sim 1$mm H2O ice spheres in the $175$K$<T<240$K temperature range and find that this increased strength only applies in a narrow temperature range plateauing near 200 K, with effective sticking speeds almost 30 times lower when the grains are cold (i.e., $<175$ K). They find that at these colder temperatures, the surface energy of H2O ice grains is about the same as the surface energy of silicate particles - suggesting that under very cold protoplanetary disk conditions, collisonal growth may not favor ice over silicates.

Based on the above findings, we consider two different values of $U_{f}$ for H2O ice: when grains are near the ice line ($\sim 5$ AU), we adopt a fragmentation velocity of $U_{f}\approx 8.9$ m/s, which we refer to as the sticky H2O ice fragmentation speed consistent with both previous studies (e.g., Wada et al., 2009), and those of Musiolik & Wurm (2019). Well outside the snow line where temperatutes are low ($\sim 30$ AU), we consider both these stronger ice particles, and also a fragmentation velocity for ice having the same value as for silicate particles ($U_{f}\approx 1.4$ m/s, Zsom et al., 2010; Güttler et al., 2010), which we refer to as the cold H2O ice fragmentation speed.

7.4. Bouncing and drift barriers

Numerous laboratory and theoretical studies have found that particle growth can be influenced by bouncing at much smaller sizes than the size that collides so energetically as to fragment the particles (Zsom et al., 2010). Though Estrada et al. (2016) derive an expression for bouncing between similar-sized particles (their Eq. 59), this size is harder to specify rigorously because it likely also depends on material properties (Güttler et al., 2010; Zsom et al., 2010). Estrada et al. (2016) emphasize that the bouncing barrier is not impermeable but merely slows growth to a great degree. For this reason, it is not as instructive to include the bouncing barrier as a constraint as we do for fragmentation. A similar situation holds for the radial drift barrier alluded to earlier; it is hard to write quantitative values for this size limit (e.g., see Eqns. 60-61, Estrada et al., 2016). Our fragmentation barrier upper limits on particle size (Eq. 26) are thus approximate. Under some conditions, particles may never get that large due to a combination of the bouncing and radial drift barriers (Birnstiel et al., 2012; Estrada et al., 2016; Sengupta et al., 2019). Under other conditions, growth can proceed somewhat beyond the fragmentation barrier if a distribution of collision velocities and mass transfer in collisions are allowed for (e.g., Windmark et al., 2012; Drazkowska et al., 2013, 2014; Estrada et al., 2016). Bouncing, drift, mass transfer, and a probability distribution for collisional velocities are included in the more realistic growth model constraints discussed in Section 7.7.

7.5. Combined limits on particle size: $\tau_{s}$.

Global numerical simulations suggest that the series of barriers to growth discussed previously are quite effective at limiting $\tau_{s}$ of the particle size that dominates the mass (e.g., Birnstiel et al., 2012; Estrada et al., 2016; Sengupta et al., 2019). Recent simulations conducted for solid particle growth over a three-order-of-magnitude range of turbulent intensities ($10^{-5}<\alpha<10^{-2}$) relevant to the first 0.5 Ma where planetesimal formation is thought to begin (see Sec. 1) further indicate that the maximum achieved Stokes numbers of the mass dominant particles151515As in some of the models discussed in Sec. 6, the particle growth models of Estrada et al. (2016, 2020) exhibit a broad size distribution and do not employ particles of a single size. However, in general, most of the particle mass is near the fragmentation size when drift is not important, or near the largest particle size in the distribution when in the drift-dominated regime - either way, defining a single particle size. range from $\sim 0.001-0.01$ over a wide range of disk models with initial global metallicities of $Z\simeq 0.01$ (Estrada et al., 2016, 2020). Table LABEL:Stokes_numbers_alpha summarizes these Stokes numbers, the corresponding particle radii, and ambient disk conditions at our three representative radial locations within the protoplanetary disk. We note that in these simulations, the local solids abundance and the nebula gas density and temperature (and thus $\delta$, as well as the pressure gradient) are evolving with time which means that $Z$ can have a range of values over the disk’s radial extent. In particular, the snow line (and various other “evaporation fronts”) evolve with time leading to both local enhancements and depletions in the amounts of solids, especially in the planet forming region. The values given in Table LABEL:Stokes_numbers_alpha are selected at locations where the solids-to-gas ratio $Z(t)\simeq 0.01$ and which correspond to the “inner disk” and “snow line”; they may thus lie at semi-major axes that are slightly different from our nominal definitions of $R=3$ and $5$ AU, respectively. The “outer disk” location reliably corresponds to 30 AU (see caption, Table LABEL:Stokes_numbers_alpha).

Despite the variable metallicity, the particle size distributions have already reached a quasi-equilibrium state at the selected times and radial locations, changing only slowly with time and depending only weakly on the instantaneous value of $Z$ (see Fig. 19 of Estrada et al., 2016, and associated discussion). The physics that limits the mass dominant particle $\tau_{s}$ in different radial regions depends on the ambient nebula conditions. In the inner disk regions where bulk composition is ice-free, $c_{s}$ and $\Omega$ are large so turbulent relative velocities are fairly large. These regions tend to be in the fragmentation regime based on previous discussion, using Eq. (26) for the fragmentation equilibrium Stokes number. However, for $\alpha\gtrsim 10^{-3}$, the values for the fragmentation $\tau_{s}$ obtained from Eq. (26) using the parameters cited in Sec. 7 are consistently smaller than the values listed in Table LABEL:Stokes_numbers_alpha, because the actual growth models include growth beyond the fragmentation barrier (see Sec. 7.4). On the other hand, for $\alpha=10^{-4}-10^{-5}$, one finds Eq. (26) gives $\tau_{s}\simeq 0.01-0.1$ which is consistently larger than the corresponding values in Table LABEL:Stokes_numbers_alpha, despite the temperatures being closer to those cited in Sec. 7, and in fact these particles have already reached the fragmentation barrier. As noted in Table LABEL:Stokes_numbers_alpha, these mass-dominant particles are in the Stokes regime, which for these relatively small $\tau_{s}$ means the relative velocities between them tend to be higher for a given Stokes number compared to the Epstein regime. This discrepancy does not depend on which flow regime the particles are in, but rather because eddy-crossing effects start to become important in weak turbulence, particularly when $\delta^{2}/\alpha\gtrsim\tau_{s}^{-1}$. Under these circumstances the relative velocities between the particles are higher, even for similar-sized particles, which means the fragmentation $\tau_{s}$ will be smaller than what Eq. (26) would predict (Ormel & Cuzzi, 2007)161616 Ormel & Cuzzi (2007) found that the transition between so-called class I and class II eddies occurs at higher values with increasing $\delta^{2}/\alpha$. In this case, class II eddies, where the fluctation times are shorter than the particle stopping time, dominate even for small particles resulting in the kinetic energy per unit mass between colliding particles of a given $\tau_{s}$ is larger. This consequently progressively lowers the estimate from Eq. (26) with decreasing $\alpha$. For more discussion, see Ormel & Cuzzi (2007)..

In the colder, ice-rich outer disk a naive application of Eq. (26) would give $\tau_{s}\gtrsim 1$ for $\alpha\lesssim 10^{-4}$ for the sticky water-ice fragmentation case. However, the simple constraint of Eq. (26) assumes that the collision velocity is dominated by turbulent relative velocities, but large values of the pressure gradient (manifested by $\delta$) in the outer nebula can drive systematic drift- and headwind-related velocities that can significantly exceed those due to turbulence for these weak values of $\alpha$. Under these circumstances, bouncing plays an even more influential role by slowing growth, enhancing the importance of the drift barrier in precluding the fragmentation barrier from being reached. Simulations of the outer disk where particles drift faster than they can grow (Figs. 12 and 13) are characterized by Stokes numbers $\tau_{s}\lesssim 0.01$ that are nearly constant or modestly varying with radius, which suggests that the decrease in particle radii with distance simply mirrors the decrease in gas density (Birnstiel et al., 2012; Estrada et al., 2016, 2020). Even when the bouncing barrier was not considered, $\tau_{s}\lesssim 0.1$ in these regions. On the other hand, for large values of $\alpha\gtrsim 10^{-2}$ one finds that even in the outer disk particles can be in the fragmentation-dominated regime even if their Stokes numbers remain comparable to the lower $\alpha$ cases, because the fragmentation barrier occurs at much smaller sizes there (see Table LABEL:Stokes_numbers_alpha).

Overall then, Eq. (26) is a handy but crude approximation in general. Estrada et al. (2016) give a more in-depth discussion and derive estimates of the limits of $\tau_{s}$ as it pertains to their models. The variation in $\tau_{s}$ and particle size over different radial locations and with time listed in Table LABEL:Stokes_numbers_alpha is subtle but secondary to our focus here. After looking at many models, we believe the values presented are representative for the purpose of constraining regions of parameter space that may permit SI to form planetesimals over a wide range of conditions (see Sec. 7.7). Additional description of the colored-symbol models will be addressed in Estrada et al. (2020). Even with their uncertainties, the message of the realistic models is that plausible, self-consistent combinations of nebula turbulence and particle size typically lie in “Zone II”, the “moderately turbulent” regime above the $\epsilon=1$ line, where SI is only “incipient”. Moreover, any degree of particle porosity leads to even smaller $\tau_{s}$ (Estrada et al., 2020).

7.6. Turbulence Constraints

It is widely believed that the region where the first planetesimals were assembled was a “Dead Zone” extending from 1AU to $\leq$ 80 AU, so-called because the temperature and transparency of protoplanetary disk material to ionizing photons are too low to allow magnetically driven (MHD) turbulence to arise (Turner et al., 2014). These areas of protoplanetary disks are now more commonly referred to as Ohmic Zones (e.g. Lyra, & Umurhan, 2019), owing to the dominance of Ohmic dissipation that suppresses self-generated activity like the linear MRI process. There remains substantial debate over the true ionization state and consequent effective magnetic resistivity of disk material, depending on the abundance of grains of sufficiently small size (Okuzumi et al., 2012; Ormel & Okuzumi , 2013; Simon et al., 2018). Even if insufficiently ionized for MRI, the disk remains susceptible to nonideal MHD processes including wind launching mechanisms (Bai , 2016; Bai et al., 2016). Nonetheless, this debate leaves some uncertainty as to whether any regions may be susceptible to the MRI, which numerical experiments show to induce very strong turbulent intensities of $\alpha\sim 10^{-3}-10^{-2}$ (Balbus & Hawley, 1998; Armitage, 2011).

Therefore, for the purposes of this study we shall assume that any turbulence that arises in the near-midplane regions of protoplanetary disks, which are of greatest interest to planetesimal formation, stems from any one of the three purely hydrodynamical mechanisms noted in section 1 that have recently been discussed in the literature – for a review see Turner et al. (2014) or Lyra, & Umurhan (2019). The operation of the three mechanisms depends upon the thermal relaxation or cooling time scale $t_{{\rm th}}$ of the particle-gas mixture, which in turn depends on the disturbance lengthscales if these lengthscales are in the optically thick regime. Specifically, the [VSI/COV/ZVI] (section 1) is expected to operate in a disk when [$\Omega t_{{\rm th}}\ll 1,\Omega t_{{\rm th}}\approx 1,\Omega t_{{\rm th}}\gg 1$] (respectively). Exactly which regions of protoplanetary disks are most susceptible to which mechanism remains under discussion (Malygin et al., 2017; Umurhan et al., 2017; Barranco et al., 2018). However, because of their complementary instability criteria, it is plausible that the full extent of protoplanetary disks that is of interest to planetesimal formation will be turbulent, due to at least one of the three mechanisms listed. Recently published numerical studies of the three processes show that the turbulent intensities arising from these mechanisms lie somewhere in the range $10^{-5}<\alpha<5\times 10^{-4}$ (Lyra, & Umurhan, 2019). Of the various published studies of the VSI, Flock et al. (2017) predict the lowest level of $\alpha\sim 4\times 10^{-5}$, and we use this as our minimum adopted value. Our understanding of the nature of cold protoplanetary disk turbulence will surely continue to evolve.

Caveat: We have assumed the mechanism driving turbulence in disks is unaffected by the degree of particle loading. This is valid a posteriori in regions covered by all of the detailed particle growth models (figures 11 - 13) which imply $\epsilon<1$, and limited or negligible particle feedback. This assumption should be regarded with caution in parameter combinations where $\epsilon>1$, but these combinations (corresponding to Zone I) are already known to lead to SI in numerical experiments.171717For example, using 2D axisymmetric simulations of a single-fluid terminal velocity approximation model (Lin & Youdin, 2017) as well as physical arguments, Lin (2019) argues that a feedback loop is set up wherein particle-loading and subsequent settling drives buoyancy oscillations that weaken the VSI (Lin & Youdin, 2015) which, in turn, further enhances particle settling. VSI stabilization by this process was argued for values of disk metallicity $Z\geq 0.01$ and appears effective for $\tau_{s}>0.003$ in the published 2D simulations. However, whether or how much particle loading would damp turbulence in full 3D, more highly resolved simulations remains to be seen, as 3D energy cascade pathways and turbulent eddy motion couplings are different (Richard et al., 2016; Lyra, & Umurhan, 2019). Additionally, whether or not this tendency toward stabilization of the VSI persists in multi-particle component 3D disk simulations remains to be seen (e.g., see discussion in Benítez-Llambay et al., 2019).

7.7. Regions of parameter space that permit the SI to operate for $Z=0.01$

Given the physical constraints outlined in the previous six subsections, we can now assess in what parts of the $\alpha$–$\tau_{s}$ parameter space the SI is likely to operate. We examine this for three locations: $R=3$AU (“inner disk”), $R=5$AU (“snowline”) and $R=30$ AU (“outer disk”), in a $Z=0.01$ disk with constant $\delta=0.05$. By implication, we imagine the following survey to be predictive for an early phase of the solar nebula with uniform Z ($=0.01$) and well before wind loss substantially evaporates the disk’s gas (say, $<2.5$ Ma; Carrera et al., 2017). It is increasingly thought that planetesimal formation had to proceed in the inner nebula before the 0.5 Ma time marker (relative to CAIs), and that even Jupiter’s core may have formed in less than 1 Ma (e.g., Kruijer et al., 2017; Simon et al., 2018). For higher, perhaps local, values of $Z$, to which particle size is insensitive, the reader is referred to Figure 4 for comparison of how Zones I and II shift.

We approach this task by first comparing the SI’s predicted growth timescales with constraining timescales based on the various disk processes and particle growth barriers discussed in sections 7.1-7.6 for our nominal disk model with $\delta=0.05$. We then delineate regimes of ($\alpha,\tau_{s}$) space deemed implausible by the best current models of growth-by-sticking (Estrada et al., 2016, 2020, sections 7.4 and 7.5) over the first 0.5 Ma, which are typically characterized by larger $\delta\gtrsim 0.05-0.1$. We declare the SI to be “incipient”, or capable of leading to some degree of enhancement, in those regions of ($\alpha,\tau_{s}$) space which are both reasonable from the standpoint of particle growth, and where the SI’s growth timescales are shorter than the aforementioned constraining timescales.

Nominal fragmentation constraints: The results for our nominal models are shown in Figures 8-10, and those with the added constraints of particle growth models are shown in Figures 11-13. All these plots overlay various excluded regions on the growth timescale plot similar to Figure 3 for the three disk zones of interest. The timescales in Figures 11-13 however were recalculated using the higher $\delta$ values associated with the particle growth models. If a zone of incipient SI is allowed, it is indicated by a red-hatched triangle in these figures.

Figure 8 shows the situation for the inner disk where the temperature is too large to support water ice particles. We find that the SI is incipient in a relatively small triangular patch of parameter space centered about a value of $\alpha\sim 5\times 10^{-5}$ and $\tau_{s}\approx 0.01-0.02$. The region of accessible parameter space is bounded on the small side of $\tau_{s}$ lines showing conditions where the growing SI pattern drifts into the star faster than it can achieve an $e$-folding level of growth (see discussion in sec. 7.2). The region is also bounded on the higher side of $\tau_{s}$ by the constraint imposed by the silicate fragmentation line. Within this SI incipient region, the growth timescales are relatively long - thousands of orbit times.

The situation at the snow line can vary significantly depending on whether one adopts the sticky or cold H2O fragmentation condition, as can be seen in Figure 9. Taking the sticky-ice constraint opens up a much larger region for incipient growth ranging from $\tau_{s}\approx 0.01-1$, and for turbulent intensities as high as $\alpha\sim 5\times 10^{-4}$ that extends into Zone III, and even Zone I for the largest $\tau_{s}$ (see Fig. 3), and growth timescales are as little as 10’s of orbits. However, the sticky ice condition may only apply to regions that are close to the water ice evaporation temperature (Musiolik & Wurm, 2019). Adopting the cold water ice fragmentation condition returns the situation to the one seen in Fig. 8 with only a tiny incipient zone. Likewise at 30 AU, the region of SI incipient growth depends strongly on how sticky water ice is. Sticky water ice, as has been adopted in most growth models to date, leads to an even larger incipient region where turbulence can be as high as $\alpha\sim 10^{-3}$, and encompasses a range from $\tau_{s}\approx 0.008-4$ again extending as far as Zones III and I. On the other hand, the cold water ice fragmentation condition (Musiolik & Wurm, 2019) restricts the incipient region to a tiny region in Zone II, centered about $\alpha\sim 7\times 10^{-5}$ and bounded by $\tau_{s}\approx 0.008-0.1$.

Modeled particle growth constraints: We now examine how the constraints imposed by the evolutionary growth model results of Estrada et al. (2016, 2020) in the epoch of interest affect the regions of SI incipient growth. In Figure 11, we return to the inner disk, but now plot the maximum achieved Stokes numbers from Table LABEL:Stokes_numbers_alpha (colored stars) which correspond to the particle sizes that carry most of the mass (see Sec. 7.5). Note that Fig. 11 contains new SI mode growth rates, recalculated for $\delta=0.07$ which represents a characteristic value for the particle growth and drift models. The main effect of this is that the radial and pattern drift constraint lines have shifted downwards, with the latter further restricting the SI incipient region to a vanishingly small range even before particle constraints are folded in. However, when one considers the growth models – even for the minimum $\alpha$ case – this region becomes inaccessible. As was discussed in Sec. 7.5, not all the plotted points fall on the representative silicate fragmentation line which assumes that turbulence dominates their relative velocities, and that eddy-crossing effects are unimportant. The points defined by $\alpha\gtrsim 10^{-4}$ broadly lie above their respective fragmentation lines (as can be determined from the temperatures listed in Table LABEL:Stokes_numbers_alpha) because there has been significant growth beyond the fragmentation barrier (section 7.5) due to both “mass transfer” (e.g., Wurm et al., 2005; Windmark et al., 2012) and “lucky particle” growth in those models (e.g., Garaud et al., 2013; Drazkowska et al., 2014), whereas for $\alpha\lesssim 10^{-4}$ the mass-dominant particles are entering a low turbulence regime in which they are subject to eddy-crossing effects so that their fragmentation Stokes numbers are smaller than what Eq. (26) would predict. Thus, all the particle evolution models discussed here are in fact in the fragmentation regime. Overall then, we find that in the silicate-dominated region illustrated by Figure 11, all plausible, self-consistent combinations of turbulent intensity and particle size lie well within Zone II, and do not overlap any of the region in $\alpha$–$\tau_{s}$ parameter space in which SI incipient growth is permissible.

For both the nominal snowline and the outer disk, Figures 12 and 13 show that the situation is also strongly limited by realistic growth models, with respect to regions of incipient growth. As was the case in Fig. 11, Figure 12 corresponds to a value of $\delta=0.07$ which represents a characteristic value across all models listed in Table LABEL:Stokes_numbers_alpha for the snow line. Like before, the higher $\delta$ further restricts the incipient growth region, though not as dramatically as for the inner disk. In this case we see that the SI appears to be incipient only within a tiny triangular patch of parameter space, again around $\alpha\sim 5\times 10^{-5}$ and $\tau_{s}\approx 0.01$ indicated in the figure by the red-hatched triangular boundary. The constraining processes of the minimum level of turbulence, and the mode pattern drift timescales, bound the region from the bottom and left. However the evolutionary growth models between $\alpha\sim 10^{-5}-10^{-4}$ cut off access to the regions of $\alpha$–$\tau_{s}$ parameter space that would be permissible from either water ice fragmentation constraint. Figure 13 shows the case for the outer nebula at 30AU, in which $\delta=0.1$. Here the mode pattern drift line moves even further downward, increasing the lower bound of $\tau_{s}\simeq 0.04$ a factor of 5 more than for the $\delta=0.05$ case (Fig. 10) and shrinking the SI incipient zone. Meanwhile, particle Stokes numbers remain about $\tau_{s}\sim 0.01$, further to the left and eliminating the small incipient zone.

In Figure 12, almost all models with $\alpha\lesssim 4\times 10^{-3}$ are in the drift-dominated regime, whereas in Figure 13 the fragmentation limit has only been reached for $\alpha=10^{-2}$ (recall the fragmentation limit for the black and yellow symbols are more closely associated with the sticky H2O ice fragmentation line). In Fig. 12, the models for $\alpha=10^{-4}-10^{-5}$ are again in a regime where eddy-crossing is starting to become important, so their fragmentation $\tau_{s}$ will be smaller than what the sticky ice line predicts, but still larger than their current values. Here, only the model for $\alpha=10^{-5}$ is in the Stokes flow regime. The remaining models in Figs. 12 and 13 (orange stars) which are associated with the cold H2O ice model, have reached the fragmentation size (though significant growth beyond it has occurred in Fig. 12 at the snow line). The curious inflection that leads to peak values for $\tau_{s}$ for $\alpha\sim 10^{-3}$ are real and due in part to enhanced growth about the snow line, whereas the trend towards larger $\tau_{s}$ seen for $\alpha\sim 10^{-2}$ at 30 AU is due to both reaching the fragmentation barrier and the lower gas surface density as a result of the more rapid viscous evolution compared to the other models. Overall though, as was the case for the silicate-dominated inner disk region, all self-consistent combinations of turbulent intensity and particle size lie within Zone II, from the standpoint of SI.

The interplay between particle drift, bouncing and fragmentation, and the manner in which particle growth proceeds, especially around evaporation fronts like the snow line, highlight the complexity of modeling particle growth and gas evolution with time, and call for more in-depth analyses of this type in the context of the theory presented herein.

8. Summary and Conclusions

Conducting high resolution numerical simulations of the interaction of gas and particles to shed light upon the streaming instability is a computationally expensive undertaking. It should be of value to have some kind of theoretical guide – however approximate – to constrain the parameter range of its validity from the standpoint of planetesimal formation. The purpose of this study is to provide a theoretical framework to address the question of how and to what extent the streaming instability might be effective for planetesimal formation under globally turbulent disk conditions.

Our model extends and generalizes behavior initially studied by YG2005 by representing the effects of local turbulence by an $\alpha$-model, and makes additional predictions that are consistent with previously reported numerical studies of the presence or absence of the streaming instability in turbulent disk simulations which include the following: Johansen et al. (2007); Balsara et al. (2009); Yang et al. (2017); Li et al. (2018); Yang et al. (2018); Gerbig et al. (2020).181818 For the one set of simulations we examined where the correspondence was weakest (i.e., Li et al., 2018) we conjecture that those simulations were not run with sufficient resolution to see the short wavelengths predicted by our model. The $\alpha$ model representation of turbulence – that characterizes its effect in the form of an enhanced isotropic turbulent viscosity and diffusion – acts locally both to stir particles and to exchange momentum. Underpinning its use here is the assumption that the processes leading to turbulence, especially in protoplanetary disk Ohmic zones, do so independent of the presence of particles with the realistic sizes and abundances treated here.

We have examined the normal mode response of the streaming instability as a function of the disk turbulence parameter $\alpha$ and particle Stokes numbers $\tau_{s}$ by identifying the wavelengths, growth times, and pattern speeds of the fastest growing modes. For given values of $\alpha$ and $\tau_{s}$, the particle to gas mass density ratio ($\epsilon=\rho_{p}/\rho_{g}$) is calculated using a turbulent dilution model (eq. 21), to represent the balance between the gravitational settling of particles toward the disk midplane and the vertical diffusion of the same particles due to turbulence (Dubrulle et al., 1995; Youdin & Lithwick, 2007; Estrada et al., 2016). We hope the theoretical framework proposed and examined in this study is useful in similar future studies.

Simple turbulence models like ours involve taking higher order moments of the equations of motion which are then truncated with some a closure relationship. In general, such resulting model equations do not conserve momentum (see extended discussion in Davidson , 2004). Tominaga et al. (2019) highlighted the lack of angular momentum conservation in a 2D non-axisymmetric disk setting with a turbulence model similar to ours, and proposed a solution to the problem. Together with the authors of Chen & Lin (2020), we have conducted a proper re-analysis of the matter in the framework of these equations (not shown here) and find that this non-conservation has negligible effect on the growth rates determined for all values of $\alpha$ pertinent for realistic protoplanetary disks. A future follow-up study explicitly detailing this result is warranted.

The study conducted here has revealed several interesting trends for a nominal (minimum mass) solar nebular model with global solids-to-gas mass ratio of $Z=0.01$ and disk opening angle $\delta=H/r=0.05$ (results are also given for values of $Z$ as large as 0.08). While the specific values summarized below primarily pertain to disk models with $Z=0.01$, the conclusions might be extended to higher $Z$ by comparison with figure 4 (although further analysis should be done to verify this assertion as well):

As turbulent intensity increases, the wavelengths of the maximally growing SI modes increase. while the growth rates of the maximally growing SI modes diminish. 2. 2.

The combination of ($\tau_{s},\alpha$) that leads to initial $\epsilon=1$ according to the turbulent dilution model traces a critical curve with important implications for SI behavior.This curve terminates at a critical point $(\tau_{c},\alpha_{c})$ corresponding to values of $\epsilon$ that monotonically increase with $Z$ for given $\delta$. For cosmic abundance $Z=0.01$, this critical point occurs at $\tau_{c}=0.45$ and $\alpha_{c}\approx 3.7\times 10^{-5}$, corresponding to $\epsilon=\epsilon_{c}=1.1$. Selected critical point values for other $Z$ are summarized in Table LABEL:Critical_Values_alpha. 3. 3.

For values of $\tau_{s}<\tau_{c}$, and for parameter combinations of $\alpha$ and $\tau_{s}$ that lead to $\epsilon=1$ according to the turbulent dilution model (eq. 21) - that is, along the critical $\epsilon=1$ curve - the least stable SI mode neither grows nor decays and the growth timescale is effectively infinite. 4. 4.

Provided $\tau_{s}<\tau_{c}$ we identify two regimes which straddle the above-mentioned critical line as being either “laminar/unstable” (Zone I) or “turbulent/saturated” (Zone II; Figure 3). 5. 5.

The spatial structure of the fastest growing mode in Zone II (the turbulent/saturated regime) typically corresponds to vertically oriented sheets with radial scale of about a pressure scale height $H$. In practice, the sheet’s vertical extent should follow the particle scale height and this appears to be consistent with the simulation results reported in Balsara et al. (2009) (see their Fig. 7). The mode structure in Zone I (the laminar/unstable regime) exhibits narrow, azimuthally oriented tubes with lengthscales much less than $H$, consistent with previous predictions and simulations (e.g., see the recent results of Carrera et al., 2015; Yang et al., 2017; Li et al., 2018). 6. 6.

The theory developed here appears to reasonably predict the onset or absence of the SI in several recently published simulations, We have plotted these correspondences in our $\alpha-$St growth timescale plots in Figures 3 and 4. For laminar disk models in which particles settling to the midplane generate their own midplane turbulence from either the from Kelvin-Helmholtz overturn and/or the SI itself, we find that our theory (i) reasonably predicts the radial spacing of emerging filaments when SI is active especially during its early onset phases well before filament-filament merging occurs and (ii) predicts when the SI does not appear. We have come to this conclusion after careful analysis of the three recent studies by Yang et al. (2017); Li et al. (2018); Gerbig et al. (2020).

In simulations where turbulence is externally driven throughout the disk model by the MRI like Johansen et al. (2007); Balsara et al. (2009); Yang et al. (2018), the theory does a reasonable job at predicting whether or not SI should or should not be present and, in some cases, predicts the radial wavelength structure and growth timescale of filaments like in Johansen et al. (2007). In particular, the theory predicts the absence of the SI in the simulations of Yang et al. (2018) for their ideal MHD model in which the MRI is rampant throughout the computational domain. On the other hand, Yang et al. (2018) also present results of Dead Zone disk models driven by waves coming from sandwiching MRI active layers. In the absence of backreaction particles accumulate along filaments likely because the waveforcing emplaces azimuthally oriented pressure fluctuations with coherent pressure maxima toward which particles naturally drift. When backreaction is turned on, the particle densities within these filaments become enhanced. Our theory does not explicitly handle SI physics in the presence of coherent structures: it implicitly assumes the scales of the system take place well inside the inertial range of the turbulent forcing where coherence is lost. This matter must be addressed in future work. 7. 7.

We believe that realistic protoplanetary disk conditions relevant to the early solar system ($\lesssim 0.5$ Ma) in which $Z=0.01$ and where typically $\delta\gtrsim 0.05-0.1$ – taking into account various barriers to particle growth, age constraints, and the disk’s likely degree of turbulence as quantified by $\alpha$ – restrict the SI to Zone II where it is “incipient”, or allowed, for a narrow range of self-consistent disk and particle properties. For example, we find at a disk location nominally representative of the ice line of Jupiter’s early core ($r=5$AU) that the incipient range of parameter space as indicated by the red-hatched triangle in Fig. 12 is restricted about a turbulent intensity $\alpha\sim 4\times 10^{-5}$, and Stokes numbers in a narrow range between $\sim 0.01<\tau_{s}<0.02$. We treat this vanishingly small allowable region as only suggestive - given uncertainties in all the constraints, it may be larger. On the other hand, no such permissible regions are found in the inner and outer disks cases for the given models. Thus achieving growth sufficient to breach the incipient regions in the inner and outer disk appears even much more challenging. 8. 8.

In Zone I, SI is robust and plausibly proceeds to planetesimal formation, as routinely observed in numerical simulations. In Zone II, we identify a new kind of behavior in which SI is only “incipient”, the growth timescales are very long ($10^{2}-10^{4}$ local orbit times) and the growth of particle density is limited. 9. 9.

It should be kept in mind that while the analysis in Section 7.5 suggest that it may be difficult for the SI to trigger sufficient particle enhancements in turbulent disks with spatially uniform values of $Z=0.01$, that this prediction is also likely to be applicable only for the most earliest phases of a model solar disk in which there is negligible particle enrichment. As a disk evolves, the disk’s gas content will evolve both globally due to wholesale loss via winds and locally variations induced by transport. Over time, then, regions will emerge with substantially enhanced values of $Z$. The regions which are “SI incipient” can, in principle, slip from being in Zone II to Zone I, wherein the growth is expected to be more rapid even in the face of turbulence. A counterargument here might be that decreasing the gas density through winds might not necessarily lead to higher values of St nor enhanced local values of $Z$ because the efficiency of radial drift can efficiently deplete a local region of its particles while particles might encounter their fragmentation barriers at smaller sizes (Birnstiel et al., 2012; Estrada et al., 2016; Carrera et al., 2017). We therefore expect the situation to be very different for conditions corresponding to the latter stages of the disk’s evolution beyond its thick gas phase, and we caution against applying our prediction of limited to no growth for the whole of the disk’s lifetime. Further follow up analysis are therefore necessary and warranted.

Other Future Work: There is much to understand about the implications of this new physics. Included in any list of future work should be (a) a better physical understanding of the critical line and critical point corresponding to the special condition $\epsilon=1$ and, especially, why this special combination of parameters leads to exactly marginal modes, (b) the role of a particle size distribution; (c) the possible role of particle loading on damping turbulence; (d) a self-consistent analytical model of particle layers and self-generated turbulence in the limit of vanishingly small global turbulence, (e) the behavior when the predicted vertical wavelength exceeds the thickness of the particle layer, and of course (f) the regime where purely hydrodynamical turbulence is operative, and its intensity. It is worthwhile to re-do this analysis within the framework of (i) a turbulent MHD model, and/or (ii) a model with multiple particle sizes as in Krapp et al. (2019), (iii) and to expand the theory to account for the presence of coherent flow structures, e.g., like those present in Yang et al. (2018).

Finally, understanding the physical mechanism of the turbulent SI remains a priority. While the mathematical description of the onset of the inviscid SI in terms of resonant drag energy exchange between gas waves and particles is satisfying (i.e., the RDI, Squire & Hopkins, 2018a), it might need modification to explain the onset of instability in this kind of turbulent model setting. The mechanistic explanation of the inviscid SI presented by Jacquet et al. (2011) offers a framework upon which an explanation for the viscous case might be built.

Acknowledgements

We thank our two reviewers who provided very thorough readings and suggestions that improved both the rigor of our analysis and the quality of this manuscript. We are indebted to Allan Treiman for organizing “Accretion: Building New Worlds Conference” held at the Lunar and Planetary Institute (August 15-18, 2017), where the motivation for this study first germinated via in-depth conversations with A. Carballido, J. Simon and D. Carrera. We thank Phil Armitage for organizing the “Theoretical and Computational Challenges in Planet Formation Workshop” at the Flatiron Institute (May 20-22, 2019) where mature results of this study were presented. We are grateful for extended conversations with K. Shariff and D. Sengupta and others at NASA Ames Research Center’s Origins Group as well as input and perspectives shared with us especially by M-K. Lin and K. Chen who developed very similar results independently of – but in parallel with – us, as well as X. Bai, H. Klahr, W. Lyra, J. Drazkowska. O.M. Umurhan acknowledges support from NASA Astrophysics Theory Program grant NNX17AK59G and The New Horizons Kuiper Belt Extended Mission. P.R. Estrada and J. N. Cuzzi acknowledge NASA Emerging Worlds grant NNX17AL60A for providing resources and support for this project. Finally, we also acknowledge NASA’s “Internal Scientist Funding Program” grant to Ames Research Center on Origins of Planetary Systems as well as NASA’s Senior NPP program.

Appendix A A. Linearized matrix equation

The linearized equations of motion are non-dimensionalized on $\Omega$ and $H$. Furthermore, in order to insure incompressibility, the radial and vertical components of the perturbation gas velocity (respectively, $u_{g}^{\prime},w_{g}^{\prime}$) is expressed in terms of the azimuthal perturbation streamfunction $\psi^{\prime}$, which is related to the gas vorticity via, $\omega_{g}^{\prime}=-K^{2}\psi^{\prime}$) to insure incompressibility. The result may be cast into the following matrix form

[TABLE]

where $\mathtt{M}$ is given by

[TABLE]

where $K^{2}\equiv k_{x}^{2}+k_{z}^{2}$ and for notational convenience we define

[TABLE]

We could approach solutions to this by solving for $\omega$ from the sixth order dispersion relation arising from

[TABLE]

where $\mathtt{I}$ is the $6\times 6$ identity matrix. However, the expression resulting from the above operation is extremely unwieldy and offers no insight toward the mechanisms operating in the instability. Instead, we choose to solve this directly by determining the eigenvalues of $\mathtt{M}$ using standard numerical techniques found in Matlab. We arrange the solutions in descending order of $\omega_{i}$. The solutions plotted throughout the text are the least stable mode. In some instances there are two unstable modes but the second mode is usually dwarfed in magnitude by the first mode. A detailed examination of the second unstable mode and its interpretation remains to be determined, but may be of interest. We are reminded that in the inviscid theory the characteristic lengthscales of the SI are $H\delta$ while lengthscales in the $\alpha$-disk model are measured on $H$. Hence, a disparity of lengthscales between the large-scale global turbulence and the SI appears with the ratio $\alpha^{1/2}/\delta$.

Appendix B B. On the Dead Zone models of Yang et al. (2018)

By contrast, the DZ model results (especially $Z=0.04,0.08$) exhibit particle accumulation along azimuthally elongated filaments which grow tighter, become increasingly coherent and achieve high densities, plausibly due to SI. However, we think that the physics in this case is more complicated than our turbulence model can address. Yang et al. (2018) depict the particle behavior in DZ models with and without backreaction (see bottom two rows of Figure 10 of Yang et al., 2018). Even in the models without backreaction (where SI is inactive), particles can be seen to accumulate into azimuthally elongated filamentary structures almost from the very beginning of the model runs. In these models the value of $\epsilon$ in the densest regions appears to exceed unity in many places. In comparison, all DZ models with backreaction show the same early accumulation as in the non-backreacting case. However, with backreaction these density enhancements continue to evolve because the SI drives further density increases.

Yang et al. (2018) characterize the vertical distribution of particles in terms of an effective particle diffusion parametrized by an analagous $\alpha$-parameter called $\alpha_{g,z}$. This diffusion parameter achieves high values ($\alpha_{g,z}\sim 4\times 10^{-3}$) similar to the values of $\alpha$ in the iMHD case, even though the value of $\alpha$ in the midplane regions of the DZ model is quite low ($\sim 2\times 10^{-4}$). Nonetheless, particles appear to concentrate into filaments which become reinforced and further enhanced as the SI takes root, this is especially apparent in the simulations involving $Z=0.02,0.04,0.08$. Such filaments are likely pressure extrema with concomitant zonal flows are known to be characteristic of MRI driven turbulence at the large scales of a simulation (Johansen et al., 2009; Simon et al., 2012).

These scales are, at best, at the top of the inertial range of the turbulent flow and are not contained inside the turbulence’s inertial range. To what extent the midplane layers in their DZ models may or may not be characterized as a uniform turbulent flow (see section 4.1 of Yang et al., 2018) and what effect external wave forcing has on inducing a priori high density azimuthally aligned filaments remains to be fully examined. By examining the development of filaments in models with backreaction against those without backreaction , it can be seen in the latter that filamentary particle enhancements emerge as a matter of course in response to the large scale wave-forcing coming from the magnetically active layers (see also Figure 8 of Yang et al., 2018).191919Similar structuring can be seen in the simulations of Johansen et al. (2007). With backreaction turned on, these particle enhancements readily continue their condensation into narrow coherent structures. It is difficult to disentangle the causal sequence of events: are the observed structures growing due solely attributed and intrinsic to the SI or does the pre-existence of overdense filaments aid in triggering the SI in these cases? It would seem that the large scale azimuthally elongated structural forcing imprinted by the waves emanating from the overlying turbulent layers predisposes the midplane particles into undergoing the SI because of the a-priori axisymmetric clumping caused by the radial pressure maxima induced by the waves. The coherent large scale structure of the non-uniform wave-forcing is a physical complexity not reflected in the simple turbulence model of our theory, nor of that of Chen & Lin (2020). We therefore view the physical structures observed in the midplane regions of their DZ model as outside the scope of the theoretical construct discussed here and could be addressed in future work. Finally, the relevance of models with such deep MRI-active layers is itself open to debate (e.g., Bai & Stone, 2010).