# Primordial Black Holes And Gravitational Waves Based On No-Scale
Supergravity

Ioanna D. Stamou  National and Kapodistrian University of Athens, Department
of Physics, Section of Nuclear & Particle Physics, GR–15784 Athens, Greece
<EMAIL_ADDRESS>
###### Abstract

In this paper we present a class of models in order to explain the production
of Primordial Black Holes (PBHs) and Gravitational Waves (GWs) in the
Universe. These models are based on no-scale theory. By breaking the
SU(2,1)/SU(2)$\times$U(1) symmetry we fix one of the two chiral fields and we
derive effective scalar potentials which are capable of generating PBHs and
GWs. As it is known in the literature there is an important unification of the
no-scale models, which leads to the Starobinsky effective scalar potential
based on the coset SU(2,1)/SU(2)$\times$U(1). We use this unification in order
to additionally explain the generation of PBHs and GWs. Concretely, we modify
well-known superpotentials, which reduce to Starobinsky-like effective scalar
potentials. Thus, we derive scalar potentials which, on the one hand, explain
the production of PBHs and GWs and, on the other hand, they conserve the
transformation laws, which yield from the parametrization of the coset
SU(2,1)/SU(2)$\times$U(1) as well as the unification between the models which
are yielded this coset. We numerically evaluate the scalar power spectra with
the effective scalar potential based on this theory. Furthermore, we evaluate
the fractional abundances of PBHs by comparing two methods the Press–Schechter
approach and the peak theory, while focusing on explaining the dark matter in
the Universe. By using the resulting scalar power spectrum we evaluate the
amount of GWs. All models are in complete consistence with Planck constraints.

##  1 Introduction

The detection of Gravitational Waves (GWs) by a binary black hole merge,
reported by LIGO/VIRGO collaboration [1, 2, 3], opens a new window in physics
of Primordial Black Holes (PBHs). As a result there are numerous recent
studies which show that the origin of PBHs can explain a fraction of Dark
Matter (DM) in the Universe. Moreover as the detection of GWs is regarded as a
milestone in Cosmology, there are signals of GWs that are expected to be
detected by future space-based interferometer, such as LISA, BBO and DECIGO. A
well-known way, in order to explain both the generation of PBHs and GWs, is
the models embedded in inflation. Specifically, an enhancement in scalar power
spectrum can be interpreted as a significant amount of both PBHs and GWs,
which are formed in the radiation dominated era with the origin to the
previous epoch of inflation.

There are numerous theoretical models, proposed in the literature for
explaining the generation of PBHs and GWs. All these models have to deal with
the observable constraints released by Planck collaboration[4] for the
prediction of spectral index $n_{s}$ and tensor-to-scalar ratio $r$. A model
which is in complete accordance with this observable bounds is the Starobinsky
potential. A theory which can lead to Starobinsky potential is the no-scale
supergravity [5, 6, 7]. For this reason we choose to work in the framework of
this theory.

Unfortunately, the Starobinsky potential cannot lead to an enhancement of the
scalar power spectrum. In the literature many mechanisms have been proposed in
order to achieve this amplification. One of the most discussed is an
inflection point in the framework of single field inflation. This inflection
point can be given as a feature in the effective scalar potential where the
first derivative and the second derivative of the potential with respect to
the field are almost zero. This mechanism is regarded as a simple way which
gives an enhancement in scalar power spectrum. However, a disadvantage of
potentials with an inflection point is that a lot of fine-tuning of the
parameters is required. In our analysis we adopt this mechanism and we analyze
the issue of fine-tuning.

This paper is organized as follows. In the section 2 we present our theory,
which is based on no-scale supergravity, in order to achieve an inflection
point in the effective scalar potential. In section 3 we show how our
modifications lead to a significant enhancement in the scalar power spectrum.
In section 4 we show how this enhancement can be interpreted as the production
of PBHs. In section 5 we calculate the energy density of induced GWs. In
section 6 we discuss about the issue of the fine-tuning and in section 7 we
present some conclusions of our work

##  2 Theory

As the theory of our analysis is based on the SU(2,1)/SU(2)$\times$U(1)
symmetry, we briefly refer to some basic aspects about symmetries in the
kinetic term in the Lagrangian. In $N=1$ supergravity the general Lagrangian
in effective field theory is given as:

$\mathcal{L}=K^{\bar{j}}_{i}\partial_{\mu}\Phi^{i}\partial^{\mu}\bar{\Phi}_{\bar{j}}-V(\Phi)$
(1)

where $\Phi_{i}$ are the chiral fields, $\bar{\Phi_{i}}$ is its conjugate and
subscripts $i,j$ running over the fields. With $K^{\bar{j}}_{i}$ we denote the
Kähler metric and $K$ is the Kähler potential. We work in Planck units
throughout this work, $M_{P}=1$. The F-term of scalar potential is given as
follows:

$V=e^{K}(D_{\Phi}WK^{\bar{\Phi}\Phi}D_{\bar{\Phi}}\bar{W}-3|W|^{2})$ (2)

where $W$ is the superpotential and the Kähler covariant derivative is given:

$D_{\Phi}W=\frac{\partial W}{\partial\Phi}+\frac{\partial K}{\partial\Phi}W.$

In order to vanish the cosmological constant, or in other words in order to
vanish the term $3|W|^{2}$ in the effective scalar potential, Eq. (2), we need
to have the following identity:

$K^{\Phi\bar{\Phi}}K_{\Phi}K_{\bar{\Phi}}=3.$ (3)

There is an elegant way to obtain this identity by considering the Kähler
potential in this form [8]:

$K=-3\ln(\Phi+\bar{\Phi}).$ (4)

This Kähler potential represent the minimal no-scale model written in term of
one chiral field. The corresponding kinetic term of the Lagrangian is:

$\mathcal{L}_{kin}=\frac{3}{(\Phi+\bar{\Phi})}\partial^{\mu}\Phi\partial_{\mu}\bar{\Phi}$.

In order to discuss the symmetries associated with the kinetic term of the
Langrangian we consider that the field transforms as $\Phi=(y+1)/(y-1)$. So
the Kähler potential in the term of the new fiels $y$ is[9]:

$K=-3\ln(1-\frac{|y|^{2}}{3}).$ (5)

In this form the corresponding kinetic term of the Lagrangian is given from:

$\mathcal{L}_{kin}=\frac{3}{(1-|y|^{2})^{2}}\partial^{\mu}y\partial_{\mu}\bar{y}$

which is invariant under the transformation of:

$y\rightarrow\frac{\alpha
y+\beta}{\bar{\beta}y+\bar{\alpha}},\quad|\alpha|^{2}-|\beta|^{2}=1.$ (6)

This invariance defines the non-compact group SU(1,1). Therefore the Kähler
potentials written in $T$ basis, from Eq. (4), and in $y$ basis, from Eq. (5),
are two forms of the non-minimal no-scale SU(1,1) model. These equivalent
forms of the Kähler potential do not lead to Starobinsky-like effective scalar
potential. So we need to consider an extension of this group, which is the
non-compact SU(2,1)/SU(2)$\times$U(1) coset. Under this coset we consider two
chiral fields, one plays the role of the inflaton and the other is the modulo
field.

As in the case of SU(1,1), there are two equivalent form of
SU(2,1)/SU(2)$\times$U(1) symmetry, which read as:

$K=-3\ln(1-\frac{|y_{1}|^{2}}{3}-\frac{|y_{2}|^{2}}{3})\quad\&\quad
K=-3\ln(T+\bar{T}-\frac{|\varphi|^{2}}{3}).$ (7)

The complex fields $(y_{1},y_{2})$ are related to $(T,\varphi)$ by the
following expressions:

$y_{1}=\Big{(}\frac{2\varphi}{1+2T}\Big{)},\quad
y_{2}=\sqrt{3}\Big{(}\frac{1-2T}{1+2T}\Big{)}$

and the inverse relations by:

$T=\frac{1}{2}\Big{(}\frac{1-y_{2}/\sqrt{3}}{1+y_{2}/\sqrt{3}}\Big{)},\quad\varphi=\Big{(}\frac{y_{1}}{1+y_{2}/\sqrt{3}}\Big{)}.$

The superpotential transforms non-trivial between these two basis.
Specifically we consider:

$W(T,\varphi)\rightarrow\bar{W}(y_{1},y_{2})=(1+y_{2}/\sqrt{3})^{3}W.$

In order to find Starobinsky-like effective scalar potential we consider two
forms of superpotentials which are written first in the more revealing basis
$(y_{1},y_{2})$ and then in basis $(T,\varphi)$:

$W_{WZ}=\Big{(}\frac{\hat{\mu}}{2}\Big{(}y_{1}^{2}+\frac{y_{1}^{2}y_{2}}{\sqrt{3}}\Big{)}-\lambda\frac{y_{1}^{3}}{3}\Big{)}\iff
W^{\prime}_{WZ}=\frac{\hat{\mu}}{2}\varphi^{2}-\frac{\lambda}{3}\varphi^{3}$
(8)

$W_{C}=m\Big{(}-y_{1}y_{2}+\frac{y_{2}y_{1}^{2}}{l\sqrt{3}}\Big{)}\iff
W^{\prime}_{C}=\sqrt{3}m\varphi\left(T-\frac{1}{2}\right).$ (9)

In the $(T,\varphi)$ basis one can recognize that Eq. (8) is the Wess-Zumino
superpotential and the Eq. (9) is the Cecotti one [5, 6]. In order to derive
Starobinsky-like effective scalar potential from the above superpotials, we
consider that the inflaton field is the $y_{1}$ and the modulo the $y_{2}$. If
we chance basis the role of inflaton for the $W^{\prime}_{WZ}$ is the
$\varphi$ and the corresponding for $W^{\prime}_{C}$ is the $T$. Moreover, we
should consider the following choices $\lambda/\mu=1/3$ and
$\hat{\mu}=\mu/\sqrt{3}$ for $W_{WZ}$ and $l=1$ for $W_{C}$. The
superpotentials in $(y_{1},y_{2})$ basis keep the transformations laws:

$y_{1}\rightarrow\alpha y_{1}+\beta y_{2},\quad
y_{2}\rightarrow-\beta^{*}y_{1}+\alpha^{*}y_{2}.$ (10)

where $\alpha,\beta$ are complex numbers and $|\alpha|^{2}+|\beta|^{2}=1$ [6,
7].

Figure 1: Figure adapted from Ref. [6].

However, Starobinsky-like effective scalar potential cannot explain the
generation of PBHs and GWs. As the detection of GWs by LIGO/VIRGO and the
future space based experiments for GWs are on the topic of the cosmology we
consider the following two superpotentials, which come from a modification of
Eqs. (8) and (9)[10]:

$W_{1}=\Big{(}\frac{\hat{\mu}}{2}\Big{(}y_{1}^{2}+\frac{y_{1}^{2}y_{2}}{\sqrt{3}}\Big{)}-\lambda\frac{y_{1}^{3}}{3}\Big{)}(1+e^{-b_{1}y_{1}^{2}}(c_{1}y_{1}^{2}+c_{2}y_{1}^{4}))$
(11)

$W_{2}=m\Big{(}-y_{1}y_{2}+\frac{y_{2}y_{1}^{2}}{l\sqrt{3}}\Big{)}(1+c_{3}e^{-b_{2}y_{1}^{2}}y_{1}^{2}).$
(12)

These superpotentials are written in the $(y_{1},y_{2})$ basis. In the
($T,\varphi$) they take the form:

$W_{1}^{\prime}=\left(\frac{\hat{\mu}}{2}\varphi^{2}-\frac{\lambda}{3}\varphi^{3}\right)\left(1+e^{-b_{1}\left(\frac{2\varphi}{1+2T}\right)^{2}}\left[c_{1}\left(\frac{2\varphi}{1+2T}\right)^{2}+c_{2}\left(\frac{2\varphi}{1+2T}\right)^{4}\right]\right)$
(13)

$W_{2}^{\prime}=\frac{\sqrt{3}}{2}m\varphi\Big{(}\frac{1}{2}-T\Big{)}\left(-1-2T-\frac{1-2T}{l}\right)\left[1+3c_{3}e^{-3b_{2}\left(\frac{1-2T}{1+2T}\right)^{2}}\left(\frac{1-2T}{1+2T}\right)^{2}\right].$
(14)

In figure 1 we notice that starting by the Kähler potential in $(y_{1},y_{2})$
basis, the effective-scalar potential can keep its form under the field
transformation colored with blue arrows. In order to find the effective scalar
potential we need to set which is the inflaton field and which is the modulo
and we need fix the non canonical kinetic term for the inflaton by a further
redefinition of the field, which is depicted in this figure with black inner
arrows. We can obtain the same effective effective scalar potential throughout
the cycle in the left side of this figure by interplay the role of inflaton
and modulo field according to this figure. Moreover, we can obtain the same
effective scalar potential if we move to the other basis $(T,\varphi)$ and
consider the other form of Kähler potential and the corresponding forms of
superpotential from Eqs. (13) and (14). This study and the importance of the
unification between the models have previously analyzed for the case of
Starobinsky scalar potential by considering Eqs. (8) and (9)[5, 6, 7]. However
these two modified superpotentials (11) and (12) (or their equivalent forms
(13) and (14)) can explain the generation of PBHs and GWs due to the
inflection point, which is presented in their effective scalar potential.

