(62 / page) 2024

# Dust growth and evolution in protoplanetary disks

Tilman Birnstiel  University Observatory, Faculty of Physics, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, 81679 Munich, Germany;
Exzellenzcluster ORIGINS, Boltzmannstr. 2, D-85748 Garching, Germany;
email<EMAIL_ADDRESS>
###### Abstract

Over the past decade, advancement of observational capabilities, specifically
the Atacama Large Millimeter/submillimeter Array (ALMA) and SPHERE instrument,
alongside theoretical innovations like pebble accretion, have reshaped our
understanding of planet formation and the physics of protoplanetary disks.
Despite this progress, mysteries persist along the winded path of micrometer-
sized dust, from the interstellar medium, through transport and growth in the
protoplanetary disk, to becoming gravitationally bound bodies. This review
outlines our current knowledge of dust evolution in circumstellar disks,
yielding the following insights:
• Theoretical and laboratory studies have accurately predicted the growth of
dust particles to sizes that are susceptible to accumulation through transport
processes like radial drift and settling. • Critical uncertainties in that
process remain the level of turbulence, the threshold collision velocities at
which dust growth stalls, and the evolution of dust porosity. • Symmetric and
asymmetric substructure are widespread. Dust traps appear to be solving
several long-standing issues in planet formation models, and they are
observationally consistent with being sites of active planetesimal formation.
• In some instances, planets have been identified as the causes behind
substructures. This underlines the need to study earlier stages of disks to
understand how planets can form so rapidly.
In the future, better probes of the physical conditions in optically thick
regions, including densities, turbulence strength, kinematics, and particle
properties will be essential for unraveling the physical processes at play.

###### keywords:

planet formation, protoplanetary disks, circumstellar matter, dust, Solar
System, accretion disks

††journal: Annu. Rev. Astron. Astrophys.

###### Contents

  1. 1 INTRODUCTION
    1. 1.1 Setting the stage: planet formation
    2. 1.2 The curtain opens: planet-forming disks observed
  2. 2 The standard disk
  3. 3 DUST DYNAMICS
    1. 3.1 Drag Forces
  4. 4 Drag Force Regimes
    1. 4.1 Systematic Drift
    2. 4.2 Turbulent Mixing
    3. 4.3 Drift / Mixing Equilibrium
    4. 4.4 Origins of pressure traps
  5. 5 DUST COLLISIONAL EVOLUTION
    1. 5.1 Introduction
    2. 5.2 Relative Velocities
  6. 6 Turbulent Relative Velocities
    1. 6.1 Cross-sections
    2. 6.2 Collisional outcomes
      1. 6.2.1 Physics of Dust Collisions
      2. 6.2.2 Sticking
      3. 6.2.3 Bouncing
      4. 6.2.4 Fragmentation, Erosion, Cratering, & Abrasion
    3. 6.3 Condensation/Deposition
    4. 6.4 Numerical methods
      1. 6.4.1 Mass grid based models
      2. 6.4.2 Monte-Carlo based models
      3. 6.4.3 Approximate Models
  7. 7 A GLOBAL PICTURE
    1. 7.1 Time scales
    2. 7.2 Growth limits
    3. 7.3 Particle size distributions
  8. 8 Particle size distributions
    1. 8.1 Global evolution of the dust mass
    2. 8.2 Global evolution of the dust surface density
  9. 9 THE PATH AHEAD
    1. 9.1 Observational Tracers
      1. 9.1.1 Dust opacities and dust spectral properties
      2. 9.1.2 Population properties
    2. 9.2 The growth towards planetesimals

##  1 INTRODUCTION

###  1.1 Setting the stage: planet formation

Today, planetary systems can be studied in great detail, from studies of our
own planet, meteoritics, and other solar system explorations, to exoplanets,
and other planetary systems in formation. It is quite surprising, that despite
this wealth of data, the question of how planets form remains unanswered on a
very fundamental level. It is generally accepted that an accretion disk forms
due to angular momentum conservation as a by-product of star formation. At
typical ISM conditions, around 1% of the mass of that disk is in sub-
micrometer dust particles (Weingartner & Draine, 2001). These particles need
to be brought together to form bodies of kilometers to thousands of kilometers
in size. However, as Youdin & Goodman (2005) write, the problems occur
embarrassingly early: are the first gravitationally bound objects, the
planetesimals, formed by gradual collisional growth or by gravitational
collapse of over-densities, or by a combination of both effects? To understand
the challenges involved and ultimately answer this question, the collisional
evolution of the solid particles and their dynamics and transport need to be
understood. Both of these aspects are intimately linked to each other: the
particle size and composition determine the aerodynamic properties of the
particle. In turn, the transport and dynamics determine collision speeds and
how particles can be locally accumulated. Furthermore, the evolution of the
particles depends on many other unknowns that can fundamentally change the
picture: this ranges from the unknown source of turbulence (or lack thereof)
to the unknown porosity evolution of particles to the unknown mechanism that
turns porous dust aggregates into igneous pieces of rock called chondrules as
they are found in meteorites (Connolly & Jones, 2016).

[] Planetesimalsbuilding blocks of planets, defined as bodies bound together
by their own gravity. Chondrulesmillimeter-sized, spherical inclusions that
make up most of the mass in a family of rocky asteroids called chondritic
meteorites. MRIMagnetorotational Instability. VSIVertical Shear Instability.

The role dust particles play in the formation of planetesimals goes far beyond
providing the material for planet formation. Dust also plays a major role in
determining the structure, evolution, and chemical composition of its parent
accretion disk. Firstly, dust particles are the primary source of continuum
opacity. They determine where starlight is absorbed or scattered and how the
disk is heated, effectively setting the hydrostatic structure of the disk
wherever dust is abundant (see, for example Gorti & Hollenbach, 2009, Woitke
et al., 2009). In the same way, it determines the UV flux inside the disk
which modulates heating, freeze-out and photo-destruction processes of gas-
phase chemical species (Jonkheid et al., 2004, Aikawa & Nomura, 2006).
Secondly, the dust particles provide the surface area on which complex surface
chemistry takes place (e.g., Garrod & Herbst, 2006). As will be described
later in this chapter, settling and radial drift of dust particles can also
act as a conveyor belt for the volatile species that are formed or frozen-out
on the dust surfaces, causing a large-scale redistribution of abundances
within the disk (e.g., Stepinski & Valageas, 1997, Cyr et al., 1998, Ciesla &
Cuzzi, 2006, Krijt et al., 2016, Stammler et al., 2017). Furthermore,
electrons can efficiently be captured on the grain surfaces, which alters the
ionization state of the disk (e.g., Okuzumi, 2009, Ivlev et al., 2016) and
influences the gas dynamics by allowing or preventing the MRI to develop
(Terquem, 2008, Okuzumi & Hirose, 2012, Delage et al., 2022). Several other
hydrodynamic instabilities such as the VSI depend on the local cooling
timescale, which is also determined to a large part by the dust properties
(e.g., Lin & Youdin, 2015, Malygin et al., 2017, Pfeil & Klahr, 2019).

Finally, starlight scattered by small dust in the disk surface as well as the
thermal continuum emission of the dust particles themselves are the most
readily available observational tracers that allow us to observe the disks
within which the planets form, as will be discussed in the following section.
All the points above underline the fact that dust, while only contributing
about 1% to the total mass of the system, is a crucial ingredient affecting
all aspects of planet formation.

###  1.2 The curtain opens: planet-forming disks observed

Figure 1: Observations for the disk around IM Lupi in scattered light
(Avenhaus et al., 2018, panel a), thermal emission (Andrews et al., 2018a,
panel b), and carbon-monoxide 2-1 peak emission (Öberg et al., 2021, Czekala
et al., 2021, panel c), and a schematic of the associated emission mechanisms
(panel d). Circles in panel a are denoting surface features discussed in
Avenhaus et al. (2018) while circles in panel b denote the rings fitted in
Huang et al. (2018) and a
$1.6\text{\,}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}$
circle denoting the outer edge of the continuum emission.

SEDSpectral Energy Distribution

Initial observations of planet-forming disks did not offer high enough
resolution to spatially resolve the circumstellar material and were therefore
limited to studying only the SEDs (for a review, see, for example Natta et
al., 2007). While extremely useful at the time, SED modeling was inherently
limited due to the degenerate nature of the problem: Not every structural
feature in the disk results in a unique spectral counterpart in the SED.
Pioneered by the Hubble Space Telescope (e.g., O’dell et al., 1993, Debes et
al., 2013), high-resolution imaging of planet-forming disks has been
revolutionizing the field. During the last decade, this was foremost driven by
two observatories: firstly, the Atacama Large Millimeter/submillimeter Array
(ALMA), a 66-antenna radio interferometer with baselines of up to 16 km,
achieving unprecedented resolution and sensitivity for observing the long
wavelength emission of disks. Around the time that the ALMA long baselines
became available, also the VLT/SPHERE instrument (Beuzit et al., 2008, 2019)
saw first light, enabling high-contrast imaging of planet-forming disks in
scattered light. Achieving comparable resolution (tens of milliarcseconds),
ALMA and VLT/SPHERE can resolve $1\text{\,}\mathrm{a}\mathrm{u}$-scale
structures in nearby disks (e.g., Andrews et al., 2016, van Boekel et al.,
2017). As an example, Figure 1 shows the same disk around the young star IM
Lup in scattered light (Avenhaus et al., 2018, H-band
$1.625\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, left) as well as dust
continuum ($1.25\text{\,}\mathrm{m}\mathrm{m}$, Andrews et al., 2018a) and gas
line emission (CO 2–1, Öberg et al., 2021, Czekala et al., 2021). In panel d
of Figure 1, the emitting regions of each mechanism are shown schematically:
the scattered light image is detecting stellar radiation that scatters off
small dust particles in the disk surface due to the high optical depth of the
disk at those wavelengths. The structure seen in the image depicts a
vertically extended disk, as shown by the non-concentric ellipses (Avenhaus et
al., 2018). In contrast, the millimeter emission shows a completely different
morphology: the radial extend of most of the emission is much smaller, spiral
features are seen in the inner parts of the disk and the fact that the circles
that trace bright and dark rings (from Huang et al., 2018 and a
$1.6\text{\,}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}$
ring) are concentric, indicates little to no vertical extend of the emission.
Comparing the scattered light image (Figure 1, left panel) to the CO peak
intensity image (right panel), it can be seen that the radial and vertical
extend of the small dust particles and the gas are comparable.

The stark differences between these images point towards significant evolution
of the solid component and can help guide the way towards a better
understanding of the early phases of planet formation. The following section
introduces the dynamics of solid particles in disks which can help explain the
features seen in the observations above. We will see that this will require
the particles to have reached sizes at least a hundred times larger than their
initial ISM size scale. Section 5 will therefore explain the physics of dust
particle growth and the combined effects of growth and transport will be
discussed in Section 7. In the final Section 9, the subsequent growth towards
planetesimals, observational methods to trace these theories, and future
directions will be discussed.

[h]

##  2 The standard disk

Within this review, unless otherwise noted, we will consider a vertically
isothermal disk with a radius-dependent temperature
$T(r)\simeq$200\text{\,}\mathrm{K}$\,(r/$\mathrm{a}\mathrm{u}$)^{-1/2}$. The
disk column density is
$\Sigma_{\mathrm{g}}(r)=\frac{M_{\mathrm{disk}}}{2\pi\,r_{\mathrm{c}}^{2}}\,\frac{r_{\mathrm{c}}}{r}\,\exp\left(-r/r_{\mathrm{c}}\right)$
which results in a vertical gas density profile of
$\rho_{\mathrm{g}}(r,z)=\frac{\Sigma_{\mathrm{g}}(r)}{\sqrt{2\pi}\,h_{\mathrm{g}}}\exp\left(-\frac{z^{2}}{2\,h_{\mathrm{g}}^{2}}\right)$,
where $z$ is the height above the mid-plane,
$h_{\mathrm{g}}=c_{\mathrm{s}}/\Omega_{\mathrm{K}}$ is the gas scale height,
$c_{\mathrm{s}}$ the isothermal sound speed and
$\Omega_{\mathrm{K}}=\sqrt{\frac{\mathrm{G}\,M_{\star}}{r^{3}}}$ the Keplerian
frequency. We will assume a stellar mass of
$M_{\star}=$0.5\text{\,}\mathrm{M}_{\odot}$$, a disk mass of
$M_{\mathrm{disk}}=$0.05\text{\,}\mathrm{M}_{\star}$$, a characteristic radius
of $r_{\mathrm{c}}=$30\text{\,}\mathrm{a}\mathrm{u}$$, and mean molecular
weight of 2.3 proton masses.

##  3 DUST DYNAMICS

The dynamical evolution of dust grains around stars can be affected by many
mechanisms including resonances, radiation pressure, solar wind interaction,
magnetic fields, and others. Inside planet-forming disks, however, small
particles are most strongly affected by their coupling to the gas via drag
forces (but see Owen & Kollmeier, 2019). In the following, the key concepts
will be introduced, and it will be shown how the drag force leads to
systematic motion of dust and how the turbulent gas can act as diffusion. For
other aspects, such as dust dynamics in disk winds, binaries, or warped disks,
the authors may refer to more in-depth literature (e.g. Sellek et al., 2021,
Aly et al., 2021, Zagaria et al., 2023, and others).

###  3.1 Drag Forces

Particles moving with respect to the surrounding gas feel a drag force that
acts to decelerate the particle until the relative velocity vanishes.
Depending on the conditions the drag forces can be linear or non-linearly
proportional to the magnitude of the relative velocity. For most regions in
the disk and particles sizes accessible to observations, the mean free path of
the gas molecules $\lambda_{\mathrm{mfp}}$ is larger than the particle radius
$a$. This ratio $\mathrm{Kn}=\frac{\lambda_{\mathrm{mfp}}}{a}$ is called the
Knudsen number. On these scales, the gas molecules are not in the hydrodynamic
limit. This means, the dust particle experiences the gas not as a flow, but as
a bombardment of individual molecules, where the impact velocities of the
molecules in the direction of motion of the particles are higher by the
relative velocity between the dust and gas. This regime of the drag force is
also called the Epstein regime after Epstein (1924) (see also the box “Drag
Force Regimes”). If the particles become larger, their Reynolds number
$\mathrm{Re}_{\mathrm{p}}$ (see box "Turbulent Relative Velocities") increases
while the Knudsen number decreases and the interaction between the dust
particle and the gas becomes of a hydrodynamic nature, where the Reynolds
number determines the flow structure, as shown in Figure 2.

By relating the momentum of the particle (relative to the gas) to the rate of
momentum change (= the drag force), we can derive the timescale on which the
particle velocity changes,

$t_{\mathrm{stop}}=\left|\frac{m\,v_{\mathrm{rel}}}{F_{\mathrm{D}}}\right|.$
(1)

In the Epstein drag regime, the linear velocity dependence cancels out, and,
if the particle is a sphere of density $\rho_{\mathrm{s}}$, the stopping time
simplifies to

$t_{\mathrm{stop}}=\frac{a\,\rho_{\mathrm{s}}}{\bar{v}_{\mathrm{th}}\,\rho_{\mathrm{g}}}.$
(2)

