approaches 1. Overplotted is the approximate value $\mu=\mu_{{\rm max}}$
corresponding to the fastest growing mode for ${\rm Ri}<1$. Note that this
corresponds to a distinct value $k_{{\rm max}}H_{p}$ that is implicitly a
function of vertical node number $m$ and other properties of the settled
particle layer like $\varepsilon_{0}$ and ${\rm Ri}$, see Eq. (58).

We therefore present here a motivated single fluid model (see details in
Appendix C) providing an analytically tractable demonstration of the SymI for
midplane disk layers. The simplified model can be thought of as the isothermal
(i.e., $\gamma\rightarrow 1$) and perfectly coupled terminal velocity limiting
form of the single-fluid model for dusty-gas mixtures (e.g., Laibe & Price,
2014; Lin & Youdin, 2017). This limiting form is the extreme limit of the
recast two-fluid equation framework developed in Laibe & Price (2014) where
the equations of motion are written without approximation instead in terms of
center of mass velocities and relative velocities of the two fluid species.
This “single fluid” limiting form emerges from this recast set of equations
after assuming both the terminal velocity approximation together with taking
the asymptotic limit of nearly zero particle stopping times.

To isolate the SymI effect from KH-roll-up, we assume there is only a purely
azimuthal mean flow state, which is the sum of the Keplerian profile plus a
departure $\tilde{v}_{{}_{0}}$ with a parabolic vertical variation given by

$\tilde{v}_{{}_{0}}=\delta
v_{{}_{00}}\left(1-\frac{1}{2}\frac{z^{2}}{H_{s}^{2}}\right),$ (45)

where $\delta v_{{}_{00}}$ characterises the mean azimuthal velocity drop
across 2-3 particle scale heights. We adopt this as a reasonable facsimile of
the center of mass azimuthal flow of the simulations, especially during the
bounce and shear-turbulent phase phase of the runs, also see Fig. 15 for more
details. We justify neglecting the vertical variation of the radial flow based
on its relatively small amplitude in comparison to the azimuthal flow – e.g.,
based on how ${\rm Ri}_{{r}}\gg{\rm Ri}_{{\phi}}$ for all of the shear phase
St = 0.04 simulations depicted in Fig. 14.

The parabolic form for $v$ adopted in Eq. (45) is a good fit for the mean
profiles that emerge in 3D simulations at various early-to-mid stages of
turbulent shear phase development (discussed at the end of Sec. 3.5) and
embodied in the assumed approximate Gaussian form expressed in Eq. (28). We
assume a mean density $\rho_{0}$ whose vertical variation is given by

$\partial_{z}\ln\rho_{0}=-\frac{z}{H_{p}^{2}},$ (46)

and, following both Chiang (2008) and Gerbig et al. (2020), subject to a
reduced gravity term acting toward the midplane, and given in magnitude by

$g=g_{{}_{\rm red}}=\frac{\varepsilon_{0}}{1+\varepsilon_{0}}\Omega_{0}^{2}z,$
(47)

where $\varepsilon_{0}$ is a constant (unlike its intepretation and usage in
the previous section), nominally representing the dust-to-gas ratio of the
midplane, in the same sense as used in Sec. 5.1. In this framework given Eqs.
45-46 together with no radial velocity, we find according to Eq. (24)’s
definition of the azimuthal Richardson number that

${\rm Ri}={\rm
Ri}_{{\phi}}=\frac{\varepsilon_{0}}{1+\varepsilon_{0}}\frac{\Omega_{0}^{2}H_{s}^{4}}{\delta
v_{{}_{00}}^{2}H_{p}^{2}}.$ (48)

We note that henceforth we consider only values of ${\rm Ri}>0$, i.e., stably
stratified flows. In order to connect to results discussed in section 3.5,
${\rm Ri}_{{\phi}}$ should be considered as being parallel to and in the same
spirit of the definition for ${\rm Ri}_{{\phi,0}}$ expressed in Eq. (29).

We assume the dynamics are incompressible, which means the radial and vertical
velocity perturbations may be written in terms of a single streamfunction
$\psi$. However, vertical density variations are advected by the perturbation
flow and give rise to buoyancy effects in the dynamics. Except for being set
in a shearing box framework, there is no other physics included including drag
exchange and fluid viscosity.

All perturbations are axisymmetric, and we assume normal mode solutions that
have the functional form for the streamfunction $\psi^{\prime}$, i.e.,
$=\hat{\psi}(z)\exp{\Big{(}-i\omega_{0}t+ikx\Big{)}}+{\rm c.c.}$, where
$\omega_{0}$ is the normal mode frequency, $k$ is the horizontal wavelength of
the disturbance, and $\hat{\psi}$ is the vertical eigenmode structure
function.

In Appendix C we analyze normal mode solutions of the above described
idealized model, which are expressed in terms of parabolic cylinder functions
${\cal D}_{m}(z/\tilde{\beta})$ (Abramowitz & Stegun, 1972). Here
$\tilde{\beta}$ is the lengthscale whose real part quantifies the Gaussian
decay lengthscale characterizing ${\cal D}_{m}(z/\tilde{\beta})$ – see Eq.
(C16). The general solution for the frequency $\omega_{0}$ is given in (C22)
which is a function of three parameters: ${\rm Ri}$ as defined in Eq. (48),
the non-negative integers $m$ signifying the number of vertical nodes in the
disturbances, and finally $\mu$, in which

$\mu\equiv\frac{(2m+1)^{2}}{2k^{2}H_{p}^{2}}\cdot\frac{\varepsilon_{0}}{1+\varepsilon_{0}},$
(49)

that characterizes a single parameter family of solutions containing $k$ and
$\varepsilon_{0}$. The general solution for the normal mode response is found
in Eq. (C22) reproduced here in slightly rewritten form:

$\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}=\sqrt{\mu\left(2\frac{{\rm Ri}-1}{{\rm
Ri}}+\mu-\frac{2i}{\nu}\sqrt{\frac{2\mu}{{\rm
Ri}}}\right)}-\mu+\frac{i}{\nu}\sqrt{\frac{2\mu}{{\rm Ri}}},$ (50)

where we have introduced $\nu\equiv 2m+1$ for notational convenience. Further
analysis found in Appendix C proves that normal modes come in growing/decaying
pairs for all finite values of Ri $\neq{\rm Ri}_{{\rm c}}(m)$, in which

${\rm Ri}_{{\rm c}}(m)\equiv 1-\frac{1}{\nu^{2}}=\frac{4m(m+1)}{(2m+1)^{2}}.$
(51)

Normal modes are marginal for the countably infinite set of ${\rm Ri}={\rm
Ri}_{c}(m)$ values. Inspection of Figure 24 indeed shows that the growth rates
become dramatically large once ${\rm Ri}$ passes below ${\rm Ri}_{c}$. Given
${\rm Ri}_{c}$’s dependence on $m$, nominally speaking we expect strong
instability when ${\rm Ri}$ passes under $1$, with instability appearing first
as $m\rightarrow\infty$.

It is instructive to showcase certain limiting cases for the above solutions
which we do in the following subsections.

####  5.2.3 $\mu\gg 1$ limiting form, a ${\rm Ri}<1$ criterion

We first start by considering the fate of modes in the limit where $\mu\gg 1$
which corresponds to large values of $m$ and/or small values of $kH_{p}$.
According to Eq. (C27), we find

$\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}\approx\left(1-\frac{{\rm Ri}_{c}}{{\rm
Ri}}\right)\left(1+i\frac{1}{\nu}\sqrt{\frac{2}{\mu{\rm Ri}}}\right)+{\cal
O}\left({\frac{1}{\mu}}\right),$ (52)

which means that up to ${\cal O}\left({1/\sqrt{\mu}}\right)$ $\omega_{0}$ is

$\displaystyle\frac{\omega_{0}}{\Omega_{0}}$ $\displaystyle=$
$\displaystyle\pm\sqrt{1-\frac{{\rm Ri}_{c}}{{\rm
Ri}}}\cdot\left(1+i\frac{1}{2\nu}\sqrt{\frac{2}{\mu{\rm Ri}}}\right)+{\cal
O}\left({\frac{1}{\mu}}\right).$ (53)

There is a dramatic change in the character of the modes when ${\rm Ri}-{\rm
Ri}_{\rm c}$ crosses zero. When ${\rm Ri}>{\rm Ri}_{c}$ the oscillating
disturbances, while being unstable, have a growth rate that is weak, and
proportional to $1/\sqrt{\mu}$. However once ${\rm Ri}<{\rm Ri}_{c}$, then the
modes show order 1 growth. This says that the character of the disturbances
undergo a stark transition at ${\rm Ri}\approx{\rm Ri}_{c}$. We also observe
that the modes with the highest values of $m$ are the most unstable when ${\rm
Ri}<{\rm Ri}_{c}$. For a layer undergoing collapse, where the value of ${\rm
Ri}$ is steadily going down, the modes with the highest values of $m$ begin to
go unstable once ${\rm Ri}<1$ since ${\rm Ri}_{c}(m\rightarrow\infty)=1$.
Outside of $m=0$, for which the theory predicts stability for all ${\rm Ri}$,
the critical values for the first few vertical nodes are ${\rm
Ri}_{c}(m=1,2,3,4,5,\cdots)=8/9,24/25,48/49,80/81,120/121$, and so on.
Therefore, for all practical reasons we view SymI to be relevant once the
particle layer satisfies ${\rm Ri}<1$ throughout most of its bulk. For the
collapsing layer solutions considered in the previous subsection, this would
appear to start taking root once the midplane particle layer achieves values
of $\varepsilon_{0}\approx 0.5$ (for example, see both panels of Fig. 23).

Figure 25: Following Eq. (50), this figure depicts predicted growth rates of
the SymI shown as a function of radial wavenumber for several vertical
parameters $m$ for a given set of values Ri, $\varepsilon_{0}$ and $H_{p}$.
Following the analysis of section 5.2.4, the asymptotic growth rate
$\sigma_{{\rm max}}(m\rightarrow\infty)$ is also shown.

####  5.2.4 $0<\mu\ll 1$ limiting form

We consider values of $\mu$ that are small, which for all practical purposes
corresponds to very large values of $kH_{p}$ and/or small values of
$\varepsilon_{0}$; however we consider the former case to be of practical use.
From Eq. (C26) we find to leading order that

$\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}\approx i\sqrt{\frac{2\mu}{{\rm
Ri}}}\left[\frac{1}{\nu}-\sqrt{1-{\rm Ri}}\right].$ (54)

The critical value Ri = Ric is easily recovered by equating to zero the
expression found within the brackets. The growth rates, while non-zero all
throughout ${\rm Ri}\neq{\rm Ri}_{c}$, are fairly muted for ${\rm Ri}>{\rm
Ri}_{c}$ but shoot up with a dramatic change in character once ${\rm Ri}<{\rm
Ri}_{c}$. The oscillation frequency also shows a strong shift around ${\rm
Ri}={\rm Ri}_{c}$.

####  5.2.5 $m\gg 1$ and $kH_{p}={\cal O}\left({m}\right)$: A maximum growth
rate

One can consider the asymptotic limiting form of Eq. (50) in the limit where
$m\gg 1$ together with maintaining $kH_{p}={\cal O}\left({m}\right)$. This
tandem limit keeps $\mu$ an ${\cal O}\left({1}\right)$ quantity according to
its definition in Eq. (49). Thus to leading order it follows that

$\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}=\sqrt{\mu^{2}-2\mu\left(\frac{1}{{\rm
Ri}}-1\right)}-\mu+{\cal O}\left({\frac{1}{m}}\right).$ (55)

An elementary analysis shows that Im$(\omega_{0})$ achieves a maximum value
when the expression underneath the squareroot operator equals zero. Therefore
we define a maximally growing value of $\mu$ to be

$\mu_{{\rm max}}=2\left(\frac{1}{{\rm Ri}}-1\right),$ (56)

whereupon it follows that – provided, of course, Ri $<1$ – the fastest growing
mode has a growth rate, $\sigma_{{\rm max}}\equiv$ Im$(\omega_{0})$,
approximately given by

$\frac{\sigma_{{\rm max}}}{\Omega_{0}}=\sqrt{\mu_{{\rm
max}}}=\sqrt{2\left(\frac{1}{{\rm Ri}}-1\right)}.$ (57)

We show overlain on Fig. 24 the approximate relationship relationship
$\mu_{{\rm max}}({\rm Ri})$. Finally, restoring the relationship between
$kH_{p}$ and $\mu$, we find a corresponding fastest growing horizontal mode,
$k_{{\rm max}}$ to be

$k_{{\rm
max}}H_{p}\approx\left(m+\frac{1}{2}\right)\sqrt{\frac{\varepsilon_{0}}{1+\varepsilon_{0}}\cdot\frac{{\rm
Ri}}{1-{\rm Ri}}}.$ (58)

Inspection of Figs. 24-25 indicates that Eq. (58) does a fairly good job at
predicting the fastest growing mode even for ${\cal O}\left({1}\right)$ values
of $m$. The asymptotic growth rate expression from Eq. (57) works well as an
upper bound. Nevertheless its utility is self-evident as it overpredicts the
low order $m$ growth rates by about 25%, at most.

###  5.3 Transition to secondary state via SI?

It is worthwhile to ask for the two St numbers simulated, if the transition
from the shear driven turbulence phase into the nonlinear pattern state is
consistent with the SI. We consider this by estimating the growth rate of the
particle scale height during this transition phase by approximating the domain
averaged time dependent scale height $\bar{H}_{p}(t)$ as

$\bar{H}_{p}(t)=H_{p,{{\alpha}}}+\delta h\exp\left(\frac{t}{t_{g}}\right)$
(59)

being a sum of a mean turbulently set particle scale height,
$H_{p,{{\alpha}}}$ that also corresponds to an estimated value of
$\alpha\approx{\rm{St}}\left(H_{p,{{\alpha}}}/H\right)^{2}$ (Dubrulle et al.,
1995), and an exponentially growing perturbation piece characterized by an
amplitude $\delta h$ and a growth time scale $t_{g}$. The aim here is to
compare the simulation derived values of $t_{g}$ with the prediction made in
the theory of turbulent SI developed in Chen & Lin (2020) and Umurhan et al.
(2020).

The left and middle panels of Fig. 26 shows a log-linear plot of the time
series for $\bar{H}_{p}$ for four simulations. Various dashed lines show our
estimated fitted values for $H_{p,{{\alpha}}}$ and corresponding estimated
slopes (dotted lines) drawn over the secondary growth phase where the
approximate values of $\bar{H}_{p}(t)\approx\tilde{H}_{p}(t)$ may be read off
of these lines for any input time $t_{i}$ falling nominally in the growth rate
time range. Based on this understanding of $\tilde{H}_{p}(t)$ it follows

