# COMETS AND PLANETESIMAL FORMATION

Jacob B. Simon Department of Physics and Astronomy, Iowa State University,
2323 Osborn Dr., Ames, IA 50011-1026, USA Jürgen Blum Institut für Geophysik
und extraterrestrische Physik, Technische Universität Braunschweig,
Mendelssohnstr. 3, D-38106 Braunschweig, Germany Til Birnstiel University
Observatory, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679
Munich, Germany Exzellenzcluster ORIGINS, Boltzmannstr. 2, D-85748 Garching,
Germany David Nesvorný Department of Space Studies, Southwest Research
Institute, 1050 Walnut St., Suite 300, Boulder, CO 80302, USA

###### Abstract

  * In this chapter, we review the processes involved in the formation of planetesimals and comets. We will start with a description of the physics of dust grain growth and how this is mediated by gas-dust interactions in planet-forming disks. We will then delve into the various models of planetesimal formation, describing how these planetesimals form as well as their resulting structure. In doing so, we focus on and compare two paradigms for planetesimal formation: the gravitational collapse of particle over-densities (which can be produced by a variety of mechanisms) and the growth of particles into planetesimals via collisional and gravitational coagulation. Finally, we compare the predictions from these models with data collected by the Rosetta and New Horizons missions and that obtained via observations of distant Kuiper Belt Objects.

##  1 INTRODUCTION

The formation of planetesimals from small sub-micron grains encompasses
approximately 10–12 orders of magnitude increase in scale and numerous
physical processes, ranging from particle-particle sticking to angular
momentum exchange between disk gas and solid particles. While there are many
unanswered questions in how exactly these physical processes work, recent
years have seen substantial progress in addressing several key questions:

  1. 1.

How do sub-micron size grains grow to the larger mm–cm solids (referred to as
pebbles; see 2.5 for a more complete definition) that have been
observationally inferred to exist in disks around young stars?

  2. 2.

How do these pebbles continue their growth to forming km-sized planetesimals?

  3. 3.

What are the properties of the resulting planetesimals and how do they compare
with observations of Solar System planetesimal populations, such as comets and
Kuiper Belt Objects (KBOs)?

Along the theoretical front, a number of new insights have been gained. In
particular, the powerful combination of analytical techniques and numerical
simulations have identified a number of instabilities (e.g., _Youdin and
Goodman_ , 2005; _Youdin_ , 2005) and mechanisms (e.g., _Cuzzi et al._ , 2008;
_Hopkins_ , 2016) that can concentrate small grains into regions of sufficient
density to gravitationally collapse into planetesimals.

These techniques have directly led to predictions that can be tested with
current and forthcoming observational campaigns, and along this front,
observational data has also been abundant. With a number missions (e.g.,
OSIRIS-REX, Rosetta, New Horizons) to small Solar System bodies (e.g., _Stern
et al._ , 2021), as well as large surveys of small body populations, such as
the Outer Solar System Origins Survey (OSSOS), a plethora of new data has
arrived on the physical characteristics of Solar System planetesimals (e.g.,
_Kavelaars et al._ , 2021; _Fraser et al._ , 2021). Furthermore, the advent of
next generation optical/IR and radio telescopes, such as the Atacama Large
Millimeter/submillimeter Array (ALMA) and the Very Large Telescope (VLT) have
revolutionized our understanding of both gas and small particle dynamics in
planet-forming disks around nearby stars (_ALMA Partnership et al._ , 2015;
_Isella et al._ , 2016; _Huang et al._ , 2018).

Finally, laboratory experiments have been very fruitful in addressing how
exactly solid particles of different compositions interact when they collide
and have thus served as a strong complement to both the observational and
theoretical/numerical studies (_Blum and Wurm_ , 2008). Moreover, the physical
properties of planetesimals formed by different hypothesized mechanisms have
also been studied in the laboratory (_Blum_ , 2018), allowing a comparison to
those of the most primitive bodies in the Solar System.

Indeed, the combination of theory, observations, and laboratory work has
proved to be very successful in improving our understanding of planetesimal
and comet formation. In this chapter, we review such work and detail the
latest understanding of how planetesimals are born. Our approach is primarily
pedagogical and as such, much of our focus will be on the physics of
planetesimal formation itself, while making connections to both observations
and laboratory experiments when appropriate.

##  2 DUST GRAIN GROWTH

In this section, we describe in detail our current understanding of how small
solids grow in protoplanetary disks (PPDs). Some discussion of dust growth is
discussed in Chapter 2 by Aikawa et al. However, here, we delve into the
details and describe the processes that limit growth of small grains beyond
mm-cm sizes, thus setting the stage for many of the mechanisms described in
Section 3.

The evolution of solids in PPDs is influenced by many effects, but they can be
broadly categorized into the two areas of 1) the evolution of their
composition, size, and morphology and 2) their transport and dynamics within
protoplanetary disks. It is important to realize that one cannot be treated
separately from the other, because these two categories affect each other in
many ways: the make-up, size and morphology of a particle affects its
aerodynamical properties, and the aerodynamic properties determine how a
particle is transported throughout the disk and how its environment changes
along the way. The environment, in turn, determines how collisions and
composition evolve: it sets the collision speeds (see, e.g., _Birnstiel et
al._ , 2016) and/or condensation/sublimation rates (e.g. _Stammler et al._ ,
2019). In the following section, we will describe how dust is transported and
how it evolves through collisions and will highlight how those aspects affect
each other.

We begin by introducing a crucial quantity to describe the evolution from
small grains to a planetesimal: the particle stopping time. It is defined as
the time scale on which the particle velocity adapts to the gas velocity
through drag forces and reads

$t_{\mathrm{stop}}=\frac{m\,\Delta
v}{F_{\mathrm{drag}}}=\frac{a\,\rho_{\mathrm{s}}}{\rho_{\mathrm{g}}v_{\mathrm{th}}},$
(1)

where $m$, $a$, $\rho_{\mathrm{s}}$, $\Delta v$, and $v_{\mathrm{th}}$ are the
particle mass, its radius, its material density, the absolute value of the
velocity difference between gas and dust, and the mean thermal velocity of the
gas (related to the gas sound speed through a factor of order unity),
respectively. The right equality in Eq. 1 assumes the drag force to be in the
Epstein regime (_Epstein_ , 1924; _Weidenschilling_ , 1977) where the particle
size is smaller than the mean-free path of the gas molecules and furthermore
that the dust particles are spherical and non-fractal. Often, the stopping
time is expressed by the Stokes number111A more accurate definition of the
Stokes number is discussed in 3.1.5. However, for our current purposes, this
definition will suffice.

$\mathrm{St}=t_{\mathrm{stop}}\cdot\Omega,$ (2)

with $\Omega$ being the Kepler frequency at the location in the disk. The
Stokes number is thus a dimensionless number specifying the aerodynamic
behavior of the particle. The Stokes number can therefore be seen as a
“dimensionless particle size” or is sometimes also called the dimensionless
stopping time as particles with $\mathrm{St}\ll 1$ are coupling much faster to
the gas flow than the orbital time scale, while particles with $\mathrm{St}\gg
1$ require many orbits to couple.

It is this coupling by drag forces that gives rise to systematic drift motion
of particles (_Whipple_ , 1972; _Nakagawa et al._ , 1986), determines how
particles are affected by turbulence (_Voelk et al._ , 1980; _Ormel and Cuzzi_
, 2007), and sets how strongly the radial gas flow drags particles along
(_Takeuchi and Lin_ , 2002). In the following, we will discuss some of these
dynamical effects in more detail.

###  2.1 Vertical Settling

A particle on an inclined orbit (and in the absence of gas) would be
oscillating vertically above and below a reference plane where the vertical
coordinate has its origin, which we will refer to henceforth as the disk “mid-
plane". This can be seen from the equation of motion of a particle near the
mid-plane (_Nakagawa et al._ , 1986),

$\ddot{z}+z\,\Omega^{2}+\frac{\dot{z}}{t_{\mathrm{stop}}}=0,$ (3)

where the first two terms describe a harmonic oscillator. The remaining term
is the velocity dependent drag term. The presence of the gas disk thus causes
a deceleration of the particle if it moves vertically up and down through the
gas, and the vertical motion of the particle can be regarded as a damped
harmonic oscillator. Particles with $\mathrm{St}\ll 1$ adapt quickly to the
gas velocity, and in this limit, the particles are not oscillating; instead
their velocity quickly tends to vertical gas velocity that we here assume to
be zero. Due to the vertical component of gravity that accelerates the
particle to the mid-plane, the particle velocity approaches a steady state
between the height-dependent gravitational acceleration and the velocity
dependent deceleration. Within a few stopping times, this equilibrium
($\ddot{z}=0$) is reached at a terminal velocity of

$v_{\mathrm{sett}}=-z\,\Omega\,\mathrm{St}\quad(\text{for }\mathrm{St}\ll 1).$
(4)

This is analogous to dropping a feather that quickly reaches a constant
velocity.

If this vertical sedimentation were to proceed unhindered, then all particles
would be accumulating at the mid-plane in a razor-thin disk, similar to
Saturn’s rings. This is however not the case, as current scattered-light
images of PPDs spectacularly depict (e.g. _Avenhaus et al._ , 2018). Some form
of gas motion (likely vertically diffusive turbulence, though bulk flows are
also possible) is believed to be at play that vertically mixes the particles
away from the mid-plane. In this case, a balance is reached when settling
leads to dust concentration gradients that in turn lead to a diffusive flux
that becomes equal, but opposing to the settling flux,

$\rho_{\mathrm{d}}\,v_{\mathrm{sett}}-D\,\rho_{\mathrm{g}}\frac{\partial\rho_{\mathrm{d}}/\rho_{\mathrm{g}}}{\partial
z}=0,$ (5)

where the diffusion flux acts to smooth out gradients in the composition as
the dust is mixed by and well coupled to the turbulent gas. The solution of
the above equation is a steady state dust distribution

$\rho_{\mathrm{d}}(z)=\rho_{\mathrm{d}}(0)\frac{\rho_{\mathrm{g}}(z)}{\rho_{\mathrm{g}}(0)}\,\exp\left(\int_{0}^{z}\frac{v_{\mathrm{sett}}}{D}\mathrm{d}z^{\prime}\right).$
(6)

The Schmidt number $\mathrm{Sc}$ is the ratio of gas and dust
diffusivities222see the discussion in (_Youdin and Johansen_ , 2007). All
particles with $\mathrm{St}\ll 1$ are well enough coupled to the turbulent
eddies of the gas, leading to a Schmidt number of effectively unity. _Youdin
and Johansen_ (2007) showed, that it can be approximated by
$\mathrm{Sc}=1+\mathrm{St}^{2}$. If the gas follows the isothermal hydrostatic
solution of a Gaussian curve with scale height $H$, then around the mid-plane,
the vertical density profile of particles of a given size also follows a
Gaussian curve but with a scale height of

$H_{\rm p}=H\sqrt{\frac{\alpha}{\alpha+\mathrm{St}}}.$ (7)

where $\alpha$ is a quantification of the turbulence. The history and various
definitions of $\alpha$ are long and not within the scope of this chapter.
Furthermore, $\alpha$ as used here should be distinguished from the so-called
“turbulent viscosity” of _Shakura and Syunyaev_ (1973). The latter quantity
depends on correlations between radial and azimuthal velocity perturbations
and includes turbulent magnetic stresses. The former is purely a dimensionless
turbulent diffusion parameter, such as has been quantified by e.g., _Johansen
and Klahr_ (2005) and _Fromang and Papaloizou_ (2006). For our purposes, it
suffices to assume that the turbulence is isotropic and define the turbulent
diffusion as $\alpha=\delta v^{2}/c_{\mathrm{s}}^{2}$,333More accurately,
$\alpha$ should include the eddy turnover time in units of $\Omega^{-1}$.
However, at the typical driving length scales for disk turbulence, the eddy
turn over time is $\sim\Omega^{-1}$, making the dimensionless turnover time
equal to unity. where $\delta v$ and $c_{\mathrm{s}}$ are the turbulent
velocity and sound speed, respectively.

The dust to gas ratio at the mid-plane is therefore enhanced for particles
with a Stokes number $\mathrm{St}>\alpha$, since

$\frac{\rho_{\mathrm{d}}(0)}{\rho_{\mathrm{g}}(0)}=\frac{\Sigma_{\mathrm{d}}}{\Sigma_{\mathrm{g}}}\sqrt{\frac{\mathrm{St}}{\alpha}+1},$
(8)

where $\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$ is the ratio of the dust and
gas column densities.

###  2.2 Radial Drift

A basic understanding of vertical settling now established, we can apply a
similar approach for the radial and azimuthal velocity of dust and gas. The
equation of motion of dust and gas are coupled by drag forces and this set of
equations can be solved for a steady state, i.e. the terminal velocities, see
_Nakagawa et al._ (1986). Summarizing these results, the gas disk in a force
balance between pressure, gravity, and centrifugal forces rotates slightly
sub-Keplerian for a typical, radially decreasing pressure profile. This leads
to the problem that the dust particle (that feels negligible pressure
acceleration) needs to be on a Keplerian orbit to reach force balance. The
resulting headwind on the particles decelerate them and remove some of their
angular momentum, which in turn causes them to spiral inwards radially. The
gas, which has gained angular momentum from the particles, moves radially
outward, and an equilibrium is reached, as derived by _Nakagawa et al._ (1986)
(see also _Whipple_ , 1972 and _Weidenschilling_ , 1977), where the dust
radial velocity is

$v_{\mathrm{r,dust}}=-\frac{2}{(1+\epsilon)^{2}\,\mathrm{St}^{-1}+\mathrm{St}}\,\eta\,v_{\mathrm{K}},$
(9)

$\epsilon=\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$ is the dust-to-gas mass ratio,

$\eta=-\frac{1}{2}\left(\frac{H}{r}\right)^{2}\frac{\partial\ln P}{\partial\ln
r},$ (10)

and $P$ is the gas pressure. The resulting time scale on which particles drift
inward is therefore (for $\epsilon,\mathrm{St}\ll 1$)

$t_{\mathrm{drift}}=\frac{1}{\mathrm{St}\,\Omega}\left(\frac{H}{r}\right)^{-2}\left|\frac{\partial\ln
P}{\partial\ln
r}\right|^{-1}\simeq\frac{64~{}\mathrm{orbits}}{\mathrm{St}}\,\left(\frac{H/r}{0.03}\right)^{-2},$
(11)

which is long compared to the orbital time scale, but short in terms of the
overall disk evolution time scale. This shows that growing particles can
travel substantial distances during their evolution.

Both radial drift and vertical settling are a result of the terminal velocity
approximation that more generally states (see _Youdin and Goodman_ , 2005)
that particles drift with respect to the gas towards higher pressure, i.e.
that the dust-gas velocity difference is $-(\nabla
P/\rho_{\mathrm{g}})\,t_{\mathrm{stop}}$.