For conditions at $1\text{\,}\mathrm{a}\mathrm{u}$ (see text box “The standard
disk”), the stopping time is only around a second for a micrometer-sized
particle, but of the order of a month for a meter sized boulder. As conditions
vary greatly throughout the disk, a more useful, dimensionless number is the
Stokes number, $\mathrm{St}=t_{\mathrm{stop}}\,\Omega_{\mathrm{K}}$, which
relates the stopping time and the orbital timescale. Again assuming the
Epstein drag regime, in the mid-plane, the Stokes number becomes

$\mathrm{St}=\frac{a\,\rho_{\mathrm{s}}}{\Sigma_{\mathrm{g}}}\,\frac{\pi}{2},$
(3)

so two particles with the same Stokes number behave aerodynamically identical,
even if their shapes or compositions are different. Eq. 3 shows, that in the
Epstein drag regime, the Stokes number is independent of velocity, and
linearly dependent on the particle radius. This changes in the Stokes drag
regime. Most relevant here is the Stokes drag at low Reynolds numbers, where
the Stokes number is still velocity-independent, but scales as $a^{2}$. For
larger Reynolds numbers, the stopping time is also velocity dependent.

[h]

##  4 Drag Force Regimes

The drag force for a dust particle of size $a$ moving at sub-sonic speed
$v_{\mathrm{rel}}$ relative to the gas can generally be written as

$\vec{F}_{\mathrm{D}}=-\frac{C_{\mathrm{D}}}{2}\,\pi\,a^{2}\,\rho_{\mathrm{g}}\,v_{\mathrm{rel}}\,\vec{}v_{\mathrm{rel}},$
(4)

where $\rho_{\mathrm{g}}$ is the gas density. $C_{\mathrm{D}}$ is the
dimensionless drag coefficient which in the Epstein drag regime
($\mathrm{Kn}\geq 4/9$) is

$C_{\mathrm{D}}^{\mathrm{Epstein}}=\frac{8\,\bar{v}_{\mathrm{th}}}{3\,v_{\mathrm{rel}}}.$

This results in a drag force that is linearly dependent on the relative speed
with respect to the gas. Beyond the Epstein drag, the Stokes drag coefficient
becomes dependent on the particle Reynolds number,
$\mathrm{Re}_{\mathrm{p}}=2\,a\,v_{\mathrm{rel}}/\nu$, where $\nu$ is the gas
viscosity. The classical parameterization of Weidenschilling (1977) for the
drag coefficient is

$C_{\mathrm{D}}\simeq\begin{cases}24\,\mathrm{Re}_{\mathrm{p}}^{-1}&\text{for
}\mathrm{Re}_{\mathrm{p}}<1\\\ 24\,\mathrm{Re}_{\mathrm{p}}^{-0.6}&\text{for
}1<\mathrm{Re}_{\mathrm{p}}<800\\\ 0.44&\text{for
}\mathrm{Re}_{\mathrm{p}}>800\end{cases}$ (5)

The first two, $\mathrm{Re}_{\mathrm{p}}$-dependent cases are often called
Stokes drag, while the last, constant value of $C_{\mathrm{D}}$ at large
$\mathrm{Re}_{\mathrm{p}}$ is sometimes called Newton drag. More recent
empirical fits that better reproduce experimental results are, for example,
given in Cheng (2009) as

$C_{\mathrm{D}}=\frac{24}{\mathrm{Re}_{\mathrm{p}}}\,\left(1+0.27\,\mathrm{Re}_{\mathrm{p}}\right)^{0.43}+0.47\,\left[1-\exp\left(-0.04\,\mathrm{Re}_{\mathrm{p}}^{0.38}\right)\right]$
(6)

and compared in Figure 2.

Figure 2: Drag coefficient and different drag regimes for subsonic conditions,
adapted from Bagheri & Bonadonna (2016). For Knudsen numbers larger than 4/9,
the Epstein regime, the drag coefficient is inversely proportional to the
dust-gas relative velocity which cancels out the velocity dependence in Eq. 4.
Here $M=v_{\mathrm{rel}}/c_{s}$ denotes the Mach number.

###  4.1 Systematic Drift

The dynamics of dust and gas are coupled via the drag forces, as described in
the previous subsection. The amount of momentum lost by a decelerated dust
grain is gained by the gas. Under nominal conditions, only about one percent
of the mass is in solids (Weingartner & Draine, 2001, Lodders, 2003), so the
dust experiences much stronger deceleration or acceleration than the gas. If
the dust is locally enhanced and approaches a dust-to-gas mass ratio of unity,
the dust and gas dynamics affect each other equally, which leads to
complicated dynamics, as discussed later in this chapter.

A dust particle orbiting a star in vacuum will follow a Keplerian orbit. In a
rest-frame that co-rotates with Keplerian velocity, this particle will
oscillate vertically around the mid-plane (if it has non-zero inclination) and
it will oscillate radially (if it has non-zero eccentricity), once per orbit.
If the particle is now embedded in a gas disk that rotates exactly Keplerian
(which is not generally the case as we will see below), the radial and
vertical oscillation of the particle will be damped by the drag forces. The
orbital oscillations will change from a harmonic oscillation to a damped
oscillation. The particle will approach an equilibrium once it has lost its
eccentricity and inclination, so it will effectively have sedimented to a
circular orbit at the mid-plane. At this point, it has reached the gas
velocity and no drag forces remain.

For most particle sizes relevant to this review, the stopping time will be
shorter than the orbital timescale ($\mathrm{St}<1$) which means that
particles approach the gas velocity before they have completed a single orbit.
For the vertical sedimentation, this corresponds to an over-damped
oscillation, similar to a falling feather that quickly approaches its terminal
velocity, where the velocity-dependent drag force becomes equal in magnitude
to the gravitational acceleration. In the case of a sedimenting dust particle,
the force balance requires
$-z\,\Omega_{\mathrm{K}}^{2}-\frac{v_{z}}{t_{\mathrm{stop}}}=0$, where the
first term is the vertical component of the stellar gravitational acceleration
and the second term is the deceleration due to vertical drag forces. This
shows, that the vertical settling speed at height $z$ above the mid-plane is

$v_{\mathrm{sett}}=-z\,\Omega_{\mathrm{K}}\,\mathrm{St},$ (7)

from which the timescale for sedimentation becomes
$t_{\mathrm{sett}}=(\mathrm{St}\,\Omega_{\mathrm{K}})^{-1}$.

Terminal velocityThe velocity reached when all acceleration terms cancel out.

A more general derivation of the dust particle velocities under the influence
of gas drag start from the Euler equations for a gas and a dust fluid of fixed
Stokes number (e.g., Youdin & Goodman, 2005). This approach allows analyzing
the system for instabilities, as discussed in Section 9.2. If we further
ignore advective contributions, as in Nakagawa et al. (1986), we arrive at

$\displaystyle\frac{\partial\vec{v}_{d}}{\partial t}$
$\displaystyle=-\frac{1}{t_{\mathrm{stop}}}(\vec{v}_{d}-\vec{v}_{g})-\frac{\mathrm{G}\,M_{\star}}{r^{3}}\,\vec{r}$
(8) $\displaystyle\frac{\partial\vec{v}_{g}}{\partial t}$
$\displaystyle=-\frac{\epsilon}{t_{\mathrm{stop}}}(\vec{v}_{g}-\vec{v}_{d})-\frac{\mathrm{G}\,M_{\star}}{r^{3}}\,\vec{r}-\frac{\nabla
P}{\rho_{\mathrm{g}}},$ (9)

where $\epsilon=\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$ is the dust-to-gas ratio.
It can be seen that a steady state is reached if

$\vec{v}_{d}-\vec{v}_{g}=\frac{t_{\mathrm{stop}}}{\rho_{\mathrm{g}}(1+\epsilon)}\,\nabla
P.$ (10)

This approximation, called the terminal velocity approximation (Youdin &
Goodman, 2005) shows that dust particles generally drift towards higher
pressure. However, there are exceptions to this rule, for example if both the
gas and the dust are on eccentric orbits. This breaks geostrophic balance and
Eq. 10 is not valid anymore. The lower velocity at apastron leads to a higher
density. Yet, as both gas and dust have the same speed along their orbit, no
drag forces are acting that accelerate or decelerate the dust particles along
the orbit (see, Hsieh & Gu, 2012).

Unlike to what was assumed above, planet-forming disks are pressure supported
with a generally negative pressure gradient. In a dust-free environment, the
force balance between gravity, centrifugal force, and pressure force leads to
a slightly sub-Keplerian gas azimuthal velocity of

$v_{g,\varphi}=v_{\mathrm{K}}\,\sqrt{1-2\,\eta}\simeq(1-\eta)\,v_{\mathrm{K}},$
(11)

where

$\eta=-\frac{1}{2}\left(\frac{h_{\mathrm{g}}}{r}\right)^{2}\,\frac{\partial\ln
P}{\partial\ln r}.$ (12)

For typical conditions, $\eta\sim\text{few}\times
10^{-3}\,\sqrt{r/$\mathrm{a}\mathrm{u}$}$ (steeper in the exponential outer
part of the disk), which means that the gas disk at
$1\text{\,}\mathrm{a}\mathrm{u}$ rotates very close to Keplerian speed, but
the difference of several tens of $\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$ still
causes significant drag on the dust particles.

The components of the gas and dust terminal velocity were derived in Nakagawa
et al. (1986, for earlier works, see ,  and , ) starting from Equations 8 and
9. The resulting deviations from Keplerian velocity in polar coordinates are
($\vec{u}=(v_{r},v_{\varphi}-v_{\mathrm{K}},0)$)

$\displaystyle v_{d,r}$
$\displaystyle=-\frac{2}{\mathrm{St}+\mathrm{St}^{-1}(1+\epsilon)^{2}}\,\eta\,v_{\mathrm{K}}$
(13) $\displaystyle v_{d,\varphi}$
$\displaystyle=-\frac{1+\epsilon}{(1+\epsilon)^{2}+\mathrm{St}^{2}}\,\eta\,v_{\mathrm{K}}$
(14) $\displaystyle v_{g,r}$
$\displaystyle=-\frac{2\,\epsilon}{\mathrm{St}+\mathrm{St}^{-1}\,(1+\epsilon)^{2}}\,\eta\,v_{\mathrm{K}}$
(15) $\displaystyle v_{g,\varphi}$
$\displaystyle=\left(\frac{\epsilon}{1+\epsilon}\cdot\frac{1}{1+\mathrm{St}^{2}\,(1+\epsilon)^{-2}}-1\right)\,\eta\,v_{\mathrm{K}}.$
(16)

It can be seen that all velocities scale in magnitude with
$\eta\,v_{\mathrm{K}}$, the amount that the gas rotates faster or slower than
the Keplerian speed. For low dust-to-gas ratio and $\mathrm{St}\ll 1$, the
dust radial speed becomes

$v_{d,r}=\mathrm{St}\left(\frac{h_{\mathrm{g}}}{r}\right)^{2}\frac{\partial\ln
P}{\partial\ln r}v_{\mathrm{K}}.$ (17)

For the values assumed here (see box “The Standard Disk”), this speed is about
$$-50\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$\simeq$-10\text{\,}\mathrm{a}\mathrm{u}\,\mathrm{k}\mathrm{y}\mathrm{r}^{-1}$$
and higher in the exponential part of the disk. This means that even particles
of Stokes numbers of $10^{-2}$ can move through the entire disk within one
million years.

Equations 13 through 16 also depend on the dust-to-gas ratio
$\epsilon=\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$. For the canonical value of 1%,
this effect is negligible and results in a substantial drift speed of the dust
and no substantial drift of the gas. As $\epsilon$ approaches unity, the
radial drift of the dust will slow down and instead the gas will start to
radially drift outward. Figure 3 shows the dust and gas velocity vectors and
components, relative to Keplerian rotation for a dust-to-gas ratio of 10% (to
amplify the effects) and for different Stokes numbers. The maximum dust radial
drift speed is reached when $\mathrm{St}=1$. The radial drift of the gas is
opposite in direction and lower by a factor of the dust-to-gas ratio,
$v_{g,r}=-\epsilon v_{d,r}$. For small $\epsilon$ and $\mathrm{St}$, both the
gas and dust remains $\eta v_{\mathrm{K}}$ slower than Keplerian. The
necessary drag therefore comes from the fact that the dust has a radial
component. In other words, the dust feels a head-wind that continuously
extracts angular momentum. This headwind is not along the azimuthal direction,
but at an angle of $\arctan(\mathrm{St}/2)$ (for $\epsilon\ll 1$), so mostly
radial for small $\mathrm{St}$ and about 26∘ for $\mathrm{St}=1$. The dust
approaches Keplerian speed only above $\mathrm{St}>1$ and the gas only
approaches Keplerian speed for dust-to-gas ratios above unity.

As discussed later, this classical picture of dust drift is complicated by
several factors. Firstly, not all particles have the same size or Stokes
number. A size distribution therefore needs to be accounted for in above
equations which results in integrals in above equations of motion,
complicating the resulting terminal velocities. See Tanaka et al. (2005),
Okuzumi et al. (2012), Dipierro et al. (2018), or Gárate et al. (2019, 2020)
for details. Secondly, the dust-to-gas ratio can vary vertically if particles
settle to the mid-plane. This will impose also a vertical stratification in
the dust and gas radial speeds and the velocities used in vertically
integrated simulations should account for that (e.g., Gárate et al., 2020). As
an example, gas pressure bumps with a high dust content would be disrupted if
the velocities are computed from the vertically averaged quantities (Taki et
al., 2016). If the vertical stratification is considered, then the impact of
the sedimented dust is even stronger, but it is spatially limited to the mid-
plane (Onishi & Sekiya, 2017) and therefore cannot disrupt the entire vertical
gas column. Similar conclusions were found for the disruption of vortices by
Lyra et al. (2018): two-dimensional considerations would predict disruption of
vortices by dust feedback effects, but in 3D, these effects are limited to the
mid-plane and do not strongly affect the vortex stability.

Figure 3: Dust and gas velocities in equilibrium (Nakagawa et al., 1986),
displayed as deviations from Keplerian rotation normalized by
$\eta\,v_{\mathrm{K}}$. The left panels show the radial (top) and azimuthal
(bottom) velocity components for the gas in green and dust in blue where solid
lines denote a case without dust affecting the gas velocity, while dashed and
dotted denote a dust-to-gas ratio of 0.1 and 1, respectively. The right panel
shows the velocities as vectors for a dust-to-gas ratio of 10%. It shows that
dust feedback is strongest for $\mathrm{St}=1$, leading to a radial outfow of
the gas which remains mostly $\eta\,v_{\mathrm{K}}$ slower than Keplerian.

###  4.2 Turbulent Mixing

The preceding section showed the systematic drift motion that particles
experience in a near-Keplerian disk. However, on top of the ordered near-
Keplerian shear flow, the gas might be turbulent. The coupling of the
particles to the turbulent gas will also cause a random velocity component of
the particles which leads to mixing/diffusion of the particle density and
relative velocities between the particles. The latter is discussed in more
detail in Section 5.2. In the following, we will discuss the effects of
turbulent mixing.

Turbulent mixing of the dust is caused by the turbulent motion of the gas that
is transferring its random motions to a certain amount to the particles (Voelk
et al., 1980). As an example, small particles that are tightly coupled to the
gas quickly adapt to the motion of the gas. As they get picked up by turbulent
eddies, they follow the motion of the eddy until it decays. Larger particles
that might not be as strongly coupled to the gas flow still feel an
acceleration due to their velocity difference with respect to the eddy
velocity and a resulting drag force. In this latter case, a turbulent eddy
induces a random kick to the particle momentum. The velocity dispersion of the
dust can therefore be thought of a superposition of these effects integrated
over all scales of the turbulent spectrum.