$\ln\left[\tilde{H}_{p}(t_{i})-H_{p,{{\alpha}}}\right]\approx\ln\delta+\frac{2\pi
t_{i}}{t_{g}},$ (60)

which, after selecting two times from the fitted lines, allows one to estimate
$t_{g}$, i.e.,

$t_{g}\approx(t_{2}-t_{1})\bigg{/}\ln\left[\frac{\tilde{H}_{p}(t_{2})-H_{p,{{\alpha}}}}{\tilde{H}_{p}(t_{1})-H_{p,{{\alpha}}}}\right].$
(61)

Taken from the theory developed in Umurhan et al. (2020) the right panel of
Fig. 26 shows the predicted growth timescales of the SI for $Z=0.01$ and
$\beta=0.05$ as a function of $\alpha$ and St as maximized over wavenumber.
Note that the growth timescales in that study are quoted in units of local
orbit times, i.e., $P_{{\rm orb}}=2\pi\Omega^{-1}$. Thus, we read off the
quoted values from that graph and multiply it by $2\pi$ and call the result
$t_{g,{{\rm th}}}$ in order to compare the predicted growth rates against the
values of $t_{g}$ derived from the simulations presented here, the latter of
which are in units of $\Omega^{-1}$. The results of this exercise are
summarized in Table 4.

Table 4: Simulations and Predicted SI growth rates maximized over wavenumber

$\displaystyle\begin{array}[]{c}{\rm Simulations}\\\ {\rm Identifier}\end{array}$ | St | $\displaystyle\frac{H_{p,\alpha}}{H}$ | $\alpha\times 10^{6}$ | $\displaystyle{t_{g}\Omega}$ | $\displaystyle{t_{g,{{\rm th}}}}\Omega$
---|---|---|---|---|---
B3D-04L | 0.04 | 0.0095 | 3.57 | 2306 | 15800
B3D-04M | 0.04 | 0.0088 | 3.10 | 59.7 | 2513
B3D-2L | 0.2 | 0.00625 | 7.83 | 10.5 | 16.5
B3D-2M | 0.2 | 0.0056 | 6.30 | 9.11 | 11.2

†Based on theoretical predictions of Umurhan et al. (2020).

Examining the right panel of Fig. 26 depicting the results of the low to
medium resolution runs with St = 0.2, we find that the level set during the
onset of the shear turbulent phase has a lower value of $H_{p,{{\alpha}}}$ (by
about 15%) for the higher resolution of the two, a trend observed in other
simulations of the SI (e.g., Yang et al., 2017; Li et al., 2018, and others).
The lower resolution run therefore operates at a higher effective $\alpha$.
Theoretically it would imply a longer growth rate for the lower resolution run
and this trend is borne out by the prediction for turbulent SI. However, a
comparison of ${t_{g,{{\rm th}}}}$ against the measured values $t_{g}$ values
found in Table 4 shows that the theory appears to systematically underpredict
growth rates.

The middle panel of Fig. 26 shows the corresponding runs for St = 0.04. The
same trends, i.e., vis-à-vis shear turbulent phase, are manifest here as well:
the lower resolution run corresponds to higher turbulence levels compared to
the higher resolution, with concomitant faster growth rates of the latter
compared to the former. However, the growth rates are vastly shorter than
those predicted for the turbulent SI according to Umurhan et al. (2020).

At this stage it is difficult to conclude as to what may be driving this
secondary growth, and we can at best only speculate. It might be possible for
these runs that the classical SI is acting in conjunction with a secondary
effect that derives from the turbulent vertical shear profile (e.g., the VSSI,
Lin, 2021), but that the effect of the classical SI is stronger than this
secondary effect in the St = 0.2 case – resulting in a better match based on
just the SI growth rates alone – while the secondary effect is stronger than
the SI in the St = 0.04 case in which the actual observed growth timescales
are much shorter than those predicted from the turbulent SI prediction. We
think the VSSI is the leading candidate mechanism to explain this effect;
however, an investigation of this hypothesis is outside the scope of this
paper but should be considered in follow-up work.

