# Architectures of Compact Multi-planet Systems:
Diversity and Uniformity

Lauren M. Weiss Department of Physics, University of Notre Dame, Notre Dame,
IN 46556, USA Sarah C. Millholland Department of Astrophysical Sciences,
Princeton University, Princeton, NJ 08544, USA Erik A. Petigura Department of
Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
Fred C. Adams Department of Physics, University of Michigan, Ann Arbor, MI
48109, USA Konstantin Batygin Division of Geological and Planetary Sciences,
California Institute of Technology, Pasadena, CA 91125, USA Anthony M. Bloch
Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA
Christoph Mordasini Physikalisches Institut, Universität Bern,
Gesellschaftstrasse 6, 3012 Bern, Switzerland

###### Abstract

“I do not believe that even in a snowflake this ordered pattern exists at
random.” – Johannes Kepler

One of the most important developments in exoplanet science in the past decade
is the discovery of multi-planet systems with sub-Neptune-sized planets
interior to 1 AU. This chapter explores the architectures of these planetary
systems, which often display a remarkable degree of uniformity: the planets
have nearly equal sizes, regular orbital spacing, low eccentricities, and
small mutual inclinations. This uniformity stands in sharp contrast to the
diverse nature of the exoplanet sample considered as a whole (as well as our
inner solar system). We begin with a critical review of the observations —
including possible biases — and find that these peas-in-a-pod planetary
systems are apparently a common outcome of the planet formation process.
Modest departures from exact uniformity suggest additional patterns, such as
the planet mass slowly increasing with semi-major axis. The star formation
process naturally produces circumstellar disks with the properties required to
produce these planetary systems, although the solid material must move inward
from its initial location. We discuss primary modes of planetary assembly, the
role of orbital migration, and post-nebular atmospheric loss. Mature planetary
systems are found to be near their minimum energy (tidal equilibrium)
configurations; this finding provides a partial explanation for their observed
properties and indicates that efficient energy dissipation must occur.
Finally, we consider population synthesis models and show that peas-in-a-pod
patterns emerge with reasonable choices for the input parameters. Nonetheless,
interesting observational and theoretical challenges remain in order to
understand how these surprisingly organized planetary systems arise from the
disorder of their formation processes.

##  1 Introduction

With thousands of extrasolar planets now detected, the inventory of the
possible architectures for planetary systems is now coming into focus. One
defining feature of this ensemble of data is that the planet formation process
can produce an enormous diversity of outcomes. The possibilities include large
planets found in tight orbits ($P\lesssim 10$ days; hot Jupiters), smaller
planets with ultra-short periods ($P\lesssim 1$ day), planets in binary star
systems, eccentric Jovian planets, planets much larger than Jupiter on distant
orbits ($a\sim 100$ AU), Earth-sized planets in potentially habitable orbits,
and even some systems roughly analogous to our own solar system.
Significantly, the collection of detected planets also includes a large number
of super-Earth and sub-Neptune planets — bodies with sizes ranging from 1 to 4
$R_{\oplus}$ and masses ranging from roughly 1 to 15 $M_{\oplus}$ — sizes and
masses not represented in our solar system.

In contrast to the diversity found in the exoplanet sample taken as a whole, a
large fraction of sub-Neptune sized planets reside in multi-planet systems
that display unexpected uniformity. These latter systems have planets with
nearly equal sizes (and, where measurements are available, masses), with
regularly-spaced orbits that are nearly circular and coplanar. Moreover,
planetary pairs in these systems are generally not found in mean-motion
resonance. These well-ordered compact planetary systems, and their degree of
intra-system uniformity, are the primary focus of this review.

To date, the majority of known multi-planet systems were discovered by the
NASA Kepler space telescope during its the prime mission (2009–2013). Kepler
discovered thousands of planets with sizes between that of Earth and Neptune,
with $\sim 40\%$ of them in systems containing more than one planet (_Batalha
et al._ 2013; for a comprehensive review of results from Kepler, see the
chapter by Lissauer, Batalha, & Borucki). The Kepler population is
particularly valuable due to the observing strategy of the mission;
$\sim$150,000 stars were monitored nearly continuously for 4 years, regardless
of whether or not they hosted transiting planets. Because the observing
strategy was straightforward and well-documented, it is possible to model how
it has shaped planet detections, and one can make projections to characterize
the underlying planet population. For example, _Mulders et al._ (2018)
estimate that $\gtrsim$42% of Sun-like stars have nearly coplanar planetary
systems with seven or more exoplanets. Although the occurrence rate is subject
to some uncertainty, such multi-planet systems represent one of the most
common outcomes of planet formation. For completeness, we note that subsequent
transit missions (the Kepler extended mission and the NASA TESS mission) and
radial velocity surveys have found additional multi-planet systems, with the
latter important for making mass estimates and detecting long-period giant
planets.

Before characterizing these compact multi-planet systems in detail, we must
specify the parameter space spanned by the systems of interest. For the sake
of definiteness, we adopt the following working definition: We define a
compact multi-planet system (‘compact multi’) to be a planetary system
containing multiple planets with radii between 0.5 and 4.0 $R_{\oplus}$ and
with orbital periods between one day and one year. In Figure 1, we show these
planets in the period-radius ($P$-$R_{p}$) plane, alongside the broader
population of known planets. Our definition is agnostic about the presence (or
absence) of giant planets and/or additional stars111The role of giant planets
in shaping the architectures of compact multis is certainly interesting, but
is beyond the scope of this review.. In spite of its utility, this
observationally-motivated definition has the following shortcomings: (1) many
single transiting systems may indeed be compact multis viewed from glancing
angles, and (2) some systems may be closely related to the compact multis from
a formation standpoint, but may fail our size and period criteria.

With the above definition in place, we note that most — but not all — of the
compact multi-planet systems contain planets with nearly equal sizes (masses)
and regularly spaced orbits. This so-called “peas-in-a-pod” pattern of intra-
system uniformity is an emergent property of compact multi-planet systems.
Note that “peas-in-a-pod” is a pattern, not a particular group of planetary
systems, although it is possible to define subsets of the compact multis that
adhere to this pattern based on a variety of metrics. The specification of
these well-ordered systems, and the mechanisms that lead to their formation,
represent a defining theme of this chapter.

For purposes of illustration, it is useful to highlight particularly
significant compact multi-planet systems according to the definitions given
above. Kepler-11 is a noteworthy and historically significant system
(_Lissauer et al._ , 2011a). This example, along with other early Kepler
discoveries, showed that the inner regions of many planetary systems differed
from that of the solar system and from the predictions of the leading
theoretical models of the day. Kepler-11 hosts five transiting planets
interior to Mercury’s orbit with masses of $M_{p}\approx$ 2–8 $M_{\oplus}$.
The total mass in these planets is $M_{T}\approx 22~{}M_{\oplus}$ compared to
only $0.05~{}M_{\oplus}$ in Mercury itself, an enhancement relative to the
solar system by a factor of $\sim$400 (_Lissauer et al._ , 2013). The planets
themselves range from 1.8 $R_{\oplus}$ to 4.2 $R_{\oplus}$ and provide early
examples of close-in Neptunes and sub-Neptunes. The planets have a 0.11 dex
dispersion in size (about 30% fractional dispersion). Many compact multi-
planet systems have an even narrower dispersion in size. Figure 1 shows three
of the most uniform Kepler systems, having size dispersion of 0.04 dex or
less. Systems like these motivate the theoretical work presented here.

On the other hand, not all compact multis have such a narrow size dispersion.
Figure 1 also highlights WASP-47 and Kepler-89 (a.k.a. KOI-94), two well-known
systems that host both sub-Neptune and Jovian-size planets (_Becker et al._ ,
2015; _Weiss et al._ , 2013). Such systems are relatively rare, but
theoretical models must allow for occasional exceptions to uniformity.

As noted earlier, some systems that meet the definition for being a compact
multi could have formed differently from the majority of such systems, whereas
others that fail the definition could have formed through similar mechanisms.
V1298 Tau is one such system (_David et al._ , 2019). With four planets with
sizes $R_{p}=6$–$10~{}R_{\oplus}$, it does not meet our compact multi
definition. However, given that the system has an age of only 20 Myr, it is
plausible that these planets are actively undergoing Kelvin-Helmholtz
contraction and/or atmospheric loss. The V1298 Tau system could thus be on its
way to becoming a compact multi in the future.

Our own solar system is another useful point of reference. Here, the planetary
orbits are mostly circular and coplanar, with orbital eccentricities $e<0.2$
and inclinations $i<7^{\circ}$. The terrestrial planets are close to the Sun
(with $a<1.6$ AU), whereas the gas and ice giants inhabit the space beyond the
snow line ($a>5$ AU). Our solar system is thus relatively well-ordered, and
has been used in the past as the ultimate example of a deterministic system or
‘clockwork’ universe (_Newton_ , 1687).222Of course, the advent of chaos
complicates this description. As outlined below, however, a large fraction of
the compact multi-planet systems of interest here are even more well-ordered
than the solar system.

This chapter is organized as follows. First, we review the prevailing physical
and orbital patterns found in compact multis (§2), including the degree of
uniformity in their masses, radii, and orbital spacing. Next, we review key
aspects of star/disk formation that set the boundary conditions for disk
properties and planet formation (§3). We then explore the specific processes
by which planets form in their disks, with a particular focus on the physics
that sets the characteristic mass scale for compact multi-planet systems (§4).
We also explore the planet-planet interactions that regulate mass growth, and
we develop a pairwise energy optimization model that naturally produces peas-
in-a-pod architectures (§5). With these results in place, we consider
population synthesis models (§6) and show that peas-in-a-pod patterns can be
produced under reasonable assumptions. The chapter concludes (§7) with a
summary of results and a discussion of how future work can move forward our
understanding of these well-ordered planetary systems.

Figure 1: Top: Census of confirmed extrasolar planets discovered by the
transit method. Top: we define a compact multi to be a system with at least
two planets with sizes and periods residing in the box. Middle: a few examples
of the many compact multis that have a narrow dispersion of size and orbital
spacing. Bottom: WASP-47 and KOI-94 are two systems that meet our definition
of compact multi but have a large dispersion in size. V1298 Tau fails our
definition, but may be evolving toward a compact multi in the distant future.
(Source: NASA Exoplanet Archive, 2021-06-03).

##  2 Properties of Compact Multi-Planet Systems

###  2.1 Homogeneous Catalogs

The discovery of several striking attributes of the compact multi-planet
systems relied on large, homogeneous catalogs of stellar and planet
properties. A particularly important advance came from the California Kepler
Survey (CKS, _Petigura et al._ 2017), in which high-resolution spectra from
the W. M. Keck Observatory were used to derive a homogeneous catalog of Kepler
stellar and planet properties for 2025 planets orbiting 1305 stars. This was
the first homogeneous catalog of over a thousand planet-hosting stars. More
recently, thanks to the ESA Gaia mission, other large catalogs of
homogeneously determined Kepler stellar and planet properties have been
constructed (e.g., _Fulton and Petigura_ , 2018; _Berger et al._ , 2020).
Their results are generally consistent with findings based on the CKS
catalog.333Note that these catalogs focus on planets from the Kepler Mission;
the K2 and TESS missions found fewer than 10 high-multiplicity (4+) planetary
systems total, and thus have little to contribute regarding the planet-planet
statistics of compact multis.

Figure 2: The distribution of transiting planet radii in the compact multis
from the California Kepler Survey (CKS). There is a significant gap at $\sim
1.8\,R_{\oplus}$. Planets larger than $4\,R_{\oplus}$, which sometimes occur
in systems that have multiple planets smaller than $4\,R_{\oplus}$, are grayed
out. Empirically motivated estimates of the planet masses, based on _Weiss and
Marcy_ (2014), are provided along the top axis. Assuming planets smaller than
$1.5\,R_{\oplus}$ are rocky, the distribution includes planets as low-mass as
$0.1\,M_{\oplus}$. The majority of planets $\sim 1.5\,R_{\oplus}$ and larger,
which likely have $\sim 5\%$ of their mass in volatile envelopes of hydrogen
and helium, are 3 to $10\,M_{\oplus}$. Figure 3: Compact multi-planet systems
with four or more transiting planets interior to 1.52 AU (data from the
California Kepler Survey, _Weiss et al._ , 2018a). The systems are ranked by
their planet radius dispersions $\sigma_{R}$, with the most uniform sizes at
the top and the greatest size diversity at the bottom. The point sizes
represent the planet radii on a logarithmic scale (see legend). The solar
system terrestrial planets are included for comparison; their $\sigma_{R}$ is
the bottom quintile of size uniformity. Although we do not know how many solar
system analogs Kepler missed (since it was insensitive to Mars-sized planets
at Mars-like orbits), the peas-in-a-pod architecture that is prevalent in the
Kepler compact multis did not emerge as strongly in the final architecture of
the solar system.

Patterns emerging from the observed architectures of compact multis, including
aspects of either diversity or uniformity, provide important clues as to how
these ubiquitous planets formed. Such patterns relate to planetary physical
properties (sizes and masses), orbital properties (spacing, eccentricities,
inclinations), or a mix of the two (e.g. size-spacing relationships).
Moreover, we can also examine whether planetary patterns are correlated with
the properties of their host stars. This section reviews the primary patterns
that have been identified in the observed sample of compact multi-planet
systems. For a preview of the full range of observations that will be
discussed, one may skip ahead to a schematic summary diagram in Figure 9.

###  2.2 Uniform Sizes and Masses

#### Radius and Mass Distributions.

The overall distribution of the Kepler planet radii is broad, with the sizes
of detected planets ranging from the size of Mars to Jupiter. However, the
size distribution is far from uniform: Earth-to-Neptune-size planets outnumber
giants by an order of magnitude. As a point of reference, there are roughly
$\sim$90 planets per 100 stars with radii $R_{p}$ = 1–4 $R_{\oplus}$ and
periods less than 1 year, compared to only $\sim$7 planets per 100 stars with
radii $R_{p}$ = 8–32 $R_{\oplus}$ in the same period range (_Petigura et al._
, 2018).

Thanks to the large number of transiting planets discovered by Kepler and
intensive follow-up efforts, fine details have been resolved within the planet
size distribution. Using the CKS sample, _Fulton et al._ (2017) identified a
paucity of planets in the size range 1.5–2.0 $R_{\oplus}$ compared to slightly
larger and smaller planets. This gap aligned with a previously observed
transition in planetary bulk density over the same size range (_Weiss and
Marcy_ , 2014). Smaller planets (known as “super-Earths”) have densities
consistent with iron/silicate compositions while larger planets (“sub-
Neptunes”) have low densities and require volatile envelopes that add
significant volume to the planets (_Rogers_ , 2015). The planet radius gap is
present in both the overall CKS sample and in the CKS compact multis (_Weiss
et al._ , 2018b).