In planet-forming disks, the turbulent diffusivity is usually parameterized in
terms of the Schmidt number $\mathrm{Sc}$. However, definitions vary between
describing the ratio of diffusivity to viscosity of the gas
$\nu/D_{\mathrm{g}}$, ratio of gas and dust diffusivity
$D_{\mathrm{d}}/D_{\mathrm{g}}$, or combinations of those, see the discussion
in Youdin & Lithwick (2007). In the following, we will assume that the gas
viscosity and diffusivity are identical (i.e. mass is diffused in the same way
as momentum), and will follow the definition of Youdin & Lithwick (2007),
where the Schmidt number is

$\mathrm{Sc}=D_{\mathrm{g}}/D_{\mathrm{d}}.$ (18)

Within this picture of a turbulent cascade, Dubrulle et al. (1995) derived a
Schmidt number of $\mathrm{Sc}=1+\mathrm{St}$ for vertical mixing. Youdin &
Lithwick (2007) generalized these considerations taking the effects of orbital
motion into account: particles in orbit around a star diffuse differently than
free bodies. The result of angular momentum conservation leads to a reduced
diffusivity of the particles of $\mathrm{Sc}=1+\mathrm{St}^{2}$. An extended
derivation and discussion of this can be found in Binkert (2023).

###  4.3 Drift / Mixing Equilibrium

The vertical evolution of a contaminant (Morfill & Voelk, 1984) under the
influence of vertical settling and diffusion is described by an advection-
diffusion equation,

$\frac{\partial\rho_{\mathrm{d}}(z)}{\partial t}+\frac{\partial}{\partial
z}\left(\rho_{\mathrm{d}}\,v_{\mathrm{sett}}-D\rho_{\mathrm{g}}\frac{\partial}{\partial
z}\left(\frac{\rho_{\mathrm{d}}}{\rho_{\mathrm{g}}}\right)\right)=0.$ (19)

A solution to this equation in steady state (and zero net vertical dust
transport) is found where the diffusive (upward) flux is exactly opposite
equal to the sedimentation (downward) flux. This solution is (see Fromang &
Nelson, 2009)

$\epsilon(z)=\epsilon(z=0)\,\exp\left(\int_{0}^{z}{\frac{v_{\mathrm{sett}}(z^{\prime})}{D(z^{\prime})}}dz^{\prime}\right),$
(20)

where $\epsilon(z)=\rho_{\mathrm{d}}(z)/\rho_{\mathrm{g}}(z)$. Using the
settling velocity (Eq. 7) and a vertical diffusivity of $D=\nu$ (i.e. assuming
$\mathrm{St}\ll 1$), the solution for the vertical distribution of the dust
density becomes

$\displaystyle\rho_{\mathrm{d}}(z)$
$\displaystyle=\rho_{\mathrm{d}}(0)\,\exp\left[-\frac{z^{2}}{2\,h_{\mathrm{g}}^{2}}-\frac{\mathrm{St}_{0}}{\alpha}\left(\exp\left(\frac{z^{2}}{2\,h_{\mathrm{g}}^{2}}\right)-1\right)\right]$
(21)
$\displaystyle\simeq\rho_{\mathrm{d}}(0)\,\exp\left[-\frac{z^{2}}{2\,h_{\mathrm{d}}^{2}}\right]$

where a vertically isothermal gas density profile was assumed, such that the
Stokes number becomes
$\mathrm{St}=\mathrm{St}_{0}\,\exp(z^{2}/(2\,h_{\mathrm{g}}^{2}))$ and the
second line assumes $z\ll h_{\mathrm{g}}$. Eq. 21 shows that the dust
distribution near the mid-plane follows a Gaussian with scale height

$h_{\mathrm{d}}=h_{\mathrm{g}}\sqrt{\frac{\alpha}{\mathrm{St}+\alpha}}$ (22)

but falls off much faster above one gas scale height. A more general
derivation of the transport equation and this equilibrium can be found in
Binkert (2023).

Eq. 21 is the dust distribution in a settling-mixing equilibrium. Similar
results can be obtained for pressure maxima in the azimuthal (Birnstiel et
al., 2013, Lyra & Lin, 2013) and radial direction (Dullemond et al., 2018).

Radial trapping of dust has been discussed in many works dating back to
Whipple (1972), including Paardekooper & Mellema (2004), Rice et al. (2006),
or Pinilla et al. (2012) and many others. It again represents a competition
between the radial advective flux (cf. Eq. 13) and radial diffusion. Radial
gradients in the temperature and density make a fully analytical solution
difficult and, unlike for the vertical structure, a radial net flux through
the pressure bump may exist: a pressure bump that is still receiving dust from
the outside and/or lets some dust pass to the inner disk. However, under the
assumption of a Gaussian pressure profile and no net radial flux, the steady
state surface density again follows a Gaussian with familiar properties

$\Sigma_{\mathrm{d}}(r)\propto\exp\left(-\frac{(r-r_{0})^{2}}{2\,w_{\mathrm{d}}^{2}}\right),$
(23)

where $w_{\mathrm{d}}=w\sqrt{\alpha/(\alpha+\mathrm{Sc}\,\mathrm{St})}$ and
$w$ is the radial standard deviation of the Gaussian gas pressure peak
(Dullemond et al., 2018). For azimuthal pressure peaks, a solution that
balances advective and diffusive fluxes can be found as (Birnstiel et al.,
2013)

$\rho_{\mathrm{d}}(y)\propto\,\rho_{\mathrm{g}}(y)\exp\left(-\frac{\mathrm{St}(y)}{\alpha}\right).$
(24)

Just like the cases of the radial and vertical directions, the solution again
depends on the ratio $\mathrm{St}/\alpha$. This ratio represents the relative
strength of advective motion (=drift) and diffusive transport and is therefore
by definition the Péclet number, apart from factors of order unity. In other
words: only particles with a Stokes number that is larger than $\alpha$ will
experience effective trapping.

The fact that the radial and vertical extent of the dust distribution is
sensitive to $\alpha/\mathrm{St}$ means that this ratio can be constrained by
observations. This has for example been pioneered in Pinte et al. (2016) for
HL Tauri, Dullemond et al. (2018) for the DSHARP sample of Andrews et al.
(2018a). Similar approaches can even constrain asymmetries in radial and
vertical direction (Doi & Kataoka, 2021), or using disk observations at
different wavelengths (Franceschi et al., 2023, Doi & Kataoka, 2023).

###  4.4 Origins of pressure traps

As seen in the previous sections, particles have a general tendency to drift
towards higher pressures and that diffusion is acting against this tendency.
The classical results are that dust sediments to the mid-plane and radially
drifts inward (Nakagawa et al., 1986). For localized pressure maxima, this
means that dust can be trapped, and the peak concentration is only limited by
diffusion. Accumulating large amounts of dust can lead to conditions where
planetesimal formation becomes possible (Youdin & Goodman, 2005, Johansen et
al., 2007). Furthermore, as discussed in Section 1.2, axisymmetric
substructure seems to be ubiquitous in planet-forming disks (Andrews et al.,
2018a, Long et al., 2018). Understanding how and where pressure maxima form is
therefore a key question in planet formation. A wide variety of mechanisms
have been proposed, too many to discuss in depth in this review, which will
only summarize the basics in the following (see Bae et al., 2023 for a
review).

Firstly, variations in the disk viscosity $\nu$ can cause pressure maxima due
to the fact that the viscous transport speed is proportional to the viscosity.
In a steady-state viscous disk,
$\dot{M}_{\mathrm{g}}\propto\Sigma_{\mathrm{g}}\,\nu$. A perturbation to the
viscosity $\nu$ (or likewise the gas pressure) will therefore cause an inverse
effect on the gas surface density $\Sigma_{\mathrm{g}}$. This principle is the
source of many flavors of pressure bump mechanisms: Dead zones (Gammie, 1996)
are regions in which non-ideal effects (mainly Ohmic resistivity) prevents the
MRI from developing. While this simple picture has been complicated by other
non-ideal effects (e.g. due to the Hall effect, see Bai, 2015, or ambipolar
diffusion, see Bai & Stone, 2011, Gressel et al., 2015), the consensus is,
that magneto-hydrodynamic effects will not create the same level of turbulent
viscosity everywhere, and such variation will imprint themselves on the
pressure profile. A related magneto-hydrodynamic effect is the formation of
zonal flows (e.g. Johansen et al., 2009, Bai & Stone, 2014, Béthune et al.,
2016), self-organizing regions of high magnetic pressure, that can be long-
lived enough to trap particles (Dittrich et al., 2013).

The leading hypothesis to form azimuthal pressure traps are vortices. Such
vortices may be caused by the Rossby Wave Instability (Lovelace et al., 1999,
Li et al., 2000, Chang et al., 2023), by the convective overstability /
subcritical baroclinic instability (e.g. Klahr & Hubbard, 2014, Lyra, 2014) or
the vertical shear instability (Nelson et al., 2013, Richard et al., 2016,
Latter & Papaloizou, 2018, Manger et al., 2020). Anticyclonic vortices in
disks have previously been shown to collect dust particles near their center
(e.g. Barge & Sommeria, 1995, Klahr & Henning, 1997). More detailed
calculations for dust trapping in vortices confirmed the amount of azimuthal
concentration in Eq. 24, see Lyra & Lin (2013). As predicted by Wolf & Klahr
(2002), ALMA detected several strongly asymmetric disks (e.g. van der Marel et
al., 2013, Casassus et al., 2013, Andrews et al., 2018a, Dong et al., 2018)
that created lots of interests in vortices as possible sites of planet
formation.

Vortices are not the only possible source of asymmetries: massive companions
can cause the disk to become eccentric (Kley & Dirksen, 2006). This
eccentricity leads to a variation in the orbital speed of the gas with a
minimum at apocenter that causes a local accumulation of gas density without
particle trapping (Hsieh & Gu, 2012, Ataiee et al., 2013). For high companion-
to-star mass ratios above 0.04, this forms a Keplerian-rotating clump with a
density contrast of up to 10 (Shi et al., 2012, Ragusa et al., 2017). This
scenario can not explain asymmetries in disks with multiple rings, such as HD
143006 (Andrews et al., 2018a, Huang et al., 2018, shown in Figure 4), or with
strong density contrasts, as in IRS 48 (van der Marel et al., 2013). van der
Marel et al. (2021) found, that the observed density contrast appears to
increase with an observationally estimated Stokes number, as expected from
azimuthal trapping (Birnstiel et al., 2013). These observations indicate that
asymmetries in the gas may be more abundant than they currently appear, based
on asymmetries in the dust emission and that more might be found when moving
to longer wavelength and comparable or higher resolutions. These findings make
a strong case for azimuthal trapping in vortices, however a small subset of
the sources are consistent with asymmetries caused by massive companions, such
as HD 142527 (Price et al., 2018) where the observed low-mass stellar
companion can account for many of the observed features of the circumbinary
disk.

Figure 4: Dust continuum emission at $1.25\text{\,}\mathrm{m}\mathrm{m}$ of
the disk HD 143006 from the DSHARP Survey (Andrews et al., 2018a). The disk
within $\sim 0.1\mbox{${}^{\prime\prime}$}{}$ appears misaligned (Pérez et
al., 2018). The substructure was proposed to be caused by a planet between the
inner disk and the outer two rings (Pérez et al., 2018, Ballabio et al.,
2021).

Overall, planets are the explanation that currently can explain the majority
of substructure: they cause pressure bumps that can trap dust into symmetric
rings (e.g. Rice et al., 2006, Pinilla et al., 2012, and many others). Those
pressure bumps can form vortices (Lovelace et al., 1999, Li et al., 2000,
Chang et al., 2023) and a single planet can cause multiple rings (e.g. Bae et
al., 2017, Dong et al., 2017), multiple vortices and spiral arms (Lobo Gomes
et al., 2015, Miranda & Rafikov, 2019). A planet can also explain the observed
perturbations in the rotation velocity around the sub-structure (Teague et
al., 2018, Izquierdo et al., 2022) and the velocity kinks in the channel maps
(Pinte et al., 2018). In rare cases, planets have even been imaged (Keppler et
al., 2018, Müller et al., 2018, Benisty et al., 2021). The conclusion is that
planets can be responsible for most observed substructures, yet it does not at
all mean that all sub-structure is caused by planets, nor that each ring
corresponds to one planet. In the future, developments in gas kinematics, and
multi-wavelength dust observations as well as more cases of more direct
evidence for or against planets will have to shed light on these questions.

##  5 DUST COLLISIONAL EVOLUTION

Dust particles in planet forming disks are thought to have initial sizes of
$\lesssim$1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$, as grain growth
in the low density pre-disk environment is expected to be inefficient (Ormel
et al., 2009, Bate, 2022). The following section will introduce the concepts
that determine the collision rates in disks (relative velocities and cross-
sections) and the collision outcome, and then discuss the contribution by
surface growth (condensation) as well as numerical implementations.

###  5.1 Introduction

The size evolution of dust can be thought of a general two-body process where
kinetics and cross-sections determine collision rates, quite similar to a
chemical network. There are however key differences: Firstly, the “reactants”
are not discrete chemical species, but are samples of a continuous
distribution. They can have a range of properties such as mass, composition,
or internal structure. For example, many mass combinations of masses $m_{1}$
and $m_{2}$ can lead to the same resulting mass $m=m_{1}+m_{2}$. Secondly, the
product can be an individual single grain (in the case of perfect sticking),
or it can produce a continuum of results, such as a distribution of fragments.

This means that the mathematical representation will be in the form of an
integro-differential equation, where the change of a quantity in time (how
many particles of a certain type exist) depends on an integral over the size
distribution itself. For the mass dimension alone, this can be expressed in a
generalized version of the Smoluchowski equation (Smoluchowski, 1916),

$\displaystyle\frac{\partial n(m)}{\partial t}=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{m_{1}}K(m,m_{1},m_{2})\,R(m_{1},m_{2})\,n(m_{1})\,n(m_{2})\mathrm{d}m_{1}\mathrm{d}m_{2}$
(25) $\displaystyle-n(m)\int_{0}^{\infty}R(m,m_{1})\,n(m_{1})\mathrm{d}m_{1},$
(26)

where we follow the notation of Stammler & Birnstiel (2022). Here, the change
of the number density distribution $n(m)$ of particles of mass $m$ has a
positive and a negative contribution. The positive contribution is caused by
two particles of masses $m_{1}$ and $m_{2}$ colliding. This collision happens
at a rate $R(m_{1},m_{2})\,n(m_{1})\,n(m_{2})$. How much of the colliding mass
$m_{1}+m_{2}$ is put into particles of mass $m$ is determined by the quantity
$K(m,m_{1},m_{2})$. For perfect sticking, this would be described by a Dirac
delta function, $K(m,m_{1},m_{2})=\delta(m-m_{1}-m_{2})$. Some aspects of
evolving this equation numerically will be discussed in Section 6.4.