Figure 26: Particle scale height versus time for several 3D simulations and a
semilog plot together with approximate baselines (dashed lines, for
$H_{p,\alpha}$ and approximate growth lines (dotted lines) described further
in text: (left panel) low to medium resolution runs for St = 0.2 , (middle
panel) low to medium resolution runs for St = 0.04. (Right panel) predicted
growth rates for the SI maximized over wavenumber (Umurhan et al., 2020) as a
function of St and $\alpha$, where the four simulations shown in the left two
panels are represented with symbols. Note all simulations appear on the
relatively active side of the SI.

##  6 Discussion

###  6.1 On turbulence in particle laden sheared midplane layers

The simulations conducted here appear to indicate that midplane settled
particle layers that are not strongly susceptible to the SI instead erupt via
other stratified axisymmetric fluid instabilities involving the vertical shear
of both the radial and azimuthal velocities. These instabilities appear to
lead to an early turbulent phase from which some secondary processes takes
root, possibly involving the SI, that leads to a sinusoidal pattern state of
varying degrees of coherence and characterized by some amount of radial drift.

The instability mechanisms leading to the early turbulent phase involve the
particle component as a collective agent, and mainly in a stabilizing role.
The flow analysis conducted indicates that the primary mechanism driving
instability depends on the St number of the particle component and although we
have examined only two values of St, we can confidently say that there are at
least two types of destabilizing processes: Based on the St = 0.04 suite of
simulations it appears that these settled layers primarily experience the
protoplanetary disk analog of the SymI, whose dynamic relies on the vertical
shear of the azimuthal velocity profile and requiring Ri $<1$. Based on the St
= 0.2 suite of numerical experiments the settled layers experience radial KH-
roll-up within layers located about 1-2 particle scale heights away from the
midplane. Further, the assessed Ri values within the settled particle layers
for St = 0.2 should support, in principle at least, azimuthal KH-roll-up as
well, but visualization of the flow fields do not present clear evidence for
it either.

During the review phase of this manuscript one of the reviewers commented on
the possibility that the non-turbulent SI could explain the early bounce phase
turbulent development. We have examined this possibility by doing a spot-check
on the predicted growth rates for the St = 0.04 case using the laminar SI
theory of YG2005. For the apparent lengthscales ($\lambda\sim 0.01H$) emerging
during the early bounce phase (e.g., see the $w_{g}$ field in Fig. 6 at
$t\Omega=6$), the predicted SI growth rate in the zero turbulent limit falls
in the range $0.05-0.1\Omega_{0}$, which is at least a factor of 5 weaker than
the corresponding growth rate of the SymI under those similar conditions. A
detailed examination of this matter, to better illustrate this case, deserves
to be done in the future.

In either case, however, the midplane layers exhibit unsteady motions –
apparently turbulent – that eventually leads into a secondary transition
resulting in a pattern forming state. We have examined whether the drifting
patterned state is an instance of the SI and the predictions also appear
mixed: for the St = 0.2 case the measured growth rates between the initial
putative turbulent phase and the final pattern state appear to be consistent
with SI under turbulent conditions if, however, systematically somewhat
faster-acting than predicted (e.g., following the theory of, Umurhan et al.,
2020; Chen & Lin, 2020). On the other hand, in the St = 0.04 case, the
corresponding measured transition phase growth rates are far more rapid than
those predicted by one to two order of magnitude (see Table 4). What is
responsible for this growth, and what is the cause for the dynamical
differences between these two St values, remains uncertain. Is it the SI
working in tandem with a secondary process, which senses the underlying
turbulent state in the mean? Might the secondary transition be an instance of
the proposed Vertical Shearing Streaming Instability (VSSI), a process
recently examined in a theoretical model by (Lin, 2021)? Furthermore, is it
possible that this emergent drifting pattern state is a numerical artifact
that goes away if the radial and/or azimuthal box scales are made larger?
These are all questions that need immediate resolution in future studies. In
Sec. 6.7 we further reflect on the possible action of KH-roll-up for Ri
$>1/4$.

Often rooted in simulations of relatively large St simulations (e.g., St
$>0.2$), it is a commonly held assumption that settled particle layers subject
to the SI are likely also susceptible to azimuthal KH-roll-up (e.g. Barranco,
2009; Lee et al., 2010a, b). Our findings suggest that there is, in fact, a
diversity of processes at play; and which actor dominates depends on the
layer’s St and local metallicity $Z$. Indeed, the St = 0.04, $Z=0.01$
simulations examined here appear to be driven entirely by the axisymmetric
SymI effect. Increasing St to 0.2 while keeping $Z$ fixed appears to lead to
both the SymI and radial KH-roll-up as the primary dynamical mechanisms.

###  6.2 On the relationship to the findings of Ishitsu et al. (2009)

In their unpublished study Ishitsu et al. (2009) examined the fate of
particle-gas setups like considered here, in which an imposed radial pressure
gradient induces a relative streaming between gas and particles. Their
simulations were 3D axisymmetric and considered two sets of St = 0.001 and
1.0. The early development of the layer in the St = 0.001 case is highly
reminiscent of the transition observed in both of our St suite of runs. In
particular, the upper right panel of Ishitsu et al. (2009)’s Figure 6 shows
off-midplane filament development where the particle layer’s vertical gradient
is greatest. This character is very similar to the dynamical structure
observed throughout the early development of both of our 3D axisymmetric and
full 3D simulations. For the St $=1.0$ simulation shown, their simulation
erupts into the kind of pattern characteristic of KH-roll-up as reported in
Barranco (2009, e.g., see Figure 7 of that work). No such KH-roll-up pattern
is discernible in the St = 0.001 run reported in Ishitsu et al. (2009).

In this sense we think these authors are justified in writing,

> “However, the instability of two-fluid shown in this work has the axis-
> symmetric unstable mode. As a result, the stabilization caused by the
> increase of the azimuthal wave number due to the radial shear is not
> effective. We expect that the instability occurs in the radial direction,
> and then the perturbation with small azimuthal wave number grows,” (Ishitsu
> et al., 2009, pg. 14).

Not only do our 3D axisymmetric simulations confirm this notion, but our full
3D investigations suggests that non-axisymmetric dynamics – e.g., like
azimuthal KH-roll-up – are not the primary instability mechanism driving
turbulence in low St midplane-settled particle layers, i.e., our findings
implicate axisymmetric modes as the primary driver of turbulence in full 3D
scenarios.

###  6.3 This study in relation to Lin (2021)

(Lin, 2021) undertake a two-fluid stability analysis of a model similar in
set-up to what has been considered here. In that study the author considers
the stability of an already settled particle layer whose particle equilibrium
is established within the framework of a local $\alpha$-disk shearing box
model often used to capture the essence of particle profiles within gas driven
turbulence (e.g., Dubrulle et al., 1995; Youdin & Lithwick, 2007; Laibe et
al., 2020) and as recently used in the analysis of the SI by Chen & Lin (2020)
and Umurhan et al. (2020). For input values of $Z$, St and $\alpha$ the
particle gas equilibrium in the radial-azimuthal mean admits solutions with
non-trivial vertical variations in all quantities on which a numerical
stability analysis is done. For conditions examined most similar to ours
(c.f., case C with the St = 0.04 runs treated here) Lin (2021) report that the
fastest growing mode has a horizontal wavenumber $kH_{g}\approx 1100$ with
growth rate $\sim 0.6\Omega$ and, most importantly, with a mode amplitude
becoming greatest at about 1-2 scale particle scale heights. Perhaps most
importantly, Lin (2021) show that the vertical shearing instability revealed
in the two-fluid setup is recovered within a single fluid framework as well.

While our findings agree in spirit with those reported in Lin (2021) there are
notable differences likely attributed to our differing steady states, since
ours are drawn directly from those exhibited by actively collapsing solutions
as opposed to those equilibria resulting from their adopted putative
turbulence model. Since our analysis includes neither particle diffusion nor a
model for turbulent viscosity, we predict that instability persists as
$m\rightarrow\infty$ with corresponding fastest growing radial wavenumber
similarly diverging, and this is similar in quality for the VSSI in the
inviscid limit as reported in Lin (2021). However, a notable difference in
predictions of the inviscid limit is that we predict growth rates that
asymptote to ${\cal O}\left({1}\right)$ values of $\Omega$ while the VSSI seem
to grow without bound as $kH\rightarrow\infty$ (e.g., see Figure 15 of Lin,
2021). Of course, viscosity due to any developing turbulence should temper
growth rates at short wavelengths, at least in reality.

Yet, the simple setup that went into motivating the SymI analysis developed in
section 5.2 captures the essence of what occurs during the settling and
transition phase of the numerical experiments conducted here, especially the
St = 0.04 simulations. Both studies demonstrate that an effectively single-
fluid process can lead to strong instability and quite likely explains the
dynamical source of turbulence in these midplane settled layers. Moreover, by
being an effectively single-fluid model, it shows that the process does not
need a relative stream between the two fluids in order to become active like
required for the SI (in this regard see recent conceptual advances of Squire &
Hopkins, 2018a, b). Despite this broad conceptual agreement, further work is
needed to reconcile these two theoretical approaches. However, it would seem
that the VSSI framework might be best applied in analyzing how an already
turbulent sublayer further develops as it already has built into it a model of
turbulence. Perhaps it might be used in explaining the secondary transition
into the drifting pattern state we have reported on here. This is grounds for
further investigation.

###  6.4 The Current work in the context of Garaud & Lin (2004)

Garaud & Lin (2004) conducted similar studies in order to investigate the
evolution of the settled dust layer in a two fluid approximation. Their work
used a mono-disperse population of small dust grains along with a gas disk
which is strictly laminar in the absence of solids. However, in their working
perturbation equations (their equations 15-19) the term containing the
coriolis force, and hence the effect of rotation in the problem is absent.
This is a major and the most important difference between their work and the
current one. The emergence of SymI in the work presented here is axisymmetric
in nature and emerges solely from the vertical variations of the azimuthal
flow velocity that has its root in the cross-velocity components in the
Coriolis term, whereas Garaud & Lin (2004) put their emphasis on the growth of
non-axisymmetric perturbations absent Corliolis influences. Hence the
evolution of the shear layer reported in their study belongs to a general
class of KH instability without any possibilities of SymI. Apart from this
point, Garaud & Lin (2004) put significant effort on the effect of cooling
time of the stability of the shear layer, whereas we adhere to a strictly
isothermal condition for all our analysis. It is noteworthy that they observe
the boundary of instability to extend beyond the classical limit of Ri =1/4
when cooling effects are included in simplified analytical modeling. Whether
these effects play a role in realistic disks remain to be explored.

###  6.5 A note on dust settling instability

Settling dust grains can also give rise to the dust settling instability
(DSI), a bi-product mechanism of the general class of resonant drag
instabilities (RDI, Squire & Hopkins, 2018a; Krapp et al., 2020). As noted in
(Squire & Hopkins, 2018a), the settling instability has two regimes: (i) for
modes in which $k_{x}w_{sx}\neq-k_{z}w_{sz}$, where ${\bf w}$ is the streaming
velocity, the growth rates are generally $\leq\Omega$ but, (ii) as
$k_{x}w_{sx}\rightarrow-k_{z}w_{sz}$ growth rates show divergence with
increasing absolute wavenumber $k$. In the latter case and so long as
$w_{sz}\neq 0$, Squire & Hopkins (2018a) predict that a rapid eruption of
activity will appear at the smallest lengthscales, with growth rates $\sim
k^{1/3}$ (for $\epsilon\ll 1)$. Krapp et al. (2020) observe the development of
fine scale structure in their highly resolved simulations centered on disk
sections centered at heights $z\sim H_{g}$ with box sizes $\sim 0.1H_{g}$.
Most of the activity of the DSI is expected to happen far from the midplane,
nominally at locations where the product combination of settling velocities
and local mean dust-densities are maximized. Moreover, the unrestricted growth
rates for increasing wavenumber suggests that resolving DSI may need more
resolution than what we have in our simulations. In our simulation results we
conjecture that the DSI may be operating, especially far from the midplane,
although, we have not notice any obvious signature in our simulation domain
spanning only $0.2~{}H_{g}$. However, in our analytical model, the primary
equations do not have any vertical velocity under equilibrium (see equations
C2 – C5) that can give rise to DSI. This lack of signs of the DSI in the
simulations reported here may be because the simulations take place close to
the midplane where the settling velocities are expected to be small since
$w_{sz}\sim\Omega H_{p}{\rm{St}}$, e.g., see Eq. (B6). However, with
sufficiently high numerical resolution the fastest growing mode ought to
appear as very short vertical wavenumber midplane parallel banded structures,
which given our current resolution capabilities, are not resolvable.
Additionally, there is as yet no Richardson Number analysis for the DSI based
on which, a more quantitative distinction between the two processes might be
made. Lastly, while the DSI can lead to disorder in the flow as the flow
settles, whether or not it can drive sustained midplane turbulence once the
main settling phase has receded is not expected.

###  6.6 On spectra and simulation convergence

The St = 0.2 kinetic energy power spectrum from Fig. 17 makes it evident that
for the 3D simulations, the power-law index of the inertial range changes when
we increase the resolution. The root cause of this effect, not being properly
pinpointed in this work, could be multifaceted. It is possible that the 3D
simulations are not yet resolved at $512^{3}$ resolution ($2560$ grids
$/H_{g}$) and that higher resolution is required. Another possibility could be
due to the number of particles used in the simulations with moderate and high
resolutions: In order to minimize the computation expense, the number of
particles used per grid point in the high resolution run (B3D-02H) is only
$0.125$ (i.e., 1 particle per 8 grids) compared to $1.0$ in the moderate
resolution run (See table 2). The gas and the particle fields in the system
communicate with each other through a drag term with a finite relaxation time
– Eqs. (2-4) – where the momentum exchange through the density fluctuations
are nonlinear, implying scale-to-scale energy transfer. So, the implementation
of a comparatively low number of super particles may influence the gas
dynamics. We disfavor this explanation because when settling is present,
particles cluster towards the disk midplane effectively reaching more than one
particle per grid. We also do not see any bottleneck effect in the gas kinetic
energy beyond the dissipation scale, so the possibility of back-scattering of
energy can be ruled out.

It is also interesting to note that for 3D-axisymmetric simulations, the
energy spectrum for both high and super-high resolution simulations are
converged with similar power-law index for the inertial range. For both St
$=0.2$ and $0.04$, the power-law index hovers between $2.1$ and $2.2$. The
small difference here for the two different St is consistent with the trend
reported by Pandey et al. (2019) where the authors found a dependence of the
energy spectrum on the overall mass loading effect, especially at the higher
wavenumbers. Note that the highest 3D axisymmetric simulations have $10240$
grids/$H_{g}$, so it is possible that for the full 3D simulations, going
beyond $512^{3}$ may show convergence in the energy spectrum. Given the
relatively minor change in power-law slope in the 3D-axisymmetric suite of
runs going from high ($512^{2}$) to super-high ($2048^{2}$) resolution runs,
we cautiously conjecture that the $512^{3}$ resolution full 3D simulation
might be close to convergence for the dynamically resolved scales.
Unfortunately, we were limited by available computational resources to carry
out simulations with resolution higher than $512^{3}$ and are unable to verify
this conjecture at this time nor characterize the turbulent behavior at high
resolution for St = 0.04.

It is also important to note that the character of the 3D axisymmetric
simulations in the inertial range are fundamentally different from the full 3D
ones, both in terms of slope of power spectrum and the location of the
integral ($k_{\rm integral}$) scale. Classic two-dimensional turbulence is
characterized by simultaneous downscale enstrophy cascade and upscale energy
cascade together with a $E_{k}\sim k^{-3}$ behavior in the inertial range
(Kraichnan & Montgomery, 1980), while in full 3D isotropic turbulence energy
cascades toward smaller scales with $E_{k}\sim k^{-5/3}$. 3D axisymmetric
scenarios in a disk might exhibit cascade properties similar to purely 2D
flows, but this is not quite certain as yet. The statistical behavior under
conditions where rotation and stratification are on equal footing is
complicated by wave dynamics, which present alternate pathways for scale-to-
scale energy exchange beyond just nonlinear velocity advection (i.e., inertial
effects). Little is known about how this unfolds under this dual influence in
both atmospheric/geophysical flows (e.g., see discussion of this for
atmospheric flows in Sec. 4.5 of Alexakis & Biferale, 2018) and much less in
dust-laden protoplanetary disk models. It is known that strongly stratified
flows, or dimensionally constrained settings can exhibit mixed/split and
anisotropic energy cascades while strong rotation can support some amount of
inverse cascade behavior under suitable forcing and/or vertical scales. In
these cases new exchange pathways emerge due to wave dynamics and produce
inertial range spectral slopes that deviate from that expected from pure 2D
turbulence ($k^{-3}$) or from full 3D Kolmogorov ($\sim k^{-5/3}$). For
example, under relatively strong stratification and moderate wavenumbers in
atmospheric models there exists the so-called Bolgiano-Obukhov scaling where
the kinetic energy behaves like $E_{k}\sim k^{-11/5}$ before it eventually
turns down to $E_{k}\sim k^{-5/3}$ at sufficiently high wavenumbers. In this
scenario energy exchange carries kinetic energy into potential energy through
wave transfer via gravity waves upon which nonlinear potential energy
advection forward spreads energy toward smaller scales. In our particle-laden
sheared setting there are nonlinearities in the gas-dust drag exchange terms.
Exactly how this dynamical pathway mitigates and/or directs the cascade of
energy – and how it ultimately shapes the resulting energy spectra – is yet to
be systematically examined in either 3D axisymmetric and/or full 3D settings.

In short: the energy content in these flows depends upon whether the equations
are being treated as full 3D or not, and that the spatial energy distribution
may not be Kolmogorov, especially at intermediate scales bridging energy
injection and the very short scales that should exhibit $k^{-5/3}$ Kolmogorov
character. Therefore, looking back upon our results: while large scale
structures look qualitatively similar in both axisymmetric and full 3D, a one-
to-one correspondence between the two may not be an accurate representation on
the small scales were particles are expected to cluster. As a result, we are
cautioned to infer very much about what happens on the small scales based on
the results of currently available high resolution axisymmetric simulations.
We expect better light on this will be shed upon the advent of high resolution
3D experiments or new simulation tools with wider dynamical scale resolution.

###  6.7 On turbulent activity where Ri $>1/4$

The persistence of activity in disks where the Richardson number exceeds the
classical critical value of $1/4$ remains enigmatic. Here, we briefly review
what is known about this feature of protoplanetary disk models and place our
findings with regards to the SymI into that context.

Gerbig et al. (2020) unequivocally demonstrated the emergence of sustained
(and probably turbulent) 3D activity in midplane settled dust layers
exhibiting Ri $>1/4$. The antecedents to this can be found in the three
studies of Gómez & Ostriker (2005), Johansen et al. (2006), and Barranco
(2009) where, in order to isolate and better understand operative physical
effects, the stability and nonlinear response of a midplane-settled particle
disk setting was examined within a restricted two-dimensional, azimuthal-
vertical slice. Initiated with vertical shears in the azimuthal velocity
field, over time it can be seen that the particle layer develops sinusoidal
undulations that eventually grow in amplitude and finally erupting to generate
strong vertical mixing (see especially Barranco, 2009). Johansen et al.
(2006), who examined this dynamic for the two-fluid model, calculate the
effective Richardson number as a function of disk height for simulations after
reaching their putative equilibrated statistical state, for which they find
that Ri $\approx 1$ within the confines of the disk containing most of the
dust. Barranco (2009) similarly report effective minimum Ri that significantly
exceed 1/4 for simulations in their well-developed stage.

Gómez & Ostriker (2005) perform a corresponding stability analysis in a
tractable one-fluid physical model and find that values of Ri as high as 5
could be linearly unstable. Similarly, Barranco (2009) report the possibility
of linear instability for values of Ri as high as 1.25. For the nonlinear
models presented in these three studies, the period of time from initiation of
a model run until the manifestation of full-fledged nonlinear layer
development ranges from 5 to 10 orbit times or, in our units, after
$t\Omega=30$-60. Similar roll-up times are reported for KH-roll-up in Johansen
et al. (2006).

Gómez & Ostriker (2005) and Barranco (2009) attribute to Coriolis effects the
tendency for such layers to be unstable even when ${\rm Ri}>1/4$. Indeed, the
classic KH-roll-up analysis involves the analysis of a two-dimenensional
vertically sheared stratified fluid. The Ri = $1/4$ criterion is technically
only appropriate to that simple setup sans rotation. The corresponding
incompressible linear stability problem, in the guise of the Taylor-Goldstein
equation (e.g., Garaud & Lin, 2004), is technically second order in time.
Introduction of Coriolis effects – i.e., in the way they appear in the disk
problem considered here and the aforementioned studies – raises the problem’s
temporal order by one, which certainly enriches the range of normal mode
behavior. As we have seen, the added complexity introduced by rotation is also
reflected in the dynamics responsible for the SymI. Indeed, the structure of
the perturbation equations for the SymI – i.e., Eqs. (C2-C5) – are similar to
the corresponding ones for KH-roll-up considered in the single-fluid setup of
Barranco (2009).

In light of these reflections, one may define a characteristic Rossby number
in terms of the amplitude and scale of the azimuthal velocity shear via

${\rm Ro}=\delta v_{00}\big{/}2\Omega H_{s},$ (62)

e.g., like those based on averaged fits to $\left<v_{g}\right>_{xy}$ or
$V_{{\rm cm}}$ discussed in Sec. 3.5 and Eq. (28) in particular. Indeed,
activity persists for Ri $>1/4$ for Ro $\gtrapprox 1$ – especially throughout
all phases of development in the St = 0.04 suite of simulations (see last
column of Table 3). This should be considered in comparison to the Rossby
number of Keplerian flow, ${\rm Ro}_{K}=3/4$. Thus it seems that for problems
of this sort the critical Richardson number should in general be a function of
the Rossby number, i.e., ${\rm Ri}_{{\rm c}}={\rm Ri}_{{\rm c}}({\rm Ro})$, in
which the classical stratified limit is recovered when rotation goes away

${\rm Ri}_{{\rm c}}\left({\rm Ro}\rightarrow\infty\right)=1/4.$ (63)

We think this is a worthwhile program – one that was instigated in Barranco
(2009) – for future clarification in the context of midplane-settled
protoplanetary disk modeling.

Nevertheless, we think the SymI dynamics unfolding in the St = 0.04 simulation
dominates the rotationally modified KH-roll-up likely simultaneously present.
While we have not done a detailed KH-roll-up analysis to complement that done
in Barranco (2009), we can see from that study that the predicted growth rates
for a configuration that most resembles our St = 0.04 simulations, with
$H_{p}\approx 0.01H_{g}$, predicts growth rates $\sim 0.1\Omega$ (see 3rd
column, 3rd row of Barranco, 2009) which is a factor of 5 times slower than
what is seen in our simulations.

We think the SymI is a part of the explanation for why such settled particle
layers go turbulent in the simulations conducted both here and in Gerbig et
al. (2020). All of our simulations manifest clear layer transition by 1
orbital period and fully developed nonlinear activity by 2-2.5 orbital periods
(i.e., $t\Omega=12$-15), which is far shorter than the time it takes non-
axisymmetric KH-roll-up to develop according to the above mentioned studies.
Therefore, we conjecture that 3D axisymmetric dynamics like the SymI and
radial KH-roll-up are the main drivers of activity, at least for the range of
St numbers considered here. A comprehensive study in this respect should start
with a theoretical single-fluid model that parses the relative importance of
KH-roll-up and the SymI in 3D-axisymmetry. This then should be followed with
understanding how these dynamics play out in the presence of two-stream
dynamics including the SI. In principle, the results of the single-fluid
theory ought to be contained in the generalized VSSI framework, particularly
within its single-fluid limit (Lin, 2021).

##  7 A Hard-boiled Summary with some final remarks

We provide summary bullet points of the major findings and conclusions of this
concerning small midplane sections of protoplanetary disks not subject to an
external source of turbulence:

  1. 1.

Midplane settled particle layers impart both upon itself (as a second fluid)
and upon the gas-fluid a complex Ekman pattern of radial and azimuthal jets in
the azimuthal-radial mean. The emergent epicyclically oscillating jet profiles
are a result of the momentum exchange between the gas and particles, and come
about while the particles are settling toward the midplane.

  2. 2.

The jet flows are subject both to classic KH roll-up as well as being unstable
to the so-called Symmetric Instability. The SymI relies on a mismatch in
density and azimuthal velocity isolines in a rotating atmosphere and is well
known to be relevant to mixing-layer dynamics in the ocean and frontogenesis
in the atmosphere, among other geophysical phenomena (e.g., Hoskins, 1974;
Bennetts & Hoskins, 1979; Thomas et al., 2013; Zeitlin, 2018; Stamper &
Taylor, 2017; Zhou et al., 2022).

  3. 3.

In the St = 0.2 simulations radial KH-roll-up appears to act in distinct
layers that are at least 2-3 $H_{p}$ away from the midplane and with
relatively  short growth rates $\sim\Omega$. This axisymmetric dynamic appears
to be present in the previously unpublished study of Ishitsu et al. (2009),
which has recently been examined in a new theoretical framework in Lin (2021).
These KH-roll-up dynamics play a primary role in driving instability of the
midplane region for St = 0.2, but appears to operate together with SymI.
Nonetheless, we note its significance in that this classical shear instability
acts on layers substantially removed from where the particles reside, contrary
to previously held assumptions or expectations.

  4. 4.

The axisymmetric SymI acts primarily within 1-2 $H_{p}$ of the midplane. Its
growth rate is relatively short at $\approx$1-2$\Omega$, it appears to be the
other main driver of midplane turbulence in the St = 0.2 simulations, and it
might be the primary –if not sole – driver of instability in the St = 0.04
runs as the Ri numbers are near 1 in this latter case. In the St $=0.04$ case,
we find that the SymI erupts into unsteady dynamics by $t\Omega=6$, which is
comparable to the KH-roll-up timescale in the St = 0.2 simulation, suggesting
that both processes are active and act on similar timescales. Thus, for the
numerical experiments conducted here we posit that while the shape of the mean
flows that develop in settling layers are St dependent, the following unstable
dynamics taking root in them are that of a single perfectly-coupled fluid
responding to a St-dependent flow profile.

  5. 5.

We note that the SymI is distinct from the SI as the latter emerges as a
result of momentum exchange between components drifting through one another,
while the former only requires strong vertical gradients in both the
perturbation azimuthal gas velocity ($v_{g}$) and mean density together with
isoline misalignment between particle density and total azimuthal gas velocity
($V_{K}+v_{g}$). In essence, the SymI is effectively a single-fluid baroclinic
dynamic.

  6. 6.

Furthermore, we have reason to suspect that the SymI we observe in these
layers is of the fundamentally same kind of mechanical process that drives the
VSI based on the connection of the latter with the “sloping convection” effect
well-known in geophysical fluid flows Yellin-Bergovoy et al. (2021).

  7. 7.

The analytical theory we developed indicates that a layer will experience
strong widespread instability in the particle layer due to the SymI once ${\rm
Ri}<1$. The fastest growing mode’s growth rate is well approximated by Eq.
(57), $\sqrt{2/{\rm Ri}-2}\Omega$, where Ri is based on the azimuthal
velocity’s vertical shear. The fastest growing wavelength depends on integer
vertical mode number $m>0$ and is given in Eq. (58).

  8. 8.

We believe this finding goes toward explaining the results reported in Gerbig
et al. (2020) where they find that settled turbulent particle layers have
effective values of ${\rm Ri}$ to be well above $1/4$ but less than $1$,
especially for simulations where $Z>0.01$ (e.g., see Figure 6 of Gerbig et
al., 2020). We therefore suspect that the SymI plays a role –and a prominent
one depending on St – in driving midplane turbulence in all numerical
computational studies of the SI in otherwise laminar disk models.

  9. 9.

Given the above findings, we therefore conclude that these particle settled
midplane layers are subject to three simultaneously acting instabilities: the
above mentioned two, and the SI. Based on the simulations with the specific
input parameters we have conducted here, we conjecture that the SI is the
weakest of the three with the slowest growth rates (Chen & Lin, 2020; Umurhan
et al., 2020). We conjecture that whenever the parameters for ${\rm{St}}$ and
$Z$ (or $Z/\Pi$ as proposed in Sekiya & Onishi, 2018) permit the SI with
relatively fast growth rates, it does so out of a turbulent state driven by
the other two shear instabilities.

  10. 10.

We have produced space-time diagrams of the azimuthally averaged particle
surface density field and have examined the azimuthally averaged particle
surface densities at late times. As has been done before, in those simulations
that admit the SI and allow for its nonlinear development and saturation
(i.e., the $Z=0.01,{\rm{St}}=0.2$ simulation), its signature is clearly
visible by the emergence and maintenence of radially drifting coherent
particle overdensities. In the high resolution simulation of
$Z=0.01,{\rm{St}}=0.04$, the spacetime diagram indicates intermittent
manifestation of overdensities – that appear to form, drift and later
dissipating only to reform and repeat this dynamic once again. At this stage
it is unclear whether this observed intermittency is a result of the SI
struggling to emerge, or if it is some other collective turbulent effect that
momentarily drives particle overdensities only to be later destroyed (e.g.,
Yang et al., 2018). This requires further analysis.

  11. 11.

We have calibrated the PENCIL code used here by simulating 3D Kolmogorov
turbulence in a non-rotating frame free of particles by forcing the simulation
at approximately 1/3 the box size. We recover the expected $n_{g}=5/3$
inertial range power law behavior (i.e., $\varepsilon_{k,g}\sim k^{-n_{g}}$)
down to a wavenumber $k_{N,D}=k_{N}/4$, where $k_{N}$ is the simulation’s
Nyquist wavenumber. At lengthscales smaller than $2\pi/k_{N,D}$ the gas
kinetic energy plummets super exponentially. It therefore means that down-
scale propagating turbulent kinetic energy does not appreciably reach scales
less than 6-8 grid points in these simulations. All subsequent analyses,
including statistical measures of underlying turbulence, are therefore here
restricted to scales larger than $2\pi/k_{N,D}$.

  12. 12.

For the particle-gas simulations done here with box sizes $=0.2H$ we
conjecture that all medium sized $256^{3}$ simulations are not statistically
converged. We have confirmed this to be the case in for the St = 0.2
simulation as the nominal inertial range power-law slopes for both the
particle and gas kinetic energies ($\varepsilon_{k,p}\sim k^{-n_{p}}$ and
$\varepsilon_{k,g}\sim k^{-n_{g}}$, respectively) steepen for the
corresponding high resolution $512^{3}$ run. Unfortunately due to limitations
of available resources we were not able to confirm this trend for a
concomitant high resolution run in the St = 0.04 case.

  13. 13.

We have run a corresponding series of high ($512^{2}$ element) and super-high
($2048^{2}$ element) resolution 3D-axisymmetric simulations for which we find
that the power-law slope of the gas kinetic energy appears to be converged at
$512^{2}$: For St = 0.2 we find $n_{g}\approx 2.1\pm 0.05$; while for St =
0.04, $n_{g}\approx 2.15\pm 0.05$. The inertial range particle kinetic
energies appears converged for the St = 0.04 case with $n_{p}\approx 0.42\pm
0.05$. However, the corresponding inertial range power-law slope shows
steepening for St = 0.2 with $n_{p}=0.92\pm 0.1$ for high resolution and
$n_{p}=1.1\pm 0.05$ at super-high resolution.

  14. 14.

Based on this 3D-axisymmetric finding we cautiously conjecture that the
inertial range power slope reported for the St = 0.2 high resolution full-3D
simulation might be converged or close to convergence with $n_{g}=1.37\pm
0.03$ and $n_{p}=0.64\pm 0.04$. In that case the inertial range appears to
emerge at $k\approx 200H_{g}^{-1}$ corresponding to about 1/6 the box size,
which is roughly 3 particle scale heights viewed symmetrically from the
midplane. What these inertial range behaviors say about the nature of
turbulent kinetic energy cascade remains to be understood.

Much of the findings reported in this study owe their illumination to the
earlier axisymmetric study of Ishitsu et al. (2009). It therefore strongly
suggests that axisymmetric instabilities primarily operate in midplane settled
particle layers in gaseous protoplanetary disks not subject to external
sources of turbulence and where the SI is at best only weakly operative. The
recent results of Gerbig et al. (2020) appears to confirm that such layers are
sufficiently active even though the SI is weakly operating. If such weakly
turbulent disk conditions are realizable in realistic protplanetary disk
models, then it remains a challenge to understand how low St number conditions
can lead to planetesimal formation. Indeed, global evolution modeling of
particle growth in turbulent disks show that disks maintain low St particles
for weak-to-moderate levels of external turbulence during the first million
years after disk formation (e.g., Estrada et al., 2016; Sengupta et al.,
2019). If such externally driven turbulence scenarios shut-off for some
reason, then such low St number particles would settle to the midplane, then
to be subjected to the particle-shear driven turbulent state (discussed thus
far) that could conceivably further act to thwart planetesimal accumulation.
We are therefore swayed by the concerns stated in the conclusions of Ishitsu
et al. (2009), where they write “Thus, even though the global turbulence is
weak in the dead zone, the turbulence due to the instability described in this
paper may play the role of avoiding planetesimal formation and floating dust
in the disk.”

##  8 ACKNOWLEDGEMENT

We are grateful to P.E. Estrada, J.N. Cuzzi and K. Shariff of NASA-ARC for
many illuminating discussions and institutional support. D.S. was supported by
the NASA Postdoctoral Program (NPP) fellowship, and partly by the NASA
Astrobiology Institute. O.M.U. acknowledges the NASA Planetary Science
Division ISFM for Planet Formation and Exoplanets Theory at NASA/Ames for
financial and computational support, and also, the NASA ROSES TCAN grant
“Dynamical instabilities in the aid of planet formation in circumstellar
disks” (20-TCAN20-0011). We thank the anonymous referee for many insightful
comments and suggestions which significantly improved the quality of this
manuscript. We are also indebted to insights shared with us by our TCAN
collaboration colleagues W. Lyra, J. Simon, C.C. Yang, and A. Youdin. All the
simulations presented in this paper are performed on the NASA Advanced
Supercomputing (NAS) facility with generous computational resources provided
through NPP and ISFM allocations.

##  Appendix A Reformulating particle fluid evolution equations in terms of
$\mu$

Our aim is to rewrite the particle-fluid evolution equations Eq.(3-4) in terms
of the variable $\hat{\mu}_{{}_{i}}\equiv\sqrt{\rho_{p}}U_{pi}$, i.e., in
terms of the full velocity field including the background Keplerian flow. We
note that within this section we use the usual Einstein index convention and,
as such, the dummy indices $i$ and $j$ used here are not to be confused with
the particle labeling indices $i$ or with the grid labeling indices $j$ used
extensively in Section. 2.1. Leaving the RHS of Eq.(4) in terms of a forcing
function $F_{{}_{i}}$ we input this definition into the total momentum
conservation and manipulate the result accordingly. The purpose of this
procedure is to have the equations appear in terms of conjugate symmetric
variables in Fourier space that make assessing the energy contained in the
particle component straightforward. We will take the result we develop in
terms of $\sqrt{\rho_{p}}U_{pi}$ and then restore the definition
$U_{pi}=V_{{}_{{\rm K}}}\delta_{{}_{i2}}+u_{pi}$ to recover the form we seek
in terms of the perturbation velocities and
$\mu_{{}_{i}}\equiv\sqrt{\rho_{p}}u_{pi}$. Starting with

$\displaystyle\partial_{{}_{t}}\rho_{p}U_{pi}+\partial_{j}\rho_{p}U_{pj}U_{pi}+\epsilon_{ijm}\rho_{p}2\Omega_{j}U_{pm}=F_{pi}$
(A1)

in which $\Omega_{i}=\delta_{{}_{i3}}\Omega_{0}$ and where $\epsilon_{ijk}$ is
the Levi-Civita symbol. We have after rewriting the above in terms of
$\hat{\mu}_{{}_{i}}$ that

$\displaystyle F_{pi}$ $\displaystyle=$
$\displaystyle\rho_{p}^{1/2}\partial_{{}_{t}}\hat{\mu}_{{}_{i}}+\small\mbox{$\frac{1}{2}$}\rho_{p}^{-1/2}\hat{\mu}_{{}_{i}}\partial_{t}\rho_{p}+\partial_{j}\hat{\mu}_{{}_{j}}\hat{\mu}_{{}_{i}}+\epsilon_{ikm}\rho_{p}^{1/2}2\Omega_{j}\hat{\mu}_{{}_{m}}.$
(A2)

We similarly re-express the dust continuity equation as

$\partial_{t}\rho_{p}+\hat{\mu}_{{}_{j}}\partial_{j}\rho_{p}^{1/2}+\rho_{p}^{1/2}\partial_{j}\hat{\mu}_{{}_{j}}=0.$
(A3)

Replacing the time derivative of $\rho_{p}$ in Eq. (A2) by the above
expression, followed by factoring out $\rho_{p}^{1/2}$ reveals

$F_{pi}=\rho_{p}^{1/2}\left[\partial_{t}\hat{\mu}_{{}_{i}}+\rho_{p}^{-1/2}\partial_{j}\hat{\mu}_{j}\hat{\mu}_{i}-\small\mbox{$\frac{1}{2}$}U_{pi}U_{pj}\partial_{j}\rho^{1/2}-\small\mbox{$\frac{1}{2}$}U_{pi}\partial_{j}\hat{\mu}_{{}_{j}}+\epsilon_{ijm}2\Omega_{j}\hat{\mu}_{{}_{m}}\right].$
(A4)

We may rewrite the third term on the RHS of the above equation through a
number of derivative by parts maneuvers to find

$\displaystyle-\small\mbox{$\frac{1}{2}$}U_{pi}U_{pj}\partial_{j}\rho_{p}^{1/2}$
$\displaystyle=$
$\displaystyle-\small\mbox{$\frac{1}{2}$}\partial_{{}_{j}}\Big{(}U_{pi}\hat{\mu}_{{}_{j}}\Big{)}+\small\mbox{$\frac{1}{2}$}\rho_{p}^{1/2}\left(U_{pi}\partial_{j}U_{pj}+U_{pj}\partial_{j}U_{pi}\right)$
(A5) $\displaystyle=$
$\displaystyle-\small\mbox{$\frac{1}{2}$}\Big{(}\hat{\mu}_{{}_{j}}\partial_{j}U_{pi}+U_{pi}\partial_{j}\hat{\mu}_{{}_{j}}\Big{)}+\small\mbox{$\frac{1}{2}$}\Big{(}\hat{\mu}_{{}_{j}}\partial_{j}U_{pi}+\hat{\mu}_{{}_{j}}\partial_{j}U_{pi}\Big{)}$
$\displaystyle=$
$\displaystyle-\small\mbox{$\frac{1}{2}$}U_{pi}\partial_{j}\hat{\mu}_{{}_{j}}+\small\mbox{$\frac{1}{2}$}\hat{\mu}_{{}_{i}}\partial_{j}U_{pj},$

while

$\rho_{p}^{-1/2}\partial_{j}\hat{\mu}_{j}\mu_{i}=U_{pj}\partial_{j}\hat{\mu}_{{}_{i}}+U_{pi}\partial_{j}\hat{\mu}_{{}_{j}}.$
(A6)

Putting these all together yields a reworked momentum evolution equation

$\partial_{t}\hat{\mu}_{{}_{i}}+U_{pj}\partial_{j}\hat{\mu}_{{}_{i}}+\small\mbox{$\frac{1}{2}$}\hat{\mu}_{{}_{i}}\partial_{j}U_{pj}+\epsilon_{ijm}2\Omega_{j}\hat{\mu}_{{}_{m}}=F_{{}_{i}}\Big{/}\sqrt{\rho_{p}}.$
(A7)

Finally, making the replacements
$\hat{\mu}_{{}_{i}}\rightarrow\mu_{{}_{i}}+\delta_{{}_{i2}}\sqrt{\rho_{p}}V_{{}_{{\rm
K}}}$ and $U_{pi}\rightarrow u_{pj}+\delta_{{}_{i2}}V_{{}_{{\rm K}}}$, we have

$\partial_{t}\mu_{{}_{i}}+u_{pj}\partial_{j}\mu_{{}_{i}}+V_{{}_{{\rm
K}}}\partial_{y}\mu_{{}_{i}}+\small\mbox{$\frac{1}{2}$}\mu_{{}_{i}}\partial_{j}v_{{}_{j}}+\biggl{(}V_{{}_{{\rm
K}}}\partial_{{}_{t}}\sqrt{\rho_{p}}+\small\mbox{$\frac{1}{2}$}\sqrt{\rho_{p}}V_{{}_{{\rm
K}}}\partial_{j}u_{pj}+u_{px}\partial_{x}\sqrt{\rho_{p}}V_{{}_{{\rm
K}}}\biggr{)}\delta_{{}_{i2}}+\epsilon_{{}_{ijm}}2\Omega_{j}\mu_{{}_{m}}-2\Omega_{0}\sqrt{\rho_{p}}V_{{}_{{\rm
K}}}\delta_{{}_{i1}}=F_{{}_{i}}\Big{/}\sqrt{\rho_{p}},$

which after making use of Eq. (A3) becomes finally

$\displaystyle\partial_{t}\mu_{{}_{i}}+u_{pj}\partial_{j}\mu_{{}_{i}}+V_{{}_{{\rm
K}}}\partial_{y}\mu_{{}_{i}}+\small\mbox{$\frac{1}{2}$}\mu_{{}_{i}}\partial_{j}u_{pj}+\epsilon_{{}_{ijm}}2\Omega_{j}\mu_{{}_{m}}$
$\displaystyle\hskip
56.9055pt+\biggl{(}u_{px}\partial_{x}\sqrt{\rho_{p}}V_{{}_{{\rm
K}}}-V_{{}_{{\rm K}}}u_{pj}\partial_{j}\sqrt{\rho_{p}}-V_{{}_{{\rm
K}}}^{2}\partial_{y}\sqrt{\rho_{p}}\biggr{)}\delta_{{}_{i2}}-2\Omega_{0}\sqrt{\rho_{p}}V_{{}_{{\rm
K}}}\delta_{{}_{i1}}=F_{{}_{i}}\Big{/}\sqrt{\rho_{p}}.$ (A8)

Energy statistics on the particle fluid may be constructed based on either of
the two formulations found in Eq. (A7) or Eq. (A8) depending upon the context
of interest; whether it be the total energies (the former) or the perturbation
energies (the latter). In particular, these expressions will be Fourier
transformed (e.g., the transform of $\mu_{i}({\bf x})$ =
$\tilde{\mu}_{{}_{k,i}}({\bf k})$) and partial energies and their scale-to-
scale transfer can be assessed within the usual prescribed narrow wavenumber
bins.

##  Appendix B Time dependent settling solutions

We consider solutions to the horizontally uniform steady settling solutions of
the system. We assume that the gas density is constant and write it here as
$\overline{\rho}_{g}$, and that the gas field has no vertical velocity
component. We use capital letters to designate these solutions, e.g., for the
horizontal gas velocity we have $u_{g}=U_{g}(z,t)$, and so forth for the other
variables. Since by assumption $W_{g}=0$, the evolution equations for the gas
are

$\displaystyle\partial_{t}U_{g}+2\Omega_{0}V_{g}$ $\displaystyle=$
$\displaystyle-\Pi_{r}-\frac{\Omega_{0}\varepsilon}{{\rm{St}}}\Big{(}U_{g}-U_{p}\Big{)},$
(B1)
$\displaystyle\partial_{t}V_{g}-\small\mbox{$\frac{1}{2}$}\Omega_{0}U_{g}$
$\displaystyle=$
$\displaystyle-\frac{\Omega_{0}\varepsilon}{{\rm{St}}}\Big{(}V_{g}-V_{p}\Big{)},$
(B2)

where the particle to gas ratio, written here as
$\varepsilon=\rho_{p}(z,t)/\overline{\rho_{g}}$, differs from the its general
definition in the text (i.e., $\epsilon$), since $\rho_{g}$ is also a
variable. In this way $\varepsilon$ is a proxy for the particle density
$\rho_{p}$. For the particle component, the particle density continuity
equation for the particle density is

$\partial_{t}{\varepsilon}+\partial_{z}\left(W_{p}\varepsilon\right)=0,$ (B3)

while for the particle momentum equations we have

$\displaystyle\partial_{t}U_{p}+W_{p}\partial_{z}U_{p}+2\Omega_{0}V_{p}$
$\displaystyle=$
$\displaystyle-\Pi_{r}-\frac{\Omega_{0}}{{\rm{St}}}\Big{(}U_{p}-U_{g}\Big{)},$
(B4)
$\displaystyle\partial_{t}V_{p}+W_{p}\partial_{z}V_{p}-\small\mbox{$\frac{1}{2}$}\Omega_{0}U_{p}$
$\displaystyle=$
$\displaystyle-\frac{\Omega_{0}}{{\rm{St}}}\Big{(}V_{p}-V_{g}\Big{)},$ (B5)
$\displaystyle\partial_{t}W_{p}+W_{p}\partial_{z}W_{p}$ $\displaystyle=$
$\displaystyle-\frac{\Omega_{0}}{{\rm{St}}}W_{p}-\Omega_{0}^{2}z.$ (B6)

We begin by analyzing Eq. (B6) by noting that $W_{p}=-\beta\Omega_{0}z$ is a
solution in which $\beta$ is the solution of the time evolution equation.

$\partial_{t}\beta=\Omega_{0}\left(\beta^{2}-\frac{\beta}{{\rm{St}}}+1\right).$
(B7)

Solutions to this equation depend upon whether or not ${\rm{St}}\leq 1/2$ Lin
(2021). We define

$\displaystyle\beta_{\pm}=\frac{1\pm\delta}{2{\rm{St}}};\quad\delta\equiv\sqrt{1-4{\rm{St}}^{2}}.$
(B8)

When ${\rm{St}}\leq 1/2$ the time-asymptotic stable settling solution is given
by $\beta=\beta_{-}$ Lin (2021) 555For ${\rm{St}}\leq 1/2$ the solution
$\beta=\beta_{+}$ corresponds to exponentially growing perturbations.. When
${\rm{St}}>1/2$ the solution is given by

$\beta=\frac{1}{2{\rm{St}}}\left[1+|\delta|\tan\left(\frac{|\delta|\Omega_{0}t}{2{\rm{St}}}\right)\right],\qquad\left|\frac{|\delta|\Omega_{0}t}{2{\rm{St}}}\right|\leq\frac{\pi}{2},$
(B9)

while noted here, we will not consider these solutions any further in this
study. With the purely linear dependence on $z$ for $W_{p}$, the solution to
Eq. (B3) has a Gaussian form

$\varepsilon=\varepsilon_{0}(t)\exp\left(-\frac{z^{2}}{2H_{p}^{2}}\right),$
(B10)

in which

$\frac{d\ln\varepsilon_{0}}{dt}=-\frac{d\ln H_{p}}{dt}=\beta.$ (B11)

The proof of this solution follows by inserting the solution form Eq. (B10)
together with $W_{p}=-\beta\Omega_{0}z$ into Eq. (B3), collecting the
resulting expression into like powers of $z$ followed by setting their
coefficients to zero, which produces the time dependencies on
$\varepsilon_{0}$ and $H_{p}$ found in Eq. (B11). We note that the total
vertically integrated mass is conserved as the product
$\displaystyle\int_{-\infty}^{\infty}\rho_{d}(z,t)dz=\sqrt{2\pi}\rho_{g}\varepsilon_{0}\cdot
H_{p}$ is always time-independent constant for $\rho_{g}$ constant. Because we
will restrict our attention to problems with ${\rm{St}}<1/2$, in all of our
subsequent analyses we will assume that $\beta$ achieves its time-asymptotic
value $\beta_{-}$, in which case

$\varepsilon_{0}=\varepsilon_{{}_{00}}e^{\displaystyle\beta\Omega_{0}t},\qquad
H_{p}=H_{p0}e^{\displaystyle-\beta\Omega_{0}t}.$ (B12)

The constant coefficients $\varepsilon_{{}_{00}}$ and $H_{p0}$ relate to one
another on the assumption that the vertically integrated particle density
$\Sigma_{p}=\sqrt{2\pi}\varepsilon_{{}_{00}}H_{p0}\rho_{g}(0)$ is constant, in
which case

$\varepsilon_{{}_{00}}=Z\left(\frac{H_{g}}{H_{p0}}\right),\qquad
Z\equiv\frac{\Sigma_{p}}{\Sigma_{g}},$ (B13)

where $Z$ is the local model disk metallicity (see text). Thus, solutions to
the remaining equation set are defined by the parameters $Z,H_{p0}$ and
${\rm{St}}$.

The solutions for the remaining flow variables, $U_{g},U_{p},V_{g},V_{p}$, may
be determined from solving equations (B1-B2) and (B4-B5) after inserting for
$\epsilon$ the solution for $\varepsilon$ found in Eq. (B10), as well as
replacing $W_{p}$ with $-\beta\Omega z$ where $\beta=\beta_{-}$, (Eq. B8).
Since the problem involves the solution of a settling layer, it proves
beneficial to go into a spatial coordinate frame that follows this descent as
the particle scale height $H_{p}(t)$ shrinks over time. For example

$U_{g}(z,t)\rightarrow
U_{g}(\varphi,t);\qquad\displaystyle\varphi=\exp\left(-\frac{z^{2}}{2H_{p}^{2}}\right),$
(B14)

and similarly for $U_{p},V_{g},V_{p}$. $\varphi$ is now an independent
variable that varies between $0$ and $1$. In these new coordinates, partial
derivatives are replaced according to their coordinate transformed forms

$\displaystyle z\partial_{z}\longrightarrow
z\left(\frac{\partial\varphi}{\partial
z}\right)\partial_{\varphi}=-2\varphi\ln{\varphi}\partial_{\varphi},$
$\displaystyle\partial_{t}\longrightarrow\partial_{t}+\left(\frac{\partial\varphi}{\partial
t}\right)\partial_{\varphi}=\partial_{t}-2(\beta\varphi\ln{\varphi})\partial_{\varphi}.$
(B15)

In these new coordinates Eqs. (B1-B2, B4-B5) are now reexpressed as

$\displaystyle\partial_{t}U_{g}-W_{\varphi}\partial_{\varphi}U_{g}-2\Omega_{0}V_{g}$
$\displaystyle=$
$\displaystyle-\Pi_{r}-\frac{\Omega_{0}}{{\rm{St}}}\varepsilon_{0}\varphi\Big{(}U_{g}-U_{p}\Big{)},$
(B16)
$\displaystyle\partial_{t}V_{g}-W_{\varphi}\partial_{\varphi}V_{g}+(1/2)\Omega_{0}U_{g}$
$\displaystyle=$
$\displaystyle\frac{\Omega_{0}\varepsilon_{0}}{{\rm{St}}}\varphi\Big{(}V_{g}-V_{p}\Big{)},$
(B17) $\displaystyle\partial_{t}U_{p}-2\Omega_{0}V_{p}$ $\displaystyle=$
$\displaystyle\frac{\Omega_{0}}{{\rm{St}}}\Big{(}U_{p}-U_{g}\Big{)},$ (B18)
$\displaystyle\partial_{t}V_{p}+(1/2)\Omega_{0}U_{p}$ $\displaystyle=$
$\displaystyle\frac{\Omega_{0}}{{\rm{St}}}\Big{(}V_{p}-V_{g}\Big{)},$ (B19)

where the speed $W_{\varphi}\equiv-2\beta\varphi\ln\varphi$ is always greater
than zero. Note how in this transformed coordinate system the vertical
advection by $W_{p}$ of the horizontal particle velocities $U_{p},V_{p}$ is
now transformed into an effective upward advection of the horizontal gas
velocities, with no more vertical advection in the particle component. This
makes sense because we have moved into a reference frame that follows the
evolving particle scale height. The solution to the above set of equations are
sought subject to the condition that
$\partial_{z}U_{g}\big{|}_{z=0}=\partial_{z}V_{g}\big{|}_{z=0}=0$. In all
cases we have determined, these gas velocity conditions automatically impose
the same conditions at the midplane for $U_{p},V_{p}$ (i.e., at $\varphi=1$).
We use the convention described in Umurhan et al. (2020) and write
$\Pi_{r}=-2\delta c_{s}\Omega_{0}$, where $\delta$ is the local disk opening
angle. In the limit as $z\rightarrow\infty$ we make sure the solution behave
according to the Nakagawa solutions in the limit $\epsilon\rightarrow 0$
(Nakagawa et al., 1986; Youdin & Goodman, 2005; Umurhan et al., 2020; Lin,
2021), i.e.,

$U_{g}(\varphi\rightarrow 0)=0,\quad V_{g}(\varphi\rightarrow 0)=-\delta
c_{s},\quad V_{p}(\varphi\rightarrow 0)=-\frac{\delta
c_{s}}{1+{\rm{St}}^{2}},\quad U_{p}(\varphi\rightarrow
0)=-\frac{2{\rm{St}}\delta c_{s}}{1+{\rm{St}}^{2}},$ (B20)

Evidently, all velocity variables are scaled by $\delta c_{s}$. Henceforth we
assume $\delta=0.05$, noting that this choice has little influence on the
qualitative solutions we demonstrate forthwith.

We note that a quick eigenvalue analysis of Eq. (B18-B19) – setting the
forcing terms $U_{g},V_{p}\rightarrow 0$ shows that the eigenvalues, $\sigma$,
are decaying epicyclic oscillations: $\sigma=\pm
i\Omega_{0}-\Omega/{\rm{St}}$. Thus, drag forcing by the gas induces
temporally decaying epicyclic motions in the particle fluid. The memory of the
forcing decays away on a time scale ${\rm{St}}/\Omega_{0}$. A cursory
inspection of Eq. (B16-B17) indicates that qualitatively similar response
occurs in the gas fluid due to the drag forcing by the particles, except that
the decay timescales depend upon the local value of $\varepsilon$, which is a
function of height.

_Numerical Method._ We use a third order correct upwind differencing scheme to
stably calculate the advection terms in (B16-B17). Given the sign of the
derivative terms, this amounts to a forward differencing scheme on a grid
$0\leq\varphi_{i}\leq 1$. To better capture the decaying behavior of the
solutions as $\varphi\rightarrow 0$, we further move into a stretched
coordinate system in which

$\varphi_{i}=\exp\left(1-\frac{1}{\zeta_{i}}\right),$ (B21)

where $\zeta_{i}$ are points on a uniform grid between 0 and 1, endpoints
included. We also note that

$\frac{\partial\zeta}{\partial\varphi}=\frac{1}{\varphi}\left(\frac{1}{1-\log\varphi}\right)^{2},$
(B22)

which is used in the actual differencing scheme described further below.

We discretize in time with time step $\Delta t$, where the time at the $n^{\rm
th}$ time step is given by $t_{n}=n\Delta t$. All dependent variables at time
step $n$ are denoted with superscripts “$n$”. Spatial values of the dependent
variables are designated by the subscript “$i$”, where $i$ runs from $1$ to
$N$, the total number of grid points in the domain. We implement an
exponential time integrator common in varied applied mathematical studies (Cox
& Matthews, 2002) including some astrophysical applications (e.g., Umurhan &
Regev, 2004; Umurhan et al., 2007). Thus for the particle component the
evolution follows

$\displaystyle U_{p,i}^{n+1}$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{\Delta
t\Omega_{0}}{{\rm{St}}}\right)\cdot\Big{[}U_{p,i}^{n}\cdot\cos\left(\Delta
t\Omega_{0}\right)+2V_{p,i}^{n}\sin\left(\Delta t\Omega_{0}\right)\Big{]}+$
(B23) $\displaystyle\ \ \ \ \ \ \ \ \frac{\Delta
t\Omega_{0}}{{\rm{St}}}\cdot\exp\left(-\frac{\Delta
t\Omega_{0}}{2{\rm{St}}}\right)\cdot\left[U_{g,i}^{n}\cdot\cos\left(\frac{\Delta
t\Omega_{0}}{2}\right)+2V_{g,i}^{n}\sin\left(\frac{\Delta
t\Omega_{0}}{2}\right)\right],$ $\displaystyle V_{p,i}^{n+1}$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{\Delta
t\Omega_{0}}{{\rm{St}}}\right)\cdot\Big{[}-\frac{1}{2}U_{p,i}^{n}\cdot\sin\left(\Delta
t\Omega_{0}\right)+V_{p,i}^{n}\cos\left(\Delta t\Omega_{0}\right)\Big{]}+$
(B24) $\displaystyle\ \ \ \ \ \ \frac{\Delta
t\Omega_{0}}{{\rm{St}}}\cdot\exp\left(-\frac{\Delta
t\Omega_{0}}{2{\rm{St}}}\right)\cdot\left[-\frac{1}{2}U_{g,i}^{n}\cdot\sin\left(\frac{\Delta
t\Omega_{0}}{2}\right)+V_{g,i}^{n}\cos\left(\frac{\Delta
t\Omega_{0}}{2}\right)\right].$