###  2.3 Collisional Growth and Destruction

When two dust grains collide, their fate is determined by the collision energy
as well as the material and makeup of the particles. In general, three
collisional outcomes are possible: sticking, bouncing, and fragmentation,
although mixtures of and transitions between these outcomes are possible too
(_Güttler et al._ , 2010, see also Figure 1).

The initial grain sizes can be estimated by comparison with observations of
meteorites, comets or interstellar dust. _Vaccaro et al._ (2015) investigated
the size distribution of matrix grains in primitive carbonaceous chondrites,
using high-resolution electron microscopy, and found geometric mean diameters
on the order of 200 nm. _Mannel et al._ (2019) determined very similar mean
particle sizes for the dust particles collected from comet 67P/Churyumov-
Gerasimenko by the Rosetta/MIDAS instrument. Depending on the particular dust
particle investigated, they found median values of the constituent grain
diameters between $\sim 100$ nm and $\sim 400$ nm. In contrast to the
collected samples of our own Solar System, dust grains in the interstellar
medium (ISM) seem not to possess a single characteristic size value. ISM
grains rather follow a power-law size-frequency distribution, with the
smallest grains in the 1 nm size range and the largest ones having diameters
of $\sim 500$ nm (_Mathis et al._ , 1977). The slope of the power law is such
that while most of the cross section resides in the smallest grains, most of
the mass is in the largest grains. Thus, the largest effect on the mass growth
will be from collisions between large grains; an appropriate range over which
to investigate the collisional growth is $0.1-1~{}\mathrm{\mu m}$.

Laboratory experiments have shown that initially, dust and ice grains in PPDs
generally stick upon mutual collisions (_Blum and Wurm_ , 2008; _Gundlach and
Blum_ , 2015), because their sizes are small (typically $\lesssim
1~{}\mathrm{\mu m}$, see above) and their collision speeds are low
(_Weidenschilling_ , 1977, 1980). Two $1~{}\mathrm{\mu m}$ solid silicate
grains stick to each other when colliding with speeds $\lesssim
1~{}\mathrm{m~{}s^{-1}}$ (_Blum and Wurm_ , 2008), whereas two
$1~{}\mathrm{\mu m}$ solid water-ice grains stick for collision velocities
$\lesssim 10~{}\mathrm{m~{}s^{-1}}$ (_Gundlach and Blum_ , 2015). The latter
result is somewhat surprising, because the surface energy of water ice at low
temperatures does not exceed that of silicates (_Gundlach et al._ , 2018;
_Musiolik and Wurm_ , 2019). Obviously, the surface energy is not the only
material parameter with relevance to the collision and sticking behavior.
_Krijt et al._ (2013) showed that the yield strength and the viscous
relaxation time are two additional parameters that determine the collisional
outcome. Recently, _Arakawa and Krijt_ (2021) showed how the earlier and
latest experimental data on the sticking threshold of $\mathrm{H_{2}O}$-ice,
$\mathrm{CO_{2}}$-ice and $\mathrm{SiO_{2}}$ grains can be reconciled with the
model by _Krijt et al._ (2013). From the various experimental and theoretical
approaches it became clear that the smaller the grains are, the higher their
threshold velocity for sticking is (_Blum and Wurm_ , 2008). A recent study
also showed that the static cohesion (a measure for the surface energy)
between grains of carbon-bearing (mostly organic) substances is heavily
material dependent, with the highest values for graphite and paraffin and the
smallest for humic acid and brown coal; the ratio in static cohesion between
these two groups is as large as a factor $1,000$ (_Bischoff et al._ , 2020).

Figure 1 shows examples of a state-of-the-art collision model for dust and ice
grains in PPDs, based on laboratory collision experiments (_Kothe_ , 2016). In
the top row, the effect of monomer-grain size on the collision outcome is
presented by a comparison between $\mathrm{SiO_{2}}$ grains with
$1.5~{}\mathrm{\mu m}$ and $0.1~{}\mathrm{\mu m}$ diameter, respectively.
Green colors denote the regions in which growth, either by direct sticking or
by mass transfer in collisions in which one of the aggregates fragments,
dominates. In the yellow (bouncing), orange (erosion, abrasion, cratering) and
red regions (fragmentation), growth is not possible. It is important to note
that the collision outcome not only depends on the sizes of the two colliding
aggregates, but also on their collision speed, which is in the example shown
determined for a minimum-mass solar nebula model (dashed contours denote
collision velocities in m/s). From this comparison, it can be seen that for
the more cohesive small grains, potential routes to planetesimal formation are
present, namely by mass transfer in collisions between large bodies and
smaller dust aggregates. However, these growth avenues are surrounded by
orange and red regions so that an easy assessment of whether growth or
destruction dominates is not possible. For large monomer-dust grains, the
pathway to planetesimals is almost absent. Further out in the PPD, where water
ice may be the dominant dust-grain material, is similar to the
$0.1~{}\mathrm{\mu m}$ $\mathrm{SiO_{2}}$ grains at 1 au, even for ice
particles of $1.5~{}\mathrm{\mu m}$ diameter (bottom row in Figure 1). This is
due to the intrinsic higher collisional stickiness of water-ice particles in
comparison to the $\mathrm{SiO_{2}}$ grains (_Gundlach and Blum_ , 2015).

Figure 1: A dust-aggregate collision model for dust and ice aggregates after
_Kothe_ (2016). Top left: for $\mathrm{SiO_{2}}$ monomer grains of
$1.5~{}\mathrm{\mu m}$ diameter at 1 au; top right: for $\mathrm{SiO_{2}}$
monomer grains of $0.1~{}\mathrm{\mu m}$ diameter at 1 au; bottom left: for
water-ice monomer grains of $1.5~{}\mathrm{\mu m}$ diameter at 3 au; bottom
right: for water-ice monomer grains of $1.5~{}\mathrm{\mu m}$ diameter at 10
au. The underlying PPD model is a minimum-mass solar nebula model. Image
credit: Stefan Kothe.

It should be noted that the collision model by _Kothe_ (2016) does not
explicitly include the possible change in porosity during the growth of dust
agglomerates. However, an earlier model by _Güttler et al._ (2010) showed that
there is an effect of porosity on the collisional outcome. Whenever the
colliding particles are in the sticking regime, this difference is small.
However, for larger aggregate masses and/or collisions speeds, compact (low-
porosity) particles enter the bouncing regions, whereas fluffy (high-porosity)
particles can still grow. The transition between fluffy and compact particles
is usually achieved when the collision energies exceed the rolling-friction
threshold, as determined by _Dominik and Tielens_ (1997). This threshold
energy is required to overcome the cohesional friction against mutual rolling
of the particles over their surfaces. It is apparent that a reliable modeling
of dust growth has to take into account the change in porosity during the
growth process. For example, _Zsom and Dullemond_ (2008) show that initially
the dust aggregates start in extremely fluffy (or even fractal, see below)
configurations, but get compacted in mutual collisions when the collision
energy becomes high enough (_Blum and Wurm_ , 2000; _Weidling et al._ , 2009).

Because of the initially slow collision speeds, sticking will be the dominant
collisional outcome for the initial stages of grain growth in PPDs (see all
four panels in Figure 1). _Dominik and Tielens_ (1997) predicted and _Blum and
Wurm_ (2000) experimentally confirmed that if the collision energy is much
less than the energy for rolling or sliding (the latter is analogous to the
rolling case and refers to the force required to slide two surfaces over each
other against cohesional friction) of the grain contacts, the collisional
adhesion will be of hit-and-stick type, i.e. the grain contact will be frozen
at the position of the first contact. In such a case, the growing aggregates
can be described by a fractal description, $m\propto s^{D_{\mathrm{f}}}$,
between their mass $m$ and their characteristic radius $s$, with the fractal
dimension $D_{\mathrm{f}}$. The role of the fractal dimension can best be
understood if one calculates the average mass density and the collision cross
section for fractal aggregates, respectively. The former can be estimated by
$\rho_{\mathrm{f}}\propto\frac{m}{s^{3}}\propto m^{1-3/D_{\mathrm{f}}}$, the
latter by $\sigma_{\mathrm{f}}\propto s^{2}\propto m^{2/D_{\mathrm{f}}}$. For
$D_{\mathrm{f}}=3$, the aggregates behave like any solid material, i.e.
$\rho_{\mathrm{f}}=\mathrm{const}$ and $\sigma_{\mathrm{f}}\propto m^{2/3}$.
These aggregates are just porous, with the same porosity throughout the
volume. For smaller fractal dimensions, however, e.g. $D_{\mathrm{f}}=2$ (see
Fig. 2 for an example), the situation changes, because now
$\rho_{\mathrm{f}}\propto m^{-1/2}$ and $\sigma_{\mathrm{f}}\propto m$. This
means that the average density decreases with increasing aggregate mass and
also from the center of mass radially outward, if one considers a specific
aggregate. On the other hand, the collision cross section increases in
proportion with mass and, thus, much faster than for higher fractal
dimensions, which can have considerable impact on the growth rate.

Figure 2: Example of a simulated dust aggregates with a fractal dimensions
$D_{\mathrm{f}}=2.0$ consisting of 1,024 monomer grains with a very narrow
size distribution.

The exact value of $D_{\mathrm{f}}$ depends on the details of the growth
process and can reach from $D_{\mathrm{f}}\approx 1.5$ in the Brownian-motion-
dominated growth regime (_Paszun and Dominik_ , 2006) to $D_{\mathrm{f}}=3$
when a large particle drifts through a cloud of small grains (_Blum_ , 2006).
Fractal aggregates with $D_{\mathrm{f}}<2$ possess stopping times or Stokes
numbers that are only weakly mass-dependent (see Eqs. 1 and 2) and thus behave
dynamically very differently than those with larger fractal dimensions. Small
fractal dimensions are equivalent to preferred sticking collisions among
similar-sized aggregates, i.e. the mass-frequency distribution of the growing
aggregates in this case is narrow at any given time. On the other hand, for
the highest $D_{\mathrm{f}}$ values the size distribution may be very wide or
even bi-disperse (_Blum_ , 2006) so that the largest aggregates preferentially
grow by collisions with much smaller particles. Detailed studies on this will
be discussed in the following sub-section.

As the collision energy of protoplanetary dust aggregates increases with
increasing mass of the particles, the criterion for hit-and-stick collisions
will no longer be fulfilled and collisions will lead to the compaction of the
aggregates (_Wada et al._ , 2009; _Güttler et al._ , 2010; _Zsom et al._ ,
2010; _Okuzumi et al._ , 2012). Thus, dust aggregates with $D_{\mathrm{f}}<2$
might be converted into those with $D_{\mathrm{f}}>2$ and ultimately may reach
$D_{\mathrm{f}}=3$. It should be mentioned that the compaction process depends
on monomer-grain size, shape and material.

For even higher collision speeds, the colliding aggregates will either bounce
or even fragment upon collision (_Blum and Wurm_ , 2008). While the former
leads to an outside-in compaction of the aggregates (_Weidling et al._ , 2009;
_Güttler et al._ , 2010), the latter changes the mass-frequency distribution
dramatically (_Güttler et al._ , 2010; _Bukhari Syed et al._ , 2017).
Laboratory studies have shown that the threshold speed for fragmentation of
cm-sized aggregates of silica grains is very close to the sticking threshold
of their constituent grains (_Blum and Wurm_ , 2008; _Gundlach and Blum_ ,
2015; _Bukhari Syed et al._ , 2017), i.e. around $1~{}\mathrm{m~{}s^{-1}}$ for
$\sim 1~{}\mathrm{\mu m}$ $\mathrm{SiO_{2}}$ grains. However, recent
experiments with cm-sized aggregates consisting of micrometer-sized water-ice
particles at low temperatures showed that their fragmentation speed is also at
$1~{}\mathrm{m~{}s^{-1}}$ (_Landeck_ , 2018) and, thus, much smaller than the
sticking threshold of their constituent particles, which is around
$1~{}\mathrm{m~{}s^{-1}}$ at low temperatures (_Gundlach and Blum_ , 2015). In
contrast to the collisional sticking process, fragmentation seems to be
dominated by the surface energy. This is consistent with the findings by
_Gundlach et al._ (2018) that the tensile strength of aggregates consisting of
$\sim 1~{}\mathrm{\mu m}$ $\mathrm{SiO_{2}}$ grains and of $\sim
1~{}\mathrm{\mu m}$ $\mathrm{H_{2}O}$ grains is very similar. In conclusion of
the limited empirical data available, the survival of aggregates in collisions
above a few $\mathrm{m~{}s^{-1}}$ is very unlikely.

In order to understand the evolution of the particle size in PPDs, given these
collisional outcomes, one needs to track all possible grain-grain collisions
and their resulting aggregates and fragments, a formidable task. Numerical
models therefore have to simplify this evolution in one way or another. Most
common approaches fall either into a mass-bin-based direct solution of the
collisional evolution (e.g. _Weidenschilling and Ruzmaikina_ , 1994;
_Dullemond and Dominik_ , 2005; _Brauer et al._ , 2008; _Birnstiel et al._ ,
2010; _Okuzumi et al._ , 2012) or into Monte-Carlo based approaches (e.g.
_Ormel and Spaans_ , 2008; _Zsom and Dullemond_ , 2008), see _Drążkowska et
al._ (2014) for a discussion of these approaches. However, many of those
detailed simulations can be understood in simpler, approximate terms following
_Birnstiel et al._ (2012): even in the most optimistic case, where every
collision is assumed to result in perfect sticking, radial drift can usually
remove particles faster than they form. A very simplified treatment of the
growth time scale can be derived from monodisperse growth (e.g. _Kornet et
al._ , 2004; _Brauer et al._ , 2008; _Birnstiel et al._ , 2012) which gives

$t_{\mathrm{grow}}\sim\frac{1}{\epsilon\,\Omega}.$ (12)

This is apart a factor of a few identical with the collisional time scale. See
_Birnstiel et al._ , 2016 for a derivation with a discussion of the
approximations under which this is derived. This time scale (see the vertical
arrows in Figure 3) is roughly independent of particle size (but see _Powell
et al._ , 2019), while the drift time scale (horizontal arrows in Figure 3)
does depend on particle size. Hence, small particles drift so slowly, that
they mainly grow in situ, while large particles can drift inward before they
substantially continue to grow. Equating these particles gives us a limiting
size, called the “drift limit” that particles approximately can reach before
drift will move them to smaller radii

$a_{\mathrm{drift}}\simeq
0.35\frac{\Sigma_{\mathrm{d}}}{\rho_{\mathrm{s}}}\left(\frac{H}{r}\right)^{-2}\,\left|\frac{\partial\ln
P}{\partial\ln r}\right|^{-1},$ (13)

which is shown as orange line in Figure 3.

