# Hydro-, Magnetohydro-, and Dust-Gas Dynamics of Protoplanetary Disks

G. Lesur1, B. Ercolano2, M. Flock3, M.-K. Lin4,5, C.-C. Yang6,
J. A. Barranco7, P. Benitez-Llambay8, J. Goodman9, A. Johansen10,11, H.
Klahr3, G. Laibe12,13, W. Lyra14, P. Marcus15, R.P. Nelson16, J. Squire17, J.
B. Simon18, N. Turner19, O.M. Umurhan20,21, A.N. Youdin22 1 Univ. Grenoble
Alpes, CNRS, IPAG, 38000 Grenoble, France
2 Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität
München, Scheinerstr. 1, 81679 München, Germany
3 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117, Heidelberg,
Germany
4 Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 10617,
Taiwan
5 Physics Division, National Center for Theoretical Sciences, Taipei 10617,
Taiwan
6 Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505
S. Maryland Parkway, Box 454002, Las Vegas, NV 89154-4002, USA
7 Dept. of Physics & Astronomy, San Francisco State University, San Francisco,
CA 94132, U.S.A.
8 Niels Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17,
DK-2100 Copenhagen Ø, Denmark
9 Department of Astrophysical Sciences, Princeton University, Princeton NJ,
U.S.A.
10 Center for Star and Planet Formation, GLOBE Institute, University of
Copenhagen, Øster Voldgade 5-7, 1350 Copenhagen, Denmark
11 Lund Observatory, Department of Astronomy and Theoretical Physics, Lund
University, Box 43, 221 00 Lund, Sweden
12 Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230, Saint-Genis,-Laval, France.
13 Institut Universitaire de France
14 Department of Astronomy, New Mexico State University, PO Box 30001, MSC
4500 Las Cruces, NM 88003, USA
15 Dept. of Mechanical Engineering, University of California, Berkeley, CA,
94720, U.S.A.
16 Department of Physics & Astronomy, Queen Mary University of London, Mile
End Road, London, E1 4NS, U.K.
17 Physics Department, University of Otago, Dunedin 9010, New Zealand
18 Department of Physics and Astronomy, Iowa State University, Ames, IA,
50010, USA
19 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Drive, Pasadena, California 91109, U.S.A.
20 SETI Institute, 339 Bernardo Ave., Suite 200, Mountain View, CA 94043,
U.S.A.
21 Space Sciences Division, NASA Ames Research Center, Moffett Field, CA
94035, U.S.A.
22 University of Arizona, Departments of Astronomy and Planetary Sciences

###### Abstract

The building of planetary systems is controlled by the gas and dust dynamics
of protoplanetary disks. While the gas is simultaneously accreted onto the
central star and dissipated away by winds, dust grains aggregate and collapse
to form planetesimals and eventually planets. This dust and gas dynamics
involves instabilities, turbulence and complex non-linear interactions which
ultimately control the observational appearance and the secular evolution of
these disks.

This chapter is dedicated to the most recent developments in our understanding
of the dynamics of gaseous and dusty disks, covering hydrodynamic and
magnetohydrodynamic turbulence, gas-dust instabilities, dust clumping and disk
winds. We show how these physical processes have been tested from observations
and highlight standing questions that should be addressed in the future.

##  1 INTRODUCTION

We focus in this chapter on the progress since PPVI in the fields of angular
momentum transport, turbulence, dust growth, and winds from planet-forming
disks. We begin by reviewing the basics of disk dynamics. The interested
reader may consult the reviews and lecture notes provided in reference for
more detailed discussion of these fundamental aspects.

###  1.1 Framework and observational constraints

The ongoing improvements in observational techniques, experiments, and
simulations are opening a new window on our understanding of protoplanetary
disks dynamics, as it becomes possible to test and validate scenarios that
were, until recently, only theoretical concepts. In parallel, dynamical models
are being refined, treating processes such as radiative heating and cooling,
non-ideal magnetohydrodynamics (MHD), and outflows, to name a few. These
efforts have led to significant revisions in our understanding of the
mechanisms responsible for mass accretion. In addition, it has been realized
that the dust/gas coupling is much richer than anticipated, leading to
recognition of new multi-phase instabilities and dynamical phenomena that can
enhance planet-formation rates. The goal of this chapter is to review these
recent advances and to point out the open questions and some links between
them.

In the following and unless stated otherwise, we consider an isolated
protoplanetary disk, perturbed neither by a binary stellar companion nor by
infalling material from a remnant envelope with mass $M_{\mathrm{disk}}$
orbiting a central young stellar object of mass $M_{\mathrm{star}}$ assumed to
be near a solar mass. This restriction allows us to make general statements
about the disk independent of environmental effects, which should be
applicable to most class II disks. We note however that disks are not always
isolated, as streamers are sometimes observed extending over thousands of
astronomical units (e.g. _Pineda et al._ , 2020). We assume the disk is low
enough in mass to be non-self-gravitating, i.e.
$M_{\mathrm{disk}}/M_{\mathrm{star}}<O(10^{-1})$ (_Armitage_ , 2011). More
precise threshold disk masses can be found in, e.g. _Kratter and Lodato_
(2016); _Haworth et al._ (2020). The disk is accreting on the star, forming
planets from dust grains, and ultimately dispersing on a timescale of a few
million years. In this introduction, we define the main quantities and
concepts used throughout the chapter relating to accretion, angular momentum
transport (1.2), and dust dynamics (1.3). In subsequent sections, we focus on
instabilities of the gas phase leading to hydrodynamical turbulence (§2). We
emphasize the conditions for existence of these instabilities which turn out
to depend mainly on thermodynamics. We next move to the coupling with dust
gains (§3), describing instabilities and dust clumping arising because of
dust-gas interactions. In §4, we discuss the roles of magnetic fields,
magneto-rotational turbulence, and winds launched by gas pressure and magnetic
forces. Finally, we link these theoretical processes to present and future
observations in §5. We conclude with a summary of the major achievements of
the field since PPVI and highlight several key questions which should be
addressed in the future.

###  1.2 Accretion, angular momentum transport, and turbulence

Protoplanetary disks are accretion disks, delivering material onto their
central stars at rates typically between $10^{-10}$ and
$10^{-7}\,M_{\odot}/\mathrm{yr}$ (_Venuti et al._ , 2014; _Hartmann et al._ ,
2016; _Manara et al._ , 2016). Their formation, evolution, and dispersal are
closely linked to their angular momentum content and, since angular momentum
is a conserved quantity, how the angular momentum is transported.

This question of angular momentum transport has been the subject of a great
deal of work over the past 50 years not only in protoplanetary disks, but also
in the disks around stellar remnants and in active galactic nuclei.
Essentially two transport mechanisms have been proposed: turbulence and
magnetized disk winds (MDW). The former has attracted most attention because
of the simplicity of treating it using the _Shakura and Sunyaev_ (1973)
pseudo-viscosity paradigm. Winds, initially proposed by _Blandford and Payne_
(1982), have regained popularity over the last ten years as a way to transport
angular momentum in otherwise ”dead” discs.

The links between accretion flow, angular momentum transport, and disks’
secular evolution can be seen in the equations of mass and angular momentum
conservation, averaged through the disk thickness (e.g. _Lesur_ 2021a):

$\displaystyle\frac{\partial\Sigma}{\partial t}+\frac{1}{2\pi
R}\frac{\partial}{\partial R}\dot{M}_{\mathrm{acc}}$
$\displaystyle=-\zeta\Sigma\Omega_{\mathrm{K}},$ (1)
$\displaystyle\frac{\dot{M}_{\mathrm{acc}}\Omega_{\mathrm{K}}}{4\pi}$
$\displaystyle=\frac{1}{R}\frac{\partial}{\partial
R}\left(R^{2}\alpha_{\mathrm{S}}P\right)$
$\displaystyle\qquad+\zeta(\lambda-1)P\left(\frac{R}{H}\right)^{2},$ (2)

where we define the disk mass accretion rate $\dot{M}_{\mathrm{acc}}$, surface
density $\Sigma$ and pressure $P$, Keplerian angular velocity
$\Omega_{\mathrm{K}}$ and geometrical thickness $H$. In these equations we
include three unknown dimensionless coefficients: $\zeta$, the mass loss
parameter due to a hypothetical outflow, $\alpha_{\mathrm{S}}$ which
corresponds to _Shakura and Sunyaev_ (1973) definition and describes radial
angular momentum transport within the disk, such as by turbulence, and finally
$\lambda$, the outflow lever arm (_Blandford and Payne_ , 1982), which
quantifies the specific angular momentum extracted vertically by an outflow.
All three of these coefficients may vary with radius and time. Once the three
coefficients are known, the secular evolution of the system can be entirely
predicted and accretion theory is complete.

If the disk accretion is driven solely by turbulence of hydrodynamical or
magnetohydrodynamical origin, then the accretion rate can be derived from (2)
and is approximately (_Hartmann et al._ , 1998)

$\displaystyle\dot{M}_{\mathrm{acc}}\sim\left(3\times
10^{-8}\,M_{\odot}/\mathrm{yr}\right)\left(\frac{\alpha_{\mathrm{S}}}{10^{-2}}\right)\left(\frac{\varepsilon}{0.1}\right)^{2}$
$\displaystyle\times\left(\frac{R}{10\,\mathrm{AU}}\right)^{1/2}\left(\frac{M_{\mathrm{star}}}{M_{\odot}}\right)^{1/2}\left(\frac{\Sigma}{10\,\mathrm{g.cm}^{-2}}\right)$
(3)

where we use the disk aspect ratio $\varepsilon\equiv H/R$ and assume
$\Sigma(R)$ follows a shallow power law (typically $\propto R^{-1/2}$). Hence
typically $\alpha_{\mathrm{S}}\sim 10^{-2}$ is required to match the observed
accretion rates by turbulence alone.

It is worth noting that gas turbulence need not always lead to outward angular
momentum transport. One can obtain $\alpha_{\mathrm{S}}\sim 0$ from turbulence
that is vigorous in the sense of large velocity fluctuations. This is
particularly important in the context of dust transport and settling: some
instabilities lead to efficient stirring of dust grains yet yield little
angular momentum transport (e.g. 2.1). Hence the $\alpha_{\mathrm{S}}$
parameter gives incomplete information on the presence and impact of
turbulence.

###  1.3 Dust transport & growth

Dust grains have porous and fractal structures (_Dominik et al._ , 2007; _Blum
and Wurm_ , 2008). They are often modeled as spheres of radii $s\lesssim 10$
cm and densities up to $\rho\sim 1$ g.cm-3 for compacted material (_Love et
al._ , 1994). Aggregates evolve as they undergo collisions. Sticking is
favored for small ice-coated grains, while collisions between large aggregates
involve redistribution of energy via elastic waves, which can trigger sliding
or breaking at points of contact between substructures. Sticking, bouncing,
compaction, abrasion or fragmentation, mass transfer, cratering, and erosion
are all possible outcomes of collisions (_Blum_ , 2018), as we discuss in §
3.1.

The drag stopping time is $t_{\rm s}\sim(\rho_{\rm m}s)/(\rho_{\rm g}c_{\rm
s})$, where $\rho_{\rm m}$ is the material density of the grain and $\rho_{\rm
g}$ is the density of the surrounding gas (_Epstein_ 1923; ; _Baines et al._
1965; _Clair et al._ 1970). Drag damps eccentricity and makes them settle and
drift radially (_Adachi et al._ , 1976). The competition between drag and the
stellar gravity is (_Safronov_ , 1969; _Whipple_ , 1972, 1973)

$\mathrm{St}\equiv\tau_{\mathrm{s}}\equiv\Omega_{\mathrm{K}}t_{\rm s}.$ (4)

for $\mathrm{St}\sim 1$. Solids concentrate strongly as they from the gas,
making for planet formation. In a vertically isothermal disk, Eq. (4)
$\mathrm{St}=((\rho_{\mathrm{m}}s)/(\sqrt{2\pi}\Sigma_{\mathrm{g}}))\exp(z^{2}/2H^{2})$,
that the surface density $\Sigma_{\rm g}$ is key for setting values of
$\mathrm{St}$ and thus, dust dynamics all through the disk. In the midplane of
a disk $\Sigma_{\rm g}\propto R^{-1}$ and $\Sigma_{\rm g}\sim 200$ g.cm-2 at
$R=1$ AU (which match the low end of the observed mass distribution, e.g.
_Andrews et al._ 2009) millimeter-sized grains have $\mathrm{St}_{0}\sim 0.04$
at $R\sim 50$ AU. , most grains have $\mathrm{St}\gtrsim 1$. Since
$\Sigma_{g}$ can not be inferred directly from observations, values of
$\mathrm{St}$ major uncertainties .

In the vertical direction, and in the absence of turbulence, grains damped
harmonic oscillators of quality factors
$\sim\mathrm{St}_{0}=\mathrm{St}\left(z=0\right)$ driven by the mean flow of
the gas at large scale. : the diffusivity $\alpha_{\mathrm{D}}$ (which is not
necessarily equal to the angular momentum transport coefficient
$\alpha_{\mathrm{S}}$; see Section 5.2) and the correlation time at eddy
scales $t_{\rm e}$ (_Youdin and Lithwick_ , 2007). A classic choice is
$\alpha_{\mathrm{D}}\sim\alpha_{\mathrm{S}}$ and $t_{\rm
e}\sim\Omega_{\mathrm{K}}^{-1}$, the typical time for vortex stretching by
differential rotation (see however § 2 and § 3.2.2). In steady state, solids
concentrate close to the midplane within a layer of height

$H_{\rm d}\simeq
H\sqrt{\frac{\alpha_{\mathrm{D}}}{\alpha_{\mathrm{D}}+\mathrm{St}}}$ (5)

(see also _Yang et al._ , 2018). As a result, larger grains .

In the radial direction, between gas and dust. for both phases. In an
inviscid, non-magnetic, planet-free Keplerian disk, , given by

$v_{\mathrm{d},R}=\frac{1}{\mathrm{St}+\mathrm{St}^{-1}\left(1+\epsilon_{\mathrm{d}}\right)^{2}}\frac{1}{\Omega_{\mathrm{K}}\rho_{\mathrm{g}}}\frac{\partial
P}{\partial R}=-\epsilon_{\mathrm{d}}^{-1}v_{\mathrm{g},R},$ (6)

where $\epsilon_{\mathrm{d}}=\rho_{\rm d}/\rho_{g}$ is the local dust-to-gas
density ratio (_Nakagawa et al._ , 1986). drifts towards pressure maxima –
either the inner or local pressure bumps – with optimal efficiency when
$\mathrm{St}\sim 1$. grains travel through the disk in a typical time $t_{\rm
drift}\sim\Omega_{\mathrm{K}}^{-1}(H/R)^{-2}(1+\mathrm{St}^{2})/\mathrm{St}$.

##  2 THERMO-HYDRODYNAMICAL INSTABILITIES

Table 1: Candidates for hydrodynamic activity in PPDs: the Vertical Shear
Instability (VSI), the Convective Overstability (COS), and the Zombie Vortex
Instability (ZVI). Figures adapted from _Pfeil and Klahr_ (2021), _Lyra_
(2014), and _Barranco et al._ (2018), respectively. The vertical component of
the vorticity is shown for each, while the Reynolds stress is also shown for
the VSI. Each instability requires a different gas cooling timescale
$\tau_{\mathrm{cool}}$, here scaled by the Keplerian orbital frequency
$\Omega_{\mathrm{K}}$. They also require different disk structures, here
expressed for a vertically isothermal disk with a radial temperature profile
$\propto R^{-q}$, midplane density profile $\propto R^{-p}$, adiabatic index
$\gamma$, the height above the midplane $z$, and the gas pressure scale height
$H$. $R$ is the cylindrical distance from the star. A Shakura-Sunyaev
viscosity, $\alpha_{\mathrm{SS}}$, can be measured from simulations as a
metric for radial angular momentum transport, but does not fully characterize
the turbulence generated by these instabilities.