The evolution of gas velocity quantities follows

$\displaystyle U_{g,i}^{n+1}$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{\Delta
t\Omega_{0}}{{\rm{St}}}\varepsilon(t_{n})\varphi_{i}\right)\cdot
U^{n}_{g,i}+\Delta t\exp\left(-\frac{\Delta
t\Omega_{0}}{2{\rm{St}}}\varepsilon(t_{n})\varphi_{i}\right)\cdot
F^{n}_{g,i},$ (B25) $\displaystyle V_{g,i}^{n+1}$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{\Delta
t\Omega_{0}}{{\rm{St}}}\varepsilon(t_{n})\varphi_{i}\right)\cdot
V^{n}_{g,i}+\Delta t\exp\left(-\frac{\Delta
t\Omega_{0}}{2{\rm{St}}}\varepsilon(t_{n})\varphi_{i}\right)\cdot
G^{n}_{g,i},$ (B26)

where

$\displaystyle F^{n}_{g,i}$ $\displaystyle=$ $\displaystyle
W_{\varphi,i}\Big{(}D_{\varphi}U^{n}_{g}\Big{)}_{i}+2\Omega_{0}V^{n}_{g,i}-\Pi_{r}+\frac{\Omega_{0}}{{\rm{St}}}\varepsilon_{{}_{0}}(t_{n})\varphi_{i}U^{n}_{p,i}$
(B27) $\displaystyle G^{n}_{g,i}$ $\displaystyle=$ $\displaystyle
W_{\varphi,i}\Big{(}D_{\varphi}V^{n}_{g}\Big{)}_{i}-(1/2)\Omega_{0}V^{n}_{g,i}+\frac{\Omega_{0}}{{\rm{St}}}\varepsilon_{{}_{0}}(t_{n})\varphi_{i}V^{n}_{p,i}.$
(B28)