The resulting particle sizes are typically below millimeters in the outer disk
but can be far beyond meters in the inner disk. This suggests that the outer
disk may be limited by these dynamical effects, while in the inner disk, the
above mentioned collisional effects (bouncing, erosion, fragmentation) are
more likely to stop particle growth. For those disk regions, a limiting
particle size can be derived if the threshold velocity for bouncing or
fragmentation is compared with the size-dependent collision speeds. For
turbulent velocities (_Ormel and Cuzzi_ , 2007), and e.g. a fragmentation
threshold velocity $v_{\mathrm{frag}}$, the fragmentation limit becomes (e.g.
_Birnstiel et al._ , 2012)

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

which shows a strong dependence on the fragmentation velocity. Similar limits
can be derived for either different sources of relative velocities (e.g. drift
velocity) or different threshold velocities (e.g. bouncing threshold
velocity). Figure 3 only considers fragmentation and displays the
fragmentation barrier as purple line. Figure 3 shows how the simulated dust
size distribution (blue contours, after _Birnstiel et al._ , 2012) is well
contained below either the drift limit in the outer disk (orange line) or the
fragmentation limit (purple line). The arrows display the approximate growth
and drift rates. Particle growth is non-local in mass space. This means any
particle size can interact with any other size and affect many other sizes,
for example in case of a distribution of fragments. Still, the evolution of
the bulk mass can often be understood by calculating trajectories in the
mass/radius space (see black dashed lines in Figure 3): small grains grow but
do not drift significantly which results in a vertical motion in the figure.
At $100\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ the particles have reached the
position marked with a cross, which shows that by that time, particles in the
outer disk have barely grown while particles in the inner disk have since long
grown to large sizes and have been lost towards the inner disk. This inward
radial drift can be stopped if pressure maxima are present (we explain this in
detail in Section 3.1.5 below, but also see e.g., _Pinilla et al._ , 2012): at
those positions, radial drift is stopped (see Eqns. 9 and 10) and dust can
accumulate.

The black-dashed trajectories in Figure 3 show that the particles tend to grow
towards and along the drift limit, and in the inner disk along the
fragmentation limit. In the latter case, continuous fragmentation and
coagulation balance each other, which maintains a constant population of both
small and large particles. Figure 4 shows the particle size distributions at
different positions in the disk. It can be seen that in all cases, most of the
mass resides near the largest particles. This, coupled with the fact that the
velocity also scales with the particle size, means that those largest sizes
carry most of the radial mass flux (_Birnstiel et al._ , 2012). The resulting
maximum particle size $a_{\mathrm{max}}$, determines the rate at which the
dust surface density is transported inward and this explains why the two-
population model of (_Birnstiel et al._ , 2012) reproduces more complex
simulations well.

This inward drifting, growing dust might feed dust pressure traps and/or
accrete onto planets and it determines the local particle sizes that may
eventually form planetesimals (_Lambrechts and Johansen_ , 2014; _Lambrechts
et al._ , 2019, see Section 3).

Figure 3: Simulated dust size distribution after _Birnstiel et al._ (2010)1,
shown as blue contours. Gray arrows show the drift and growth rates. The
growth limits are displayed as brown/orange (drift limit) and purple
(fragmentation limit) lines. Gray lines denote the corresponding Stokes
numbers of the particles. Monodisperse growth/drift trajectories are shown as
black dashed lines, ending at $700\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ and
denoting the position at $100\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$ with a
cross. The vertical lines correspond to the size distributions shown in Figure
4. Figure adapted from _Birnstiel et al._ (2016). Figure 4: Dust size density
distributions (surface density per logarithmic mass) measured at 3, 30, and
$100\text{\,}\mathrm{a}\mathrm{u}$ from the simulation in Figure 3. The
corresponding slope of the MRN size distribution as well as a top-heavier
exponent ($\propto a^{-2.5}$) are shown as references. Figure after _Birnstiel
et al._ (2012).

###  2.4 Condensational growth

Collisional growth is not the only possible mode of grain growth in
protoplanetary disks. Dust particles can also grow by deposition of vapor
species on their surface. The key difference is, however, that condensation
stops once the vapor is consumed, i.e. continuous condensation needs
continuous resupply of condensible vapor. In contrast, collisional growth does
not need to be sustained by adding mass. This means that significant growth by
orders of magnitude in size can be easily achieved by collisions (e.g. _Ormel
et al._ , 2007; _Brauer et al._ , 2008; _Zsom and Dullemond_ , 2008; _Okuzumi
et al._ , 2009; _Birnstiel et al._ , 2010), but via condensation would require
orders of magnitude of resupplied mass. A possible way around this limitation
would be continuously repeating events of sublimation and re-condensation.
However, condensation/deposition is a surface effect and small particles
usually dominate the dust surface area of a size distribution even if they
constitute only a negligible fraction of the total dust mass. This means, that
condensational growth as described can only increase the sizes of small
particles efficiently (_Hirashita_ , 2012). This was also seen in _Stammler et
al._ (2017) where both coagulation and sublimation/re-condensation were
simulated around the CO snow line. More recently, _Ros et al._ (2019) pointed
out that heterogeneous ice nucleation on small silicate grains may be
hindered. This, however could at most offset the problem to larger sizes:
whenever there is a wide distribution of particle sizes, the smallest
particles will grow more efficiently than the largest particles.

So while condensation cannot explain growth from ISM sizes to mm-cm sizes, it
is still thought to be an important effect near the water snow line. As the
most abundant volatile, water sublimation at the water snow line can deposit
significant amounts of water vapor (_Cuzzi and Zahnle_ , 2004; _Ciesla and
Cuzzi_ , 2006; _Gárate et al._ , 2020, e.g.,). Outward diffusion of this vapor
can lead to deposition of water ice on the inward-drifting icy particles. This
can locally increase both the particle size as well as the dust surface
density and lead to conditions supporting the streaming instability
(_Drążkowska et al._ , 2016; _Schoonenberg and Ormel_ , 2017). The exact
location where this takes place with respect to the snow line, depends however
also on the effects of back-reaction of the dust onto the gas density
distribution (_Hyodo et al._ , 2019; _Gárate et al._ , 2020).

###  2.5 Grain Size Distributions

Starting from an idealized mono-disperse distribution of individual dust
grains and assuming a high sticking probability in every occurring collision,
the resulting mass-frequency distribution of the growing dust aggregates
depend on (i) the velocity field among the particles and (ii) the fractal
dimension of the aggregates. If we approximate the mass dependency of the
velocity field by a power law $v(m)=v_{0}\left(m/m_{0}\right)^{\gamma}$, with
$v_{0}$ and $m_{0}$ being the velocity and mass of a monomer grain, and the
collision cross section of the fractal aggregates by
$\sigma=\sigma_{0}\left(m/m_{0}\right)^{2/D_{\mathrm{f}}}$, with $\sigma_{0}$
being the collision cross section of the monomer grains, the approximate
growth equation reads

$\frac{\mathrm{d}\left(m/m_{0}\right)(t)}{\mathrm{d}t}=\frac{1}{\tau}~{}\left(m/m_{0}\right)^{\gamma+(2/{D_{\mathrm{f}}})},$
(15)

with $\tau=1/(n_{0}\sigma_{0}v_{0})$ and $n_{0}$ being the collision timescale
and the initial number density of monomer grains, respectively (_Blum_ ,
2006). For $\gamma+(2/D_{\mathrm{f}})<1$, we are in the regime of orderly
growth, for which the average mass grows with a power law of time and the
mass-frequency distribution at any given time is quasi-mono-disperse (i.e.
narrow). For $\gamma+(2/D_{\mathrm{f}})>1$, the growth is characterized by the
so-called runaway process in which the mass-frequency distribution becomes bi-
disperse and the aggregates in the upper mass peak grow infinitely large in a
finite time.

Thus, the growth behavior depends massively on the two parameters $\gamma$ and
$D_{\mathrm{f}}$. The prior one can obtain values between $\gamma=-1/2$ for
Brownian-motion-dominated growth and $\gamma=+1/3$ for drift-dominated
velocities. The fractal dimension can obtain values between
$D_{\mathrm{f}}=1.5-2$ for Brownian-motion-driven growth and
$D_{\mathrm{f}}=3$ for either the pre-compaction stage (see above) or in the
runaway-growth regime. Unfortunately, these two parameters are not independent
and cannot easily be predicted ab initio. Thus, detailed numerical studies are
required for the prediction of the growth speed, the morphology and the mass-
frequency distribution of the aggregates.

As pointed out in Sect. 2.3, growth stops at the bouncing barrier or at the
fragmentation barrier (_Güttler et al._ , 2010; _Zsom et al._ , 2010). As
these barriers typically occur after the sticking-with-compaction stage has
been reached, the aggregates have $D=3$ and an internal porosity of $\sim
60-70\%$ (_Weidling et al._ , 2009; _Güttler et al._ , 2010; _Zsom et al._ ,
2010). Their final size depends on the size and material of the monomer
dust/ice grains as well as on the location in the PPD and can reach values of
$\sim 1$ cm at 1 au (_Zsom et al._ , 2010) and $\sim 1$ mm at 30-50 au (_Lorek
et al._ , 2018). These aggregates have been termed "pebbles".

While growth can proceed orderly (as a well defined peak), or in a run-away
fashion, the outcome of many growth (and transport) simulations tend to reach
a quasi-steady-state size distribution, the shape of which is determined by
the collisional outcome.

Firstly, if particles grow up to the fragmentation barrier, mass is
effectively redistributed to all sizes smaller than the fragmenting particle.
In that case dust mass is transported to larger sizes by growing collisions
and to smaller sizes by fragmenting or eroding collisions. This quickly leads
to a steady state where continuous growth and fragmentation balance each
other. An example of this can be seen in Figure 4 at
$3\text{\,}\mathrm{a}\mathrm{u}$, which shows slices of the dust distrivbution
of Figure 3. The resulting size distribution can depend on the way the
particles fragment, particularly whether most fragmented mass is in the
smallest or the largest fragments. Analytic and semi-analytic models of these
scenarios were derived in _Birnstiel et al._ (2011). Interestingly, for the
most common regime of turbulent collision velocities and fragment
distributions, the resulting size distribution is very close to $n(a)\propto
a^{-3.5}$ which is the result for a fragmentation cascade, despite the
different physics involved. Particles below roughly a micrometer in size,
however are even more strongly depleted, such that the particle surface area
is not dominated by the smallest particles present, but instead by particles
of about one micron (_Birnstiel et al._ , 2018).

Secondly, if particles do not fragment, but instead bounce off each other,
most particle sizes tend to be around the bouncing barrier, leading to a very
narrow size distribution with strongly reduced amounts of fine grains (_Zsom
et al._ , 2010; _Windmark et al._ , 2012b).

Thirdly, if growth is not halted by a collisional barrier, but instead by the
dynamical removal of larger particles in the drift limit, then again the
particle size distribution becomes “top-heavy” with most mass in a narrow
range around the maximum particle size, since small dust is not efficiently
reproduced and instead swept up by the largest particles. The particle size
distribution was measured to follow approximately $n(a)\propto a^{-2.5}$ in
the drift limited case (see _Birnstiel et al._ , 2012, or Figure 4). The
shallower slopes towards the smallest grains is caused either by radial mixing
from the inner disk, where small grains are produced by fragmenting collisions
(see the $30\text{\,}\mathrm{a}\mathrm{u}$ case in Figure 4 and _Birnstiel et
al._ , 2015). In the $100\text{\,}\mathrm{a}\mathrm{u}$ case, the small grains
are remnants of the initial condition as many have not yet grown to larger
sizes.

###  2.6 Constraints from Observations

In recent years, much progress has been made towards testing the theoretical
expectations above through observations. Traditionally, this was done by
measuring the spectral index at millimeter wavelengths
$F_{\nu}\propto\nu^{\alpha_{\mathrm{mm}}}$, as $\alpha_{\mathrm{mm}}$ (in the
optically thin Rayleigh-Jeans limit) is linked to the grain opacity
$\kappa_{\mathrm{abs}}\propto\nu^{\beta}$ via $\alpha_{\mathrm{mm}}=\beta+2$.
$\beta$ in turn depends on the maximum particle size and to a lesser extent on
the particle size distribution. While particles much smaller than the
wavelength typically present values of $\beta\sim 1.7$, values of $\beta<1$
tend to be only reached if the maximum particle size is larger than the
observational wavelength (e.g. _Ricci et al._ , 2010; _Testi et al._ , 2014,
and references therein).

Initial disk-integrated measurements already indicated low $\beta$ values at
millimeter wavelengths, indicative of millimeter sized or larger particles,
which are in agreement with the predicted maximum particle sizes of
theoretical models. A large caveat, however, is that the lifetime of particles
of such sizes, should however be much shorter than the disk age. Theoretical
models produced the right sizes, but only for a short time of a few
${10}^{5}\text{\,}\mathrm{y}\mathrm{r}$ (_Brauer et al._ , 2007, 2008;
_Birnstiel et al._ , 2011). One possible solution, already envisioned by
_Whipple_ (1972) was the introduction of pressure perturbations that have a
vanishing pressure gradient, thus locally halting radial drift. Strong enough
pressure perturbations could locally trap the large particles long enough, as
required by observations (_Pinilla et al._ , 2012).

The theoretically predicted size sorting of larger grains at smaller radii
(_Birnstiel et al._ , 2012) was also confirmed by modeling the radially
resolved spectral indices in _Pérez et al._ (2012); _Tazzari et al._ (2016).
However, it was not until ALMA’s large baseline observations became available
that sub-structures were imaged _ALMA Partnership et al._ (2015); _Andrews et
al._ (2016); _Huang et al._ (2018) that seemed to agree with the idea of
pressure traps. These observations also revived the idea of _Ricci et al._
(2012) that a significant fraction of the flux could be optically thick – then
the observed low spectral index is not caused by large grains, but instead by
optically thick emission (_Tripathi et al._ , 2017).

ALMA has also opened the way to using polarization to measure particle sizes
via self-scattering (_Kataoka et al._ , 2015), which indicated particles not
much larger than about $100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, a
clear conflict with spectral-index based measurements that indicated
millimeter or centimeter sized particles. While this conflict continues to
linger (see _Lin et al._ , 2020, for a possible solution), the width of
several of the imaged rings of DSHARP (_Andrews et al._ , 2018; _Huang et al._
, 2018) is so narrow that dust trapping, i.e. $\mathrm{St}\gtrsim\alpha$ is
most likely required (_Dullemond et al._ , 2018), putting a constraint on the
particle sizes. Furthermore, these rings might be providing the right
conditions for planetesimal formation via the streaming instability and
gravitational collapse (see Section 3.1.3), as also argued for by _Stammler et
al._ (2019) as this would naturally lead to the observed optical depths of
around 0.4 in all of the well resolved sub-structures. The formation of
planetesimals in pressure traps is briefly discussed further in Section 3.1.5.

