# Effects of Dust Evolution on the Vertical Shear Instability in the Outer
Regions of Protoplanetary Disks

Yuya Fukuhara Satoshi Okuzumi Tomohiro Ono Department of Earth and Planetary
Sciences, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan Yuya
Fukuhara<EMAIL_ADDRESS>
###### Abstract

The vertical shear instability (VSI) is a hydrodynamical instability that
requires rapid gas cooling and has been suggested to operate in outer regions
of protoplanetary disks. The VSI drives turbulence with strong vertical
motions, which could regulate the dust growth and settling. However, dust
growth and settling can regulate the VSI because dust depletion makes gas
cooling inefficient in outer disk regions that are optically thin to their own
thermal emission. In this study, we quantify this potentially stabilizing
effects of dust evolution on the VSI based on the linear analysis. We
construct a model for calculating the cooling timescale, taking into account
dust growth beyond micron sizes and size-dependent settling. Combining the
model with the linear stability analysis, we map the region where the VSI
operates, which we call the VSI zone, and estimate the maximum growth rate at
each radial position. We find that dust growth as well as settling makes the
VSI zone more confined around the midplane. This causes a decrease in the
growth rate because the vertical shear of the rotation velocity, which is the
source of the instability, is weaker at lower altitude. In our default disk
model with 0.01 solar masses, dust growth from $10~{}\micron$ to $1{\rm~{}mm}$
causes a decrease in the growth rate by a factor of more than 10. The
suppression of VSI-driven turbulence by dust evolution may promote further
dust evolution in the outer regions and also explain a high degree of dust
settling observed in the disk around HL Tau.

protoplanetary disks — hydrodynamics — instabilities

††journal: ApJ

##  1 Introduction

Planet formation begins with the evolution of dust grains in protoplanetary
disks into kilometer-sized planetesimals. This first stage is initially driven
by the growth of dust grains through mutual sticking and condensation (Chokshi
et al., 1993; Dominik & Tielens, 1997). Large dust particles settle to the
midplane and may experience the streaming instability (Youdin & Goodman, 2005;
Johansen & Youdin, 2007) and gravitational instabilities (Goldreich & Ward,
1973; Youdin, 2011; Takahashi & Inutsuka, 2014; Tominaga et al., 2018, 2019,
2020; Pierens, 2021), which concentrate the dust particles in a runaway
fashion and thereby form planetesimals (e.g., Johansen et al., 2009; Carrera
et al., 2015; Yang et al., 2017). The dust particles may also grow directly
into planetesimals if the particles are sticky enough (e.g., Okuzumi et al.,
2012; Windmark et al., 2012; Kataoka et al., 2013).

Dust evolution in protoplanetary disks depends on gas disk turbulence in many
ways. Turbulence enhances the relative velocity of solid aggregates and may
prevent them from sticking together through collisions (e.g., Brauer et al.,
2008; Okuzumi & Hirose, 2012). Turbulence may also inhibit dust settling
toward the disk midplane and formation of a dense dust layer at the midplane
(e.g., Dubrulle et al., 1995). Therefore, constraining the level of turbulence
in real protoplanetary disks is essential for fully understanding how the dust
in the disks evolves into planetesimals.

Recent radio interferometric observations with the Atacama Large Millimeter-
submillimeter Array have provided us detailed information of dust evolution
and gas turbulence in the outer part of protoplanetary disks. The observations
have provided us with ample evidence that massive and large disks commonly
have rings and gaps of dust (e.g., ALMA Partnership et al., 2015; Andrews et
al., 2018; Long et al., 2018; van der Marel et al., 2019). Although there are
a number of potential mechanisms that provide such substructures (for a
review, see Andrews, 2020), many of them assume that the dust particles
comprising the rings have already grown to $0.1$–10 mm in size so that they
can concentrate radially under the influence of gas drag. Furthermore, the
well separated morphology of the dust rings in the disk around HL Tau (ALMA
Partnership et al., 2015) indicates that the dust particles comprising the
rings have already settled significantly, with a dust scale height being ten
times smaller than the gas scale height (Pinte et al., 2016). Assuming that
the dust rings are indeed dominated by millimeter-sized particles, the high
degree of settling also points to a low level of turbulence near the midplane.
Molecular line emission observations suggest that turbulence in the upper
layers of the outer disk regions is also weak (Flaherty et al., 2015, 2017,
2018, 2020).

The absence of strong turbulence in the outer disk regions is consistent with
the theoretical expectation that the magnetorotational instability (MRI,
Balbus & Hawley, 1991) in the outer regions is suppressed by ambipolar
diffusion (Simon et al., 2013a, b; Bai, 2015; Riols & Lesur, 2018). However,
the suppression of the MRI is not enough to explain the absence of strong
turbulence because purely hydrodynamic disk instabilities can also drive
turbulence (for a review, see Lyra & Umurhan, 2019). Among them, the most
robust one in outer regions of protoplanetry disks is the vertical shear
instability (VSI, Urpin & Brandenburg, 1998; Nelson et al., 2013; Lin &
Youdin, 2015). The VSI is an instability caused by a vertical gradient in the
gas rotation velocity together with a cooling timescale much shorter than the
orbital timescale (Urpin, 2003; Nelson et al., 2013; Lin & Youdin, 2015). The
requirement of rapid gas cooling tends to be met in outer disk regions where
the optical depth is low (Malygin et al., 2017; Pfeil & Klahr, 2019). Once the
VSI operates, it produces turbulence with predominant vertical motion (e.g.,
Stoll & Kley, 2014), which would efficiently prevent vertical dust settling
(Flock et al., 2017, 2020).

The question then is what can suppress the VSI in outer regions of
protoplanetary disks. One candidate is strong magnetic field as suggested by
Nelson et al. (2013) and Cui & Bai (2020). In this paper, we explore the
potential role of dust evolution, i.e., dust growth and settling, in
suppressing the VSI. It is natural to expect that dust evolution should affect
the VSI as it is the dust that is responsible for disk cooling. Malygin et al.
(2017) already noted that a depletion of small grains slows down cooling in
disk regions that are optically thin to their own thermal emission. This
implies that both dust evolution growth and settling should lead to
suppression of the VSI, although no quantitative assessment of the effects has
been made so far.

The goal of this study is to clarify the influences of dust growth and
settling on the stability of outer (5–100 au) protoplanetary disk regions
against the VSI. We calculate the thermal relaxation (cooling) timescale in
the outer disk regions using a parametrized dust model in which the maximum
particle size and dust vertical diffusion coefficient are given as free
parameters. Using the two-dimensional maps of the thermal relaxation timescale
together with linear stability analysis, we study how the location of the VSI-
active region and the VSI growth rate in the region vary as dust grows and
settles.

This paper is organized as follows. In Section 2, we review the basic
properties of the VSI, deriving the linear dispersion relation that gives the
growth rate of the VSI at each location in a disk. We then describe our model
in Section 3, present the main results in Section 4, and discuss limitations
and implications of our study in Section 5. Section 6 presents a summary.

##  2 The VSI

In this section, we review the basic properties of the VSI and derive the
dispersion relation that is used in the following section. The VSI is a type
of the Goldreich–Schubert–Fricke instability (Goldreich & Schubert, 1967;
Fricke, 1968) known in the context of differentially rotating stars. The
presence of a vertical gradient in gas angular velocity $\Omega$ is one of the
necessary conditions for the VSI to operate. The radial and vertical force
balances give

$\Omega^{2}=\frac{GM_{\ast}}{\left(R^{2}+z^{2}\right)^{3/2}}+\frac{1}{R\rho_{g}}\frac{\partial
P}{\partial R},$ (1) $\frac{\partial\left(R\Omega\right)}{\partial
z}=\frac{1}{2\Omega\rho_{g}^{2}}\left(\frac{\partial P}{\partial
z}\frac{\partial\rho_{g}}{\partial R}-\frac{\partial P}{\partial
R}\frac{\partial\rho_{g}}{\partial z}\right),$ (2)