Vertical Shear Instability | Convective Overstability | Zombie Vortex Instability
---|---|---
|  |
$\tau_{\mathrm{cool}}\Omega_{\mathrm{K}}\ll 1$ | $\tau_{\mathrm{cool}}\Omega_{\mathrm{K}}\sim 1$ | $\tau_{\mathrm{cool}}\Omega_{\mathrm{K}}\gg 1$
$q\neq 0$ | $-1<p/q<1/(\gamma-1)$ | $|z|\gtrsim\sqrt{\gamma/(\gamma-1)}H$
$\alpha_{\mathrm{SS}}\sim 10^{-4}$ | $\alpha_{\mathrm{SS}}\sim 10^{-3}$ |  $\alpha_{\mathrm{SS}}\sim(10^{-5}$–$10^{-4})$$\dagger$
Outcome: turbulence & vortices

  * $\dagger$

Based on incompressible simulations (_Barranco and Marcus_ , 2005).

PPDs are poorly ionized and weakly coupled to magnetic fields (see §4 for
details). The lack of magnetically-driven turbulence in some regions has led
to a surge in the interest in purely hydrodynamic origins of turbulence. To
this end, three new hydrodynamic instabilities have been discovered:  1) the
Vertical Shear Instability (VSI);  2) the Convective Overstability (COS); and
3) the Zombie Vortex Instability (ZVI).  Key features of these instabilities
are summarized in Table 1. Which of these operate depends on the disk’s
structure and the gas thermal timescales. However, taken together they likely
pertain to large swaths of PPDs (see §2.4) and are expected to produce weak to
moderate turbulence, transport, and large-scale structures such as vortices,
which are all relevant to interpreting disk observations as well as
planetesimal formation (§3).

To anchor the following discussion, we consider a three-dimensional (3D) PPD
orbiting a central star of mass $M_{*}$. Cylindrical $(R,\phi,z)$ coordinates
are centered on the star with the disk plane located at $z=0$. The equilibrium
disk is steady and axisymmetric with a purely azimuthal flow
$\bm{v}_{\mathrm{g}}=R\Omega(R,z)\hat{\bm{\phi}}$, where $\Omega$ is the gas’
angular velocity, which is set by centrifugal balance and is close to
Keplerian rotation but not exactly. Note that we do not consider warped or
eccentric disks. The gas density $\rho_{\mathrm{g}}$ is set by vertical
hydrostatic balance and we take a power-law midplane density profile
$\rho_{\mathrm{g0}}(R)\propto R^{-p}$. For an ideal gas its temperature $T$ is
defined through $P=\mathcal{R}\rho T/\mu$, where $\mathcal{R}$ is the gas
constant and $\mu$ is the mean molecular weight. The gas’ specific entropy is
$S\equiv c_{P}\ln\left(P^{1/\gamma}/\rho_{\mathrm{g}}\right)$, where $c_{P}$
is the specific heat capacity at constant pressure and $\gamma$ is the
constant adiabatic index.

We consider vertically isothermal disks with power-law equilibrium temperature
profiles $T\propto R^{-q}$, which are often used to model the disk bulk or
interior in the outer regions of irradiated PPDs (_Hubeny_ , 1990; _Chiang and
Goldreich_ , 1997). The isothermal sound-speed is
$c_{s}\equiv\sqrt{\mathcal{R}T/\mu}$ and $H=c_{s}/\Omega_{\mathrm{K}}$. Then
the gas surface density $\Sigma_{\mathrm{g}}\propto R^{-\sigma}$ with
$\sigma=p+q/2-3/2$. PPDs are thin with aspect-ratios $\varepsilon\equiv H/R\ll
1$. Finally, we assume gravitationally stable disks with Toomre parameters
$Q_{g}=c_{s}\Omega/\pi G\Sigma_{g}\gg 1$, see _Kratter and Lodato_ (2016) for
a discussion of gravitational instabilities and the Chapters by Bae et al. and
Pinte et al. for their potential to explain spiral structures observed in
PPDs.

Prior to PPVI, PPDs were mostly considered hydrodynamically inactive, at least
to small-scale, local, infinitesimal perturbations. This notion stems from the
classic Solberg-Høiland (SH) stability criteria (e.g. _Tassoul_ , 1978), which
can be written as $\kappa^{2}+N_{R}^{2}+N_{z}^{2}>0$ and
$-\partial_{z}P\left(\kappa^{2}\partial_{z}S-R\partial_{z}\Omega^{2}\partial_{R}S\right)>0$.
Here, $\kappa^{2}=R^{-3}\partial_{R}\left(R^{4}\Omega^{2}\right)$ is the
square of the radial epicycle frequency; and
$N_{R,z}^{2}=-\left(c_{P}\rho_{\mathrm{g}}\right)^{-1}\partial_{R,z}P\partial_{R,z}S$
are the squares of the radial and vertical contributions to the buoyancy (or
Brunt-Väisälä) frequency, respectively. For radially smooth, vertically thin
disks, the SH stability criteria are generally satisfied (_Lin and Youdin_ ,
2015) because usually $N_{z}^{2}\geq 0$ (so the disk is stable against
vertical convection as well) and $|N_{R}^{2}|\ll\kappa^{2}$, while the second
criterion reduces to $\gamma>1+O(\varepsilon^{2})$, which is satisfied for
typical values of $\gamma=5/3$ or $7/5$.

However, the SH stability criteria assume the gas is adiabatic, i.e. it does
not experience heat gains or losses, and disturbances are axisymmetric and
infinitesimal. The key to the new class of thermo-hydrodynamic instabilities
in Table 1 is the violation of these assumptions (see discussion in _Fromang
and Lesur_ 2019). Both the VSI and COS require finite thermal losses. The ZVI
applies to adiabatic gas, but requires finite amplitude, non-axisymmetric
perturbations. These conditions can be expected in PPDs and thus the SH
criteria are not applicable. Furthermore, parts of PPDs may have a non-
standard radial structure such that $N_{R}^{2}<0$, which is necessary for the
COS.

The revived interest in hydrodynamic instabilities in PPDs thus resulted from
more careful treatments of the disk’s thermal structure and evolution. To this
end, we write the gas energy equation as
$\partial_{t}S+\bm{v}_{\mathrm{g}}\cdot\nabla S=-\Lambda.$ Here, $\Lambda$
represents all heating and cooling processes, but for simplicity we will refer
to it as ‘cooling’. Physically, cooling is mediated by radiation and its form
depends on the gas and dust properties (e.g. _Malygin et al._ , 2017;
_Barranco et al._ , 2018; _Pfeil and Klahr_ , 2019, see also §2.4). However,
for discussion purposes it is sufficient to consider the Newtonian cooling
prescription,

$\displaystyle\Lambda=\frac{c_{P}}{\gamma
T}\frac{\left(T-T_{\mathrm{ref}}\right)}{\tau_{\mathrm{cool}}},$ (7)

where the reference temperature $T_{\mathrm{ref}}$ is usually taken to be the
initial temperature field; and $\tau_{\mathrm{cool}}$ is the cooling timescale
over which the local temperature relaxes back to $T_{\mathrm{ref}}$. We define
the dimensionless cooling time
$\beta_{\mathrm{cool}}\equiv\tau_{\mathrm{cool}}\Omega_{\mathrm{K}}$.

The beta cooling prescription encapsulates all possible gas thermodynamic
responses. This enables the three hydrodynamic instabilities to be considered
by varying $\beta_{\mathrm{cool}}$. The limit $\beta_{\mathrm{cool}}\to\infty$
corresponds to adiabatic evolution, where the specific entropy is materially
conserved, which enables the ZVI provided $|N_{z}|$ is sufficiently large. In
the isothermal limit, $\beta_{\mathrm{cool}}\to 0$, or instant cooling, $T\to
T_{\mathrm{ref}}$ and the gas temperature does not evolve, which enables the
VSI provided there is a radial temperature gradient to produce vertical shear
(see below). In the intermediate case, $\beta_{\mathrm{cool}}\sim 1$, the disk
cools on the orbital timescale, which enables the COS if $N_{R}^{2}<0$ also.

The onset of these hydrodynamic instabilities thus has both structural and
thermodynamical requirements (Table 1). If these are met, they can lead to
hydrodynamic turbulence and vortex formation. In simulations, it is
straightforward to measure the conventional Shakura-Sunyaev
$\alpha_{\mathrm{S}}$ to describe radial angular momentum transport mediated
by said turbulence. However, it is important to note that
$\alpha_{\mathrm{S}}$ alone is insufficient to fully characterize the ensuing
turbulence. The fundamental measure of turbulence is the root-mean-squared
velocities or Mach numbers $\mathcal{M}$ at various scales. For planetesimal
formation, one is often interested in the resulting turbulent diffusion of
particles, $\alpha_{\mathrm{D}}$ (see §1.3 and §3). One often assumes
$\mathcal{M}^{2}$, $\alpha_{\mathrm{S}}$, and $\alpha_{\mathrm{D}}$ are equal
and isotropic, but this is not necessarily the case, especially for the
anisotropic and inhomogeneous turbulence associated with the three
instabilities. Observationally, this means that the turbulence parameters
associated with mass accretion, dust distributions, and line broadening can
all be different, even when they have the same physical origin.

###  2.1 Vertical shear instability

The VSI (_Arlt and Urpin_ , 2004; _Nelson et al._ , 2013) is a linear,
axisymmetric instability that is a disk-variation of the Goldreich-Schubert-
Fricke instability in differentially rotating stars (GSFI, _Goldreich and
Schubert_ , 1967; _Fricke_ , 1968). The VSI can be considered as the
$\mathrm{Pr}\rightarrow 0$ limit of the GSFI (_Knobloch and Spruit_ , 1982),
where the Prandtl number $\mathrm{Pr}$ is the ratio of viscosity to thermal
diffusivity.

The VSI is powered by the disk’s vertical differential rotation, i.e.
$\partial_{z}\Omega\neq 0$, which is exhibited by any baroclinic equilibrium
wherein surfaces of constant density and pressure are misaligned
($\nabla\rho_{\mathrm{g}}\times\nabla P\neq\bm{0}$). This includes vertically
isothermal disks with a radial temperature gradient ($q\neq 0$), for which one
finds

$\displaystyle R\frac{\partial\Omega^{2}}{\partial z}\simeq
q\varepsilon\frac{z}{H}\Omega_{\mathrm{K}}^{2}.$ (8)

Although vertical shear is weak
($\left|R\partial_{z}\Omega\right|\sim\varepsilon\Omega_{\mathrm{K}}\ll\Omega_{\mathrm{K}}$),
the associated free energy can be released through radially-short, vertically-
elongated disturbances (_Umurhan et al._ , 2013) provided that the stabilizing
effect of vertical buoyancy is either absent ($N_{z}^{2}=0$) or eliminated by
rapid cooling when $N_{z}^{2}>0$ (_Nelson et al._ , 2013), as is also the case
for the GSFI when Pr$\rightarrow 0$ (_Knobloch and Spruit_ , 1982). In the
case $N_{z}^{2}>0$, _Lin and Youdin_ (2015) find the criterion

$\displaystyle\tau_{\mathrm{cool}}\lesssim\frac{|q|\varepsilon}{\gamma-1}\Omega_{\mathrm{K}}^{-1},$
(9)

or $\beta_{\mathrm{cool}}\lesssim\varepsilon\ll 1$ for typical disk
parameters. That is, cooling timescales must be much shorter than the
dynamical timescale.

Linear VSI theory was initially developed in the local approximation for
axisymmetric disturbances (_Urpin and Brandenburg_ , 1998; _Urpin_ , 2003).
More recently, non-axisymmetric (_Volponi_ , 2014; _Latter and Papaloizou_ ,
2018) and radially-local, vertically global models have been developed
(_Nelson et al._ , 2013; _Barker and Latter_ , 2015; _McNally and Pessah_ ,
2015; _Lin and Youdin_ , 2015; _Umurhan et al._ , 2016b). The linear
instability mechanism can be explained in terms of angular momentum,
vorticity, and energetic considerations (_Yellin-Bergovoy et al._ , 2021).
Weakly non-linear theories have been developed by _Latter and Papaloizou_
(2018) and _Shtemler and Mond_ (2019, 2020).

Semi-global linear analyses find:  1) destabilized inertial waves or ‘body
modes’ that permeate the disk column, and  2) ‘surface modes’ concentrated
near the vertical boundaries.  Body modes tend to dominate nonlinear
simulations, which are typically characterized by narrow bands of large-scale
vertical motions (_Nelson et al._ , 2013), similar to that exhibited by low-
order body modes with widths $\sim\varepsilon H$. Interestingly, local, linear
VSI modes are also nonlinear solutions, which partly explains why they can
persist with large amplitudes (_Latter and Papaloizou_ , 2018).

Nevertheless, linear VSI modes are subject to secondary, Kelvin-Helmholtz-type
instabilities, which are often associated with vortex formation (_Nelson et
al._ , 2013; _Latter and Papaloizou_ , 2018). Indeed, vortices – non-
axisymmetric flow features – are readily observed in full 3D simulations
(_Richard et al._ , 2016; _Manger and Klahr_ , 2018), although there it has
been attributed to the Rossby Wave Instability (RWI _Lovelace et al._ , 1999;
_Li et al._ , 2000, 2001) of narrow vorticity rings directly formed by the VSI
or pressure bumps that arise from radial variations in the rate of angular
momentum transport by the VSI. These vortices may be subsequently amplified by
the Sub-critical Baroclinic Instability (_Klahr and Bodenheimer_ , 2003;
_Petersen et al._ , 2007a, b; _Lesur and Papaloizou_ , 2010, see also §2.2)
for steep radial temperature profiles and moderate cooling timescales; or
attacked by 3D Elliptic Instabilities (_Kerswell_ , 2002; _Lesur and
Papaloizou_ , 2009), especially if their aspect-ratios $\chi$ are small
($\lesssim 4$). Recent simulations with beta and radiative cooling find both
compact, short-lived (a few orbits) vortices and large-scale, elongated
($\chi\gtrsim 8$) vortices that survive for hundreds of orbits (_Manger et
al._ , 2020; _Flock et al._ , 2020). The latter can arise from pressure bumps
at the boundary of ‘VSI dead zones’ (_Flock et al._ , 2020), which develops
from the strong dependence of the VSI on the cooling time that leads to a
situation where the VSI does not operate at small radii ($\lesssim 10$ au)
while being active at larger radii. High resolution simulations also find
meridional vortices (_Flores-Rivera et al._ , 2020).

VSI-active disks are turbulent with moderate outwards angular momentum
transport: the radial stress-to-pressure ratio, or an alpha parameter,
$\alpha_{\mathrm{S}}$ of order $10^{-4}$ is found for $\varepsilon=0.05$ and
is an increasing function of $\varepsilon$ (_Manger et al._ , 2020). However,
VSI-turbulence is strongly anisotropic because typically $|v_{z}/v_{R}|\gg 1$.
This results in a vertical stress-to-pressure ratio of
$\alpha_{z}\equiv\langle\rho_{\mathrm{g}}v_{z}\delta v_{\phi}\rangle/P\sim
10^{-2}$ (here $\delta$ denotes the deviation from equilibrium), i.e. two
orders of magnitude larger than radial transport (_Stoll et al._ , 2017a).
Furthermore, both $\alpha_{\mathrm{S}}$ and $\alpha_{z}$ are non-uniform.
These are important distinctions from standard, constant-alpha disks as the
anisotropy of the VSI implies that the gas’ and dust’s radial and vertical
evolution can differ significantly.

The basic properties of the VSI can be captured by locally isothermal disks. A
first approximation to modeling the VSI in more realistic PPDs is to use the
beta cooling prescription with an estimate of physical cooling timescales (see
§2.4) in linear theory (_Lin and Youdin_ , 2015; _Fukuhara et al._ , 2021) or
numerical simulations (_Pfeil and Klahr_ , 2021). Beyond this, one can
explicitly account for radiative losses and stellar irradiation by solving the
radiative hydrodynamic equations (_Stoll and Kley_ , 2014; _Flock et al._ ,
2017b; _Flock et al._ , 2020). Such models show that the VSI can indeed be
expected in the outer regions of PPDs with an $\alpha_{\mathrm{S}}\sim
10^{-4}$ and flow structures largely consistent with locally isothermal
models. Turbulent velocities on the order of $10^{-2}{c_{s}}$ are found in the
midplane regions ; whereas in the coronal regions ($|z|\gtrsim 0.2R$, or
roughly beyond two pressure scale heights) these increase to $\sim 0.06c_{s}$
and can reach $O(10^{-1}c_{s})$ at some points (_Flock et al._ , 2017b). These
are consistent with the upper limits obtained from ALMA observations of
turbulent line widths (§5.2). However, identifying VSI turbulence this way
appears difficult owing to the VSI’s strong anisotropy; instead, its vertical
motions can be detected by spatially resolved observations of the CO line
emission (_Barraza-Alfaro et al._ , 2021).