##  3 DUST TO PLANETESIMALS

Having described the formation and evolution as well as properties of small
solids in PPDs, we now turn to a description of how pebbles transition to the
next phase of planet formation: planetesimals. Currently, there are two
schools of thought for how planetesimals form. The first relies on pebbles
being concentrated somehow to the point where their mutual gravitational
attraction overpower competing effects, such as stellar tides and turbulent
diffusion. As we will see, this approach has provided the largest number of
potential mechanisms, which we refer to as “gravitational collapse models".

The second school of thought takes the opposing view of bottom-up growth.
These models, which we refer to as “coagulation models", rely on the
imperfection of the the various growth barriers described above; as described
further below, there are routes toward sticking mm–cm size pebbles together to
further grow in size, such as the collision-induced mass transfer from smaller
aggregates to larger bodies. Thus, coagulation models describe a route to
planetesimal formation through purely collisional growth.

Most of the literature to date has focused on the former set of models, and
consequently, we devote most of the remainder of this chapter to their
description. However, in Section 3.2, we describe in detail the coagulation
models.

###  3.1 Gravitational Collapse Models

As mentioned above, the gravitational collapse models all provide a route for
pebbles to be concentrated sufficiently well that their mass density
$\rho_{\mathrm{d}}$ exceeds that which is required to counteract the
destructive tidal forces from the central star, the so-called Roche density,

$\rho_{\rm R}\equiv 3.5\frac{M_{*}}{r^{3}}$ (16)

where $M_{*}$ is the mass of the central star, and $r$ is the radial distance
of the self-gravitating object from the star (for a derivation see _Armitage_
2007).

In the remainder of this section, we describe a number of different mechanisms
that all provide a route toward $\rho_{\mathrm{d}}\gtrsim\rho_{\rm R}$
locally, and thus for planetesimals to be born.

####  3.1.1 The Goldreich-Ward Mechanism

The first such model put forth, and in some ways, the simplest one is that of
a gravitationally unstable collisionless particle disk. This mechanism, often
referred to as the Goldreich-Ward mechanism (_Goldreich and Ward_ ,
1973)444Similar work can be found in (_Safronov_ , 1972)., is no longer
thought to work for reasons we will address shortly. However, we include it
here as a starting point for the discussion of other instabilities, especially
as some of the ideas described here carry over to other mechanisms.

Following _Binney and Tremaine_ (1987), our starting point with this
discussion will be the consideration of a razor thin gaseous disk that extends
infinitely, has a constant surface density $\Sigma$, has sound speed
$c_{\mathrm{s}}$, and is rotating uniformly at a rate $\Omega$. By introducing
small perturbations of the form ${\rm exp}\left[{i\mkern
1.0mu}(k_{r}r+k_{\phi}\phi-\omega t)\right]$ (where $\omega$ is the temporal
frequency of the perturbed wave mode and $k_{r}$ and $k_{\phi}$ are the wave
vector components in the radial and azimuthal directions, respectively) to the
governing equations and keeping only leading order (linear) terms (i.e.,
carrying out a “linear perturbation analysis")555We encourage the reader to
read _Binney and Tremaine_ (1987) for more details., one can show that the
gravitational stability of this system is characterized by the so-called
(_Toomre_ , 1964) $Q$ parameter,

$Q_{\rm g}\equiv\frac{c_{\mathrm{s}}\Omega}{\pi G\Sigma}.$ (17)

where the subscript on $Q$ denotes that it is relevant to the disk gas.

This key parameter quantifies the ratio of stabilizing to destabilizing
parameters. More specifically, the gas disk is stabilized on small scales by
thermal motions, hence the presence of $c_{\mathrm{s}}$ in the numerator,
whereas on larger scales, shear acts to stabilize the fluid, as shown by
$\Omega$ (also in the numerator).666This dependence of stability on scales can
more directly be seen by examining the dispersion relation, as shown in
_Binney and Tremaine_ (1987). For more massive disks (i.e., larger $\Sigma$ in
the denominator), the effect of the stabilizing parameters is diminished and
the disk is more subject to gravitational instability.

More quantitatively, if $Q_{\rm g}<Q_{\rm g,crit}$, the system is
gravitationally unstable, and for the setup considered here $Q_{\rm
g,crit}=0.5$ (see _Binney and Tremaine_ 1987). However, gas is a collisional
system, which is different than our collisionless solid particle disk. Yet,
even in this system, $Q$ remains a very useful quantity. If we define $\sigma$
as the one-dimensional velocity dispersion of particles (e.g., due to
turbulent stirring by the gas) and replace $\Sigma$ with $\Sigma_{\rm p}$, the
solid particle surface density, our “particle" Toomre parameter becomes

$Q_{\rm p}\equiv\frac{\sigma\Omega}{\pi G\Sigma_{\rm p}}.$ (18)

The precise value of the critical $Q$ for a differentially rotating Keplerian
disk of particles (the system under consideration here), $Q_{\rm p,crit}$,
will be different than the gaseous disk considered above. However, as in the
gas case, $Q_{\rm p,crit}\approx 1$ and $Q_{\rm p}\lesssim 1$ implies
gravitational instability of the particle layer.

This is the crux of the Goldreich-Ward mechanism; for sufficiently large
$\Sigma_{\rm p}$ and sufficiently low $\sigma$, $Q<1$ and the particle layer
fragments into bound objects.

However, how viable is this mechanism given reasonable assumptions for disk
parameters? Assuming the minimum mass solar nebula model (which is arguably
not the most accurate model, but serves our purposes for now; _Hayashi_ 1981),
the gas surface density at 1 AU is $\Sigma=1700\mathrm{g/cm^{2}}$. With a
standard dust-to-gas ratio of 1%, $\Sigma_{\rm p}\approx 10\mathrm{g/cm^{2}}$
at this location. Further assuming a solar mass star at the center of the
disk, $Q=1$ corresponds to $\sigma\approx 10~{}\mathrm{cm/s}$. Since $H_{\rm
p}=\sigma/\Omega$ and $c_{\mathrm{s}}\sim 10^{5}~{}\mathrm{cm/s}$, we find
that for gravitational instability, the ratio of particle to gas scale heights
is $H_{\rm p}/H_{\rm g}\sim 10^{-4}$.

For cm-sized particles at 1 AU, $\mathrm{St}\sim 10^{-3}$. Using Eq. 7, in
order to maintain such a thin particle layer, $\alpha\sim 10^{-11}$ or less,
whereas even very weak turbulent stirring often leads to a minimum of
$\alpha\sim 10^{-6}$–$10^{-5}$! At this point, the problem with the Goldreich-
Ward mechanism should be clear; turbulence produced by any number of
magnetohydrodynamic or purely hydrodynamic instabilities (see _Balbus and
Hawley_ 1998 and _Lyra and Umurhan_ 2019), would far exceed this tiny number.
Even absent these processes, the velocity shear induced by the particles
pushing back on the gas near the mid-plane and the slower velocity away from
the mid-plane (again, due to the radial pressure gradient) leads to the
Kelvin-Helmholtz instability, which easily produces turbulence that prevents
such a thin particle layer from forming (_Cuzzi et al._ , 1993).

####  3.1.2 Secular Gravitational Instability

While the Goldreich-Ward mechanism is very unlikely to produce planetesimals,
there is a related phenomenon that acts similarly and can get around this thin
layer problem. This process, known as the secular gravitational instability
(SGI) (_Ward_ , 1976) comes about by considering the gravitational instability
problem we just described, but with gas drag included.

The physical mechanism of this process is as follows (see also _Goodman and
Pindor_ 2000). Upon a steady state particle layer with surface density
$\Sigma_{\rm p}$, we impose an axisymmetric particle overdensity. The width of
this overdensity $\delta r$ satisfies $H_{\mathrm{p}}\ll\delta r\ll
H_{\mathrm{g}}$. Due to the effects of particle self-gravity, this overdensity
will contract. However, as this contraction proceeds, the outer edge will move
inward but will maintain its original angular momentum (i.e., there is no
torque introduced by the self-gravity), thus particles increase their orbital
velocity and resist the pull of the self-gravity. Up until now, what we have
described is the stability of the system to pure gravitational instability,
which can be overtaken if the mass of the overdensity is sufficiently strong.
However, another route to overpower this stability is gas drag. Since $\delta
r\ll H_{\mathrm{g}}$, the gas’s orbital speed is much less affected; thus, as
the solid particles move inward, they feel an increased headwind, which
removes their angular momentum and allows self-gravity to win and produce an
even larger overdensity; from here, the whole process accelerates. Similar
arguments can be applied to the inner part of the overdensity, and the result
is run away growth to eventually produce gravitationally bound objects
(though, there are caveats as we will see shortly)

A number of works have examined different aspects of this mechanism, ranging
from its linear growth (_Ward_ , 2000; _Youdin_ , 2005) to the behavior in the
fully non-linear (i.e., no small perturbations assumed) regime (e.g.,
_Pierens_ 2021). _Takahashi and Inutsuka_ (2014) and _Takahashi and Inutsuka_
(2016) carried out a linear analysis of the SGI but with particle feedback on
the gas included, finding that this slightly modified SGI could be responsible
for producing dust rings in protoplanetary systems, such as those observed in
HL Tau (_ALMA Partnership et al._ , 2015).

Other works further explored this linear regime. _Youdin_ (2011) and
_Michikoshi et al._ (2012) found that for small ($\leq$ mm) particles, radial
drift timescales can limit the effectiveness of the SGI to regions with weak
turbulent diffusion, namely $\alpha\lesssim 10^{-4}$. _Shadmehri et al._
(2019) found that non-axisymmetric modes of the SGI are more robust than
axisymmetric modes, though are still limited by turbulent diffusion and are
absent for $\alpha\geq 10^{-3}.$

More recently, a number of simulations of the non-linear SGI have been carried
out. _Tominaga et al._ (2018) ran 2D global simulations in the $r\phi$ plane
and found that indeed the SGI can form dust rings, as was hypothesized by
earlier work. More recently, _Pierens_ (2021) carried out similar simulations
and found that such dust rings can trigger the Rossby Wave Instability (RWI;
_Lovelace et al._ 1999), which then gives rise to gravitationally bound
clusters of particles via trapping in the RWI-induced vortices (see Section
3.1.5 for more details on vortices trapping pebbles). _Tominaga et al._ (2020)
followed up on their previous work by including radial drift and turbulent
diffusion in their simulations, though in the more recent work, they assumed
axisymmetry (i.e., their simulations were 1D); they found that the SGI
maintained its robustness in the nonlinear regime, even with these additional
effects included.

Finally, _Abod et al._ (2019) carried out a number of local, shearing box
simulations designed to study the streaming instability (see below), but one
of their simulations produced gravitationally bound clumps in the presence of
gas drag but with no radial drift (which requires that the clumps form via the
SGI). The formed clumps had a initial mass function very similar to the
streaming instability calculations, though with a slightly steeper power law
towards low masses. Furthermore, the maximum clumps that formed had masses
similar that produced by the streaming instability (see below) and are thus
roughly consistent with (though actually larger than) the largest Solar System
planetesimals. However, as pointed out in _Abod et al._ (2019), modeling the
SGI in the setup that they used posed challenges, and as such, caution is
warranted in interrupting these results.

Ultimately, studies of the SGI, while promising, are in their infancy and more
work is needed to test the viability of this mechanism, particularly compared
to others such as the streaming instability. Indeed, since streaming-
instability clumping is now known to be triggered at relatively low dust
abundances (see below and _Li and Youdin_ 2021) more work is needed to
understand whether SGI could even be triggered before streaming instability
clumping takes over. If the SGI can past these tests of robustness, then the
next stage will be to quantify the properties of any planetesimals formed via
this process.

####  3.1.3 Streaming Instability

The streaming instability (hereafter SI) is the most well-studied of the
mechanisms to produce gravitational collapse. As such, we will devote
considerable attention to its description and the implications for
planetesimal formation.

We use the term “streaming instability" (or “SI") here to refer to a broad
class of behaviors that have been observed in simulations that includes at the
minimum two key ingredients: a background pressure gradient (and thus radial
drift; see below) and momentum feedback from solids to the gas. If these
criterion are included, both linear analyses (e.g., _Youdin and Goodman_ ,
2005) and non-linear simulations (e.g., _Johansen and Youdin_ , 2007) alike
have demonstrated the unstable growth of perturbations that ultimately, as
seen in the non-linear simulations, lead to clumping of solids.

However, as we describe in detail below, there has been considerable follow up
work to this, which has expanded our understanding of the behavior of this
coupled solid-gas system both in the linear and nonlinear state with
additional affects, such as vertical gravity, turbulence, and particle self-
gravity. Some of these additional effects have led to the emergence of new
routes to instability, such as the vertical shearing streaming instability of
_Lin_ (2021). We touch upon some of these extensions below, but for simplicity
(and to be consistent with the literature), we collectively refer to these
processes as the SI, deviating only when appropriate.

Figure 5: Basic mechanism for the streaming instability. A small dust
overdensity (Panel A) brings gas inwards via drag. This inward gas motion is
deflected azimuthally via the Coriolis force (Panel B). Dust already present
downstream of the overdensity then feels an azimuthal boost from the deflected
gas via drag, which then causes this downstream dust to move to larger radii
(as it now has increased angular momentum), adding to the original dust
overdensity. Note that this entire process is in a frame corresponding to the
average radial drift of the dust. (Panel C). Modified from _Squire and
Hopkins_ (2020).

Simplest Model: Linear, Unstratified, & No Self-gravity
Our starting point is the simplest model consisting of a uniform distribution
of solid particles in a uniform gas, both subject to the radial component of
the central star’s gravity (but not vertical!) and with angular momentum
exchange between solids and gas present. By carrying out a linear perturbation
analysis on the governing equations, _Youdin and Goodman_ (2005) and _Youdin
and Johansen_ (2007) showed that this system is unstable, with growth rates
and unstable wavelengths that depend on $\epsilon,\mathrm{St}$, and $\eta$ (as
defined previously).

While a physically intuitive explanation for this instability remained elusive
for many years, recent work by _Squire and Hopkins_ (2018, 2020) has provided
some clarity. This work led to the realization that in the limit of
$\epsilon<1$, the SI is actually one of a broader class of “resonant drag
instabilies”. These instabilities arise when the drift velocity of solids
equal the phase velocity of a specific wave mode in the gas. Such a
configuration equates to a resonance in the frame of the solids. In the case
of the SI, this resonance occurs between radial drift and epicyclic
oscillations of the gas, which lead to slow growth of the dust concentration
over time.