where $R$ is the cylindrical distance from the central star, $z$ is the height
from the midplane, $M_{\ast}$ is the mass of the central star, $P$ is the gas
pressure, $\rho_{g}$ is the gas density, and $G$ is the gravitational
constant. The vertical shear $\partial(R\Omega)/\partial z$ is nonzero if the
temperature gradient exists in the radial direction (Urpin, 2003, see also
Equation (17) in Section 3.1).

###  2.1 The Thermal Criterion for Instability

However, buoyant forces can stabilize the VSI when the entropy increases in
the direction of decreasing gas pressure (Lin & Youdin, 2015). The
Brunt–Väisälä frequency $N_{z}$ is given by

$N_{z}^{2}\equiv-\frac{1}{\rho_{g}C_{P}}\cdot\frac{\partial P}{\partial
z}\cdot\frac{\partial s}{\partial z},$ (3)

where $C_{P}$ and $s$ are the specific heat at constant pressure and the
specific entropy, respectively. The specific entropy is given by
$s=C_{V}\log{\left(P/\rho_{g}^{\gamma}\right)}$, where $C_{V}$ is the specific
heat at constant volume and $\gamma$ is the heat capacity ratio. The vertical
buoyancy is stabilizing if $N_{z}^{2}>0$. In protoplanetary disks, the
pressure decreases with distance $z$ from the midplane. In outer disk regions
where the temperature is determined by stellar irradiation (e.g., Chiang &
Goldreich, 1997), the entropy increases with $z$, thus stabilizing the VSI.

Therefore, the VSI requires fast thermal relaxation that reduces buoyancy
(Nelson et al., 2013). This requirement can be expressed as

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

Here, $\tau_{{\rm crit}}$ is the critical thermal relaxation timescale defined
by (Lin & Youdin, 2015)

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

where $\Omega_{\rm K}$ is the Keplerian frequency, $H_{g}$ is the gas scale
height, and $q$ is the radial gradient of the temperature. Strictly speaking,
high-$k_{x}$ unstable modes persist for $\tau_{\rm relax}\gtrsim\tau_{\rm
crit}$, but we neglect these modes because their growth rates are much smaller
than the maximum growth rates for $\tau_{\rm relax}\lesssim\tau_{\rm crit}$.
We describe more details of the thermal relaxation timescale $\tau_{{\rm
relax}}$ in Section 3.2.1.

In realistic protoplanetary disks, Equation (4) is fulfilled in regions around
the midplane (Malygin et al., 2017; Pfeil & Klahr, 2019, see also Section 4 of
this paper). In this study, we refer to such a region as a VSI zone.

###  2.2 Linear Analysis

A number of previous studies already derived dispersion relations with and
without gas cooling and vertical stratification (Urpin & Brandenburg, 1998;
Urpin, 2003; Arlt & Urpin, 2004; Nelson et al., 2013; Barker & Latter, 2015;
Lin & Youdin, 2015, 2017; Latter & Papaloizou, 2018). Here, we follow Nelson
et al. (2013, 2016) and derive a local dispersion relation.

We consider a VSI zone with the finite vertical extent and focus on linear
modes whose wavelengths are short enough to fit into the zone (see Section
3.2.2 for the allowed range of wavenumbers). We apply a locally isothermal
equation of state to the VSI zone and regard the modes within the zone as
radially and vertically local. We note that Lin & Youdin (2015) presents a
dispersion relation for vertically global VSI modes including the effect of
finite thermal relaxation. However, we do not use this here because it is not
strictly applicable to disks with vertically varying $\tau_{\rm relax}$.

Under the local shearing box approximation, the equation of continuity can be
written as

$\frac{\partial\rho_{g}}{\partial t}+\nabla\cdot\left(\rho_{g}{v}\right)=0,$
(6)

where $\rho_{g}$ is the gas density and ${v}=(v_{x},~{}v_{y},~{}v_{z})$ is the
gas velocity with three components for the radial, azimuthal and vertical
velocities in the Cartesian shearing box. The equations of motion of gas in
the shearing box are

$\frac{\partial v_{x}}{\partial
t}+\left({v}\cdot\nabla\right)v_{x}=-\frac{1}{\rho_{g}}\frac{\partial
P}{\partial x}+2\Omega_{0}v_{y}+3\Omega_{0}^{2}x,$ (7) $\frac{\partial
v_{y}}{\partial
t}+\left({v}\cdot\nabla\right)v_{y}=-\frac{1}{\rho_{g}}\frac{\partial
P}{\partial y}-2\Omega_{0}v_{x},$ (8) $\frac{\partial v_{z}}{\partial
t}+\left({v}\cdot\nabla\right)v_{z}=-\frac{1}{\rho_{g}}\frac{\partial
P}{\partial z}-g,$ (9)

where $\Omega_{0}=\mathrm{constant}$ is the angular velocity in the shearing
box, and $g=\Omega_{0}^{2}z$ is the vertical component of stellar gravity.
Following Goldreich & Schubert (1967), we consider the vicinity of a given
height $z=z_{0}$ and assume that $g$ and the isothermal sound speed
$c_{s}=\sqrt{P/\rho_{g}}$ are approximately constant on this small vertical
scale.

We assume axisymmetry and consider perturbations of the form $\propto
e^{-i\omega t+ik_{x}x+i(k_{z}-ig/(2c_{s}^{2}))z}$ on the steady background.
Here, $\omega$ is the angular frequency and $k_{x}$ and $k_{z}$ are the radial
and vertical wavenumbers, respectively. The factor $e^{gz/(2c_{s}^{2})}$
accounts for the variation of the background density in the vertical direction
(Nelson et al., 2013), The dispersion relation for the perturbations is (see
Nelson et al. (2013) and Nelson et al. (2016) for a derivation)

$\displaystyle\omega^{4}-\left[c_{s}^{2}\left(k_{x}^{2}+k_{z}^{2}\right)+\kappa_{0}^{2}+{\frac{g^{2}}{4c_{s}^{2}}}\right]\omega^{2}$
$\displaystyle-2\Omega_{0}c_{s}^{2}k_{x}k_{z}\frac{\partial v_{y0}}{\partial
z}+i\Omega_{0}gk_{x}\frac{\partial v_{y0}}{\partial z}+$
$\displaystyle\kappa_{0}^{2}\left(c_{s}^{2}k_{z}^{2}+{\frac{g^{2}}{4c_{s}^{2}}}\right)=0,$
(10)

where $v_{y0}$ and $\kappa_{0}$ are the rotation velocity deviation and
epicyclic frequency, respectively, for the background; the latter one can be
written as

$\kappa_{0}^{2}=2\Omega_{0}\left(\frac{\partial v_{y0}}{\partial
x}+2\Omega_{0}\right).$ (11)

For an unstable perturbation that satisfies Equation (2.2), the growth rate is
given by

$\Gamma_{{\rm VSI}}=\mathrm{Im}\left(\omega\right).$ (12)

##  3 Model

We study the VSI linear stability of protoplanetary disks taking into account
dust growth and settling. We present a gas disk model in Section 3.1, a model
for calculating the disk relaxation time and VSI growth rate in Section 3.2, a
dust model used to calculate the thermal relaxation time in Section 3.3, and
our parameter choice in Secrion 3.4.

###  3.1 Gas Disk Model

Figure 1: Local rotational angular velocity $\Omega$ (upper panel) and
vertical shear $\partial\left(R\Omega\right)/\partial z$ (lower panel), as a
function of $R$ and $z/R$. The dashed lines represent $H_{g}$, $2H_{g}$, and
$3H_{g}$ in height from the midplane.

We consider an axisymmetric disk around a solar-mass star. The gas surface
density is given by

$\Sigma_{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]},$
(13)

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. Equation (13) is
motivated by the similarity solution of the viscous accretion disk model
(Lynden-Bell & Pringle, 1974; Hartmann et al., 1998) although we do not
consider disk evolution in this study. Since we are primarily interested in
the effects of dust growth and settling on the VSI, we fix the gas disk
parameters to $M_{{\rm disk}}=0.01M_{\sun}$, $R_{c}=100~{}{\rm au}$, and
$\beta_{\Sigma}=1$ except in Section 4.3, where we show that disks of higher
disk dust masses have more extended VSI zones. The heat capacity ratio
$\gamma$ and mean molecular mass $m_{g}$ of the gas is taken to be 1.4 and
$2.3m_{p}$, respectively, where $m_{p}$ is the proton mass.