The collision rate $R(m_{1},m_{2})$ is the product of the collision cross-
section and the relative velocity of the colliding particles. It is
complicated by two effects: Firstly, one needs to distinguish different
outcomes of the collisions, such as perfect sticking, fragmentation, or
erosion, and chose the adequate distribution $K(m,m_{1},m_{2})$. Secondly, the
probability of each outcome is velocity dependent. A common choice is to
assume an average speed of the particles to compute the collision rate. This
might neglect rare but important outcomes because higher- or lower-than-
average collision speeds are ignored. It was proposed that a highly unlikely
series of events become possible for a very small fraction of the particles,
given that the numbers of particles involved in building planets are extremely
large (around ${10}^{39}\text{\,}\mathrm{m}$icrometer-sized particles are
needed to form an Earth mass). This way, few lucky particles might be able to
continue to grow in an environment where growth would otherwise be hindered
(Windmark et al., 2012b, Garaud et al., 2013). In case mass is the only
particle property of interest and if only mean collision speeds are assumed,
the collision rate factor $R_{C}(m_{1},m_{2})$ of a given collision outcome
$C$ (e.g. sticking) can be simply written as

$R_{C}(m_{1},m_{2})=\sigma_{12}\,\overline{\Delta
v}_{12}\,p_{C}(\overline{\Delta v}_{12}),$ (27)

with the collisional cross-section $\sigma_{12}$, the mean relative velocity
$\overline{\Delta v}_{12}$, and the probability of the collision outcome
$p_{C}$. These terms will be discussed further below. If the distribution of
collision velocities is taken into account, the collision rate factor becomes
an integral over this distribution,

$R_{C}(m_{1},m_{2})=\int_{0}^{\infty}\sigma_{12}\,\Delta v_{12}\,p_{C}(\Delta
v_{12})\,p_{v}(\Delta v_{12},v_{\mathrm{rms}})\,\mathrm{d}\Delta v_{12},$ (28)

where $p_{v}(\Delta v_{12},v_{\mathrm{rms}})$ is the distribution of relative
velocities for a given RMS speed. For the case of a Maxwellian velocity
distribution and simple step-functions of the collision probability, analytic
solutions can be found (e.g. Stammler & Birnstiel, 2022).

###  5.2 Relative Velocities

For particle collisions to occur, particles need to move relative to each
other. Relative motion can be caused by the systematic drift motion of the
particles in radial, vertical, or azimuthal direction. As these are size-
dependent, particles of different sizes move at different velocities, which
results in a relative velocity. Following Stammler & Birnstiel (2022) (see
also earlier works by Weidenschilling, 1984, Dullemond & Dominik, 2005, Tanaka
et al., 2005, Ormel et al., 2007, Brauer et al., 2008a), we can write these
contributions as

$\begin{split}\Delta
v_{ij}^{\mathrm{rad}}&=\left|v_{\mathrm{d,r}}(\mathrm{St}_{i})-v_{\mathrm{d,r}}(\mathrm{St}_{j})\right|,\\\
\Delta
v_{ij}^{\mathrm{azi}}&=\left|v_{\mathrm{d,\varphi}}(\mathrm{St}_{i})-v_{\mathrm{d,\varphi}}(\mathrm{St}_{j})\right|,\\\
\Delta
v_{ij}^{\mathrm{vert}}&=\left|h_{\mathrm{d},i}\min\left(\mathrm{St}_{i},\frac{1}{2}\right)-h_{\mathrm{d},j}\min\left(\mathrm{St}_{j},\frac{1}{2}\right)\right|\,\Omega_{\mathrm{K}},\end{split}$
(29)

for the radial, azimuthal, and vertical relative velocity, respectively, where
$h_{\mathrm{d},i}$ is computed as in Eq. 22, for a Stokes number
$\mathrm{St}_{i}$.

In this picture, equal sized particles (or better, particles of the same
Stokes number) would not experience collisions. This is circumvented by the
fact that the size distribution is never exactly monodisperse, which leads to
a dispersion in relative velocities. Furthermore, there are random
contributions to the particle velocity which give rise to random velocity for
all particles: For very small particles, these are caused by Brownian motion,

$\Delta
v_{ij}^{\mathrm{BM}}=\sqrt{\frac{8\,k_{\mathrm{B}}\,T\,(m_{i}+m_{j})}{\pi\,m_{i}\,m_{j}}}.$
(30)

Larger particles experience random velocities caused by the drag forces of the
turbulent gas. Closed form expressions for turbulent relative motion, based on
the model of Voelk et al. (1980), were introduced by Ormel & Cuzzi (2007) and
are now widely used (however, see Pan et al., 2014 and following works in that
series for more recent results). Figure 5 shows the different contributions to
relative velocities as well as the total relative velocities $\Delta
v_{ij}=\sqrt{\sum_{k}\Delta v_{ij,k}^{2}}$ which is a sum over all
contributions. Several trends are visible:

  * •

All contributions except Brownian motion tend to increase with the particle
size. The maximum is reached at a Stokes number of untity (around
$43\text{\,}\mathrm{c}\mathrm{m}$ for this example). Only azimuthal drift
plateaus at high speeds beyond $\mathrm{St}=1$, other contributions decrease
at these sizes.

  * •

As explained above, only Brownian and turbulent motion have non-zero terms on
the diagonal (i.e. equal sized particle collisions).

  * •

In all cases we find that unequal sized particles (e.g. a large and a smaller
grain) tend to collide faster than equal-sized collision pairs.

  * •

In this example, $\alpha=${10}^{-3}\text{\,}$$. As turbulent velocities scale
with the sound speed, turbulence is the key driver of relative velocities in
the inner disk. In the colder outer disk, radial drift can become the dominant
contribution.

  * •

The transition between Brownian Motion and turbulent relative velocities
happens around $1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$.

Figure 5: Relative velocity contributions in the standard disk at
$1\text{\,}\mathrm{a}\mathrm{u}$ assuming
$\alpha=$5\text{\times}{10}^{-4}\text{\,}$$. The maximum relative velocity is
approximately $62\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$.

[h]

##  6 Turbulent Relative Velocities

Turbulent relative velocities are not as straight-forward to compute as the
relative velocity contributions in Eq. 29 (see Ormel & Cuzzi, 2007, Pan et
al., 2014). While closed form solutions exist (Ormel & Cuzzi, 2007), is can be
helpful and instructive to use approximations that allow order-of-magnitude
estimates (see also Birnstiel et al., 2011, Powell et al., 2019). Following
the approximation of Ormel & Cuzzi (2007), collision speeds between particles
of Stokes number $\mathrm{St}$ and monomers of Stokes number $\mathrm{St}_{0}$
are

$\Delta
v^{\mathrm{mon}}_{0}=u_{\mathrm{rms}}\,\mathrm{Re}^{1/4}(\mathrm{St}-\mathrm{St}_{0})$

if the largest grain has a size below
$a_{12}={2\,\Sigma_{\mathrm{g}}}/(y_{a}\,\pi\,\rho_{\mathrm{s}}\,\sqrt{\mathrm{Re}}),$
where
$\mathrm{Re}=\frac{\alpha\Sigma_{\mathrm{g}}\,\sigma_{H_{2}}}{2\mu\,m_{\mathrm{p}}}$
is the turbulent Reynolds number, $y_{a}=1.6$ and
$u_{\mathrm{rms}}=c_{\mathrm{s}}\sqrt{\frac{3}{2}\alpha}$ is the gas turbulent
velocity. If the larger particle radius $a>=a_{12}$ (intermediate turbulent
regime in Ormel & Cuzzi, 2007),

$\Delta
v^{\mathrm{mon}}_{1}=u_{\mathrm{rms}}\,\sqrt{\frac{3}{\mathrm{St}+1/\mathrm{St}}}$

A good approximation therefore is

$\Delta
v^{\mathrm{mon}}=u_{\mathrm{rms}}\,\min\left(\mathrm{Re}^{1/4}(\mathrm{St}-\mathrm{St}_{0}),\sqrt{\frac{3}{\mathrm{St}+1/\mathrm{St}}}\right)$

while collisions among equal-sized particles can be approximated by

$\Delta v^{\mathrm{eq}}=\sqrt{\frac{2}{3}}\,\begin{cases}0&\text{if
}a<a_{12}\\\ \Delta v^{\mathrm{mon}}&\text{if }a\geq a_{12}.\end{cases}$

###  6.1 Cross-sections

The next ingredient to the collision kernel is the cross-section of the
respective collision. In the simplest case for spherical particles, this
amounts to the geometric cross-section

$\sigma_{12}=\pi\left(a_{1}+a_{2}\right)^{2}.$

For spherical, compact grains (this means uniform porosity or fractal
dimension $f_{\mathrm{D}}=3$), the particle radius is simply related to the
mass as $m=4\pi\rho_{\mathrm{s}}\,m^{3}/3$. The situation gets more complex if
the grains become fractal. This means that the material density varies with
radius such that the fractal dimension $f_{\mathrm{D}}$, defined as $m\propto
a^{f_{\mathrm{D}}}$, is less than 3. In such cases, different definitions are
used, such as the characteristic particle radius
$a_{\mathrm{c}}=\sqrt{\frac{5}{3\,N}\sum_{i=1}^{N}\vec{r}_{i}^{2}}$, where
$\vec{r}_{i}$ is the distance vector from the center of mass of the $N$
monomers that make up the aggregate. This definition is following Mukai et al.
(1992) and was used in Okuzumi et al. (2009) and following works. Earlier
works used the area-equivalent radius which is the radius of a circle with the
same area as the projected area of the aggregate.

Further modifications of the cross-section can stem from particle charging. In
the surfaces of planet-forming disks, dust grains can be positively charged
due to photoelectric charging (Pedersen & Gómez de Castro, 2011, Akimkin,
2015). However, in the denser parts of the disk, below the photon-dominated
region, grains collide more often with electrons than with ions since the
electron velocity is higher than the velocity of the more massive ions. This
causes a net negative charge on the grains (Spitzer, 1941, Okuzumi, 2009,
Ivlev et al., 2016). The grains will therefore feel a repulsive potential and
the collision cross-section needs to be modified with a multiplicative factor,

$f_{\mathrm{ch}}=\begin{cases}1-\frac{E_{\mathrm{el}}}{E_{\mathrm{kin}}}&\text{if
}E_{\mathrm{el}}\leq E_{\mathrm{kin}}\\\ 0&\text{else},\end{cases}$

where $E_{\mathrm{kin}}=\frac{1}{2}m_{ij}\Delta v_{ij}^{2}$ is the kinetic
energy of the collision (with $m_{ij}$ being the reduced mass) and
$E_{\mathrm{el}}=Q_{i}Q_{j}e^{2}/(a_{i}+a_{j})$ is the Coulomb energy barrier
at contact between the two grains of charges $Q_{i}$ and $Q_{j}$ times the
elementary charge $e$.

[] Charging Barriersize limit imposed when electrostatic repulsion overcomes
the kinetic energy of a collision.

Okuzumi (2009) showed that this effect can give rise to a charging barrier
that would prevent small particles from colliding. It was pointed out that
this barrier does not hold for collisions with larger particles that, due to
their larger mass and their larger relative velocity possess much more kinetic
energy. The growth of larger particles can happen in regions where the
charging barrier is not relevant (or affecting only at larger sizes), so
possibly in the disk surface, or the inner and outermost parts of the disk
Okuzumi et al. (2011). Transport processes like vertical or radial drift or
turbulent diffusion could move those seeds into regions where small grains are
prevented from equal size collisions, but could then readily collide with
those larger aggregates.

###  6.2 Collisional outcomes

After having discussed the collision rates and collision speeds, the missing
ingredient at this point is the outcome of the collision, often termed the
collision model. Such a model would need to predict the properties of the
resulting particles, depending on the parameters of the collision and the
colliding particles. With impact parameter, impact speed, particle
composition, particle structure, charges, and more, this is a tremendously
large parameter space. Collision models therefore mostly rely on a limited
subset of this parameter space, where experiments exist and on interpolating
or extrapolating the observed behavior. Experiments can either be done in the
laboratory (see, e.g. Blum & Wurm, 2008) where analogues of cosmic dust
aggregates are produced and collided, mostly under micro-gravity conditions,
or the experiments are of numerical nature, where molecular dynamics codes are
used to simulate aggregate collisions (e.g. Paszun & Dominik, 2009, Wada et
al., 2009, Seizinger & Kley, 2013, Hasegawa et al., 2021, and many others).
Both approaches come with the main drawback that we do not know the sizes and
properties of the real cosmic dust particles well, the experimentally produced
or numerically simulated particles might therefore not represent those well.
Furthermore, the microphysics of the collisions is not fully understood, which
means that the numeric and experimental results do not always agree (Krijt et
al., 2015). This section will summarize the main collisional outcomes and some
aspects of current experimental results (for more details see for example
Blum, 2018 and references therein).

####  6.2.1 Physics of Dust Collisions

The outcome of two colliding dust aggregates is modulated on smallest scales
by the surface between the constituent building blocks of the aggregates (see
Dominik & Tielens, 1997 who build on the elastic theory of Johnson et al.,
1971). Two monomers in contact share a contact area. Deformation of the
aggregate requires that these contracts are either broken, or the monomers are
rolled, twisted, or slid with respect to each other. Dominik & Tielens (1997)
argue that contact breaking and rolling require the least energy and are
therefore the dominant processes. The resulting collisional behavior of
aggregates consisting of many monomers typically involve molecular dynamics
simulations or is based on laboratory work. The outcome (see Figure 6) is
generally divided into growth-positive collisions (sticking, fragmentation
with mass deposition), growth-neutral collisions (aggregates bouncing off each
other, usually with some restructuring involved) and growth negative
collisions (fragmentation or erosion). Wada et al. (2013) find that the
critical collision velocity above which collisions change from growth positive
to negative happens at

$v_{\mathrm{crit}}=C\left(\frac{a_{0}}{$0.1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$}\right)^{-5/6},$
(31)

where the constant $C$ is around $80\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$ for
water ice and $8\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$ for silicates. These
findings are in good agreement with Güttler et al. (2010) and Schräpler & Blum
(2011) who used approximately micrometer-sized silicate monomers. The higher
erosion threshold velocity for ices was experimentally confirmed by Gundlach &
Blum (2015). However, recent laboratory experiments investigated the
temperature dependence of the material constants. Musiolik & Wurm (2019) found
that the critical sticking and rolling forces of ice aggregates drop at low
temperatures to values comparable to silicates. This would imply that the
enhanced stickiness of ices would only apply in a narrow region where the
temperature is high enough, but not as high as to sublimate the ices and that
the pressure and temperature conditions are crucial in determining the
collisional behavior (Gärtner et al., 2017).

Figure 6: Simplified classification of collisional outcomes. Red and green
colors of the collision partners denote net mass loss or gain respectively.

####  6.2.2 Sticking

As outlined above, collisions that happen below the critical velocity can lead
at most to restructuring of the aggregates. Low-speed collisions are therefore
mostly of the category sticking depicted in Figure 6. This outcome can be
further subdivided (see Güttler et al., 2010) into hit-and-stick at lowest
speeds (sticking at first point of contact without relevant structural
changes), or, for increasing impact speed, sticking with surface deformations
or sticking with penetration.

####  6.2.3 Bouncing