##  3 Evaluating Scalar Power Spectrum

In this section we present our analysis in order to derive numerically the
scalar power spectrum. The equation of motion for the inflaton field is given
in efold times by the expression:

$\chi^{\prime\prime}+3\chi^{\prime}-\frac{1}{2}\chi^{\prime
3}+\left(3-\frac{1}{2}\chi^{\prime 2}\right)\frac{d\ln V(\chi)}{d\chi}=0$ (15)

where with primes we denote the derivative in efold time and as the equations,
which follow, are valid for canonical kinetic term in the Lagrangian. In
figure 1, one can notice how the canonical $\chi$ field is connected to the
$(y_{1},y_{2})$ fields and the $(T,\varphi)$ each time.

If we consider that the perturbation of the field is given as
$\chi+\delta\chi$ we have for the field’s perturbation the following
expression:

$\delta\chi^{\prime\prime}=-\left(3-\frac{1}{2}\chi^{\prime
2}\right)\delta\chi^{\prime}-\frac{1}{H^{2}}\frac{d^{2}V}{d\chi^{2}}\delta\chi-\frac{k^{2}}{a^{2}H^{2}}\delta\chi+4\Psi^{\prime}\chi^{\prime}-\frac{2\Psi}{H^{2}}\frac{dV}{d\chi}$
(16)

and for the Bardeen potential $\Psi$ we have the equation:

$\Psi^{\prime\prime}=-\left(7-\frac{1}{2}\chi^{\prime
2}\right)\Psi^{\prime}-\left(2\frac{V}{H^{2}}+\frac{k^{2}}{a^{2}H^{2}}\right)\Psi-\frac{1}{H^{2}}\frac{dV}{d\chi}\delta\chi$
(17)

where $k$ is the comoving wavenumber and $H$ is the Hubble parameter.

After integrating Eqs.(15), (16) and (17) for each comoving wavenumber $k$, we
can evaluate the scalar power spectrum from the following expression:

$P_{R}=\frac{k^{3}}{2\pi^{2}}\left|R_{k}\right|^{2},$ (18)

where $R_{k}$ is the comoving curvature perturbation:

$R_{k}=\Psi+\frac{\delta\phi}{\phi^{\prime}}.$ (19)

An analysis of the numerical treatment can be found in Ref.[11, 12].

Figure 2:  Scalar power spectrum for the superpotenrial $W_{1}$ and $W_{2}$
(or $W^{\prime}_{1}$ and $W^{\prime}_{2}$). Parameters: $c_{1}=1.7$,
$c_{2}=-0.65$, $c_{3}=14$, $b_{1}=3.701266$, $b_{2}=8.995722$, $l=1.0002$,
$\lambda/\mu=1/3$.

In figure 2 we show the scalar power spectrum for the superpotentials $W_{1}$
and $W_{2}$, given from Eqs. (11) and (12), with red and blue respectively. We
can notice that we have an significant enhancement of power at around
$10^{13}$ $Mpc^{-1}$. This enhancement can be interpreted as the production of
PBH and GWs in radiation dominated epoch with the origin to the previous epoch
of inflation. Moreover it can be interpreted as a significant fraction of the
DM in the Universe. The results for their equivalent forms $W_{1}^{\prime}$
and $W_{2}^{\prime}$, given from Eqs. (13) and (14) are the same. For initial
condtions of the field we choose $\chi_{ic}=4.95$ for $W_{1}$ and
$\chi_{ic}=4.87$ for $W_{2}$.

##  4 Production of PBHs

The present abundance of PBHs is given by the integral:

$f_{PBH}=\int d\ln M\frac{\Omega_{PBH}}{\Omega_{DM}}$ (20)

where

$\frac{\Omega_{PBH}}{\Omega_{DM}}=\frac{\beta(M(k))}{8\times
10^{-16}}\left(\frac{\gamma}{0.2}\right)^{3/2}\left(\frac{g(T_{f})}{106.75}\right)^{-1/4}\left(\frac{M(k)}{10^{-18}\textsl{g}}\right)^{-1/2}.$
(21)

The mass is given as a function of k mode:

$M_{PBH}=10^{18}\left(\frac{\gamma}{0.2}\right)\left(\frac{g}{106.75}\right)^{-1/6}\left(\frac{k}{7\times
10^{13}Mpc^{-1}}\right)^{-2}\textsl{g}$ (22)

where $g$ is the effective number of degree of freedom and $\gamma$ is the
correction factor which depends on the gravitational collapse. The mass
fraction of the Universe collapse in PBH of $M_{PBH}$, denoted as $\beta$, is
calculated with two methods, the Press-Schecter approach (PS) and the peak
theory (PT).

In the PS approach the mass fraction $\beta_{PS}$ is given by the probability
that the overdensity $\delta$ is above a certain threshold of collapse,
denoted as $\delta_{c}$:

$\beta_{PS}(M_{PBH})=\frac{1}{\sqrt{2\pi\sigma^{2}(M)}}\int^{\infty}_{\delta_{c}}d\delta~{}e^{\frac{-\delta^{2}}{2\sigma^{2}(M)}}=\frac{\Gamma\left(\frac{1}{2},\frac{\delta_{c}^{2}}{2\sigma^{2}}\right)}{2\sqrt{\pi}}.$
(23)

The variance of the curvature perturbation $\sigma$ is related to the power
spectrum:

$\sigma^{2}(M_{PBH}(k))=\frac{4(1+\omega)^{2}}{(5+3\omega)^{2}}\int\frac{dk^{\prime}}{k^{\prime}}\left(\frac{k^{\prime}}{k}\right)^{4}P_{R}(k^{\prime})\tilde{W}^{2}\left(\frac{k^{\prime}}{k}\right)$
(24)

where $\tilde{W}$ is a window function. Throughout this work we consider the
Gaussian distribution for this function. With $\omega$ we denote the equation
of state and as we work in radiation dominated epoch we consider that
$\omega=1/3$.

In the PT, the number density of peaks above a threshold is given by the
expression:

$n_{peaks}=\frac{1}{(2\pi)^{2}}\left(\frac{\left\langle{{k^{2}}}\right\rangle}{3}\right)^{3/2}\left(\left(\frac{\delta_{c}}{\sigma}\right)^{2}-1\right)exp\left(-\frac{\delta_{c}^{2}}{2\sigma^{2}}\right)$
(25)

where $\left\langle{{k^{2}}}\right\rangle$ reads as:

$\left\langle{{k^{2}}}\right\rangle=\frac{1}{\sigma^{2}}\int^{\infty}_{0}\frac{dk^{\prime}}{k^{\prime}}{k^{\prime}}^{2}\tilde{W}\left(\frac{k^{\prime}}{k}\right)^{2}P_{\Delta}(k^{\prime}).$
(26)

The density power spectrum is defined as:

$P_{\Delta}=\frac{4(1+\omega)^{2}}{(5+3\omega)^{2}}\left(\frac{k^{\prime}}{k}\right)^{4}P_{R}(k^{\prime}).$
(27)

The mass fraction for a flat Universe assuming a Gaussian distribution is
given by the expression:

$\beta_{PT}=n_{peaks}(2\pi)^{3/2}(1/k)^{3}.$ (28)

Figure 3:  Fractional abundance of PBH for $W_{1}$ and $W_{2}$ (or
$W^{\prime}_{1}$ and $W^{\prime}_{2}$).

In figure 3 we present the abundances derived by the superpotentials $W_{1}$
(red) and $W_{2}$ (blue), or their equivalent forms $W_{1}^{\prime}$ and
$W_{2}^{\prime}$. Solid lines correspond to the PT and dashed line to PS. We
notice that there is a difference between these two methods. However this
difference can be compensated by the fact that the result are sensitive to the
value of $\delta_{c}$. In this figure we consider that $\delta_{c}=0.45$. The
choice of parameters are the same with the figure 2. In this figure we also
depict for comparison reason the observable bounds that come from
extragalactic gamma ray from PBH evaporation ($EG\gamma$), microlensing for
Subaru (HSC), Eros/Macho/Ogle and Kepler dynamical heating of ultra faint
dwarf (UFD), CMB measurements and radio observation.

In table LABEL:tabu1 we show the fractional abundances for the two proposed
models $W_{1}$ and $W_{2}$. We have notice that the values for fractional
abundance using PT are bigger than those of PS, as we previously observe by
the figure 3. In this table we show the value of the peak of the power
spectrum and the value for the observable constraints of $n_{s}$ and $r$. It
is important to remark here that our proposed models are in complete
consistence with the observable constraints released by Planck collaboration
2018[4]. The $n_{s}$ and $r$ are calculated at leading order in the slow-roll
expansion from the following expressions:

$n_{s}\simeq 1+2\eta_{V}-6\varepsilon_{V},\quad r\simeq 16\varepsilon_{V},$
(29)

where the slow-roll $\varepsilon_{V}$ and $\eta_{V}$ parameters are given by
the following expression:

$\varepsilon_{V}=\frac{1}{2}\left(\frac{V^{\prime}(\chi)}{V(\chi)}\right)^{2},\quad\eta_{V}=\frac{V^{\prime\prime}(\chi)}{V(\chi)}.$
(30)

Table 1: The value for the $n_{s}$, $r$ | $n_{s}$ | $r$ | $P_{R}^{peak}$ | $f_{PBH}^{PS}$ | $f_{PBH}^{PT}$
---|---|---|---|---|---
1 | $0.9634$ | $0.0099$ | $0.036$ | $0.00239$ | $0.752$
2 | $0.9605$ | $0.0127$ | $0.034$ | $0.00201$ | $0.651$

##  5 Production of GWs

The energy density of the GWs in terms of scalar power spectrum is given by
thr integral:

${\Omega_{GW}(k)}=\frac{c_{g}\Omega_{r}}{36}\int^{\frac{1}{\sqrt{3}}}_{0}\mathrm{d}d\int^{\infty}_{\frac{1}{\sqrt{3}}}\mathrm{d}s\left[\frac{(s^{2}-1/3)(d^{2}-1/3)}{s^{2}+d^{2}}\right]^{2}\times
P_{R}(kx)P_{R}(ky)(I_{c}^{2}+I_{s}^{2})$ (31)

where the radiation density is $\Omega_{r}\approx 5.4\times 10^{-5}$. The
variables $x$ and $y$ are defined as:

$x=\frac{\sqrt{3}}{2}(s+d),\quad y=\frac{\sqrt{3}}{2}(s-d).$ (32)

Finally, the functions $I_{c}$ and $I_{s}$ are given by the equations:

$\displaystyle
I_{c}=-36\pi\frac{(s^{2}+d^{2}-2)^{2}}{(s^{2}-d^{2})^{3}}\theta(s-1)$ (33)
$\displaystyle
I_{s}=-36\frac{(s^{2}+d^{2}-2)^{2}}{(s^{2}-d^{2})^{2}}\left[\frac{(s^{2}+d^{2}-2)}{(s^{2}-d^{2})}log\left|\frac{d^{2}-1}{s^{2}-1}\right|+2\right].$
(34)

Figure 4:  Energy density of GWs for $W_{1}$ and $W_{2}$ (or $W^{\prime}_{1}$
and $W^{\prime}_{2}$).

In figure 4 we present the energy density for GWs for the superpotentials
$W_{1}$ (red) and $W_{2}$ (blue), or for their equivalent forms
$W_{1}^{\prime}$ and $W_{2}^{\prime}$. The choice of parameters are the same
with figure 2, as in the study of PBHs. We plot them in respect to the
frequency in order to have a comparison with future space-based experiments,
such as LISA, DECIGO and BBO. We can notice that both proposed models $W_{1}$
and $W_{2}$ give us significant results for the amount of GWs.

##  6 Fine-tuning

Finally, the fine-tuning is regarded as a drawback in models with an
inflection point. In our analysis the parameters that need fine tuning are
$b_{1,2}$. In figure 5 we depict how the value of $b_{2}$, for instance, which
comes from $W_{2}$, depends on the enhancement of the power spectrum. In the
case of the production of PBHs a peak at least $10^{-2}$ is required.
Specifically, the value of the parameters $b_{1,2}$ which are related of the
value of the peak need to take such values in order to use the acceptable
threshold $\delta c$. As we see in the right panel of figure 5, if we want to
use the value of $\delta c$ at a range between $[0.4-0.6]$, we should fix the
value of $b_{2}$ according to this plot. This plot is made assuming that the
fractional abundance of PBHs is $f_{PBH}=0.1$. We numerically evaluate the
value of $b_{2}$ by repeating many times the calculation for PBHs production.
Specifically, we define a value for $\delta c$ and we stop the running when
the $f_{PBH}$ reaches the value $0.1$. The dashed line corresponds to PS
approach and the solid to the PT. We can observe again a difference between
these two approaches.

As a next step in order to have a measurement for the required fine tuning we
evaluate the quantity of $\Delta$, which is the maximum value of the following
logarithmic derivative:

$\Delta=\mathrm{Max}\left|\frac{\partial ln(P_{R}^{PEAK})}{\partial
ln(b_{i})}\right|.$ (35)