We focus on the outer region of the disk where the temperature is determined
by stellar irradiation. Assuming that the disk is optically thick to stellar
radiation 111 A disk region that is optically thin to its own thermal emission
can be optically thick to the radiation from the central star because the
opacity in the visible is higher than that in the infrared and, more
importantly, because the radial optical depth is $\sim R/H_{g}$ ($\sim
10$–100) times larger than the vertical optical depth (Chiang & Goldreich,
1997). and that the stellar luminosity is equal to the solar luminosity, the
temperature of the disk interior is given by

$T(R)=130\left(\frac{R}{1{\rm~{}au}}\right)^{q}{\rm~{}K}$ (14)

with $q=-3/7$ (Chiang & Goldreich, 1997). Viscous heating is negligible as
long as we focus on $R\gtrsim 10~{}\rm au$ (e.g., Bitsch et al., 2015). We
assume that the disk interior is vertically isothermal, neglecting warmer
surface layers that are optically thin to the starlight (Chiang & Goldreich,
1997).

From vertical hydrostatic equilibrium, the gas density is given by

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

with $H_{g}=c_{s}/\Omega_{\rm K}$, $\Omega_{\rm K}=\sqrt{GM_{\ast}/R^{3}}$,
and $M_{*}=1M_{\sun}$.

The vertical shear of the gas rotation velocity, $\partial v_{y0}/\partial
z=\partial(R\Omega)/\partial z$, characterizes the strength of the VSI.
Because we assume a radially varying temperature profile,
$\partial(R\Omega)/\partial z$ is nonzero at all height except at the
midplane. This can be analytically shown for $R\ll R_{c}$, where
$\Sigma_{g}\propto R^{-\beta_{\Sigma}}$. In this region, $\Omega$ and
$\partial(R\Omega)/\partial z$ can be approximated as (Takeuchi & Lin, 2002)

$\Omega\left(R,~{}z\right)\approx\Omega_{\rm
K}\left[1+\frac{1}{2}\left(\frac{H_{g}}{R}\right)^{2}\left(-\beta_{\Sigma}{+\frac{q-3}{2}}+\frac{q}{2}\frac{z^{2}}{H_{g}^{2}}\right)\right],$
(16)

and

$\frac{\partial\left(R\Omega\right)}{\partial
z}\approx\frac{q}{2}\frac{z}{R}\Omega_{\rm K},$ (17)

respectively. Equation (17) proves that a nonzero radial temperature gradient
is the source of the vertical shear (e.g., Lin & Youdin, 2015). Figure 1 shows
the gas angular velocity $\Omega(R,z)$ and the vertical shear of the angular
velocity $\partial(R\Omega)/\partial z$ in the disk, including the region
$R\gtrsim R_{c}$. From the lower panel of Figure 1, and also from Equation
(17), the vertical shear increases with $z$, suggesting that the VSI is
stronger at higher altitude as long as the cooling criterion (Equation (4)) is
fulfilled.

We assume that the disk is weakly turbulent and express the turbulent
diffusion coefficient as $\alpha_{D}c_{s}H_{g}$, where $\alpha_{D}$ is the
dimensionless parameter characterizing the level of turbulent diffusion.
Turbulent diffusion controls the maximum wavenumber of the VSI modes (Section
3.2.2) and the vertical scale height of dust particles (Setion 3.3). In
principle, $\alpha_{D}$ in outer disk regions should depend on the strength of
VSI-driven turbulence, and hence on the VSI growth rate $\Gamma_{{\rm VSI}}$,
which is the output of our model. Therefore, a self-consistent determination
of $\alpha_{D}$ and $\Gamma_{{\rm VSI}}$ requires a model that predicts the
former as a function of the latter. Lacking such a model, we opt for taking
$\alpha_{D}$ as a free parameter, although we do discuss potential feedback of
dust settling on the level of VSI-driven turbulence in Section 5.2.

###  3.2 Instability Analysis

We analyze the linear stability of the model disk against the VSI in the
following three steps. In the first step, we use the thermal relaxation
criterion (Equation (4)) to search for the disk region, which we call the VSI
zone, where the linear VSI operates (Section 3.2.1). In the second step, we
compute the range of wavenumbers for the VSI modes that fit into the active
zone (Section 3.2.2). In the third step, we calculate the maximum growth rate
of the VSI at each radial location (Section 3.2.3).

####  3.2.1 Defining the VSI Zone

We compute the radial and vertical extent of the VSI zone by applying the
thermal relaxation criterion (Equation (4)) to each point ($R$, $z$) in the
disk. In protoplanetary disks, local thermal relaxation is regulated either by
collisional heat transfer from gas to dust or by radiative cooling (Malygin et
al., 2017). The area around the boundary of the VSI zone can be regarded as
optically thin (Malygin et al., 2017). Following Malygin et al. (2017) and
Pfeil & Klahr (2019), we estimate the local thermal relaxation timescale
$\tau_{{\rm relax}}$ as

$\tau_{{\rm relax}}=\max{\left(\tau_{{\rm coll}},\tau_{{\rm emit}}\right)},$
(18)

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

The timescale of collisional heat transfer is given by

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

where $\ell_{{\rm gd}}$ is the mean travel length of gas molecules before
colliding with dust particles and $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}}\approx\sqrt{\frac{8k_{\rm B}T}{\pi m_{g}}},$ (20)

where $k_{\rm B}$ is the Boltzmann constant. 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)

where $\kappa_{\rm P}(T)$ is the Planck mean opacity per unit gas mass and
$\sigma_{{\rm SB}}$ is the Stefan–Boltzmann constant. Both $\ell_{{\rm gd}}$
and $\kappa_{\rm P}$ depend on the local size distribution of the dust
particles. Our dust model is described in Section 3.3.

In general, collisional heat transfer regulates the cooling timescale at high
altitude where the dust density is low (Malygin et al., 2017). Because the
dust density decreases monotonically with $|z|$, there is the height $|z|=L$
above which the VSI is stable. In other words, the VSI zone refers to the
region where $|z|<L$. Note that $L$ generally depends on $R$ and becomes zero
where the VSI is stable at all heights.

####  3.2.2 Wavenumber Range of the VSI Modes

The VSI modes generally have $k_{x}k_{z}<0$ (Arlt & Urpin, 2004), so we
restrict $k_{x}>0$ and $k_{z}<0$. Below we further restrict the wavelength
range permitted for VSI modes by accounting for viscous damping and the finite
thickness of the VSI zone.

Viscous damping erases unstable modes of short wavelengths, giving upper
limits on $k_{x}$ and $|k_{z}|$ for the true VSI modes. We particularly focus
on the upper limit on the radial wavenumbers because VSI modes typical have
$k_{x}\gg|k_{z}|$ (Arlt & Urpin, 2004). The maximum wavenumber set by viscous
damping can be estimated as (Lin & Youdin, 2015)

$(k_{x,{\rm max}}H_{g})^{2}\approx\frac{|q|}{\alpha_{D}}\frac{H_{g}}{R}$ (22)

or equivalently,

$k_{x,{\rm max}}\approx\sqrt{\frac{|q|}{\alpha_{D}RH_{g}}}.$ (23)

Because the VSI zone has a finite vertical extent, only modes whose vertical
wavelengths are short enough to “fit” into the zone can be unstable.
Specifically, we require the VSI modes at each $R$ to have half wavelengths
shorter than the vertical zone width $2L$. In other words, the vertical
wavenumbers of the VSI modes must satisfy $|k_{z}|>|k_{z}|_{\rm min}$, where

$|k_{z}|_{\rm min}\approx\frac{\pi}{2L}.$ (24)

####  3.2.3 The Maximum Growth Rate

The local dispersion relation (Equation (2.2)) formally gives the growth rate
of the VSI at each ($R$, $z$). However, the most unstable VSI modes typically
have vertical wavelengths comparable to $H_{g}$ (e.g., Nelson et al., 2013,
see also Section 4.2). Physically, such long-wavelength modes should be
regarded as extending over the whole vertical extent of the VSI zone.