Bouncing is a process that happens mainly at impact speeds beyond sticking and
is mainly understood as a collision that does not involve any relevant mass
transfer between the colliding bodies, although subtypes of bouncing with mass
transfer are discussed in Güttler et al. (2010). Bouncing involves some
amounts of inelastic interaction mainly through compaction around the contact
area of the collision. Initially, bouncing was seen in many experimental
studies (e.g. Weidling et al., 2009, Güttler et al., 2010) but could not be
reproduced in simulations for similar parameters (Suyama et al., 2008, Paszun
& Dominik, 2009). Using molecular dynamics simulations, Wada et al. (2011) and
Seizinger & Kley (2013) showed that the exact method of preparation of the
aggregate matters. Their results depend on the coordination number and equally
on the filling factor, while only the latter is easily accessible in
laboratory experiments. They conclude that bouncing of particles below
$100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ is rare, but can happen
for compact particles (filling factor above 50%) and low impact velocities
($<$0.1\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$$).

Coordination numberAverage number of contacts of a monomer to its neighbors in
an aggregate. Filling factor, PorosityThe fraction of the aggregate volume
filled with material, as opposed to vacuum/air. Porosity is the volume
fraction of vacuum/air.

####  6.2.4 Fragmentation, Erosion, Cratering, & Abrasion

Even higher impact velocities tend to lead to mass loss of at least one of the
collision partners. The outcome depends on the mass ratio: mass ratios much
different from unity tend to lead to destruction of the smaller particle, but
often are not energetic enough to fragment the larger particle. Such cases are
typically called erosion as depicted in Figure 6. For less extreme mass
ratios, when more mass is excavated, the process is sometimes called cratering
(Blum, 2018) which appears to gradually transition to complete fragmentation
(right panel in Figure 6). Another effect in same family of outcomes is
abrasion (Kothe, 2016, Blum, 2018) which are small relative mass losses (of
the order of a permill in mass) between similar sized aggregates at low
speeds.

A special case in this family of outcomes is fragmentation with mass transfer,
as observed in experiments of Wurm et al. (2005), Teiser & Wurm (2009), Kothe
et al. (2010), and several following works (see Hasegawa et al., 2023, for a
numerical study). This effect occurs when a small projectile hits a much
larger target at high speed. The impact fragments the projectile, but deposits
some mass of the projectile onto the target, similar to a snowball sticking to
a wall. As impact speeds tend to increase with particle size (see Section
5.2), this provides a potential pathway for continued growth as a larger
particle that gets continuously bombarded with smaller grains can slowly but
continuously grow as long as it does not experience equal-sized collisions.
This possibility will be discussed further in Section 9.2.

###  6.3 Condensation/Deposition

An alternative way of growing dust particles is not by mutual collisions
(often called coagulation in the context of planet formation), but by
deposition from the gas phase to the dust particle surface. This process is
often called condensation111Condensation, strictly speaking, is the transition
from vapor to liquid state. Under disk conditions the liquid phase usually
does not exist, so the vapor is turned to a solid instead. The proper term
therefore should be re- or desublimation, or deposition. However, the term
condensation is still commonly used in the field of planet formation. In the
interstellar medium context, this process if often called accretion.,
deposition or re-sublimation. The process can be modeled as a flux of atoms or
molecules impinging onto the particle surface
$F_{\mathrm{gain}}=\bar{v}_{\mathrm{th}}\,n_{i}/4$ (Hollenbach et al., 2009),
with the vapor number density $n_{i}$. At the same time, this effect is
counteracted by a thermal desorption rate of atoms/molecules leaving the
surface, which is given by the Polanyi-Wigner equation (e.g. Minissale et al.,
2022), $F_{\mathrm{loss}}=-N_{i}\,\nu_{i}\,\exp(-T_{i}/T_{\mathrm{d}})$. Here
$N_{i}$ is the number of desorption sites per area (typically around
${10}^{15}\text{\,}\mathrm{c}\mathrm{m}^{-2}$), and $\nu_{i}$ is the
vibrational frequency of the species $i$ with dimensionless molecular weight
$\mu_{i}$ and a temperature $T_{i}$ that corresponds to the binding energy
$E_{i}=k_{\mathrm{B}}\,T_{i}$. The sum of these fluxes determines whether the
species sublimates or resublimates and an equilibrium is reached if the vapor
pressure reaches the saturation pressure. The desorption rate is complicated
by the detailed microphysics of the substrate and by the surface curvature.
Often, empirically derived values are used, see Minissale et al. (2022) for
details and recommended values of $\nu_{i}$ and $E_{i}$.

The sublimation rate above has an exponential temperature dependence. This
means that sublimation can be very fast once the temperature rises over the
sublimation temperature. The opposite process, resublimation, however, is
limited by the collision with the dust surface area. Typically, small
particles $\leq$1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ dominate the
available dust surface area. The deposition timescale
$t_{\mathrm{resub}}\simeq\frac{\rho_{\mathrm{s}}\,a}{3\,\Sigma_{\mathrm{d}}}\,t_{\mathrm{orb}}$
is then a small fraction of the orbital timescale as long as conditions are
far away from equilibrium.

This already leads to the two main issues of particle growth by condensation:
Firstly, it is a surface effect, meaning equal amounts of mass are deposited
on equal areas. Particle size distributions tend to contain most mass in the
largest particles, however most surface area in the smallest particles (as
long as $q>3$ in $n(a)\propto a^{-q}$). This means that most vapor is
deposited on the smallest grains available instead of significantly growing
the largest particles (although surface curvature might affect those rates as
well). Secondly, while the timescales can be short, this process quickly runs
out of fuel. If vapor is not resupplied, the mass growth stops after a few
$t_{\mathrm{resub}}$. It has been proposed that continuous cycles of
sublimation and redeposition can lead to significant growth (Ros & Johansen,
2013), however this neglected the contribution of small silicate particles.
Stammler et al. (2017) simulated dust growth and fragmentation together with
sublimation and redeposition of carbonmonoxide and found no significant
modifications to the overall size distribution, but the CO ice dominated the
mass for small particles. A preference for homogeneous over heterogeneous
deposition, as proposed by Ros et al. (2019) likely does not change this
picture unless the particle size distribution is very narrow. In any case the
required continuous resupply of vapor limits the relevant effects of
deposition on dust particles to sublimation fronts. Sublimation front /
lineRadial distance from the star where the radial temperature profile crosses
the sublimation temperature of a given volatile species. Also called ice
lines, snow lines, or condensation fronts.

Even if resublimation might not cause strong effects on the particle growth,
the effects of sublimation and resublimation around sublimation fronts, mainly
the water snow line, have a wide range of consequences for planet formation.
If particles become less sticky upon losing water ice, their reduced drift
speed can lead to a traffic jam of solids as shown in Birnstiel et al. (2010)
and Saito & Sirono (2011). The recondensing water vapor can enhance the solid
surface density which can help trigger planetesimal formation (Drążkowska &
Alibert, 2017, Schoonenberg & Ormel, 2017, Lichtenberg et al., 2021, however,
back reactions can hamper these effects, see Gárate et al., 2020).

###  6.4 Numerical methods

Analytic treatment of coagulation and fragmentation is generally limited to
few simple cases that are of limited applicability apart from validating
numerical schemes. They include the constant and linear kernel (Silk &
Takahashi, 1979, Wetherill, 1990) and Brownian motion (Lai et al., 1972). In
some cases a hybrid analytical-numerical scheme can be applied (e.g. Marchand
et al., 2021), but most cases need to be treated numerically. Numeric
solutions to the coagulation equation are challenging, firstly due to
resolution constraints that require large numbers of particles or mass bins
and secondly by the limited precision of floating point numbers. A full
discussion of the numerics of dust evolution is beyond the scope of this
review and only a few references exemplary to the different methods will be
discussed briefly. Numerical dust evolution models in protoplanetary disks and
beyond can generally be classified into the three categories of mass-grid
based, Monte-Carlo based, and approximate models, discussed in the following.

####  6.4.1 Mass grid based models

The particle mass dimension is represented by a grid of size $N$. Dust growth
and fragmentation is then based on computing the source and loss terms. Early
works using such methods include Kovetz & Olund (1969), Weidenschilling
(1980), Nakagawa et al. (1981), Ossenkopf (1993), Lee (2000), or Dullemond &
Dominik (2005). It was shown by Ohtsuki et al. (1990), that realistic growth
can only modeled if more than 7 bins per mass decade are used, otherwise
diffusive effects cause excessive underprediction of the growth time scales
(while some aspects require much higher resolution, see Drążkowska et al.,
2014). Examples of global models that treat dust collisional evolution along
with radial transport are Brauer et al. (2008a), Birnstiel et al. (2010),
Okuzumi et al. (2012) or the open source code dustpy (Stammler & Birnstiel,
2022).

For the case of fragmenting particles, this grid-based approach has
computational costs that scale with $N^{3}$ as there are $N^{2}$ collisions
between the $N$ bins and each collision can affect up to $N$ of the bins
through a distribution of fragments. Rafikov et al. (2020) recently showed
that for self-similar fragment distributions, the costs can be reduced to
$N^{2}$. Attempts to model dust collisional evolution with fewer numbers of
bins were presented in Liu et al. (2019) and Lombart & Laibe (2021) which look
very promising, but have yet to be implemented in multi-dimensional transport
or hydrodynamic codes.

The advantage of grid-based integration of the coagulation equation is that it
can be quite accurate and, unlike Monte-Carlo methods, it is neither subject
to numerical noise nor are some size ranges possibly under-represented.
Implicit integration also allows significant speed-up near steady-state
solutions. Furthermore, the grid-based treatment of dust mass means that there
is no conceptual challenge in implementing it in hydrodynamic grid codes as
they are commonly used in the field. The challenges, however, arise due to the
numerical diffusion and computational costs. First two-dimensional coagulation
simulations coupled to hydrodynamics were shown in Drążkowska et al. (2019).

A major downside of grid-based models is the challenge of adding additional
dimensions beyond the particle mass as this increases both the numeric
complexity and the computational costs. While this was done in few cases (e.g.
Ossenkopf, 1993), recent works used moment approaches to decouple the
properties which results in separate coagulation properties, one for the mass
and another for the addition properties. This was derived and demonstrated for
particle porosity in Okuzumi et al. (2009) and used to treat the
carbonmonoxide abundance in solids in Stammler et al. (2017).

####  6.4.2 Monte-Carlo based models

Monte-Carlo methods for coagulation have the advantage, that the
dimensionality of the model does not increase if additional particle
properties are tracked. In all cases, there is a number of particles $N$ that
are interacting in a randomized way that sample the true collisional
statistics of the model, although larger numbers of particles might be needed
to well sample the larger parameter space. Early Monte-Carlo models were
presented in Gillespie (1975), however, when many small particles are joined
together to form few larger ones, the number of particles is not conserved and
the sampling is reduced. For significant growth, such direct Monte-Carlo
methods are therefore usually not suited and some form particle
spawning/grouping or representative particles need to be used. Examples
include Ormel et al. (2007), Ormel & Spaans (2008), Zsom & Dullemond (2008).
Some models (e.g. Zsom & Dullemond, 2008) sample mass fractions, which usually
means that after substantial growth, no information about the small grain
population is retained, once most mass is in the largest bodies. The grouping
method of Ormel & Spaans (2008) can alleviate this issue, however it comes at
additional computational complexity. Monte-Carlo methods are a straight-
forward way to sample the true statistical variation in the collision outcome,
but they also suffer from the numerical noise from limited sampling
statistics. Drążkowska et al. (2014) compared both grid- and Monte-Carlo-based
methods for particularly challenging situations of low-numbers of particles.
Monte-Carlo based coagulation can also be extended to include particle
transport described by a stochastic equation of motion (Ciesla, 2010, Ormel &
Liu, 2018), as was done in Zsom et al. (2011), Drążkowska et al. (2013) and
others. A coupling to hydrodynamic codes is however challenging as large
numbers of particles are required per grid cell in order to well sample the
dust distribution. Furthermore, the time step of Monte-Carlo based schemes is
bound to the collisional timescale which can be restrictively short. Recently
significant advances of representative particle Monte-Carlo methods stem from
bucketing methods (Beutel & Dullemond, 2023).

####  6.4.3 Approximate Models

Due to the numerical challenges and computational costs of a full treatment of
coagulation, simplified or approximate models have often been preferred. This
often allows a speed-up of many order of magnitude, but comes at the cost of
now knowing if the approximations are valid. Monodisperse growth is one
example, (e.g. Stepinski & Valageas, 1997, Kornet et al., 2001), that often
allows accurate estimates of the overall particle growth timescale (Dullemond
& Dominik, 2005, Birnstiel et al., 2010). More accurate are models that assume
a given size distribution and approximate the coagulation process with moments
(e.g. Garaud, 2007, Estrada & Cuzzi, 2008, Sato et al., 2016), other models
estimate the processes based on the time scales involved (e.g. Ciesla & Cuzzi,
2006, Birnstiel et al., 2012, Vorobyov et al., 2018). Modern approaches
involve neural networks to speed-up computations and/or to make those
approximate treatments more accurate (e.g Pfeil et al., 2022). In the
following section, we will follow the treatment of Birnstiel et al. (2012) due
to the fact that it is simple, yet well reproduces most parts of global dust
evolution models.

##  7 A GLOBAL PICTURE

Having introduced the transport processes in Section 3 and the collisional
evolution in Section 5, this section aims at synthesizing a global picture of
dust evolution. It is important to realize the difference: simulating the
transport of a distribution of small and large grains looks very different
from a simulation that evolves this distribution during transport. On global
scales dust collisional evolution cannot be neglected and even on small
scales, the coagulation timescale can approach the dynamical timescale if the
dust-to-gas ratio is increased (see below). In this section we will discuss
the relevant time scales, the particle size distribution, and how this
knowledge can be used to understand the global behavior of dust in disks on
global scales.

###  7.1 Time scales

To develop an approximate understanding we turn to the timescales involved. In
addition to diffusion and advection time scales, we can derive the collisional
timescale by assuming monodisperse growth such that a particle doubles its
mass on every collision. The resulting growth rate is then
$\dot{a}=\rho_{\mathrm{d}}\Delta v/\rho_{\mathrm{s}}$. The actual growth rate
thus depends linearly on the relative velocity. As shown in Figure 5,
turbulent velocities are often the largest contribution. If mid-plane
conditions, Epstein drag, $\alpha<\mathrm{St}<1$ are assumed and the relative
velocity is approximated with $\Delta
v\sim\sqrt{3\,\alpha\,\mathrm{St}}\,c_{\mathrm{s}}$, one finds
$t_{\mathrm{grow}}\sim 1/(Z\,\Omega_{\mathrm{K}})$. This approximation works
surprisingly well (Birnstiel et al., 2012), however for small particle sizes,
low $\alpha$, or different stellar masses, deviations are to be expected
because other relative velocities apply (e.g. Powell et al., 2019). To
summarize, the important time scales of the system are (with assumptions as in
Eq. 17)

$\displaystyle t_{\mathrm{grow}}$
$\displaystyle\simeq\frac{1}{Z}\,\frac{1}{\Omega_{\mathrm{K}}}$ (growth
timescale) (32) $\displaystyle t_{\mathrm{drift}}$
$\displaystyle\simeq\frac{1}{\mathrm{St}\,\gamma}\frac{\Delta
x}{r}\,\left(\frac{h_{\mathrm{g}}}{r}\right)^{-2}\frac{1}{\Omega_{\mathrm{K}}}$
(drift timescale) (33) $\displaystyle t_{\mathrm{diff}}$
$\displaystyle\simeq\frac{1}{\alpha}\,\left(\frac{\Delta
x}{h_{\mathrm{g}}}\right)^{2}\,\frac{1}{\Omega_{\mathrm{K}}}$ (diffusion
timescale) (34) $\displaystyle t_{\mathrm{sett}}$
$\displaystyle\simeq\frac{1}{\mathrm{St}\,\Omega_{\mathrm{K}}}$ (settling
timescale) (35)