If the max value of this logarithmic derivative is big, we need a lot of fine
tuning. In case of studying PBHs the value of $\Delta$ is up to $10^{6}$.
Hence a lot of fine tuning is required. This is due to the value of
$\delta_{c}$, which should be $\delta_{c}=0.4-0.6$.

Figure 5:  How $b_{2}$ releated to the enchantment of power spectrum and to
the $\delta_{c}$.

However in the study of GWs the required fine-tuning can be decreased because
there is not the restriction of $\delta c$ and because the peak of scalar
power spectrum should not be at least $10^{-2}$. Therefore, if we demand a
lower peak, this quantity $\Delta$ can be decreased by two or three order of
magnitude.

##  7 Conclusions

In this work we provide a class of superpotential in order to generate PBHs
and GWs, which is derived, by no-scale supergravity. We study superpotentials
using the non-compact SU(2,1)/SU(2)$\times$ U(1) symmetry. By a proper
modification in the superpotentials, which in unmodified case lead to
Starobinsky-like effective scalar potential, we derive potentials which are
capable not only to generate PBHs and GWs, but also to obtain the same
transformation laws in $(y_{1},y_{2})$ and keep the form of effective scalar
potential in both basis $(T,\varphi)$ and $(y_{1},y_{2})$. The machanism in
order to achieve enhancement in scalar power spectrum and explain the
production of PBHs and GWs is an inflection point in the effective scalar
potential.

As a next step, we evaluate the production of PBHs by using two approaches the
PS approach and the PT. Furthermore, we evaluate the abundances GWs by using
the scalar power spectra. Our proposed models give us significant results in
order to expain the production of PBHs and the DM in the Universe as well the
generation of GWs. For both the production of PBHs and GWs we assume that they
are formed in radiation dominated epoch.

Finally, we discuss the issue of fine-tuning of the parameters of the
effective scalar potential. We conclude that the apparent drawback of models
with an inflection point is that it is required a lot of fine-tuning.
Therefore, other theoretical approaches, such as multi-field inflation, should
be taken into consideration in the study of PBHs and GWs.

## Acknowledgments

The research work was supported by the Hellenic Foundation for Research and
Innovation (H.F.R.I.) under the First Call for H.F.R.I. Research Projects to
support Faculty members and Researchers and the procurement of high-cost
research equipment grant (Project Number: 824). I.D.S. would like to thank
V.C. Spanos for useful discussions.

## References

## References

  * [1] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116 (2016) no.6, 061102 doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * [2] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116 (2016) no.24, 241103 doi:10.1103/PhysRevLett.116.241103 [arXiv:1606.04855 [gr-qc]].
  * [3] B. P. Abbott et al. [LIGO Scientific and VIRGO], Phys. Rev. Lett. 118 (2017) no.22, 221101 [erratum: Phys. Rev. Lett. 121 (2018) no.12, 129901] doi:10.1103/PhysRevLett.118.221101 [arXiv:1706.01812 [gr-qc]].
  * [4] Y. Akrami et al. [Planck], Astron. Astrophys. 641 (2020), A10 doi:10.1051/0004-6361/201833887 [arXiv:1807.06211 [astro-ph.CO]].
  * [5] J. Ellis, D. V. Nanopoulos and K. A. Olive, Phys. Rev. Lett. 111 (2013), 111301 [erratum: Phys. Rev. Lett. 111 (2013) no.12, 129902] doi:10.1103/PhysRevLett.111.111301 [arXiv:1305.1247 [hep-th]].
  * [6] J. Ellis, D. V. Nanopoulos, K. A. Olive and S. Verner, JHEP 03 (2019), 099 doi:10.1007/JHEP03(2019)099 [arXiv:1812.02192 [hep-th]].
  * [7] J. Ellis, D. V. Nanopoulos and K. A. Olive, JCAP 10 (2013), 009 doi:10.1088/1475-7516/2013/10/009 [arXiv:1307.3537 [hep-th]].
  * [8] E. Cremmer, S. Ferrara, C. Kounnas and D. V. Nanopoulos, Phys. Lett. B 133 (1983), 61 doi:10.1016/0370-2693(83)90106-5
  * [9] J. R. Ellis, C. Kounnas and D. V. Nanopoulos, Nucl. Phys. B 241 (1984), 406-428 doi:10.1016/0550-3213(84)90054-3
  * [10] I. D. Stamou, Phys. Rev. D 103 (2021) no.8, 083512 doi:10.1103/PhysRevD.103.083512 [arXiv:2104.08654 [hep-ph]].
  * [11] C. Ringeval, Lect. Notes Phys. 738 (2008), 243-273 doi:10.1007/978-3-540-74353-8_7 [arXiv:astro-ph/0703486 [astro-ph]].
  * [12] D. V. Nanopoulos, V. C. Spanos and I. D. Stamou, Phys. Rev. D 102 (2020) no.8, 083536 doi:10.1103/PhysRevD.102.083536 [arXiv:2008.01457 [astro-ph.CO]].


$\langle$reception date$\rangle$ $\langle$acception date$\rangle$
$\langle$publication date$\rangle$

protoplanetary disks — hydrodynamics — instabilities — turbulence

# A self-consistent model for dust settling and the vertical shear instability
in protoplanetary disks

Yuya Fukuhara11affiliation: Department of Earth and Planetary Sciences, Tokyo
Institute of Technology, Meguro, Tokyo 152-8551, Japan Satoshi
Okuzumi11affiliationmark<EMAIL_ADDRESS>
###### Abstract

The spatial distribution of dust particles in protoplanetary disks affects
dust evolution and planetesimal formation processes. The vertical shear
instability (VSI) is one of the candidate hydrodynamic mechanisms that can
generate turbulence in the outer disk region and affect dust diffusion.
Turbulence driven by the VSI has a predominant vertical motion that can
prevent dust settling. On the other hand, the dust distribution controls the
spatial distribution of the gas cooling rate, thereby affecting the strength
of VSI-driven turbulence. Here, we present a semi-analytic model that
determines the vertical dust distribution and the strength of VSI-driven
turbulence in a self-consistent manner. The model uses an empirical formula
for the vertical diffusion coefficient in VSI-driven turbulence obtained from
our recent hydrodynamical simulations. The formula returns the vertical
diffusion coefficient as a function of the vertical profile of the cooling
rate, which is determined by the vertical dust distribution. We use this model
to search for an equilibrium vertical dust profile where settling balances
with turbulent diffusion for a given maximum grain size. We find that if the
grains are sufficiently small, there exists a stable equilibrium dust
distribution where VSI-driven turbulence is sustained at a level of
$\alpha_{z}\sim 10^{-3}$, where $\alpha_{z}$ is the dimensionless vertical
diffusion coefficient. However, as the maximum grain size increases, the
equilibrium solution vanishes because the VSI can no longer stop the settling
of the grains. This runaway settling may explain highly settled dust rings
found in the outer part of some protoplanetary disks.

##  1 Introduction

The initial stage of planet formation is the formation of kilometer-sized
planetesimals from micron-sized dust grains (for reviews, Johansen et al.
2014; Drazkowska et al. 2023). These dust growth and planetesimal formation
depend on the vertical distribution of dust particles. If the dust density at
the midplane well exceeds the gas density, the streaming and gravitational
instabilities set in (e.g., Goldreich & Ward 1973; Sekiya 1998; Youdin & Shu
2002; Youdin & Goodman 2005), leading to planetesimal formation. Understanding
what processes determine the dust vertical profile in disks is also essential
for interpreting recent high-resolution radio observations showing that the
dust at different radial locations or in different disks is settled to
different extents (Pinte et al. 2016; Doi & Kataoka 2021; Villenave et al.
2022, 2023; Pizzati et al. 2023, and for a review, Miotello et al. 2023).

In protoplanetary disks, the vertical profile of dust particles critically
depends on the intensity of gas disk turbulence causing dust diffusion. The
question then is what mechanisms generate disk turbulence. Recent theoretical
studies have shown that not only the magnetorotational instability (MRI;
Balbus & Hawley 1991) but also some thermo-hydrodynamical instabilities are
important turbulence-driving mechanisms (for reviews, Lyra & Umurhan 2019;
Lesur et al. 2023). Among them, the vertical shear instability (VSI; Urpin &
Brandenburg 1998; Arlt & Urpin 2004; Nelson et al. 2013; Lin & Youdin 2015) is
thought to be the leading mechanism for driving turbulence in outer disk
regions where the current radio observations best constrain the degree of dust
settling. The VSI requires rapid cooling of disk gas in addition to a
vertically varying gas orbital velocity (Urpin, 2003; Nelson et al., 2013; Lin
& Youdin, 2015; Manger et al., 2021). Therefore, the VSI tends to operate in
the outer disk region with low optical depths (Malygin et al., 2017; Pfeil &
Klahr, 2019; Fukuhara et al., 2021; Melon Fuksman et al., 2024a), where it may
dominate over the MRI (e.g., Cui & Bai 2022). The VSI generates turbulence
with a predominant vertical gas motion (e.g., Nelson et al. 2013; Stoll & Kley
2014) that can prevent dust settling (Stoll & Kley, 2016; Flock et al., 2017b,
2020; Dullemond et al., 2022). This turbulence can also stir dust particles
significantly (Stoll & Kley, 2016; Flock et al., 2017b) and thus induce their
collisional velocities, leading to suppression of planetesimal formation
through coagulation (e.g, Ormel & Cuzzi 2007; Brauer et al. 2008; Okuzumi &
Hirose 2012). On the other hand, VSI-driven turbulence can produce both small
short-lived and azimuthally large long-lived vortices (Richard et al., 2016;
Manger & Klahr, 2018; Flock et al., 2020; Pfeil & Klahr, 2021; Melon Fuksman
et al., 2024b). These vortices may lead to dust concentration and subsequent
planetesimal formation through gravitational collapse (e.g., Barge & Sommeria
1995; Raettig et al. 2021; Lehmann & Lin 2022).

Importantly, dust particles are the dominant opacity source and determine the
local cooling rates of protoplanetary disks (Malygin et al., 2017; Barranco et
al., 2018). Therefore, their spatial distribution controls where in the disks
the VSI operates (Pfeil & Klahr, 2019; Fukuhara et al., 2021) and can even
affect the dust vertical diffusivity at the midplane (Pfeil & Klahr, 2021;
Fukuhara et al., 2023; Pfeil et al., 2023). The consequence of this dust–VSI
thermal interaction for dust settling remains unclear.

In this paper, we model the above-mentioned interaction between dust and the
VSI to study how dust setting and diffusion balance in VSI-driven turbulence.
To this end, we present a semi-analytic model that determines the vertical
dust distribution and the strength of VSI-driven turbulence in a self-
consistent manner. The model uses an empirical formula for the vertical
diffusion coefficient in VSI-driven turbulence obtained from our recent
hydrodynamical simulations (Fukuhara et al., 2023). The formula returns the
vertical diffusion coefficient as a function of the vertical profile of the
cooling rate, which is determined by the vertical dust distribution. We use
this model to search for an equilibrium vertical dust profile where settling
balances with turbulent diffusion for a given grain size.

This paper is organized as follows. In section 2, we describe our self-
consistent model. We present the main results in section 3 and discuss the
implication of our study in section 4. Section 5 presents a summary.

##  2 Method

Figure 1: Overview showing the self-consistent model determining the vertical
dust profile and VSI-driven turbulence intensity in this study. The model
considers a protoplanetary disk consisting of gas and dust (section 2.1).
Assuming dust particle size and vertical diffusion coefficient $\alpha_{z,\rm
assume}$, we compute the vertical profile of dust particles for the single-
sized model and power-law size distribution model (section 2.2). The dust
vertical profile determines the vertical profile of the cooling timescale
(section 2.3) and thereby localizes the region where the VSI operates (VSI-
unstable layer), using the linear VSI criterion (section 2.4). Based on the
empirical formula derived from the hydrodynamical simulations, we estimate the
vertical diffusion coefficient of VSI-driven turbulence $\alpha_{z,\rm VSI}$
(section 2.5). Comparing estimated coefficient $\alpha_{z,\rm VSI}$ with
assumed coefficient $\alpha_{z,\rm assume}$, we search an equilibrium vertical
dust profile where settling balances with turbulent diffusion for a given
grain size (section 2.6). Our parameters are the grain size for the single-
sized model or the maximum grain size for the power-law size distribution
model, and the distance from the central star (section 2.7).

In this section, we describe our self-consistent model that determines the
dust vertical distribution and VSI-driven turbulence intensity simultaneously
(see figure 1 for an overview). We assume a protoplanetary disk consisting of
gas (section 2.1) and dust (section 2.2). The dust grains dominate the disk’s
opacity (section 2.2) and their distribution determines the disk’s local
cooling rate distribution (section 2.3). Because the VSI requires rapid gas
cooling, the cooling rate distribution in turn determines the location where
the linear VSI operates (section 2.4), which we call the VSI-unstable
layer111Fukuhara et al. (2021) called these the VSI zones.. The VSI produces
turbulent gas motion propagating across the boundary between the VSI-unstable
and stable regions. We calculate the VSI-driven vertical diffusion coefficient
from the thicknesses of the VSI-unstable and -stable layers using an empirical
formula based on our previous hydrodynamical simulations (section 2.5). The
vertical diffusion coefficients can be used to update the vertical dust
distribution. Iterating these calculating steps, we determine the vertical
dust distribution and the strength of VSI-driven turbulence in a self-
consistent manner (section 2.6). In the following subsections, we describe the
model in more detail.