Recent studies have also begun to incorporate the VSI in models of
planetesimal formation and planet evolution. Axisymmetric vorticity
perturbations induced by the VSI, and discrete vortices in the disk plane,
have been shown to act as effective traps of mm dust (_Stoll and Kley_ , 2016;
_Flock et al._ , 2017b; _Flock et al._ , 2020). The former gives rise to
prominent ring-like structures, while the latter produces localized
enhancements of dust. Large-scale vertical gas motions induced by the VSI can
effectively lift up passive dust grains to the disk atmosphere (_Stoll and
Kley_ , 2016; _Flock et al._ , 2017b; _Flock et al._ , 2020). When dust-on-gas
feedback is accounted for, solids can settle against VSI-stirring if the
overall dust abundance is sufficiently super-solar (_Lin_ , 2019). The dense
midplane dust layer can further undergo the Streaming Instability (_Youdin and
Goodman_ , 2005, see also §3), in which case the VSI can enhance radial dust
concentrations (_Schäfer et al._ , 2020). Following planet formation, disk-
planet interaction and pebble accretion with VSI turbulence have been
investigated by _Stoll et al._ (2017b) and _Picogna et al._ (2018),
respectively. For this particular problem they find results similar to viscous
disk models with an equivalent $\alpha_{\mathrm{S}}$ as measured in VSI runs,
but in the case of pebble accretion one also needs to include stochastic kicks
onto the solids.

Finally, while the VSI is fundamentally a hydrodynamic instability, PPDs are
magnetized. In ideal MHD, a strong field with a plasma beta parameter $\gtrsim
1/\varepsilon$ can stabilize the VSI (_Latter and Papaloizou_ , 2018; _Cui and
Lin_ , 2021), which may apply to the well-ionized disk atmospheres. On the
other hand, in the disk bulk where non-ideal effects dominate, the VSI
persists and can co-exist with magnetized disk winds (_Cui and Bai_ , 2020).
The dynamics of dust grains under such conditions remains to be investigated.

###  2.2 Convective overstability and baroclinic vortex amplification

The COS is a linear, axisymmetric instability that can be interpreted as
radial convection mediated by cooling (_Klahr and Hubbard_ , 2014; _Lyra_ ,
2014). Unstratified models show that the COS requires  1) a radial buoyancy
such that $N_{R}^{2}<0$, and  2) cooling on dynamical timescales,
$\beta_{\mathrm{cool}}\sim 1$.  For power-law disks the first criterion
becomes

$\displaystyle-1<\frac{p}{q}<\frac{1}{\gamma-1}\quad\text{for $N_{R}^{2}<0$}.$
(10)

This analytical criterion applies to the disk midplane where $N_{z}^{2}=0$.
For $q>0$, Eq. 10 translates to
$-(q+3)/2<\sigma<[q(\gamma+1)/(\gamma-1)-3]/2$, implying a shallow or even
rising surface density profile. For example, in irradiated disks with $q=0.5$,
one requires $-1.75<\sigma<0$ for $\gamma=7/5$. Such atypical profiles may be
realized at special locations such as gap edges or the outer edge of dead
zones (e.g. _Delage et al._ , 2022).

However, even if $N_{R}^{2}>0$ at $z=0$, in typical disk models the fact that
$\varepsilon$ increases with $R$ can drive $N_{R}^{2}<0$ away from the
midplane, where $N_{z}^{2}>0$. Unlike the VSI, though, vertical buoyancy is
not expected to suppress COS modes, which have little vertical motion (_Lyra_
, 2014). A similar situation applies to the symmetric instability in the
Earth’s atmosphere (_Schultz and Schumacher_ , 1999). Nevertheless, the
existence and saturation of the COS in vertically stratified disks off the
midplane still need to be demonstrated explicitly, which is not trivial as
potentially also VSI can be operative in this region.

The COS corresponds to growing epicyclic oscillations. In a thin disk,
rotation is strongly stabilizing. This means that, even if a gas parcel is
radially buoyant ($N_{R}^{2}<0$), it will still perform radial oscillations on
the epicyclic frequency. Indeed, in the absence of cooling the SH conditions
imply stability as $\kappa^{2}\gg|N_{R}^{2}|$ since typically
$|N_{R}|\sim\varepsilon\Omega_{\mathrm{K}}$. However, if cooling occurs on the
timescale on the order of the orbital period, the gas parcel can exchange heat
with the surrounding such that it experiences a buoyant acceleration, thereby
increasing its oscillation amplitude (_Latter_ , 2016). On the other hand, if
cooling is too rapid then the parcel experiences no buoyancy and epicycles are
again stable. Ultimately, the COS is mainly powered by the unstable entropy
gradient (cf. vertical shear for the VSI). Also unlike the VSI, the COS
prefers vertically short, radially extended length scales.

The linear analysis for the COS was possible as there are, just like for the
VSI, axisymmetric modes which are unstable (_Klahr and Hubbard_ , 2014; _Lyra_
, 2014; _Latter_ , 2016; _Volponi_ , 2016). Weakly nonlinear theory (_Latter_
, 2016) suggested that COS modes might be limited to small amplitudes due to
secondary parasitic instabilities. However, recent axisymmetric numerical
simulations (_Teed and Latter_ , 2021) showed a progression in the nonlinear
regime with increased Reynolds number as the system transitions from a weakly
nonlinear state with relatively ordered waves, to wave turbulence, and finally
to the formation of intermittent and then persistent zonal flows. This process
is likely what spawned vortices in 3D, non-axisymmetric models of the COS
(_Lyra_ , 2014). The axisymmetric simulations of _Teed and Latter_ (2021) also
showed the existence of ”elevator modes”, in which a mode of uniform vertical
velocity dominates the vertical extent of the (unstratified) box. Whether or
not this leads to a large-scale meridional circulation in a stratified, global
disk needs to be clarified.

Yet, well before the discovery of the axisymmetric, linear COS, the
amplification of vortices had already been observed in simulations of radially
stratified disks with cooling (_Klahr and Bodenheimer_ , 2003; _Petersen et
al._ , 2007a, b; _Lesur and Papaloizou_ , 2010; _Lyra and Klahr_ , 2011).
These authors find that such amplification of the vortices also requires
$N_{R}^{2}<0$ and cooling that is neither too fast nor slow. _Lesur and
Papaloizou_ described this effect as a Sub-critical Baroclinic Instability
(SBI), where ‘sub-critical’ refers to the fact that a vortex is not an
infinitesimal perturbation in the flow, but needs a finite initial size to
have the amplification be stronger than viscosity and shear. That is, SBI is a
non-linear phenomenon. In addition, the SBI operates in 2D, thin disks, unlike
the linear, axisymmetric COS that requires the vertical dimension.

Ultimately, however, the criteria for COS and SBI are qualitatively the same,
suggesting that both instabilities share the same physical driving, i.e.
radial buoyancy and cooling rates on their individual radial oscillation
frequencies. Whereas for the linear COS, this frequency is the epicyclic
frequency ($\kappa$); vortices typically rotate much slower than the disk
rotation, with a turnover frequency of order $\Omega_{\mathrm{K}}/\chi$
(recall $\chi$ is the vortex aspect ratio). Thus the larger the vortex aspect
ratio, the longer the rotation time, which gives a wide range of cooling times
that will lead to the amplification of a range of vortex shapes (_Raettig et
al._ , 2013).

Global numerical simulations of the SBI in 2D and 3D disks readily show the
development of large-scale vortices that can persist for hundreds of orbits
(_Klahr and Bodenheimer_ , 2003; _Petersen et al._ , 2007a, b; _Lyra and
Klahr_ , 2011; _Barge et al._ , 2016). Simulations also find that COS/SBI
leads to outwards angular momentum transport at a typical level of
$\alpha_{\mathrm{SS}}\sim 10^{-3}$, but can vary by an order of magnitude
depending on buoyancy and cooling parameters (_Raettig et al._ , 2013). Local
2D models suggest that vortices can be amplified even when radial buoyancy is
only weakly unstable (_Raettig et al._ , 2013). High resolution, local 3D
simulations show that SBI growth is balanced by secondary, 3D Elliptic
Instabilities (_Lesur and Papaloizou_ , 2009), resulting in a vortex with a
turbulent core (_Lesur and Papaloizou_ , 2010).

Nevertheless, the ability for COS-induced or SBI-amplified vortices to act as
efficient dust traps have been demonstrated by _Raettig et al._ (2015) in 2D,
and _Lyra et al._ (2018) and _Raettig et al._ (2021) in 3D. On the other hand,
_Lobo Gomes et al._ (2015) find the interaction between massive planets and
radially buoyantly unstable, cooling disks renders $N_{R}^{2}>0$, in which
case planet-induced gaps become smoother and leads to weaker vortices that are
associated with the RWI of the gap edge, rather than COS or SBI. However,
models of planetesimal formation and disk-planet interaction that explicitly
incorporate COS/SBI-induced turbulence and vortices are still lacking (cf. as
being developed for the VSI), which should be pursued in the future,
especially given the relevance of vortices for planet formation, and the
predicted occurrence of the COS in the 1-10 AU region.

###  2.3 Zombie vortex instability

Unlike the VSI and the COV, the ZVI is a non-axisymmetric, non-linear
instability, the latter meaning that it requires finite-amplitude
perturbations. The ZVI can develop if there is sufficient vertical buoyancy.
In practice, this translates to  1) $N_{z}^{2}\gtrsim\Omega_{\mathrm{K}}^{2}$,
which is generally the case for $|z|\gtrsim 1.5H$; and  2)
$\beta_{\mathrm{cool}}\gg 1$, otherwise buoyancy is diminished by cooling.

_Barranco and Marcus_ (2005) first observed coherent, long-lived vortices
naturally form in the strongly stratified regions above and below the PPD
midplane in their local, shearing box simulations with a spectral method. At
first, they hypothesized that the vortices were created by breaking internal
gravity waves, but _Marcus et al._ (2013), _Marcus et al._ (2015), and _Marcus
et al._ (2016) demonstrated the existence of a new type of purely hydrodynamic
instability associated with “baroclinic critical layers”: narrow structures in
a stratified shear flow where a wave’s phase speed matches the shear velocity
plus/minus the vertical buoyancy frequency divided by the wavenumber. They
named the instability the “Zombie Vortex Instability” (ZVI) not only because
it may occur in dead zones, but also because of the way one zombie vortex
“infects” neighboring baroclinic critical layers, spawning new zombie
vortices, which “infect” farther critical layers, and so on, filling a region
of the disk with vigorous turbulence. ZVI is not an artifact of the numerical
method as it was computed with spectral codes and finite-volume codes, with
Boussinesq, anelastic, and fully compressible treatments of the continuity
condition, and with and without the shearing box. One may ask, if it is so
robust, how was it missed in previous numerical calculations? Prior work often
lacked one or more of the crucial ingredients: ZVI requires strong vertical
stratification, high resolution to resolve the narrow critical layers
($\gtrsim 256$ cells or spectral modes per $H$), a broad spectrum of
perturbations (i.e., Kolmogorov, but not Gaussian-peaked, so that the
vorticity peaks on the small scales), and enough simulation time to allow the
critical layers to amplify perturbations.

_Lesur and Latter_ (2016) confirmed the existence of the instability with
their own spectral simulations, but also pointed out the sensitivity of the
ZVI to dissipation processes and small scale physics: fast cooling and/or
sufficient viscosity suppresses the ZVI in numerical simulations (with
Reynolds number $\mathrm{Re}\lesssim 10^{6}$). However, the relevant Reynolds
number in PPDs is $\mathrm{Re}\sim 10^{14}$, so viscosity may not inhibit the
ZVI in PPDs, but does make it difficult to produce ZVI in a laboratory (_Wang_
, 2016). _Barranco et al._ (2018) investigated ZVI under a wider range of
realistic scenarios, including with nonuniform vertical stratification and
radiative cooling, and demonstrated that the off-midplane regions throughout
the planet-forming regions of protoplanetary disks may be susceptible to ZVI
so long as the cooling time is longer than a few orbital periods, which may be
likely with modest grain settling and/or growth. _Umurhan et al._ (2016a) and
_Wang and Balmforth_ (2020) have undertaken theoretical investigations into
the earliest phases of the instability considering both linear growth and the
nonlinear dynamics within the baroclinic critical layers. Most recently, _Wang
and Balmforth_ (2021) developed a reduced model for the nonlinear dynamics
within forced baroclinic critical layers and have demonstrated that
perturbations grow secularly, generating jet-like defects within the shear
which are then later susceptible to secondary instabilities that yield
coherent vortical structures.

Turbulence from ZVI has some unique characteristics that may be highly
relevant to dust transport. At late times, ZVI turbulence results in
significant vertical mixing that tends to homogenize the background
stratification resulting in a nearly uniform $N_{z}$. Similar behavior is seen
in atmospheric and oceanic flows in which the breaking of internal gravity
waves creates step-like or staircase patterns in stratification (_Orlanski and
Bryan_ , 1969; _Phillips_ , 1972; _Pelegrí and Sangrã_ , 1998). While the
region in the immediate vicinity of the midplane lacks the requisite
stratification for the excitation of baroclinic critical layers, _Barranco et
al._ (2018) observed that zombie turbulence from the ZVI susceptible regions
can penetrate into the midplane, albeit with a smaller magnitude. At late
times, ZVI turbulence resulted in the creation of azimuthal quasi-steady-state
zonal flows. The zonal flows consisted of 5-6 pairs of dipolar vortex layers
within a radial extent of $8H$ (see the third figure in Table 1). The radial
thickness of the cyclonic layers is approximately one-third the thickness of
the anticyclonic layers, but this is compensated by the fact that the cyclonic
layers are roughly three times more intense than the anticyclonic layers.
Fully-developed zombie turbulence shows intermittency where the flow cycles
through near-laminar phases of zonal flow punctuated by chaotic bursts of new
zombie vortices. In some simulations, the bursting is quasi-periodic in time
(with periods between 100-150 orbits), whereas in other cases, it appeared
more stochastically.

The ultimate impact of ZVI on PPD dynamics is still to be demonstrated. Future
topics include: (1) how efficient is ZVI at transporting angular momentum in
compressible simulations ($\alpha_{\mathrm{S}}$); (2) the properties of ZVI
turbulence as a function of the Reynolds number, especially beyond
$O(10^{7})$; (3) how does ZVI turbulence transport and mix dust and how does
dust affect the ZVI, as well its mutual interaction with planetesimals and
planets; (4) what does the ZVI look like in global simulations; and (5) what
are the observational signatures of ZVI turbulence, e.g., broadening of
molecular lines, and can that be distinguished from turbulence from other
purely hydrodynamic instabilities.

###  2.4 Occurrence in physical disk models

As Table 1 summarizes, the onset of the above hydrodynamic instabilities
depends on both the disk structure and thermal relaxation or cooling times
($\tau_{{\rm cool}}$), which are intimately interdependent on one another. As
discussed in §2.1 – 2.3, the structural requirements for the VSI,
$\partial_{z}\Omega\neq 0$, can be met by a non-zero radial temperature
gradient ($q\neq 0$), that for the COV where the radial buoyancy
$N_{R}^{2}<0$, and that for the ZVI where the vertical buoyancy
$|N_{z}|>\Omega_{\mathrm{K}}$. As for their thermodynamic requirements, with a
broad brush stroke, we can say that VSI requires cooling times significantly
shorter than an orbital period, COS needs cooling times of order the orbital
period, whereas ZVI is operable when the cooling time is longer than a few
orbital periods. In this section, we describe the ingredients for estimating
cooling times and discuss recent efforts to do so in order to assess the
occurrence of the above instabilities. We will see that at the current stage,
results are still model-dependent as they are subject to a number of
uncertainties.