where the magnitude of the logarithmic pressure gradient is
$\gamma=\left|\partial\ln P/\partial\ln r\right|$ which is 2.75 in the
standard disk (but larger in the exponential part of the gas density). The
timescales correspond to e-folding of the particle size (Eq. 32), drifting a
length $\Delta x$ (Eq. 33), diffusion over length scale $\Delta x$ (Eq. 34),
and sedimentation towards mid-plane (Eq. 35).

###  7.2 Growth limits

As discussed in Section 5, dust collisions tend to become growth-neutral or
even destructive as the impact speed increases. Section 3 showed that the
collision speeds tend to increase with the particle sizes. As a result, we can
expect a maximum size that particles can reach before they stop growing due to
bouncing, erosion, or fragmentation. A widely used case is the turbulent
fragmentation barrier where the fragmentation threshold velocity
$v_{\mathrm{frag}}$ is equated with an approximate turbulent collision speed
which results in a Stokes number of

$\mathrm{St}_{\mathrm{frag}}=\frac{1}{3\,\alpha}\left(\frac{v_{\mathrm{frag}}}{c_{\mathrm{s}}}\right)^{2},$
(36)

for Epstein drag and mid-plane conditions, the corresponding particle size is

$a_{\mathrm{frag}}=\frac{2\,\Sigma_{\mathrm{g}}}{3\pi\,\rho_{\mathrm{s}}\,\alpha}\left(\frac{v_{\mathrm{frag}}}{c_{\mathrm{s}}}\right)^{2}\propto\frac{v_{\mathrm{frag}}^{2}}{\alpha\,T}.$
(37)

For very low values of the turbulence parameter $\alpha$, the maximum
turbulent relative velocity $\sim\sqrt{3\alpha/2}\,c_{\mathrm{s}}$ can be too
small to cause fragmentation. Particles then potentially grow to sizes where
the radial drift speed becomes dominant. A drift-induced fragmentation barrier
can be derived in the same way using the relative drift speed (see Birnstiel
et al., 2012),

$\mathrm{St}_{\mathrm{df}}=\frac{2\,v_{\mathrm{frag}}\,v_{\mathrm{K}}}{\gamma\,c_{\mathrm{s}}^{2}},$
(38)

unless the maximum drift speed is also too slow to cause fragmentation. In
that case, the maximum impact speed is reached by azimuthal drift (see Figure
5), where particles with $\mathrm{St}\gg 1$ experience small dust impacting at
the full sub-Keplerian speed $\eta\,v_{\mathrm{K}}$, typically several tens of
meters per second and likely causing erosion (Krijt et al., 2015). It is worth
mentioning that in a pressure maximum, all systematic drift speeds are zero.
In this case, only Brownian motion and turbulence remain as source of relative
velocity.

For bouncing, a size limit can be derived analogously to Eq. 37, however, the
bouncing threshold velocity was found to depend on the particle mass itself.
For threshold velocities reported by Weidling et al. (2012) the bouncing
barrier has a very weak radial dependency (see for example Stammler et al.,
2023, Figure 1). In the simulations shown in Figure 7, it would be
approximately
$$2\text{\times}{10}^{-3}\text{\,}\mathrm{c}\mathrm{m}$\,(r/$10\text{\,}\mathrm{a}\mathrm{u}$)^{-3/16}$.

The charging barrier (see Okuzumi et al., 2011, and Section 6.1) applies to
sizes of few micrometers at most, but mixing and drift is thought to alleviate
the problem. Furthermore, the particle size measurements (see Section 9.1)
appear to show sizes well beyond that barrier, indicating that charging might
at most stall, but not prevent particle growth.

The barriers above are all growth barriers, meaning that they halt particle
growth or make it inefficient. A further particle size barrier that is
limiting particle sizes without directly limiting the growth process is the
radial drift barrier (Birnstiel et al., 2012). It works through transporting
larger particles away, faster than particle growth can resupply them. This is
not unlike hail stones that cannot remain lifted in the cloud once they reach
a certain size. In protoplanetary disks, this limit can be approximated by
equating the drift and the growth timescale: with the growth timescale being
approximately constant with size this finds the particle size that is removed
through drift more efficiently than growth can increase its size. This limit,
expressed in terms of Stokes number and in terms of particle sizes, with
$Z=\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$ is

$\displaystyle\mathrm{St}_{\mathrm{dr}}$
$\displaystyle=\frac{Z}{\gamma}\,\left(\frac{h_{\mathrm{g}}}{r}\right)^{-2}$
(39) $\displaystyle a_{\mathrm{dr}}$
$\displaystyle=\frac{2\,\Sigma_{\mathrm{d}}}{\pi\,\rho_{\mathrm{s}}\,\gamma}\left(\frac{h_{\mathrm{g}}}{r}\right)^{-2}.$
(40)

Figure 7: Snapshots of the dust particle size distribution throughout the disk
at $\sim$${10}^{6}\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$
(left panels) and slices of these distributions at 3, 30, and
$100\text{\,}\mathrm{a}\mathrm{u}$ (right panels). The lower simulation is
identical to the upper simulation but includes a deep gap, as expected for a
Jupiter-mass planet at $20\text{\,}\mathrm{a}\mathrm{u}$. Arrows in the left
panels visualize the growth and drift speeds while the red dashed lines
indicate the trajectories of a particle under the influences of radial drift
and monodisperse growth, where numbers mark the time the particle would need
to reach those points, for example
${10}^{5}\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$ are
indicated with the number 5. Purple lines show the size limits imposed by
radial drift (solid) and fragmentation (dashed). Gray contours indicate the
Stokes number of particles at a given size and radius. The quantity
$\sigma_{\mathrm{d}}(r,a)$ is the dust surface density per decade in size, as
discussed in the Box of Section 8.

###  7.3 Particle size distributions

Another crucial piece of information for any accurate approximation is the
particle size distribution. Depending on the physical process, different
aspects of the distribution might matter. For chemistry, or vapor deposition
the surface area matters, while for the evolution of the dust surface density,
knowing the particle sizes that cause the highest dust flux is important. Dust
particles in the interstellar medium and also particles in collisional
cascades such as in debris disks tend to a size distribution of $n(a)\propto
a^{-3.5}$ (Dohnanyi, 1969, Mathis et al., 1977, Tanaka et al., 1996). In
planet-forming disks, however, not only fragmentation is at play, but also
coagulation. The case of an equilibrium where particles coagulate up to a
maximum size at which they break up (fragmentation-coagulation equilibrium)
was discussed in Birnstiel et al. (2011) where it was found that the
equilibrium size distribution can depend on the size distribution of
fragments, unlike in a classical collisional cascade (but see also Gáspár et
al., 2012, Pan & Schlichting, 2012). The resulting steady-state size
distributions in fragmentation-coagulation equilibrium can well be understood
as a series of analytical approximations, where the exponent of each interval
of the size distribution depends on the dominant source of relative velocities
in that interval.

Figure 7 shows size distributions at different positions in dust evolution
simulations carried out with dustpy (Stammler & Birnstiel, 2022). These
simulations have a simple collisional model with a sharp threshold velocity
that divides sticking and fragmentation (but consider a distribution of impact
velocities). If a bouncing barrier were included, this would typically lie
between the sticking and the fragmentation regime. Bouncing tends to stall the
growth of the largest grains, while small particles might still stick to the
largest grains. In such cases, the size distribution can become almost mono-
disperse (Zsom et al., 2010, Windmark et al., 2012a, b, Stammler et al.,
2023).

The distribution in Figure 7 can be divided into two separate types: the
fragmentation-limited and the drift-limited distributions. The panels on the
left show the size distribution at $1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$
overlaid with the fragmentation and drift size limits. In both cases the mass
is dominated by the largest particles. The size exponent is not constant over
all sizes, but the fragmentation-limited distributions (the
$3\text{\,}\mathrm{a}\mathrm{u}$ case up to
$0.1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ in both simulations and the
$30\text{\,}\mathrm{a}\mathrm{u}$ case in the lower panels at all times) have
a size exponent close to -3.5.

The fragmentation barrier $a_{\mathrm{frag}}$ depends mostly on variables that
are constant or weakly time-dependent and so is the barrier itself. The drift
limit, however, has a linear dependence on the dust surface density
$\Sigma_{\mathrm{d}}$. As drift moves particles inward, the dust surface
density decreases with time and therefore also the size limit.

Dust size distributions that are limited by radial drift tend to be more ‘top-
heavy’, where the upper part of the distribution often follows an exponent of
-2.5 (see Figure 7). This comes about from the fact that particles are only
sticking and growing. Hence, small dust particles are swept-up by the bigger
ones and, without replenishment, get depleted. The shape of the top-end of the
distribution depends on the level of turbulence: while drift and growth tend
to push particles to the drift limit, radial diffusion mixes particles from
radially outside the drift limit (i.e. usually smaller sizes) to the inside
and, vice versa, particles from inside the drift limit (usually larger
particles) outside, hence causing a dispersion in particle sizes. For low
levels of turbulence the dispersion can become very narrow. The size
distribution in the drift limit for small grains is mostly set by the small
grains that are diffused outward from the inner regions of the disk, where
fragmentation retains abundant amounts of small fragments (Birnstiel et al.,
2015).

[h]

##  8 Particle size distributions

Particle size distributions are often stated in terms of $n(a)\propto a^{-q}$
(with varying sign conventions of the exponent), where a canonical choice is
$q=3.5$ as found in the ISM (Mathis et al., 1977) or debris disks and
fragmentation cascades (Dohnanyi, 1969, Tanaka et al., 1996). In this
definition, $n(a)$ describes how many particles per infinitesimal size-
interval $\mathrm{d}a$ exist per unit of volume. The units are therefore
$\text{\,}\mathrm{c}\mathrm{m}^{-3}\,\mathrm{c}\mathrm{m}^{-1}$. In the
context of protoplanetary disk and planet formation, the available dust mass
(not number of particles) is more relevant and due to the disk geometry, the
column or surface density is the density of choice. Vertical integration gives
a column-number density $N(a)=\int_{-\infty}^{\infty}n(a,z)\mathrm{d}z$, where
the vertical dependency of the density was discussed in Section 4.3. The total
dust surface density can then be computed as
$\Sigma_{\mathrm{d}}=\int_{0}^{\infty}m(a)\,N(a)\,\mathrm{d}a$. In case the
particle size distribution is of interest, not just the total dust surface
density, one could define a size distribution in units of surface density
$m(a)\,N(a)$, however this would still not well represent the fact that masses
and sizes of larger particles are often many orders of magnitude larger than
the smallest grains and thus contribute more to the mass integral. It is
therefore common practice to define the quantity on a logarithmic scale,
$\sigma_{\mathrm{d}}(a)=N(a)\,m(a)\,a,$ where the mass integral is understood
as being on the logarithmic axis,
$\Sigma_{\mathrm{d}}=\int_{0}^{\infty}\sigma_{\mathrm{d}}(a)\,\mathrm{dln}a$.
This has the same units as a surface density,
$\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-2}$. This quantity, displayed in
Figure 7, therefore directly visualizes on a logarithmic plot, where most of
the mass is located, as it corresponds to the mass per decade in size.

###  8.1 Global evolution of the dust mass

Equipped with information on the time scales, the sizes particles can reach,
and their distributions, we can draw a global picture of dust evolution. The
growth and drift time scales are visualized in Figure 7 by the arrows, where
short arrows indicate long time scales. The growth timescale (Eq. 32) shows
that particles in the inner disk grow the fastest. Particle growth takes only
few decades to proceed. The inward speed of small particles is very slow,
except at large radii, where the gas density drops off exponentially. Radial
drift is negligible up to a size where the Stokes number (depicted as gray
contours in Figure 7) exceeds $\alpha$ (within factors of a few, comparing gas
and dust radial speeds). Above those sizes, particles start to move
significantly relative to the gas. The fact that larger particles drift faster
can be seen by the longer inward-pointing arrows in Figure 7. This is further
visualized by the red dashed lines in Figure 7 that depict the trajectories of
mono-disperse growth, starting at the lower end of the size axis at the time
of the snapshot. The numbers denote the time these trajectories take to reach
those positions (e.g. 5 denotes
${10}^{5}\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$). The
trajectories show clearly: 1) how dust growth moves initially vertically (to
larger sizes) turns inward (to smaller radii) when the growth and drift time
scales become equal at the drift limit and 2) how particle growth and drift
happen much faster in the inner disk.

Therefore, radial drift acts akin to an inside-out collapse (Shu, 1977): dust
at small radii grow quickly, but they also drift inward quickly. On longer and
longer time scales, regions further and further out have reached sizes where
grains start to drift and supply themselves towards the inner regions.

In the inner parts of the disk, collision speeds are often so high that
fragmentation is limiting particle growth. The speed at which the largest
grains move is therefore the radial drift speed of particles at the
fragmentation size. This speed is not negligible, but considerable smaller
than the maximum drift speed at $\mathrm{St}=1$. The drift speed of particles
(Eq. 17) at the drift limit (for $\mathrm{St}\ll 1$, $Z\ll 1$) can be written
as

$v_{r}(a_{\mathrm{drift}})=Z\,v_{\mathrm{K}},$ (41)

which shows that drift quickly removes the first order of magnitude of dust
mass. The next order of magnitude in mass is then removed on a 10 times longer
timescale (determined by the growth timescale). Globally, the dust-to-gas
ratio should therefore be set by the growth timescale of the outermost parts
of the disk that still contain a significant dust mass, the leaky reservoir as
termed by Garaud (2007). In those regions, the growth and drift time scales
become comparable to the age of the disk after about one or two orders of
magnitude reduction of the dust-to-gas ratio. Where exactly this reservoir is
located is less clear, as the cold outer parts of the disk may be of low
surface brightness but could still contain relevant dust masses (e.g. Ilee et
al., 2022). As time progresses there will be a smaller and smaller total mass
of dust at lower and lower dust-to-gas ratio with longer and longer
evolutionary times that supply a smaller and smaller amount of dust flux to
the disk inward of it. An outer limit to the dust distribution is imprinted
early in the regions where the gas density drops off significantly. There,
already the ISM-sized, (sub-)$\mu$m grains would be rapidly drifting towards
higher densities. This initial phase of drift can form a sharp outer edge in
the dust-to-gas ratio within the first
$\sim$0.1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$$ outside
$\sim$100\text{\,}\mathrm{a}\mathrm{u}$$ to a few hundred
$\text{\,}\mathrm{a}\mathrm{u}$ (Birnstiel & Andrews, 2014).

Figure 8: Evolution of the total dust mass (solid) and gas mass (dashed) of
the two simulations from Figure 7. Orange lines denote the simulation without
a pressure trap, yellow lines the simulation with substructure, but no other
effects that might reduce the observable dust mass, such as planetesimal
formation or pebble accretion. The vertical line shows an appropriate time
that particles in the outer disk require to grow and start drifting. Drift-
limited growth tends to a time dependency of $\propto t^{-1}$. On timescales
longer than the viscous timescale of the disks
($\sim$5\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$$ in this simulation), viscous
effects should steepen it towards $\propto t^{-1.5}$.