For this reason, we opt for a semi-local approach in which we calculate the
maximum VSI growth rate $\Gamma_{\rm VSI,maax}$ at each $R$ using the linear
dispersion relation, but limiting the range of vertical wavenumbers to
$|k_{z}|\geqslant|k_{z}|_{\rm min}$ to account for the finite vertical extent
of the VSI zone. The vertically local quantities $g(=-\Omega_{\rm K}^{2}z)$
and $v_{y0}(=R\Omega(R,z))$ involved in the local dispersion relation are
evaluated at the vertical boundary of the VSI zone, $|z|=L$, because the VSI
is generally more vigorous at higher altitude unless buoyancy suppresses it
(Nelson et al., 2013). We also limit radial wavenumbers to $k_{x}\leqslant
k_{x,{\rm max}}$ to account for viscous damping.

Besides, the local analysis in this study and the vertically global analysis
by Lin & Youdin (2015) give similar predictions for the body mode in a
vertically wide VSI zone (see Section 5.1).

###  3.3 Dust Model

We here describe the dust model we use to calculate $\ell_{\rm dg}$ and
$\kappa_{\rm P}$. We consider spherical, icy dust particles of internal
density $\rho_{{\rm int}}=1{\rm~{}g~{}cm^{-3}}$. The ratio between the dust
surface density $\Sigma_{d}$ and $\Sigma_{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. As we discuss in Section 4.3, the VSI
zone shrinks as $\Sigma_{d}$ decreases.

The particle size distribution is assumed to follow a power law

$\displaystyle\frac{dN_{d}(a)}{da}=\left\\{\begin{array}[]{ll}\frac{3\Sigma_{d}}{8\pi\rho_{\rm
int}\left(\sqrt{a_{{\rm max}}}-\sqrt{a_{{\rm min}}}\right)}a^{-7/2},&a_{{\rm
min}}<a<a_{{\rm max}},\\\ 0,&{\rm otherwise},\end{array}\right.$ (27)

where $dN_{d}(a)/da$ is the number surface density per unit particle size $a$,
$\Sigma_{d}$ is the total dust mass surface density, and $a_{{\rm min}}$ and
$a_{{\rm max}}$ are the minimum and maximum particle sizes, respectively. The
size distribution given by Equation (27) satisfies the normalization

$\Sigma_{d}=\int_{a_{{\rm min}}}^{a_{{\rm max}}}m\frac{dN_{d}(a)}{da}da,$ (28)

where $m=(4\pi/3)\rho_{\rm int}a^{3}$ is the particle mass. The power-law
slope assumed in Equation (27) is simply taken from the interstellar particle
size distribution (Mathis et al., 1977). We note, however, that the size
distribution can be somewhat shallower or steeper than assumed in Equation
(27) depending on the details of collisional growth and fragmentation
(Birnstiel et al., 2011).

We take the maximum particle size $a_{{\rm max}}$ as a free parameter to study
the impact of dust growth on the VSI. The minimum particle size is less well
defined, but we may crudely taken it to be $\sim$ 0.1–1 $\micron$ because
particles smaller than these sizes grow quickly through Brownian motion
(Birnstiel et al., 2011). We fix $a_{{\rm min}}=1~{}\micron$ throughout this
study.

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

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

where $dn_{d}(a,z)/da$ is the particle number density per unit radius at
height $z$, $\mathrm{St}_{{\rm mid}}$ is the Stokes number of the particles at
the midplane, and $C_{d}(a)$ is the normalized constant determined by the
condition

$\frac{dN_{d}(a)}{da}=\int\frac{dn_{d}(a,z)}{da}dz.$ (30)

The Stokes number is the product of the stopping time and Keplerian frequency.
Assuming that the particle radius are 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}}=\frac{\pi\rho_{{\rm int}}a}{2\Sigma_{g}}.$ (31)

To evaluate $C_{d}$, we note that the integral in Equation (30) is dominated
by the region $z\lesssim H_{g}$, for which the Equation (29) can be
approximated as ${dn_{d}(a,z)}/{da}\approx C_{d}(a)\exp(-z^{2}/(2H_{d}^{2}))$,
where

$H_{d}=\left(1+\frac{\mathrm{St}_{{\rm mid}}}{\alpha_{D}}\right)^{-1/2}H_{g}$
(32)

represents the scale height of particles with size $a$ (Dubrulle et al., 1995;
Youdin & Lithwick, 2007). This approximation gives

$C_{d}(a)=\frac{1}{\sqrt{2\pi}H_{d}}\frac{dN_{d}(a)}{da}.$ (33)

The vertical–size distribution $dn_{d}(a,z)/da$ gives $\ell_{{\rm gd}}$ and
$\kappa_{\rm P}$ 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_{a_{{\rm min}}}^{a_{{\rm max}}}\pi
a^{2}\frac{dn_{d}}{da}da\right)^{-1}.$ (34)

The Planck mean opacity $\kappa_{\rm P}$ per unit gas mass can be written as

$\kappa_{\rm P}(T)=\frac{1}{\rho_{g}}\int_{a_{{\rm min}}}^{a_{{\rm
max}}}\sigma_{\rm abs,P}(a,T)\frac{dn_{d}}{da}da,$ (35)

where $\sigma_{\rm abs,P}(a,T)$ is the Planck mean absorption cross section of
the particles. In this study, we approximate $\sigma_{\rm abs,P}$ with the
monochromatic absorption cross section $\sigma_{\rm abs}(a,\lambda)$ at
wavelength $\lambda=\lambda_{\rm peak}(T)$, where

$\lambda_{{\rm peak}}(T)=10{\rm~{}\mu m}\left(\frac{300{\rm~{}K}}{T}\right)$
(36)

is the peak wavelength of the Planck function. Furthermore, we crudely
approximate $\sigma_{\rm abs}(a,\lambda)$ as

$\sigma_{{\rm abs}}(a,\lambda)=\pi a^{2}\min\left(1,\frac{2\pi
a}{\lambda}\right),$ (37)

where the factor $2\pi a/\lambda$ is called the size parameter in Mie
scattering theory (see, e.g., Bohren & Huffman, 1983). This approximate
expression satisfies the asymptotic properties of Mie scattering that
$\sigma_{{\rm abs}}$ approaches the geometric cross section $\pi a^{2}$ in the
short wavelength limit ($\lambda\ll 2\pi a$) and that the dust mass opacity
$\sigma_{{\rm abs}}/m$ is independent of $a$ in the long wavelength limit
($\lambda\gg 2\pi a$).

We neglect the radial inward drift of the dust particles due to gas drag
(Whipple, 1972; Adachi et al., 1976; Weidenschilling, 1977) and assume that
the dust surface density profile is constant in time. This approach is valid
if the VSI growth timescale is shorter than the dust drift timescale. For
${\rm St}_{\rm mid}<1$, which is the case for the dust particles considered in
our model, the drift timescale can be estimated as $\sim(R/H_{g})^{2}{\rm
St}_{\rm mid}^{-1}\Omega_{\rm K}^{-1}$ (see the references listed above). In
the limit of short relaxation timescales, the VSI grows on a timescale of
$\sim(R/H_{g})\Omega_{\rm K}^{-1}$ (e.g., Urpin & Brandenburg, 1998; Nelson et
al., 2013), which is shorter than the drift timescale as long as ${\rm
St}_{\rm mid}<1$. However, radial dust drift may not be negligible when a
finite relaxation time suppresses the growth of the VSI. In this case, dust
depletion due to the radial inward drift would further suppress the VSI as we
demonstrate in Section 4.3.

###  3.4 Computational Domain and Parameter Choices

We consider a disk region defined by $5{\rm~{}au}<R<100{\rm~{}au}$ and
$0<z/R<0.2$ and divide the region into 1000 logarithmically spaced radial
grids and 1000 linearly spaced vertical grids. The dust size distribution is
divided into logarithmic bins of 10 grids per decade in $a$. The main free
parameters in our model are the vertical diffusion coefficient $\alpha_{D}$
and maximum dust particle size $a_{{\rm max}}$. We take
$\alpha_{D}=10^{-3},10^{-4}$, and $10^{-5}$ and $a_{{\rm
max}}=10~{}\micron,100~{}\micron$, and $1{\rm~{}mm}$.

##  4 Results