###  2.1 Disk model

We consider an axisymmetric disk around a solar-mass star. We adopt the
cylindrical coordinate system $(R,~{}z)$, where $R$ and $z$ are the distance
from the central star and the height from the midplane, respectively. The gas
surface density is given by

$\Sigma_{\rm g}(R)=\frac{\left({2-\beta_{\Sigma}}\right)M_{{\rm disk}}}{2\pi
R_{c}^{2}}\left({\frac{R}{R_{c}}}\right)^{-\beta_{\Sigma}}\exp{\left[-\left({\frac{R}{R_{c}}}\right)^{2-\beta_{\Sigma}}\right]},$
(1)

where $M_{{\rm disk}}$ is the total mass of the gas disk, $R_{c}$ is the
characteristic radius, and $\beta_{\Sigma}$ is a dimensionless number
characterizing the radial slope of the gas surface density. We fix the gas
disk parameters to $M_{{\rm disk}}=0.01\MO$, $R_{c}=100~{}{\rm au}$, and
$\beta_{\Sigma}=1$.

Assuming that the disk is optically thick to stellar radiation and that the
stellar luminosity is equal to the solar luminosity, the temperature of the
disk is given by

$T(R)=T_{0}\left({\frac{R}{1{\rm~{}au}}}\right)^{q},$ (2)

with $T_{0}=130{\rm~{}K}$ and $q=-3/7$ (Chiang & Goldreich, 1997). We assume
that the disk is vertically isothermal, ignoring warmer surface layers that
are optically thin to the starlight (Chiang & Goldreich, 1997).

From vertical hydrostatic equilibrium, the gas density is given by

$\rho_{\rm g}(R,~{}z)=\frac{\Sigma_{\rm g}}{\sqrt{2\pi}H_{\rm
g}}\exp{\left(-\frac{z^{2}}{2H_{\rm g}^{2}}\right)}.$ (3)

where $H_{\rm g}=c_{\rm s}/\Omega_{\rm K}$ is the gas scale height with
$c_{\rm s}$ and $\Omega_{\rm K}$ being the isothermal sound speed and
Keplerian frequency, respectively. The isothermal sound speed is given by
$c_{\rm s}=\sqrt{k_{\rm B}T/m_{\rm g}}$, where $k_{\rm B}$ is the Boltzmann
constant and $m_{\rm g}$ is the mean molecular mass of the gas. The Keplerian
frequency is given by $\Omega_{\rm K}=\sqrt{GM_{*}/R^{3}}$, where $G$ is the
gravitational constant and $M_{*}$ is the mass of the central star. In this
study, $m_{\rm g}$ and $M_{*}$ are taken to be $2.3m_{\rm p}$ and $1\MO$,
respectively, where $m_{\rm p}$ is the proton mass.

###  2.2 Dust model

We here describe the dust model that we use to calculate the cooling timescale
(see section 2.3). The ratio between the dust surface density $\Sigma_{\rm d}$
and $\Sigma_{\rm g}$ is fixed to the interstellar dust abundance of $1\%$,
whereas the local dust-to-gas ratio is allowed to vary with $z$ considering
dust settling.

We consider two grain size-distribution models. The first model is the single-
sized model where all dust particles are assumed to have equal size $a_{1}$.
Formally, the size distribution for the single-sized model can be written as

$\frac{dN_{\rm d}(a)}{da}=\frac{3\Sigma_{\rm d}}{4\pi\rho_{\rm
int}a_{1}^{3}}\delta(a-a_{1})$ (4)

where $dN_{\rm d}(a)/da$ is the number surface density per unit particle size
$a$, $\Sigma_{\rm d}$ is the total dust mass surface density, $\rho_{\rm int}$
is the grains’ internal density, and $\delta$ is the delta function. Equation
(4) satisfies the normalization

$\Sigma_{\rm d}=\int m_{\rm d}\frac{dN_{\rm d}(a)}{da}da,$ (5)

where $m_{\rm d}=(4\pi/3)\rho_{\rm int}a^{3}$ is the particle mass. The second
model is the power-law model where the grain size distribution is given by