Although masses have been measured for only $\sim$5% of compact multis, we can
obtain mass estimates using a mass-radius relationship. Here we adopt the
piecewise mass-radius relationship of _Weiss and Marcy_ (2014), but note that
there is an astrophysical dispersion about this relationship of $\sim$4.2
$M_{\oplus}$.444We do not model mass dispersion in this chapter. Using this
empirical mass-radius function, we find that the planets in compact multi-
planet systems have masses ranging from 0.1 to 20 $M_{\oplus}$, with the
majority of the detected planets having $M_{p}$ $\approx$ 3–10 $M_{\oplus}$.
Figure 2 shows the full distributions for both radius and mass estimates for
the sample. Although most of the planets are super-Earths and sub-Neptunes,
the distribution has a tail extending below 1 $M_{\oplus}$ (containing
$\sim$8.6% of the planets).

#### Uniformity in Planet Sizes.

Various studies have found that planets within a given multi-planet system
have similar sizes to their neighbors. The size similarity was first noted in
(_Lissauer et al._ , 2011b), and was more robustly confirmed of 909 planets
homogeneously characterized planets in 355 Kepler multi-planet systems (_Weiss
et al._ , 2018a). Most of these planetary systems fit the definition of
compact multis provided in §1 (833 planets in 322 systems). The peas-in-a-pod
pattern is visible by eye in the high-multiplicity (4+ planet) systems shown
in Figure 3. In the most uniform systems, the planets appear to have a
characteristic size that is well-correlated with their neighbors.

One might expect the size of each detected planet to be consistent with a
random draw from a single underlying distribution, with a combination of
detection and selection biases sculpting the pattern. However, the size
similarity found among the compact multis is inconsistent with random draws.
_Weiss et al._ (2018a) and _Weiss and Petigura_ (2020) performed a simple
controlled experiment in which they held fixed the stellar, orbital, and noise
properties of the observed planets, but drew synthetic planet radii at random
from both the observed and log-normal radius distributions. As shown in Figure
4, the observed planetary systems show more uniformity in radius than the
systems produced through random sampling. Several studies have assumed more
complicated underlying radius distributions (including covariance between
planet radius and other properties) in an attempt to simultaneously reproduce
a variety of observed distributions of the compact multis (_Mulders et al._ ,
2018; _He et al._ , 2019; _Sandford et al._ , 2019; _He et al._ , 2020). These
studies also found that the planet sizes are more correlated within a given
planetary system than would be produced through random draws.

Figure 4: Top left: The relationship between the radii of adjacent transiting
planets, $j$ and $j+1$, from the CKS compact multi-planet systems. Top right:
The same, but for planet radii drawn at random with replacement from a log-
normal distribution that survive a mock Kepler detection pipeline and yield
synthetic compact multis. Bottom left: The relationship between the period
ratios of adjacent pairs of transiting planets from the CKS compact multis.
Note the absence of period ratios $<1.2$ (dotted lines), which is likely due
to instability. Bottom right: The same, but for period ratios drawn at random
from the observed distribution, with replacement (i.e., bootstrap). The random
draws of planet radius and period ratio, when subjected to detection bias,
selection criteria, and the stability requirement $\mathcal{P}>1.2$, do not
resemble the observed planetary systems.

Just how similar are the planet sizes in a given system? A useful way to
approach this question is to define a metric for size similarity. The metric
should be independent of the typical planet size. For a given system of $N$
detected planets indexed by increasing orbital distance ($j=1,2,\dots N$), we
consider the fractional dispersion of the planet radii,

$\sigma^{2}_{\mathcal{R}}={\rm
Variance}\left\\{\mathrm{log_{10}}(R_{j}/R_{\oplus})\right\\}\,,$ (1)

where we take the logarithm of planet radius to ensure that the dispersion is
fractional and we use $R_{\oplus}$ as a reference scale to make the variable
dimensionless. The choice of log-base 10 is arbitrary. Note that the size
dispersion is only defined for systems with three or more planets; a similar
metric for the dispersion in planet spacing (presented below) requires four or
more planets. In this review, we compute the size and spacing variances for
systems with four or more transiting planets. For completeness, we note that
the literature contains many additional metrics for size diversity (e.g.,
_Millholland et al._ , 2017; _Gilbert and Fabrycky_ , 2020).

In the sample of 4+ planet systems plotted in Figure 3, the systems are ranked
in order of increasing $\sigma_{\mathcal{R}}$, with the systems with the
lowest size dispersions at the top. The median value of the fractional
dispersion is $\sigma_{\mathcal{R}}=0.14$ dex. This value is lower than the
median size dispersion in systems with radii drawn at random (0.18 dex). As
another point of comparison, the inner solar system (Mercury, Venus, Earth,
and Mars) has $\sigma_{\mathcal{R}}=0.20$ dex, placing it in the most
disordered quintile of high-multiplicity systems. Several multis with more
diverse sizes than the terrestrial planets have a giant planet within 1 AU
that is partially responsible for the size dispersion (e.g., KOI-834; KOI-94,
_Weiss et al._ 2013; and KOI-351 _Schmitt et al._ 2014; _Cabrera et al._
2014), although some systems have sub-Neptunes with diverse sizes (e.g.,
KOI-70, _Buchhave et al._ 2016; KOI-82, _Marcy et al._ 2014).

Although there is some size dispersion in all systems, the variations in
planet radii tend to be orderly. For example, _Ciardi et al._ (2013) found
that the outer planet is larger in $68\pm 6\%$ of pairs of transiting planets
(where this trend is limited to cases where the smaller planet would be
detectable with the orbital period of either planet). For comparison, this
pattern holds for 66% of pairs in the inner solar system (with Mars breaking
the pattern). _Kipping_ (2018) reproduced this result and also found that size
ordering in high-multiplicity systems tends to be monotonic, indicating that
the planet sizes are in a low-entropy state.

#### Uniformity in Planet Masses.

The intra-system similarity found for planetary radii extends to the planet
masses as well. This finding suggests that the bulk compositions of planets,
including their envelope-to-core fractions, are also likely to be relatively
uniform. _Millholland et al._ (2017) first demonstrated mass uniformity using
a sample of 37 systems with masses determined from Transit-Timing Variations
(TTVs), as measured by _Hadden and Lithwick_ (2017). They found that masses of
planets were more similar to the masses of neighboring planets than would be
expected from random draws. _Millholland et al._ (2017) also identified radius
uniformity, radius ordering, and mass ordering within this sample. _Wang_
(2017) later studied a sample of multi-planet systems with mass measurements
from radial velocities and showed that this sample also exhibited intra-system
mass uniformity. The tendency for uniformity in both radii and masses implies
that the significant scatter in the mass-radius relation (the $M_{p}$-$R_{p}$
plane) is dominated by system-to-system variance rather than intra-system
variance.

For completeness, we note that the current determinations of mass uniformity
rely primarily on mass estimates made from TTVs (with the exception of _Wang_
2017). TTVs are most readily detected for large planets with ultra-compact
(period ratios $\lesssim 2$) or near-resonant orbital spacing (_Mills and
Mazeh_ , 2017), which might necessitate that they have low masses to maintain
stable orbits. Thus, it remains possible that the systems that display TTVs
are systematically different from those systems that do not (_Weiss and Marcy_
, 2014; _Lee and Chiang_ , 2016). _Millholland_ (2019), for instance, showed
that planets wide of resonance have systematically larger radii than non-
resonant planets, suggesting some degree of physical uniqueness of the TTV
planet sample. Exploring whether the similarity in planet masses and
compositions robustly extends beyond the TTV systems will be a major avenue of
investigation in the era of sub-meter-per-second RVs.

Despite this caveat, there is ample reason to suggest that intra-system mass
uniformity is a widespread trend. First, the intra-system uniformity of planet
radii exists both for systems with larger planets (which are likely to have
envelopes) and for the 72 planets with radii smaller than Earth, including a
tail of Mars-sized planets. At an orbital distance of 0.1 AU, a Mars-sized
planet is unlikely to hold a volatile atmosphere. We can thus assume that the
Mars-sized planets are rocky, and the size similarity in systems such as
KOI-3158 (see Figures 1 and 3) is indicative of similarity in mass (in
addition to radius). Moreover, even among planets with envelopes, atmospheric
mass loss (see §4.4) tends to reduce planet radii down to $R_{p}\approx 2\
R_{\mathrm{core}}$ (_Owen and Wu_ , 2017), such that their uniform radii are
indicative of uniform core sizes (and thus masses, assuming compositional
similarity). Lastly, mass uniformity within systems is a fundamental outcome
of the planet formation process (as discussed in §3–6).

###  2.3 Uniform, Non-Resonant Orbital Spacing

Figure 5: Left: Planet occurrence as a function of orbital period for the
following size classes: super-Earths ($R_{p}=1.0{-}1.7\ R_{\oplus}$), sub-
Neptunes ($R_{p}=1.7{-}4.0\ R_{\oplus}$), sub-Saturns ($R_{p}=4{-}8\
R_{\oplus}$), and Jupiters ($R_{p}=8{-}24\ R_{\oplus}$). The Figure is from
_Petigura et al._ (2018) who removed Kepler’s observational biases when
computing the distributions. Right: Period ratios of adjacent pairs of planets
discovered by Kepler. Top panel: The full period ratio distribution of Kepler
DR25 compact multis (containing planets with $0.5\ R_{\oplus}<R_{p}<4\
R_{\oplus}$ and $P<1$ year, as defined in §1). The smooth red line is a fit of
a shifted log-normal distribution. Kepler’s observational biases have not been
removed. Pairs with period ratios within 0.06 of the three strongest first-
order MMRs are shaded. While there are overabundances wide of MMRs, such pairs
make up a small fraction of the total distribution. Bottom panel: Period
ratios of pairs of planets within four example 4+ planet systems,
demonstrating that individual systems tend to have period ratios that cluster
within the period ratio distribution.

Kepler’s four-year observational baseline and high-quality photometry enabled
precise measurements of orbital periods, with typical fractional errors of
$\mathcal{O}(10^{-5})$. Consequently, there is a rich literature not only on
planet periods themselves, but also on the pairwise relationships between
periods (i.e., orbital spacing) and higher-order relationships (i.e., orbital
spacing dispersion).

#### Period and Period Ratio Distributions.

Figure 1 shows that compact multis span a broad range of orbital separation.
However, the density of points in this figure does not visually convey the
intrinsic period distribution, because of the declining sensitivity with
period due to falling transit probability and signal-to-noise. _Petigura et
al._ (2018) accounted for these effects and found that the period distribution
of super-Earths and sub-Neptunes rises over the range $P\approx 1$–10 days
according to $dN/d\log P\propto P^{\alpha}$ with $\alpha\approx 2$. For longer
periods, the occurrence function is roughly log-uniform, i.e. $\alpha\approx
0$ (see Figure 5).

Beyond the overall period distribution, further insight can be gained from the
the distribution of period ratios, which reveals pairwise relationships.
Figure 5 shows the distribution of period ratios of pairs of adjacent planets
in compact multis. Here, we concentrate on the tightly-packed region with
$\mathcal{P}_{j}\equiv P_{j+1}/P_{j}<3.5$, leaving out $\sim 20\%$ of planet
pairs. This region is well-described by a shifted log-normal distribution of
the form

$f(\mathcal{P};\delta,\mu,\sigma)=\frac{1}{(\mathcal{P}-\delta)\sigma\sqrt{2\pi}}\exp\left[\frac{-(\ln(\mathcal{P}-\delta)-\mu)^{2}}{2\sigma^{2}}\right],$
(2)

with best-fit parameters given by $\delta=0.61$, $\mu=0.25$, and
$\sigma=0.39$. The mode of the distribution falls at $\mathcal{P}\sim 1.7$.

Two primary features of the period ratio distribution stand out. First, the
distribution is smooth and broad overall, as evidenced by its close agreement
with a log-normal. This indicates that the compact multis have predominantly
non-resonant orbital architectures. A second feature, however, is the
existence of narrow peaks just outside (within about 5%) of the first-order
mean-motion resonances and deficits just inside of the nominal resonance
locations. These near-resonance features were first noticed through studies of
the Kepler multi-planet population (_Lissauer et al._ , 2011b; _Fabrycky et
al._ , 2014), and they have been shown to be statistically significant (e.g.,
_Steffen and Hwang_ , 2015). Both the paucity of resonances and the near-
resonance features offer important challenges for planet formation theories.
We will thus defer theoretical discussion to §4.3, where we will review
efforts to understand their origins.

#### Uniform Orbital Spacing.

While the period ratio distribution provides a population-level summary of
pairwise orbital spacing, another level of sophistication comes through the
study of spacing variations from system-to-system. _Weiss et al._ (2018a)
found that compact multis tend to have regular orbital spacing.555 _Weiss et
al._ (2018a) considered planets with $P_{j+1}/P_{j}<4$ to mitigate the effects
of detection bias mainly because the probability of a long-period planet
transiting is low. That is, planets in a given system are roughly evenly
spaced in terms of log-orbital period in a manner that is inconsistent with
random draws from the ensemble period ratio distribution (Figure 4). Because
planets in such systems orbit a common host star, regular log-period spacing
also implies regular separations in terms of log-semi-major-axis. It is
important to note that the tendency for uniform orbital spacing is not in
conflict with the smoothness of the period ratio distribution. Rather, the
broader distribution is composed of individual systems with period ratios that
are clustered within sections of the distribution. This effect is illustrated
with four example systems shown in the bottom right panel of Figure 5.

How regular is the typical planet spacing? A slight modification of our radius
dispersion metric, equation (1), yields the definition

$\sigma^{2}_{\mathcal{P}}={\rm
Variance}\left\\{\mathrm{log_{10}}(\mathcal{P}_{j})\right\\}\,,$ (3)

where we take the logarithm of each period ratio before computing the
dispersion to keep the comparison fractional. Note that there are only $N-1$
adjacent period ratios in an $N$-planet system. For completeness, we note that
other spacing metrics exist, such as the gap complexity defined in _Gilbert
and Fabrycky_ (2020). In the CKS compact multis with four or more transiting
planets, the median spacing dispersion is 0.05 dex (0.04 dex in systems with
$\sigma_{R}<0.1$). In synthetic 4-planet systems where $\mathcal{P}$ is drawn
at random from the ensemble of observed period ratios with $\mathcal{P}<4$,
the median value of $\sigma_{\mathcal{P}}$ is 0.10 dex.