Figure 2: Collisional timescale $\tau_{{\rm coll}}$ normalized by radiative
cooling timescale $\tau_{{\rm emit}}$ (upper panel) and critical thermal
relaxation timescale $\tau_{{\rm crit}}$ (lower panel) as a function of $R$
and $z/R$ for $\alpha_{D}=10^{-4}$ and $a_{{\rm max}}=10~{}\micron$. The
thermal relaxation timescale is given by the larger of the two (see Equation
(18)). The dashed lines represent $z=H_{g}$, $2H_{g}$, and $3H_{g}$. The solid
line marks $\tau_{{\rm coll}}=\tau_{{\rm crit}}$.

In this section, we use the model presented in Section 3 to study how dust
growth and settling affects the VSI stability of protoplanetary disks. We map
the VSI zones in Section 4.1 and then compute the radial distribution of the
growth rate and wavenumbers of the most unstable VSI model in Section 4.2. We
study the dependence of the VSI zone size on the dust and gas surface
densities in Section 4.3.

###  4.1 Radial and Vertical Extent of the VSI Zone

Figure 3: Location of the VSI zone (the shaded area) for different values of
$\alpha_{D}$ with $a_{{\rm max}}=100~{}\micron$. The dashed lines represent
$H_{g}$, $2H_{g}$, and $3H_{g}$ in height from the midplane. Figure 4: Mean
travel length of gas molecules before colliding with dust particles,
$\ell_{{\rm gd}}$, as a function of $R$ and $z/R$ for different values of
$\alpha_{D}$ with $a_{{\rm max}}=100~{}\micron$. The dashed lines represent
$H_{g}$, $2H_{g}$, and $3H_{g}$ in height from the midplane. The solid lines
correspond to the boundary of the VSI zone.

As described in Section 2, the VSI is active where the thermal relaxation time
$\tau_{{\rm relax}}$ (Equation (18)) is shorter than the critical timescale
$\tau_{\rm crit}$ (Equation (5)). We find that $\tau_{{\rm relax}}$ (Equation
(18)) is determined by the collisional cooling timescale $\tau_{{\rm coll}}$
for all parameters and all regions explored in this study. Figure 2 plots
$\tau_{\rm coll}/\tau_{\rm coll}$ and $\tau_{\rm coll}/\tau_{\rm crit}$ as a
function of $R$ and $z$ for $\alpha_{D}=10^{-4}$ and $a_{{\rm
max}}=10~{}\micron$, showing that $\tau_{{\rm coll}}$ is about an order-of-
magnitude larger than $\tau_{{\rm emit}}$ at all locations. Our result is
consistent with that of Malygin et al. (2017), who showed that $\tau_{{\rm
coll}}>\tau_{{\rm emit}}$ in an optically thin region away from the central
star (see their Figure 3). The maps of $\tau_{\rm relax}$ (in units of
$\Omega_{\rm K}^{-1}$) for all parameter sets are shown in Figures 11 and 12
in Appendix A.

The solid line in the lower panel of Figure 2 indicates the boundary of the
the VSI zone; below this line, one has $\tau_{{\rm relax}}(=\tau_{{\rm
coll}})<\tau_{{\rm crit}}$ and the VSI operates.222 To be precise, the VSI
zone does not include the midplane ($z=0$), where the vertical shear vanishes.
In this example, the VSI zone extends to $z\sim$0.5–$2H_{\rm g}$ at $5~{}{\rm
au}<R<30~{}{\rm au}$. The vertical width of the VSI zone diminishes as $R$
increases. The vertical optical depth from infinity to the VSI zone boundary
is $O(1)$ at $R\sim 10~{}\rm au$ and decreases with increasing $R$, so using
the thermal relaxation timescale for the optically thin limit (Section 3.2.1)
is marginally justified.

Figure 5: Location of the VSI zone (the shaded area) for different values of
$a_{{\rm max}}$ with $\alpha_{D}=10^{-4}$. The dashed lines represent $H_{g}$,
$2H_{g}$, and $3H_{g}$ in height from the midplane. Figure 6: Mean travel
length of gas molecules before colliding with dust particles, $\ell_{{\rm
gd}}$, as a function of $R$ and $z/R$ for different values of $a_{{\rm max}}$
with $\alpha_{D}=10^{-4}$. The dashed lines represent $H_{g}$, $2H_{g}$, and
$3H_{g}$ in height from the midplane. The solid lines correspond to the
boundary of the VSI zone.

Figure 3 indicates the location of the VSI zone for $a_{\rm max}=100~{}\rm\mu
m$ but with different values of $\alpha_{D}$. This figure illustrates how dust
settling affects the extent of the VSI zone; in our model, the dust scale
height decreases with decreasing $\alpha_{D}$ (see Equation (32)). We find
that dust settling leads to a VSI zone that is more confined to the midplane
region and more extended to larger radial distances. This is because the
settling causes dust depletion and dust concentration, which increase and
decrease the thermal relaxation timescale $\tau_{\rm relax}=\tau_{\rm
coll}\propto\ell_{\rm gd}$, at high and low altitude, respectively (see also
Figure 11). See Figure 4 for the variation of $\ell_{\rm gd}$ at different
locations with $\alpha_{D}$.

The stabilizing effect of dust setting on the VSI at the midplane can also be
confirmed by looking at how $\ell_{\rm gd}$ at the midplane depends on
$\alpha_{\rm D}$. For simplicity, we assume $\alpha_{D}<{\rm St}_{{\rm mid}}$,
which holds in our model at sufficiently large $R$ ($R\gtrsim 50~{}{\rm au}$
for $\alpha_{D}=10^{-4}$). With this assumption, we approximate
$H_{d}\sim(\alpha_{D}/\mathrm{St}_{{\rm mid}})^{1/2}H_{g}$ and
$dn_{d}(z=0)/da\propto({\rm St}_{\rm mid}/\alpha_{D})^{1/2}dN_{d}/da$. Using
this, the integration in Equation (34) can be analytically performed,
resulting in

$\ell_{{\rm gd}}(z=0)\propto\frac{\Sigma_{\rm g}^{1/2}}{\Sigma_{\rm
d}}\frac{a_{{\rm max}}^{1/2}}{\log\left(a_{{\rm max}}\right)}\alpha_{D}^{1/2}$
(38)

for $a_{{\rm max}}\gg a_{{\rm min}}$. Equation (38) confirms that $\ell_{{\rm
gd}}$ at the midplane decreases with decreasing $\alpha_{D}$ (see Figure 4).

Figure 5 indicates the location of the VSI zone for different values of
$a_{\rm max}$, illustrating how the VSI zone evolves with dust growth. The
figure shows that the VSI zone shrinks toward the midplane and toward the
central star as $a_{{\rm max}}$ increases. This is because increasing $a_{\rm
max}$ decreases dust particles’ total surface area and thus increases
$\ell_{\rm gd}$ as shown in Figure 6. In fact, increasing $a_{\rm max}$ also
promotes dust settling, which acts to decrease $\ell_{\rm gd}$ at the
midplane. However, we find that this effect is minor compared to the increase
in the midplane $\ell_{\rm gd}$ due to local dust growth. This can also be
confirmed from Equation (38) implying that $\ell_{\rm gd}$ at the midplane
increases with $a_{{\rm max}}$. Well above the midplane, both dust settling
and local dust growth increase $\ell_{{\rm gd}}$, stabilizing the VSI. A
closer inspection shows that local dust growth dominates the increase of
$\ell_{\rm gd}$ at $z\lesssim 1$–$2H_{g}$.

###  4.2 The Maximum Growth Rate and Wavenumbers

Figure 7:  VSI growth rate $\Gamma_{{\rm VSI}}$ as a function of $k_{x}$ and
$|k_{z}|$ at the top of the VSI zone at $R=60~{}{\rm au}$ for $a_{{\rm
max}}=10~{}\micron$ and $\alpha_{D}=10^{-4}$. The vertical and horizontal
dotted lines indicate $k_{x}=k_{x,{\rm max}}$ and $k_{z}=|k_{z}|_{\rm min}$,
respectively. The star symbol marks the most unstable VSI mode in the range
$k_{x}<k_{x,\max}$ and $|k_{z}|>|k_{z}|_{\rm min}$, which is limited by
viscous damping and the finite vertical extent of the VSI zone. For
comparison, the filled circle marks the mode that would be most unstable if
the ranges of $k_{x}$ and $k_{z}$ were not limited.