Hence, we can expect that dust in the inner regions of a disk ($r\ll
r_{\mathrm{c}}$), if unaffected by substructure, will grow and drain quickly.
These regions do contain only a negligible fraction of the total dust mass.
Regions further outside contain most of the mass, but their growth lags
behind, so during this initial growth phase, we expect an almost constant dust
mass in the disk, as seen in Figure 8. Particles in the outer disk require a
few size-doubling time scales to reach the drift limit, so we can estimate
this time of constant mass by a few growth time scales (Eq. 32) at distances
of around 1-2 characteristic radii.

After this initial growth phase (and depending on the collision model),
particles can reach large sizes where fragmentation becomes the limiting
factor. At this constant size limit (Eq. 37), also the depletion timescale is
constant which leads to an initially steep decrease in the dust mass (around
$0.1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}0.3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$
in Figure 8). This decrease in dust mass (and dust-to-gas ratio) moves the
drift limit towards smaller particle sizes until it becomes the limiting
growth barrier. This is seen in the top-left panel of Figure 7, where
fragmentation only happens inside of $1\text{\,}\mathrm{a}\mathrm{u}$ while
the entire outer disk follows the drift limit. Under these conditions, we can
expect the majority of the dust mass to drain with a speed of Eq. 41, which
therefore tends to a $t^{-1}$ dependency, as shown in Figure 8 after about
$0.5\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. If this timescale approaches the
viscous timescale, it will tend towards $t^{-1.5}$ due to the additional mass
loss through viscous transport.

If significant sub-structure is present, this process can look drastically
different, as shown in the lower panels in Figure 7 and yellow lines in Figure
8: the displayed simulation is identical apart from the fact that a
modification in the gas viscosity imposes a gap (using the gap profile of
Duffell, 2020 for Jupiter mass planet). The gap at
$20\text{\,}\mathrm{a}\mathrm{u}$ causes a pressure maximum close to
$30\text{\,}\mathrm{a}\mathrm{u}$. With the pressure gradient approaching zero
towards the maximum, radial drift as well as radial and azimuthal relative
velocities vanish. This means that only the fragmentation barrier (or a
bouncing barrier if included in the collision model) remains as growth barrier
as long as the gas turbulent velocity is high enough for fragmenting
collisions. The size limits in the lower left panel of Figure 7 indicate that
particles in the pressure maximum are fragmentation limited, while the entire
rest of the disk is already in the drift limit. The radial directions of the
arrows show that this pressure maximum moves the particles between 20 and
$30\text{\,}\mathrm{a}\mathrm{u}$ radially outward towards the pressure bump
while the dust outside the trap and inward of the gap keeps drifting radially
inward. The dust inside $20\text{\,}\mathrm{a}\mathrm{u}$ is therefore
drained, while practically all dust outside $20\text{\,}\mathrm{a}\mathrm{u}$
is steadily transported into the pressure trap. Figure 8 shows that in this
simulation, the mass evolution levels off at approximately the initial dust
mass outside the gap.

###  8.2 Global evolution of the dust surface density

The size distributions discussed previously were shown to be top-heavy, i.e.
contain most total mass in the largest grains. At the same time, the largest
grains tend to have the highest drift speeds. Therefore, the total dust flux
(sometimes called the pebble accretion rate) is dominated by the upper end of
the size distribution. This upper end can be simply approximated: initially
all particles are small, but they grow at a speed of
$\dot{a}\simeq\rho_{\mathrm{d}}\,\Delta v/(2\,\rho_{\mathrm{s}})$ (see Section
7.1), which is often dominated by turbulent relative velocities and then
becomes $\dot{a}\simeq a\,Z\,\Omega_{\mathrm{K}}/2$ (see Eq. 32). After a few
local growth time scales, the particles reach one of the size limits described
in Section 7.2, whichever is more constraining. At this point, the dust mass
flux $\dot{M}_{\mathrm{d}}=2\,\pi\,r\,\bar{v}_{d,r}\,\Sigma_{\mathrm{d}}$, is
largely set by the mass and drift speed of the largest grains, or stated
differently, the mass averaged velocity
$\bar{v}_{d,r}=(\int_{0}^{\infty}v(a)\,\sigma_{\mathrm{d}}\,\mathrm{dln}a)/\Sigma_{\mathrm{d}}$
will be within a factor of order unity of the drift speed of the largest
grains. Birnstiel et al. (2012) therefore proposed a two-population model
where the entire dust distribution is represented by only two populations of
dust: small dust that effectively follows the gas radial speed (and contains
an often negligible fraction of the mass) and a population of large dust that
moves with the drift speed of the largest grains, multiplied by a fudge factor
that is calibrated by numerical simulations.

If we ignore this order-of-unity factor and the small contribution of small
dust grains, we can (slightly over-)estimate the total dust flux as
$\dot{M}_{d}=2\pi\,r\,\bar{v}_{d,r}(a_{\mathrm{max}})\,\Sigma_{\mathrm{d}}$.
The previous subsection showed that the outermost regions of the disk retain
their dust the longest, so they provide a slowly decreasing reservoir of dust
mass and set the value of $\dot{M}$. The regions inward will mostly just
transport this dust further inside with the velocity given by the largest
grains. This allows solving the mass flux equation for the surface density,

$\Sigma_{\mathrm{d}}\simeq\frac{\dot{M}_{d}}{2\pi\,r\,\bar{v}_{d,r}(a_{\mathrm{max}})},$
(42)

where we can now insert which ever size limit sets $a_{\mathrm{max}}$. As
shown in Birnstiel et al. (2012), this would predict a dust surface density of
$\Sigma_{\mathrm{d}}\propto r^{-3/4}$ for a gas density of
$\Sigma_{\mathrm{g}}\propto r^{-1}$ if drift is the limiting growth barrier.
For the fragmentation barrier (and constant fragmentation velocity and
$\alpha$), this would always result in a profile of
$\Sigma_{\mathrm{d}}\propto r^{-3/2}$. Eq. 42 also shows that variation in the
fragmentation velocity cause the maximum particle size to vary (see Eq. 37)
and opposite variations in the dust surface density. As shown in Birnstiel et
al. (2010), lower fragmentation speeds inside the snow-line would cause an
increase of the surface density. The same mechanism was applied in Okuzumi et
al. (2016) where sintering was proposed to lower the fragmentation velocity in
the regions outside sublimation fronts and this way causing local enhancements
in the dust surface density.

##  9 THE PATH AHEAD

The previous sections drew a broad picture of the processes and timescales
involved in particle growth and transport, the earliest stages of forming
planets. The next sections will look beyond those processes towards testing
the predictions made via observational tracers (Section 9.1), towards growth
beyond the growth barriers (Section 9.2) and then conclude with some future
directions of the field.

###  9.1 Observational Tracers

In Section 1.2, some observational trends have already been introduced: Within
the last decade, it became possible to observe planet-forming disks at high
angular resolution and sensitivity, thanks in particular to the ALMA
Observatory and the SPHERE instrument at the VLT. ALMA covers a wavelength
range from a fraction of a millimeter to a few millimeters and both ALMA and
SPHERE have capabilities to observe polarized light. These capabilities enable
not just imaging the disks, but also to develop observational probes of the
underlying physics. In the following, the most important observational tracers
and trends and their impact on our picture of early planet formation will be
discussed. For more details, the reader may refer to the review in this series
by Andrews (2020).