The metric $\sigma_{\mathcal{P}}$ is especially efficient at identifying high-
multiplicity systems that mostly follow the peas-in-a-pod pattern, but have
one or more planets missing. For instance, systems KOI-2433 and KOI-1306 in
Figure 3 have some of the most regular planet sizes, but each has a clear gap
that could be filled by a planet. These two systems emerge as outliers in the
$\sigma_{R}$-$\sigma_{\mathcal{P}}$ plane (Figure 6).

The solar system has long been known to have regular planet spacing (_Bode_ ,
1768), with ${a_{j+1}/a_{j}}\approx 1.73$ or $\mathcal{P}\approx 2.27$
(_Blagg_ , 1913). The inner solar system (Mercury-Mars) has
$\sigma_{\mathcal{P}}=0.1$ dex. Many of the high-multiplicity CKS compact
multis have lower size dispersion (80%) and lower spacing dispersion (68%)
than the inner solar system. Thus, the Titius-Bode rule apperas to be a
special case–and yet not the strongest example–of the more general empirical
rule that planets tend to have similar sizes and regular orbital spacing.

Within the solar system, the slight deviations to the regular spacing appear
to have some structure that might be astrophysical in origin (_Blagg_ , 1913;
_Nieto_ , 1970). The majority of the compact multis have too few planets to
investigate whether the deviations from regularity are orderly, but such an
investigation will be valuable when higher-multiplicity exoplanet systems are
discovered.

Figure 6: The typical size vs. spacing dispersions, $\sigma_{R}$ and
$\sigma_{\mathcal{P}}$, of the CKS compact multis with 4 or more planets. The
inner solar system (terrestrial planets only) is not atypical, but is more
disorderly in both its size dispersion and spacing dispersion than the
majority of compact multis. Two systems with low size dispersions but high
spacing dispersions are highlighted; these systems have clear gaps in Figure 3
where an additional planet likely resides.

Figure 7: Left: Distribution of orbital separations between 812 pairs of
adjacent planets, plotted in units of mutual Hill radii. Planet masses were
estimated using the empirical mass-radius relations from _Weiss and Marcy_
(2014). Planet pairs with $\Delta>60$ (21) are not shown. The high-
multiplicity systems (4+ planets) tend to have closer packing than lower-
multiplicity systems. Note that planets with apparent large separations might
have an undetected planet in between. Right: The average planet-to-star mass
ratio vs. the average orbital period ratio in the CKS peas-in-a-pod (systems
with 4 or more transiting planets with $\sigma_{R}<0.1$). Both the system-wide
geometric means (squares) and individual planet pairs (circles) are shown. The
grayscale represents their separations in units of mutual Hill radii.

###  2.4 Ordered Distributions of Spacing and Size

One notable attribute of the high-multiplicity (4+ planet) systems is that the
typical period ratio is related to the typical planet size, especially for the
planets that have the least radius dispersion ($\sigma_{\mathcal{R}}<0.1$). We
convert the planet radii to masses through a nominal mass-radius relationship
(_Weiss and Marcy_ , 2014) so that we can compute the planet-to-star mass
ratios, which are dynamically relevant via the mutual Hill radius,

$R_{H}=\left({M_{p,j}+M_{p,j+1}\over
3M_{\star}}\right)^{1/3}\left(\frac{a_{j}+a_{j+1}}{2}\right).$ (4)

The separation between two planets in units of mutual Hill radii is

$\Delta=(a_{j+1}-a_{j})/{R_{H}}.$ (5)

Histograms of the pairwise separations in mutual Hill radii are shown in the
left panel of Figure 7. The typical separations are $\sim 20$ mutual Hill
radii. Systems with more transiting planets tend to have smaller separations.
However, an interesting architectural feature emerges among the highest
multiplicity (4+ planets) systems when we plot the planet-to-star mass ratios
versus the period ratios, representing the typical mass and spacing in each
system with the geometric mean. The geometric mean Hill separation in each
system, computed via equation 5, is indicated in color.

The right panel of Figure 7 displays several striking features. First, there
are no systems with typical period ratios $<1.2$. The absence of planets
closer than this minimum period ratio spacing is likely because resonance
overlap becomes important at $\mathcal{P}<1.2$ for $\sim 1\,M_{\oplus}$
planets, leading to chaotic orbital perturbations and eventual instability
(_Deck et al._ , 2013). Second, the smallest planets (with $R_{p}\lesssim
1\,R_{\oplus}$ and $M_{p}/M_{\star}\lesssim 5\times 10^{-6}$) have orbital
period ratios confined to $1.2<\mathcal{P}<2$, although the planets would have
been both stable and detectable at larger period ratios (_Weiss et al._ ,
2018a). We will explore possible origins for the tight spacing, and also
relationships between planets spacing and size, in §3-6.

###  2.5 Low Eccentricities

As discussed in the previous two sections, the architectural uniformity of the
compact multis extends beyond physical planet properties (masses and radii) to
orbital properties such as semi-major axis spacing. This finding immediately
raises the question of whether additional orbital properties, such as
eccentricities and inclinations, also exhibit orderly structure. The answer
appears to be affirmative; the orbits of the compact multis are nearly
circular (and nearly coplanar — see the following section). Moreover, the
orbital configurations are consistent with being the dynamically coolest
extreme of a spectrum of compact multi-planet systems. In this section and the
next, we review the observations that have led to this conclusion, focusing
first on eccentricities and then proceeding to inclinations.

Eccentricities are difficult to measure for individual compact multi-planet
systems. The multis have RV signals that are small and complex, and their host
stars are generally faint, thus prohibiting strong RV-derived eccentricity
constraints. The majority of individual measurements have been made using
Transit Timing Variations (TTVs). _Hadden and Lithwick_ (2017) performed a
systematic analysis of 145 Kepler planets with detectable TTVs and found
predominantly small but nonzero eccentricities, with a median value of $e\sim
0.02$.

Beyond TTV analyses in individual systems, further eccentricity constraints
have been made using population-level techniques. One approach is to exploit
the planetary transit duration $T$, which depends on eccentricity along with
period, mean stellar density $\rho_{\star}$, argument of periapse $\omega$,
and impact parameter $b$ according to (e.g. _Winn_ , 2010):

$T=\left(\frac{3}{\pi^{2}}\frac{P(1-b^{2})^{3/2}}{G\rho_{\star}}\right)^{1/3}\frac{\sqrt{1-e^{2}}}{(1+e\sin\omega)}.$
(6)

For individual planets, measurements of the duration yield weak constraints on
eccentricity due to uncertainties in $b$. However, one may constrain the
distribution of eccentricities in a sample of transiting planets with measured
$\rho_{\star}$ (_Ford et al._ , 2008), since the distribution of $b$ is known.

This analysis has been performed for a variety of sub-populations of Kepler
systems. _Van Eylen and Albrecht_ (2015) and _Van Eylen et al._ (2019)
analyzed transit durations of $\sim$100 planets with asteroseismically
constrained $\rho_{\star}$ and modeled the eccentricity distribution with a
Rayleigh form with a scale parameter (mode) of $\sigma_{e}$. Taken together,
these two studies found that mulits have lower typical eccentricities than
single planet systems: $\sigma_{e}\sim 0.05$ versus $\sigma_{e}\sim 0.25$.
Similar trends have been observed in larger samples that used different
techniques to measure $\rho_{\star}$. In a sample of $\sim$700 planets with
spectroscopic determination of $\rho_{\star}$, _Xie et al._ (2016) found
$\sigma_{e}\sim 0.03$ for multis and $\sigma_{e}\sim 0.25$ for singles.
Similarly, in a sample of $\sim$1300 planets with hosts with
spectroscopic/astrometric determination of $\rho_{\star}$, _Mills et al._
(2019) found $\sigma_{e}\sim 0.04$ for multis and $\sigma_{e}\sim 0.17$ for
singles.

Taken together, ensemble studies of Kepler planets have found typical
eccentricities of $e\approx 0.15$–0.25 among singles and $e\approx 0.02$–0.05
among multis. In some sense, the low eccentricities of the compact multis are
not surprising given their tight orbital spacings. However, the eccentricities
are even lower than that required by orbital stability. Among a sample of TTV-
active systems with 3 or more planets, _Yee et al._ (2021) found that their
TTV-derived eccentricities were smaller than the maximum eccentricities
allowed by long-term orbital stability by factors of 2–10, suggesting that
compact multis attained their final orbital configurations in the presence of
efficient dissipation. We will revisit the role of dissipation when we discuss
the origin of intra-system uniformity from the perspective of energy
optimization in §5.2.

###  2.6 Low Mutual Inclinations

In addition to being nearly circular, the orbits of the compact multis are
approximately coplanar, in line with the expectations of dynamical
equipartition. This can be seen merely from the fact that the planets are co-
transiting, which requires approximate coplanarity (unless there is a rare
chance alignment of the nodes). For individual systems, quantitative
measurements of mutual inclinations are even more challenging than
eccentricities, and thus almost all constraints are from population-level
analyses.

In the previous section, we described how the distribution of transit
durations provides population-level constraints on eccentricities. Similarly,
the distribution of transit duration ratios encodes information about mutual
inclinations. The ratio of the transit chord lengths,
$\xi\equiv(T_{\mathrm{in}}/T_{\mathrm{out}})(P_{\mathrm{in}}/P_{\mathrm{out}})^{-1/3}$,
between pairs of transiting planets in the same system is sensitive to the
planets’ mutual inclinations through their relative impact parameters
(_Steffen et al._ , 2010). Early analyses of these ratios showed that compact
multis typically have low mutual inclinations (_Fang and Margot_ , 2012;
_Fabrycky et al._ , 2014). They described inclination dispersion as Rayliegh
distributions with scale parameters of $\sigma_{i}\sim 1^{\circ}-2^{\circ}$.

A natural question is whether the small $\sigma_{i}$ of the co-transiting
compact multis is fully representative of the mutual inclination distribution
of the underlying (intrinsic) population of multi-planet systems (including
both transiting and non-transiting planets). Studying this broader, “parent”
population gives context to the compact multis. The observed transiting
multiplicity distribution (the number of systems with $1,2,\dots N$ transiting
planets) depends on both the underlying mutual inclination and multiplicity
distributions (_Tremaine and Dong_ , 2012).

Early in the Kepler mission, _Lissauer et al._ (2011b) noted that a single
mutual inclination distribution of $\sigma_{i}\sim 1^{\circ}-2^{\circ}$ and a
constant number of planets per system did not reproduce the transiting
multiplicity distribution; there was an overabundance of singles relative to
multis, a discrepancy now known as the “Kepler dichotomy” (e.g., _Johansen et
al._ , 2012; _Ballard and Johnson_ , 2016). The dichotomy is likely
exaggerated (although not entirely explained) by biases due to planet
detection order (_Zink et al._ , 2019). The leading interpretation of the
dichotomy is that it signals the existence of a sub-population of systems with
mutual inclinations larger than a few degrees.

Several authors have constrained the joint inclination-multiplicity
distributions through forward modeling (e.g., _Mulders et al._ , 2018; _He et
al._ , 2019, 2020). Here, one generates synthetic populations of planets from
parameterized size, period, inclination, and multiplicity distributions and
then “observes” these samples with a simulated Kepler mission. The population
parameters are adjusted until the observed population agrees with the actual
Kepler census. The latest models indicate that the transiting multiplicity
distribution can be well-described by a mixture of high- and low-inclination
sub-populations with $\sigma_{i,\mathrm{low}}\approx 1^{\circ}-2^{\circ}$ and
$\sigma_{i,\mathrm{high}}\approx 30^{\circ}-65^{\circ}$ making up $\sim 60\%$
and $\sim 40\%$ of systems, respectively (_He et al._ 2019; see also _Mulders
et al._ 2018). However, the data are also consistent with a continuous (i.e.,
non-dichotomous) distribution of relatively low ($i\lesssim
5^{\circ}-10^{\circ}$) mutual inclinations with a dispersion that is
inversely-correlated with the system’s intrinsic multiplicity,
$\sigma_{i}\propto n^{\alpha}$ (where $n$ is the multiplicity and $\alpha<0$;
_He et al._ 2020; see also _Zhu et al._ 2018).

These dichotomous and continuous models reproduce the transiting multiplicity
distribution roughly equally well, on account of the degeneracy between the
intrinsic multiplicity and mutual inclinations. However, this degeneracy can
be broken with further observational inputs, making it possible to distinguish
the models. Using the statistics of Transit-Duration Variations (TDVs) of the
Kepler planet population as an additional constraint, _Millholland et al._
(2021) showed that the continuous distribution of low ($i\lesssim
5^{\circ}-10^{\circ}$) mutual inclinations (_He et al._ , 2020) is most
favored by the data. Specifically, long-term TDV signals are driven by orbital
precession between inclined planets, and they increase in magnitude with
mutual inclination. The Kepler planet population cannot contain a large
fraction of systems with large mutual inclinations, since there would have
been a larger number of TDV detections than actually observed.

Cumulatively, the data show that systems with higher intrinsic multiplicity
are dynamically cooler in both eccentricities and inclinations. The high
multiplicity compact multis appear to be the “coolest” extreme of a dynamical
continuum of close-in, multi-planet systems. In addition, if we consider the
$e=0$ and $i=0$ configurations as the lowest energy state, then the excitation
is roughly comparable in the two dynamical variables. Specifically, the degree
of inclination excitation (expressed in radians) is given by
$\sigma_{i}\approx 0.02-0.04$, which compares well with the degree of
eccentricity excitation $\sigma_{e}\approx 0.03-0.05$. In other words,
existing data are roughly consistent with equipartition in energy excitation
as measured by $(e,i)$. This finding poses an interesting avenue for future
exploration: If one considers the degree of mass uniformity, size uniformity,
or orbital spacing uniformity as additional “dynamical variables” one could
investigate whether or not excitation in these variables is correlated with
excitation in eccentricity and inclination.

###  2.7 Weak Correlation Between Stellar Properties and Architectures of
Compact Multis

As we review in §3, the planet formation process begins with the collapse of a
star and the formation of a circumstellar (or protoplanetary) accretion disk.
For nearly all known exoplanets, the circumstellar accretion disks in which
they formed are long gone. Only the host stars remain. Nonetheless, the
present-day host star properties offer a window into the environment at the
epoch of formation. The exoplanet literature contains many studies that have
investigated correlations between star and planet properties. Such studies
treated Jovian planets first because they were discovered first.

As a point of comparison for the small planets, it is noteworthy that giant
planets with $a=0.1{-}1.0$ AU are more common around massive and metal-rich
stars, and they are more eccentric than their counterparts orbiting sub-solar
hosts (_Gonzalez_ , 1997; _Santos et al._ , 2004; _Fischer and Valenti_ ,
2005; _Johnson et al._ , 2007; _Dawson and Murray-Clay_ , 2013). We expect
that both stellar mass and metallicity correlate with the total inventory of
solids in the protoplanetry disk and conclude that solid-rich disks produce
giant planets more efficiently. The prevalence of planets smaller than Neptune
only slightly increases with stellar metallicity (see, e.g., _Buchhave et al._
2012; _Petigura et al._ 2018) and decreases with stellar mass (see, e.g.,
_Howard et al._ 2012; _Mulders et al._ 2015). These trends suggest a different
relationship between disk solid inventory and the formation of small planets.