$\frac{dN_{\rm
d}(a)}{da}=\left\\{\begin{array}[]{ll}{\displaystyle\frac{(12+3p)\Sigma_{\rm
d}}{4\pi\rho_{\rm int}\left(a_{\rm max}^{4+p}-a_{\rm
min}^{4+p}\right)}a^{p},}&a_{\rm min}<a<a_{\rm max},\\\ 0,&{\rm
otherwise},\end{array}\right.$ (6)

where $p~{}(\neq-4)$ is the slope of size distribution, and $a_{{\rm min}}$
and $a_{{\rm max}}$ are the minimum and maximum particle sizes, respectively.
We fix $p=-3.5$ and $a_{{\rm min}}=1~{}\micron$. The maximum particle size
serves as a free parameter in this study (see also section 2.7). Equation (6)
also fulfills equation (5).

Dust particles settle toward the midplane owing to stellar gravity and diffuse
away from the midplane owing to turbulence. We assume that the turbulent
diffusion coefficient is constant in the vertical direction. This approach is
valid if the VSI is the dominant source of disk turbulence and determines dust
vertical diffusion. Even if there is the linear-VSI stable layer near the
midplane, VSI-driven turbulence can penetrate the stable midplane layer and
thereby produce nearly constant turbulent vertical intensity in the vertical
direction (Pfeil & Klahr, 2021; Fukuhara et al., 2023). This situation can be
realized when the stable midplane layer thickness is thinner than two gas
scale heights, or the thickness of the unstable layer above the stable
midplane layer is thicker than a few gas scale heights (Fukuhara et al.,
2023). In this study, the results for all parameter ranges satisfy these
conditions. We note that, however, the diffusion coefficient may not be
uniform in the vertical direction when these conditions are broken.

Assuming the balance between settling and diffusion, the vertical distribution
of the particles can be written as (Takeuchi & Lin, 2002)

$\frac{dn_{\rm d}(a,z)}{da}=C_{d}(a)\exp{\left[-\frac{z^{2}}{2H_{\rm
g}^{2}}-\frac{{\rm St}_{\rm
mid}(a)}{\alpha_{z}}\left({\exp\frac{z^{2}}{2H_{\rm
g}^{2}}-1}\right)\right]},$ (7)

where $dn_{\rm d}(a,z)/da$ is the particle number density per unit radius at
height $z$, $\mathrm{St}_{{\rm mid}}(a)$ is the Stokes number of the particles
at the midplane, $\alpha_{z}$ is a dimensionless parameter that characterizes
the level of the dust vertical diffusion caused by turbulence, and $C_{d}(a)$
is the normalized constant determined by the condition $dN_{\rm
d}(a)/da=\int(dn_{\rm d}(a,z)/da)dz$. Because most dust particles lie at the
region of $z\ll H_{\rm g}$, the exponential factor in equation (7) can be
approximated as $\exp[-z^{2}/(2H_{\rm d}^{2})]$, yielding (Fukuhara et al.,
2021)

$C_{d}(a)=\frac{1}{\sqrt{2\pi}H_{\rm d}(a,\alpha_{z})}\frac{dN_{\rm
d}(a)}{da}.$ (8)

Here, $H_{\rm d}(a,\alpha_{z})$ is the scale height of particles with size $a$
given by (Dubrulle et al., 1995; Youdin & Lithwick, 2007)

$H_{\rm d}(a,\alpha_{z})=\left[1+\frac{\mathrm{St}_{{\rm
mid}}(a)}{\alpha_{z}}\right]^{-1/2}H_{\rm g}.$ (9)

The Stokes number is the product of the stopping time and Keplerian frequency.
Assuming that the particle radius is smaller than the mean free path of the
disk gas molecules, gas drag onto the particles follows Epstein’s law, which
gives (see, e.g., Birnstiel et al., 2010)

$\mathrm{St}_{{\rm mid}}(a)=\frac{\pi\rho_{{\rm int}}a}{2\Sigma_{\rm g}}.$
(10)

The vertical–size distribution $dn_{\rm d}(a,z)/da$ gives the collisional heat
transfer and opacity, which control the cooling timescale, as a function of
$z$. The mean travel length of gas molecules colliding with dust particles
$\ell_{{\rm gd}}$ is given by

$\ell_{{\rm gd}}=\left(\int\pi a^{2}\frac{dn_{\rm d}}{da}da\right)^{-1}.$ (11)

This determines the timescale of collisional heat transfer (see section 2.3).

To calculate the cooling time, we use the Plank mean opacity and Rosseland
mean opacity defined by

$\kappa_{\rm P}=\frac{\displaystyle\int_{0}^{\infty}\kappa_{{\rm
g},\lambda}B_{\lambda}(T)d\lambda}{\displaystyle\int_{0}^{\infty}B_{\lambda}(T)d\lambda},$
(12)

and

$\frac{1}{\kappa_{\rm
R}}=\frac{\displaystyle\int_{0}^{\infty}\frac{1}{\kappa_{{\rm
g},\lambda}}\frac{\partial B_{\lambda}(T)}{\partial
T}d\lambda}{\displaystyle\int_{0}^{\infty}\frac{\partial
B_{\lambda}(T)}{\partial T}d\lambda},$ (13)

respectively, where $\kappa_{{\rm g},\lambda}$ is the wavelength-dependent
opacity per gas mass, $\lambda$ is the wavelength, and $B_{\lambda}(T)$ is the
Planck function.

The opacity per gas mass is related to the size distribution of the dust
particles as (e.g., Kondo et al. 2023)

$\kappa_{{\rm g},\lambda}=\frac{1}{\rho_{\rm g}}\int\kappa_{{\rm
d},\lambda}(a)m_{\rm d}\frac{dn_{\rm d}}{da}da,$ (14)

where $\kappa_{{\rm d},\lambda}(a)$ is the opacity per dust mass. Assuming
that a dust particle is an uniform sphere, $\kappa_{{\rm d},\lambda}(a)$ can
be written as (Kataoka et al., 2014)

$\kappa_{{\rm d},\lambda}(a)=\frac{\pi a^{2}}{m_{\rm
d}}\times\left\\{\begin{array}[]{ll}{\displaystyle\frac{24nkx}{\left(n^{2}+2\right)^{2}},}&x\leq
1,\\\
{\displaystyle\min\left\\{\frac{8kx}{3n}\left(n^{3}-\left(n^{2}-1\right)^{3/2}\right),0.9\right\\},}&x>1,\end{array}\right.$
(15)

where $x=2\pi a/\lambda$ is the size parameter and $n$ and $k$ are the real
and imaginary parts of the complex refractive index, respectively. We
calculate $n$ and $k$ at each wavelength using the Bruggeman mixing rule,
assuming that the dust is a mixture of silicate and ice with a mass mixing
ratio of 1:1 with no porosity. The values of $n$ and $k$ for silicate and ice
are taken from Draine (2003) and Warren & Brandt (2008), respectively. The
internal density of the dust grains is then $\rho_{\rm int}=1.46~{}{\rm
g~{}cm^{-3}}$.

###  2.3 Cooling model

The spatial profile of the gas disk cooling rate depends on the size and
spatial distribution of the dust particles (Malygin et al., 2017; Barranco et
al., 2018; Fukuhara et al., 2021). Using the dust model described in section
2.2, we approximate the local cooling (thermal relaxation) timescale
$\tau_{\rm relax}$ as (Malygin et al., 2017; Pfeil & Klahr, 2019)

$\tau_{\rm relax}(z)=\max{\\{\tau_{\rm diff}(z),~{}\tau_{\rm
coll}(z),~{}\tau_{\rm emit}(z)\\}},$ (16)

where $\tau_{\rm diff}$, $\tau_{\rm coll}$, and $\tau_{\rm emit}$ are the
timescales of radiative diffusion, collisional heat transfer, and radiative
cooling, respectively.

The cooling time in the optically thick regime is dominated by the radiative
diffusion timescale $\tau_{\rm diff}$ as (Malygin et al., 2017)

$\tau_{\rm diff}=\frac{1}{\bar{D}k^{2}},$ (17)

where $\bar{D}$ and $k$ are the effective energy diffusion coefficient and the
wavenumber of the perturbation, respectively. The effective energy diffusion
coefficient is given by

$\bar{D}=\frac{\lambda_{\rm r}c}{\kappa_{\rm R}(T)\rho_{\rm
g}}\frac{4\eta}{1+3\eta}$ (18)

where $\lambda_{\rm r}$ is the flux limiter, $c$ is the light speed, and
$\eta$ is the ratio between radiation energy density $E_{\rm r}$ to combined
radiation and internal energy density $E_{\rm int}$. Here, $\eta$ is given by
$\eta=E_{\rm r}/(E_{\rm r}+E_{\rm int})$ with $E_{\rm r}=4\sigma_{\rm
SB}T^{4}/c$ and $E_{\rm int}=\rho_{\rm g}C_{V}T$, where $\sigma_{{\rm SB}}$ is
the Stefan–Boltzmann constant and $C_{V}=5k_{\rm B}/2m_{\rm g}$ is the
specific heat at constant volume. The flux limiter is fixed to $\lambda_{\rm
r}=1/3$, which is the value for the optically thick limit (Levermore &
Pomraning, 1981). For $k$, we assume that the thermal perturbation is
determined by the wavenumber of the VSI-driven turbulence structure. The VSI
unstable modes emerge when a radial perturbation wavenumber is larger than a
vertical one (e.g., Arlt & Urpin 2004), and typically takes $\sim 20/H_{\rm
g}$ by previous simulations of VSI-driven turbulence222 Note that the radial
wavenumber of the VSI-driven turbulence structure can depend on the density
and opacity. The previous hydrodynamical simulations including radiative
transport by Stoll & Kley (2014) have shown that an increase in the density
leads to a smaller radial wavelength of turbulence structure. (see appendix of
Pfeil & Klahr 2021). Therefore, we set $k=20/H_{\rm g}$.

The timescale of collisional heat transfer is given by (Fukuhara et al., 2021)

$\tau_{{\rm coll}}=\frac{\ell_{{\rm gd}}}{v_{{\rm th}}},$ (19)

where $v_{{\rm th}}$ is the mean relative velocity between the gas molecules
and dust particles. The relative velocity $v_{{\rm th}}$ can be approximated
as the mean thermal speed of the molecules,

$v_{{\rm th}}=\sqrt{\frac{8k_{\rm B}T}{\pi m_{\rm g}}}.$ (20)

We note that this collisional timescale is an approximation to the actual
thermal accommodation timescale of the gas molecules and dust particles, which
differs by a factor of $\gamma/(\gamma-1)$ (Burke & Hollenbach, 1983; Barranco
et al., 2018; Pfeil et al., 2023), where $\gamma$ is the heat capacity ratio.

The radiative cooling timescale $\tau_{{\rm emit}}$ in the optically thin
limit is given by (Malygin et al., 2017)

$\tau_{{\rm emit}}=\frac{C_{V}}{16\kappa_{\rm P}(T)\sigma_{{\rm SB}}T^{3}}.$
(21)

Both $\ell_{{\rm gd}}$ and $\kappa_{\rm P}$ depend on the local size
distribution of the dust particles (see section 2.2).

###  2.4 Defining the VSI-unstable layer

The cooling timescale determines the VSI-unstable layer because rapid cooling
reduces buoyancy preventing instability growth (e.g., Nelson et al. 2013).
Following Lin & Youdin (2015), the criterion for the VSI can be expressed as

$\tau_{\rm relax}(z)\lesssim\tau_{\rm crit}.$ (22)

Here, $\tau_{\rm crit}$ is the vertically global critical cooling timescale
defined by

$\tau_{\rm crit}=\frac{H_{\rm g}}{R}\frac{|q|}{\gamma-1}\Omega_{\rm K}^{-1},$
(23)

with $\gamma=1.4$. By applying the VSI criterion of equation (22) to each
point in the disk, we search for the VSI-unstable and -stable layers. We apply
this vertically global criterion instead of the local criterion [equation (4)
in Lin & Youdin 2015]. This is because the intensity of VSI-driven turbulence
is tightly correlated with the thicknesses of the VSI-unstable and -stable
layers predicted by the global criterion (Fukuhara et al. 2023, see also
section 2.5). The unstable and stable layers exist on the upper ($z>0$) and
lower ($z<0$) halves of the disk, symmetrically concerning the midplane
($z=0$, see the middle panel of figure 1). For the upper halves, we determine
the height of the unstable layer’s upper boundary $z_{\rm u}$ and the height
of the midplane stable layer’s upper boundary $z_{\rm s}$. When the midplane
stable layer is absent, we set $z_{\rm s}=0$. Following Fukuhara et al.
(2023), we also define the thicknesses of the linearly stable and unstable
layers as $\Delta L_{\rm s}=2z_{\rm s}$ and $\Delta L_{\rm u}=2z_{\rm
u}-\Delta L_{\rm s}$, respectively (see figure 1).

###  2.5 Estimating dust vertical diffusion coefficient

The thicknesses of the VSI-unstable and -stable layers can be used to
determine the vertical diffusion of dust particles (Fukuhara et al., 2023).
VSI-driven turbulence generally exhibits the vertical gas motion that is
uniform in the vertical direction (e.g., Pfeil & Klahr 2021; Fukuhara et al.
2023). Therefore, we estimate the dimensionless dust vertical diffusion
coefficient of VSI-driven turbulence $\alpha_{z,\rm VSI}$ as

$\alpha_{z,\rm VSI}=\frac{\langle v_{z}^{2}\rangle}{c_{\rm
s}^{2}}\cdot\tau_{\rm corr}\Omega_{\rm K},$ (24)

where $\langle v_{z}^{2}\rangle$ is the time-averaged squared gas vertical
velocity at the midplane and $\tau_{\rm corr}$ is the correlation time of VSI-
driven turbulence. The correlation time of VSI-driven turbulence is still
indeterminate, ranging from $\tau_{\rm corr}\Omega_{\rm K}\sim 0.2$ to $\sim
20$ (Stoll & Kley 2016; Flock et al. 2020, see also section 4.1 of Fukuhara et
al. 2023). In this study, we fix $\tau_{\rm corr}\Omega_{\rm K}=1.0$.

Figure 2: Dimensionless dust vertical diffusion coefficient of VSI-driven
turbulence $\alpha_{z,\rm VSI}$ [equation (24)] as a function of $\Delta
L_{\rm u}$ and $\Delta L_{\rm s}$. The dashed lines show $\alpha_{z,\rm
VSI}=10^{-3}$, $10^{-4}$, and $10^{-5}$.

Global axisymmetric simulations by Fukuhara et al. (2023) show that the mean
squared gas vertical velocity $\langle v_{z}^{2}\rangle$ in VSI-driven
turbulence is tightly correlated with the thicknesses of the VSI-unstable and
-stable layers. We estimate $\langle v_{z}^{2}\rangle$ from $\Delta L_{\rm u}$
and $\Delta L_{\rm s}$, using an empirical formula provided by Fukuhara et al.
(2023),

$\frac{\langle v_{z}^{2}\rangle}{c_{\rm s}^{2}}=f_{\rm T}(\Delta L_{\rm
u},~{}\Delta L_{\rm s})+f_{\rm pT}(\Delta L_{\rm u},~{}\Delta L_{\rm s}),$
(25)

where $f_{\rm T}$ and $f_{\rm pT}$ represent the dependence of $\Delta L_{\rm
u}$ and $\Delta L_{\rm s}$ on $\langle v_{z}^{2}\rangle$. For $f_{\rm T}$ and
$f_{\rm pT}$, we use equations (17) and (19) of Fukuhara et al. (2023), which
represent the sharp decrease of $\langle v_{z}^{2}\rangle$ in $\Delta L_{\rm
u}\lesssim 2H_{\rm g}$ and $\Delta L_{\rm s}\gtrsim 2H_{\rm g}$. In figure 2,
We plot $\alpha_{z,\rm VSI}$ [equation (24)] as a function of $\Delta L_{\rm
u}$ and $\Delta L_{\rm s}$. When $\Delta L_{\rm u}\gtrsim 2H_{\rm g}$ and
$\Delta L_{\rm s}\lesssim 2H_{\rm g}$, equation (24) predicts $\alpha_{z,\rm
VSI}\approx 2\times 10^{-3}$; otherwise, $\alpha_{z,\rm VSI}$ approaches zero.

###  2.6 Calculation procedure

The vertical dust distribution depends on the turbulent diffusion coefficient
$\alpha_{z}$ [equation (7)]. On the other hand, if we assume that the VSI is
the main driver of disk turbulence, its diffusion strength $\alpha_{z,\rm
VSI}$ depends on the disk’s cooling structure [equation (24)] and in turn on
vertical dust distribution. In general, $\alpha_{z,\rm VSI}$ produced by VSI-
driven turbulence under a given vertical dust distribution does not
necessarily match $\alpha_{z}$ required to maintain the dust distribution. We
search for an equilibrium state where the two turbulence strengths match in
the following steps.

  1. 1.

For a given grain size distribution and radial position $R$, use equation (7)
to generate vertical dust profiles for various values of trial turbulence
strengths $\alpha_{z}=\alpha_{z,\rm assume}$. The trial values are generated
by dividing the range $10^{-8}\leq\alpha_{z,\rm assume}\leq 10^{-2}$ into
grids of an equal logarithmic interval of 0.05.

  2. 2.

For each value of $\alpha_{z,\rm assume}$, use the cooling model [equation
(16)], VSI criterion [equation (22)], and the empirical formula for VSI-driven
turbulence strength [equation (24)] to evaluate $\alpha_{z,\rm VSI}$.

  3. 3.

Interpolate $\alpha_{z,\rm VSI}$ as a smooth function of $\alpha_{z,\rm
assume}$, and search for self-consistent equilibrium solutions where the
predicted VSI-driven turbulence strength $\alpha_{z,\rm VSI}$ equals the
assume turbulence strength $\alpha_{z,\rm assume}$. Below, we denote the
vertical diffusion coefficient for an equilibrium solution by $\alpha_{z,\rm
equi}$. If $\alpha_{z,\rm VSI}<\alpha_{z,\rm assume}$ for all values of
$\alpha_{z,\rm assume}$, VSI-driven turbulence cannot be strong enough to
sustain the assumed dust vertical distribution, meaning that it cannot stop
grain settling. We call this runaway dust settling.

For every set of grain size and $R$, we repeat this procedure and search for
equilibrium solutions.

We ignore that VSI-driven turbulence revives before the dust settles toward
the midplane in the runaway fashion. This assumption is valid if the dust-
settling timescale is longer than the growth timescale of turbulence. The
settling timescale can be estimated as $\sim z/|v_{z}|\approx{\rm St}_{\rm
mid}^{-1}\Omega_{\rm K}^{-1}$ (e.g., Dubrulle et al. 1995; Youdin & Lithwick
2007), where $v_{z}\approx-{\rm St}_{\rm mid}\Omega_{\rm K}z$ is the vertical
velocity of dust particles for the terminal velocity approximation. The
typical growth timescale of VSI-driven turbulence is $\sim 10^{2}\Omega_{\rm
K}^{-1}$ (e.g., Nelson et al. 2013), which is longer than the settling
timescale for large dust grains with ${\rm St}_{\rm mid}\gtrsim 10^{-2}$.
However, for small grains with ${\rm St}_{\rm mid}\lesssim 10^{-2}$, the
growth timescale is shorter than the settling timescale, suggesting re-
development of VSI-driven turbulence. In this case, weak turbulence may
operate before the dust settles.

###  2.7 Parameter choices

Our model involves two main parameters: the size of dust grains $a_{1}$ for
the single-sized model or the maximum size of dust grains $a_{\rm max}$ for
the power-law size distribution model, and the distance from the central star
$R$. We divide the parameter ranges $1~{}\micron<a<1{\rm~{}cm}$ or
$1~{}\micron<a_{\rm max}<1{\rm~{}cm}$ and $1{\rm~{}au}<R<100{\rm~{}au}$ into
grids of an equal logarithmic interval of 0.1.

##  3 Results

In this section, we use the model described in section 2 to search for
equilibrium solutions where dust settling balances with VSI-driven turbulent
diffusion. In section 3.1, we describe the properties of the equilibrium
solutions and their dependence on the dust grain size and radial distance for
the single-sized model. In section 3.2, we show the equilibrium solutions for
the power-law size distribution model.

###  3.1 Equilibrium solutions for the single-sized model

Figure 3: Equilibrium vertical diffusion coefficient $\alpha_{z.\rm equi}$
(upper panel) and its dust scale height (lower panel), for different values of
dust grain’s size at $R=10~{}{\rm au}$ in the single-sized model. The symbols
indicate stable (circles) and unstable (crosses) equilibrium solutions. The
dashed lines show a constant value with $\alpha_{z}=2\times 10^{-3}$. The
dotted lines in the lower panel represent dust scale height with
$\alpha_{z}=10^{-2},~{}10^{-3},~{}10^{-4},~{}10^{-5},~{}10^{-6},$ and
$10^{-7}$.

First, we use the single-sized model at $R=10{\rm~{}au}$ to illustrate the
properties of the equilibrium solutions. Following the procedure described in
section 2.6, we search for the equilibrium vertical diffusion coefficient
$\alpha_{z,\rm equi}$ for different values of single grain size $a_{1}$. The
upper panel of figure 3 shows $\alpha_{z,\rm equi}$ for each $a_{1}$. We find
that when $a_{1}\lesssim 10{\rm~{}\micron}$, there exist two equilibrium
solutions. The equilibrium solutions vanish when the grain size exceeds $\sim
10{\rm~{}\micron}$. We also plot in the lower panel of figure 3 the dust scale
height with equilibrium solutions defined by $H_{\rm d,equi}(a)\equiv H_{\rm
d}(a,\alpha_{z,\rm equi})$ as a function of $a_{1}$.

Figure 4 illustrates two cases at $R=10~{}\rm au$ with
$a_{1}=10~{}{\rm\micron}$ and $100~{}{\rm\micron}$, representing the cases
with two and no equilibrium solutions, respectively. The upper panel plots
$\alpha_{z,\rm VSI}$ as a function of $\alpha_{z,\rm assume}$. For
$a_{1}=10~{}{\rm\micron}$, $\alpha_{z,\rm VSI}$ drops sharply from $\approx
2\times 10^{-3}$ to $\ll 10^{-3}$ as $\alpha_{z,\rm assume}$ falls below $\sim
10^{-4}$. A decrease in $\alpha_{z,\rm assume}$ leads to dust depletion at
high altitudes, yielding a decrease in the thickness of the VSI-unstable
layer. This decrease in the unstable layer’s thickness results in a strong
suppression of VSI-driven turbulence around $\Delta L_{\rm u}\approx 2H_{\rm
g}$ (see figure 2). For $a_{1}=100~{}{\rm\micron}$, $\alpha_{z,\rm VSI}$ takes
lower values of $\lesssim 10^{-10}$. This is because the VSI-unstable layer is
thin ($\Delta L_{\rm u}\lesssim 2H_{\rm g}$) for any values of $\alpha_{z,\rm
assume}$ due to dust vertical settling, resulting in suppression of VSI-driven
turbulence.

Equilibrium solutions correspond to the points where $\alpha_{z,\rm
VSI}=\alpha_{z,\rm assume}$. For $a_{1}=10{\rm~{}\micron}$, two equilibrium
solutions exist at $\alpha_{z,\rm equi}\approx 2\times 10^{-3}$ and $2\times
10^{-5}$. For $a_{1}=100{\rm~{}\micron}$, no equilibrium solution exists
because $\alpha_{z,\rm VSI}<\alpha_{z,\rm assume}$ for all values of
$\alpha_{z,\rm assume}$. This implies that the vertical diffusion by VSI-
driven turbulence alone cannot sustain dust vertical distribution of any scale
height in the latter case.

To understand the results presented in the upper panel of figure 3 in terms of
the competition between vertical settling and turbulent diffusion, we consider
a time-dependent toy model where the vertical distribution of dust grains
evolves through settling and diffusion. The equation that describes the
evolution of the squared mean of the grains’ vertical positions, $\langle
z^{2}\rangle$, is given by (for a derivation, see appendix A)

$\frac{1}{2}\frac{d\langle z^{2}\rangle}{dt}=-{\rm St}_{\rm mid}\Omega_{\rm
K}\langle z^{2}\rangle+D_{z}(\alpha_{z,\rm VSI})\left(1-\frac{\langle
z^{2}\rangle}{H_{\rm g}^{2}}\right),$ (26)

where $D_{z}(\alpha_{z})=\alpha_{z}c_{\rm s}H_{\rm g}$ is the dust vertical
diffusion coefficient. Below we simply call $\langle z^{2}\rangle^{1/2}$ the
dust scale height. In the right-hand side of equation (26), the first term
represents settling toward the midplane at the terminal vertical velocity
$-{\rm St}_{\rm mid}\Omega_{\rm K}z$. The second term corresponds to the
vertical diffusion by VSI-driven turbulence in a disk of gas scale height
$H_{\rm g}$, with the factor $(1-{\langle z^{2}\rangle}/{H_{\rm g}^{2}})$
guaranteeing that the dust scale height never exceeds $H_{\rm g}$ in the limit
of strong diffusion (Ciesla, 2010). It can be shown from equation (9) that our
single-sized model assumes

$\langle z^{2}\rangle=\left[H_{\rm d}(a,\alpha_{z,\rm
assume})\right]^{2}=\left[1+\frac{{\rm St}_{\rm mid}(a)}{\alpha_{z,\rm
assume}}\right]^{-1}H_{\rm g}^{2}$ (27)

for an arbitrary value of $\alpha_{z,\rm assume}$. Substituting this into
equation (26), we find that $d\langle z^{2}\rangle/dt=0$ if $\alpha_{z,\rm
assume}=\alpha_{z,\rm VSI}$, meaning that the dust settling and diffusion
indeed balance in the equilibrium solutions defined in section 2.6. This toy
model also allows us to predict what would occur when $\alpha_{z,\rm assume}$
and $\alpha_{z,\rm VSI}$ are unequal. For example, if $\alpha_{z,\rm
assume}>\alpha_{z,\rm VSI}$ (i.e., if VSI-driven turbulence is not strong
enough to sustain the vertical dust distribution), equation (26) shows that
settling dominates over diffusion and the dust scale height decreases.

Figure 4: Upper panel: VSI-driven vertical diffusion coefficient
$\alpha_{z,\rm VSI}$ at $R=10~{}\rm au$ as a function of $\alpha_{z,\rm
assume}$ from the single-sized model with $a_{1}=10~{}\micron$ and
$100~{}\micron$. The dotted line represents $\alpha_{z,\rm VSI}=\alpha_{z,\rm
assume}$. The circle and cross symbols mark stable and unstable equilibrium
solutions, respectively. Lower panel: $d\langle z^{2}\rangle/dt$ as a function
of $\langle z^{2}\rangle$ for the two cases presented in the upper panel. The
dashed lines show cases of a constant vertical diffusion coefficient as
$\alpha_{z}=2\times 10^{-3}$.

We now use equation (26) to explain why either two or no equilibrium solutions
emerge depending on the value of $a_{1}$. The lower panel of figure 4 shows
$d\langle z^{2}\rangle/dt$ as a function of $\langle z^{2}\rangle$ for
$a_{1}=10$ and $100{\rm~{}\micron}$. In general, there exists one equilibrium
solution for a constant value of $D_{z}$ because $d\langle z^{2}\rangle/dt$ is
then a monotonically decreasing function of $\langle z^{2}\rangle$, i.e.,
diffusion (settling) tends to dominate at small (large) dust scale heights.
This is illustrated by the dashed curves, which show $d\langle
z^{2}\rangle/dt$ in the case where turbulence with $\alpha_{z,\rm VSI}\approx
2\times 10^{-3}$ is present independently of $\langle z^{2}\rangle$.
Therefore, multiple or no equilibrium solution occurs only when $D_{z}$ is a
function of $\langle z^{2}\rangle$. For $a_{1}=10{\rm~{}\micron}$, the two
equilibrium solutions have dust scale heights of $\langle
z^{2}\rangle^{1/2}\sim H_{\rm g}$ and $0.3H_{\rm g}$. Comparison with the
dashed curve shows that the solution with higher $\langle z^{2}\rangle^{1/2}$
corresponds to the equilibrium solution for constant diffusion coefficient
with $\alpha_{z,\rm VSI}\approx 2\times 10^{-3}$. The solution with lower
$\langle z^{2}\rangle^{1/2}$ emerges because VSI-driven turbulent diffusion is
suppressed for small dust scale heights, i.e., for small $\alpha_{z,\rm
assume}$, as mentioned earlier. For $a_{1}=100~{}{\rm\micron}$, $d\langle
z^{2}\rangle/dt$ is negative for all dust scale heights, indicating runaway
dust settling.

Equation (26) can also be used to predict whether the equilibrium solutions
are stable or not against a small perturbation in the dust scale height. In
general, equilibrium solutions can be classified into stable and unstable
ones. In this context, stable (unstable) solutions are the ones for which
$d\langle z^{2}\rangle/dt$ becomes negative (positive) as we slightly increase
$\langle z^{2}\rangle$ from the equilibrium value, so that the perturbation
decays (diverges). In other words, stable (unstable) solutions have a negative
(positive) slope of $d\langle z^{2}\rangle/dt$ as a function of $\langle
z^{2}\rangle$. In the case of $a_{1}=10~{}\rm\mu m$, the equilibrium solution
with the larger $\langle z^{2}\rangle$ is stable, whereas that with the
smaller $\langle z^{2}\rangle$ is unstable (see figure 4).

Figure 5: Vertical diffusion coefficients $\alpha_{z,\rm equi}$ for the stable
equilibrium solution from the single-sized model as a function of radial
distance $R$ and grain size $a_{1}$. The black area indicates the parameter
space where no equilibrium solution exists. In this area, $\alpha_{z,\rm VSI}$
is lower than $\alpha_{z,\rm assume}$ for all $\alpha_{z,\rm assume}$,
implying runaway dust settling. The dotted lines mark, from top to bottom,
${\rm St}_{\rm mid}(a_{1})=10^{-1},~{}10^{-2},~{}10^{-3},$ and $10^{-4}$.
Figure 6: Radial profiles of the height of the unstable layer’s upper boundary
(circles), $z_{\rm u}$, and the height of the midplane stable layer’s upper
boundary (crosses), $z_{\rm s}$, for $a_{1}=10{\rm~{}\micron}$ with
$\alpha_{z,\rm assume}=2\times 10^{-3}$ from the single-sized model. The
regions of $z_{\rm s}<z<z_{\rm u}$ and $0<z<z_{\rm s}$ indicate the unstable
and midplane stable layers, respectively. The left and right vertical dashed
lines mark $\Delta L_{\rm s}(\equiv 2z_{\rm s})\approx 2H_{\rm g}$ and $\Delta
L_{\rm u}(\equiv 2z_{\rm u}-\Delta L_{\rm s})\approx 2H_{\rm g}$,
respectively.

Next, we describe how the stable equilibrium solution varies with radial
distance $R$. Figure 5 shows $\alpha_{z,\rm equi}$ for the stable solution as
a function of $R$ and $a_{1}$ for the single-sized model. The black area in
the $R$–$a_{1}$ plane indicates the parameter space where no equilibrium
solution exists and runaway dust settling is expected. The figure shows that
equilibrium solutions only exist beyond $R\sim 3~{}\rm au$. At $R\lesssim
3~{}\rm au$, the midplane region is significantly optically thick to its own
thermal emission, yielding a thick VSI-stable layer ($\Delta L_{\rm s}\gtrsim
2H_{\rm g}$) around the midplane. We map in figure 6 the unstable and stable
layers for $a_{1}=10{\rm~{}\micron}$ with $\alpha_{z,\rm assume}=2\times
10^{-3}$ from the single-sized model at the $R$–$z$ plane. This figure
indicates that the thickness of the stable midplane layer increases sharply as
$R$ decreases from $\sim 3{\rm~{}au}$. This thick stable layer prevents VSI-
driven turbulence from developing ($\alpha_{z,\rm VSI}\ll 10^{-3}$; see figure
2), resulting in runaway dust settling ($\alpha_{z,\rm VSI}<\alpha_{z,\rm
assume}$ for all values of $\alpha_{z,\rm assume}$). For the cases of large
dust, optical depth decreases and the VSI-stable layer becomes thin. However,
because the dust settles vertically, the VSI-unstable layer is also thin
($\Delta L_{\rm u}\lesssim 2H_{\rm g}$), leading to suppression of VSI-driven
turbulence. The stable equilibrium solutions are generally accompanied by
fully developed VSI-driven turbulence with $\alpha_{z,\rm VSI}\approx 2\times
10^{-3}$. This is reasonable because any vulnerable turbulence whose strength
decreases steeply with decreasing dust scale height would lead to an unstable
solution.

At $R\gtrsim 3~{}\rm au$, equilibrium solutions exist for sufficiently small
$a_{1}$ as we already showed in the particular case of $R=10~{}\rm au$. The
maximum grain size for equilibrium decreases with increasing $R$. This is
because the dust number density decreases as $R$ increases: decreasing dust
density leads to a longer cooling time, which makes the VSI-unstable layers
vertically thinner (Malygin et al. 2017; Pfeil & Klahr 2019; Fukuhara et al.
2021).

Figure 7: Radial profiles of the dust scale height from the stable equilibrium
solution for the single-sized dust model with different values of $a_{1}$. The
dotted lines show the dust scale height for the fixed vertical diffusion
coefficient of $\alpha_{z}=2\times 10^{-3}$.

Figure 7 shows the radial profiles of the dust scale height from the stable
equilibrium solutions, $H_{\rm d,equi}(a_{1})$, for different values of
$a_{1}$. We set $H_{\rm d,equi}(a_{1})$ to zero if runaway settling occurs. As
mentioned above, the stable solutions are generally accompanied by fully
developed VSI-driven turbulence with $\alpha_{z,\rm equi}\approx 2\times
10^{-3}$. For reference, the dotted lines in the figure show the dust scale
height for $\alpha_{z,\rm equi}=2\times 10^{-3}$. Because the stable solutions
require small $a_{1}$, they always lead to a thick dust disk with $H_{\rm
d,equi}(a_{1})\approx H_{\rm g}$.

Because the dust scale height is a function of ${\rm St}_{\rm mid}(a_{1})$, we
can expect that the value of ${\rm St}_{\rm mid}(a_{1})$ rather than $a_{1}$
more directly determines whether there exist equilibrium solutions or not. To
confirm this expectation, we overplot in figure 5 contours of constant ${\rm
St}_{\rm mid}(a_{1})$. We find that the equilibrium solution at $R\gtrsim
3~{}\rm au$ vanishes for ${\rm St}_{\rm mid}(a_{1})\gtrsim 3\times 10^{-4}$.

###  3.2 Equilibrium solutions for the power-law size distribution model

The grain size distribution affects the cooling timescale profile. Of the
three timescales determining the cooling timescale (see section 2.3), the
timescales of collisional heat transfer $\tau_{\rm coll}$ and radiative
diffusion $\tau_{\rm diff}$ determine the thicknesses of the unstable and
midplane stable layers, respectively. For the power-law size distribution with
a slope of $p=-3.5$, the largest grains dominate the dust mass budget, whereas
the smallest grains dominate the total geometric cross-section and thereby
control collisional heat transfer. Therefore, even if large grains exist,
small grains can still sustain the unstable layer.

Figure 8: Same as figure 5, but from the power-law size distribution model.
The dotted lines mark, from top to bottom, ${\rm St}_{\rm mid}(a_{\rm
max})=10^{-1},~{}10^{-2},~{}10^{-3},$ and $10^{-4}$.

Figure 8 illustrates how the small grains in the power-law size distribution
model extend the parameter space where strong VSI-driven turbulence is
sustained. This figure shows $\alpha_{z,\rm equi}$ for the stable equilibrium
solutions as a function of $R$ and $a_{\rm max}$ for the power-law size
distribution model. At $R=10{\rm~{}au}$, the equilibrium solutions exist up to
$a_{\rm max}=100~{}\rm\mu m$. This is in contrast to the single-sized model,
which only allows equilibrium solutions up to $10{\rm~{}\micron}$ at the same
radial position (see figure 5). In terms of the Stokes number, the power-law
size distribution model allows the equilibrium solutions up to ${\rm St}_{\rm
mid}(a_{\rm max})\approx 10^{-3}$–$10^{-2}$ depending on $R$. As in the
single-sized model, the stable solution is accompanied by fully developed VSI-
driven turbulence with $\alpha_{z,\rm equi}\approx 2\times 10^{-3}$. At
$R\lesssim 3~{}{\rm au}$ or at sufficiently large $a_{\rm max}$, the
equilibrium solution vanishes and runaway settling ($\alpha_{z,\rm
assume}>\alpha_{z,\rm VSI}$) occurs for the same reason as in the single-sized
model.

Figure 9: Same as figure 7, but from the power-law size distribution model.

Figure 9 plots the dust scale height for the stable equilibrium solution as a
function of $R$ from the power-law size distribution model with different
values of $a_{\rm max}$. As in figure 7, we set $H_{\rm d,equi}=0$ for runaway
settling. For $a_{\rm max}=10~{}{\rm\micron}$, VSI-driven turbulence maintains
a thick dust layer of $H_{\rm d,equi}\gtrsim 0.8H_{\rm g}$ at
$2{\rm~{}au}\lesssim R\lesssim 40{\rm~{}au}$. As $a_{\rm max}$ increases, the
region with the stable equilibrium solution shrinks, and the equilibrium dust
scale height within that region decreases because the settling velocity
increases while keeping $\alpha_{z,\rm equi}\approx 2\times 10^{-3}$. For
$a_{\rm max}>1~{}{\rm mm}$, the equilibrium solution vanishes at all $R$.

##  4 Discussion

###  4.1 Implications for disk observations

In section 3, we have shown that the dust grain size and radial distance
determine whether the VSI can sustain high dust diffusion. VSI-driven
turbulence can maintain a thick dust layer of $H_{\rm d}\approx H_{\rm g}$
within an annular region with a moderate optical depth and small grains.
However, when both the dust size and radial distance are large, VSI-driven
turbulence is suppressed, resulting in runaway dust settling leading to
$H_{\rm d}\approx 0$.

These differences in dust settling levels may provide an explanation for
varying degrees of dust settling as inferred from recent radio interferometric
observations of some protoplanetary disks. The disks around HL Tau and Oph
163131 exhibit well-defined dust gaps, limiting the dust scale height to
$\lesssim 0.1H_{\rm g}$ at $R\gtrsim 10{\rm~{}au}$ (Pinte et al., 2016) and at
$R\approx 100{\rm~{}au}$ (Villenave et al., 2022), respectively. These
estimated dust scale heights are in line with runaway dust settling in the
outer disk regions predicted in this study. In contrast, the disk around HD
163296 has two major dust rings at $70$ and $100{\rm~{}au}$ that show high and
low degrees of dust settling with $H_{\rm d}/H_{\rm g}\gtrsim 0.8$ and
$\lesssim 0.1$, respectively (Doi & Kataoka, 2021). This difference in dust
diffusion by radial distance is consistent with the trend in the radial
profile of dust scale heights predicted by our model (see figures 7 and 9).
Therefore, we hypothesize that VSI-driven turbulence dominates vertical dust
diffusion in these disks, and that this turbulence is suppressed in the outer
disk regions. Testing this hypothesis requires detailed modeling of gas, dust,
and cooling rate profiles in these disks.

###  4.2 Implications for disk and dust evolution

In section 3, we have shown that there exists a parameter space where VSI-
driven turbulence is too weak to stop dust settling. This occurs either when
the VSI-stable layer at the midplane is too thick or when the VSI-unstable
layer is too thin for the VSI to operate. The former condition is met at small
orbital radii ($R\lesssim$ 3 au), whereas the latter condition is met when the
grains grow beyond a certain size.

This runaway settling is beneficial for dust growth and planetesimal
formation. A high degree of dust settling generally promotes planetesimal
formation through the streaming and gravitational instabilities (e.g., Sekiya
1998; Youdin & Shu 2002; Johansen et al. 2009; Gole et al. 2020; Umurhan et
al. 2020; Chen & Lin 2020). Weak turbulence leading to high dust concentration
is also conducive to dust growth through coagulation without collisional
fragmentation and erosion (e.g., Brauer et al. 2008; Okuzumi & Hirose 2012).

Figures 7 and 9 imply that strong dust diffusion by the VSI tends to occur in
an annular region where the optical depth is moderate (see also figure 6). The
vertically extended dust in this annulus may block the radiation of the
central star and thereby cast a shadow beyond the annulus (e.g., Dullemond et
al., 2001; Dullemond & Dominik, 2004). The shadowing generally results in a
significant drop in temperature and therefore can affect the chemical
evolution in that region (Ohno & Ueda, 2021; Notsu et al., 2022).

The inner and outer edges of the VSI-turbulence zone also serve as potential
sites for planetesimal formation. A steep radial change in turbulence
viscosity at the edges may trigger the Rossby wave instability (Lovelace et
al., 1999; Li et al., 2000, 2001) and create a long-lived vortex, promoting
dust concentration and subsequent planetesimal formation through gravitational
collapse (e.g., Barge & Sommeria, 1995). The sharp drop in turbulent
diffusivity at the inner edge of the VSI-driven turbulence zone can also
trigger a runaway pile-up of dust grains (Hyodo et al., 2021, 2022) because
the increase of dust-to-gas ratio at the midplane reduces the radial drift
velocity of dust. Moreover, if the turbulent viscosity caused by the VSI
dominates gas disk accretion, the outer edge of VSI-driven turbulence may
create a local maximum in the radial profile of the gas pressure because the
turbulent viscosity decreases sharply toward the outside. This is similar to
the mechanism by which a pressure maximum is generated near the dead-zone
inner edge of MRI (e.g., Dzyurkevich et al. 2010; Flock et al. 2016, 2017a).
The pressure maximum can trap dust particles (Whipple, 1972; Adachi et al.,
1976; Weidenschilling, 1977). All these dust concentration processes could
lead to planetesimal formation via the streaming and the gravitational
instabilities (e.g., Youdin & Goodman 2005; Johansen & Youdin 2007; Johansen
et al. 2009; Carrera et al. 2015; Yang et al. 2017). We plan to quantify the
impact of these dust concentration processes on planetesimal formation by
conducting hydrodynamical simulations around the edge of the VSI-driven
turbulence zone.

###  4.3 Effects of minimum grain size, grain size distribution, and disk mass

The results presented in section 3.2 depend on the minimum grain size and
slope of the size distribution because the small grains dominate the cooling
rate. In this study, we have fixed the minimum size of grains and slope of
grain size distribution to $a_{\rm min}=1~{}{\rm\micron}$ and $p=-3.5$,
respectively. Smaller dust grows quickly through Brownian motion (Birnstiel et
al., 2011), but its size limit is uncertain, approximately ranging from $0.1$
to $1{\rm~{}\micron}$. As $a_{\rm min}$ decreases, the cooling timescale
decreases, leading to a vertically more extended VSI-unstable layer. The grain
size distribution limited by the radial drift can also become significantly
steeper (e.g., Birnstiel et al. 2011; Stammler & Birnstiel 2022; Birnstiel
2023), with $p$ reaching approximately $-2.5$. As $p$ increases, cooling would
be less efficient because the number of the smallest grains is smaller,
leading to more suppressed the VSI.

Furthermore, variations in the disk mass can alter the cooling rate profile.
In this study, we have fixed the disk mass to $M_{\rm disk}=0.01M_{\odot}$.
The disk mass can depend on disk age (Cazzoletti et al., 2019; Testi et al.,
2022); in particular, the disk mass of the young disk around HL Tau can be
estimated as $M_{\rm disk}\sim 0.1M_{\odot}$ (Kwon et al., 2015). The massive
disk would extend the parameter space where the equilibrium solutions exist
because the larger disk mass can make cooling more efficient, leading to the
radial and vertical expansion of the VSI-unstable layers (see section 4.3 and
figure 10 of Fukuhara et al., 2021). However, as the disk mass increases, the
region with no equilibrium solutions due to the thick stable layer around the
midplane, corresponding to the regions with $R\lesssim 3~{}{\rm au}$ for
figures 5 and 8, would be more extended to larger radial distance. This is
because as the dust density increases, optical depth increases, leading to
inefficient cooling.

###  4.4 Limitations of the model

So far we have assumed the local model in the radial direction, meaning that
gas and dust remain stationary radially. However, VSI-driven turbulence can
diffuse dust radially, altering its spatial distribution (Stoll & Kley, 2016;
Flock et al., 2020; Dullemond et al., 2022). This effect implies that the
edges of the VSI-driven turbulence zone evolve with dust. Quantifying this
effect should be studied in hydrodynamical simulations that include the
dynamic and thermal coupling between gas and dust.

In the region where dust settles in a runaway fashion, other mechanisms that
drive disk turbulence may also contribute to the vertical diffusion of dust
grains. The MRI can drive gas turbulence in both the inner and outer regions
of protoplanetary disks because of the high ionization of gas (for reviews,
Turner et al. 2014; Lesur et al. 2023). The streaming instability caused by
strong dust settling can also trigger turbulent gas motion (Johansen & Youdin,
2007; Yang & Zhu, 2021). Turbulence driven by these mechanisms can maintain a
vertical dust distribution that is balanced between turbulent diffusion and
settling. Moreover, the dust vertical diffusion due to these other mechanisms
may revive VSI-driven turbulence through changes in the vertical profiles of
dust particles and the cooling rate. To understand this, we should investigate
which mechanism dominates turbulence generation in each region of
protoplanetary disks.

Furthermore, this study ignores the effects of dust, magnetic field, and
vertical thermal structure. Dust would increase the effective buoyancy
frequency of the gas (Lin & Youdin, 2017) that prevents the growth of the
linear VSI. The gas–dust drag force can also suppress the VSI and dust
vertical diffusion (Lin, 2019; Lehmann & Lin, 2022, 2023). Moreover, magnetic
fields threading the global disk may suppress the VSI either directly through
magnetic tension or indirectly through MRI turbulence (Nelson et al., 2013;
Latter & Papaloizou, 2018; Cui & Bai, 2020). The roles of magnetic fields in
the VSI suppression can be positive or negative depending on non-ideal
magnetohydrodynamical effects (ambipolar diffusion, Ohmic resistivity, and
Hall effect; Cui & Bai 2020, 2022; Cui & Lin 2021; Latter & Kunz 2022).
Additionally, Zhang et al. (2024) recently found that vertical thermal
stratification with a colder interior and a hotter surface can suppress VSI-
driven turbulence around the midplane. They may change the levels of the
equilibrium vertical dust profile.

##  5 Summary

We have searched for the equilibrium vertical dust profile where settling
balances with diffusion caused by VSI-driven turbulence. We construct the
semi-analytic model that determines the vertical profile of dust grains and
the intensity of VSI-driven turbulence in a self-consistent manner (figure 1).
Our key findings are summarized as follows.

  1. 1.

We find that there exist equilibrium solutions where dust settling balances
with VSI-driven turbulent diffusion for small grains (figure 3). If we assume
that all grains have equal size, there exist two equilibrium solutions when
the single grain size is smaller than $10{\rm~{}\micron}$ at $10{\rm~{}au}$.
If the grain size exceeds $\sim 10{\rm~{}\micron}$, the equilibrium solutions
vanish.

  2. 2.

For the cases of small grains, two equilibrium solutions are classified into
stable and unstable ones [equation (26) and figure 4]. The stable ones
correspond to the dust scale height of $\alpha_{\rm z,equi}\approx 2\times
10^{-3}$, where $\alpha_{\rm z,equi}$ is the equilibrium dimensionless
vertical diffusion coefficient. For the cases of large grains, no equilibrium
solutions indicate runaway dust settling.

  3. 3.

The existence of the equilibrium solutions depends on the radial distance $R$
as well as dust size (figure 5). The equilibrium solutions only exist beyond
$R\sim 3{\rm~{}au}$ because the midplane region at $R\lesssim 3{\rm~{}au}$ is
optically thick, yielding a thick VSI-stable layer (figure 6). The maximum
grain size that allows for the equilibrium solutions also decreases with
increasing $R$. Because the equilibrium solutions require a small grain size,
they lead to a thick dust disk with $H_{\rm d,equi}\approx 0.8H_{\rm g}$
(figure 7), where $H_{\rm d,equi}$ is the dust scale height with the
equilibrium solutions.

  4. 4.

If the particle size distribution is assumed to follow a power law, the small
grains extend the parameter space where strong VSI-driven turbulence is
sustained (figure 8). At $10~{}{\rm au}$, the equilibrium solutions exist up
to $a_{\rm max}=100{\rm~{}\micron}$, where $a_{\rm max}$ is the maximum size
of dust grains. For $a_{\rm max}=10~{}{\rm\micron}$ in the entire disk, VSI-
driven turbulence maintains a thick dust layer of $H_{\rm d,equi}\gtrsim
0.8H_{\rm g}$ at $2{\rm~{}au}\lesssim R\lesssim 50{\rm~{}au}$ (figure 9).

Our results suggest that dust diffusion by VSI-driven turbulence has different
levels depending on the radial distance. This variation may explain the
different degrees of dust settling inferred from observations of some
protoplanetary disks. This implies that VSI-driven turbulence plays a dominant
role in vertical dust diffusion within these disks. Testing this hypothesis
requires a more quantitative investigation of these disks’ cooling rate
structure.

We thank Tomohiro Ono for discussions that motivated this project. We also
thank Akimasa Kataoka, Kiyoaki Doi, Hidekazu Tanaka, Mario Flock, and Takahiro
Ueda for the useful discussions of applications to observed protoplanetary
disks. We appreciate the anonymous referee for comments that greatly helped
improve the manuscript. This work was supported by JSPS KAKENHI Grant Numbers
JP20H01948, JP20H00182, JP22KJ1337, JP23H01227, and JP23K25923.

## References

  * Adachi et al. (1976) Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756
  * Arlt & Urpin (2004) Arlt, R., & Urpin, V. 2004, A&A, 426, 755
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  * Barge & Sommeria (1995) Barge, P., & Sommeria, J. 1995, A&A, 295, L1
  * Barranco et al. (2018) Barranco, J. A., Pei, S., & Marcus, P. S. 2018, ApJ, 869, 127
  * Birnstiel (2023) Birnstiel, T. 2023, arXiv e-prints, arXiv:2312.13287
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79
  * Birnstiel et al. (2011) Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11
  * Brauer et al. (2008) Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859
  * Burke & Hollenbach (1983) Burke, J. R., & Hollenbach, D. J. 1983, ApJ, 265, 223
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Cazzoletti et al. (2019) Cazzoletti, P., Manara, C. F., Liu, H. B., et al. 2019, A&A, 626, A11
  * Chen & Lin (2020) Chen, K., & Lin, M.-K. 2020, ApJ, 891, 132
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
  * Ciesla (2010) Ciesla, F. J. 2010, ApJ, 723, 514
  * Cui & Bai (2020) Cui, C., & Bai, X.-N. 2020, ApJ, 891, 30
  * Cui & Bai (2022) —. 2022, MNRAS, 516, 4660
  * Cui & Lin (2021) Cui, C., & Lin, M.-K. 2021, MNRAS, 505, 2983
  * Doi & Kataoka (2021) Doi, K., & Kataoka, A. 2021, ApJ, 912, 164
  * Draine (2003) Draine, B. T. 2003, ApJ, 598, 1026
  * Drazkowska et al. (2023) Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 717
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  * Dullemond & Dominik (2004) Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075
  * Dullemond et al. (2001) Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957
  * Dullemond et al. (2022) Dullemond, C. P., Ziampras, A., Ostertag, D., & Dominik, C. 2022, A&A, 668, A105
  * Dzyurkevich et al. (2010) Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70
  * Flock et al. (2016) Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ApJ, 827, 144
  * Flock et al. (2017a) —. 2017a, ApJ, 835, 230
  * Flock et al. (2017b) Flock, M., Nelson, R. P., Turner, N. J., et al. 2017b, ApJ, 850, 131
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155
  * Fromang & Papaloizou (2006) Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751
  * Fukuhara et al. (2021) Fukuhara, Y., Okuzumi, S., & Ono, T. 2021, ApJ, 914, 132
  * Fukuhara et al. (2023) —. 2023, PASJ, 75, 233
  * Goldreich & Ward (1973) Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132
  * Hyodo et al. (2021) Hyodo, R., Ida, S., & Guillot, T. 2021, A&A, 645, L9
  * Hyodo et al. (2022) —. 2022, A&A, 660, A117
  * Johansen et al. (2014) Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 547
  * Johansen & Youdin (2007) Johansen, A., & Youdin, A. 2007, ApJ, 662, 627
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75
  * Kataoka et al. (2014) Kataoka, A., Okuzumi, S., Tanaka, H., & Nomura, H. 2014, A&A, 568, A42
  * Kondo et al. (2023) Kondo, K., Okuzumi, S., & Mori, S. 2023, ApJ, 949, 119
  * Kwon et al. (2015) Kwon, W., Looney, L. W., Mundy, L. G., & Welch, W. J. 2015, ApJ, 808, 102
  * Latter & Kunz (2022) Latter, H. N., & Kunz, M. W. 2022, MNRAS, 511, 1182
  * Latter & Papaloizou (2018) Latter, H. N., & Papaloizou, J. 2018, MNRAS, 474, 3110
  * Lehmann & Lin (2022) Lehmann, M., & Lin, M. K. 2022, A&A, 658, A156
  * Lehmann & Lin (2023) Lehmann, M., & Lin, M.-K. 2023, MNRAS, 522, 5892
  * Lesur et al. (2023) Lesur, G., Flock, M., Ercolano, B., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 465
  * Levermore & Pomraning (1981) Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321
  * Li et al. (2001) Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874
  * Li et al. (2000) Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023
  * Lin (2019) Lin, M.-K. 2019, MNRAS, 485, 5221
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17
  * Lin & Youdin (2017) —. 2017, ApJ, 849, 129
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805
  * Lyra & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001
  * Malygin et al. (2017) Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dullemond, C. P. 2017, A&A, 605, A30
  * Manger & Klahr (2018) Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125
  * Manger et al. (2021) Manger, N., Pfeil, T., & Klahr, H. 2021, MNRAS, 508, 5402
  * Melon Fuksman et al. (2024a) Melon Fuksman, J. D., Flock, M., & Klahr, H. 2024a, A&A, 682, A139
  * Melon Fuksman et al. (2024b) —. 2024b, A&A, 682, A140
  * Miotello et al. (2023) Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 501
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Notsu et al. (2022) Notsu, S., Ohno, K., Ueda, T., et al. 2022, ApJ, 936, 188
  * Ohno & Ueda (2021) Ohno, K., & Ueda, T. 2021, A&A, 651, L2
  * Okuzumi & Hirose (2012) Okuzumi, S., & Hirose, S. 2012, ApJ, 753, L8
  * Ormel & Cuzzi (2007) Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413
  * Pfeil et al. (2023) Pfeil, T., Birnstiel, T., & Klahr, H. 2023, ApJ, 959, 121
  * Pfeil & Klahr (2019) Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150
  * Pfeil & Klahr (2021) —. 2021, ApJ, 915, 130
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25
  * Pizzati et al. (2023) Pizzati, E., Rosotti, G. P., & Tabone, B. 2023, MNRAS, 524, 3184
  * Raettig et al. (2021) Raettig, N., Lyra, W., & Klahr, H. 2021, ApJ, 913, 92
  * Richard et al. (2016) Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571
  * Sekiya (1998) Sekiya, M. 1998, Icarus, 133, 298
  * Stammler & Birnstiel (2022) Stammler, S. M., & Birnstiel, T. 2022, arXiv e-prints, arXiv:2207.00322
  * Stoll & Kley (2014) Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77
  * Stoll & Kley (2016) —. 2016, A&A, 594, A57
  * Takeuchi & Lin (2002) Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344
  * Testi et al. (2022) Testi, L., Natta, A., Manara, C. F., et al. 2022, A&A, 663, A98
  * Turner et al. (2014) Turner, N. J., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 411–432
  * Umurhan et al. (2020) Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4
  * Urpin (2003) Urpin, V. 2003, A&A, 404, 397
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399
  * Villenave et al. (2022) Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11
  * Villenave et al. (2023) Villenave, M., Podio, L., Duchêne, G., et al. 2023, ApJ, 946, 70
  * Warren & Brandt (2008) Warren, S. G., & Brandt, R. E. 2008, Journal of Geophysical Research (Atmospheres), 113, D14220
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Yang et al. (2017) Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80
  * Yang & Zhu (2021) Yang, C.-C., & Zhu, Z. 2021, MNRAS, 508, 5538
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
  * Youdin & Shu (2002) Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494
  * Zhang et al. (2024) Zhang, S., Zhu, Z., & Jiang, Y.-F. 2024, arXiv e-prints, arXiv:2404.05608

##  Appendix A Vertical diffusion of dust particles

The evolution equation of the squared mean of the grains’ vertical positions
[equation (26)] is derived by considering the flux of dust in the vertical
direction. The vertical diffusion of dust density $\rho_{\rm d}$ in a
turbulent protoplanetary disk is given by (Dubrulle et al., 1995; Fromang &
Papaloizou, 2006)

$\frac{\partial\rho_{\rm d}}{\partial t}=-\frac{\partial}{\partial
z}\left(\rho_{\rm d}v_{z}\right)+\frac{\partial}{\partial z}\left[\rho_{\rm
g}D_{z}\frac{\partial}{\partial z}\left(\frac{\rho_{\rm d}}{\rho_{\rm
g}}\right)\right],$ (28)

where $v_{z}$ is the vertical velocity of dust and $D_{z}$ is the vertical
turbulent diffusion coefficient of dust. Expanding the second term of the
right side in equation (28) and assuming vertical hydrostatic equilibrium for
gas [equation (3)] yield (Ciesla, 2010)

$\frac{\partial\rho_{\rm d}}{\partial t}=-\frac{\partial}{\partial
z}\left(\rho_{\rm d}v_{z}\right)+\frac{\partial}{\partial
z}\left(D_{z}\frac{\partial\rho_{\rm d}}{\partial
z}\right)+\frac{\partial}{\partial z}\left(D_{z}\frac{z}{H_{\rm
g}^{2}}\rho_{\rm d}\right).$ (29)

Next, we multiply both sides of this equation by $z^{2}$ and integrate it from
$z=-\infty$ to $z=\infty$. The partial integrals are performed repeatedly and
surface terms are dropped, resulting in

$\frac{1}{2}\frac{\partial}{\partial t}\int_{-\infty}^{\infty}\rho_{\rm
d}z^{2}dz=\int_{-\infty}^{\infty}\rho_{\rm
d}v_{z}zdz+D_{z}\left(\int_{-\infty}^{\infty}\rho_{\rm d}dz-\frac{1}{H_{\rm
g}^{2}}\int_{-\infty}^{\infty}\rho_{\rm d}z^{2}dz\right).$ (30)

Assuming the terminal vertical velocity $v_{z}=-{\rm St}_{\rm mid}\Omega_{\rm
K}z$ and a constant dust surface density defined by $\Sigma_{\rm
d}\equiv\int\rho_{\rm d}dz$, equation (30) can be written as

$\frac{1}{2}\frac{d\langle z^{2}\rangle}{dt}=-{\rm St}_{\rm mid}\Omega_{\rm
K}\langle z^{2}\rangle+D_{z}\left(1-\frac{\langle z^{2}\rangle}{H_{\rm
g}^{2}}\right),$ (31)

where $\langle z^{2}\rangle$ is the ensemble average of $z^{2}$ for dust
particles defined by

$\langle z^{2}\rangle\equiv\frac{\displaystyle\int_{-\infty}^{\infty}\rho_{\rm
d}z^{2}dz}{\displaystyle\int_{-\infty}^{\infty}\rho_{\rm d}dz}.$ (32)

If we assume the single-sized model and $\rho_{\rm
d}\propto\exp{\left(-z^{2}/2H_{\rm d}^{2}\right)}$, equation (32) results in
$\langle z^{2}\rangle^{1/2}=H_{\rm d}$, where $H_{\rm d}$ is the dust scale
height.

We use the Stokes number at the midplane ${\rm St}_{\rm mid}$ to calculate the
vertical velocity. In reality, the Stokes number increases as $z$ increases,
causing the small grains to settle toward the midplane faster. However, in
this study, the Stokes numbers of the smallest grains exceed unity only at
$z\gtrsim 4H_{\rm g}$. This height is higher than the height of the VSI-
unstable layer ($z\lesssim 3H_{\rm g}$, see also figure 6 in section 3.1).
Therefore, the effect of ignoring the regime for larger Stokes numbers than
unity has little impact on our results.