The derivative stencil is given by a third order forward differencing scheme
in $\varphi$. Because we go into the stretched coordinates defined in
equations (B21-B22), this derivative operation is written out as

$\Big{(}D_{\varphi}U^{n}_{g}\Big{)}_{i}\longrightarrow\left(\frac{\partial\zeta}{\partial\varphi}\right)_{i}\left(-\frac{11}{6}U^{n}_{g,i}+3U^{n}_{g,i+1}-\frac{3}{2}U^{n}_{g,i+2}+\frac{1}{3}U^{n}_{g,i+3}\right),$
(B29)

and similarly for $\displaystyle\Big{(}D_{\varphi}V^{n}_{g}\Big{)}_{i}$. All
numerical solutions for $U_{g},V_{g},U_{p},V_{p}$ are initiated with the
spatially uniform solutions found in Eq. (B20). Also, all solutions displayed
in the text were derived on a uniform grid $0\leq\zeta_{i}\leq 1$ either
$N=512$ or $N=1024$ points (the latter usually to check for convergence) and
all time steps were in the range by $0.01{\rm{St}}<\Delta t<0.025{\rm{St}}$.
And finally,
$\displaystyle\varepsilon_{{}_{0}}(t_{n})=\varepsilon_{{}_{00}}\exp\big{(}\beta
t_{n}\big{)}$.