The connection between stellar metallicity and fundamental properties for
smaller planets is less clear due to the challenge of obtaining an accurate
inventory of the small planets, particularly in systems that appear to have
just one transiting planet. An analysis of transit durations by _Mills et al._
(2019) suggested a positive correlation between eccentricity and stellar
metallicity for small planets. _Dong et al._ (2018) found a separate line of
evidence for dynamically excited small planets orbiting metal-rich stars. Sub-
Neptunes with $P<10$ days are preferentially found around stars with super-
solar metallicities and are rarely members of multi-transiting systems,
suggesting a correlation between stellar metallicity and mutual inclination.
No strong correlation has been found between host star properties and the
radius and period distributions of small planets as a function of planet
multiplicity (_Xie et al._ , 2016; _Weiss et al._ , 2018b).

Here, we investigate the relationship between host star properties and the
fundamental properties of the compact multis, including a subset of systems
with very low size dispersion ($\sigma_{R}<0.1$, Figure 8). The correlations
between host star properties and the typical planet size for the most uniform
systems are weaker than for the population as a whole. This is perhaps because
we have selected a smaller sample of planets, or because the compact multis
with the most uniform sizes are slightly smaller on average than the typical
sub-Neptune. Nonetheless, the lack of strong relationships between the host
star properties and the emergent peas-in-a-pod architecture suggests that the
host star itself has little influence on the planetary patterns that form.

Note that the host stars occupy a narrow range of mass and metallicity
($0.6<M_{\star}/M_{\odot}<1.3$ and $-0.5<\mathrm{[Fe/H]}<0.5$). However, if we
were to include, e.g., TRAPPIST-1, which has a typical planet size of $\sim$1
$R_{\oplus}$ for a $0.1\,M_{\odot}$ star, this would only contribute to the
flatness of the planet radius vs. stellar mass panel, underscoring that
stellar mass is not the primary driver of the planet radii. This flatness of
planet size as a function of stellar mass was also noted in _Dai et al._
(2020), who interpreted it as evidence that planet formation is more efficient
in disks around low-mass stars.

Figure 8: Left: Mean planet radius vs. stellar mass in the CKS compact multis
with four or more planets, including a subset of planets with low size
dispersion ($\sigma_{R}<0.1$ dex). The error bars represent $\sigma_{R}$ in
units of Earth radii. Right: Same as left, but as a function of host star iron
abundance. Although the typical planet size is marginally correlated with
stellar mass and [Fe/H] (Pearson-R and p-values indicated), the scatter in
typical planet size greatly exceeds the intra-system radius dispersion,
indicating that stellar mass and metallicity are not the dominant factors in
setting the planet sizes in the peas-in-a-pod pattern.

###  2.8 Potential Biases

Before transitioning to theoretical discussions of the peas-in-a-pod patterns,
let us pause to consider the potential influence of observational biases and
selection effects. To our knowledge, there are two recent papers that have
interpreted the observed intra-system uniformity in planet sizes as the result
of observational bias.

In one study, _Zhu_ (2020) considered an alternative randomization procedure
aimed at assessing the significance of intra-system uniformity. They generated
random multi-planet systems by starting with the _Weiss et al._ (2018a) sample
of planet sizes and rescaling them according to
$R_{p}^{\prime}=R_{p,\mathrm{obs}}\,(\mathrm{SNR}_{\mathrm{rand}}/\mathrm{SNR}_{\mathrm{obs}})^{1/2}$
where $R_{p,\mathrm{obs}}$ is the observed planet radius,
$\mathrm{SNR}_{\mathrm{obs}}$ is the observed signal-to-noise ratio of the
transit and $\mathrm{SNR}_{\mathrm{rand}}$ is drawn at random from the
ensemble distribution of $\mathrm{SNR}_{\mathrm{obs}}$. _Zhu_ (2020) observed
an intra-system uniformity in radius, albeit with a weaker correlation than
reported in _Weiss et al._ (2018a), and interpreted this uniformity as
evidence that the entire peas-in-a-pod pattern was due to SNR effects.
However, _Weiss and Petigura_ (2020) showed that this procedure does not prove
that uniform radii are due to selection effects alone. If the radii are indeed
correlated, as shown in Figure 4, then adding random multiplicative scalings
will blur — but not entirely erase — the correlation. As a result, the
correlation observed by _Zhu_ (2020) does not falsify an astrophysical
explanation for the peas-in-a-pod pattern.

In another study, _Murchikova and Tremaine_ (2020) found that if the
underlying radius distribution is sufficiently steep, $dN/dR_{p}\propto
R_{p}^{-4}$, then random draws can reproduce the observed size similarity,
since planets just above the detection limit are much more common than
slightly larger planets. Although this approach provides, in principle, a
plausible way for selection effects to produce size uniformity, the assumed
power-law is inconsistent with the true distribution of planet sizes.
Measuring the size distribution of Kepler planets while accounting for
selection effects has been the subject of many independent analyses (see,
e.g., _Howard et al._ 2012, _Fressin et al._ 2013, _Petigura et al._ 2013,
_Fulton et al._ 2017, or the review by _Winn_ 2018). Smoothing over the fine
structure observed by _Fulton et al._ (2017), the radius distribution is
roughly log-uniform from $R_{p}$ = 1–3 $R_{\oplus}$, i.e., $dN/dR_{p}\propto
R_{p}^{-1}$ instead of $dN/dR_{p}\propto R_{p}^{-4}$. In addition, the
_Murchikova and Tremaine_ (2020) model predicts that most planets will be
detected near the minimum detectable SNR, which is $\sim$ 10 (_Christiansen et
al._ , 2015). However, $\sim$70% of the _Weiss et al._ (2018a) multis have SNR
$>$ 20\. In summary, the _Murchikova and Tremaine_ (2020) model can only
produce the peas-in-a-pod pattern through selection effects by invoking a
planet size distribution and a planet signal-to-noise distribution that are
inconsistent with the Kepler data.

_Murchikova and Tremaine_ (2020) also considered a second model, with a less
steep radius dependence, but with 9 free parameters, including a prescription
that the planets in a given system tend to be either $1.5$ or
$2.4\,R_{\oplus}$ (based on the peaks of the planet radius distribution from
_Fulton et al._ 2017), with a tuneable size dispersion for planets of either
radius type. This model somewhat reproduces the size similarity observed in
CKS, but it does not reproduce the other patterns described above, including
the spacing similarity and the size-spacing relation. Also, the prescription
that planets in the same system tend to be centered around $1.5$ or
$2.4\,R_{\oplus}$ is not a null hypothesis because it assumes that planets
within the same system have similar sizes.

The role of detection biases in shaping our understanding of patterns in
planetary systems is indeed important, but the case for the peas-in-a-pod
pattern having an astrophysical origin remains strong. Many different groups —
using independent methods and tools — find evidence for patterns of uniformity
(e.g. _Weiss et al._ , 2018a; _He et al._ , 2019; _Gilbert and Fabrycky_ ,
2020; _He et al._ , 2020; _Mishra et al._ , 2021; _Millholland and Winn_ ,
2021; _Otegi et al._ , 2021). Nonetheless, the question of the veracity of the
peas-in-a-pod phenomenon will ultimately be decided by the continued
accumulation of observational data.

Figure 9: Summary of physical and orbital properties of compact multis
reviewed in §2. For each planet observable we summarize the observations and
provide a sketch of a key figure from the following references: size— _Fulton
et al._ (2017) Figure 7, size dispersion— _Weiss et al._ (2018a) Figure 2,
orbital distance— _Petigura et al._ (2018) Figure 7, spacing— _Fabrycky et
al._ (2014) Figure 4, Spacing dispersion— _Weiss et al._ (2018a) Figure 7,
size-spacing—this work Figure 7, size-mass dispersion— _Millholland et al._
(2017) Figure 2, eccentricity— _Van Eylen et al._ (2019) Figure 5,
inclination— _He et al._ (2020) Figure 7.

###  2.9 Summary of Observations: Peas-in-a-Pod

Let us take stock of the physical and orbital properties of compact multi-
planet systems described in this section. Figure 9 lists a number of planet
observables and summarizes the current state of the observations. For each
item in the list, we provide a reference to a key figure in the literature. To
summarize: The prevalence of planets increases with decreasing size, with a
gap in the distribution at $R_{p}=1.5-2.0$ $R_{\oplus}$. Compact multis have
correlated sizes and masses. The occurrence of small planets is roughly log-
uniform down to $P\approx 10$ days, and below that value it falls
precipitously. Planetary pairs have a relatively broad distribution of orbital
spacing with slight over-abundances wide of first-order mean-motion
resonances. Despite the overall range in observed spacing parameters, however,
the spacings within a given system are more highly correlated. The orbital
spacing distribution within a planetary system depends on the size of the
constituent planets, the smallest planets preferring tighter spacing. Compact
multis have low eccentricities and inclinations. Finally, the host star masses
and metallicities do not strongly correlate with the peas-in-a-pod
architectural features.

Taken together, these findings indicate that planets in compact multi-planet
systems are dynamically cool. The planets have somehow coordinated their
sizes, masses, and orbital spacings in a manner that depends only weakly on
the properties of the host star. In the following sections, we consider the
physical mechanisms that can account for these features. We explore how the
formation of the star-disk system (§3) and the emergence of planets in the
disk (§4) sets the characteristic mass of the planets. In the final stages of
planet formation and evolution, planet-planet interactions dominate (§5), and
we consider how the dynamics in mature planetary systems sculpt both the
regularity of the peas-in-a-pod pattern and the conditions for departures from
that pattern. Figure 10 provides a schematic roadmap for our review of these
physical processes.

##  3 Connection to Star Formation and Circumstellar Disk Properties

Planets form within circumstellar disks, which are naturally produced
alongside their host stars. As a result, the star formation process places
important constraints on the subsequent process of planet formation, including
the compact multi-planet systems of interest here. More specifically, the
properties of disks, which provide the initial conditions of planet formation,
are sculpted by the earlier action of star formation. Here we briefly review
these constraints — including the mass, radius, and surface density
distributions of the disks — with a focus on properties relevant to the
emergence of the peas-in-a-pod pattern and other observed properties of the
compact multis (for a more comprehensive treatment, see the Chapters by
Miotello et al. and by Manara et al.). We consider these properties at the
epoch corresponding to the end of star formation (§3.1) and over the
subsequent epochs when planets form (§3.2, §3.3). This span of time thus
covers the transition from protostars to planets.

###  3.1 Constraints from Star Formation

Star formation itself occurs within molecular clouds that can include millions
of solar masses of material (_Shu et al._ , 1987). However, individual star-
forming events take place within molecular clouds cores, which are much
smaller sub-units containing several times the mass of the stars that they
produce. At intermediate scales, stars form within embedded clusters that
contain 100s to 1000s of stars (_Lada and Lada_ , 2003). While these cluster
systems can sculpt planetary systems after their formation (_Adams_ , 2010),
this background environment is subdominant during the collapse phase that
actually produces a star/disk system.

Figure 10: Schematic diagram showing the various physical processes involved
in the formation of compact multi-planet systems and the possible emergence of
peas-in-a-pod architectures. These processes can be broadly divided into three
categories: disk properties, planet-disk effects, and planet-planet effects.
Disk properties are set by the formation of the star-disk system, and they
include the global characteristics and evolution of the accretion disk.
Planet-disk effects include the conglomeration of dust and rocky material in
larger bodies, accretion of nebular gas onto the planet (and subsequent gas
loss), and migration driven by planet-disk interactions. Lastly, planet-planet
effects include any pairwise interactions, such as gravitational effects
leading to enhanced long-term stability. These three stages cover the full
planet formation process roughly temporally, forming the roadmap for §3, §4,
and §5, respectively. Finally, all of these stages will be explored together
in the population synthesis models of §6.

We consider a simple model for the collapse of the molecular cloud core, from
which a few key physical scales of the planet formation environment emerge. We
can thus describe the molecular cloud core with two variables: the sound speed
$v_{s}$ sets the degree of thermal support and the rotation rate $\Omega$
determines the amount of angular momentum (e.g., see the review of _Shu et
al._ 1987). From these quantities, we can specify the maximum mass available

$M_{\rm max}={v_{s}^{3}\over G\Omega}\,,$ (7)

the size of the region for an individual star formation event

$R_{\rm core}=v_{s}/\Omega\,,$ (8)

and size of the disk $R_{\rm disk}$ as determined by the centrifugal barrier
$R_{C}$ in the collapse flow

$R_{\rm disk}\sim R_{C}={G^{3}M^{3}\Omega^{2}\over 16v_{s}^{8}}\,.$ (9)

This radius should be considered as the starting radius for a disk when a
total mass $M$ has fallen to form the star/disk system. Under the action of
viscosity, subsequent disk evolution will spread out the disk, so that the
outer disk radius is expected to grow.

The collapse conditions for star formation define a mass infall rate
${\dot{M}}\sim v_{s}^{3}/G$. As a result, the formation of a star with mass
$M_{\star}$ must take place over a time scale $t_{\star}\sim
M_{\star}/{\dot{M}}\sim GM_{\star}/v_{s}^{3}\sim 0.1-0.3$ Myr.

In the collapse flow that produces stars, nearly all of the mass falls to
radii much larger than a stellar radius, i.e., material fall first onto the
disk. Gravitational instabilities in the disk (along with other mechanisms)
act to transfer angular momentum and facilitate accretion of disk material
onto the star. Gravitational instabilities grow on a dynamical time scale when
the parameter $Q=v_{s}\Omega/\pi G\Sigma\sim 1$ (_Toomre_ , 1964), which
results when the disk mass is comparable to the stellar mass. In order for the
disk instability time scale to become longer than the collapse time of a cloud
core (about 0.1 Myr), the disk needs $Q\sim 10$ which corresponds to a disk
mass

$M_{\rm disk}\sim{1\over 10}M_{\star}\,.$ (10)

Note that this mass scale for the disk is the starting disk mass and/or the
maximum disk mass at birth. Observed disk masses have an upper envelope
bounded by this value (e.g., _Hartmann_ 2008). In other words, the data show a
range of disk masses for a given stellar mass, but the maximum observed disk
mass scales according to equation (10). More specifically, the range of disk
masses for a given stellar mass spans an order of magnitude (_Andrews et al._
, 2013).