When $\epsilon>1$, on the other hand, the feedback from the dust overwhelms
the gas, which accelerates the dust concentration drastically. Both regimes,
however, operate on the same basic idea, as shown in Fig. 5: a dust
overdensity drags gas inwards (Panel A in the figure). This inward motion is
deflected into the azimuthal direction via the Coriolis force (Panel B). Dust
already present downstream of the overdensity then feels a boost from the
deflected gas via drag, which then causes this downstream dust to move outward
radially, adding to the original dust overdensity (Panel C). There are a
number of subtleties associated with this picture, including the difference
between the $\epsilon<1$ and $\epsilon>1$ regimes, which are described in
detail in _Squire and Hopkins_ (2020). However, for the purposes of this
chapter, this picture captures the essence of the SI and demonstrates how at
its core, it is a concentration-inducing instability.777While the original
form of the SI, as derived by _Youdin and Goodman_ (2005) is 2D in the radial-
vertical plane, this diagram sacrifices information about the vertical
dimension and includes the azimuthal dimension in order to provide a more
intuitive understanding of the mechanism. Furthermore, the mechanism as
depicted here applies to the axisymmetric case as well.

Nonlinear Simulations
The linear phase of the SI elucidated, we now turn to a more in-depth
discussion of the non-linear regime, probed with numerical simulations. While
unstratified simulations have been run (_Johansen and Youdin_ , 2007; _Bai and
Stone_ , 2010b), they have largely been used to test the employed code’s
ability to reproduce the linear growth phase.888For brevity, we will neglect a
description of unstratified simulations beyond reproducing linear results.

Most SI simulations have included the vertical component to the central star’s
gravity (the so-called “stratified” setup). Furthermore, with a few exceptions
(see _Kowalik et al._ 2013; _Mignone et al._ 2019; _Schäfer et al._ 2020;
_Flock and Mignone_ 2021), these simulations have all employed the local,
shearing box approximation. The shearing box treats a local, co-rotating patch
of the disk of sufficiently small size to expand the relevant equations into
Cartesian coordinates (see e.g., _Hawley et al._ 1995). This approach is
generally done when small scales ($\ll H_{\mathrm{g}}$) must be adequately
resolved.

Within the context of the SI, such simulations have been done both in 2D
(assuming axisymmetry) and full 3D. The primary goal of these simulations has
been to study the properties of and parameters responsible for “strong
clumping”, defined here as clumping sufficiently strong so as to locally reach
the Roche density. Perhaps the most impactful of such studies have been a
series of papers started by _Carrera et al._ (2015); the authors of that work
ran a series of 2D axisymmetric SI simulations finding an approximate
“U-shaped" region of the $\mathrm{St}$ vs. $Z$ (where $Z\equiv\Sigma_{\rm
p}/\Sigma$ is the dust concentration) parameter space where strong clumping
occurs. To cover such a large parameter space, for each $\mathrm{St}$, $Z$ was
slowly increased by removing gas from the domain.

_Yang et al._ (2017) expanded on this work to find that the SI can produce
strong clumping at even smaller $\mathrm{St}$ than was determined by _Carrera
et al._ (2015) if the numerical resolution is high enough and the simulation
allowed to run for sufficiently long time. They also verified that for a
subset of their parameters, the results are similar between 2D and 3D.

More recently, _Li and Youdin_ (2021) re-examined the problem again, over a
larger range of $\mathrm{St}$ values than explored by _Yang et al._ (2017) but
also running 2D and 3D calculations for very long integration times. They
found that the strong clumping criterion is even less restrictive than
originally thought, allowing such clumping to occur for $Z=0.004$, well below
the standard value assumed to be inherited from the ISM. Figure 6 summarizes
the results from these three works and indicates our current understanding of
the parameters for which strong clumping is allowed.

Beyond $\mathrm{St}$ and $\epsilon$, the latter of which can be related to $Z$
via $Z=\epsilon\left(H_{g}/H_{p}\right)$, the third parameter that controls
the linear SI growth is the pressure gradient associated with stellar
illumination.999The gas pressure decreases with radius to partially balance
the disk against the radial component of gravity. The remainder of the balance
comes from azimuthal flow (i.e., centrifugal balance). That pressure aids
rotation in balancing gravity is the reason for the gas to be sub-Keplerian,
albeit slightly. This pressure gradient can be written in dimensionless form
as

$\Pi\equiv\frac{\Delta
v}{c_{\mathrm{s}}}=-\frac{1}{2}\left(\frac{c_{\mathrm{s}}}{v_{\rm
k}}\right)\frac{\partial\ln P}{\partial\ln r}$ (19)

where $v_{\rm k}$ is the Keplerian velocity, $P$ is the gas pressure, and
$\Delta v$ is the headwind experienced by solid particles

$\Delta v\equiv v_{\rm k}-u_{\phi}=\eta v_{\rm k}$ (20)

_Carrera et al._ (2015) found that lower $\Pi$ decreased the threshold for
strong clumping, consistent with earlier work by _Bai and Stone_ (2010c).
While this third parameter may make further parameter studies of the SI
expensive, recent results (_Sekiya and Onishi_ , 2018) have demonstrated that
$Z/\Pi$ is actually the more relevant parameter for determining the clumping
properties of the SI. That is, if one changes $Z$ and $\Pi$ separately but
their ratio remain constant, the SI clumping appears to be very similar
(_Sekiya and Onishi_ , 2018).

This critical result allows us to expand the parameter surveys described above
to define the critical concentration parameter $Z_{\rm crit}$ for any value of
$\Pi$ as (see equation 9 in _Li and Youdin_ 2021)

$Z_{\mathrm{crit}}\left(\Pi\right)=\frac{\Pi}{0.05}Z_{\rm
crit}\left(\Pi=0.05\right).$ (21)

This result has a clear physical interpretation when considering the result
that larger $\Pi$ is associated with more turbulence induced by the SI
(_Carrera et al._ , 2015; _Abod et al._ , 2019). Higher turbulence equates to
a larger $H_{p}$, which for a given $Z$ will decrease $\epsilon$ (again
because $Z=\epsilon H_{g}/H_{p}$), which as has been shown from unstratified
simulations is a control parameter for the fully non-linear SI.