[] Population propertiesObservable properties of the stars or disks that can
be measured for large numbers of sources, including disk integrated fluxes,
spectral indices, disk sizes, dust/gas masses, or stellar masses, ages, and
accretion rates. Effective radiusThe radius encompassing a defined fraction of
the total flux, often 68% or 99%, depending on application. AlbedoThe (single-
scattering) albedo
$\omega=\bar{\kappa}_{\nu,\mathrm{sca}}/(\bar{\kappa}_{\nu,abs}+\bar{\kappa}_{\nu,sca}$
is the probability of a photon being scattered instead of being absorbed when
interacting with a dust grain.

Figure 9: Mass absorption and scattering opacity (top row, at
$1.25\text{\,}\mathrm{m}\mathrm{m}$), spectral indices (middle row, measured
between 1.25 and $3\text{\,}\mathrm{m}\mathrm{m}$), and scattering properties
of dust particles as function of maximum particle size $a_{\mathrm{max}}$,
averaged over a size distribution with size exponent of $q=2.5$ (left) and
$q=3.5$ (right). Lines correspond to default DSHARP (green), DSHARP with
different carbonaceous material and computational method (orange, labeled
DSHARP-ZBED) and default DIANA (purple). Solid and dashed lines in the top row
correspond to absorption and scattering, respectively. The bottom panels show
the albedo $\omega$ (dashed), the probability of perpendicular scattering
$P_{90}$ (dotted), and the product of both (solid) at
$1.25\text{\,}\mathrm{m}\mathrm{m}$. Gray lines show how these values are
shifted for a wavelength of $3\text{\,}\mathrm{m}\mathrm{m}$.

####  9.1.1 Dust opacities and dust spectral properties

In the context of dust continuum emission, the spectral index is defined as
the exponent of the spectral energy distribution in the (sub-)millimeter
range, $F_{\nu}\propto\nu^{\alpha}_{\mathrm{mm}}$, sometimes given with a
specific wavelength range over which it was measured, such as
$\alpha_{\mathrm{1-3mm}}$. This should not be confused with $\alpha$ (without
subscript), that denotes the turbulence parameter in this review. In a simple
slab model, neglecting the effects of scattering or background, the emitted
intensity is

$I_{\nu}=B_{\nu}(1-\exp(-\tau_{\nu}))=B_{\nu}\,\times\begin{cases}1&\text{for
}\tau_{\nu}\gg 1\\\ \tau_{\nu}&\text{for }\tau_{\nu}\ll 1,\end{cases}$ (43)

where $B_{\nu}$ and $\tau_{\nu}$ are the Planck function and the absorption
optical depth along the line of sight. The latter is often well approximated
by $\tau_{\nu}=\Sigma_{\mathrm{d}}\,\bar{\kappa}_{\nu}/\cos i$, where
$\bar{\kappa}_{\nu}$ is the opacity (absorption cross-section per gram of
dust, units of $\text{\,}\mathrm{c}\mathrm{m}^{2}\,\mathrm{g}^{-1}$), averaged
over the dust size distribution
($\bar{\kappa}_{\nu}=\int\sigma_{\mathrm{d}}(a)\,\kappa_{\nu}(a)\mathrm{dln}a/\Sigma_{\mathrm{d}}$),
and $i$ the inclination.

If the opacity follows at least locally a power-law,
$\bar{\kappa}_{\nu}\propto\nu^{\beta}$, the frequency dependency of $I_{\nu}$
in the Rayleigh-Jeans limit becomes $\alpha_{\mathrm{mm}}=2+\beta$ for
optically thin emission and $\alpha_{\mathrm{mm}}=2$ for optically thick
emission (where the number 2 might be smaller if the emission is not deep
enough in the Rayleigh-Jeans limit). Figure 9 shows the optical properties of
two particle size distributions (with size exponents $q=2.5$ and $q=3.5$) of
different materials and employing different computational methods. DSHARP
denotes the standard choice of Birnstiel et al. (2018) (mostly backward-
compatible with D’Alessio et al., 2001), DSHARP-ZBED uses the same mass
fractions, but uses different optical constants for the carbonaceous
materials, Zubko et al. (1996, sample ‘BE’), instead of the ones by Henning &
Stognienko (1996) and furthermore uses DHS (Min et al., 2005, with their
parameter $f=0.8$) to compute the opacities, instead of classical Mie theory.
DIANA denotes the standard opacities of Woitke et al. (2016). This shows that
the choice of carbonaceous material, and to some extent also the computational
method, can strongly affect the shape and strength of the absorption opacity
with very limited effect on the scattering opacity (see top row in Figure 9).
The particle size distribution exponent $q$ affects particularly the opacity
spectral index $\beta$ for large maximum particle sizes (comparing the left
and right columns of Figure 9). For the chosen opacities, only the top-heavy
particle size distribution with $q=2.5$ (which resembles the radial drift
limited simulations shown in Figure 7) reaches very low values of $\beta$. For
a given (uncertain) assumption of the opacities, measurements of the spectral
index can thus in principle be used to constrain the maximum particle sizes in
the sub-millimeter to centimeter range (Ricci et al., 2010, Isella et al.,
2010, Pérez et al., 2012, e.g). This approach is further complicated by the
unknown optical depth, which can mimic a low dust opacity spectral index via
high optical depth as both lead to $\alpha_{\mathrm{mm}}\sim 2$ (Ricci et al.,
2012). The measured low values of $\alpha_{\mathrm{mm}}$ (see next section)
are indicative of significant grain growth, high optical depth, or likely a
combination of both.

If scattering is considered (Miyake & Nakagawa, 1993, Birnstiel et al., 2018),
the resulting intensity in the optically thin region is increased (a photon
takes a longer path and appears more thermalized with the medium), while the
optically thick emission is reduced (see Figure 7 of Miotello et al., 2023, or
also Zhu et al., 2019). The amount of reduction in the latter case depends on
the albedo
$\omega=\bar{\kappa}_{\nu,\mathrm{sca}}/(\bar{\kappa}_{\nu,\mathrm{abs}}+\bar{\kappa}_{\nu,\mathrm{sca}})$.
At wavelengths, where the albedo is strongly wavelength-dependent (see bottom
row in Figure 9, near $0.2\text{\,}\mathrm{m}\mathrm{m}$), this might lead to
spectral indices even below $\alpha_{\mathrm{mm}}=2$ (Liu, 2019).

For increasing maximum particle size around
$a_{\mathrm{max}}\sim\lambda/(2\pi)$ the probability of scattering
($\omega_{\nu}$) increase while at the same time scattering becomes more
directed into the forward direction. Perpendicular scattering produces the
most polarized intensity, but the probability of perpendicular scattering,
$P_{90}$, decreases with $a_{\mathrm{max}}$. The combination of scattering
becoming more likely, but less isotropic leads to a sharp peak in the degree
of polarized intensity (Kataoka et al., 2015), as shown in the bottom row of
Figure 9. This polarization by self-scattered dust emission at different
wavelength allows setting tight constraints on the dust grain size at around
$0.1\text{\,}\mathrm{m}\mathrm{m}0.2\text{\,}\mathrm{m}\mathrm{m}$ (Kataoka et
al., 2016, 2017), although this picture is complicated by other mechanisms
that can cause polarization at radio wavelengths (Kataoka et al., 2017, 2019).

At the same time, multi-wavelength analysis of the continuum emission can be
used to constrain the maximum particle size, dust surface density, and
temperature profiles of disks, (e.g. Carrasco-González et al., 2019, Ueda et
al., 2020, Sierra et al., 2021, Guidi et al., 2022) often finding slightly
larger particle sizes around $1\text{\,}\mathrm{m}\mathrm{m}$. This
inconsistency between measuring particle sizes with (1) the polarization
signal and (2) the continuum spectral index might be solved by the
contributions of optically thick emission, as proposed in Lin et al. (2020)
for HD 163296. Modeling the rings in this disk, Doi & Kataoka (2023) computed
much larger values, however, they also show how much the absolute values of
$a_{\mathrm{max}}$ depend on the choice of opacity. In the future, it will be
instructive to bring different methods of measuring particle sizes in
agreement when applied to the same disk. Doi & Kataoka (2023) also derive
constraints on the turbulence parameter, finding relatively high values of
$\alpha\sim$6\text{\times}{10}^{-2}\text{\,}$$ in the vertical direction,
while similar approaches for other disks indicate extreme levels of settling,
consistent with $\alpha<10^{-4}$ (Villenave et al., 2022). This shows the
diversity of disk physical properties which will require studying not only
individual sources, but also population statistics of disks in large numbers,
as discussed in the next section.

####  9.1.2 Population properties

Population properties such as integrated fluxes or radii were accessible
already pre-ALMA, but population studies benefitted strongly from the high
sensitivity of ALMA. Studies from Andrews et al. (2013, for the Taurus star
forming region), Ansdell et al. (2016, for Lupus), Barenfeld et al. (2016, for
Upper Scorpius), Pascucci et al. (2016, for Chamaeleon I), Ansdell et al.
(2018, $\sigma$ Orionis), and Tychoniec et al. (2018, of Class 0/I sources in
Perseus) and many others brought the numbers of observed disks close to 1000
(see Miotello et al., 2023 and Manara et al., 2023 for recent reviews). These
studies allow moving beyond studies of individual disks and enable probing
general correlations or age trends that any theory of disk evolution needs to
be able to explain. Tripathi et al. (2017) for example discovered a disk-size-
flux relation, often called size-luminosity-relation (SLR, see also Andrews et
al., 2018b) in which the integrated dust continuum flux scales very closely
with $r_{\mathrm{eff}}^{2}$. This is surprising since disks do not display a
radially constant surface brightness. Rosotti et al. (2019a) and Zormpas et
al. (2022) showed that drift-limited particle growth could naturally explain
this trend. However, this seems at odds with the fact that significant
substructure is deemed necessary to explain (1) high-resolution images (e.g.
Andrews et al., 2018a, Dullemond et al., 2018) as well as (2) the generally
low spectral indices of disks, shown in the left panel of Figure 10. For
comparison, two evolutionary tracks of two simulations are shown overlaid on
an observed population of disks (from Testi et al., 2014, Tazzari et al.,
2021, Andrews et al., 2018b). The green line depicts a simulation with a
strong pressure trap (as caused by a Jupiter-mass planet, based on the profile
of Duffell, 2020, corresponding to the lower row of Figure 7). The orange line
shows the identical simulation but without a pressure trap (corresponding to
the upper row in Figure 7). Since radial drift preferentially removes the
large grains, the spectral index in the latter case increases as the dust mass
(and therefore also the millimeter flux) decreases. It can be seen that the
simulation with pressure trap is broadly consistent with the observed values.
However, large numbers of simulations are required to draw firm conclusions
and lead to the expectation that sub-structured disks would follow a shallower
size-luminosity trend as observed (Zormpas et al., 2022).

It is remarkable that none of the observed disks shows high spectral indices
towards lower-flux-ranges. This is indicating that dust grains remain large
and/or optically thick throughout their later evolution, at least at the flux
levels where spectral indices are currently available. The very low spectral
indices below $\alpha_{\mathrm{mm}}=2$ are partially attributed to errors, but
some can be caused by the effects of scattering (see Zhu et al., 2019, Liu,
2019, and the discussion above).

While population synthesis of disks is still in the infancy (see Rosotti et
al., 2019b, Somigliana et al., 2020, Zagaria et al., 2022, Zormpas et al.,
2022, for some recent work), future efforts are needed to find the initial
conditions and the disk evoluationary processes that can explain the large
ranges of population properties available today (gas and dust mass
constraints, multi-wavelength-dust and gas disk sizes, spectral indices,
accretion rates, and more). A more in-depth discussion of disk populations and
the measurements of disk fundamental properties can be found in Manara et al.
(2023) and Miotello et al. (2023).

Figure 10: Left: Distribution of millimeter fluxes and spectral indices
compared to simulations of dust evolution simulations. Black and gray dots are
the samples of Testi et al. (2014) and Tazzari et al. (2021), respectively.
Their combined kernel density estimate distribution is shown in color scale
where the solid/dashed black contour encompasses 25 and 68%. Right: Size-
Luminosity correlation (gray lines) of Tripathi et al. (2017) and Andrews et
al. (2018b) and the samples from which the correlation was derived (points and
kernel density estimate). In both panels, colored lines depict the
evolutionary tracks of the simulations shown in Figure 7 and described in the
text box “the standard disk”, where numbers denote the simulation time in
millions of years. The green line is a simulation with an efficient dust trap,
the orange a simulation without a trap.

###  9.2 The growth towards planetesimals

Given the size limits set by bouncing, fragmentation, erosion, or radial
drift, it might seem as growth towards gravitationally bound bodies, called
planetesimals is impossible. There are however proposed pathways of growth
beyond the size barriers.

Firstly, _lucky particles_ were proposed by Windmark et al. (2012a), Windmark
et al. (2012b), and Garaud (2007). As discussed in Section 5, this relies on
the fact that small particles colliding at tens of meters per second may
deposit mass onto much larger target particles (e.g. Teiser & Wurm, 2009), so
a single larger particle might be able to grow in an environment where all
other particles remain limited to small sizes by a barrier. Producing such a
single larger particle might be possible given the randomness in the collision
speed (or other collision parameters), such that a particle could happen to be
always be colliding with particularly beneficial outcome, giving rise to the
term _lucky_ particles. Such a scenario is very improbable, but given the
large numbers of particles, conceivable. However, there are two main caveats
of this idea. Firstly, the radial drift barrier cannot be circumvented by
lucky collisions as it is not a collisional, but a dynamical barrier, it would
have to be prevented by a pressure trap. Secondly, even if the large particle
is able to grow, its growth timescale would be very long, too long to
efficiently produce the planetessimals needed for the formation of planets
(e.g. Estrada et al., 2016).

Secondly, _fluffy growth_ was proposed as a way to grow quickly enough to
bypass the radial drift barrier. When including porosity evolution in models
of particle growth and transport (e.g. Okuzumi et al., 2009), it was found
that large particles in the inner disk can grow through the drift barrier
(Okuzumi et al., 2012). This is enabled by the fact that particles in the
inner disk, if gas densities are high enough, can experience Stokes drag
instead of the Epstein drag (see Section 3.1). In that case, the growth
timescale is not independent of size (as in Eq. 32), but inversely
proportional to the particle size. Therefore, particles with a low fractal
dimension and consequently large cross-section can reach very short growth
time scales (Okuzumi et al., 2012, Kataoka et al., 2014). This, however,
requires that the particles do not fragment, even at tens of meters per second
(see Krijt et al., 2015). Furthermore, highly porous particles currently do
not seem to be able to reproduce observed disk populations (Zormpas et al.,
2022) due to their low opacity (Kataoka et al., 2014). Still, more work is
required to draw firm conclusions about particle porosity.

The currently favored scenario of planetesimal formation is gravitational
collapse of a dust cloud, similar to the original works of Safronov (1972) and
Goldreich & Ward (1973). The challenge in this scenario are the many orders of
magnitude difference between the initial dust volume density and the density
that is required for the dust cloud to be gravitationally bound, which happens
at around the Hill density $\rho_{\mathrm{H}}\simeq 9M_{\star}/(4\pi\,r^{3})$.
At this point, the gravitational collapse can proceed until material densities
of order unity in CGS units are reached, unless effects on smaller scales
intervene (Klahr & Schreiber, 2021). Figure 11 displays the dust, gas and Hill
densities in the mid-plane of a disk at
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. The dust, at the canonical dust-
to-gas ratio of 1% needs to be accumulated by more than four orders of
magnitude. As discussed in Section 4.3, dust is sedimenting towards higher
pressure, so in a ring of high pressure, dust is accumulated in the radial and
vertical direction, as seen at around $20\text{\,}\mathrm{a}\mathrm{u}$ in
Figure 11. This can raise the dust-to-gas ratio towards unity. At this point,
the dust mass and momentum are not negligible compared to the gas and
therefore, the dust starts to strongly affect the gas dynamics (e.g. Fu et
al., 2014, Taki et al., 2016, Onishi & Sekiya, 2017, Dipierro et al., 2018,
Raettig et al., 2021). Initially, it was thought that those dynamical affects
are detrimental to planet formation (Cuzzi et al., 1993) as the dust-gas shear
motion would trigger turbulence which through turbulent mixing inhibits
further concentration of the dust. The other face of turbulence is, however,
that it does not act as a homogeneous diffusion on all scales, but that it can
also form localized pockets of high and low densities. Under the right
conditions (see below), one effect coming into play is the Streaming
Instability (SI) (Youdin & Goodman, 2005, Johansen et al., 2007) which is an
instability from a family of resonant drag instabilities (Squire & Hopkins,
2018). Its basic mechanism works as follows: radially drifting dust drags gas
radially along, but angular momentum conservation causes the gas to increase
its azimuthal velocity. This increased orbital gas speed drags dust forward in
azimuth, accelerating it in the azimuthal direction and thus pushing it
outward. This enhances the original dust enhancement and becomes a run-away
effect. A detailed review of the SI related literature is beyond the scope of
this review (see Simon et al., 2022 for a recent review) and only the key
findings are summarized in the following: While this linear phase of the
classical streaming instability is reasonably well understood (Squire &
Hopkins, 2020), the open questions lie in its interplay with other effects
(external turbulence, particle size distributions) and in the clumping in the
non-linear phase of the instability.

It was found that the SI linear growth rate (Auffinger & Laibe, 2018, Umurhan
et al., 2020) and strong clumping (Gole et al., 2020) can be significantly
reduced by viscosity. In contrast, Yang et al. (2018) showed that strong
clumping can occur despite significant stirring in an MRI dead zone. This was
confirmed by Xu & Bai (2022) where moderate MRI turbulence levels even helped
trigger the SI. Other mechanisms to drive turbulence, such as the VSI were
shown to be able to coexist with the SI (Schäfer et al., 2020).

When particle size distributions are included, Krapp et al. (2019) found that
the SI can be suppressed (see also Bai & Stone, 2010, Zhu & Yang, 2021 and
Paardekooper et al., 2020). In contrast, Schaffer et al. (2021) find that dust
size distributions, although changing the dynamics of the system, do not
suppress strong clumping, which the authors contribute to the fact that
vertical settling leads in all cases to mid-plane concentrations where also
linear analysis finds sufficiently high growth rates.

Our most current understanding (despite the caveats mentioned before, which
clearly show that more work is required in the future), still remains that
_high_ dust-to-gas ratios and _large_ particles are required to trigger strong
clumping. What exactly _high_ and _large_ means was consolidated in Simon et
al. (2022), based on previous works by Carrera et al. (2015), Yang et al.
(2017), and Li & Youdin (2021), which is shown in the right panel of Figure
11. Dust to gas surface density ratios of
$\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$ above 3% and/or Stokes numbers above
0.01 are generally susceptible to strong clumping and dust evolution models
allow the production of such sizes in pressure traps. In addition to that, the
conditions for SI and strong clumping need to be retained long enough: Carrera
& Simon (2022), for example, showed that a wide pressure bump in their
simulation was not able to trigger strong clumping despite far exceeding the
required criteria. As the authors state, this is likely due to the fact that
the particles drift through the region of high concentration on time scales
shorter than the growth timescale of the instability.

There do exist alternative proposals to concentrate particles to the onset of
gravitational collapse: on small scales turbulent eddies can centrifugally
eject particles, causing concentrations _between_ eddies, opposite to vortex
trapping discussed in Section 4.3. This idea, proposed in Cuzzi et al. (2008),
was criticized as unlikely in more realistic models of turbulence by Pan et
al. (2011), but it remains the subject of ongoing discussions (Hartlep &
Cuzzi, 2020, Klahr & Schreiber, 2020). A second alternative is called the
secular gravitational instability, where the self-gravity of a dust ring is
aided by gas drag which can damp the dust eccentricy and enable continued
growth of the initial perturbation (e.g. Ward, 2000, Takahashi & Inutsuka,
2014). Recent works on this mechanism includes Tominaga et al. (2020) and Abod
et al. (2019) which show that the boundaries between the mechanisms get
blurred in a pressure bump and that on small scales, the collapse is likely
governed by the diffusing effects of turbulence (see Klahr & Schreiber, 2020,
Gerbig et al., 2020).

In summary, most models agree that accumulations of partially decoupled
particles can trigger planetesimal formation. This means that some amount of
particle growth is required, in order to provide particles large enough to
participate in these mechanisms. Global models of dust evolution routinely
produce particles significantly larger than the initially well coupled
monomers. It is therefore likely that coagulation models (that mostly fail to
produce planetesimals) and models of gravitational collapse of pebble-
overdensities (that would not work with small dust grains) go hand-in-hand. In
order to model planetesimal formation starting from small dust, therefore
requires recipes that predict when the conditions for planetesimal formation
(that is strong clumping) are reached. In doing so, a wide variety of models
exist, see the discussion in Drążkowska et al. (2023). The right panel in
Figure 11 presents the combined criteria of Li & Youdin (2021) and Simon et
al. (2022) in comparison to a widely-used criterion which requires the mid-
plane dust-to-gas ratio to reach unity (Drążkowska et al., 2016, Schoonenberg
& Ormel, 2017, and others). At the bottom, the ranges of Stokes numbers are
displayed that are reached in fragmentation-limited growth for fragmentation
velocities between 1 and $10\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$. It can be
seen that for a high metallicity, there always exists a range of Stokes
numbers that can contribute to planetesimal formation. Such high metallicities
$Z$ are however only reached through some form of dust-concentrating
mechanism. This means that the planetesimal formation process can bee seen as
a three stage-process:

  * •

Dust needs to grow via collisions to sizes that are partially decoupled from
the gas, often called pebbles. There is no fixed size associated with this
process, it is rather understood as particles with a Stokes number
$\mathrm{St}>\alpha$.

  * •

Pebbles are susceptible to radial, vertical, and azimuthal drift and can
therefore be efficiently trapped, or through radial drift efficiently be
supplied to the inner disk regions. Some mechanism accumulates the dust beyond
the percent-level. A range of potential mechanisms exist such as radial drift
(e.g. Youdin & Shu, 2002), pressure traps (e.g. Whipple, 1972, Klahr &
Henning, 1997, Rice et al., 2006, Pinilla et al., 2012, and many others) or
dead-zone boundaries (e.g. Kretke & Lin, 2007, Lyra et al., 2008, Brauer et
al., 2008b, Dzyurkevich et al., 2010, Pinilla et al., 2016, Charnoz et al.,
2021), zonal flows (Johansen et al., 2009, Bai & Stone, 2014), the secular
gravitational instability (e.g. Takahashi & Inutsuka, 2014, Tominaga et al.,
2020), or traffic-jam effects at the water snow line (Birnstiel et al., 2010,
Saito & Sirono, 2011) or outside of condensation fronts through sintering
(Okuzumi et al., 2016) or aided by dust back-reactions (Drążkowska & Alibert,
2017, Schoonenberg & Ormel, 2017, Gárate et al., 2020).

  * •

If the dust-to-gas ratio near the mid-plane approaches unity, and conditions
for _strong clumping_ are fulfilled, pockets of high dust density can exceed
the Hill density and further collapse to planetesimals is possible. This
closes the gap between the high dust-to-gas ratio in the pressure trap (see
Figure 11) and the even higher Hill density that is required for gravitational
collapse.

At this point is worth revisiting the observations of disk sub-structure in
disks: the optical depth at millimeter wavelength appears universally within
the order of $\sim 0.4$, (Cazzoletti et al., 2018, Isella et al., 2018,