For completeness, we note that the disk must have an approximate minimum mass
for a given disk radius. The initial molecular cloud cores that form stars
contain more angular momentum than can be carried by a star (the star would
have to rotate much faster than its breakup angular velocity to contain that
angular momentum – this finding represents the well-known angular momentum
problem of star formation). As a result, the disk must initially contain most
of the angular momentum of the initial cloud core (note that subsequent
outflows can remove angular momentum). In approximate terms, the disk angular
momentum is given by $J\sim M_{\rm disk}\sqrt{GM_{\star}R_{\rm disk}}$, where
$R_{\rm disk}$ is the outer disk radius and where we ignore dimensionless
factors of order unity. For typical stellar masses $M_{\star}\sim 1M_{\odot}$
and observed disk radii $R_{\rm disk}\sim 100$ AU, the disk must contain
roughly $\sim 0.01M_{\star}$ in order to carry the angular momentum. Again,
these values represent the disk properties at the end of star formation, and
the beginning of planet formation (coincident with the T Tauri phase of
stellar evolution). As the system evolves, the disk mass decreases due to
accretion, evaporation, and planet formation. Accretion causes the disk radius
to grow, evaporation carries away some of the angular momentum with the flow,
and the remaining angular momentum is locked up in the planets (primarily in
their orbits, with some accounting for planetary rotation).

The maximum and minimum disk masses discussed above indicate that the starting
mass available for planet formation falls in the range $0.01-0.1M_{\star}$. If
we consider typical stars with mass $M_{\star}=M_{\odot}$ and typical
metallicity $Z=0.01$ (near solar), the available mass in solids thus falls in
the range $M_{\rm solid}=30-300$ $M_{\oplus}$ (e.g., see _Tychoniec et al._
2018 for observational constraints). Given that compact multis typically
contain at least 4 planets with $\sim 10M_{\oplus}$ of rocky material, the
heavy metal inventory must be greater than about 0.1 of the solar value to
produce such typical systems. In other words, systems with sufficiently low
$Z$ are predicted to not have multiple $\sim 10M_{\oplus}$ planets. Systems
with moderately lower $Z$ can in principle produce systems with lower mass
planets, but the planet masses must vanish in the limit $Z\to 0$.

In addition to the total mass of solids in a disk, the radial distribution of
solids is relevant for planet formation. When disks are built from the infall-
collapse flow that produces their stellar hosts, the surface density
distribution has a nearly power-law form

$\Sigma(r)=\Sigma_{0}\left({r\over r_{0}}\right)^{-p}\,,$ (11)

where $r_{0}$ is a reference orbital radius and the coefficient $\Sigma_{0}$
is determined by the total mass of the disk. The power-law index depends on
the manner in which the disk forms. If the starting angular momentum profile
of the core corresponds to an isothermal sphere with uniform rotation, and if
the material does not change radius as it become incorporated into the disk,
then the power-law index $p=5/3$. However, the incoming material does not, in
general, have the correct azimuthal velocity to become part of a Keplerian
disk. If one assumes that parcels of gas lose energy, at constant angular
momentum, as they join the disk, then the power-law index $p\approx 7/4$
(_Cassen and Moosman_ 1981; _Terebey et al._ 1984; see also _Adams and Shu_
1986). For comparison, the well known Minimum Mass Solar Nebula (MMSN) also
has a power-law form (starting with _Hayashi_ 1981). This profile, obtained by
augmenting the observed masses of solar system planets and by using their
current orbits, has power-law index $p=3/2$.

Observations of circumstellar disks tell a similar story, but suggest some
differences. A full review of the data is beyond the scope of this
contribution (again, see the Chapters by Miotello et al. and Manara et al.),
but the bottom line is that observed surface density estimates also show a
power-law form, in agreement with theoretical considerations and estimates
made from the exoplanet sample. On the other hand, observed disks tend to have
shallower density profiles – smaller indices – where current SMA and ALMA maps
suggest that $p\approx 1$ (e.g., _Andrews et al._ 2009).

Finally, we note that various attempts have been made to use the sample of
observed exoplanets to construct a minimum-mass extrasolar nebula. These
analyses also generally give power-law forms with varying indices and
normalizations. Using all of the planets available at the time, _Chiang and
Laughlin_ (2013a) find a power-law index $p=1.6$ for both the gaseous disk
component and the solids. Although this form is widely used, it is
inconsistent with the peas-in-a-pod configuration (note that a purely peas-in-
a-pod system, which has equal mass in equal units of $\mathrm{dln}a$, has
$p=2$, whereas the energy optimized configuration of _Adams_ 2019 has $p$ =
11/6; see §5.2). Earlier work (_Kuchner_ , 2004) find similar results with
$p\sim 3/2$, consistent with the minimum mass solar nebula, but with large
uncertainties (due to the smaller amount of data available at the time).

The above collection of power-law indices for circumstellar disks can be
summarized as follows:

Source | $p$
---|---
ALMA observations | $\sim 1$
Minimum Mass Solar Nebula | $3/2$
Minimum Mass Extrasolar Nebula | $\sim 3/2$
Star formation theory | $5/3-7/4$
Pair-wise energy optimization | $11/6$
Pure peas-in-a-pod configuration | $2$

Both circumstellar disks and planetary systems are observed to have nearly
power-law surface density distributions, with indices falling in the range
$p=1-2$. Disk configurations of this general form are produced naturally from
the collapse flow that forms stars. Moreover, these results are consistent
with known (expected) properties of the angular momentum profiles of the
molecular cloud cores that provide the initial conditions. In spite of the
apparent consistency, one important discrepancy remains. The peas-in-a-pod
pattern has a steeper surface density distribution than the disks that gave
rise to it (as well as other planetary systems). This difference suggests that
either the raw materials for planet formation, or the planets themselves, must
move relative to the positions indicated by the initial disk profiles.

With the disk masses, outer radii, and surface density profiles specified, we
now consider the inner boundary condition. The inner disk edge is controlled
by the presence of strong magnetic fields generated within the star. The
circumstellar disks are truncated at the radius $R_{X}$ where the inward
pressure due to the accretion flow is balanced by the outward pressure from
stellar magnetic fields. For dipole field structures, the radius $R_{X}$ can
be written in the form

$R_{X}=\zeta_{X}\left({B_{\star}^{4}R_{\star}^{12}\over
GM_{\star}{\dot{M}}^{2}}\right)^{1/7}\,,$ (12)

where ${\dot{M}}$ is the mass accretion rate through the disk and $B_{\star}$
is the magnetic field strength on the stellar surface. The dimensionless
parameter $\zeta_{X}$ is expected to be of order unity and depends on the
model assumptions (_Ghosh and Lamb_ , 1978; _Blandford and Payne_ , 1982). The
magnetic field strengths have been measured for a collection of T Tauri stars
(_Johns-Krull_ , 2007) and fall in the range $B_{\star}\approx 1-2$ kilogauss.
As a result, the truncation radii that define the inner disk edge have typical
values $R_{X}\sim 0.1$ AU. Note that this size scale is comparable to the
semimajor axes of the observed compact multis (where the latter are observed
at much later ages). Notice also that as the mass accretion rate decreases,
the truncation radius $R_{X}$ increases, so that this coincidence only holds
at early epochs of disk evolution (ages of few Myr).

One thus obtains the ordering of radial scales

$R_{\rm core}\gg R_{\rm disk}\sim R_{C}\gg R_{X}>R_{\star}$ (13)

and the corresponding ordering of mass scales

$M_{\rm core}>M_{\star}\gg M_{\rm disk}\gg M_{p}\,,$ (14)

where $M_{p}$ is the mass scale of planets ($\sim 10$ $M_{\oplus}$ for the
typical planet in a compact multi, of interest here). Another defining
characteristic of the disks is that they are geometrically thin, so that the
scale height of the disk much smaller than the radius, i.e., $h/r\ll 1$. Since
the scale height ratio is equivalent to the ratio of speeds, $h/r=v_{s}/v_{\rm
orb}$, the disks are cold. These scales and their ordering sets the stage for
the subsequent process of planet formation.666In spite of the disks being
geometrically thin, their vertical extent is far greater than the radii of
planets, $h\gg R_{p}$.

###  3.2 Evolution of Disks, Stars, and Clusters

The disk properties outlined above result from the process of star formation
and correspond to the time when the star is first formed, i.e., the epoch when
the star/disk system has gathered the majority of its final mass. The time
scale for this process, and hence the system age when the disk has the
aforementioned properties, is of order 0.1 Myr. Significantly, this collapse
time scale is much shorter than several other scales of interest. Estimates
for the expected time required for planet formation vary widely, but extend up
to 100 Myr for the formation of Earth in our solar system. The observed
lifetimes of circumstellar disks typically fall in the range 1 – 3 Myr. The
pre-main-sequence contraction time for stars of solar mass and smaller are
also of order 10 Myr. Finally, the lifetimes of embedded stellar clusters are
$\sim 10$ Myr, with the more robust open clusters lasting $\sim 100$ Myr. All
of these astrophysical processes taking place in the background can influence
planet formation, as outlined below.

Perhaps most importantly, lifetimes for circumstellar disks are relatively
short. Several observational surveys have studied disk signatures, such as
infrared excess emission and veiled emission lines, in star forming regions
with measured ages (_Hernández et al._ , 2007; _Mamajek_ , 2009; _Fedele et
al._ , 2010). The fraction of stars with disk signatures is found to be a
decreasing function of time with an approximately exponential form. Roughly
half of the disks (or at least their observational signatures) are gone by an
age of 3 Myr, with only $\sim 10\%$ of the systems retaining disks at 10 Myr
(_Meyer et al._ , 2007; _Williams and Cieza_ , 2011; _Hartmann et al._ , 2016;
_Manzo-Martínez et al._ , 2020). These observations correspond to the presence
of gas and accretion processes in the disks, and thus indicate that gas is
only available for planet formation over 3 – 10 Myr. On the other hand, for
the production of smaller rocky planets, with little or no gas, this time
constraint is less problematic. The compact multis fall in an intermediate
regime. Although they are thought to have a primarily rocky composition, some
of the planets retain a hydrogen/helium atmosphere (_Weiss and Marcy_ , 2014;
_Rogers_ , 2015). Although the gaseous mass is small, only a few percent of
the total, its presence indicates that the planets are forming before the disk
fully dissipates. On the other hand, if large rocky planets form while the
disk retains a large gas supply, the planets can readily accrete thicker
atmospheres than indicated by observations (for further discussion and
possible solutions, see _Lee et al._ 2014; _Ormel et al._ 2015; _Lambrechts
and Lega_ 2017).

Disks evolve through a variety of processes, including viscous accretion,
evaporation, and disk winds. The latter can be driven by radiation fields from
the central star itself (_Owen and Jackson_ , 2012) and/or by external sources
(_Adams et al._ , 2004). In addition, some mass becomes locked up within
forming planets. All of these effects combine to make the disk mass decrease
with time. In addition, viscosity acts as a diffusive process (_Shu_ , 1992)
so that the disk surface density is expected to spread out with time. In other
words, the natural evolution of the disk is to become less centrally condensed
with time, whereas the peas-in-a-pod pattern corresponds to a state of greater
concentration. This mismatch indicates that the rocky material that makes up
planets must move (migrate) relative to the gas that comprises most the disk
mass. How this redistribution takes place represents an important unresolved
issue. Note that this movement can take place at different stage of planetary
development, from early times when the rocky material is still small (rocks)
to late times after the planets have already finished fully forming.
Nonetheless, by an age of $\sim 5$ Myr, the outskirts of disks are depleted in
dust mass, suggesting that the mass has already been consolidated into pebbles
or planetary-sized bodies or migrated inward (_Ansdell et al._ , 2015).

As the disks spread out and lose their mass, and planets form within them, the
host star evolves on a similar time scale. During the time in which planetary
systems are forged, the central star contracts in radius and decreases its
luminosity. Over the same time, it slows down its rotation rate and decreases
its quadrupole moment $J_{2}$, thereby changing the manner in which the star
couples dynamically to the disk and/or planets.

###  3.3 Global Angular Momentum Constraints

Circumstellar disk properties are determined to a large extent by conservation
of angular momentum, as discussed in §3.1. The resulting configurations have
roughly solar system sizes ($R_{\rm disk}\sim 100$ AU) and moderate masses
(starting at $M_{\rm disk}\sim 0.1M_{\star}$). To constrain how such star/disk
systems subsequently evolve, consider the following question: What is the
lowest energy state accessible to a star/disk system subject to conservation
of angular momentum? The optimum energy is achieved by placing essentially all
of the mass in the central object, while leaving one small particle with a
large orbit to carry the angular momentum. Note that this optimum energy state
is nearly realized by our solar system: Essentially all of the mass is
contained in the Sun, with the majority of the angular momentum carried by
Jupiter in its orbit. Since any mechanism for energy dissipation causes
star/disk systems to evolve toward a similar state, exoplanetary systems are
expected to have similar properties.

Disks can lose angular momentum through the action of photoevaporative winds,
which carry away gaseous material. If the solids remain in small entities,
they can be swept away with the gas. If, instead, the solids grow to size
scales large enough to decouple them from the gas, they remain behind as the
gas evaporates. As a result, stars must have ‘solar systems’ consisting of
some rocky material. Disks typically contain a mass $M_{\rm
solid}=30-300M_{\oplus}$ of solid material (§3.1). Suppose that a fraction
${\cal F}$ of the original solid content of the disk is left behind. Even if
all of the gas can be successfully evaporated, the disk (planetary system)
will contain a mass ${\cal F}M_{\rm solid}$. If the evaporating gas does not
carry away the specific angular momentum of the rocky material, the remaining
system will have characteristic size $R\sim R_{C}/{\cal F}^{2}$ (where $R_{C}$
is given by equation [9]).

Compact multis typically consist of 3 or 4 detected planets, each with mass
$\sim 10M_{\oplus}$ and orbital periods less than $\sim 100$ days. The total
mass in rocky material is thus a significant fraction of the expected mass
$M_{\rm solid}$ from the initial state. Since the semimajor axes are small,
however, the angular momentum of the system is only a fraction
$\sim(a/R_{C})\sim 1/10$ of the total expected value. This discrepancy implies
that the compact multis must either have additional material in orbit at
larger radii, or have been formed with highly efficient removal of angular
momentum.

##  4 Formation of Sub-Neptune-Sized Planets

The formation of the star-disk system sets the global conditions in which
planet formation occurs, as detailed in the previous section. Within this
environment, dust and rocky material coalesces into progressively larger
bodies, eventually generating proto-planets (recall the schematic outline in
Figure 10). The specific processes by which this growth occurs is a topic of
vigorous research. Nevertheless, numerous new ideas pertaining to the
formation of sub-Neptune-sized planets offer insights into how planets might
emerge at a variety of physical scales. Here we briefly review the detailed
formation processes with an eye towards quantifying the characteristic planet
mass and size scale (for a more general treatment, see the Chapter by
Drazkowska et al.). We proceed temporally, focusing first on the growth of
plantesimals (§4.1) and then planets (§4.2), before considering the role of
orbital migration (§4.3) and atmospheric loss (§4.4).