Figs. 27 and 20 display a sample set of generated solutions of the horizontal
velocities and associated vorticities. A detailed discussion of these
solutions may be found in the main body of the text. Note, by the Howard
Semicircle Theorem, the magnitude of the azimuthal vorticity denote upper
bounds of the growth rates of shear induced instabilities when, of course, the
Richardson criterion is met.

Figure 27: Like Fig. 20; settling solutions for St = 0.04, $Z=0.01$, and
initial particle scale height $H_{p0}=0.05H_{g}$, at three selected times.

.

##  Appendix C A motivated single fluid model problem for the disk-analog of
the symmetric instability

We consider the linear normal mode response of a simplified single fluid “toy”
problem that exemplifies the essential processes we believe to be responsible
for the layer overturn during the bounce-phase of the simulations. As observed
in the simulations, settling dust imparts momentum upon the gas with a
vertical dependence. The resulting collective motion of the gas and dust shows
a net vertical variation in the azimuthal flow in the center of mass frame,
like depicted in Fig. 22 as well shown in Fig. 20 and Fig. 27. In the terminal
velocity and zero stopping time limit of the recast two-fluid equations, the
momentum exchange forcing will manifest as an effective non-zero vertically
varying external force acting in both the radial and azimuthal directions. In
the following examination we consider only the radial component of this force,
$f_{r}(z)$, as present and we say that it is not dynamically active. This
immediately leads to a new effective steady state azimuthal velocity
$\tilde{v}_{0}$ that adjusts to the force, i.e.,

$-2\Omega_{0}\tilde{v}_{0}(z)=f_{r}(z),$ (C1)

the analysis considered henceforth assumes single-fluid perturbations atop
this basic state.

We view the instability as being primarily axisymmetric, driven by the shear
in that part of the mean azimuthal flow that departs from the basic Keplerian
($=-(3/2)\Omega_{0}x$), which is driven by the background pressure gradient
and the particle stream, defined in the text as $\tilde{v}_{0}(z)$. We treat
buoyancy effects driven by the particle layer in the Boussinesq approximation,
namely that density fluctuations – whether they are driven by advective
motions or are a result of weak compressibility dependencies on temperature –
are dynamically significant only when coupled to the external gravitational
field (Spiegel & Veronis, 1960). The background density field of the model
single fluid is represented by $\rho_{0}(z)$. It is also assumed that the gas
fluid is incompressible. Thus the simple model for the perturbation quantities
$u^{\prime},v^{\prime},w^{\prime},\Delta^{\prime}$

$\displaystyle\partial_{t}u^{\prime}-2\Omega_{0}v^{\prime}$ $\displaystyle=$
$\displaystyle-(1/\rho_{00})\partial_{x}p^{\prime},$ (C2)
$\displaystyle\partial_{t}v^{\prime}+(1/2)\Omega_{0}u^{\prime}+w^{\prime}\partial_{z}\tilde{v}_{0}$
$\displaystyle=$ $\displaystyle 0,$ (C3) $\displaystyle\partial_{t}w^{\prime}$
$\displaystyle=$
$\displaystyle-(1/\rho_{00})\partial_{z}p^{\prime}-g(z)\Delta^{\prime},$ (C4)
$\displaystyle\partial_{t}\Delta^{\prime}+w^{\prime}\rho_{0}^{-1}\partial_{z}\rho_{0}$
$\displaystyle=$ $\displaystyle 0.$ (C5)

where $\rho_{00}$ is taken as a reference constant scale measure of the
density. These equations are supplemented with the statement of
incompressiblity,

$\partial_{x}u^{\prime}+\partial_{z}w^{\prime}=0.$ (C6)

The definition found in Eq. (C6) motivates the definition of a stream function
($\psi$) and vorticity ($\omega^{\prime}$) where

$u^{\prime}=\partial_{z}\psi^{\prime},\quad
w^{\prime}=-\partial_{x}\psi^{\prime},\qquad\omega^{\prime}\equiv\partial_{z}u^{\prime}-\partial_{x}w^{\prime}=\Big{(}\partial_{x}^{2}+\partial_{z}^{2}\Big{)}\psi^{\prime}.$
(C7)

The vertical component of gravity is taken to be

$g(z)=-\frac{\varepsilon_{0}}{1+\varepsilon_{0}}\Omega_{0}^{2}z$ (C8)

where we have included the factor $\varepsilon_{0}/(1+\varepsilon_{0})$, where
$\varepsilon_{0}$ is the midplane value of $\epsilon$, to represent an
effective reduced gravity in a single fluid formulation of the gas-particle
fluid (e.g. Chiang, 2008). We note that as per our observations in the
simulations, the mean gas density varies little and its dynamics are largely
incompressible, especially on the scales on which the particles cluster about
the midplane. We therefore represent the mean density layer profile and its
gradient via

$\rho_{0}^{-1}\partial_{z}\rho_{0}=\partial_{z}\ln\rho_{0}=-\frac{z}{H_{p}^{2}}.$
(C9)

The set of Eqs. (C2-C9) may be viewed as the linearized perturbations of the
dual terminal velocity (i.e., “strong-drag/small-grain”) and isothermal limit
of the one-fluid equivalent representation of a two-fluid system, as motivated
by both Laibe & Price (2014) and Lin & Youdin (2017). We note that this limit
contains no two-stream effects.

We rewrite all perturbation quantities as normal modes, e.g.,
$\psi^{\prime}\rightarrow\psi_{{}_{k}}(z)\exp\left[{ikx-i\omega_{0}t}\right]$,
and similarly for all the other quantities. We derive an equation for the
perturbation vorticity by operating on Eq. (C2) by $\partial_{z}$ and
subtracting from it the result of operating on Eq. (C4) by $-\partial_{x}$.
After inserting the assumed normal mode form into the result, as well as into
Eq. (C3) and Eq. (C5) we find the following three equations

$\displaystyle-i\omega_{0}\Big{(}\partial_{z}^{2}-k^{2}\Big{)}\psi_{{}_{k}}-2\Omega_{0}\partial_{z}v_{{}_{k}}$
$\displaystyle=$ $\displaystyle ikg(z)\Delta_{{}_{k}},$ (C10)
$\displaystyle-i\omega_{0}v_{{}_{k}}-ik\psi_{{}_{k}}\partial_{z}\tilde{v}_{0}+(\Omega_{0}/2)\partial_{z}\psi_{{}_{k}}$
$\displaystyle=$ $\displaystyle 0,$ (C11)
$\displaystyle-i\omega_{0}\Delta_{{}_{k}}-ik\psi_{{}_{k}}\cdot\partial_{z}\ln\rho_{0}$
$\displaystyle=$ $\displaystyle 0.$ (C12)

The above three equations may be further combined into a single equation for
$\psi_{k}$:

${\cal
L}\psi_{{}_{k}}=\Big{(}\Omega_{0}^{2}-\omega_{0}^{2}\Big{)}\partial_{z}^{2}\psi_{{}_{k}}+2\Omega_{0}ik\partial_{z}\Big{(}\psi_{{}_{k}}\partial_{z}\tilde{v}_{{}_{0}}\Big{)}+\omega_{0}^{2}k^{2}\psi_{{}_{k}}+k^{2}\psi_{{}_{k}}\Big{(}g\partial_{z}\ln\rho_{0}\Big{)}=0.$
(C13)

Our aim is to assess the growth rate, Im$(\omega_{0})$, and frequency
response, Re$(\omega_{0})$, to perturbuations to Eq. (C13) subject to the
condition that all perturbation quantities (i.e.,
$\psi_{{}_{k}},\partial_{z}\psi_{{}_{k}}$) decay as $z\rightarrow\pm\infty$.
To facilitate analytical treatment, we consider a model for the background
functions that permit us to seek solutions of Eq. (C13) in terms of parabolic
cylinder functions ${\cal D}_{m}(\zeta)$. We therefore adopt the following

$g\partial_{z}\ln\rho_{0}=-\frac{\varepsilon_{0}\Omega_{0}^{2}}{1+\varepsilon_{0}}\frac{z^{2}}{H_{p}^{2}},\qquad\tilde{v}_{{}_{0}}=v_{{\rm
cm},\infty}+\delta
v_{{}_{00}}\left(1-\frac{1}{2}\frac{z^{2}}{H_{s}^{2}}\right),$ (C14)