Establishing the relevant cooling times for the mechanism in question depends
upon the spatial lengthscale of the fastest growing disturbance $\ell_{m}$ and
how $\ell_{m}$ compares to the photon mean-free path $\ell_{{\rm
ph}}=1\big{/}(\kappa_{\mathrm{op}}\rho)$, where $\kappa_{\mathrm{op}}$ is the
material’s frequency integrated opacity. This places the dynamics in either
the optically thin ($\ell_{{\rm m}}\ll\ell_{{\rm ph}}$) or thick ($\ell_{{\rm
m}}\gg\ell_{{\rm ph}}$) regime (e.g., see the discussion in _Barranco et al._
, 2018). Since H2 and other molecules are inefficient radiators compared to
dust particles (however molecular line cooling may play a role in the disk’s
coldest regions, e.g., _Woitke et al._ , 2016), the primary cooling pathway
involves the collisional transfer of thermal energy from gas to particles, the
latter of which efficiently radiates (e.g., see discussions in _Barranco et
al._ , 2018; _Pfeil and Klahr_ , 2019; _Lyra and Umurhan_ , 2019).

In practice, $\tau_{{\rm cool}}$ can be thought of as the longer timescale of
the particle-gas collisional energy exchange timescale $\tau_{{\rm gas}}^{{\rm
col}}$ and the radiative timescale $\tau_{{}_{\rm cool}}^{{\rm rad}}$, i.e.,
$\tau_{{\rm cool}}={{\rm max}}\left(\tau_{{\rm cool}}^{{\rm rad}},\tau_{{\rm
gas}}^{{\rm col}}\right)$. Furthermore, $\tau_{{\rm cool}}^{{\rm rad}}$
depends upon whether the lengthscale of interest is in the optically thick or
thin regimes embodied in the expression (_Spiegel_ , 1957)

$\frac{1}{\tau_{{\rm cool}}^{{\rm rad}}}=\frac{\mu}{\tau_{{\rm cool}}^{{\rm
thin}}},\qquad\mu\equiv\left[1-\frac{\tan^{-1}\left(2\pi\ell_{\mathrm{ph}}/\ell_{m}\right)}{2\pi\ell_{\mathrm{ph}}/\ell_{m}}\right];$
(11)

where the frequency integrated optically thin cooling timescale is given by
$\tau_{{\rm cool}}^{{\rm thin}}\equiv{\ell_{{\rm
ph}}\rho_{g}C_{P}}\big{/}{16\sigma_{\mathrm{B}}T^{3}},$ and
$\sigma_{\mathrm{B}}$ is the Stefan-Boltzmann constant. A further complication
is how to estimate the frequency-integrated opacity $\kappa_{\mathrm{op}}$ as
a function of disk location and epoch that, in turn, depends on the disk’s
turbulence-mitigated spatiotemporal particle distribution $n(a)$, where $a$ is
particle size (e.g. _Birnstiel et al._ , 2012; _Estrada et al._ , 2016). This
reveals a causality dilemma as the key determinant driving disk hydrodynamic
instabilities, i.e. $\tau_{{\rm cool}}\left(\kappa_{{\rm
op}},n(a),\cdots\right)$ and associated disk structure, depends on the
turbulence generated by the very same instabilities. These, in turn, are
sensitive to the assumed disk mass and maximum particle size as a function of
disk epoch – two quantities that are currently only weakly constrained by
observations.

The currently adopted strategy in confronting the problem of determining
$\kappa_{{\rm op}}$, and subsequently $\tau_{{\rm cool}}$, is to assume an
alpha-type viscous disk model in which vertical particle settling is balanced
by upward turbulent transport (_Dubrulle et al._ , 1995). Together with
evolutionary models of dust mass and size distributions – and their resulting
frequency-integrated Rosseland or Planck mean opacities (c.f., _Semenov et
al._ , 2003; _Woitke et al._ , 2016; _Cuzzi et al._ , 2014) – several recent
studies have made preliminary forays into mapping out disk hydrodynamic
instabilities using the derived thermal relaxation times based on the
aforementioned opacity models (_Malygin et al._ , 2017; _Barranco et al._ ,
2018; _Pfeil and Klahr_ , 2019; _Fukuhara et al._ , 2021). Each of these
studies extracts various dynamical and thermodynamical profiles resulting from
and/or are inputs to various 1+1 global disk evolutionary models – e.g.,
quantities like density, pressure, total disk mass ($M_{{\rm disk}}$),
particle and gas surface densities ($\Sigma_{g},\Sigma_{d}$ respectively),
particle size distribution, the turbulent $\alpha_{\mathrm{S}}$ for gas
accretion, the turbulent $\alpha_{\mathrm{D}}$ for dust diffusion, etc., all
as functions of position and disk age. One can then assess whether the
structure and thermodynamic criteria for the onset of the instabilities
discussed in §2.1–2.3 are met, at least marginally.

_Pfeil and Klahr_ (2019), following the framework described in _Malygin et
al._ (2017), assume a submicron-sized monodisperse population of grains and
primarily analyze disks with mass $M_{\mathrm{disk}}=0.1M_{\odot}$. Fig 1a is
a corresponding map where their results predict that the VSI is widespread
across the vast majority of the disk (1AU$<R<$100AU), with some coincidental
potential for active COS for 1AU$<R<$10AU in regions containing the midplane
and for regions extending away from the midplane for $R>$10AU. This study
implies that $\tau_{{\rm cool}}$ is limited by $\tau_{{\rm cool}}^{{\rm rad}}$
as the grain number density and disk mass is sufficiently high that the
collision time is considerably short, i.e., $\tau_{{\rm cool}}^{{\rm
rad}}\gg\tau_{{\rm gas}}^{{\rm col}}$.

Figure 1: Occurrence of hydrodynamic instabilities in different PPD models.
(a) From _Pfeil and Klahr_ (2019). This monodisperse ($a\approx 1\mu$m) grain
population model considers $M_{{\rm disk}}=0.1M_{\odot}$ and
$\alpha_{\mathrm{S}}=10^{-3}$. Black lines delineate successive pressure scale
heights. In their models, a Vertical Convective Instability (VCI) can develop
because $N_{z}^{2}<0$, but this is not considered in the main text. (b) From
_Fukuhara et al._ (2021). Panels (i)–(iii) show VSI activity for
$M_{\mathrm{disk}}/M_{\sun}=10^{-2},\,10^{-3}$,
$\Sigma_{d}/\Sigma_{g}=10^{-3},10^{-2}$, and polydisperse grains with $a_{{\rm
max}}=1$mm, $s=3.5$, $\alpha_{\mathrm{D}}=10^{-4}$. Panel (iv) shows
$\tau_{{}_{\rm cool}}^{{\rm col}}/\tau_{{}_{\rm cool}}^{{\rm rad}}$ for
$a_{{\rm max}}=10\mu$m, $M_{{\rm disk}}=0.01M_{\odot}$, and
$\alpha_{\mathrm{D}}=10^{-4}$. Dashed black lines delineate successive
pressure scale heights. (c) From _Barranco et al._ (2018). ZVI-active regions
are shaded. Black contours are $\log(\Omega\tau_{{\rm cool}})$ while dotted
magenta contours signify $\log(\Omega\tau_{{}_{\rm cool}}^{{\rm rad}})$. Here
$M_{{\rm disk}}\approx 0.04M_{\odot}$, $\Sigma_{d}/\Sigma_{g}=0.01$, and
grains have $a_{{\rm max}}=10$ cm and $s=3.25$. Panels correspond to different
$\alpha_{\mathrm{D}}$ (labeled as $\alpha$). Simulations with vertical extent
and at the radii of the red lines develop ZVI turbulence, while those of the
blue dashed lines do not. Note that in (a) $\alpha_{\mathrm{S}}$ parameterizes
mass transport, while in (b) and (c) $\alpha_{\mathrm{D}}$ parameterizes the
dust settling.

_Fukuhara et al._ (2021) treat disk models with a polydisperse population of
grains with a number density size distribution $n(a)\propto a^{-s}$ with a
constant power-law index $s=3.5$ between a minimum and maximum grain sizes
$a_{{\rm min}}$ and $a_{{\rm max}}$, respectively. From this they construct an
effective mean molecule travel length $\ell_{gd}$, from which they estimate
the gas-grain collisional energy exchange timescale in terms of a simpler
collisional timescale, with $\tau_{{\rm gas}}^{{\rm
col}}\equiv\ell_{gd}/v_{{\rm th}}$, where $v_{{\rm th}}$ is a typical thermal
velocity ($\propto T^{1/2}$). Their analysis reveals the useful approximate
form $\ell_{gd}\propto a_{{\rm
max}}^{1/2}\Sigma_{g}^{-1/2}(\Sigma_{d}/\Sigma_{g})^{-1}\propto M_{{\rm
disk}}^{-1/2}(\Sigma_{d}/\Sigma_{g})^{-1}$. As noted earlier in _Barranco et
al._ (2018), as grains get incorporated into larger particles, the total
effective collisional cross-section of the grain population decreases. Thus, a
molecule must travel further before the next collision, which leads to
increasingly inefficient cooling. However, this reduction in the local
medium’s effective collisional cross-section may be halted if particles grow
as porous fractal aggregates instead of as zero-porosity spheres (_Okuzumi et
al._ , 2012).

For a range of $M_{{\rm disk}}$ values, _Fukuhara et al._ (2021) examine the
onset of the VSI, finding that it is absent in low mass disks ($M_{{\rm
disk}}=0.001M_{\odot}$), while for relatively high disk masses ($M_{{\rm
disk}}=0.1M_{\odot}$) it can potentially take root in fairly narrow zones
sandwiching the midplane and out to the disk’s outer edge ($|z|/R\lessapprox
0.05-0.1$, see also Fig 1b). They show how the suppression of the VSI,
especially in the low disk mass models, is due to the inability of gas
molecules to efficiently communicate with the dust particles causing
$\beta_{\mathrm{cool}}$ to greatly exceed the critical value of 0.1. However
their results would suggest that the increased cooling times could make their
disks susceptible to the COS or ZVI instead.

In a prior study, _Barranco et al._ (2018) take a more sophisticated approach
to estimating cooling times by explicitly considering the finite time for
energy to be exchanged between gas molecules and dust grains via collisions
(e.g., _Hollenbach and McKee_ , 1979; _Burke and Hollenbach_ , 1983;
_Glassgold et al._ , 2004). They allow for a variable power-law index
($3<s<4$) in the particle size distribution and find $\tau_{{\rm gas}}^{{\rm
col}}\propto a_{{\rm S}}\rho_{g}^{-1}T^{-1/2}(\Sigma_{d}/\Sigma_{g})^{-1}$,
where $a_{{\rm S}}$ is the Sauter-mean radius of grains (i.e. the radius of a
monodisperse population of grains that would have the same total surface area
and total volume of a given polydisperse population of grains). Their study
makes predictions for a Minimum Mass Solar Nebula disk model with $M_{{\rm
disk}}\approx 0.042M_{\odot}$ (_Cuzzi et al._ , 1993) for a variety of
$\alpha_{\mathrm{SS}}$ and $a_{{\rm max}}$ values, finding under these
parameter conditions that $\tau_{{\rm cool}}$ is generally controlled by
$\tau_{{\rm gas}}^{{\rm col}}$ . Fig 1c exhibits a set of results for large
values of $a_{{\rm max}}$ ($\sim 10$ cm) in which the VSI is entirely ruled
out as $\beta_{\mathrm{cool}}\geq 0.1$ everywhere in the model, while allowing
for the feasibility of the ZVI generally away from the midplane where
$\beta_{\mathrm{cool}}\gtrapprox 20$, and possibly the COS in midplane layers
where $\beta_{\mathrm{cool}}\sim 1$. For smaller $a_{{\rm max}}$ ($\sim 1$ cm)
their results indicate that the VSI is feasible closer to the midplane and
closer to the star ($<$ 10 AU) where $\beta_{\mathrm{cool}}\leq 0.1$.

The general trends reported in _Fukuhara et al._ (2021) – especially the
limiting behavior posed by collisional exchange timescales (e.g., see bottom
right panel of Fig. 1b) – are largely consistent with the findings of
_Barranco et al._ (2018); while their mutual differences are likely due to the
use of different disk model parameters as well as in the approaches applied in
estimating $\tau_{{\rm gas}}^{{\rm col}}$. These are matters that need future
resolution.

In summary: a systematic examination is required to determine what type of
instability is operative when and where in a PPD. The strategies used in the
aforementioned studies ought to be applied toward more realistic global
evolution models that improve upon simplistic power-law formalisms. The maps
in Fig. 1 should not be taken as definitive but are meant to illustrate the
diverse outcomes depending on model assumptions. However, some trends may be
inferred from such recent studies: the picture of the disk during its earliest
stages – i.e., when the $a_{{\rm max}}$ are still small and $M_{{\rm disk}}$
still relatively high – may be more like Fig. 1a, while as the disk evolves
with $M_{{\rm disk}}$ decreasing and $a_{{\rm max}}$ inexorably growing the
turbulent state of the disk may be driven by processes designated like those
shown in Figs. 1b-c. In any case, given the possible range of $\beta_{{\rm
cool}}$ values likely relevant across the disk bulk, some type of
hydrodynamical instability that leads to turbulence – whether it be ZVI and/or
COS, and/or VSI (or perhaps something yet to be discovered) – is likely active
in PPDs.

As a final remark, we emphasize that alpha-type viscous disk models are
generally no replacement for direct numerical simulations for modeling the
turbulence generated by these hydrodynamic instabilities. It should be kept in
mind that PPDs are essentially inviscid but computational limitations prohibit
one to probe realistic Reynolds numbers. An important future direction is to
develop detailed models of the hydrodynamic turbulence generated by these
instabilities that are suitable for use in large-scale, long-term global
simulations of PPDs. Until then, alpha-type viscous disk models – while
convenient – should be interpreted with care.

##  3 MULTI-PHASE INSTABILITIES WITH DUST-GAS INTERACTIONS

###  3.1 Barriers to planetesimal formation in turbulent disks

Bridging the gap between dust grains and solid planetary cores in the core
accretion paradigm involves intermediate gravitationally-bound bodies of
typical sizes $\sim$0.1–100 km called planetesimals . Conversion of pebbles
into planetesimals requires for solids to overcome several physical barriers
as they grow, and is constrained by spatially resolved observations of young
disks. Strong concentrations of solids must be driven and retained in the
disk. Privileged places for large dust enrichments are long-lived pressure
maxima, where dust grains tend to drift into. Moreover, dust grains settle
towards the mid-plane of the disk. Turbulent of order
$\alpha_{\mathrm{D}}\sim\alpha_{\mathrm{S}}\gtrsim 10^{-5}$–$10^{-4}$ solid-
to-gas density ratio $\epsilon_{\mathrm{d}}\lesssim 1$ in the mid-plane may
prevent the gravitational collapse of the dust layer . This turbulence can be
hydrodynamical, MHD, or even of dust origin since the dust layer may itself be
Kelvin-Helmholtz unstable .

####  3.1.1 The radial-drift barrier

Disks need to somehow retain their solids as the solids grow in size and reach
, at which their radial drift onto the star becomes (Eq. (6)). Pure growth
could in principle assist grains for decoupling from the gas before they get
accreted, but bouncing or fragmentation probably limit this possibility in
practice (_Brauer et al._ 2008; _Birnstiel et al._ 2009; _Zsom et al._ 2010;
see also §3.1.2). Local dust traps offer powerful alternatives to concentrate
grains (see chapter by Bae et al. in the same book). Some traps are
axisymmetric. Among suggestions: disk boundaries, dead zones, snow-lines,
zonal flows, and edges of planetary gaps. the dust layer, , modifies the
transport of gas according to

$\displaystyle\frac{\partial\Sigma}{\partial
t}-\frac{1}{R}\frac{\partial}{\partial R}$
$\displaystyle\left[aZ\frac{\partial\left(c_{\mathrm{s}}^{2}\Sigma\right)}{\partial
R}\right]$ $\displaystyle+\frac{1}{R}\frac{\partial}{\partial
R}\left(\frac{1+Zb}{1+Z}R\,v_{\rm visc}\right)=0,$

where $a=1/\left[\mathrm{St}+(1+Z)^{2}\mathrm{St}^{-1}\right]$,
$b=\mathrm{St}^{2}/\left[\left(1+Z\right)^{2}+\mathrm{St}^{2}\right]$, and
vertical integration being restrained to the dust layer (_Gonzalez et al._ ,
2017). Relative intensities of are of order $\sim aZ/\alpha$ in smooth disks,
meaning that for $10^{-2}\lesssim\mathrm{St}\lesssim 10^{2}$ and a large
$Z\gtrsim 0.1$, gas dynamics can be dominated by viscosity, favoring the
formation of self-induced dust traps (_Gonzalez et al._ , 2017).