We search for the most unstable VSI mode at each $R$ using the procedure
described in Section 3.2.3. For the entire parameter space explored in this
study, we find that the must VSI mode always lies at $k_{x}=k_{x,\rm max}$ and
$|k_{z}|=|k_{z}|_{\rm min}$. This is illustrated in Figure 7, where we plot
the VSI growth rate $\Gamma_{{\rm VSI}}$ as a function of $k_{x}$ and $k_{z}$
at the top of the VSI zone at $R=60~{}{\rm au}$ in the case of $a_{{\rm
max}}=10~{}\micron$ and $\alpha_{D}=10^{-4}$. In this example, we have
$k_{x,\rm max}\approx 18/H_{g}$ and $|k_{z}|_{\rm min}\approx 1.7/H_{g}$
($L\approx 0.9H_{g}$). In the range of $k_{x}\leqslant k_{x,\rm max}$ and
$|k_{z}|\geqslant|k_{z}|_{\rm min}$, $\Gamma_{\rm VSI}$ reaches a maximum of
$0.003\Omega_{\rm K}$ at $k_{x}=k_{x,\rm max}$ and $|k_{z}|=|k_{z}|_{\rm min}$
(marked by the star point in Figure 7). If we did not limit $k_{x}$ and
$|k_{z}|$, a higher maximum growth rate of $0.015\Omega_{\rm K}$ would be
reached at the point $k_{x}\approx 100/H_{g}$ and $|k_{z}|\approx 1.5/H_{g}$
($L\approx H_{g}$) marked by the filled circle in Figure 7.

Figure 8:  Maximum radial wavenumber $k_{x,\max}$ (dotted lines; Equation
(23)) and minimum vertical wavenumber $|k_{z}|_{\rm min}$ (solid lines;
Equation (24)) as a function $R$ for different values of $\alpha_{D}$ and
$a_{{\rm max}}$. These are equal to the wavenumbers of the most unstable VSI
mode (see Section 4.2).

Figure 8 shows the wavenumbers of the most unstable VSI modes, $k_{x,\max}$
and $|k_{z}|_{\rm min}$, as a function of $R$ for various values of
$\alpha_{D}$ and $a_{\rm max}$. Note that $k_{x,\max}$ depends only on
$\alpha_{D}$ whereas $|k_{z}|_{\rm min}$ depends on both $\alpha_{D}$ and
$a_{\rm max}$. Being inversely proportional to the VSI zone vertical thickness
$L$, $|k_{z}|_{\rm min}$ increases as $a_{{\rm max}}$ increases or
$\alpha_{D}$ decreases. Our assumption $|k_{z}|\ll k_{x}$ breaks down near the
outer edge of the VSI zone where $|k_{z}|_{\rm min}$ diverges. However, this
region is narrow compared to the VSI zone itself.

Figure 9: Maximum VSI growth rate $\Gamma_{{\rm VSI,max}}$ as a function of
$R$ for different values of $\alpha_{D}$ with $a_{{\rm max}}=10~{}\micron$,
$100~{}\micron$, and $1{\rm~{}mm}$.

Figure 9 shows the maximum VSI growth rate $\Gamma_{\rm VSI,max}$ as a
function of $R$ for various values of $\alpha_{\rm D}$ and $a_{\rm max}$.
Overall, $\Gamma_{\rm VSI,max}$ decreases with increasing $a_{\rm max}$,
reflecting the fact that the VSI zone shrinks as $a_{\rm max}$ increases. In
this default disk model of $M_{\rm disk}=0.01M_{\sun}$, an increase in $a_{\rm
max}$ from 10 $\micron$ to 100 $\micron$ causes a decrease in $\Gamma_{\rm
VSI,max}$ by a factor of more than 10 at all $R\gtrsim 5~{}\rm au$, and
increasing $a_{\rm max}$ to 1 mm completely stabilizes the VSI at $R>$ 10 au
(but see Section 4.3 for the dependence on the disk mass). The maximum growth
rate decreases mainly because the vertical shear $\partial(R\Omega)/\partial
z$ at the top of the VSI zone, where we evaluate $\Gamma_{\rm VSI,max}$,
decreases as the VSI zone shrinks vertically (see Figure 1(c)). Decreasing
$\alpha_{\rm D}$ also causes a decrease in $\Gamma_{\rm VSI,max}$, but this
effect is minor compared to the effect of varying $a_{\rm max}$. Furthermore,
decreasing $\alpha_{D}$ broadens the range of $k_{x}$ to the short wavelength
side (see Equation (23)) and promotes instability. However, this effect plays
a minor role in the variation of $\Gamma_{\rm VSI,max}$ with $\alpha_{D}$.

Since our dispersion relation assumes zero viscosity and infinitesimally short
cooling times, it is likely to underestimate the growth rates of modes with
$k_{x}\sim k_{x,\rm max}$ and $|k_{z}|\sim|k_{z}|_{\rm min}$, respectively.
Because the most unstable modes in our analysis have both $k_{x}=k_{x,\rm
max}$ and $|k_{z}|=|k_{z}|_{\rm min}$, the maximum growth rate of the VSI
predicted from our analysis should be regarded as an upper limit.

###  4.3 Variation of the VSI Zone Size with the Dust-to-Gas Ratio and Disk
Mass

Figure 10: Location of the VSI zone (shaded area) for different values of
$\Sigma_{d}/\Sigma_{g}$ and $M_{\rm disk}$ with $a_{{\rm max}}=1~{}{\rm mm}$
and $\alpha_{D}=10^{-4}$. The dashed lines represent $H_{g}$, $2H_{g}$, and
$3H_{g}$ in height from the midplane.

So far we have fixed the dust-to-gas mass ratio and disk mass to
$\Sigma_{d}/\Sigma_{g}=0.01$ and $M_{\rm disk}=0.01M_{\sun}$. The dust-to-gas
ratio can decrease with time as dust particles tend to drift toward the
central star owing to disk gas drag (Brauer et al., 2008). The disk mass can
also depend on disk age; in particular, very young disks like the one around
HL Tau can be as massive as $M_{\rm disk}\sim 0.1M_{\sun}$ (Kwon et al.,
2015). As highlighted in Section 4.1, the size of the VSI zone is controlled
by $\ell_{{\rm gd}}$, with larger $\ell_{{\rm gd}}$ leading to smaller VSI
zones. Because $\ell_{\rm gd}(z=0)\propto\Sigma_{\rm g}^{1/2}/\Sigma_{\rm
d}\propto(\Sigma_{\rm d}/\Sigma_{\rm g})^{-1}M_{\rm disk}^{-1/2}$ (see
Equation (38)), the VSI zone shrinks and expands as $\Sigma_{\rm
d}/\Sigma_{\rm g}$ decreases and $M_{\rm disk}$ increases, respectively.

We illustrate this in Figure 10, where we show the maps of the VSI zone for a
disk with a lower dust-to-gas ratio of $\Sigma_{d}/\Sigma_{g}=0.001$ and for a
disk with a higher disk mass of $M_{\rm disk}=0.1~{}M_{\sun}$ (panels (b) and
(c), respectively), both with $a_{{\rm max}}=1{\rm~{}mm}$ and
$\alpha_{D}=10^{-4}$. We find that the VSI zone is completely removed from
$R>5~{}\rm au$ for $\Sigma_{d}/\Sigma_{g}=0.001$, whereas it extends out to
100 au for $M_{\rm disk}=0.1M_{\odot}$ even with $a_{\rm max}=1~{}\rm mm$. The
radial extent of the VSI zone for the latter case is almost the same as that
for $M_{\rm disk}=0.01M_{\sun}$ and $a_{\rm max}=10~{}\micron$ (see Figure
5(a)), consistent with the scaling $\ell_{\rm gd}(z=0)\propto a_{\rm
max}^{1/2}(\Sigma_{\rm d}/\Sigma_{\rm g})^{-1}M_{\rm disk}^{-1/2}$ from
Equation (38).

##  5 Discussion