###  4.1 From Dust to Planetesimals

The first step towards quantifying how planetary objects coalesce lies in
understanding the genesis of $\sim 10-100\,$km planetary building blocks,
i.e., the planetesimals within the protoplanetary disk. Rocky and icy grains
within the nebulae originate as $\sim\mu$m dust grains that are delivered to
the circumstellar disks along with the infalling gas (§3). Like snowflakes in
the Earth’s atmosphere, these grains gradually grow to become $\sim$mm-cm
“pebbles". Their growth, however, does not continue beyond this characteristic
size scale due to the so-called fragmentation barrier. Both modeling and
experimental studies have shown that beyond a characteristic particle size on
the order of a centimeter (millimeter), even low-velocity collisions among icy
(rocky) particles lead to breakup or mass-transfer among impactors (_Blum and
Münch_ , 1993; _Windmark et al._ , 2012). Accordingly, the planetesimal
formation process must boil down to the recurrent and rapid conversion of
innumerable number of pebbles into massive asteroid-like bodies.

A number of empirical lines of evidence pertaining to the surviving small-body
populations of the solar system – including the size-distribution of the
asteroid belt (_Morbidelli et al._ , 2009) as well as the inclination
distribution of binaries within the Kuiper belt (_Nesvorný et al._ , 2019) –
indicate that the physical machinery that drives planetesimal formation is
driven by direct gravitational collapse. In other words, pebbles within the
protoplanetary disk concentrate into vast clouds that grow massive enough to
become gravitationally unstable. In turn, concentration of pebbles into these
clouds is facilitated by a distinct hydrodynamical process – the streaming
instability (_Youdin and Goodman_ , 2005; _Johansen et al._ , 2007).

Although the detailed picture of the streaming instability (and resonant drag
instabilities in general) is intricate (_Squire and Hopkins_ , 2018), its
basic physics can be understood in a straightforward manner. As a starting
point, recall that due to internal pressure support, the orbital velocity of
gas within the circumstellar disk lags the Keplerian speed by a small margin
$\Delta v=\eta_{w}\,v_{\rm{kep}}\sim(h/r)^{2}\,v_{\rm{kep}}$. Due to
interactions with a sub-Keplerian nebular flow, dust grains within
protoplanetary disks gradually drift inwards. Among other factors, their drift
rate is controlled by the ratio of the particle cross-sectional area to the
mass. If we consider the evolution of two nearby dust grains in a
protoplanetary disk, their cumulative drift rate can be diminished if one
particle “hides" in the hydrodynamic wake of the other, such that the mass of
the two-particle aggregate is doubled while leaving its effective cross-
sectional area unchanged. The differential drift among solids created by this
process allows the two-particle aggregate to encounter and capture a third
particle, and so on. Extending this thought experiment to an entire cloud of
pebbles, one can understand how a sufficiently massive peloton of icy and
rocky grains can locally accelerate the gas towards the Keplerian value,
thereby creating a broadly favorable environment for continued capture of
drifting particles and facilitating its own continued growth.

The purely hydrodynamic accumulation of solids into massive filaments
facilitated by the streaming instability cannot continue indefinitely, and the
nonlinear evolution of pebble clouds within protoplanetary nebulae culminates
in their gravitational collapse. Drawing upon the pronounced analogy of this
process with the formation of stars, _Klahr and Schreiber_ (2021) derived an
effective Jeans mass for planetesimals by considering the competition between
the timescales for pebble turbulent diffusion and gravitational collapse. This
relation sets the characteristic mass scale for bodies generated from self-
gravitational collapse:

$\displaystyle m\approx
100\,\bigg{(}\frac{\alpha}{\mathcal{Z}_{\rm{Hill}}}\bigg{)}^{3/2}\,\bigg{(}\frac{h}{r}\bigg{)}^{3}\,M_{\star},$
(15)

where $\alpha\sim 10^{-4}-10^{-2}$ denotes the turbulent diffusivity
parameter, $\mathcal{Z}_{\rm{Hill}}\sim 10-100$ is the mid-plane dust-to-gas
ratio upon reaching Hill density, and $h/r\gtrsim 10^{-2}$ is the disk’s
geometric aspect ratio.

_Liu et al._ (2020) took a different approach to computing $m$, and obtained a
characteristic planetesimal mass by extrapolating the results of numerical
streaming instabiity simulations. For their preferred fit parameters, the
relevant expression takes the form:

$\displaystyle m\approx
200\,\big{(}h\,r\,\Sigma\big{)}^{3/2}\,\sqrt{\frac{Z}{M_{\star}}}\,\frac{M_{\oplus}}{M_{\odot}},$
(16)

Given nominal disk properties, the above expressions yield characteristic
planetesimal radii on the order of one to a few hundred kilometers, with the
analytical estimate systematically overestimating the numerical results.
Nonetheless, the emerging paradigm suggests that the formation of planetary
building blocks stems from a near-universal gravity-hydrodynamic mechanism and
instills a pronounced degree of mass uniformity onto newly born planetesimals.

###  4.2 From Planetesimals to Planets

Planetesimals generated through gravitational instabilities of pebble clouds
can perpetuate their growth through two distinct pathways: pairwise collisions
with other planetesimals, and/or continued accumulation of small rocky bodies.
At a basic level, both of these processes can be understood from analytic
considerations, and the characteristic rate of mass-accretion can be derived
from simple $n$-$\sigma$-$v$ type relations. We begin by reviewing the former
process.

#### Pairwise Collisions.

To leading order, the collisional mass accretion rate can be approximated as
the product of the rate of collisions $\Gamma$ and the impactor mass $m$ such
that $\dot{M}\sim m\,\Gamma$. Replacing the number density with the
planetesimal surface density $n\,m\sim\Sigma_{\rm{pl}}/h$, and approximating
the velocity dispersion as $\langle v\rangle\sim(h/r)\,v_{\rm{kep}}$, the mass
accretion rate takes the form (_Lissauer_ , 1993; _Kokubo and Ida_ , 1996):

$\displaystyle\dot{M}=\Sigma_{\rm{pl}}\,\pi\,R^{2}\,(1+\Theta)\,\Omega.$ (17)

Note that in the above expression, the physical cross-section for collisions
is augmented by the gravitational focusing factor $(1+\Theta)$, where
$\Theta=(v_{\rm{esc}}/\langle v\rangle)^{2}$ is the Safronov number. This
factor accounts for the fact that in the limit of very low velocity
dispersion, the Safronov number approaches $\Theta\rightarrow\infty$, and
accretion ensues at an accelerated rate due to enhanced gravitational
focusing. Alternatively, if the velocity dispersion of the planetesimal swarm
becomes comparable to the escape velocity of the constituent planetesimals,
the collision cross-section becomes comparable to the physical radius.

For a given value of $\Theta$ (with the caveat that $\Theta$ does not exceed
unity by a large margin), it is convenient to express $\dot{M}$ as
$4\,\pi\,\rho_{\star}\,R^{2}\,dR/dt$, and note that the dependence on the
physical radius cancels, yielding a constant $dR/dt$ (equivalently, $M\propto
t^{3}$). Taking advantage of this simplification, we can rewrite equation (17)
as an accretion timescale. The time required for the radius to evolve from
$R_{0}$ to $R$ is then:

$\displaystyle\mathcal{T}_{\rm{accr}}=\frac{4\,\rho_{\star}\,(R-R_{0})}{\Sigma_{\rm{pl}}\,\Omega\,(1+\Theta)}.$
(18)

Before evaluating $\mathcal{T}_{\rm{accr}}$, we remark that collisional growth
cannot persist indefinitely, and stalls when the planetary mass reaches the
isolation threshold (_Lissauer_ , 1993; _Kokubo and Ida_ , 1996). This
argument hinges on the fact that the planetesimal feeding zone available to
the planet only extends over a few Hill spheres in either radial direction,
i.e., $\Delta a\sim 3\,a\,(M/(3M_{\star}))^{1/3}$. Consequently, isolation is
reached when the cumulative planetesimal mass contained within an an annulus
of radial extent $2\,\Delta a$ matches the planetary mass itself:
$M=4\,\pi\,a\,\Sigma_{\rm{pl}}\,\Delta a$. Rearranging this expression for the
isolation mass, we have:

$\displaystyle
M_{\rm{iso}}=\xi\,\pi^{3/2}\,\frac{\Sigma_{\rm{pl}}^{3/2}}{\sqrt{M_{\star}}}\,a^{3},$
(19)

where $\xi\approx 30$ is a dimensionless constant. Figure 11 shows the
isolation mass $M_{\rm{iso}}$ as a function of radial distance for a series of
surface density profiles, normalized to an MMSN-like value of the solid
surface density $\Sigma_{0}=20\,$g/cm2 at 1 AU. An immediately evident result
is that the isolation mass itself can be rather uniform if the solid surface
density power law index approaches $p=2$. As importantly, if we assume that
planetesimal formation can efficiently operate significantly interior to
$1\,$AU, then the timescale needed to reach isolation in the inner-most
regions of the nebula is on the order of $\mathcal{T}_{\rm{accr}}\sim
10^{5}\,$years, although this value increases to $\sim 10\,$Myr at $1\,$AU.
More strikingly, however, the value of $M_{\rm{iso}}$ is on the order of the
mass of Mars, and is thus significantly smaller than the typical mass-scale of
short-period extrasolar planets.

The above argument indicates that pairwise-collisional mode of accretion
within a MMSN-like profile of solids does not naturally yield multi-Earth-mass
objects. Instead, such a scenario necessitates a post-isolation phase of giant
impacts to consolidate the emergent embryos into bonafide planets. Of course,
a Mars-like value of $M_{\rm{iso}}$ is not universal: one may invoke a
substantially more solid-rich (e.g., MMEN-type; _Chiang and Laughlin_ 2013b)
disk, and rely on the $M_{\rm{iso}}\propto\Sigma_{\rm{pl}}^{3/2}$
proportionality to boost the isolation scale (thereby bringing the planet
formation paradigm closer to the standard model of giant-planet satellite
formation; _Canup and Ward_ 2002; _Batygin and Morbidelli_ 2020). Yet another
alternative is that of a disk with a markedly non-uniform distribution of
solids, where $\Sigma_{\rm{pl}}$ can achieve high values at discrete locations
while maintaining a low overall solid mass-budget within the disk. To this
end, recent work has shown how radial concentration of solids within
protoplanetary nebulae can lead to the formation of planetesimals within dense
narrow annuli (_Drążkowska et al._ , 2016; _Morbidelli et al._ , 2022).
Although such a scenario has gained traction within the solar system formation
literature (see e.g., _Izidoro et al._ 2021 and the references therein), the
viability of collisional accretion of mass-uniform sub-Jovian planets within
rings of solid debris remains to be demonstrated quantitatively.

#### Pebble Accretion.

A different mode of planetary growth proceeds through the capture of mm- to
cm-size particles by a growing body — a regime of conglomeration known as
pebble accretion (_Ormel and Klahr_ , 2010; _Lambrechts and Johansen_ , 2012).
Counter-intuitively, under certain conditions this process can lead to a mass
accretion rate that is considerably faster than that given by pairwise
collisions, even though the mass of each captured particle is very small in
comparison with the mass of a planetesimal.

The most important distinction between planetesimal accretion and pebble
accretion arises from the accretion cross-section. Particularly, in the case
of pebble accretion, the physical radius (augmented by gravitational focusing)
present in equation (17) is replaced by a characteristic length-scale over
which a drifting particle experiences a large-angle deflection by the
planetary gravity. For low planetary masses, this length scale is given by the
modified Bondi radius (that accounts for the fact that smaller particles are
more tightly coupled to the gas and are therefore difficult to capture; see
_Lambrechts and Johansen_ 2012). For larger planets, the critical impact
parameter for accretion becomes comparable to the Hill radius.

If vertical settling of solids is efficient (such that pebbles form a thin
sub-disk), accretion can proceed in 2D, and in analogy with equation (17), the
rate of mass accretion in the Hill regime can be written as

$\displaystyle\dot{M}_{\mathrm{2D}}\sim
2\,\Sigma_{\rm{pebble}}\,R_{\rm{Hill}}^{2}\,\tau^{2/3}\,\Omega,$ (20)

where $\tau=\Omega\,(\rho_{\rm{pebble}}/\rho_{\rm{gas}})(s/v_{\rm{th}})$ is
the Stokes number (dimensionless stopping time) of the particles. Given that
the Hill radius exceeds the planetary physical radius by a factor of order
$\sim a/R_{\star}$, it is easy to see how under the favorable assumption of
$\tau\sim 1$, the process of pebble accretion can facilitate the formation of
multi-Earth-mass bodies on a timescale that is much shorter than the age of
the protoplanetary disk (e.g., $\sim 10^{4}$ years; _Lambrechts and Johansen_
2012; _Morbidelli and Nesvorny_ 2012). Moreover, unlike the case of
planetesimal accretion where the local mass budget plays a limiting role,
pebble accretion is facilitated by the radial drift of solids, and can
therefore “access" the mass reservoirs that initially reside in the outer disk
region and continuously drift inward to cross the planet’s orbit.

Although the 2D regime of gravito-hydrodynamic dust capture can in principle
be staggeringly efficient, the advantages of pebble accretion largely fade if
the solid component of the disk is vertically well-mixed with the gas. In this
case, accretion proceeds in the 3D regime, at an exponential, but slow rate
(_Ormel_ , 2017):

$\displaystyle\dot{M}_{\mathrm{3D}}\sim
6\,\pi\,\rho_{\rm{pebble}}\,R_{\rm{Hill}}^{3}\,\tau\,\Omega,$ (21)