spirals, vortices and shocks induced by a companion (see chapter by Bae et al.
in the same book). Studying the resilience of these traps is key for
understanding planetesimal formation. The existence and the morphology of
these traps combined with the location of the dust outer radius, and/or the
remaining dust flux outside the observed traps provide indirect but
degenerated constraints on the parameters $\alpha_{\mathrm{D}}$, $\rm{St}$ and
$\epsilon_{\mathrm{d}}$ (_Birnstiel and Andrews_ , 2014; _Rosotti et al._ ,
2019). Models of dust traps should be compatible with the constraints provided
by structural data. Observation of silicate reveals crystallinity in the
grains in the cold outer regions , although it requires high temperatures to
develop. Explanations invoke outward drift mechanisms, redistribution of
material by , without any major consensus so far . Isotopic analysis reveals
that the early Solar System has likely been split into two parts in less than
$\sim$3 Myrs (_Kruijer et al._ , 2017). This puts a drastic constraint on the
existence of a dust trap in its early life, with subsequent consequences for
the formation of terrestrial and giant planets (_Desch et al._ , 2018;
_Kruijer et al._ , 2020).

####  3.1.2 The growth barriers

Lab experiments suggest that make dust aggregates grow up to a few centimeters
in size by hit-and-stick collisions. occur without restructuring at first, but
. Above a few centimeters in size, collisions are found to be on average
either non-adhesive, or destructive (_Blum_ , 2018), which is referred as the
bouncing, erosion or fragmentation barriers. Fragmentation occurs for
collision velocities of order $\sim$1–10 m.s-1, when monomer bonds are broken.
The fragmentation thresholds appear to for icy particles, but this is not in
general true and depends on temperature and surface coating (_Homma et al._ ,
2019; _Musiolik and Wurm_ , 2019; _Steinpilz et al._ , 2019; _Bischoff et al._
, 2020; _Arakawa and Krijt_ , 2021) The lab experiments appear to be in
agreement with a maximum size of order $\sim$1 mm for the size distribution of
dust grains analyzed on the Comet 67P/Churyumov-Gerasimenko (_Blum et al._ ,
2017a), and with a maximum size of $\sim$1 cm for meteoritic inclusions.

Besides the dust-gas instabilities discussed in the following sections, there
may be other possible pathways to overcome these growth barriers. One involves
rare low-velocity events that seed the growth of large bodies (_Windmark et
al._ , 2012; _Garaud et al._ , 2013; _Booth et al._ , 2018). However, the
growth timescales may be too long for solids to resist both the radial-drift
barrier and erosion (_Schräpler et al._ , 2018). Porous aggregates made of
small ultra-sticky icy monomers have also been invoked (_Potapov et al._ ,
2020). The robustness of this scenario remains to be validated experimentally
with matrix monomers of size 1–10 $\mu$m. Electrostatic charges such as
resulting from photoelectric and plasma charging were thought to inhibit
growth efficiency (_Okuzumi_ , 2009; _Akimkin et al._ , 2020). Triboelectric
charging may relieve this constraint, given collision velocities $\lesssim$0.1
m.s-1 (_Steinpilz et al._ , 2020b, a; _Teiser et al._ , 2021).

In any case, the lithospheric pressure of comets is only compatible with
material of low tensile strength. Experiments reveal that such tensile
strength can not originate from collisional growth through mass transfer, but
fits with a gentle collection of solids via a local gravitational collapse
(_Blum et al._ , 2017a). So far, observations have probed continuum emission
of solids up to a few centimetres (_Lommen et al._ , 2009; _Casassus et al._ ,
2019). Although the forthcoming Square Kilometre Array will allow to probe
larger wavelengths, it is not clear whether continuum emission of decimeter-
sized grains will be lost into free-free electron emission or not (_Dewdney et
al._ , 2009).

Another possibility to overcome the growth barriers is to concentrate grains
aerodynamically by means of dust-gas instabilities into the form of dusty
clouds, up to the stage where the local mass of solids is such that it
collapses gravitationally (_Ward_ , 2000; _Shariff and Cuzzi_ , 2011; _Youdin_
, 2011; _Shi and Chiang_ , 2013; _Latter and Rosca_ , 2017). We hereby
dedicate the rest of the section to this pathway.

###  3.2 Streaming instability & resonant drag instabilities

####  3.2.1 Linear phase

A linear analysis of the SI proceeds from the equations of motion for gas and
dust under the local-shearing-box approximation, treating dust as a
pressureless fluid that drifts inwards (see, e.g., _Youdin and Johansen_ ,
2007, Equations (1)–(4)). Considering a small, homogenous patch at $z=0$
within the midplane dust layer, one linearizes the equations to obtain
algebraic equations for each Faster-growing, larger-scale modes are considered
likely to have a stronger nonlinear influence, although the link between
linear-SI growth rates and planetesimal formation remains highly uncertain at
this time.

Figure 2: The basic ingredients involved in the mechanism of the streaming
instability. The feedback cycle (steps (1)–(4)) shows how an asymmetric clump
of dust can dynamically drive gas and dust flows that subsequently enhance the
density of the clump. We show the radial-azimuthal plane with dust density
shown using brown dots and the motions of gas and dust relative to the
background drifts shown with pink and brown arrows, respectively. During this
process, the dust continues its bulk inwards drift; the pictured deflections
are on top of this drift. For this feedback to operate, the perturbation must
have a vertical structure, otherwise radial gas motions in step (1) are halted
by the gas pressure because they are compressive. While the key physical
principles shown here are valid in both the low- and
high-$\epsilon_{\mathrm{d}}$ regimes, it is necessarily overly simplified;
more detailed versions are given in _Squire and Hopkins_ (2020).

##### Ideal streaming instability

Motivated by the study of _Goodman and Pindor_ (2000), the analysis described
above was first carried out by _Youdin and Goodman_ (2005), who discovered the
SI and suggested its importance to planetesimal formation. They showed that
the key parameters governing the instability are the dimensionless stopping
time $\tau_{\mathrm{s}}$ and the solid-to-gas density ratio
$\epsilon_{\mathrm{d}}$,111The pressure support parameter $\eta$ (defined by
_Nakagawa et al._ , 1986) sets the velocity scale via the differential speed
$\eta v_{\mathrm{K}}$ between gas and dust, but otherwise does not change the
properties of the instability. demonstrated that the motions induced by the SI
act to clump the dust, and argued that in most relevant regimes the growth
rates are large enough for the instability to reach its nonlinear stages well
before the radial dust drift significantly modifies the equilibrium. A key
feature uncovered in their analysis was the sudden increase in the SI’s growth
rate for $\epsilon_{\mathrm{d}}\gtrsim 1$.

Most of the analysis of _Youdin and Goodman_ (2005) solved for the eigenvalues
numerically, which can be inconvenient for making simple estimates across a
wide parameter range. We see that the SI grows faster – although at smaller
scales – for near vertical modes ($k_{z}\gg k_{r}$, $\theta_{k}\approx 0$),
and grows more slowly for smaller grains (smaller $\tau_{\mathrm{s}}$). At
$\epsilon_{\mathrm{d}}\gtrsim 1$, likely the more relevant regime for
planetesimal formation, _Squire and Hopkins_ (2018a) showed that the SI
changes its character and mechanism, . For $\epsilon_{\mathrm{d}}\gg 1$ and
$\tau_{\mathrm{s}}\ll 1$, its maximum growth rate is ${\rm
Im}(\omega)/\Omega\approx\sqrt{\epsilon_{\mathrm{d}}-1}$ at radial wavenumber
$k_{r}\eta r\approx 0.4\epsilon_{\mathrm{d}}^{19/8}\tau_{\mathrm{s}}^{-5/4}$.
This growth is very rapid (${\rm Im}(\omega)>\Omega$), with a rate that
increases monotonically with $\epsilon_{\mathrm{d}}$, although it operates at
very small scales for small $\tau_{\mathrm{s}}$.

Even in its most basic form, the mechanism for the streaming instability is
more complex than most astrophysical fluid instabilities. Various works have
presented simplified analyses, exploring reduced equations and the behavior of
the roots of the dispersion relation to understand how the instability
operates (see, e.g., _Jacquet et al._ , 2011; _Squire and Hopkins_ , 2018a;
_Zhuravlev_ , 2019; _Jaupart and Laibe_ , 2020; _Pan_ , 2020; _Squire and
Hopkins_ , 2020). These have shown that rotation, dust drift, and two-
dimensional structure are key to the SI, although the background shear flow is
not. The basic mechanism of its operation is sketched in Fig. 2, which shows
how a radial clump of dust can drive a gas flow that – due to the Coriolis
force – drives a radial drift of dust outwards (relative to the background
dust drift), thus enhancing the original clump. In order for the associated
gas flow to be incompressible, and thus not strongly resisted by gas pressure,
the flows must have vertical structure. . Note also that the gas pressure
perturbations, which are driven in step (2) , are not in phase with the dust-
density perturbations in most regimes, complicating the picture beyond dust
trapping in pressure bumps. This is consistent with _Lin and Youdin_ (2017)’s
finding that an isothermal gas embedded with small grains is equivalent to a
single fluid with a modified cooling, in which case a phase lag between
pressure and dust-density perturbations is a necessary feature of growing
oscillations, such as the SI. Interestingly, this feature seems to carry over
to the nonlinear regime, where an anti-correlation between dust and gas
density has been observed (see figure 7 of _Yang and Johansen_ 2014 as well as
_Li et al._ 2018). Thus, the conventional notion of dust-trapping at pressure
maxima does not necessarily apply under dynamical (i.e., non-steady)
conditions such as instabilities.

##### Extensions: polydisperse grains, turbulence, and other drag
instabilities

Any real disk will involve a wide spectrum of grain sizes
($\tau_{\mathrm{s}}$) at a given radial location, which is not taken into
account in the above analysis. _Krapp et al._ (2019) first explored this
“polydisperse” streaming instability, uncovering the unexpected result that in
regimes with $\epsilon_{\mathrm{d}}\lesssim 1$ the maximal growth rate of the
SI decreases monotonically as the number of discrete grain species is
increased. Very low growth rates, or only upper bounds in some cases, are
found with a true continuous distribution. Their results, which have been
further explored in _Zhu and Yang_ (2021) and _Paardekooper et al._ (2020,
2021), suggest that the polydisperse SI is primarily controlled by
$\epsilon_{\mathrm{d}}$ and the maximum $\tau_{\mathrm{s}}$ of the
distribution ($\tau_{{\rm s,max}}$). Specifically, _Zhu and Yang_ (2021) found
that the instability is robust (converged growth rates $\gtrsim 0.1\Omega$)
only above a sharp, distinct boundary in $\epsilon_{\mathrm{d}}$-$\tau_{{\rm
s,max}}$ space, with $\epsilon_{\mathrm{d}}\gtrsim 1$ or $\tau_{{\rm
s,max}}\gtrsim 1$ (see also §3.2.2). The results also suggest that
calculations involving polydisperse grains must carefully check convergence in
the number of grain species.

The presence of background turbulence, even at very low levels, may also
adversely affect the SI by diffusing small-scale motions and dust density
variations. Simple estimates can be obtained by modeling the effect of
turbulence as causing a gas viscosity, dust pressure, and dust density
diffusion, each parameterized by the disk $\alpha_{\mathrm{S}}$. The approach,
explored by _Umurhan et al._ (2020) and _Chen and Lin_ (2020), suggests that,
because turbulence decreases its growth rate and increases its spatial scale,
the SI can operate only relatively far out in the disk beyond $\simeq\\!10{\rm
AU}$, and only for very low (though realistic) levels of turbulence
$\alpha_{\mathrm{S}}\lesssim 10^{-4}$ and relatively large grains (see also
_Squire and Hopkins_ , 2018a; _Jaupart and Laibe_ , 2020; _Zhuravlev_ , 2020,
for further mathematical details). On the other hand, turbulence can also act
to clump grains itself, an effect that cannot be captured by diffusive
treatments (e.g., _Cuzzi et al._ , 2001; _Pan et al._ , 2011; _Hartlep et al._
, 2017; _Hartlep and Cuzzi_ , 2020).

Finally, it is worth noting several extensions to the standard SI analysis,
which have revealed related instabilities. _Pan and Yu_ (2020) considered non-
axisymmetric motions in the radial-azimuthal plane (perturbations without
vertical structure across the dust layer). They found a similar instability,
as expected from numerical simulations (_Schreiber and Klahr_ , 2018),
although there are some unique features and different scalings from the
axisymmetric SI. _Auffinger and Laibe_ (2018) found an interesting related
instability when the background profile is modified by a pressure bump,
although its details and relationship to the standard SI remain unclear. The
vertically global SI study of _Lin_ (2021) revealed a faster-growing
instability driven by the vertical gradient of the dusty gas’ rotation
velocity that is not captured by standard analyses (see also the unpublished
study of _Ishitsu et al._ , 2009). This may dominate the SI and control the
dust-layer thickness, showing interesting similarities to results of nonlinear
simulations; more investigation is needed. Finally, the RDI method of _Squire
and Hopkins_ (2018a) also uncovered a related instability, termed the
“settling instability,” which occurs as grains settle towards the midplane of
the disk (see also _Lambrechts et al._ , 2016). Unlike standard SI, the
instability grows rapidly even for the smallest grains at low
$\epsilon_{\mathrm{d}}$, suggesting the possibility of early seeding of grain
growth through clumping before grains reach the midplane. It is also more
robust than SI to polydisperse grains, although the nonlinear clumping it
causes may be rather weak (_Krapp et al._ , 2020).

####  3.2.2 Nonlinear saturation and dust concentration

The discussion in §3.2.1 considers how a system of coupled gas and dust in a
Keplerian disk responds to small perturbations. As the instability drives
these perturbations to grow exponentially with time, they should ultimately
become nonlinear, leading to a saturated, potentially turbulent, state.

_Johansen and Youdin_ (2007) first studied the nonlinear saturation of the
streaming instability and the resulting turbulence with an unstratified local
shearing box and a single dust species. Depending on the dust size, the system
could reach two distinctly different saturation states, which continues to be
phenomenologically accurate for more complex systems. When the dimensionless
stopping time $\tau_{\mathrm{s}}\gtrsim 1$, the dust undergoes traffic jams
and organizes into axisymmetric dust filaments. When
$\tau_{\mathrm{s}}\lesssim 0.1$, by contrast, the system generates numerous
dust-gas vortices vertically, collecting dust in between. The former could
lead to a $\sim$103 enhancement in dust concentration, while the latter a much
weaker $O(10)$ enhancement. The gas remains incompressible to a high degree
($\lesssim$0.1%). These saturated states are qualitatively consistent across
various simulations with different numerical methods, although some issues in
numerical convergence at high-end tail of the particle density distribution
might exist (_Bai and Stone_ , 2010; _Yang and Johansen_ , 2016; _Benítez-
Llambay et al._ , 2019).

As mentioned in §3.2.1, two distinct regimes occur in the linear growth of the
streaming instability with polydisperse dust grains – fast and slow growth –,
and it appears that this dichotomy also carries over to their nonlinear
saturation in unstratified disks (_Yang and Zhu_ , 2021). When
$\epsilon_{\mathrm{d}}\gtrsim 1$ or $\tau_{\mathrm{s,max}}\gtrsim 1$, the
dust-gas dynamics at the saturation state is similar to that driven by
monodisperse dust grains with $\tau_{\mathrm{s}}\sim\tau_{\mathrm{s,max}}$
discussed above (see also _Schaffer et al._ , 2021). On the contrary, when
$\epsilon_{\mathrm{d}}\lesssim 1$ and $\tau_{\mathrm{s,max}}\lesssim 1$, the
system at the saturation state appears to be close to laminar. Interestingly,
_Yang and Zhu_ (2021) found that when the saturation state is turbulent (i.e.,
$\epsilon_{\mathrm{d}}\gtrsim 1$ or $\tau_{\mathrm{s,max}}\gtrsim 1$),
significant dust segregation by size occurs and the mean radial drift of dust
grains of different size is noticeably altered from the drag-force equilibrium
of _Nakagawa et al._ (1986).