where $v_{{\rm cm},\infty}$ is the asymptotic constant value far from the
particle layer and $H_{s}$ is the scale factor associated with the vertical
variation of the azimuthal perturbation flow field as uncovered and discussed
in the text. Typically it is slightly larger than the particle scale height
$H_{p}$, but for our purposes here it is immaterial to the final outcome
below. It is important that the variation in $v$ be parabolic for what
follows. Thus, Eq. (C13) may be rewritten in terms of a standard parabolic
cylinder equation form if we introduce

$\psi_{k}=e^{-az^{2}/4}\Psi_{k};\qquad a=-\frac{2ik\delta
v_{{}_{00}}\Omega_{0}}{H_{s}^{2}(\omega_{0}^{2}-\Omega_{0}^{2})}.$ (C15)

and the vertically stretched coordinate $z=\tilde{\beta}\zeta$, where the
scale $\tilde{\beta}$, satisfies

$\frac{\tilde{\beta}^{4}\gamma^{2}}{\Omega_{0}^{2}-\omega_{0}^{2}}\left({\rm
Ri}-\frac{\Omega_{0}^{2}}{\Omega_{0}^{2}-\omega_{0}^{2}}\right)=\frac{1}{4};\qquad\gamma^{2}(k)\equiv\frac{\delta
v_{{}_{00}}^{2}k^{2}}{H_{s}^{4}}.$ (C16)

We note that $\tilde{\beta}$ may be complex and that we have re-expressed
constants in terms of the azimuthal form of the Richardson number, ${\rm
Ri}_{{\phi}}$ defined in Eq. (24), i.e.,

${\rm Ri}\longleftrightarrow{\rm
Ri}_{{\phi}}=\displaystyle{-\frac{\displaystyle{g\partial_{z}\ln\rho_{0}}}{\displaystyle{\big{(}\partial_{z}\tilde{v}_{0}}\big{)}^{2}}=\frac{\displaystyle{\frac{\varepsilon_{0}\Omega_{0}^{2}}{1+\varepsilon_{0}}\frac{z^{2}}{H_{p}^{2}}}}{\displaystyle{\frac{\delta
v_{{}_{00}}^{2}z^{2}}{H_{s}^{4}}}}}=\frac{\varepsilon_{0}}{1+\varepsilon_{0}}\frac{H_{s}^{4}}{H_{p}^{2}}\frac{\Omega_{0}^{2}}{\delta
v_{{}_{00}}^{2}}.$ (C17)

Eq. (C13) now simplifies into the canonical form in terms of the non-
dimensional independent variable $\zeta$

$\partial_{\zeta}^{2}\Psi_{k}+\left[\frac{\omega_{0}^{2}k^{2}-i\Omega_{0}\gamma}{\Omega_{0}^{2}-\omega_{0}^{2}}\tilde{\beta}^{2}-\frac{1}{4}\zeta^{2}\right]\Psi_{{}_{k}}=0.$
(C18)

The solution to the above are the standard Parabolic Cylinder functions ${\cal
D}_{m}(z/\tilde{\beta})$ (Abramowitz & Stegun, 1972), provided certain
conditions ensuring that either $\psi_{{}_{k}}$ rapidly decays as
$z\rightarrow\pm\infty$ or the quantization condition is satisfied: (i) first,
the quantization condition is

$\displaystyle\frac{\omega_{0}^{2}k^{2}-i\Omega_{0}\gamma}{\Omega_{0}^{2}-\omega_{0}^{2}}\tilde{\beta}^{2}=m+\frac{1}{2},$
(C19)

which must be satisfied for non-negative integer values of $m$; (ii) the
asymptotic behavior of $\psi_{k}$ for $|z|$ large is given by

$\psi_{k}\Big{(}|z|\rightarrow\infty\Big{)}\sim\exp\left[-\left(a+\frac{1}{\tilde{\beta}^{2}}\right)\frac{z^{2}}{4}\right],$
(C20)

which means that we require that

${\rm Re}\left(a+\frac{1}{\tilde{\beta}^{2}}\right)>0.$ (C21)

The quantization condition Eq. (C19) together with Eq. (C16) straightforwardly
leads to solutions for $\omega_{0}^{2}$ given by

$\displaystyle\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}$ $\displaystyle=$
$\displaystyle\sqrt{\mu\left(2\frac{{\rm Ri}-1}{{\rm
Ri}}+\mu-2i\frac{\sqrt{2\mu/\nu^{2}}}{\sqrt{{\rm
Ri}}}\right)}-\mu+i\frac{\sqrt{2\mu/\nu^{2}}}{\sqrt{{\rm Ri}}},$ (C22)

in which

$\nu\equiv
2m+1,\qquad\mu\equiv\frac{1}{2}\frac{\varepsilon_{0}}{1+\varepsilon_{0}}\frac{\nu^{2}}{k^{2}H_{p}^{2}}.$
(C23)

We have deliberately rewritten the above in order to express all vertical
velocity gradient information in terms of Ri.

The above form for $\omega_{0}^{2}$ satisfies the condition that solutions
exhibit Gaussian decay as $|z|\rightarrow\infty$. (The algebraic procedure
also admits solutions for $\omega_{0}^{2}$ where the sign of the first
(squareroot) term on the RHS of Eq. (C22) is negative. However, we find that
these solutions violate the large $|z|$ asymptotic decay criterion in both the
Ri $\rightarrow\infty$ and Ri $\rightarrow 0$ limits. There still remains the
possibility a second branch of solutions associated with this possible root
might be viable in some intermediate Ri limit, but this is not yet determined
at this stage.) The solutions are characterized by three parameters, the
vertical node $m$ (through $\nu$), the parameter $\mu$ which contains
information about the stratification as well as the horizontal wavenumber, and
Ri containing the relative measure of shear to stratification.

We note several features. Most importantly we see that $\omega_{0}=0$
solutions exist for specific values of ${\rm Ri}$. Setting $\omega_{0}$ to
zero in both Eqs. (C16) and (C19) reveals this is possible if both of the
following relationships are simultaneously satisfied

$\tilde{\beta}^{4}\gamma^{2}\big{(}{\rm Ri}-1\big{)}=\frac{1}{4},\qquad
i\tilde{\beta}^{2}\gamma=-\frac{\nu}{2}.$ (C24)

Upon eliminating $\tilde{\beta}$ in both expressions shows that this is
possible only when ${\rm Ri}={\rm Ri}_{c}$ where,

${\rm Ri}_{c}\equiv 1-\frac{1}{\nu^{2}}.$ (C25)

This condition signifies a strong change in the character of these solutions.
Furthermore, a series expansion of Eq. (C22) for $0<\mu\ll 1$ shows

$\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}\approx i\sqrt{\frac{2\mu}{{\rm
Ri}}}\left[\frac{1}{\nu}-\sqrt{1-{\rm Ri}}\right]+{\cal O}\left({\mu}\right),$
(C26)

which is written by taking into account the branch cut at ${\rm Ri}=1$ in the
first squareroot term on the RHS of Eq. (C22). The sign of the square
bracketed term changes sign when ${\rm Ri}$ crosses ${\rm Ri}_{c}$. In the
other limit where $\mu\rightarrow\infty$ we find

$\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}\approx\left(1-\frac{{\rm Ri}_{c}}{{\rm
Ri}}\right)\left(1+i\frac{1}{\nu}\sqrt{\frac{2}{\mu{\rm Ri}}}\right)+{\cal
O}\left({\frac{1}{\mu}}\right).$ (C27)

Once again, what stands out is the gross change in character when ${\rm Ri}$
passes ${\rm Ri}_{c}$. The large $\mu$ limit (high $m$ or small horizontal
wavenumbers $k$) shows that the leading order behavior of $\omega_{0}$ is that
of weakly growing oscillations with growth rates ${\cal
O}\left({\mu^{-1/2}}\right)$ when ${\rm Ri}>{\rm Ri}_{c}$ to very strong
${\cal O}\left({1}\right)$ widespread growth when ${\rm Ri}<{\rm Ri}_{c}$.
While ${\rm Ri}=1$ designates an important boundary indicating a stark
transition in character, we note that all normal modes show growth for all
non-infinite values of ${\rm Ri}$.

The limiting form for when the shear is nearly absent recovers at leading
order the known oscillatory behavior of disk intertia-gravity oscillations. A
series expansion of Eq. (C22) in inverse powers of Ri reveals that

$\frac{\omega_{0}^{2}}{\Omega_{0}^{2}}=\sqrt{\mu(2+\mu)}-\mu+\frac{i\sqrt{2}}{\nu\sqrt{{\rm
Ri}}}\left(\sqrt{\mu}-\mu\sqrt{\frac{1}{2+\mu}}\right)+{\cal
O}\left({\frac{1}{{\rm Ri}}}\right)$ (C28)

If we write
$\displaystyle\omega_{0}=\omega_{{}_{00}}+i\omega_{{}_{i,0}}+\cdots$ and
assume that $\omega_{{}_{i,0}}$ scales like ${\rm Ri}^{-1/2}$, then to leading
order we find

$\frac{\omega_{{}_{00}}}{\Omega_{0}}=\pm\sqrt{\sqrt{\mu(2+\mu)}-\mu},\qquad\frac{\omega_{{}_{i,0}}}{\Omega_{0}}=\frac{1}{\sqrt{2}\nu\omega_{00}}\left(\sqrt{\mu}-\mu\sqrt{\frac{1}{2+\mu}}\right)$
(C29)

The leading order expression for $\tilde{\beta}^{2}$ (denoted by
$\beta_{{}_{00}}^{2}$) is given as

$\frac{\tilde{\beta}_{{}_{00}}^{2}}{H_{p}^{2}}=\frac{1}{\nu\sqrt{2}}\left(1+\frac{1}{\varepsilon_{0}}\right)\sqrt{\left[\mu\left(1+\mu-\sqrt{\mu(2+\mu)}\right)\right]}.$
(C30)

We observe that the growth rate $\omega_{{}_{0,i}}$ has a $\mu$ dependence,
predicting growth for all values of $\mu$, with a corrsponding maximum value
at $\mu\equiv\mu_{{\rm max}}=2/3$. In this large Ri limit, we define the
maximum growth rate $\omega_{{}_{{\rm m}}}\equiv\omega_{{}_{i}}(\mu_{{}_{{\rm
max}}})$ and find that it is given by

$\frac{\omega_{{}_{{\rm
m}}}}{\Omega_{0}}=\frac{\sqrt{1+\frac{1}{\varepsilon_{0}}}}{2\sqrt{2}(2m+1)}\left(\frac{\delta
v_{{}_{00}}}{2\Omega_{0}H_{p}}\right)=\frac{1/\sqrt{{\rm
Ri}}}{2\sqrt{2}(2m+1)}.$ (C31)

We can similarly identify the fastest growing radial wavenumber based on the
definition of $\mu$ found in Eq. (C23). We define a corresponding radial
wavelength for the fastest growing mode ($\lambda_{{}_{{\rm m}}}$), related to
$k_{{\rm m}}=k(\mu=\mu_{{}_{{\rm max}}})$ via $k_{{}_{{\rm m}}}\equiv
2\pi/\lambda_{{}_{{\rm m}}}$. This leads to

$\displaystyle\frac{\lambda_{{}_{{\rm
m}}}}{H_{p}}=\frac{2\pi/\sqrt{3}}{m+\small\mbox{$\frac{1}{2}$}}\sqrt{1+\frac{1}{\varepsilon_{0}}}.$
(C32)

Finally the vertical scale scale characterizing variations in the parabolic
cylinder function is set by $\beta$ and we estimate its value at $k_{{}_{{\rm
m}}}$ based on $\beta$’s leading order behavior found in Eq. (C30). We
therefore say $\tilde{\beta}_{{}_{{\rm max}}}$ where

$\tilde{\beta}_{{}_{{\rm
m}}}\approx\tilde{\beta}_{{}_{00}}\Big{(}\mu=\mu_{{}_{{\rm
m}}}\Big{)}=\left(\frac{1}{6\sqrt{3}}\right)^{1/4}\left[\frac{\varepsilon_{0}+1}{\varepsilon_{0}(m+{1}/{2})}\right]^{1/2}H_{p}.$
(C33)

In order to help generate an analytic solution to discuss the physical content
of the simulations we have developed in the text, we have chosen a vertical
profile for $v$ that extends parabolically to infinity. However we know that
in the mean time-instantaneous profiles found in section 5.1 approach constant
profiles in that limit, but show a parabolic profile in the region around the
midplane. Furthermore, since the function ${\cal D}_{m}(z/\beta)$ show
Gaussian decay on a length scale defined by $\beta$, it is reasonable to
suppose that these analytical solutions are representative of the more
realistic circumstances found section 5.1 so long as the decay scale is less
than the particle scale height, i.e., $\beta<H_{p}$. Additionally, the first
few zeros of ${\cal D}_{m}(\zeta)$ for given values of $m$, occur for values
of $\zeta=\zeta_{0i}<1$. This feature ensures that the modes have multiple
nodes across the vertical extent. Based on this reasoning, we expect that the
fastest growing of the analytic modes we have developed here are
representative so long as $\tilde{\beta}_{{}_{{\rm m}}}<1$. Thus by Eq. (C33)
we require

$\left(\frac{1}{6\sqrt{3}}\right)^{1/4}\left(\frac{\varepsilon_{0}+1}{\varepsilon_{0}(m+{1}/{2})}\right)^{1/2}\leq
1,$ (C34)

and this expression places the following constraint on $m$,

$m+\frac{1}{2}>\left(\frac{1}{6\sqrt{3}}\right)^{1/2}\left(\frac{1+\varepsilon_{0}}{\varepsilon_{0}}\right),$
(C35)

in order to plausibly apply these approximate – though analytically derived –
solutions to interpret the flow transitions considered in this study.

Let $\pm\zeta_{{}_{m,i}}$ denote the $i^{{\rm th}}$ zero of ${\cal
D}_{m}(\zeta)$. For each order $m$ there are $m$ such zeros which are also
those of the the Hermite polynomials ${\cal H}_{m}(\zeta)$. The values of
these are: $\zeta_{{}_{2,\pm 1}}=\pm\sqrt{2}/2$, $\zeta_{{}_{3,i=\\{0,\pm
1\\}}}\approx 0,\pm 1.22474$, $\zeta_{{}_{4,i=\\{\pm 1,\pm 2\\}}}\approx\pm
0.52464,\pm 1.65068$, $\zeta_{{}_{4,i=\\{0,\pm 1,\pm 2\\}}}\approx 0,\pm
0.95857,\pm 2.02018$ and so on (Abramowitz & Stegun, 1972). Because of the
vertical scaling associated with the fastest growing mode, $\beta_{m}$ defined
in Eq. (C33), we have therefore the distance between successive nodes for
solution index $m$ to be