###  5.1 Correspondence between Vertically Local and Global Analyses

As mentioned in Section 3, we have employed the vertically local approximation
to treat the modes within a VSI zone. For vertically wide ($L\sim H_{g}$) VSI
zones, our local modes should correspond to some modes in the vertically
global linear analysis.

Below, we show that the modes with the longest vertical wavelength
($|k_{z}|=|k_{z}|_{\rm min}$, Equation (24)) indeed correspond to the
fundamental corrugation modes in the global analysis. The fundamental
corrugation modes are the modes with uniform vertical motion and are known to
dominate the nonlinear phase of the VSI (Nelson et al., 2013). In the limit of
short cooling times, the vertically global dispersion relation for the
fundamental corrugation modes can be written as (Lin & Youdin, 2015)

$\omega^{2}=\frac{1+ihq\hat{k}_{x}}{1+\hat{k}_{x}^{2}}\Omega_{\rm K}^{2},$
(39)

where $h=H_{g}/R$ and $\hat{k}_{x}=k_{x}H_{g}$. To show the correspondence
between these vertically uniform modes and our local VSI modes with
$|k_{z}|=|k_{z}|_{\rm min}$, we focus on local inertial modes and neglect the
$\omega^{4}$ in our local dispersion relation (Equation (2.2)).We also use
$\partial v_{y0}/\partial z\sim(qH_{g}/R)\Omega_{\rm K}$,
$N_{z}\ll\kappa_{0}\approx\Omega_{\rm K}$, $g\sim\Omega_{\rm K}^{2}H_{g}$, and
$k_{z}=-|k_{z}|_{\rm min}\sim-1/H_{g}$ for $L\sim H_{g}$. Our dispersion
relation then reduces to

$\omega^{2}\sim\frac{1-2\hat{k}_{x}|q|h+ihq\hat{k}_{x}}{2+\hat{k}_{x}^{2}}\Omega_{\rm
K}^{2},$ (40)

Since $\hat{k}_{x}|q|h\leq\hat{k}_{x,\rm max}|q|h\lesssim 1$ for
$\alpha_{D}\gtrsim 10^{-5}$ (see Equation (23)), Equation (40) agrees with
Equation (39) to within a factor of order unity. Because the most unstable
modes in our analysis have $|k_{z}|=|k_{z}|_{\rm min}$, our $\Gamma_{\rm
VSI,\rm max}$ serves as a good estimate for the growth rate of the fundamental
corrugation modes when $L\sim H_{g}$.

###  5.2 Can VSI-driven Turbulence Stably Sustain Vertical Dust Distribution?

The main limitation of our model is that it has to assume the strength of disk
turbulence to calculate the vertical dust distribution. In reality, in outer
disk regions where ambipolar diffusion suppresses MRI, the VSI itself can be
the dominant source of disk turbulence. If this is the case, our model
effectively assumes that the strength of vertical diffusion caused by VSI-
driven turbulence matches the diffusion strength $\alpha_{D}$ required to
sustain the vertical dust distribution. The system would evolve until two
diffusion coefficients match, but we cannot tell if such equilibrium states
would exist because our current model relying on linear stability analysis
does not predict the strength of VSI-driven turbulence.

Moreover, even if there exists an equilibrium state, the state can be unstable
against perturbations to the vertical dust distribution. Unstable equilibrium
is expected if, for instance, a small decrease in the the dust scale height
causes a large decrease in the VSI turbulence strength, in which case dust
settling would proceed in a runaway fashion . Dust settling also introduces
effective vertical buoyancy that further stabilizes the VSI in the settled
dust layers (Lin, 2019; Schäfer et al., 2020). The stability of the system
should be studied in future hydrodynamical simulations that include both the
thermal and frictional coupling between the gas and dust.

###  5.3 Implications for Dust Growth, Settling, and Planetesimal Formation in
Outer Disk Regions

We have shown that dust growth can substantially suppresses the VSI beyond 10
au. This suggests that disk turbulence in this outer disk region would become
weaker as dust grows. This may provide positive feedback to dust growth and
also to planetesimal formation. Weaker turbulence would suppress collisional
fragmentation of the dust particles and thereby further promote dust growth
(e.g., Brauer et al., 2008; Birnstiel et al., 2010; Okuzumi & Hirose, 2012).
Weak turbulence are also preferred for planetesimal formation via the
streaming and gravitational instabilities, both of which require substantial
dust settling toward the midplane (Sekiya, 1998; Youdin & Shu, 2002; Johansen
et al., 2009). These positive feedback effects are potentially important for
understanding planet formation and dust ring/gap formation in outer regions of
protoplanetary disks.

Suppression of the VSI at large radial distances due to dust growth may also
explain the high degree of dust settling in the HL Tau disk (Pinte et al.,
2016). However, quantitative estimates for the strength of VSI-driven
turbulence are needed to test this hypothesis because, as we have seen in
Section 4.3, the VSI tends to be vigorous in massive disks like the HL Tau
disk ($M_{\rm disk}\sim 0.1M_{\sun}$; Kwon et al. 2015). Very recently, Doi &
Kataoka (2021) measured the vertical thicknesses of two dust rings at 70 and
100 au in the massive disk around HD 163296, and showed that the outer ring is
much thinner than the gas disk but the inner ring is as thick as the gas disk.
This may suggest that the VSI is active at $\sim 70~{}\rm au$ in this disk
(see Bi et al. 2021 and Binkert et al. 2021 for another potential
interpretation). We plan to address these open issues in future work.

##  6 Conclusions

We have investigated the impacts of dust growth and settling on the VSI in the
outer regions of protoplanetary disks using a model based on linear stability
analysis. Our key findings are summarized as follows.

  1. 1.

For fixed dust particle size distribution, a higher degree of dust settling
(corresponding to a lower value of turbulence strength $\alpha_{D}$) leads to
a VSI zone that is more confined to the midplane and more extended to larger
radial distances (Figure 3). This is because in outer disk regions of low
optical depths, dust settling causes a decrease and an increase in the
timescale of thermal conduction from gas to dust ($\propto\ell_{\rm gd}$;
Equation (19)) at the midplane and well above the midplane, respectively
(Figure 4).

  2. 2.

For fixed turbulence strength, the VSI zone shrinks toward the midplane and
also toward the central star as dust particles grow (Figure 5). Dust growth
also causes dust settling, but the decrease in the total surface area of the
particles due to dust growth dominates the change in the thermal relaxation
timescale (Figure 6).

  3. 3.

The maximum growth rate of the VSI decreases as dust particles grow. In our
default disk model assuming $M_{\rm disk}=0.01M_{\sun}$, dust growth to 1 mm
in size completely stabilizes the VSI exterior to 10 au (Figure 9). On the
other hand, the VSI is more vigorous in disks with larger disk masses (Section
4.3).

Our results suggest that dust evolution, in particular dust growth, should
lead to suppression of VSI-driven turbulence. This effect may enable further
dust coagulation and settling and may potentially promote planetesimal
formation in outer disk regions. The effect may also explain the high degree
of dust settling observed in the dust rings around HL Tau, but testing this
hypothesis requires more qualitative investigation of the interplay between
dust evolution and the nonlinear development of the VSI.

We thank Shoji Mori for the discussions that motivated this project. We also
thank the referee for a helpful report that motivated us to discuss the
correspondence between local and global linear analyses. This work was
supported by JSPS KAKENHI Grant Numbers JP20H01948, JP20H00182, JP19K03926,
JP18H05438 and JP20J01376.