When vertical gravity from the central star is considered, the saturated
turbulence driven by the streaming instability sustains a vertically
stratified layer of dust with a finite scale height. The same dust-gas
vortices could be seen near the mid-plane in the stratified simulations of
_Yang et al._ (2017) with single dust species. With multiple species, _Bai and
Stone_ (2010) and _Schaffer et al._ (2018) demonstrated clear scale separation
between different species (see also _Yang and Zhu_ , 2021). If the initial
solid loading is not sufficiently high, this equilibrated stratified dust
layer can be maintained for more than thousands of orbital periods, when
simulations ended.

If the solid loading, quantified by the dust-to-gas _column_ density ratio
$Z\equiv\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$, is above some critical value
$Z_{\mathrm{c}}$, strong concentration of solids occurs (_Johansen et al._ ,
2009). This leads to nearly _axisymmetric_ , dense, narrow filamentary
structures of solids (_Yang and Johansen_ 2014; _Li et al._ 2018; see also
_Flock and Mignone_ 2021), which may have observational consequences
(_Scardoni et al._ , 2021). Without external turbulence, typical radial
separation between adjacent filaments is about 20% of the gas scale height
$H_{\mathrm{g}}$. With external turbulence such as non-ideal MHD, the
separation can be on the order of about $H_{\mathrm{g}}$ (_Yang et al._ ,
2018). In any case, the separation appears to decrease with increasing $Z$
(_Yang et al._ , 2017, 2018). As mentioned in §3.2.1, these dense filaments
are _not_ correlated with local pressure maxima (_Yang and Johansen_ , 2014;
_Li et al._ , 2018) and hence the latter may not be responsible for driving
the concentration of solids in this scenario. Furthermore, the boundary
between strong clumping of solids when $Z\gtrsim Z_{\mathrm{c}}$ and no strong
clumping when $Z\lesssim Z_{\mathrm{c}}$ appears to be sharp (see also below).
However, the relationship of this boundary to the linear streaming instability
discussed in §3.2.1 (e.g., its change at $\epsilon_{\mathrm{d}}\simeq 1$) is
tenuous (see _Li and Youdin_ , 2021, §3.5). It remains to be understood
whether this dichotomy arises because other linear instabilities are important
(e.g., that of _Lin_ 2021 with vertical stratification), or whether clumping
to Roche densities is inherently nonlinear. Further study of this important,
but subtle, issue is needed.

Figure 3: Critical solid abundance $Z_{\mathrm{c}}$ as a function of
dimensionless stopping time $\tau_{\mathrm{s}}$. Each line is the best fit to
the numerical simulations conducted in different work, and above the line may
strong clumping of solids occur and trigger the formation of planetesimals.
All works compiled here assumed a dimensionless radial pressure gradient of
$\Pi=0.05$.

The value of $Z_{\mathrm{c}}$ should depend on several parameters, including
the dimensionless stopping time $\tau_{\mathrm{s}}$ (single or multiple
species), the radial pressure gradient, as well as the dynamics of external
disturbances. The radial pressure gradient is often quantified by the
dimensionless number $\Pi\equiv\Delta u_{\phi}/c_{\mathrm{s}}$ (_Bai and
Stone_ , 2010) or $\eta\equiv\Delta u_{\phi}/v_{\mathrm{K}}$ (_Nakagawa et
al._ , 1986), where $\Delta u_{\phi}$ is the reduction in the equilibrium
azimuthal velocity of the gas , and $c_{\mathrm{s}}$ and $v_{\mathrm{K}}$ are
the local speed of sound and Keplerian speed, respectively. At a given and
representative radial pressure gradient $\Pi=0.05$, without external
disturbances, and with single species, _Carrera et al._ (2015) first measured
$Z_{\mathrm{c}}$ as a function of $\tau_{\mathrm{s}}$ by gradually removing
the gas over the course of the simulations. By fixing $Z$ instead in each
simulation, the measurement for $\tau_{\mathrm{s}}<0.1$ was later modified by
_Yang et al._ (2017), and _Li and Youdin_ (2021) further found significantly
lower threshold for $\tau_{\mathrm{s}}\gtrsim 0.015$. The comparison between
the three works is shown in Fig. 3. To trigger strong solid concentration for
small particles with $\tau_{\mathrm{s}}\lesssim 10^{-2}$, a value of $Z$
greater than a few percent is in general required. The timescale for forming
dense solid filaments when $Z>Z_{\mathrm{c}}$ appears reciprocal to
$\tau_{\mathrm{s}}$, i.e., proportional to the timescales of radial drift and
vertical sedimentation, and in fact a smaller dust species may lead to a much
stronger concentration than a larger one (_Yang et al._ , 2017). When varying
radial pressure gradient, it appears that the stronger the gradient, the more
difficult to drive strong concentration of solids (_Bai and Stone_ , 2010;
_Abod et al._ , 2019), and $Z_{\mathrm{c}}$ may simply scale linearly with
$\Pi$ (_Sekiya and Onishi_ , 2018). Finally, for a system with a dust-size
distribution, it remains unclear how to best quantify the critical condition
(_Bai and Stone_ , 2010; _Schaffer et al._ , 2021).

The evaluation of this critical solid abundance $Z_{\mathrm{c}}$ becomes
complicated when external disturbances to dust-gas dynamics driven by
(magneto-)hydrodynamical instabilities are present (see §2 and §4). For large
particles with $\tau_{\mathrm{s}}\sim 1$, it was found that the MHD turbulence
driven by ideal MRI in general assists local concentration of solids and hence
triggers their strong clumping with a lower $Z_{\mathrm{c}}\simeq 1\%$ (as
compared to $Z_{\mathrm{c}}\simeq 2$–3% with a similar numerical setup but
without MHD turbulence) (_Johansen et al._ , 2007, 2011). However, it appears
that under such an ideal MHD turbulence, strong clumping of particles with
$\tau_{\mathrm{s}}\sim 0.1$ is significantly more difficult to trigger,
potentially making $Z_{\mathrm{c}}\gtrsim 8\%$ (_Yang et al._ , 2018).

the _anisotropic_ turbulence driven by the non-ideal MHD or the VSI might
alleviate this difficulty for planetesimal formation. _Yang et al._ (2018)
found that in a layered accretion disk with an Ohmic dead zone,
$Z_{\mathrm{c}}\simeq 2\%$ for particles of $\tau_{\mathrm{s}}\sim 0.1$.
Similar $Z_{\mathrm{c}}$ was also found in disks dominated by ambipolar
diffusion (_Xu and Bai_ , 2021). _Schäfer et al._ (2020) found that in some
cases, the interaction between the VSI and the streaming instability makes the
clumping of solids even stronger as compared to the latter alone. In both
scenarios, the turbulence near the mid-plane is highly anisotropic, with the
vertical stirring of particles much stronger than the radial one (see also
_Zhu et al._ , 2015; _Stoll and Kley_ , 2016; _Riols and Lesur_ , 2018). This
also raises the question of how sensitively the triggering of strong clumping
depends on vertical sedimentation of solids and hence the local dust-to-gas
density ratio $\epsilon_{\mathrm{d}}$ near the mid-plane.

###  3.3 Secular gravitational instabilities

The idea that gravitational instabilities (GI) in the dust layer cause solids
to collapse into planetesimals has a long history (_Safronov_ , 1969;
_Goldreich and Ward_ , 1973; _Youdin and Shu_ , 2002; reviewed in more detail
by _Chiang and Youdin_ , 2010). Several Solar System observations support the
hypothesis of planetesimal formation by gravitational collapse (_Morbidelli et
al._ , 2009; _Nesvorný et al._ , 2010, 2019; _Blum et al._ , 2017b; _McKinnon
et al._ , 2020) .

The gaseous component of protoplanetary disks is at best marginally
gravitationally unstable, and the dust component has less self-gravitating
mass at early times. If the gas could be ignored, then collisional damping of
particle motions ensures gravitational fragmentation (_Tanga et al._ , 2004;
_Michikoshi et al._ , 2007). Gas cannot be ignored in the early stages of
planet formation, which introduces both stabilizing and destabilizing effects
on a self-gravitating dust layer.

Gas drag can also lead to strong particle clumping by the SI or in pressure
bumps (§3.1 and §3.2.2; see also _Pinilla and Youdin_ , 2017). Strong particle
clumping clearly facilitates GI. The conditions for gravitational collapse in
a dynamically clumping medium is addressed below. Dust gravitational
instabilities in an otherwise smooth disk of gas and dust strong clumping,
e.g. outside pressure bumps and below SI clumping thresholds (Fig. 3; _Li and
Youdin_ , 2021).

Dust that is perfectly coupled to the gas gives a limiting case of GI.
_Sekiya_ (1983) analyzed a midplane layer in this limit and found that GI
requires a total density, $\rho_{0}>0.17\rho_{\rm R}$, where $\rho_{\rm
R}=3.53M_{\ast}/r^{3}$ is the Roche limit below which a hydrostatic fluid body
is tidally disrupted (see _Chiang and Youdin_ , 2010 for discussion). If the
midplane layer were dust-free, then applying this _Sekiya_ (1983) condition to
the midplane of a gas disk (that is vertically isothermal with $\gamma=7/5$)
translates to $Q_{\rm g}=c_{s}\Omega/(\pi G\Sigma_{\rm g})\lesssim 0.25$,
stricter than usual due to finite thickness, the assumed incompressibility and
other details of the sublayer analysis. With dust (more relevantly),
instability occurs for midplane dust-gas ratios $\epsilon_{\mathrm{d}}\gtrsim
4Q_{\rm g}-1$, which is large for gravitationally stable gas disks.

For larger grains that slip relative to gas, a new mode of secular
gravitational instability (SGI) exists (_Ward_ , 1976, 2000; _Youdin_ , 2005;
_Michikoshi et al._ , 2012). Due to the transfer of angular momentum from dust
to gas, SGI modes can have wavelengths longer than the standard Toomre (dust)
upper limit of $\lambda_{\rm T,d}=4\pi^{2}G\Sigma_{\rm d}/\Omega^{2}$. The
possibility of SGI forming wide dust rings has been proposed as an explanation
for observed structures in protoplanetary disks (_Takahashi and Inutsuka_ ,
2016) and — after the rings fragment into many planetesimals — for the radial
zonation of asteroid spectral classes (_Youdin_ , 2011).

However, the viability of the SGI in real protoplanetary disks is unclear.
Radial turbulent diffusion was shown to have a strong stabilizing effect on
SGI, especially for smaller solids (_Youdin_ , 2011; _Shariff and Cuzzi_ ,
2011). _Takahashi and Inutsuka_ (2014) analyzed “two fluid” models that
include drag feedback (similar to the linear SI, but with self-gravity and
without the vertical motions) and showed that sufficiently long wavelength
modes are stabilized. Their stability criterion, valid for
$\alpha_{\mathrm{D}}\ll\tau_{\rm s}\ll 1$, i.e. weak turbulent diffusion and
small solids, is

$\displaystyle\epsilon_{\mathrm{d}}(1+\epsilon_{\mathrm{d}})$
$\displaystyle>Q_{\rm g}^{2}$ (12)

where we use $Q_{\rm g}=Q_{\rm g,sub}H_{\rm p}/H_{\rm g}=Q_{\rm
g,sub}\sqrt{\alpha_{\mathrm{D}}/\tau_{\rm s}}$ to replace the Toomre $Q_{\rm
g,sub}$ of the gas subdisk that appears in the original works (see also
_Latter and Rosca_ , 2017; _Tominaga et al._ , 2020). If we require the gas
disk to be gravitationally stable, $Q_{\rm g}\gtrsim 2$, then the instability
condition becomes $\epsilon_{\mathrm{d}}\gtrsim Q_{\rm g}$, i.e. similar to
the Roche-like criteria of _Sekiya_ (1983).

In general, an effective viscosity that transports momentum is known to
facilitate GI (_Lin and Kratter_ , 2016). _Tominaga et al._ (2019) discovered
a viscous GI in a gas-dust sublayer with an approximate stability criterion

$\displaystyle\frac{\epsilon_{\mathrm{d}}}{\sqrt{1+\epsilon_{\mathrm{d}}}}$
$\displaystyle>Q_{\rm g}\sqrt{\frac{3\alpha}{\tau_{\rm s}}}\,,$ (13)

where
$\alpha\sim\alpha_{\mathrm{D}}\sim\alpha_{\mathrm{S}}\ll\tau_{\mathrm{s}}$ is
assumed. Since $\alpha\ll\tau_{\rm s}$ is again required, this criterion is
easier to satisfy than Equation (12). Simulations with turbulence that could
act as an effective viscosity are thus strongly motivated.

In general, investigations that go beyond height-integrated models of the
isolated dust sublayer are needed to better determine the relevance of the SGI
and related viscous GI for observed disk structures and planetesimal
formation.

###  3.4 Formation of planetesimals and their expected properties

Since PPVI, it has become feasible to include the self-gravity of pebbles in
the simulations of the streaming instability to trigger gravitational collapse
of local pebble concentrations into a large number of bound objects and hence
planetesimals, leading to a statistically meaningful mass function of newborn
planetesimals. At similar conditions, the result appears insensitive to
numerical methods, or the use of sink particles (_Johansen et al._ , 2015) or
clump finding (_Simon et al._ , 2016) to identify the bound objects, although
some variability may exist between different realizations of the same system
(_Rucska and Wadsley_ , 2021).

Excluding the high-mass cutoff, it was found that the differential mass
function could be represented by a power law $\mathrm{d}N/\mathrm{d}M\propto
M^{-q}$ with $q\simeq 1.6$, where $M$ is the mass of a planetesimal.
Therefore, small planetesimals dominate in number, while a few largest
planetesimals dominate the total mass. This function seems not particularly
sensitive to the pebble size ($\tau_{\mathrm{s}}$) or the solid abundance
($Z$) (_Simon et al._ , 2017). On the other hand, the _efficiency_ of
planetesimal formation does appear to depend on the radial pressure gradient
(_Abod et al._ , 2019) and the strength of turbulence (_Gole et al._ 2020; see
also the linear analyses by _Chen and Lin_ 2020 and _Umurhan et al._ 2020;
§3.2.1). The higher the pressure gradient or the stronger the turbulence, the
less efficient the formation of planetesimals is. In fact, the former could be
understood by the more turbulent saturation state of the streaming instability
(_Abod et al._ , 2019), and hence indicating that the region near a pressure
bump may be a favorable location for planetesimal formation (_Carrera et al._
, 2021). Moreover, there may exist a maximum turbulent forcing at
$\alpha_{\mathrm{S}}\sim 10^{-3}$ above which planetesimal formation is
inhibited (_Gole et al._ , 2020).

In contrast to the differential mass function of planetesimals, the cumulative
mass function showed that the largest few planetesimals are situated with a
much steeper slope, indicating a high-mass cutoff. It was convenient to
express this cutoff as an exponential tapering (_Johansen et al._ , 2015),
while it is perhaps more accurate to describe the whole mass function by a
broken power law (_Li et al._ , 2019). In any case, the “knee” of the mass
function defines the characteristic mass $M_{\mathrm{c}}$ of newborn
planetesimals. _Schäfer et al._ (2017) found that $M_{\mathrm{c}}$ correlates
with the separation between axisymmetric dense filaments of pebbles driven by
the streaming instability, and hence that $M_{\mathrm{c}}$ is directly
proportional to the total mass reservoir inside each filament. Moreover,
$M_{\mathrm{c}}$ could also depend on the balance between gravitational
contraction and turbulent diffusion of a collapsing pebble cloud (_Schreiber
and Klahr_ , 2018; _Gerbig et al._ , 2020; _Klahr and Schreiber_ , 2020,
2021). In an attempt to consolidate the simulations in the literature, _Liu et
al._ (2020) fit and gave an expression for $M_{\mathrm{c}}$ as a function of
$Z$, $\Pi$, disk aspect ratio, as well as strength of the disk gravity.