This result also led _Li and Youdin_ (2021) to write an equation for the
critical $\epsilon$ responsible for strong clumping by combining the _Sekiya
and Onishi_ (2018) result with their simulation suite. _Li and Youdin_ (2021)
also included the effect of externally driven turbulence (e.g., from the MRI)
in their equation (though external turbulence was not included in their
simulations through the traditional $\alpha$ parameter (defined above), as
this also sets $H_{p}$ and thus $\epsilon$. While we will return to the effect
of externally driven turbulence shortly, it is worth writing this critical
$\epsilon$ here as it is our most current understanding of the condition for
strong clumping due to the SI. From _Li and Youdin_ (2021),

${\rm log_{10}}\,\epsilon_{\mathrm{crit}}\simeq A^{\prime}\left({\rm
log_{10}}\,\mathrm{St}\right)^{2}+B^{\prime}\left({\rm
log_{10}}\,\mathrm{St}\right)+C^{\prime}$ (22)

where

$\begin{cases}A^{\prime}=0,B^{\prime}=0,C^{\prime}={\rm
log_{10}(2.5)}&\text{if \, }\mathrm{St}<0.015\\\
A^{\prime}=0.48,B^{\prime}=0.87,C^{\prime}=-0.11&\text{if \,
}\mathrm{St}>0.015\end{cases}$

This can then be included into a new equation for $Z_{\mathrm{crit}}$ that
takes into account the effect of $\Pi$ (see Equation 21) and externally driven
turbulence parameterized via $\alpha$,

$Z_{\mathrm{crit}}\simeq\epsilon_{\mathrm{crit}}(\mathrm{St})\sqrt{\left(\frac{\Pi}{5}\right)^{2}+\frac{\alpha}{(\alpha+\mathrm{St})}}$
(23)

While the above set of equations represents our current best understanding of
where in physical parameter space, strong clumping occurs and planetesimal
formation is likely, numerical effects must also be quantified when employing
simulations. _Yang and Johansen_ (2014) studied the effect of domain size,
which is particularly important as the small scales over which the SI operate
($\ll H_{\rm g}$) restrict simulations to small domain sizes, typically
$\approx 0.2H_{\rm g}$ on a side. The authors found an average separation of
$\sim 0.2H_{\rm g}$ between SI-induced clumps, suggesting that at least this
domain size should be used in SI simulations. The authors also found that the
maximum degree of clumping (measured as the maximum particle density
throughout the simulation domain) increases with increasing resolution,
consistent with prior work (_Johansen et al._ , 2007, 2012; _Bai and Stone_ ,
2010b), though _Yang and Johansen_ (2014) did find tentative evidence for
convergence in this quantity at $160$–$320$ grid zones per $H_{g}$.

_Li et al._ (2018) followed up on this work, re-examining the effect of box
size and resolution, while also exploring different vertical boundary
conditions. Given the small domain size typically used in such simulations,
previously used boundary conditions were either periodic or reflecting. _Li et
al._ (2018) examined the effect of a modified open boundary (first developed
for MHD simulations in _Simon et al._ 2011 and then modified by _Li et al._
2018), finding that such a boundary led to better convergence of results with
box size, a smaller $H_{p}$ (due to the absence of reflecting or periodic
vertical gas flows), but no appreciable change in the maximum particle
density.

More recently, _Li and Youdin_ (2021) found that larger domains in the radial
direction may be necessary for strong clumping when the relevant parameters
lie close to the clumping boundary (the solid black line in Fig. 6).
Furthermore, _Rucska and Wadsley_ (2021) found evidence that scales larger
than $0.2H_{\rm g}$ may be important in determining the properties of
planetesimal born from the streaming instability.

Figure 6: Regions in $\mathrm{St}$–$Z$ space for which strong clumping does or
does not occur, based on the work by _Carrera et al._ (2015); _Yang et al._
(2017); _Li and Youdin_ (2021) and assuming $\Pi=0.05$. The pink region
corresponds to no strong clumping, whereas both green regions correspond to
strong clumping. The darker green region, outlined by the blue dashed line
indicates the work by _Carrera et al._ (2015) in which strong clumping was
more rapid due to the continuous removal of gas, amounting to a 50% reduction
of the remaining gas every 50 orbits of the local shearing box.

Clearly, a picture of what parameters lead to strong clumping is emerging.
However, there still remains more work to be done, exploring both physical
parameters and numerical effects to fully flesh out the nature of strong
clumping. For example, a major unanswered question is: why does strong
clumping occur for the parameter choices that it does? As pointed out by _Li
and Youdin_ (2021), large linear growth rates do not necessarily lead to
strong clumping. Not only does this motivate a deeper understanding of how
strong clumping works, but it also warrants caution in over-interpreting the
results of linear analysis. Furthermore, improving upon other conditions used
in previous work, such as assuming axisymmetry for most of the $\mathrm{St}$,
$Z$ parameter exploration, as well as neglecting particle self-gravity and its
influence on strong clumping, may prove to change (yet again) the parameters
for which strong clumping occurs.

Figure 7: Log of the pebble surface density (relative to the average value) in
a series of SI simulation snapshots. The simulation domain is a small patch of
a disk, with the vertical axis being azimuthal and the horizontal axis radial.
As the SI proceeds, pebble over-densities begin to form, ultimately creating
azimuthally elongated filaments. Eventually (right-most panel) the mass
density in these filaments exceeds Roche, leading to gravitationally bound
clouds of pebbles. Modified from _Abod et al._ (2019). $\Omega$ as shown on
the figure refers to the rotation rate at the center of the disk patch.

Extensions: Turbulence, Multiple Particle Species
There have been a number of extensions applied to the SI problem, in the
linear and nonlinear regimes alike. For instance, externally driven turbulence
has been shown to significantly reduce the growth rate of the SI in the linear
regime (_Umurhan et al._ 2020; _Chen and Lin_ 2020 and see preliminary
comments on the matter in _Youdin and Goodman_ 2005). However, as we just
emphasized (and also mentioned in other works; see _Li and Youdin_ 2021), one
should be cautious in interpreting the linear growth rates as having a direct
influence on the degree of strong clumping.

Indeed results from nonlinear simulations appear mixed. Simulations by _Gole
et al._ (2020), for which turbulence was injected into a domain of
sufficiently high resolution and small scale to resolve the SI, shows that
strong clumping (and planetesimal formation; they included particle self-
gravity) can be significantly reduced by turbulence.

A number of other works, on the other hand, (_Johansen et al._ , 2007; _Yang
et al._ , 2018; _Schäfer et al._ , 2020; _Xu and Bai_ , 2022a) included the
turbulence self-consistently (e.g., included magnetic fields to activate the
MRI). These studies found that the sources of turbulence explored led to solid
particle concentration on scales $\sim 1$–10$H_{\mathrm{g}}$, pushing the
local concentration into the streaming-unstable regime. While the obvious
implication is that turbulence may actually aid in planetesimal formation,
many of these studies were carried out at modest resolution. Given that the
rate at which over-densities are diffused via turbulence scales with $k^{2}$,
where $k=2\pi/\lambda$ and $\lambda$ is the typical size scale of an
overdensity, resolving smaller scales may enhance the turbulent diffusion of
small-scale clumping and thus reduce or quench the production of
planetesimals. Ultimately, more work remains to answer the question of whether
turbulence hinders or helps SI-induced planetesimal formation; such studies
are already underway.

Another extension is the inclusion of multiple particle sizes. Of course, a
non monodisperse set of particles naturally result from collisional
growth/fragmentation models, as outlined above (see Sec. 2.5). However, for
simplicity, many previous SI studies have assumed just a single size for all
particles. Recently, a linear analysis carried out by _Krapp et al._ (2019)
has shown that for $\epsilon\lesssim 1$, the fastest growth rate of the SI
decreases with increasing number of discrete particle sizes without any sign
of convergence (i.e., the growth rate approaches zero as the continuum limit
of particle sizes is reached). This surprising result was further explored by
_Zhu and Yang_ (2021) and _Paardekooper et al._ (2020), which found that the
multi-species SI is largely controlled by $\epsilon$ and the maximum
$\mathrm{St}$ of the distribution, $\mathrm{St_{max}}$. Furthermore, for
$\epsilon\gtrsim 1$ and $\mathrm{St_{max}}\gtrsim 1$, the SI growth rates are
of order the dynamical time, regardless of the number of particle sizes.

More encouragingly, non-linear, vertically stratified SI simulations with
multiple particle sizes have also demonstrated the presence of the SI (_Bai
and Stone_ , 2010a; _Schaffer et al._ , 2018) and even the convergence of
particle concentration via the SI with increasing number of particle sizes
(_Schaffer et al._ , 2021). The latter work suggests that their results
converge because the conditions near the disk mid-plane where the SI is active
satisfy the convergent region of parameter space in the linear regime. While
this result is encouraging, these considerations imply that careful tests of
convergence should be carried out when running multiple-species SI
simulations.

Formation of Gravitationally Bound Clouds
Once strong clumping is achieved and the Roche density surpassed,
gravitational forces overpower tidal and turbulent forces (see, e.g., _Gerbig
et al._ 2020) and gravitationally bound clouds form101010We use the term cloud
here to distinguish these diffuse, spheroidal structures from the SI-induced
clumps, which are not necessarily gravitationally bound.. This has been
demonstrated numerically with a number of calculations that include particle-
particle self-gravitational forces.

One such example is shown in Fig. 7; the SI produces elongated (nearly
axisymmetric) filaments that continue to increase in concentration.
Eventually, the Roche density is reached locally, and gravitationally bound
clouds of particles spawn off of the SI-induced filaments. While the figure
shows just one such simulation (from _Abod et al._ 2019), a number of works
have demonstrated that this is a robust process (e.g., _Johansen et al._ ,
2007, 2015; _Simon et al._ , 2016, 2017; _Abod et al._ , 2019; _Li et al._ ,
2019). Furthermore, these works have shown that the largest clouds have masses
on roughly the same order as the largest Solar System planetesimals and dwarf
planets, such as Ceres (though as we will see shortly, this agreement is at
best very approximate).

A more in-depth comparison of these clouds with observations of Solar System
planetesimals is addressed in Section 4. However, at this point, a technical
issue associated with many of the above calculations must be discussed. The
algorithm to solve Poisson’s equation for the particle gravitational
potential, the particle-mesh (PM) method (which has been the only method
employed in calculations of SI-induced planetesimal formation to date), limits
the further collapse of such clouds to a factor of $10^{3}$–$10^{4}$ larger
than the scales of the largest observed planetesimals ($\approx$ 100 km). The
specific details of why this happens are beyond the scope of this chapter.
However, this is a purely numerical limitation, and in some cases, a work-
around has already been achieved. For instance, recent work with the Pencil
code has included “sink particles” (_Johansen et al._ , 2015). Another work
around is an effective zoom-in on individual bound clouds with an improved
treatment of gravity, as we now describe in more detail.111111We should note
that simulations including the entire range of scales from birth of
planetesimals to final collapse are under development.

####  3.1.4 Final Stage of Collapse: Planetesimals are Born

Recently, _Nesvorný et al._ (2021) carried out a series of N-body calculations
with a modified version of Pkdgrav that solves particle gravity with a tree
algorithm (see _Barnes and Hut_ 1986), includes collisions between particles,
and properly accounts for energy dissipation via a new superparticle method
(see _Nesvorný et al._ 2020). Each of these calculations focused on individual
particle clouds extracted from SI calculations. Furthermore, the particle
collisions are expected to be inelastic: they extract the kinetic energy from
the collapsing cloud, which allows it to collapse further towards denser and
denser objects.

Figure 8: Gravitational collapse of a particle cloud. The four snapshots show
different stages of collapse at simulation time $t=0$, 6, 13 and 50 yr (from
left to right). The simulation was performed with the superparticle module of
the PKDGRAV $N$-body code (Nesvorný et al. 2020). The view is projected down
the clump angular momentum axis. Each frame is 450,000 km across. The equal-
size binary in the right-hand panel has primary diameter $D_{1}\approx 110$
km, secondary diameter $D_{2}\approx 78$ km, semimajor axis $a_{\rm b}=10400$
km, and eccentricity $e_{\rm b}=0.19$. The binary represents 98% of the
initial mass and 58% of the initial angular momentum of its parent cloud.

_Nesvorný et al._ (2021) found that the particle clouds rapidly collapse into
short-lived disk structures from which planetesimals form (see Fig. 8). The
planetesimal properties depend on the cloud’s scaled angular momentum,
$l=L/(MR_{\rm H}^{2}\Omega)$, where $L$ and $M$ are the angular momentum and
mass, $R_{\rm H}$ is the Hill radius (i.e., the radius at which the gravity of
the planetary body or planetesimal dominates over that from the Sun), and
$\Omega$ is the orbital frequency. Low-$l$ particle clouds produce tight (or
contact) binaries and single planetesimals. Compact high-$l$ clouds give birth
to binary planetesimals with attributes that closely resemble the equally
sized binaries found in the Kuiper Belt (see Section 4). Furthermore, 50%-100%
of the original mass of a particle cloud is converted to the final largest
planetesimal or planetesimal binary. The remaining $<50$% of the original mass
remains in particles and small planetesimals (_Robinson et al._ , 2020).

We will return shortly to binary production in collapsing particle clouds,
particularly within the context of comparison with observations. However, at
this point, it is worth describing studies of interior planetesimal
structures, as this will also be directly relevant to the examination of
actual Solar System bodies.

Planetesimal Structures - Predictions from Theory
_Wahlberg Jansson and Johansen_ (2014), and the subsequent works of these
authors (see below) considered gravitational collapse of pebble clouds, though
their setup was generic (i.e., not taken from an SI calculation) and had
$L\sim 0$ (i.e., no global rotation). These authors developed statistical
schemes to characterize the effects of collisions on the pebble size
distribution, with different collisional regimes being considered, including
sticking, bouncing, mass transfer, and fragmentation of pebbles (_Güttler et
al._ , 2010). They found that the initial mass of a pebble cloud determines
the interior structure of the resulting planetesimal.

The pebble collision speeds in low-mass clouds are below the threshold for
fragmentation ($\sim 1-10~{}\mathrm{m~{}s^{-1}}$, see Sect. 2.3), forming
pebble-pile planetesimals consisting of the primordial pebbles from the PPD.
Planetesimals above 100 km in radius should consist of mixtures of dust
(pebble fragments) and pebbles which have undergone substantial collisions
with dust and other pebbles. If comet-sized planetesimals form by
gravitational collapse of small pebble clouds, this model would predict that
the recently visited Solar System bodies 67P/Churyumov-Gerasimenko and
Arrokoth (by Rosetta and New Horizons, respectively) are pebble-pile
planetesimals consisting of primordial pebbles from the solar nebula. If,
instead, comets are collisional fragments of much larger planetesimals, the
memory of original pebbles may have been lost.

These results were extended (_Wahlberg Jansson and Johansen_ , 2017) by
allowing the individual pebble sub-clouds to contract at different rates and
including the effect of gas drag (on the contraction speed and energy
dissipation). Their results yield comets that are porous pebble-piles with
particle sizes varying with depth. The interior should consist of primordial
pebbles with a narrow size distribution, yielding higher porosity. The surface
layers should be a mixture of primordial pebbles and pebble fragments. The
detailed predictions depend on the applied fragmentation model (_Wahlberg
Jansson et al._ , 2017).

The physical properties of primordial pebble piles (ignoring pebble
fragmentation during collapse) were calculated by _Skorov and Blum_ (2012) and
_Blum_ (2018) and compared to known properties of comets (see Sect. 4).
_Skorov and Blum_ (2012) calculated the tensile strength of the surface layers
of pebble piles, which is many orders of magnitude lower than that of non-
hierarchical, porous bodies consisting of (sub-)micrometer-sized dust/ice
grains. The reason for the low strength is that the number of grain-grain
contacts, which determines the cohesion of granular matter under low-gravity
conditions, is greatly reduced. This in turn results from the contacts
occurring only in the pebble-pebble contact area, which is much smaller than
the pebble cross section. Other major differences to the non-hierarchical
makeup of planetesimals are the vastly different heat conductivity and gas
permeability, which are important properties to understand the thermal
evolution of planetesimals and the activity of comets (see below).

Due to the hierarchical setup of these kinds of planetesimals (dust/ice grains
– pebbles – planetesimal), a number of characteristics are present:
Low mass density / high porosity. Due to the effect of bouncing after the
pebbles had been formed in the PPD, their internal volume filling factor was
estimated to be in the range $\Phi_{\mathrm{int}}\approx 0.33-0.58$ (_Güttler
et al._ , 2009; _O’Rourke et al._ , 2020). The pebbles themselves form a
pebble pile in the gravitational collapse that has a filling factor of
$\Phi_{\mathrm{pp}}\approx 0.55-0.64$ (_O’Rourke et al._ , 2020). Thus, small
planetesimals ($\lesssim 100$ km) as a whole should possess a volume filling
factor of
$\Phi_{\mathrm{pl}}=\Phi_{\mathrm{int}}\cdot\Phi_{\mathrm{pp}}\approx
0.18-0.37$ and a mass density of
$\rho_{\mathrm{pl}}=\Phi_{\mathrm{pl}}\cdot\rho_{\mathrm{solid}}$, with
$\rho_{\mathrm{solid}}$ being the average mass density of the solid monomer
grains. For large planetesimals ($\gtrsim 100$ km), their self-gravity is so
large that the pebbles collapse and fill the macroscopic void spaces. Thus,
these planetesimals have reduced porosities and increased densities and also
possess other very different physical properties (for details see _Blum_ ,
2018).
Low tensile, compressive and fragmentation strengths. Due to the macroscopic
size of the pebbles, the inner cohesion of planetesimals formed in a gentle
gravitational collapse of a pebble cloud is extremely weak (_Skorov and Blum_
, 2012). For mm- to cm-sized pebbles, tensile strength values for the
planetesimals of $\lesssim 1$ Pa are expected. Similar low values are expected
for the compressive strength (_O’Rourke et al._ , 2020). In contrast, a
homogeneous (i.e. non-hierarchical) packing of microscopic solid grains leads
to tensile strengths $\gtrsim 1$ kPa (_Blum_ , 2018).

For the dynamical evolution of planetesimals, it is essential to consider
their collisional behavior. A characteristic value often used in collisional
evolution studies is the fragmentation strength, which is the collision energy
per unit mass required to yield a largest fragment with half the mass of the
original body. _Katsuragi and Blum_ (2017) and _San Sebastián et al._ (2020)
showed that the fragmentation strength scales with the static tensile
strength. Thus, the energy required for the fragmentation of pebble-pile
planetesimals can be many orders of magnitude smaller than in the non-
hierarchical case (_Krivov et al._ , 2018; _San Sebastián et al._ , 2020).
High gas permeability. The macroscopic void spaces between the pebbles in a
planetesimal that was formed by a gentle gravitational collapse of a pebble
cloud play an important role in vapor transport (see _Lichtenberg et al._ ,
2016, for the influence of porosity on heating and gas transport). Besides
current Solar System bodies exhibiting gas and/or dust activity, gas transport
might be important in the first few million years after their formation if the
abundance of short-lived radiogenic nuclei is high enough to vaporize volatile
(e.g. water ice) or super-volatile (e.g. CO ice) materials (_Lichtenberg et
al._ , 2021; _Golabek and Jutzi_ , 2021; _Malamud et al._ , 2022). Due to a
negative temperature gradient in the radial direction (in this context, radial
means away from the center of the planetesimal), volatiles could be driven
radially outward and could lead to their enhancement close to the surface of
the planetesimal.
Low solid-state and high radiative thermal conductivity. Besides the gas
permeability and the mechanical strengths, the heat conductivity is heavily
dependent on the size of the building blocks of the planetesimal. For a
homogeneous body consisting of microscopic dust or ice grains, the heat
transport occurs through the network of grains. Thus, the heat conductivity
decreases for increasing porosity and increasing grain size. In contrast, a
pebble-pile makeup of the planetesimals may lead to an increased heat
transport, due to the influence of radiative heat transfer, in particular at
higher temperatures and for larger pebble sizes (_Bischoff et al._ , 2021).
The strong temperature dependency of the heat conductivity in the pebble case
can lead to interesting diurnal effects, with high energy transport from the
surface into the sub-surface regions during periods of high-intensity
insolation and low energy transport at night times (_Bischoff et al._ , 2021).

As with the other subsections in Section 3, we will return to
planetesimal/comet structures shortly as we make a comparison with
observations. However, for now, we conclude this section with a brief
discussion of another route to producing planetesimals via gravitational
collapse.

####  3.1.5 Concentration in Gas Features

The final gravitational collapse mechanisms we will discuss are the broad
group of “concentration mechanisms". These mechanisms provide a route
(alternative to the SI) toward concentrating pebbles to the point that their
local density surpasses the Roche density. In what follows, we only briefly
describe these mechanisms, as the literature is considerably less extensive
compared with that of the SI.

The simplest concentration mechanism is an axisymmetric gas pressure
enhancement. These can be induced by any number of things, including
magnetohydrodynamical effects (e.g., _Johansen et al._ , 2009; _Simon and
Armitage_ , 2014; _Bai and Stone_ , 2014; _Béthune et al._ , 2017),
hydrodynamic structures from instabilities (e.g, the vertical shear
instability; _Schäfer et al._ 2020), and planets themselves, which produce
pressure maxima at the outer edge of planet-induced gaps (e.g., _Picogna and
Kley_ , 2015; _Fedele et al._ , 2017; _Hendler et al._ , 2018; _Baruteau et
al._ , 2019; _Pérez et al._ , 2019; _Veronesi et al._ , 2020). Indeed, ALMA
has observed the effect of dust concentration in such rings in a number of
planet-forming disks (_Huang et al._ , 2018), though their origin has not yet
been well constrained by such observations.

The process by which such concentration occurs is as follows. Pebbles on the
outer part of a pressure enhancement drift towards increasing pressure due to
the headwind they feel from the partially pressure supported sub-Keplerian gas
(see Section 2.2). However, on the other side of the pressure enhancement
(i.e., at smaller radii), the opposite is true. In order for the pressure
enhancement to be long lived, the Coriolis force and gas pressure gradient
must balance (this is known as geostrophic balance). Thus, inward of the
enhancement, gas is moving super Keplerian, which provides a tailwind to
pebbles, increasing their angular momentum; the pebbles move towards pressure
maxima and are trapped in the case where ${\rm d}P/{\rm d}r\geq 0$.121212In
the case in which a local maximum is not produced, the inward radial drift of
pebbles continues, but their drift slows down as they enter the pressure
enhancement causing a “traffic jam”.

A related process is the concentration of pebbles in vortices. Gas vortices
are another form of geostrophic balance, but without axisymmetry. In the non-
inertial co-rotating shearing box frame, cyclonic vorticies (that is with
${\mbox{\boldmath{$\nabla$}}}\times{\mbox{\boldmath{$v$}}}>0$), are sheared
out by the anti-cyclonic shear flow. Only anti-cyclonic vortices survive, and
in this configuration a high pressure region exists within the vortex that
balances the Coriolis force; geostrophy is maintained.

As the vortex orbits the central body (at Keplerian velocities), it sweeps up
pebbles that are inwardly drifting from larger radii. These pebbles are
captured in the outer region of the vortex and are swept along with the vortex
streamlines. However, since they experience drag but no pressure gradient,
they move within the vortex at a higher velocity relative to the gas, which
causes them to slowly spiral towards the center of the vortex. An apt analogy
is to think of the vortex as a miniature disk with the role of radial gravity
being replaced by the Coriolis force. Thus the pebbles orbit the center of the
vortex as they would a disk, at a higher velocity than the gas. The resulting
headwind on the pebbles causes them to spiral inward to the center of the
vortex.

Indeed, there is already evidence that dust can become entrapped by these
vortices. For instance, early dust continuum observations with ALMA (e.g.,
_van der Marel et al._ , 2013; _Pérez et al._ , 2014) depicted the presence of
large-scale asymmetric concentrations of dust grains, which has been
interpreted as trapping of these grains in vortices.

A third mechanism, known as turbulent clustering is the concentration of
pebbles due to turbulent gas in planet-forming disks (_Klahr and Henning_ ,
1997; _Cuzzi et al._ , 2008; _Pan et al._ , 2011; _Hopkins_ , 2016; _Hartlep
et al._ , 2017; _Hartlep and Cuzzi_ , 2020), namely through the interaction
between pebbles and turbulent eddies. Within this mechanism, pebbles that are
coupled to the turbulent gas are centrifugally “flung" outward away from the
centers of turbulent vorticies/eddies and are concentrated between the eddies
to the point of gravitational collapse. However, for this to work effectively,
several conditions must be satisfied:

  1. 1.

The pebble clump is of sufficiently high density to resist the destruction by
fluid motion within the turbulent eddy as the clump and eddy interact

  2. 2.

The pebble clump is of sufficiently high density that it can resist disruption
from “ram pressure" (i.e., feeling a headwind as it moves through the gas)

  3. 3.

The pebble clump is of sufficiently low density that the feedback from pebbles
does not destroy the turbulent eddies.

More quantitatively, if one considers pebbles with a specific stopping time
$t_{\rm stop}$, the Stokes number at length-scale $\ell$ is defined
as131313The ratio of stopping time to eddy turnover time is the true
definition of the Stokes number and the quantity St defined in Equation 2 is
the Stokes number at scale $H$ in which $t_{{\rm eddy},H}=\Omega^{-1}$. When
planetary scientists and astronomers use the term Stokes number, they most
often to refer to the stopping time in units of $\Omega^{-1}$, which is
equivalent to the Stokes number at the driving scale of the turbulence.

${\rm St}_{\ell}\equiv\frac{t_{\rm stop}}{t_{{\rm eddy},\ell}}$ (24)

where $t_{{\rm eddy},\ell}$ is the eddy turnover time at scale $\ell$. For
scales smaller than $H$, $t_{{\rm eddy},\ell}<\Omega^{-1}$.

In the case of a turbulent cascade of energy toward smaller scales (e.g.,
Kolmogorov turbulence; _Kolmogorov_ 1941), $t_{\rm eddy}$ decreases with
decreasing spatial scale $\ell$. Thus, for a given pebble size, ${\rm
St}_{\ell}$ increases with decreasing $\ell$.

The necessary conditions for collapse, as enumerated above, are best reached
at scales where ${\rm St}_{\ell}\sim 0.3$ (_Hartlep et al._ , 2017). At larger
scales, the Stokes number is sufficiently small that the pebbles are strongly
coupled to turbulent eddies and thus cannot be concentrated to meet the above
collapse criteria. On smaller scales, the pebbles are only very weakly coupled
to the gas and thus feel the turbulent motions as kicks that do not
appreciably change the pebble density. The “sweet spot" for planetesimal
formation happens at intermediate values, namely ${\rm St}_{\ell}\sim 0.3$
(see _Hartlep et al._ 2017 and _Hartlep and Cuzzi_ 2020).

We will return to this mechanism later in this chapter within the context of
predictions for planetesimal sizes. For now, however, it is worth mentioning
that while this mechanism is promising, more work is needed to study its role
in planetesimal formation and in particular in comparing predictions made by
this model with Solar System constraints.

Finally, it is essential to point out that these processes are not mutually
exclusive with the SI. Indeed, it has recently been shown that pebble
concentration in axisymmetric pressure bumps (_Taki et al._ , 2016; _Onishi
and Sekiya_ , 2017; _Carrera et al._ , 2021, 2022; _Xu and Bai_ , 2022a, b)
can provide a sufficient increase in localized dust-to-gas ratio such that the
SI becomes active near the peak of the pressure bump. However, while this
process may be robust for some parameters (e.g., the strength of the force
generating the bumps; _Carrera et al._ 2022), the simulations by _Carrera and
Simon_ (2022) suggest that for sufficiently small pebbles, planetesimals are
only formed when these pebbles are completely trapped in the pressure bump. In
this case, which requires a bump of sufficiently high amplitude that it
becomes unstable to the RWI (based on the conditions in _Ono et al._ 2016),
planetesimals are only formed by concentrating the pebbles beyond the Roche
density; planetesimals were formed from direct gravitational collapse with no
help from the SI.

Similarly, vortices can also provide sufficient concentration to kickstart the
SI. In particular, _Raettig et al._ (2015) and _Raettig et al._ (2021) carried
out local simulations of pebble concentration in vortices and found that the
SI can be activated for dust-to-gas ratios much less than 1%.

Ultimately, however, these additional mechanisms coupled with the SI deserve
further study. Encouragingly, such efforts are already underway.

###  3.2 Coagulation Models

Although the pebble-collapse model has had a great impact in the scientific
community and can explain many observational facts (see below), there are
alternative planetesimal-formation models that still merit study. In this
section, we discuss the “coagulation models", of which there are two different
types. The first involves the growth (via sticking) of particles beyond the
various barriers discussed above, whereas the second involves the
gravitational growth of small planetesimals into larger ones. We discuss both
here for completeness, though given the limited literature in this area
compared to collapse models, this section is necessarily short.

####  3.2.1 Growth Past the Barriers

The obstacles of the bouncing and fragmentation barriers discussed in Section
2.3 are not absolute. The physical principles to overcome these barriers have
been discussed before (see, e.g., _Blum_ , 2018) and will only be briefly
reviewed here.

The first possibility is to avoid reaching these barriers altogether. _Okuzumi
et al._ (2012) and _Kataoka et al._ (2013) pointed out that extremely small
(i.e., nanometer-sized) monomer grains possess a very high sticking-bouncing
threshold (see Fig. 13 in _Blum and Wurm_ , 2008) and a rather high tensile
(cohesive) strength (see Fig. 4 in _Gundlach et al._ , 2018). This means that
the growing aggregates may remain in the fractal growth regime (see Sect. 2.5)
long enough to grow to macroscopic sizes. _Okuzumi et al._ (2012) and _Kataoka
et al._ (2013) show that aggregates can grow to cm-sizes with a low fractal
dimension $D_{\mathrm{f}}\lesssim 2$ so that their gas-grain response time
does not change by much and radial drift is slow. Although the fractal
aggregates are very rigid, their mechanical stiffness is finite. Thus,
collisions among cm-sized fractal aggregates lead to a compaction without
changing the sticking efficiency. _Kataoka et al._ (2013) also took into
account compaction through gas drag and self-gravity. Beyond the collisional-
compaction size, the aggregates continue to grow with only slightly increasing
densities (from $\sim 10^{-5}-10^{-4}~{}\mathrm{g~{}cm^{-3}}$ for $\sim$cm-
sized aggregates to $\sim 10^{-3}~{}\mathrm{g~{}cm^{-3}}$ for
$\sim$100-m-sized bodies). Self-gravity is responsible for a rather steep
increase in density for bodies above $\sim$100 m in size. Ultimately, the
resulting planetesimals have sizes of $\sim$10 km and densities in the range
$\sim 0.1-1~{}\mathrm{g~{}cm^{-3}}$. With this growth pattern, the radial-
drift barrier (see Sect. 2.2) can be avoided. Recently, _Kobayashi and Tanaka_
(2021) showed that based on such a scenario even the cores of the giant
planets can form in a relatively short time. Although the basic assumptions of
high sticking threshold velocity and high tensile strength is certainly
undoubted for very small monomer grains, the collision behavior of the
extremely fluffy bodies above $\sim$1 cm in size cannot easily be tested in
the lab and thus remains uncertain. In particular, ignoring fragmentation as a
possible collision outcome, deserves a re-visit in future lab experiments.

A second possibility is that of mass transfer during collisions. Although
fragmentation in collisions among dust aggregates in principle leads to the
mass loss of the colliding bodies, there are situations for which this is not
the case for one of the collision partners. If the mass ratio of the two
aggregates is large and the impact energy is moderate, only the smaller
aggregate will fragment and may transfer some of its mass to the intact target
aggregate (_Meisner et al._ , 2013). Based upon the mass-transfer process,
_Windmark et al._ (2012a, b, c); _Garaud et al._ (2013) showed that
planetesimal formation is feasible, in spite of the dominance of fragmentation
in collisions among similar-sized bodies. _Windmark et al._ (2012a, b, c)
based their dust-evolution model on the collision model by _Güttler et al._
(2010) and took into account a variety of collisional outcomes. In order to
overcome the bouncing barrier and enter the mass-transfer regime, _Windmark et
al._ (2012a) artificially inserted 1-cm-sized seed particles, large enough to
start the mass-transfer process. In contrast, _Windmark et al._ (2012b, c)
used velocity-probability distribution functions to initiate the mass-transfer
process. While _Windmark et al._ (2012a) needed $10^{6}$ years to grow
planetesimals with 100 m size, _Windmark et al._ (2012b, c) achieved this
result in $\sim 5\cdot 10^{4}$ years. _Garaud et al._ (2013) used a similar
approach as _Windmark et al._ (2012b, c) and achieved aggregate sizes on the
order of several hundred meters within $10^{4}$ years at 1 au, whereas further
out, the aggregates remained much smaller. Of critical importance for the
credibility of the models is under which conditions the mass-transfer growth
process dominates over mass-loss effects, such as fragmentation, cratering, or
erosion. As reviewed by _Blum_ (2018) and based upon the dust-aggregate
collision model by _Güttler et al._ (2010) with later supplements from lab
experiments, fragmentation is important as long as the size ratio between the
colliding bodies is $\lesssim 10$ under solar-nebula conditions; cratering
becomes the dominating process when the projectile sizes are $\gtrsim 1$ cm,
because in that case the impact energies are too large for mass transfer to
occur; finally, erosion happens when the projectiles are $\lesssim
100~{}\mathrm{\mu m}$ in size (_Schräpler et al._ , 2018). This leaves only a
very small size window for the projectiles to reach the mass-transfer realm.
Whether fragmentation, cratering, and erosion can completely prevent the
formation of planetesimals by the mass-transfer process needs still to be
clarified.

In conclusion, more empirical data on the collision behavior of fluffy
aggregates are required to determine whether protoplanetary growth of solid
bodies by coagulation inevitably stop at pebble sizes or continue towards
planetesimal scales.

####  3.2.2 Formation of 100 km Planetesimals

The second type of coagulation model has been examined in a number of works
(e.g., _Stern and Colwell_ , 1997; _Weidenschilling_ , 1997; _Kenyon and
Bromley_ , 2004, 2012; _Schlichting and Sari_ , 2011; _Schlichting et al._ ,
2013) and is built on the assumption that smaller grains have already grown to
some distribution of initial sizes, often $\sim$ km. These smaller
planetesimals then collisionally grow to larger 10-100 km size bodies.

The basic physics of how this works is relatively straightforward. At the
scales associated with the initial planetesimal populations (e.g., km scales),
$\mathrm{St}\ggg 1$, and drag forces are negligible. Instead, small bodies
grow larger when they collide with one another, “sticking" together primarily
due to gravitational forces. Gravitational focusing also plays a role, though
is significantly more relevant for larger bodies as the square of the impact
parameter is enhanced by a number that is linearly proportional to the mass
(assuming equal mass bodies; see _Armitage_ 2007).

Again, this mechanism for planetesimal formation141414Though, one could argue
that because these models often start with planetesimals, albeit small ones,
that it is more a mechanism to change the size distribution of planetesimals
than it is a process to form them in the first place. has not been as
extensively studied as the gravitational collapse models. However, later in
this chapter, we will come back to it and the predictions that it makes
compared with the SI model.

##  4 PLANETESIMALS IN THE OUTER SOLAR SYSTEM

Having discussed the various mechanisms thought to produce planetesimals, we
now turn to a specific discussion of comets and Kuiper Belt Objects (KBOs). We
focus in particular on comparing observed and/or measured characteristics of
KBOs and comets to those predicted by models.

###  4.1 Mass and Size Distributions

The size distribution of planetesimals in the Solar System has been one of the
primary diagnostics employed to test planetesimal formation theory. In this
section, we first describe the size distribution inferred from planetesimal
formation calculations, followed by a discussion of the impact of collisional
evolution, ending with a comparison with observational data (a brief
discussion on this can also be found in the chapter by Fraser et al.).

####  4.1.1 Primordial Mass Distribution

As discussed previously, a number of works have explored the characteristics
of SI-generated planetesimals, generally agreeing that the largest
planetesimals have masses roughly consistent with large Solar System
planetesimals. In this section, we dive deeper into this topic: the initial
mass and (equivalently) size distributions of planetesimals. As described
above, many SI simulations prevent collapse below certain scales due to
numerical limitations. Thus, in what follows, the size we reference for
gravitationally bound pebble clouds is that which the total mass bound within
a cloud would have if it were collapsed to material densities, generally
assumed to be $\sim 1{\rm~{}g~{}cm^{-3}}$.

_Johansen et al._ (2015) carried out high resolution simulations of SI-induced
planetesimal formation in order to allow a large range of planetesimal sizes
to form. Modeling the mass distribution of planetesimals as a simple power
law,

$N(>M_{p})\propto M_{p}^{-p+1},$ (25)

they found that $p\approx 1.6$. More recently, a number of works followed up
on this and found that $p\approx 1.6$ independent of relative strength of
gravity (effectively the value of $Q$), numerical resolution (_Simon et al._ ,
2016), $\mathrm{St}$, $Z$ (_Simon et al._ , 2017), and $\Pi$ (_Abod et al._ ,
2019); that is, with this simple power law shape to the mass and size
distributions, the so-called “slope" appeared to be nearly universal.

While _Abod et al._ (2019) did fit a simple power law to their planetesimal
masses, they were able to increase the number of detected bound clouds with a
tree-based algorithm based on halo finders in cosmological simulations (the
Plan code; _Li_ 2019). They found that an exponentially tapered distribution,

$N(>M_{p})\propto M_{p}^{-p+1}\exp\left[-M_{p}/M_{c}\right],$ (26)

where $M_{c}$ is a characteristic “cut-off” mass, fit the planetesimal masses
significantly better than the simple power law. In their fit, $p\approx 1.3$,
approximately independent of $\Pi$ (other parameters were not explored in this
work).

This work was roughly consistent with _Johansen et al._ (2015) and _Schäfer et
al._ (2017), though both studies fit a slightly modified version of the
exponentially tapered distribution, introducing a third parameter, $\beta$:

$N(>M_{p})\propto
M_{p}^{-p+1}\exp\left[\left(-M_{p}/M_{c}\right)^{\beta}\right].$ (27)

Given their resolution limits, _Schäfer et al._ (2017) were unable to
constrain $p$ but found that $\beta=0.3$–0.4, independent of their domain
size. _Johansen et al._ (2015) carried out simulations in a smaller domain and
at higher resolution (in terms of grid cells per $H_{\mathrm{g}}$). They found
best fit values $p=1.6$ and $\beta=4/3$; the slope of the power law component
was consistent with _Simon et al._ (2016) and _Simon et al._ (2017) but the
extent of exponential tapering was significantly larger than _Schäfer et al._
(2017).151515We reserve discussion of the remaining parameter, $M_{c}$, to
Section 4.1.3 below.

Most recently, _Li et al._ (2019) used Plan and the highest resolution SI
simulations to date to further quantify the mass distribution. Analyzing two
simulations with different $\mathrm{St}$ and $Z$ values, they found that the
best fit model differed for each simulation. Furthemore, when comparing
identical models (e.g., the exponentially tapered power law), the best fit
parameters were different at a statistically significant level.

Taken together, these results call into question the originally hypothesized
universality of the initial mass distribution. However, further work is
required to determine what trends emerge with varying physical parameters.
Furthermore, a more rigorous analysis is called for in which the same tools
(e.g., advanced clump-finding with Plan and improved statistics via Markov-
Chain Monte Carlo techniques) are applied consistently across a large region
of parameter space.

Despite this complexity, there are some consistent features that emerge from
these calculations. The slope of the mass distribution towards smaller masses
generally falls within the range 1–2 (though not in all models; see _Li et
al._ 2019) and thus the mass function is dominated in number by small
planetesimals but in mass by large planetesimals (see Figure 4 of _Li et al._
2019). Furthermore, the mass distributions are in the broadest sense composed
of a low mass region (or regions) with shallow slopes and higher mass regions
that drop more rapidly with mass, either in the form of an expoential cut-off
or a steeper power law.

Of course, the ultimate question is what do these results tell us, if
anything, when compared to observations? This is a point to which we will
return shortly. However, as these results only address the initial
planetesimal mass distribution, we must briefly discuss the role of post-
formation collisional effects and how these may affect our understanding of
the mass and size distributions of planetesimal populations.

####  4.1.2 Collisional Evolution

There are two main stages to the collisional evolution of comets and KBOs.
During the first stage, bodies presumably resided in a massive ($\sim 20$
$M_{\oplus}$) trans-Neptunian disk (see Chapter 4 by Kaib & Volk, Chapter 5 by
Fraser et al. and _Morbidelli and Nesvorný_ 2020 for further discussion of KBO
dynamics) and experienced strong collisional grinding over the lifetime of the
disk of solids ($t_{\rm s,disk}$). They were subsequently implanted into the
Kuiper belt, scattered disk, and Oort cloud when Neptune migrated to 30 au
(e.g., _Dones et al._ 2015). The collisional evolution during the subsequent
4-plus Gyr, the second stage of collisional evolution, is thought to be less
significant; at these times, these populations were distant and/or did not
contain much mass (_Singer et al._ , 2019; _Greenstreet et al._ , 2019;
_Morbidelli et al._ , 2021).

The collisional evolution during the massive disk stage was studied in
_Nesvorný et al._ (2018). They showed that the Patroclus-Menoetius (P-M)
binary in the Jupiter Trojan population poses an important constraint on the
disk lifetime. This is because the longer the P-M binary stays in the disk,
the greater is the likelihood that its components will be stripped from each
other (by impacts). This indicates $t_{\rm s,disk}<100$ Myr (_Clement et al._
, 2018; _de Sousa et al._ , 2020; _Morgan et al._ , 2021).

Different studies assumed different initial size distributions and impact
scaling laws to model the collisional evolution of KBOs/comets (_Pan and Sari_
, 2005; _Fraser_ , 2009; _Campo Bagatin and Benavidez_ , 2012). For example,
_Nesvorný and Vokrouhlický_ (2019) adopted the initial size distribution
inspired by the SI simulations (e.g., _Simon et al._ 2017), and impact scaling
laws for weak ice, to demonstrate that the collisional grinding would remove
very large bodies and produce the KBO size-distribution break at $D\sim 100$
km (_Bernstein et al._ , 2004; _Fraser_ , 2009). If that is the case, most
comets with $D<10$ km must be fragments of large planetesimals since in
collisional grinding, the collisional lifetime increases with the size of the
body (_Morbidelli and Rickman_ , 2015). The small survivors would have
suffered multiple reshaping collisions (_Jutzi et al._ , 2017). Bi-lobed
comets should either be the product of sub-catastrophic collisions (_Jutzi and
Benz_ , 2017) or fragments of catastrophic collisions (_Schwartz et al._ ,
2018).

Jupiter Trojans have a well defined cumulative power index $-2.1$ from below
10 to 100 km (_Emery et al._ , 2015), which is very close to the equilibrium
slope expected for the gravity regime ($D>0.5$ km; _O’Brien and Greenberg_
2003). This represents a very important constraint on the collisional
evolution. Specifically, if the Jupiter Trojan size distribution was
established during the massive disk stage, the disk must have been dynamically
cold (large impact velocities would produced waves in the size distribution
that are not observed; e.g., _Fraser_ 2009; _Kenyon and Bromley_ 2020). It is
also possible, however, that $t_{\rm s,disk}=0$, which corresponds to
immediately after gas dispersal, in which case the observed size distribution
of Jupiter Trojans and KBOs above $\sim$10 km would be primordial (and an
important constraint on the formation processes). Even in this case, however,
it is difficult to avoid the collisional disruption of comet-size bodies
($D<10$ km; e.g., _Jutzi et al._ 2017; _Nesvorný and Vokrouhlický_ 2019),
unless comets are jocularly stronger that we think.

The first stage of evolution in the massive planetesimal disk, when most
collisional grinding happened, only applies to bodies that formed in the
massive disk below 30 au. This presumably includes most bodies in the
dynamically hot populations of the Kuiper belt (hot classicals, resonant
objects, etc.), and comets, making a direct comparison with SI planetesimal
formation models difficult (as the SI simulations do not include collisional
evolution). However, cold-classical KBOs are thought to have formed in a low-
mass planetesimal disk at $>40$ au (_Parker and Kavelaars_ , 2010; _Batygin et
al._ , 2011; _Dawson and Murray-Clay_ , 2012), not significantly collisionally
evolved during stage two (stage one does not apply to them), and thus likely
represent the most pristine planetesimals in the Solar System. As such, we
restrict our comparison to these objects in what follows.

####  4.1.3 Comparison Between Models and Cold Classical KBOs

To best compare numerical models with observations of the cold-classicals (a
discussion of the other KBO populations is covered in Chapter 5 by Fraser et
al.), we can rewrite the relevant parameters of our mass functions as
equivalent size distribution parameters. Namely, the cumulative size
distribution for bodies of diameter $D$ in the three-parameter exponentially
tapered power law (which we choose based on its common use in the literature
thus far) is

$N(>D)\propto D^{-q+1}\exp\left[\left(-D/D_{c}\right)^{\gamma}\right]$ (28)

where $q=3p-2$, $D=2\left(M_{c}/4\pi\rho_{\mathrm{s}}\right)^{1/3}$, and
$\gamma=3\beta$.

From the discussion in Section 4.1.1 and using the above relation between $q$
and $p$, $q$ ranges from 1.9 in the tapered models of _Abod et al._ (2019) to
2.8 for the simple power law (_Simon et al._ , 2016, 2017) and the tapered
model of _Johansen et al._ (2015).

On the observational side, _Fraser et al._ (2014) analyzed the H band
magnitude distribution of KBOs by combining the results of a number of
independent KBO surveys. In fitting their data to two separate power laws
(joined at a “break point") of the form

$N(<H)\propto 10^{\alpha H}$ (29)

they were able to determine $\alpha\approx 1.5$ for the largest planetesimals,
translating to $q=5\alpha+1=8.5$, and $\alpha\approx 0.38$ for small
planetesimals, equating to $q=2.9$. The magnitude at the break corresponds to
a KBO diameter of 140 km.

More recently, _Kavelaars et al._ (2021) re-examined this problem with new
data collected from OSSOS as well as previous observations. They found that
the cumulative magnitude distribution is best fit with an exponential cut-off
corresponding to diameters ranging from 80 km to 130 km (assuming an albedo of
0.15) and a slope for small objects (the data for which were collected from
other sources) very similar to the _Fraser et al._ (2014) value.

While the $q\approx 2.8$ from the numerical models is temptingly close to the
value inferred from these observations, we must reiterate that $q=2.8$ is only
the best fit for simpler power laws (which are clearly not an accurate
representation of planetesimal size distributions). In fact, with the
exception of the one simulation in _Johansen et al._ (2015),161616While this
simulation has a slope that agrees with observations at small sizes, there is
significant mismatch between the simulations and observations at larger sizes
(not shown in Fig. 9). the models that are fit with an exponentially tapered
power law (which at least agrees with the shape of the _Kavelaars et al._
(2021) fit in principle, though see below) have much smaller $q$ values.

Beyond the slope, one can also compare characteristic planetesimal sizes. A
fit to a characteristic diameter was not analyzed explicitly in _Simon et al._
(2016) and _Simon et al._ (2017), but the effective diameters of their largest
planetesimals were $\sim 200$–800km. Similarly, _Schäfer et al._ (2017) found
large planetesimals with diameters $\approx$ 600–1200 km (when scaling their
units to the Kuiper Belt), and the maximum planetesimal diameters of _Li et
al._ (2019) ranged from 320 km to 1200 km. On the other hand, _Abod et al._
(2019) found a narrower range of characteristic diameters: $\sim 200$–400 km,
which is more consistent with observations, though still slightly larger than
those inferred from observations.

Finally, one should consider the total number of planetesimals, as the
debiased data in _Kavelaars et al._ (2021) put a constraint on the total
number of cold-classical KBOs. When adjusting for the fact that local,
shearing boxes represent only a very small fraction of the Cold-Classical
belt, SI calculations generally produce too many planetesimals. For instance,
the $\mathrm{St}=0.3,Z=0.02$ model of (_Simon et al._ , 2017) produces a
factor of $\sim$ 10 more planetesimals (at least at the resolution probed by
the simulations) compared to the OSSOS data.

To summarize the take away points of this comparison, Figure 9 shows the
cumulative size distribution from three SI simulations representing different
regions of numerical and physical parameter space as well as a fit to the
observed data as described in _Kavelaars et al._ (2021). Clearly, there
remains a mismatch between the size distribution shape, the total number of
planetesimals, and the maximum (or characteristic) planetesimal size. Though,
as shown in the right plot of Fig. 9, the slope of the SI simulations with
$\mathrm{St}=2,Z=0.1$ do match the observations reasonably well when the
theoretical curves have been manually rescaled to lie on top of the
observational data.

Figure 9: A fit to the observed cumulative size distribution of KBOs from the
Outer Solar System Origins Survey (solid blue; _Kavelaars et al._ 2021)
compared with theoretical predictions made by SI simulations with different
parameters (orange pluses, green stars, red triangles) and from the
gravitational coagulation model of _Schlichting et al._ (2013) (purple
circles). In the left plot, both the SI and coagulation calculations have been
adjusted to account for the full extent of the Cold-Classical belt (e.g., we
multiplied the data from _Schlichting et al._ (2013) by the projected area of
the belt on the sky, which we assumed to be $180^{\circ}\times 5^{\circ}$).
The right plot shows the same data but with the models arbitrarily re-scaled
so that at least some of the predicted distribution lines up with the data.
For the SI simulations, $\mathrm{St}$, $Z$, and the resolution are varied.
Neither the SI simulations, nor the coagulation model, match the observations,
both in terms of numbers of planetesimals and the shapes of the distributions.

These discrepancies between theory and observation are not yet resolved. One
possibility regarding planetesimal sizes and the shape of their size
distribution is that due to pebble clouds not collapsing below a numerically
limited spatial scale in the Athena simulations (see above), such clouds may
accrete pebbles at an unrealistically large rate or suffer interactions with
each other (e.g., tidal stripping, merging). However, this issue is not
present with the sink particle runs of Pencil, and there are still
discrepancies between those simulations and KBO observations (e.g., _Schäfer
et al._ 2017).

Another possibility is that since the continued collapse of the pebble clouds
leads to the formation of multiple smaller planetesimals, as discussed in
Section 3.1.4 and _Nesvorný et al._ (2021), the initial mass/size
distributions will change. While the fact that many of these clouds put most
of their mass into a binary seems to suggest that the distributions will not
change appreciably, more work is required to resolve this issue. In
particular, a larger exploration of parameter space (both for the SI and for
the relevant parameters in the Pkdgrav calculations) is required, as well as a
direct calculation of the size distribution from the final collapsed
planetesimals.

It is also possible that a larger volume of the SI parameter space requires
exploration. Even with the many publications on the SI mass function, there
remain many parameter combinations that have not been explored. For the sake
of computational expediency, most of the SI simulations have been carried out
far from the border between strong clumping and no strong clumping (see Fig.
6). In reality though, as pebbles grow from smaller grains, regions of the
disk may move from left to right in the $Z$-$\mathrm{St}$ parameter space, and
thus the relevant $\mathrm{St}$ values may be smaller than what has been
studied so far. While, to the best of our knowledge, small $\mathrm{St}$
values have only been explored within the context of the size distribution in
_Simon et al._ (2017), a less sophisticated clump-detection method in that
work led to poorer sampling of planetesimals, ultimately resulting in a simple
power law fit to the data (as discussed above). Similarly, the $Z$ values may
be closer to the border and not as high as $Z=0.1$, unless there is a sudden
rapid increase in the dust-to-gas ratio in the disk such that the system is
rapidly placed well within the strong clumping regime.

Indeed, simulations along this border have shown some indication of marginal
behavior (e.g., transient excursions into the strong clumping regime; _Li and
Youdin_ 2021). If true, then this would be a natural route towards reducing
the number of planetesimals formed by the SI by e.g., limiting the amount of
time or spatial regions in which strong clumping occurs. How exactly this
marginal behavior would affect the size distribution and maximum planetesimal
sizes is less clear, but certainly remains a question worthy of study.

Finally, it is also possible that including a range of particle sizes will
change the results, given the role that multiple particle sizes play in the
linear regime of the SI (_Krapp et al._ , 2019; _Paardekooper et al._ , 2020;
_Zhu and Yang_ , 2021). However, as discussed above, simulations of the non-
linear state of the SI that include the vertical component of stellar gravity
suggest that the nature of strong clumping does not change drastically change
when multiple particle sizes, at least not to the extent suggested by linear
analyses (_Bai and Stone_ , 2010a; _Schaffer et al._ , 2018, 2021).

The coagulation models for planetesimal formation also make specific
predictions for the size distribution of these bodies. _Schlichting and Sari_