where
$\rho_{\rm{pebble}}\approx\Sigma_{\rm{pebble}}/(\sqrt{2\,\pi}\,h_{\rm{pebble}})$
and $h_{\rm{pebble}}$ is the dust sub-disk’s scale-height. If we adopt a dust
particle radius on the order of a millimeter (in accordance with the
fragmentation barrier; _Windmark et al._ 2012), the Stokes number evaluates to
$\tau\lesssim 0.001$ in the inner disk. Given that $h_{\rm{pebble}}\sim
h\,\sqrt{\alpha/\tau}$ (_Chiang and Youdin_ , 2010), a nebular turbulence
parameter that is substantially lower than $\alpha\ll 0.001$ is required to
break the system out of the comparatively inefficient 3D regime. Thus,
adopting an oft-quoted radial pebble flux of $F_{\rm{pebble}}\sim
10^{-4}\,M_{\oplus}/$year (which translates to a surface density of
$\Sigma_{\rm{pebble}}\sim
F_{\rm{pebble}}/(4\,\pi\,r\,v_{\rm{kep}}\,\eta_{w}\,\tau)$; _Lambrechts et
al._ 2019) and setting $h_{\rm{pebble}}\sim h$, we obtain a mass $e$-folding
timescale of
$(M/\dot{M}_{\rm{3D}})\sim\sqrt{8\,\pi}\,(M_{\star}/F_{\rm{pebble}})\,(h/r)^{3}\sim
6\times 10^{8}\,$years. This timescale is much longer than the typical
lifetime of a protoplanetary disk, and unlike 2D accretion – which is expected
to ensue beyond the ice-line where particles can attain much larger Stokes
numbers (_Batygin and Morbidelli_ , 2022) – 3D accretion in the inner disk is
unlikely to drive significant planetary growth.

Independent of the accretion mode, an important feature of the pebble
accretion paradigm is that the protoplanetary agglomeration process itself is
self-limiting (_Lambrechts et al._ , 2014). This aspect of the theoretical
picture can be understood by considering the gravitational back-reaction of
the growing proto-planet upon the gaseous disk. In the vicinity of a
sufficiently massive planet, planetary gravity can locally accelerate the gas
above the Keplerian velocity, reversing the action of gas-drag. When this
happens, pebbles circumvent the planet, and the very effect that facilitates
pebble accretion subsides.

To estimate the characteristic planetary mass at which accretion subsides we
equate the Hill radius to the disk scale-height to obtain (_Ormel_ , 2017):

$\displaystyle\frac{M_{\mathrm{iso}}}{M_{\star}}=\zeta\bigg{(}\frac{h}{r}\bigg{)}^{3},$
(22)

where $\zeta$ is a numerical factor of order unity. Importantly, numerical
experiments of _Lambrechts et al._ (2014) and _Bitsch et al._ (2015) recover
the cubic scaling of the mass ratio upon the disk’s geometric aspect ratio and
sharpen the above estimate. In particular, for a flared disk with an aspect
ratio that scales as $h/r\propto r^{1/4}$, _Lambrechts et al._ (2014) argue
for an isolation mass of $M_{\mathrm{iso}}\sim
20\,M_{\oplus}\,(a/5\,\mathrm{AU})^{3/4}$. This estimate is shown on Figure 11
as a black line.

A steep dependence on $h/r$ – a quantity that can vary by a factor of $\sim$2
from system to system – implies that pebble accretion can generate a broad
range of planetary properties, while preserving a uniform mass-scale within a
single system. Within the same framework, the stellar mass plays a sub-
dominant role. Nevertheless, equation (22) clearly predicts that despite
inherent variability, planets that orbit less massive stars should also be
less massive on average.

Taken at face-value, this result is attractive as an explanation for the
architecture of close-in extrasolar planets: beyond reproducing the
appropriate mass-scale for the known super-Earth/sub-Neptune population, it
further illuminates that the pebble accretion model of planet formation
predicts a distinct pattern of uniformity among the generated planets. On the
other hand, inefficient settling of solids within the inner region of
protoplanetary disks insinuates that pebble accretion only operates vigorously
beyond the ice-line, thereby preferentially generating water-rich planets – an
expectation that is confirmed by detailed numerical simulations (_Izidoro et
al._ , 2017). Therefore, recent determination of the broad prevalence of
silicate-rich short-period planets (_Rogers and Owen_ , 2021; _Zeng et al._ ,
2019) fosters substantial compositional tension between the pebble accretion
model and observations, casting doubt on pebble accretion as the primary
process responsible for the genesis of sub-Neptune-sized extrasolar planets.

Figure 11: Characteristic mass scales of planet formation. Planetesimals form
within the protoplanetary disk through direct gravitational collapse. The
planetesimal mass estimates are shown in the figure with a dotted lines.
Growth beyond planetesimals occurs through collisions with other planetesimals
as well as gravito-hydrodynamic capture of pebbles. Planetesimal accretion
subsides when protoplanetary embryos reach the planetesimal isolation mass –
shown on the figure with gray lines for a variety of surface density profiles.
Beyond this scale, planetary growth can proceed via giant impacts. Pebble
accretion terminates at the considerably greater pebble isolation mass, shown
on the figure as a black line. Once planets become sufficiently massive to
raise significant wakes in their natal nebulae, orbital migration sets in.

###  4.3 Orbital Migration

An inescapable consequence of planetary growth within a massive gaseous disk
is the onset of planet-disk interactions. Distinct from aerodynamic drag
(which is negligible for planetary-mass objects) the primary mode of angular
momentum exchange between planets and their surrounding gas is gravitational
(_Goldreich and Tremaine_ , 1979; _Lin and Papaloizou_ , 1979; _Ward_ , 1997).
In particular, planetary perturbation of the gas flow in its local
neighborhood manifests as a spiral density wave – or “wake" – that trails the
planet’s orbital phase exterior to its orbit and leads the planet at interior
orbital radii (_Ogilvie and Lubow_ , 2002). It is the gravitational back-
reaction of this wake upon the planet that drives migration.

Somewhat counter-intuitively, gravitational interactions between the planet
and its wake are effectively repulsive: the exterior, trailing part of the
wave continuously pulls back on the planet, sapping it of orbital energy and
causing it to fall deeper into the stellar potential well. The reverse effect
manifests as a consequence of interaction between the planet and the leading
wake. Consequently, orbital changes themselves arise from a subtle imbalance
of the torques exerted by the inner and outer arms of the spiral density wave.

Because the wake itself is a response to planetary gravity and planetary
migration is a response to the wake, the characteristic rate of orbital
transport scales as the product of the effective disk mass,
$\sim\Sigma\,r^{2}$, and the mass of the planet, $M_{p}$. The relevant
expression further exhibits a non-trivial inverse-square dependence on the
speed of sound within the disk (via the disk aspect ratio
$h/r=c_{\rm{s}}/v_{\rm{kep}}$), which stems in part from the fact that waves
are only excited where the gaseous flow past the planet becomes supersonic, as
well as the fact that the torque difference in leading and trailing arms
itself depends on the disk scale height (_Tanaka et al._ , 2002).
Cumulatively, a general expression for the rate of migration has the form:

$\displaystyle\frac{\dot{r}}{r}=\Upsilon\,\bigg{(}\frac{M_{p}}{M_{\star}}\bigg{)}\bigg{(}\frac{\Sigma\,r^{2}}{M_{\star}}\bigg{)}\,\bigg{(}\frac{r}{h}\bigg{)}^{2}\,\Omega,$
(23)

where $\Upsilon$ is a dimensionless coefficient that depends on specific disk
properties.

Substituting the characteristic values for the pebble isolation mass (20) into
this expression, for nominal MMSN-like disk parameters, we obtain
characteristic orbital migration rates that are undeniably fast. For example,
the characteristic migration timescale of a $20\,M_{\oplus}$ planet at $5\,$AU
is $\Upsilon\,r/\dot{r}\lesssim 10^{5}\,$years. Migration rates only
accelerate for shorter-period orbits. In fact, even a Mars-mass embryo in the
inner regions of the disk will be characterized by a nominal migration rate on
the order of $\sim 10^{4}\,$years. The same flavor of planet-disk interactions
additionally leads to decay of planetary eccentricities and inclinations on a
timescale that is even more rapid (by a factor of $\sim(r/h)^{2}$; _Tanaka and
Ward_ 2004; _Kley and Nelson_ 2012). Consequently, orbital migration is almost
certainly an important and active process within protoplanetary disks that
plays a crucial role in regulating the terminal architectures of planetary
systems.

Despite the pronounced importance of planetary migration, its direct
consequences are difficult to calculate in a robustly predictive manner. Being
dictated by a difference in two large numbers (i.e., inner and outer spiral
arm torques), both the direction and magnitude of disk-forced orbital
evolution exhibit sensitive dependence on poorly constrained disk parameters –
a complexity that is captured in the highly non-trivial form of the
coefficient $\Upsilon$. To this end, _Paardekooper et al._ (2010) have
obtained the following form for this coefficient:

$\displaystyle\Upsilon=\gamma\,(p/10-17\beta/10-5/2),$ (24)

where $p$ is the surface density power-law index as before, $\beta$ is the
power-law index of the temperature profile, and $\gamma$ is the ratio of
specific heats (see also _Paardekooper et al._ 2011 for a discussion on the
corotation torque and its saturation). More generally, hydrodynamic
simulations reveal complex mass-period “migration maps" that are flecked with
intricate regions of both inward and outward migration, with variable rates
(_Bitsch et al._ , 2015; _Coleman and Nelson_ , 2014, 2016). Indeed, the
picture is rendered even more complex by consideration of disk turbulence –
that introduce a random component into the migration process (_Nelson_ , 2005;
_Adams et al._ , 2008; _Rein and Papaloizou_ , 2009) – as well as fluctuations
in the nebular entropy gradient, which are themselves controlled in part by
the chemical structure of the disk. All of these intriguing aspects of
migration theory remain an active area of research.

Returning to the specific consideration of close-in sub-Neptunes, let us
highlight two key results pertinent to orbital migration. First and foremost,
the disk’s inner edge, which is presumably set by the magnetic truncation
radius (equation 12), appears to act as a trap to inward-migrating planets
(_Masset et al._ , 2006; _Liu et al._ , 2017; _Romanova et al._ , 2019). This
is not because the planets are pushed into the cavity where the gas density is
tenuous, but instead because the sharp gradient in the surface density profile
dramatically enhances the corotation torque, which causes the inner edge to
repel the planets outward (_Paardekooper and Johansen_ , 2018). On the other
hand, irregularities in the migration portraits of typical disks tend to
diminish over the disk’s lifetime, leading to an expectation of an overall
tendency towards inward migration in mature disks (_Lyra et al._ , 2010).
Thus, the architectures of close-in exoplanetary systems should be broadly
consistent with a formation scenario that is shaped by sustained inward decay
of planetary orbits that terminates at the disk’s inner edge. However, some
challenges still remain for the migration paradigm, particularly the question
of whether it can explain the largely non-resonant orbital architectures of
the compact multis (§2.3) as well as the radial mass-ordering of the planets
(see e.g., _Ogihara et al._ 2015; _Coleman et al._ 2019). We now revisit this
observation in the context of its implications for understanding sub-Neptune
formation and migration.

#### Paucity of Resonances and Near-Resonance Features.

Convergent orbital migration naturally shepherds planets into mean-motion
resonances (e.g., _Mustill and Wyatt_ 2011; _Batygin_ 2015). Thus, the
question of why resonances are intrinsically rare among the compact multis
(recall §2.3, Figure 5) remains an active area of research. Multiple ideas
have been put forth to explain this puzzle over the last decade and a half.
One suggestion is that, in a sufficiently turbulent protoplanetary disk,
resonant capture can be impeded through stochastic fluctuations of the disk’s
gravitational potential (_Adams et al._ , 2008; _Rein and Papaloizou_ , 2009).
Alternatively, _Goldreich and Schlichting_ (2014) argued that a specific
choice of planet-disk interaction parameters can render resonances metastable
within the nebulae. Finally, _McNally et al._ (2019) have argued that in
sufficiently laminar disks, emergence of vortices can prevent resonant
trapping of planets. While these models can indeed provide a resonance
disruption mechanism under certain conditions (_Deck and Batygin_ , 2015;
_Batygin and Adams_ , 2017), a scenario which has had the most success in
matching the period ratio distribution on a quantitative level is one where
primordial resonances are disrupted through a phase of post-nebular
instabilities (_Izidoro et al._ , 2017, 2019; _Matsumoto and Ogihara_ , 2020).

Beyond reproducing the typical orbital architectures themselves, the resonant-
chain instability scenario is broadly consistent with the observed pattern of
uniformity. To this end, the numerical experiments reported in _Goldberg and
Batygin_ (2022) show that if the majority of planetary systems originate as
resonant chains that relax through a transient period of planet-planet
scattering, the process of collisional consolidation that ensues during the
dynamical instabilities does not degrade the mass uniformity of the individual
systems beyond that exhibited by the data. Numerous ideas have been put
forward to explain how sub-Neptunes can enter a phase of planet-planet
scattering (e.g., _Johansen et al._ 2012; _Spalding and Batygin_ 2016;
_Spalding et al._ 2018; _Pichierri and Morbidelli_ 2020; _Petit et al._ 2020).
However, a specific understanding of how dynamical instabilities are triggered
in the proportion needed to explain the period ratio distribution remains
elusive.

While resonances are uncommon overall, they are not completely absent from the
population. In §2.3, we showed that the period ratio distribution contains
small but significant overabundances of planet pairs wide of resonance (see
Figure 5). There have been a wide range of efforts to understand the origin of
these near-resonance features. Toward this end, many authors (_Papaloizou and
Terquem_ , 2010; _Lithwick and Wu_ , 2012; _Batygin and Morbidelli_ , 2013;
_Baruteau and Papaloizou_ , 2013a; _Delisle and Laskar_ , 2014; _Chatterjee
and Ford_ , 2015) have demonstrated that these features can be accounted for
by slow divergence of initially resonant orbits, facilitated by sustained
eccentricity damping. Pertinently, _Millholland and Laughlin_ (2019) have
shown that obliquity tides provide a natural source of this dissipation (see
also _Choksi and Chiang_ 2020; _Goldberg and Batygin_ 2021). Other proposed
ideas involve various (higher order) planet-disk interactions effects (e.g.
_Baruteau and Papaloizou_ , 2013b; _Migaszewski_ , 2015; _Ramos et al._ ,
2017). Nevertheless, even though the dominant sculptors of these near-
resonance features are still uncertain, any comprehensive theory of sub-
Neptune formation must somehow account for them. We will return to this point
when reviewing population synthesis models in §6.

###  4.4 Atmospheric Escape and Sculpting

Regardless of their orbital architectures, planets that form within H/He
nebulae capture gaseous atmospheres. The planets in the compact multis,
however, do not always retain the atmospheres they accrete during the nebular
phase. Instead, their highly-irradiated orbits render these planets
susceptible to thermally-driven atmospheric mass loss. Planets with
sufficiently low masses and close orbits can lose their initial envelopes
entirely, and thus end up as stripped cores. Atmospheric escape sculpts the
radii of the compact multis and weakens their size uniformity from its
primordial state.