Even though the mass function of planetesimals could be produced and analyzed
with hydrodynamical simulations, these simulations could not resolve the
process of gravitational collapse of the concentrated pebble clouds. To study
the internal structure of a planetesimal resulting from such a process, these
pre-collapse pebble clouds are taken as initial conditions and evolved using
statistical or $N$-body methods. Factors under consideration included
planetesimal mass, initial pebble size or size distribution, pebble porosity,
pebble collision, and gas drag. In general, the timescale for the collapse
decreases with increasing planetesimal mass and is limited by the free-fall
time when the planetesimal radius $R\gtrsim 100$ km (_Wahlberg Jansson and
Johansen_ , 2014). At such a large size $R\gtrsim 100$ km, though, the
collisions between pebbles become fragmenting and much less primordial pebbles
are preserved (_Wahlberg Jansson et al._ , 2017). When $R\lesssim 10$ km, by
contrast, most of the collisions may be bouncing and the resulting
planetesimal could be a rubble pile. Nevertheless, the pebbles are compressed
at bouncing collisions. Icy pebbles may be compacted up to a volume filling
factor of $\sim$20%, while silicate pebbles $\sim$40% (_Lorek et al._ , 2016).
The porosity of a planetesimal then depends on the final size distribution of
pebbles; the more small dust, especially resulting from fragmenting
collisions, the less porous the planetesimal is. Finally, gas drag causes size
sorting of pebbles in a collapsing pebble cloud (_Wahlberg Jansson and
Johansen_ , 2017; _Visser et al._ , 2021). While small pebbles have slow
terminal velocities, large pebbles have long deceleration times. The
combination of both effects along with the nebula conditions determine the
pebble size distribution as a function of radial distance to the planetesimal
center.

Another essential factor is the initial angular momentum of a collapsing
pebble cloud. It appears that almost all pebble clouds in nonlinear
simulations of the streaming instability with self-gravity contained too much
angular momentum and should fragment into planetesimal binaries (_Nesvorný et
al._ , 2019). The distribution of the angular momentum vector showed that
about 80% of the binaries were prograde. By simulating the gravitational
collapse of these pebble clouds using $N$-body methods, these clouds tend to
form binaries with similar-sized objects and appreciable eccentricity
$e\lesssim 0.6$ (_Nesvorný et al._ , 2010; _Robinson et al._ , 2020; _Nesvorný
et al._ , 2021). Clouds with sufficiently low angular momentum may result in
low contact speeds between the two components of the binaries and hence form a
pathway to contact binaries (_Robinson et al._ 2020; _Nesvorný et al._ 2021;
see also _Lyra et al._ 2021). On the other hand, clouds with high mass or
angular momentum may lead to a hierarchical system or multiple binaries. The
orbital inclination of the binary should be similar to that of the initial
pebble cloud, which is limited by the vertical scale height of the pebbles and
hence should be small.

##  4 MAGNETOHYDRODYNAMICAL PROCESSES & WINDS

###  4.1 Framework

####  4.1.1 Ionization

Protoplanetary disks are observed to have both inflows (i.e., accretion) and
outflows (winds). While purely hydrodynamic and thermal processes might be
responsible for both phenomena, various arguments suggest that
magnetohydrodynamic (MHD) processes are important. For example, known
hydrodynamic instabilities that might give rise to a turbulent viscosity are
weak (§2), except gravitational instabilities of disks that are quite massive
and therefore probably quite young (Class 0). Also, accelerating molecular
outflows to speeds of several $\mathrm{km\,s^{-1}}$ is difficult to do purely
thermally without heating the gas to such a temperature that the molecules
dissociate, whereas magnetocentrifugal winds can in principle expel gas
without heating it at all, at least in ideal MHD. MHD, however, requires the
gas to support electrical currents, which in turn requires mobile charge
carriers: free electrons, atomic and molecular ions, and even small charged
grains or polycyclic aromatic hydrocarbons (PAH). In a typical low-mass Class
II system, temperatures are high enough for (partial) thermal ionization only
at small radii ($\lesssim 0.1\textsc{au}$), where only a small fraction of the
mass of the disk resides. The thermally ionized region becomes larger both for
systems with high luminosity (e.g. for a Herbig Star with 55 solar luminosity
this region can extend to 1 au (_Flock et al._ , 2016)), and for systems with
high accretion, as here the MRI heating can heat up and so ionize the disk
further outward (_Mohanty et al._ , 2018). Elsewhere in the disk, as in much
of the interstellar medium (ISM), ionization relies on nonthermal sources.
These include X-ray (_Glassgold et al._ , 2017) and far ultraviolet (_Perez-
Becker and Chiang_ , 2011) photons emitted by the protostar’s corona and
accretion layers, cosmic rays (_Rab et al._ , 2017; _Seifert et al._ , 2021),
radionuclides such as $\mathrm{{}^{26}Al}$ (_Umebayashi and Nakano_ , 2009),
and perhaps strong electric fields associated with MRI (_Inutsuka and Sano_ ,
2005; _Okuzumi et al._ , 2019) or even lightning (_Whipple_ , 1966; _Desch and
Cuzzi_ , 2000). _Güdel_ (2015) gives a pedagogical review of many of these
processes.

The degree of ionization of protoplanetary disks is uncertain because of the
uncertain strengths of the various ionizing processes listed above, but also
because of the uncertain abundance and size distribution of grains, which tend
to trap free charges and promote recombination (_Ivlev et al._ , 2016; _Thi et
al._ , 2019, and references therein). Close to the star where midplane
temperatures are 500-1500 K, however, thermionic emission from grains may
actually increase the abundance of free electrons (_Desch and Turner_ , 2015;
_Jankovic et al._ , 2022); at still hotter temperatures where the grains
sublime, the main source is thermal ionization of alkali metals (K, Na,
$\ldots$), which have low first-ionization potentials.

####  4.1.2 Plasma conductivity

At magnetic field strengths and ionization levels sufficient to drive
accretion at observed rates, the electrical conductivity behaves as a tensor
rather than a scalar (_Wardle and Ng_ , 1999). That is to say, the electric
field $\bm{E^{\prime}}$ in the local rest frame of the gas and the current
density $\bm{J}$ are generally not parallel, because the mobile charges are
deflected by the magnetic field to varying degrees depending on their charge-
to-mass ratios and collision rates with the neutral gas. When the mass density
in the charged species is negligible compared to that of the neutrals, the
tensorial Ohm’s Law $\bm{E^{\prime}}=\bm{\upsigma^{-1}\cdot J}$ reduces to
three terms (_Balbus and Terquem_ , 2001; _Bai_ , 2014)

$\bm{E^{\prime}}=\eta_{\mathrm{O}}\bm{J}+\eta_{\mathrm{H}}\bm{J\times\hat{B}}+\eta_{A}\underbrace{\bm{\hat{B}\times}(\bm{J\times\hat{B}})}_{\bm{J}_{\perp}},$
(14)

in which $\bm{\hat{B}}=\bm{B}/|\bm{B}|$ is a unit vector parallel to the
magnetic field.

The Ohmic diffusivity $\eta_{\mathrm{O}}$ is independent of $|\bm{B}|$ and
therefore dominates when the field is weak compared to the pressure of the gas
and the ionization fraction is very low; these conditions often prevail at the
disk midplane. $\eta_{\mathrm{O}}$ is also the only one of the three
diffusivities that contributes to $\bm{E^{\prime}\cdot B}$, and therefore the
only one that can alter magnetic flux; the others preserve flux but allow it
to drift relative to the neutral gas. Insofar as the charged particles collide
with the neutrals more often than one another, $\eta_{\mathrm{O}}$ is
inversely proportional to the fractional ionization, which typically needs to
be $n_{e}/n_{\textsc{H}}\gtrsim 10^{-13}$ to support MRI in protoplanetary
disks.

The ambipolar diffusivity $\eta_{A}\propto|\bm{B}|^{2}$ dominates the slippage
between the charged and neutral components when the field is strong and the
ionization fraction is well above $10^{-13}$ (but still $\ll 1$), as tends to
be the case at high altitudes in the disk and in their MHD winds. Although it
preserves magnetic flux, ambipolar diffusion contributes to
$\bm{E^{\prime}\cdot J}$, and therefore dissipates magnetic energy (not flux),
which heats the gas.

Finally, the Hall diffusivity $\eta_{\mathrm{H}}$ scales linearly with field
strength. It can be important when the dominant charge carriers of one sign
are magnetized—meaning that their collision rate with the neutrals is less
than their rate of gyration around field lines—while those of the other sign
are not. At inferred field strengths, this tends to be the case at
intermediate altitudes in the disk: i.e. above the midplane, but below the FUV
dissociation front where disk winds may be launched (_Wardle_ , 2007), and
$\eta_{\mathrm{H}}$ may be important even where subdominant (_Wurster_ ,
2021). Unlike the other two diffusivities, the Hall term contributes neither
to $\bm{E^{\prime}\cdot B}$ nor to $\bm{E^{\prime}\cdot J}$ and is therefore
not directly dissipative. But it does contribute to flux migration, yielding
an $\bm{E^{\prime}\times B}$ drift along $-\bm{J}_{\perp}$.

####  4.1.3 Magnetic field

Once the gas is coupled to magnetic fields, the magnetic field strength and
topology become a key control parameters to the dynamics of the system. It is
usual to quantify the field strength with the dimensionless plasma
$\beta\equiv 8\pi P/B^{2}$ parameter. While it is in principle possible to
generate poloidal field loops in situ by dynamo action, such a dynamo remains
difficult to sustain in the regions $R>1\,\mathrm{AU}$ because of the strong
magnetic diffusion (_Simon et al._ , 2013; _Riols et al._ , 2021). Hence, the
large-scale poloidal field threading the disk $B_{\mathrm{p}}$ (that is to say
averaging out the turbulent part) is conserved: i.e it can only be accreted or
diffused away from the central star during the system lifetime. As a result,
it is often the key control parameter of MHD processes in disks, quantified in
terms of $\beta_{\mathrm{p}}\equiv 8\pi P_{\mathrm{mid}}/B_{p}^{2}$, where we
have used the mid-plane pressure $P_{\mathrm{mid}}$. While it is a useful
parameter from a theoretical standpoint, it is also useful to connect it to
physical values. Using a “typical” T-tauri disk with
$\Sigma=200\,\mathrm{g.cm}^{-2}R_{\mathrm{AU}}^{-1}$, one gets $B_{p}\simeq
10\,R_{\mathrm{AU}}^{-11/8}\beta_{\mathrm{p}}^{-1/2}\,\mathrm{G}$ (_Lesur_ ,
2021a, eq. 2.4) . Hence, the typical value expected in disks with
$\beta_{\mathrm{p}}\sim 10^{4}$ is $B_{p}\sim 5\,\mathrm{mG}$ at
10$\,\mathrm{AU}$. If the large-scale field dominates angular momentum
transport, and hence a MHD wind is launched (§4.3.2), it can be shown that the
total field $B$ is tightly linked to the accretion rate. It is then possible
to deduce a minimum field strength for a given accretion rate: $B\gtrsim
0.31\,R_{\mathrm{AU}}^{-5/4}(\dot{M}_{\mathrm{acc}}/10^{-7}M_{\odot}.\mathrm{yr}^{-1})^{1/2}\,\mathrm{G}$
(_Bai and Goodman_ , 2009). Note that this lower limit should be taken with
care given that in most MHD winds, $B$ is dominated by its toroidal component,
which changes sign on both sides of the disk, and hence is relatively
difficult to measure directly from optically thin tracers.

###  4.2 MRI-driven turbulence

The magneto-rotational instability (_Balbus and Hawley_ , 1991) plays an
important role for the accretion processes and overall for the protoplanetary
disk evolution. Since the last 30 years, many works have investigated the MRI
in the context of protoplanetary disks, using theoretical models (_Turner et
al._ , 2014) and lab experiments (_Rüdiger et al._ , 2018). From our current
understanding, we believe that in the inner (T $>$ 1000 K) and outer
protoplanetary disk regions ($\Sigma\leq 15\rm\,gcm^{-2}$), the conditions for
the MRI are met, which are: a weak magnetic field $\beta_{\mathrm{p}}\geq 1$
(_Lesur_ , 2021a, eq. 6.64), shear and high enough ionization (§4.1).

In this section we highlight the most important achievements since the last
PPVI, more specifically, the convergence of the MRI under a net-vertical
field, the MRI with realistic thermodynamics, and the MRI in the regime from
ideal to non-ideal MHD using local and global simulations.

####  4.2.1 Quasi-Ideal regime

In the quasi-ideal MHD regime (that is to say when diffusive length-scales are
much smaller than $H$) and with a net-flux vertical magnetic field present,
the saturation of the MRI is reached by the super-exponential growth of
parasitic instabilities, halting the further growth of the MRI (_Murphy and
Pessah_ , 2015; _Hirai et al._ , 2018). Such conditions of saturated MRI
activity are present in the inner and outer regions of protoplanetary disks,
and many global simulations confirmed the efficiency of the MRI using a net-
vertical magnetic field (_Parkin_ , 2014; _Suzuki and Inutsuka_ , 2014; _Flock
et al._ , 2017a; _Zhu and Stone_ , 2018). Fully saturated MRI turbulence shows
an angular momentum transport coefficient between $\alpha_{\mathrm{S}}=0.01$
and $1.0$ depending on the strength and topology of the initial magnetic
field. For instance, _Salvesen et al._ (2016) found $\alpha_{\mathrm{S}}\simeq
11\beta_{\mathrm{p}}^{-0.53}$ for $10<\beta_{\mathrm{p}}<10^{5}$ in stratified
shearing box simulations threaded by a vertical field, indicating that
$\alpha_{\mathrm{S}}$ is roughly proportional to $B_{p}$. While this result
holds in ”ideal” MHD, in protoplanetary disks we expect the magnetic Prandtl
number $\mathrm{Pm}$, the ratio between the viscous and the magnetic diffusion
rate, to be very small, even in the regions where ideal MHD applies. This
small $\mathrm{Pm}$ regime was investigated using non-stratified high-
resolution simulations, reaching for the first time convergence (_Meheut et
al._ , 2015) with a value of $\alpha_{\mathrm{S}}=3.2\times 10^{-2}$ and a net
horizontal magnetic flux with $\beta=10^{3}$. In contrast, zero-net flux MRI
simulations in the low $\mathrm{Pm}$ limit show turbulent decay
(_Mamatsashvili et al._ , 2020) and no sustained turbulence. Further studies
are needed to understand the characteristics of the MRI at different levels of
$Pm$ and different magnetic configurations.

####  4.2.2 Non-ideal regime

Moving radially outward in the protoplanetary disk, the non-ideal terms of
Ohmic, Hall and Ambipolar diffusion become more important (§4.1), efficiently
reducing the level of MRI turbulence (_Bai_ , 2015; _Simon et al._ , 2018)
down to levels of several $\sim 10^{-4}$ for $\alpha_{\mathrm{S}}$. For these
non-ideal terms, many disk parameters play now an important role, since
magnetic diffusivities (14) depend on the magnetic field strength, local
density of gas and dust and the local high-energetic radiation field of
X-Ray,Cosmic ray and even FUV (_Gole and Simon_ , 2018). These non-ideal MHD
regions are prone to magnetic flux concentrations, triggering zonal flows
which appear as axisymmetric long-lived perturbation and which appear in the
Ohmic (_Bai and Stone_ , 2014), Hall (_Kunz and Lesur_ , 2013) and Ambipolar
regimes (_Simon and Armitage_ , 2014). At the transition between dead and
active zones, these perturbations and flux concentrations are enhanced (_Flock
et al._ , 2015) halting the radial drift of pebbles (_Ruge et al._ , 2016).
The study of the non-ideal MHD evolution of the disk has become more important
in the recent years, especially due to the ability to cause variations in the
surface density. Future studies should extend the parameter regime of the non-
ideal terms and investigate their importance for the magnetic evolution of the
disk.

####  4.2.3 Thermodynamics