$\Delta z_{{}_{m,i}}=\tilde{\beta}_{{}_{\rm
m}}\Big{(}\zeta_{{}_{m,i}}-\zeta_{{}_{m,i-1}}\Big{)}.$ (C36)

Due to symmetry considerations it is enough to quantify $\Delta z_{{}_{m,i}}$
by restricting consideration to values of $1\leq i\leq m/2$ for $m$ even, and
$1\leq i\leq(m-1)/2$ for $m$ odd.

##  Appendix D Hyperdiffusion and hyperviscosity scheme in PENCIL code

The flux of viscous momentum can be cast to be proportional to the rate of
strain tensor, which can be written as

$S_{mn}=\frac{1}{2}\left(\frac{\partial u_{m}}{\partial x_{n}}+\frac{\partial
u_{n}}{\partial x_{m}}-\frac{1}{3}\delta_{mn}\nabla\cdot\bf{{\bf
u}_{g}}\right),$ (D1)

which, in the incompressible limit, along with constant dynamical viscosity
$\mu=\nu\rho$, translates to a viscous acceleration as $\nu\nabla^{2}{\bf
u}_{g}$, where $\nu$ is the kinematic viscosity. Note that, in Eq. (D1), $p$
and $q$ are used as dummy indices in order to denote the gas velocity
components. The viscous acceleration $f_{{\text{visc}}}$, in its most general
form can be written as

$f_{{\text{visc}}}=\frac{1}{\rho}\nabla\cdot\left[2\rho\nu\bm{S}\right].$ (D2)

With the PENCIL code, however, we are using a higher order (sixth) hyper-
dissipation, and hence, the rate-of-strain tensor gets replaced by a higher-
order version as below:

$f_{{\text{visc}}}^{({\text{hyper})}}=\frac{1}{\rho}\nabla\cdot\left[2\rho\nu_{n}{\bm{S}}^{(n)}\right].$
(D3)

Here, $n=3$ which corresponds to the sixth order hyperdissipation. With this,
a simple hyperviscosity is applied in the following form:

$f_{{\text{visc}}}=\begin{cases}\nu_{3}\nabla^{6}{\bf u}_{g},&\text{if $\mu$ =
constant;}\\\ \nu_{3}\left[\nabla^{6}{\bf
u}_{g}+2\bm{S}^{(3)}\cdot\nabla\ln\rho_{g}\right],&\text{if $\nu$ =
constant.}\end{cases}$ (D4)

The spectral range over which the hyperviscosity operates is small, and its
used only as a high frequency filter in the numerical setup. Moreover, given
its artificial nature, the strict requirement of momentum conservation is
dropped (Lyra et al., 2017, also see PENCIL code manual). It is important to
note that the simplified expression from Eq. (D4) can be written as the
divergence of a rate-of-strain tensor as

$S_{mn}^{(3)}=\frac{\partial^{5}u_{m}}{\partial x_{n}^{5}}.$ (D5)

For more details on the hyperdissipation scheme used in the code, see the
PENCIL CODE manual and Lyra et al. (2017).

## References

  * Abramowitz & Stegun (1972) Abramowitz, M. & Stegun, I. A. 1972, Handbook of Mathematical Functions, New York: Dover, 1972
  * Abarbanel et al. (1984) Abarbanel, H. D. I., Holm, D. D., Marsden, J. E., et al. 1984, Phys. Rev. Lett., 52, 2352. doi:10.1103/PhysRevLett.52.2352
  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192. doi:10.3847/1538-4357/ab40a3
  * Alexakis & Biferale (2018) Alexakis, A. & Biferale, L. 2018, Phys. Rep., 767, 1. doi:10.1016/j.physrep.2018.08.001
  * Barranco (2009) Barranco, J.-A. 2009, ApJ, 691, 907
  * Bennetts & Hoskins (1979) Bennetts, D. A. & Hoskins, B. J. 1979, Quarterly Journal of the Royal Meteorological Society, 105, 945. doi:10.1002/qj.49710544615
  * Brandenburg (2001) Brandenburg, A. 2001, ApJ, 550, 824. doi:10.1086/319783
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Chandrasekhar (1961) Chandrasekhar, S. 1961, International Series of Monographs on Physics, Oxford: Clarendon, 1961
  * Chen & Lin (2020) Chen, K. & Lin, M.-K. 2020, ApJ, 891, 132. doi:10.3847/1538-4357/ab76ca
  * Chambers (2010) Chambers, J. E. 2010, Icarus, 208, 505. doi:10.1016/j.icarus.2010.03.004
  * Chiang (2008) Chiang, E. 2008, ApJ, 675, 1549
  * Cox & Matthews (2002) Cox, S. M. & Matthews, P. C. 2002, Journal of Computational Physics, 176, 430. doi:10.1006/jcph.2002.6995
  * Cuzzi et al. (1993) Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
  * Davidson (2004) Davidson, P. A. 2004, Turbulence : an introduction for scientists and engineers
  * Dobrovolskis et al. (1999) Dobrovolskis, A. R., Dacles-Mariani, J. S., & Cuzzi, J. N. 1999, J. Geophys. Res., 104, 30805
  * Drazin & Reid (2004) Drazin, P. G. & Reid, W. H. 2004, Hydrodynamic Stability, by P. G. Drazin and W. H. Reid, pp. 626. ISBN 0521525411. Cambridge, UK: Cambridge University Press, September 2004., 626
  * Drazkowska et al. (2022) Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2022, arXiv:2203.09759
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  * Estrada et al. (2016) Estrada, P. R., Cuzzi, J. N., & Morgan, D. A. 2016, ApJ, 818, 200
  * Garaud & Lin (2004) Garaud, P. & Lin, D. N. C. 2004, ApJ, 608, 1050. doi:10.1086/420839
  * Gerbig et al. (2020) Gerbig, K., Murray-Clay, R. A., Klahr, H., et al. 2020, ApJ, 895, 91. doi:10.3847/1538-4357/ab8d37
  * Goldreich & Lynden-Bell (1965) Goldreich, P. & Lynden-Bell, D. 1965, MNRAS, 130, 125. doi:10.1093/mnras/130.2.125
  * Gómez & Ostriker (2005) Gómez, G. C. & Ostriker, E. C. 2005, ApJ, 630, 1093. doi:10.1086/432086
  * Hasegawa & Tsuribe (2014) Hasegawa, Y. & Tsuribe, T. 2014, PASJ, 66, L2. doi:10.1093/pasj/psu013
  * Haugen & Brandenburg (2004) Haugen, N. E. L. & Brandenburg, A. 2004, Phys. Rev. E, 70, 026405. doi:10.1103/PhysRevE.70.026405
  * Hockney & Eastwood (1981) Hockney, R. W. & Eastwood, J. W. 1981, Computer Simulation Using Particles, New York: McGraw-Hill, 1981
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., et al. 2020, ApJ, 904, 132. doi:10.3847/1538-4357/abc334
  * Hartlep & Cuzzi (2020) Hartlep, T. & Cuzzi, J. N. 2020, ApJ, 892, 120. doi:10.3847/1538-4357/ab76c3
  * Hoskins (1974) Hoskins, B. J. 1974, Quarterly Journal of the Royal Meteorological Society, 100, 480. doi:10.1002/qj.49710042520
  * Howard (1961) Howard, L. N. 1961, Journal of Fluid Mechanics, 10, 509. doi:10.1017/S0022112061000317
  * Ishitsu et al. (2009) Ishitsu, N., Inutsuka, S.-.I., & Sekiya, M. 2009, arXiv:0905.4404
  * Johansen et al. (2006) Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219. doi:10.1086/502968
  * Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022
  * Kraichnan & Montgomery (1980) Kraichnan, R. H. & Montgomery, D. 1980, Reports on Progress in Physics, 43, 547. doi:10.1088/0034-4885/43/5/001
  * Krapp et al. (2020) Krapp, L., Youdin, A. N., Kratter, K. M., et al. 2020, MNRAS, 497, 2715. doi:10.1093/mnras/staa1854
  * Laibe & Price (2014) Laibe, G. & Price, D. J. 2014, MNRAS, 440, 2136. doi:10.1093/mnras/stu355
  * Laibe et al. (2020) Laibe, G., Bréhier, C.-E., & Lombart, M. 2020, MNRAS, 494, 5134. doi:10.1093/mnras/staa994
  * Latter & Papaloizou (2017) Latter, H. N. & Papaloizou, J. 2017, MNRAS, 472, 1432. doi:10.1093/mnras/stx2038
  * Lee et al. (2010a) Lee, A. T., Chiang, E., Asay-Davis, X., et al. 2010, ApJ, 718, 1367. doi:10.1088/0004-637X/718/2/1367
  * Lee et al. (2010b) Lee, A. T., Chiang, E., Asay-Davis, X., et al. 2010, ApJ, 725, 1938. doi:10.1088/0004-637X/725/2/1938
  * Li et al. (2018) Li, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14
  * Li & Youdin (2021) Li, R. & Youdin, A. N. 2021, ApJ, 919, 107. doi:10.3847/1538-4357/ac0e9f
  * Lin & Youdin (2017) Lin, M.-K., & Youdin, A. N. 2017, ApJ, 849, 129
  * Lin (2021) Lin, M.-K. 2021, ApJ, 907, 64. doi:10.3847/1538-4357/abcd9b
  * Lyra et al. (2017) Lyra, W., McNally, C. P., Heinemann, T., et al. 2017, AJ, 154, 146. doi:10.3847/1538-3881/aa8811
  * Lyra et al. (2018) Lyra, W., Raettig, N., & Klahr, H. 2018, Research Notes of the American Astronomical Society, 2, 195. doi:10.3847/2515-5172/aaeac9
  * Lyra, & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001
  * Matthaeus et al. (1991) Matthaeus, W. H., Stribling, W. T., Martinez, D., et al. 1991, Phys. Rev. Lett., 66, 2731. doi:10.1103/PhysRevLett.66.2731
  * Miles (1961) Miles, J. W. 1961, Journal of Fluid Mechanics, 10, 496. doi:10.1017/S0022112061000305
  * Miles (1986) Miles, J. 1986, Physics of Fluids, 29, 3470. doi:10.1063/1.865812
  * Mkhinini et al. (2013) Mkhinini, N., Dubos, T., & Drobinski, P. 2013, Journal of Fluid Mechanics, 728, 29. doi:10.1017/jfm.2013.250
  * Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Nesvorný et al. (2019) Nesvorný, D., Li, R., Youdin, A. N., et al. 2019, Nature Astronomy, 3, 808
  * Nimmo et al. (2017) Nimmo, F., Umurhan, O., Lisse, C. M., et al. 2017, Icarus, 287, 12. doi:10.1016/j.icarus.2016.06.027
  * Pandey et al. (2019) Pandey, V., Perlekar, P., & Mitra, D. 2019, Phys. Rev. E, 100, 013114. doi:10.1103/PhysRevE.100.013114
  * Raettig et al. (2021) Raettig, N., Lyra, W., & Klahr, H. 2021, ApJ, 913, 92. doi:10.3847/1538-4357/abf739
  * Regev et al. (2016) Regev, O., Umurhan, O. M., & Yecko, P. A. 2016, Modern Fluid Dynamics for Physics and Astrophysics. Series: Graduate Texts in Physics, ISBN: 978-1-4939-3163-7. Springer New York
  * Richard et al. (2016) Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571
  * Schäfer et al. (2020) Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190. doi:10.1051/0004-6361/201937371
  * Sekiya (1998) Sekiya, M. 1998, Icarus, 133, 298
  * Sekiya & Onishi (2018) Sekiya, M. & Onishi, I. K. 2018, ApJ, 860, 140. doi:10.3847/1538-4357/aac4a7
  * Sengupta et al. (2019) Sengupta, D., Dodson-Robinson, S. E., Hasegawa, Y., et al. 2019, ApJ, 874, 26. doi:10.3847/1538-4357/aafc36
  * Simon et al. (2016) Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55
  * Simon et al. (2017) Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, 847, L12
  * Spiegel & Veronis (1960) Spiegel, E. A. & Veronis, G. 1960, ApJ, 131, 442
  * Squire & Hopkins (2018a) Squire, J., & Hopkins, P. F. 2018, ApJ, 856, L15
  * Squire & Hopkins (2018b) Squire, J., & Hopkins, P. F. 2018, MNRAS, 477, 5011
  * Stamper & Taylor (2017) Stamper, M. A. & Taylor, J. R. 2017, Ocean Dynamics, 67, 65. doi:10.1007/s10236-016-1011-6
  * Stone (1966) Stone, P. H. 1966, Journal of Atmospheric Sciences, 23, 390. doi:10.1175/1520-0469(1966)023¡0390:ONGBS¿2.0.CO;2
  * Thomas et al. (2013) Thomas, L. N., Taylor, J. R., Ferrari, R., et al. 2013, Deep Sea Research Part II: Topical Studies in Oceanography, 91, 96. doi:10.1016/j.dsr2.2013.02.025
  * Umurhan & Regev (2004) Umurhan, O. M., & Regev, O. 2004, A&A, 427, 855
  * Umurhan et al. (2007) Umurhan, O. M., Menou, K., & Regev, O. 2007, Phys. Rev. Lett., 98, 034501. doi:10.1103/PhysRevLett.98.034501
  * Umurhan et al. (2020) Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4. doi:10.3847/1538-4357/ab899d
  * Youdin & Johansen (2007) Youdin, A. & Johansen, A. 2007, ApJ, 662, 613. doi:10.1086/516729
  * Vanneste (1993) Vanneste, J. 1993, The Kelvin-Helmholtz instability in a nongeostrophic baroclinic unstable flow. Math Comput Model 17(1):149–154
  * Verma & Kumar (2015) Verma, M. & Kumar, A. 2015, APS Division of Fluid Dynamics Meeting Abstracts
  * Weidenschilling (1980) Weidenschilling, S. J. 1980, Icarus, 44, 172
  * Yang et al. (2017) Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27
  * Yellin-Bergovoy et al. (2021) Yellin-Bergovoy, R., Umurhan, O. M., & Heifetz, E. 2021, Geophysical & Astrophysical Fluid Dynamics, 895, 4. DOI: 10.1080/03091929.2021.1941921
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, A&A, 620, 459 (YG2005)
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
  * Zeitlin (2018) Zeitlin, V. 2018, Physics of Fluids, 30, 061701. doi:10.1063/1.5031099
  * Zhou et al. (2022) Zhou, H., Dewar, W., Yang, W., et al. 2022, Communications Earth and Environment, 3, 28. doi:10.1038/s43247-022-00362-4