At the population-level, evidence for atmospheric escape can be seen in the
distribution of radii of $1-4\ R_{\oplus}$ planets. The distribution is
bimodal, containing a dearth of planets with radii between $1.5-2R_{\oplus}$
that is known as the “radius valley”. The existence of the radius valley was
predicted by models of atmospheric escape (_Owen and Wu_ , 2013; _Lopez and
Fortney_ , 2013; _Jin et al._ , 2014), which tend to produce planets that
either retain their $\sim 0.5\%-1\%$ H/He-dominated envelopes (resulting in
sub-Neptunes) or lose them entirely (resulting in super-Earths). The radius
valley was later confirmed through observations from CKS (_Fulton et al._
2017; see §2.2).

Two types of thermally-driven atmospheric mass loss are influential for
planets in compact multis: photoevaporation and core-powered mass loss. They
differ primarily in their heat source. With photoevaporation, the energy
source is the EUV and X-ray flux from the host star, which photoionizes
hydrogen and heats the upper atmosphere to high temperatures ($T\sim 10^{4}$
K), producing a hydrodynamic outflow that is strongest in the first $\sim 100$
Myr (_Lammer et al._ , 2003; _Baraffe et al._ , 2004; _Murray-Clay et al._ ,
2009; _Owen and Wu_ , 2013; _Lopez and Fortney_ , 2013). With core-powered
mass loss, the energy source is a combination of the stellar bolometric flux
and remnant thermal energy from formation, which is slowly released from the
core to the atmosphere, producing a cooler hydrodynamic outflow that is
sustained over billions of years (_Ginzburg et al._ , 2016, 2018; _Gupta and
Schlichting_ , 2019).

Both the photoevaporation and core-powered mass loss models reproduce
demographic features of the close-in super-Earths/sub-Neptunes as a whole and
the compact multis in particular. They predict similar slopes of the radius
valley as a function of period, incident flux, and stellar mass, and these
agree with the data (_Owen and Wu_ , 2017; _Jin and Mordasini_ , 2018; _Gupta
and Schlichting_ , 2019, 2020). They also predict a breakdown of perfectly-
uniform peas-in-a-pod patterns, given that the least massive and most
irradiated planets in these systems will lose their atmospheres. The data show
a degree of size diversity that is consistent with atmospheric escape; systems
exhibit a tendency towards size ordering (_Ciardi et al._ , 2013; _Millholland
et al._ , 2017; _Kipping_ , 2018; _Weiss et al._ , 2018a; _Millholland and
Winn_ , 2021), with systems much more frequently hosting super-Earths interior
to sub-Neptunes than vice versa. For completeness, we note that mass loss can
be enhanced by large radioactive abundances in close-in planets. The magnetic
truncation radius given by equation (12) is coincident with the semimajor axes
of planets in compact multis. In the region where the disk is magnetically
truncated, the fields continually wrap up, short out, and reconnect. This
activity leads to the acceleration of particle radiation (cosmic rays) that
can drive spallation processes in the reconnection region (_Lee et al._ ,
1998; _Shu et al._ , 1997). A number of short-lived radionuclides (SLRs) can
be generated via spallation, and the predicted cross sections are large enough
that the compact multis could be enhanced in SLR abundances (_Adams_ , 2021).
Energy released from SLR decay can, in turn, amplify mass loss mechanisms from
young and forming planets, as well as remove volatile components (e.g.,
water).

Finally, it is important to note that thermally-driven mass loss is not the
only driver of intra-system size diversity via atmospheric loss. Late-stage
giant impacts, particularly those occurring after disk dispersal, can
significantly reduce the H/He envelopes of sub-Neptunes and sometimes strip
the cores entirely (e.g. _Inamdar and Schlichting_ , 2015, 2016; _Liu et al._
, 2015; _Biersteker and Schlichting_ , 2019; _Kegerreis et al._ , 2020).
Despite their efficiency, the frequency of atmosphere-stripping giant impacts
is still poorly-constrained. The observed degree of intra-system uniformity of
planet sizes and masses could feasibly provide a constraint on their
prevalence; quantitative efforts towards this goal remain to be seen.

##  5 Planet-Planet Interactions

The planet formation process, as detailed in the previous sections, is
exceedingly efficient at producing close-in planets with a characteristic mass
set by a series of planet-disk interactions. So far, however, we have not
comprehensively addressed the interactions between planets that occur during
their growth and after they are fully formed (recall the schematic outline in
Figure 10). These planet-planet interactions dominate at the later stages of
formation, and they can also play an important role in establishing intra-
system uniformity of planet masses and orbital spacings. This section
discusses the influences of dynamical interactions, in the form of
instabilities (§5.1) and pairwise energy exchange via dissipation (§5.2). We
also review how these dynamics may sculpt the planetary mass function (§5.3).

###  5.1 Dynamical Stability Constraints

Once planets have fully formed, with given values for their masses and well-
defined orbital elements, the resulting planetary systems must be dynamically
stable. In general, planetary orbits cannot be spaced too closely without
rendering the system unstable, where the minimum spacing depends on the planet
masses (and orbital eccentricities). This constraint can be written in terms
of the ratio of semi-major axes of adjacent planets,

${a_{2}\over a_{1}}\equiv\Lambda>1+\Delta{R_{H}\over a}\,,$ (25)

where $R_{H}/a$ is determined by the mutual Hill radius (equation 4). For two
planet systems, the minimum value of $\Delta$ can be calculated
$\Delta=2\sqrt{3}$ (_Gladman_ , 1993). For systems of three or more planets,
the minimum value of $\Delta$ required for stability is larger, where
numerical simulations typically find $\Delta=10$ (e.g., _Pu and Wu_ 2015; see
also _Petit et al._ 2020; _Pichierri and Morbidelli_ 2020; _Tamayo et al._
2021 for detailed theoretical inquiries). Note that one way to characterize
planetary systems is by their distributions of spacing parameters $\Lambda$.
The key features of compact multis are that [1] the spacing parameters are
large enough to imply dynamical stability, [2] the distribution of spacing
parameters for the entire sample shows a broad peak at values just larger than
those required for stability, and [3] the variation of spacing parameters for
pairs of planets within a given system is smaller (tighter) than for the
distribution as a whole (see §2.3).

###  5.2 Pairwise Energy Optimization

Many of the observed properties of the peas-in-a-pod architectures can be
understood by considering the tidal equilibrium states for a pair of planets
(_Adams_ , 2019; _Adams et al._ , 2020). Suppose that two planets are forming
from an annulus in a circumstellar disk. Under a range of circumstances, we
expect that the total mass $M_{T}$ available to make planets will be fixed,
and that the total angular momentum $L$ will be conserved. Let us further
constrain the spacing parameter $\Lambda=a_{2}/a_{1}$ of the planetary orbits
to also be determined (but different values of $\Lambda$ can be considered
after the fact). The properties of the system can be specified by the planet
masses $M_{p,1}$ and $M_{p,2}$, the semi-major axes of the orbits $a_{1}$ and
$a_{2}$, the orbital eccentricities $e_{1}$ and $e_{2}$, and the mutual
inclination $i$ between the orbits. The tidal equilibrium state corresponds to
the lowest energy state available to the system subject to the constraints. In
other words, we can find the values of the system parameters
$(M_{p,1},M_{p,2},a_{1},a_{2},e_{1},e_{2},i)$ that provide the lowest energy
state for given fixed values of total mass, angular momentum, and spacing
$(M_{T},L,\Lambda)$.

The solution to this optimization problem shows that the minimum energy state
corresponds to co-planar and circular orbits ($i=e_{1}=e_{2}=0)$ with nearly
equal masses. The optimized value of the mass fraction is given by

$f={M_{p,1}\over M_{T}}={\Lambda+\sqrt{\Lambda}-2\over
3(\Lambda-1)}={\sqrt{\Lambda}+2\over 3(\sqrt{\Lambda}+1)}\,.$ (26)

The corresponding ratio of masses $\eta=M_{p,2}/M_{p,1}$ can be written in the
form

$\eta={2\sqrt{\Lambda}+1\over\sqrt{\Lambda}+2}\,.$ (27)

In the limit of close separations $\Lambda\to 1$, the mass fraction $f\to 1/2$
and the mass ratio $\eta\to 1$. In other words, for close orbital spacing, the
minimum energy state accessible for a pair of forming planets has exactly
equal masses. The mass fraction $f$ decreases slowly with increasing
$\Lambda$, but only falls to $f=1/3$ ($\eta\to 2$) in the limit
$\Lambda\to\infty$. For any reasonable finite value of the spacing parameter,
energy optimization predicts nearly equal mass planets, as suggested by
observations (which typically measure nearly equal radius planets).

The above considerations apply directly to pairs of forming planets. For the
case of three or more planets, energy optimization can be carried out in two
conceptually different ways. Consider a three planet system: In the first
case, one can find the global energy minimum for all three planets. In the
second case, one can find the minimum energy state for both pairs of planets
and then match boundary conditions (such that the middle planet is the outer
planet of the inner pair, as well as the inner planet of the outer pair). The
second approach leads to systems like those observed. Specifically, in their
optimal state, the planets all have zero eccentricity and inclination, and
have nearly equal masses. The planet mass grows slowly with increasing
semimajor axis (by a factor of $\eta$ given by equation (27)). Note that this
slowly increasing progression of planetary masses corresponds to an underlying
surface density distribution of $\Sigma\sim r^{-11/6}$ (see §3.1).

The energy optimization approach has another feature that also conforms to
properties of observed planetary systems. The discussion thus far has not
included the self-gravity of the planets in the energy budget. If we add this
complication to the analysis (_Adams et al._ , 2020), the tendency for
planetary pairs to have nearly equal masses continues to hold when the total
mass (of the pair) is less than a threshold value given by

$M_{C}=M_{\star}\left({R_{p}\over
a_{2}}\right){(\sqrt{\Lambda}+2)(\Lambda-1)^{2}\over
4\alpha_{g}\sqrt{\Lambda}}\,.$ (28)

Here, $\alpha_{g}$ is a dimensionless constant of order unity that depends on
the internal structure of the planet. The mass scale $M_{C}$ is a bifurcation
parameter and has the value $M_{C}\sim 40M_{\oplus}$ for the compact multis of
interest. For low mass planetary pairs, with total mass $M_{T}<M_{C}$, energy
optimization leads to nearly equal mass bodies. For the high mass case,
$M_{T}>M_{C}$, it becomes energetically favorable for one planet to experience
runaway growth and consume the majority of the mass. This trend is also seen
in the observed planetary sample: Systems that contain large planets (e.g,
those of Jovian mass) do not display the same peas-in-a-pod behavior as
systems with lower mass planets (e.g., _Wang_ 2017).

Finally, we note that the bifurcation mass scale (28) depends on the location
of the planets, and varies inversely with the semimajor axis. As a result,
planetary pairs forming in the outer regions of their systems will have much
smaller values of $M_{C}$, which implies that most planetary pairs will exceed
the critical mass scale. This trend implies a clean prediction: Planets with
larger orbits ($a\gg 1$ AU) are much less likely to exhibit peas-in-a-pod
properties.

We note that there are forms of dissipation that can decrease energy while
ensuring angular momentum conservation. This type of mechanism provides a path
to the optimal energy state while preserving angular momentum. One can also
show that if the system is not at a minimum of the energy-momentum, then such
dissipation destabilizes any other type of critical point. This property is
important for understanding stability in this context for the following
reason: finding a Lyapunov function for an equilibrium implies stability.
However, the failure of a Lyapunov function to be definite does not preclude
stability – another such function may exist, for example, or one may be able
to prove stability in some other fashion. The dissipation result implies that
for practical purposes failure of the energy momentum analysis to prove
stability implies instability. Angular preserving dissipation often has a
classic form related to the Hamiltonian dynamics of the system. A canonical
and instructive case is that for rigid body dynamics where the dissipative
angular momentum preserving dynamics take the form:

$I\dot{\omega}=I\omega\times\omega+\alpha
I\omega\times(I\omega\times\omega)\,,$ (29)

where $I$ is the moment of inertia matrix and $\omega$ is the body angular
velocity vector and $\alpha$ is a positive constant which determines the
dissipation rate. Similar forms can be written down for multi-body systems and
fluids; it is also possible to write down effective angular momentum
preserving integrators that generalize symplectic integrators. For a history
of the energy momentum method and such forms of dissipation in various
contexts, including astrophysical settings, see _Bloch et al._ (1996).

###  5.3 The Planetary Mass Function

The distribution of masses of the fundamental objects of interest represents
an important issue for any astronomical study. For the multi-planet systems
considered here, their mass distribution is markedly different from the mass
distribution of the entire exoplanet sample. It is thus useful to consider how
the distributions of planetary masses for the compact multis are related to
those of the entire sample of exoplanets.777Note that this discussion is by
necessity preliminary: The planetary mass function for the entire sample is
still being determined, and many observational selection effects and biases
remain. In addition, for the compact multis, we generally measure the
planetary radius and infer the planetary mass, and this limitation leads to
additional uncertainties.

Mass distributions often display power-law forms (e.g., the stellar initial
mass function; _Salpeter_ 1955). Current data (e.g., _Cumming et al._ 2008;
_Rosenthal et al._ 2021) indicate that the planetary mass function (PMF) has
the approximate form

${df\over dM_{p}}\propto M_{p}^{-1.3}$ (30)

over the mass range $0.1M_{J}\leq{M_{p}}\leq 10M_{J}$, or, equivalently
$30M_{\oplus}\leq{M_{p}}\leq 3000M_{\oplus}$. In contrast, the finding
emerging from observations of compact multis is that, to leading order, the
planets have nearly the same mass $M_{p}\sim 10M_{\oplus}$ (_Millholland et
al._ , 2017). More specifically, the inferred masses of members of these
multi-planet systems show a peaked distribution near this mass scale, which
lies in the regime of ‘super-Earths’. Note that the lower end of applicability
for the power-law mass distribution (30 – 40 $M_{\oplus}$) is roughly
coincident with the bifurcation mass scale of equation (28).

The typical masses of the multi-planet members thus fall just below the masses
found in the global sample of larger planets that display a power-law PMF.
According to the current theory of planet formation (see the following section
and references therein), most planets in this mass range consist of a rocky
core surrounded by an envelope of gas. The mass scale of the rocky core is of
order 10 $M_{\oplus}$, comparable to masses inferred for compact multis and
comparable to the cores of Jupiter and Saturn.888For completeness, we note
that Jupiter might not have a solid core per se. Nonetheless, the chemical
composition of Jupiter is consistent with it having a rocky core of order
$\sim 10M_{\oplus}$ surrounded by a gaseous envelope (which can also have an
enhanced metallicity).

We can thus organize the current observational picture as follows:
Circumstellar disks readily form rocky bodies with super-Earth masses, of
order 10 $M_{\oplus}$. In many systems, including those that ultimately
produce peas-in-a-pod planets, these rocky bodies acquire relatively little
additional mass. In particular, they accrete little gas. In other systems, the