In recent years, local and global simulations investigated the effect of
combined radiative transfer and magnetohydrodynamics, capturing the important
heating and cooling processes together with the MRI using radiation MHD
simulations. In optical thin regimes of protoplanetary disks, the temperature
is mostly defined by the irradiation heating, and the MRI strength adapts to
the local profile of temperature and pressure, as shown by (_Flock et al._ ,
2013) using global radiation magnetohydrodynamical simulations of
protoplanetary disks including irradiation. In the optically thick regime
(_Hirose et al._ , 2014; _Scepi et al._ , 2018) the heating by the MRI leads
to an increase of the disk temperature and so to the disk scale height which
again enhances the turbulence and leads to convection. Most of the MHD heating
is taking place in the magnetic reconnection layers (_Ross and Latter_ , 2018)
for which in general high resolution is required (_McNally et al._ , 2014) to
resolve the heating rate. In protoplanetary disks the inner gas disk is
believed to remain mostly optically thin for its own thermal radiation (_Flock
et al._ , 2017a) however a detailed set of gas opacities remain under
investigation. Especially the inner gas disk, between the silicate sublimation
line and the co-rotation radius, depends heavily on the heating and cooling
processes and so the gas opacity includes various molecular lines (_Malygin et
al._ , 2014). At the same time, MRI activity in this region controls the last
stage of the accretion process, transporting matter onto the star while at the
same time feeding the magnetically driven winds in this inner region. Future
studies of the thermodynamics and the MRI should improve the gas opacity to
show the importance of the heating and cooling processes, especially for the
evolution of the inner, MRI-active, gas disk.

####  4.2.4 Outlook

The saturation level, the strength and the appearance of MRI turbulence
dependents on many disk parameters. Most important here are the local degree
of ionization and the large-scale magnetic field $B_{p}$ threading the disk,
as they determine the appearance and the strength of MRI activity, both in the
ideal and non-ideal regimes. Future studies should focus on the non-ideal
terms, the long-term magnetic flux evolution and the proper thermodynamics to
study the importance of the MRI for the disk evolution. Solving for detailed
thermodynamics is crucial not only to obtain realistic levels of turbulence
but also to enable synthetic observations. Heating by the non-ideal MHD terms
(_Mori et al._ , 2019) is significant and should be captured in future
studies.

###  4.3 Outflows

Outflows have been proposed since the 1990s to explain the atomic jets
observed in several YSOs (_Frank et al._ , 2014). The distinction between jets
and winds is relatively unclear in the literature, but it is usually
implicitly linked to the outflow collimation and velocity: while jets are
assumed to be fast (100-1000km/s) and well collimated, and therefore appear
narrow at large distances, winds are usually seen as slow (1-30km/s) conical
outflows. These slow outflows are sometimes called “molecular“ since they are
detected in molecular lines. Here, we make the distinction between two kinds
of outflows: thermally driven outflows (photoevaporation), for which the
energy source is primarily an external radiation field (which could be due to
the central star or the immediate environement), and magnetised outflows,
which are driven by accretion energy.

####  4.3.1 Thermal outflows

Thermal (photoevaporative) winds are covered in detail in the Chapter by
Pascucci et al., where current efforts to detect these winds in gas line
observations are also reviewed. While significant progress has been done on
theoretical modeling of these winds since PPVI , the broad-brush picture of
how thermal winds work and their effect on the surface density evolution of
disks has been largely confirmed. Recent calculations, however, have
highlighted some important discrepancies in the prediction of wind properties,
which will need to be resolved before the models can be used to quantitatively
interpret observations, as detailed in the Chapter by Pascucci et al.

In the simple picture gas can only be thermally unbound, if it is heated to
temperatures that approach the local escape temperature, $T_{\rm
esc}=GM_{*}\bar{m}/rk_{\textsc{b}}$, where $\bar{m}$ is the mean mass per gas
particle; this is $\approx 10^{4}$ K for $M_{*}=M_{\odot}$, $r=1\textsc{au}$,
and $\bar{m}=0.68m_{p}$ (ionized hydrogen at solar abundance). Depending on
the irradiating spectrum, gas is heated to different temperatures. As an
example, Extreme-Ultraviolet (EUV) radiation ($13.6eV<E<100eV$) heats the gas
mainly by photoionization of hydrogen and helium, generally yielding an almost
isothermal gas with temperature around $10^{4}\,\mathrm{K}$. For geometrically
thin, non self-gravitating disks around a 1$M_{\odot}$ star, the radius at
which the escape temperature of the gas equals $10^{4}\,\mathrm{K}$ is
approximately $1\,\mathrm{AU}$. The mass loss profile, in this case, is peaked
around the so-called gravitational radius of the system, which is the
cylindrical radius at which the sound speed of the heated gas equals the
Keplerian orbital speed (_Alexander et al._ , 2014). and has a total rate of
about 10${}^{-10}M_{\odot}/yr$ assuming an EUV flux of 10${}^{41}phot/sec$
(_Alexander et al._ , 2006a, b), as demonstrated by isothermal hydrodynamical
simulations. Soft X-ray radiation ($100eV<E<1keV$) is not efficient at
ionizing hydrogen directly, but it ejects inner-shell electrons from abundant
metals (e.g. Carbon and Oxygen); the ejected photoelectrons have suprathermal
energies and result in secondary ionizations. This yields a more weakly
ionized gas of a range of temperatures from a few thousand to
$10^{4}\,\mathrm{K}$. Two-dimensional radiation-hydrodynamic models yield mass
loss profiles that are more extended with respect to the EUV case and total
mass loss rates of order 10${}^{-8}M_{\odot}/yr$, assuming a soft X-ray flux
of 10${}^{30}erg/sec$ (_Owen et al._ , 2010; _Picogna et al._ , 2019), or even
higher for Carbon-depleted disks (_Wölfer et al._ , 2019). The dominant gas
heating mechanism for FUV photons ($E<13.6eV)$) is photoelectric heating from
dust grains, which is proportional to the grain surface. Thus in the presence
of abundant small grains or PAHs in the disk atmosphere, FUVs can efficiently
heat the gas and may also contribute to drive a thermal wind, as suggested by
hydrostatic 1+1D models (_Gorti et al._ , 2009, e.g.). However, uncertainties
in the atmospheric abundance of PAHs, which are rarely observed in T-Tauri
disks (_Seok and Li_ , 2017), and the variable nature of the FUV flux, which
is linked to the accretion rate onto the central object, make it difficult to
assess the role of FUV-driven winds on the final disk dispersal.

Recent calculations of thermal winds have combined two-dimensional
(axisymmetric) hydrodynamics with explicit, though simplified, treatments of
the chemical and thermal state of the gas under the influence of these
radiative processes (_Wang and Goodman_ , 2017; _Nakatani et al._ , 2018a, b).
These works find that X-rays alone are less effective at driving
photoevaporation than EUV (_Wang and Goodman_ , 2017) or FUV (_Nakatani et
al._ , 2018b), except at very low gas metallicities or if some cooling
processes are suppressed, contrary to _Owen et al._ (2010). However, the
former assumed harder X-ray spectra than the latter: hard ($E\geq 1$ keV)
X-rays penetrate more deeply than softer ones, thereby warming a larger volume
of gas to a lower temperature, possibly $<T_{\rm esc}$, for the same
luminosity. _Wang and Goodman_ (2017) and _Nakatani et al._ (2018b) also
differ among themselves as to the relative importance of EUV and FUV. Clearly,
more work is needed, but it might be helpful to propose a fiducial set of
radiation fields and disk structures against which different groups could test
their codes.

####  4.3.2 MHD outflows

The lack of any systematic process to drive turbulence in protoplanetary disks
has revived the interest in the magnetised outflow (MO) scenario, especially
in the weakly ionized regions, following an idea initially proposed by _Wardle
and Königl_ (1993) and subsequently developed by _Königl et al._ (2010) and
_Salmeron et al._ (2011). These MOs launched from weakly ionized disks
received more support from numerical simulation, first in shearing box models.
_Bai and Stone_ (2013) first found that a magnetically dead disk dominated by
Ohmic and ambipolar diffusion could still emit a MO thanks to the ionized
layer localized at the disk surface. This result was quickly generalized to
the outer, more turbulent, regions ($R>30$ AU, _Simon et al._ 2013) and to
models including the Hall effect (_Lesur et al._ , 2014). Despite their
importance, these early shearing box models suffer from an oversimplification
of the gravitational potential at large $z$ (_Turner et al._ , 2014) which
makes them unreliable to study outflows. Hence, since PPVI, a lot of effort
has been devoted to confirming these results in global models, which do not
suffer from these drawbacks.

The first generation of global simulations including non-ideal MHD effects
(_Gressel et al._ , 2015; _Béthune et al._ , 2017; _Bai_ , 2017) demonstrated
that these outflows were not a mere artifact of the shearing box model, and
could indeed explain accretion rates of the order of
$10^{-8}\,M_{\odot}/\mathrm{yrs}$, assuming a disk surface density
$\Sigma\simeq 10\,\mathrm{g.cm}^{-2}$ at $R=10\mathrm{AU}$ with a disk
magnetization $\beta_{\mathrm{p}}\simeq 10^{5}$. It was also found that the
mass loss rate in the wind
$\dot{M}_{\mathrm{wind}}=2\pi\int\,R\mathrm{d}R\zeta\Sigma\Omega_{K}$ (see eq.
1) was of the order of the mass accretion rate $\dot{M}_{\mathrm{acc}}$. A way
to interpret these high mass ejection rates is to use the locally defined
ejection efficiency

$\displaystyle\xi=\frac{1}{\dot{M}_{\mathrm{acc}}}\frac{\mathrm{d}\dot{M}_{\mathrm{wind}}}{\mathrm{d}\log
R}$

to compare ejection to accretion rates. In steady state, and assuming a
constant $\xi$, it is easy to show that

$\displaystyle\dot{M}_{\mathrm{wind}}=\dot{M}_{\mathrm{acc}}(R_{\mathrm{in}})\left(\left(\frac{R_{\mathrm{out}}}{R_{\mathrm{in}}}\right)^{\xi}-1\right)$

for a wind-emitting disk with inner radius $R_{\mathrm{int}}$ and outer radius
$R_{\mathrm{out}}$ ( these two radii can in principle be arbitrarily chosen in
the wind-emitting region, but we will assume here that they correspond to the
radial extension of the wind-emitting region). It can be shown that if
accretion is the only energy source for the outflow, then energy conservation
requires $\xi<1$ (_Ferreira and Pelletier_ , 1995; _Lesur_ , 2021a, eq.
11.40). This limit is usually exceeded in the first generation of simulations
($\xi\sim 1$ in _Béthune et al._ 2017, $\xi=1.5$ in _Bai_ 2017), indicating
that these outflows cannot be driven solely by accretion energy, and that
another energy source must contribute to the outflow. Since all of these early
solutions included some sort of prescribed coronal heating, it was quickly
realized that these outflows could be of magneto-thermal origin (i.e. the
outflow is a mixture of photo-evaporation and MHD wind).

####  4.3.3 Magneto-thermal winds

The very high mass loading in these simulations triggered the search for more
accurate models of the wind thermodynamics, including some of the processes
involved in photo-evaporation (4.3.1). Using this type of approach, _Wang et
al._ (2019) found ejection indices slightly lower than in the first global
models, with $\xi=0.5-1$, a result also recovered by _Gressel et al._ (2020).
Hence, the thermal driving turned out to be less important than anticipated in
earlier models. _Rodenkirch et al._ (2020) showed solutions of a hot magneto-
thermal wind with the transition to a photoevaporative dominated wind for very
small magnetic fields for $\beta>10^{7}$. Using numerical self-similar
solutions, _Lesur_ (2021b) also found solutions with $\xi\sim 0.5$ without any
coronal heating, confirming that massive outflows with $\xi\lesssim 1$ are
mostly a result of the low magnetizations (high $\beta_{p}$) required to match
the accretion rates expected in YSOs than a thermal effect. Overall, since
$\xi\sim 1$ in all of the MO simulations published to date with
$\beta_{\mathrm{p}}\gg 1$, one expects
$\dot{M}_{\mathrm{wind}}\gtrsim\dot{M}_{\mathrm{acc}}$ if
$R_{\mathrm{out}}/R_{\mathrm{in}}\gtrsim 2$ (i.e. if the wind is launched from
a significant fraction of the disk surface). Hence, very massive outflows
should be the rule, more than the exception, from a theoretical standpoint.

The fact that these outflows are very massive also has an impact on their
kinematics and observational signature. For instance, it is customary to
measure the amount of angular momentum extracted by a MO by its lever arm
$\lambda$ (see eq.2). Theoretically, $\lambda$ can be related to $\xi$ through
$2\xi\simeq 1/(\lambda-1)$ (_Lesur_ , 2021a, eq. 11.40). Hence, the outflows
emitted in magneto-thermal models all have $\lambda\simeq 1.5$. Because
$\lambda$ is a direct measure of the angular momentum ”excess”, it can also be
deduced by measuring the azimuthal and poloidal velocity components of a given
outflow (_Anderson et al._ , 2003) which gives a relatively direct test of
these models on observed jets and molecular winds (ref to I. Pascucci’s
chapter)

Because of these very low $\lambda$ and high $\xi$, these magneto-thermal
outflows are relatively far in parameter space from the magneto-centrifugal
mechanism of _Blandford and Payne_ (1982) from which previous accretion-
ejection solutions have been derived (_Wardle and Königl_ , 1993; _Ferreira_ ,
1997). Indeed, these early solutions all had $\beta_{\mathrm{p}}\simeq 1$,
with $\lambda\gg 1$ and $\xi\ll 1$ and lead to sonic accretion speeds which
makes them incompatible with the known accretion rate and surface density of
PPDs. As pointed out by _Bai et al._ (2016), the new $\beta_{\mathrm{p}}\gg 1$
solutions are closer to the magnetic tower scenario of _Lynden-Bell_ (2003)
since they are mostly driven by the warping of the toroidal magnetic field at
the disk ionized surface. However, it is worth noting that MOs exist at all
$\beta_{\mathrm{p}}$ and make a continuum between these two extreme limits
(_Jacquemin-Ide et al._ , 2019; _Lesur_ , 2021b), hence MHD-driven outflows
should be much more general in accreting systems than thought in the past.

Outlook Finally, if MOs are key in the angular momentum transport of PPDs, it
also implies that our understanding of their secular evolution should be
revised. First and foremost, if the mass accretion rate
$\dot{M}_{\mathrm{acc}}$ is driven by the outflow stress (eq. 2), then
$\dot{M}_{\mathrm{acc}}$ is not proportional to the density gradient as in the
viscous disk model. As a result, such a disk might not spread radially as
efficiently as viscous disks (though see _Yang and Bai_ 2021) nor fill density
gaps spontaneously. Related to this question, there are hints that planet
migration (rate and direction) is probably affected by MOs (_McNally et al._ ,
2020; _Kimmig et al._ , 2020) though these early results have to be confirmed
by models including self-consistent MOs. This will likely have measurable
consequences on planet population synthesis models and therefore predicted
planet populations (see chapter led by S-J Paardekopper).

Since MOs are structures driven by the large-scale field threading the disk,
the strength of this field becomes crucial to predict the evolution of the
system. Indeed, all of the models to date show that the mass loss $\zeta$ and
lever arm $\lambda$ are decreasing functions of $\beta_{\mathrm{p}}$. As a
result, the mass accretion and ejection rates are partially controlled by the
local field strength, which in principle can evolve with time. As a result,
the secular evolution of PPDs might be dictated by the evolution of their
magnetic flux. One of the consequences is that it becomes possible to create
large cavities in a disk as a result of magnetic field concentration in the
inner disk regions, leading for instance to an alternative model of transition
disks (see §5.3).

These prospects are tightly linked to the question of the large-scale magnetic
field strength and transport in these disks, which is still largely unknown
and poorly described. From a theoretical point of view, there is no consensus
on how this large-scale field should evolve. Global simulations in the non-
ideal MHD regime (_Bai and Stone_ , 2017; _Gressel et al._ , 2020; _Lesur_ ,
2021b) suggest that the large-scale field tend to diffuse outwards for
$R\gtrsim 1$ AU, with a velocity increasing with the field strength, but some
theoretical models suggest the field can be transported inwards (_Leung and
Ogilvie_ , 2019). In the ideal MHD regime, which should be valid for
$R\lesssim 1AU$ in PPDs (§4.1), the magnetic field is inversely transported
inwards (_Zhu and Stone_ , 2018; _Jacquemin-Ide et al._ , 2021), at a speed
that also increases with the field strength. Hence, the behavior of disks in
the non-ideal and ideal regimes regarding field transport appears to be
opposite (fig. 4). If this picture is confirmed, it would suggest the
formation of an inner strongly magnetized disk with an outer disk that slowly
loses its initial magnetic field. There is little doubt that such a scenario
would in turn impact planet formation mechanisms (e.g. _Suzuki et al._ 2016),
but it still needs to be confirmed.