## References

  * Adachi et al. (1976) Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756, doi: 10.1143/PTP.56.1756
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
  * Andrews (2020) Andrews, S. M. 2020, ARA&A, 58, 483, doi: 10.1146/annurev-astro-031220-010302
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Arlt & Urpin (2004) Arlt, R., & Urpin, V. 2004, A&A, 426, 755, doi: 10.1051/0004-6361:20035896
  * Bai (2015) Bai, X.-N. 2015, ApJ, 798, 84, doi: 10.1088/0004-637X/798/2/84
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Barker & Latter (2015) Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21, doi: 10.1093/mnras/stv640
  * Bi et al. (2021) Bi, J., Lin, M.-K., & Dong, R. 2021, arXiv e-prints, arXiv:2103.09254.  https://arxiv.org/abs/2103.09254
  * Binkert et al. (2021) Binkert, F., Szulágyi, J., & Birnstiel, T. 2021, arXiv e-prints, arXiv:2103.10177.  https://arxiv.org/abs/2103.10177
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79, doi: 10.1051/0004-6361/200913731
  * Birnstiel et al. (2011) Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11, doi: 10.1051/0004-6361/201015228
  * Bitsch et al. (2015) Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28, doi: 10.1051/0004-6361/201424964
  * Bohren & Huffman (1983) Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles
  * Brauer et al. (2008) Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859, doi: 10.1051/0004-6361:20077759
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43, doi: 10.1051/0004-6361/201425120
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368, doi: 10.1086/304869
  * Chokshi et al. (1993) Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806, doi: 10.1086/172562
  * Cui & Bai (2020) Cui, C., & Bai, X.-N. 2020, ApJ, 891, 30, doi: 10.3847/1538-4357/ab7194
  * Doi & Kataoka (2021) Doi, K., & Kataoka, A. 2021, arXiv e-prints, arXiv:2102.06209.  https://arxiv.org/abs/2102.06209
  * Dominik & Tielens (1997) Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647, doi: 10.1086/303996
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237, doi: 10.1006/icar.1995.1058
  * Flaherty et al. (2020) Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109, doi: 10.3847/1538-4357/ab8cc5
  * Flaherty et al. (2015) Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99, doi: 10.1088/0004-637X/813/2/99
  * Flaherty et al. (2018) Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117, doi: 10.3847/1538-4357/aab615
  * Flaherty et al. (2017) Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150, doi: 10.3847/1538-4357/aa79f9
  * Flock et al. (2017) Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131, doi: 10.3847/1538-4357/aa943f
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155, doi: 10.3847/1538-4357/ab9641
  * Fricke (1968) Fricke, K. 1968, ZAp, 68, 317
  * Goldreich & Schubert (1967) Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571, doi: 10.1086/149360
  * Goldreich & Ward (1973) Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051, doi: 10.1086/152291
  * Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385, doi: 10.1086/305277
  * Johansen & Youdin (2007) Johansen, A., & Youdin, A. 2007, ApJ, 662, 627, doi: 10.1086/516730
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75, doi: 10.1088/0004-637X/704/2/L75
  * Kataoka et al. (2013) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4, doi: 10.1051/0004-6361/201322151
  * Kwon et al. (2015) Kwon, W., Looney, L. W., Mundy, L. G., & Welch, W. J. 2015, ApJ, 808, 102, doi: 10.1088/0004-637X/808/1/102
  * Latter & Papaloizou (2018) Latter, H. N., & Papaloizou, J. 2018, MNRAS, 474, 3110, doi: 10.1093/mnras/stx3031
  * Lin (2019) Lin, M.-K. 2019, MNRAS, 485, 5221, doi: 10.1093/mnras/stz701
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17, doi: 10.1088/0004-637X/811/1/17
  * Lin & Youdin (2017) —. 2017, ApJ, 849, 129, doi: 10.3847/1538-4357/aa92cd
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * Lyra & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001, doi: 10.1088/1538-3873/aaf5ff
  * Malygin et al. (2017) Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dullemond, C. P. 2017, A&A, 605, A30, doi: 10.1051/0004-6361/201629933
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610, doi: 10.1093/mnras/stt1475
  * Nelson et al. (2016) —. 2016, MNRAS, 456, 239, doi: 10.1093/mnras/stv2440
  * Okuzumi & Hirose (2012) Okuzumi, S., & Hirose, S. 2012, ApJ, 753, L8, doi: 10.1088/2041-8205/753/1/L8
  * Okuzumi et al. (2012) Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106, doi: 10.1088/0004-637X/752/2/106
  * Pfeil & Klahr (2019) Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150, doi: 10.3847/1538-4357/aaf962
  * Pierens (2021) Pierens, A. 2021, MNRAS, doi: 10.1093/mnras/stab183
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25, doi: 10.3847/0004-637X/816/1/25
  * Riols & Lesur (2018) Riols, A., & Lesur, G. 2018, A&A, 617, A117, doi: 10.1051/0004-6361/201833212
  * Schäfer et al. (2020) Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190, doi: 10.1051/0004-6361/201937371
  * Sekiya (1998) Sekiya, M. 1998, Icarus, 133, 298, doi: 10.1006/icar.1998.5933
  * Simon et al. (2013a) Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013a, ApJ, 775, 73, doi: 10.1088/0004-637X/775/1/73
  * Simon et al. (2013b) Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013b, ApJ, 764, 66, doi: 10.1088/0004-637X/764/1/66
  * Stoll & Kley (2014) Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77, doi: 10.1051/0004-6361/201424114
  * Takahashi & Inutsuka (2014) Takahashi, S. Z., & Inutsuka, S.-i. 2014, ApJ, 794, 55, doi: 10.1088/0004-637X/794/1/55
  * Takeuchi & Lin (2002) Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344, doi: 10.1086/344437
  * Tominaga et al. (2018) Tominaga, R. T., Inutsuka, S.-i., & Takahashi, S. Z. 2018, PASJ, 70, 3, doi: 10.1093/pasj/psx143
  * Tominaga et al. (2019) Tominaga, R. T., Takahashi, S. Z., & Inutsuka, S.-i. 2019, ApJ, 881, 53, doi: 10.3847/1538-4357/ab25ea
  * Tominaga et al. (2020) —. 2020, ApJ, 900, 182, doi: 10.3847/1538-4357/abad36
  * Urpin (2003) Urpin, V. 2003, A&A, 404, 397, doi: 10.1051/0004-6361:20030513
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399, doi: 10.1046/j.1365-8711.1998.01118.x
  * van der Marel et al. (2019) van der Marel, N., Dong, R., di Francesco, J., Williams, J. P., & Tobin, J. 2019, ApJ, 872, 112, doi: 10.3847/1538-4357/aafd31
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57, doi: 10.1093/mnras/180.2.57
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Windmark et al. (2012) Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012, A&A, 544, L16, doi: 10.1051/0004-6361/201220004
  * Yang et al. (2017) Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80, doi: 10.1051/0004-6361/201630106
  * Youdin (2011) Youdin, A. N. 2011, ApJ, 731, 99, doi: 10.1088/0004-637X/731/2/99
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459, doi: 10.1086/426895
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588, doi: 10.1016/j.icarus.2007.07.012
  * Youdin & Shu (2002) Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494, doi: 10.1086/343109

##  Appendix A Thermal Relaxation Timescale

Figure 11:  Thermal relaxation timescale $\tau_{{\rm relax}}$ normalized by
$\Omega_{\rm K}^{-1}$ as a function of $R$ and $z/R$ for different values of
$\alpha_{D}$ with $a_{\rm max}=100{\rm~{}\micron}$. The dashed lines represent
$z=H_{g}$, $2H_{g}$, and $3H_{g}$. The solid line marks $\tau_{{\rm
relax}}=\tau_{{\rm crit}}$.  Figure 12:  Thermal relaxation timescale
$\tau_{{\rm relax}}$ normalized by $\Omega_{\rm K}^{-1}$ as a function of $R$
and $z/R$ for different values of $a_{\rm max}$ with $\alpha_{D}=10^{-4}$. The
dashed lines represent $z=H_{g}$, $2H_{g}$, and $3H_{g}$. The solid line marks
$\tau_{{\rm relax}}=\tau_{{\rm crit}}$.

Figures 11 and Figure 12 show $\tau_{\rm relax}$ ($=\tau_{\rm coll}$)
normalized by $\Omega_{\rm K}^{-1}$ as a function of $R$ and $z$ for all
parameter sets considered in this study. A smaller $\alpha_{D}$ produces a
larger $\tau_{\rm relax}$ well above the midplane (because $\tau_{\rm
relax}=\tau_{\rm coll}\propto\ell_{\rm gd}$; see Figure 4), resulting in a VSI
zone that is more confined to the midplane (see the boundaries of the VSI
zones marked by the solid lines in Figure 11). A larger $a_{\rm max}$ causes
an increase in $\tau_{\rm relax}$ at all heights and hence leads to a narrower
VSI zone.