# Debiasing the Minimum-Mass Extrasolar Nebula: On the Diversity of Solid Disk
Profiles

Matthias Y. He Department of Physics & Astronomy, 225 Nieuwland Science Hall,
The University of Notre Dame, Notre Dame, IN 46556, USA  Department of
Astronomy & Astrophysics, 525 Davey Laboratory, The Pennsylvania State
University, University Park, PA 16802, USA  Center for Exoplanets & Habitable
Worlds, 525 Davey Laboratory, The Pennsylvania State University, University
Park, PA 16802, USA  Center for Astrostatistics, 525 Davey Laboratory, The
Pennsylvania State University, University Park, PA 16802, USA  Institute for
Computational & Data Sciences, 525 Davey Laboratory, The Pennsylvania State
University, University Park, PA 16802, USA  Eric B. Ford Department of
Astronomy & Astrophysics, 525 Davey Laboratory, The Pennsylvania State
University, University Park, PA 16802, USA  Center for Exoplanets & Habitable
Worlds, 525 Davey Laboratory, The Pennsylvania State University, University
Park, PA 16802, USA  Center for Astrostatistics, 525 Davey Laboratory, The
Pennsylvania State University, University Park, PA 16802, USA  Institute for
Computational & Data Sciences, 525 Davey Laboratory, The Pennsylvania State
University, University Park, PA 16802, USA  Institute for Advanced Study, 1
Einstein Drive, Princeton, NJ 08540, USA  Center for Computational
Astrophysics, Flatiron Institute, New York, NY 10010, USA

###### Abstract

A foundational idea in the theory of in situ planet formation is the “minimum
mass extrasolar nebula” (MMEN), a surface density profile ($\Sigma$) of disk
solids that is necessary to form the planets in their present locations. While
most previous studies have fit a single power-law to all exoplanets in an
observed ensemble, it is unclear whether most exoplanetary systems form from a
universal disk template. We use an advanced statistical model for the
underlying architectures of multi-planet systems to reconstruct the MMEN. The
simulated physical and Kepler-observed catalogs allow us to directly assess
the role of detection biases, and in particular the effect of non-transiting
or otherwise undetected planets, in altering the inferred MMEN. We find that
fitting a power-law of the form $\Sigma=\Sigma_{0}^{*}(a/a_{0})^{\beta}$ to
each multi-planet system results in a broad distribution of disk profiles;
$\Sigma_{0}^{*}=336_{-291}^{+727}$ g/cm2 and $\beta=-1.98_{-1.52}^{+1.55}$
encompass the 16th-84th percentiles of the marginal distributions in an
underlying population, where $\Sigma_{0}^{*}$ is the normalization at
$a_{0}=0.3$ AU. Around half of the inner planet-forming disks have minimum
solid masses of $\gtrsim 40M_{\oplus}$ within 1 AU. While transit observations
do not tend to bias the median $\beta$, they can lead to both significantly
over- and under-estimated $\Sigma_{0}^{*}$ and thus broaden the inferred
distribution of disk masses. Nevertheless, detection biases cannot account for
the full variance in the observed disk profiles; there is no universal MMEN if
all planets formed in situ. The great diversity of solid disk profiles
suggests that a substantial fraction ($\gtrsim 23\%$) of planetary systems
experienced a history of migration.

Exoplanet systems (484); Exoplanet formation (492); Exoplanets (498);
Extrasolar rocky planets (511); Planetary system formation (1257);
Protoplanetary disks (1300)

††software: NumPy (Harris et al., 2020), Matplotlib (Hunter, 2007),
ExoplanetsSysSim (Ford et al., 2022), SysSimData (Ford, 2019),
SysSimExClusters (He & Ford, 2022), SysSimPyPlots (He, 2022b), SysSimPyMMEN
(He, 2022a)

##  1 Introduction

Well before the discovery of extrasolar planets, the idea of a “minimum mass”
solar nebula (MMSN) was introduced to posit the in situ formation of the
planets in our solar system (Weidenschilling, 1977; Hayashi et al., 1985). The
MMSN framework is intuitively simple: in this view, the initial proto-
planetary disk must have enough material to form the final planets in their
present locations, and in particular must contain enough solids for the
planets to locally accrete up to their core masses within their feeding zones.
This implies that one can work backwards to infer the requisite (solid)
surface density, $\Sigma$, local to each planet by spreading its mass in an
annulus centered around its semi-major axis ($a$). By assuming a smooth disk
profile, a power-law for $\Sigma$ as a function of $a$ is typically fit to
then infer the radial distribution of disk solids.

Since the discovery of thousands of exoplanets, largely due to the
transformational success of NASA’s Kepler mission, numerous studies have
applied the MMSN template to these extrasolar worlds in order to form an
analogous minimum mass extrasolar nebula (MMEN; Chiang & Laughlin 2013;
Raymond & Cossou 2014; Schlichting 2014; Dai et al. 2020). Yet, it continues
to be debated whether most exoplanetary systems conform to a universal disk
profile. Most previous studies have fit a single disk profile, in the form of
a single power-law for $\Sigma(a)$, to all the exoplanet candidates observed
by Kepler simultaneously (Chiang & Laughlin, 2013; Schlichting, 2014; Dai et
al., 2020). This construction has multiple shortcomings: (1) it washes out any
potential system-level correlations, (2) the resulting “MMEN” does not
represent the properties of any single true/physical disk, and (3) it does not
account for planet multiplicity, and specifically the effect of missing
(undetected) planets in each system. One notable exception111A pioneering
study by Kuchner (2004) also fit power-laws to individual planetary systems,
but with exoplanets detected by the radial velocity method and were limited to
a much smaller sample of systems. is Raymond & Cossou (2014) (hereafter RC14),
who fit a power-law to each individual system with three or more observed
planets. They showed that this produces a broad diversity of minimum-mass
disks with profiles ranging from $\Sigma\propto a^{-3.2}$ to $a^{0.5}$, thus
retaining the variances across individual systems and illustrating the
inconsistency of a universal disk profile.

Previous studies have also relied on simplistic treatments for detection
biases and the exoplanet mass-radius relationships to construct the MMEN from
the Kepler planet catalog, typically by applying a correction factor for the
transit geometric and detection probability (i.e., a form of inverse detection
efficiency) of each planet in an attempt to “debias” the observed sample
(Chiang & Laughlin, 2013; Dai et al., 2020). While this approach effectively
weights the planets such that longer period and smaller sized planets are
compensated for their reduced detectability by transits, it is a clear
oversimplification of the Kepler detection pipeline (Christiansen et al.,
2020) and does not correct for missing planets in systems with known
planet(s). The recent development of detailed forward models for the Kepler
mission has enabled unprecedented inferences on the intrinsic population of
inner planetary systems from the observed population (e.g., He et al. 2019),
leading to advanced statistical models such as the “maximum AMD model” that
captures the underlying architectures and correlations in multi-planet systems
(He et al., 2020). The SysSim simulated catalogs, comprised of physical and
observed catalog pairs, provide a way to directly address the impact of non-
transiting or otherwise unseen planets in Kepler-like systems.

In this article, we use simulated catalogs (physical and observed) to assess
how detection biases can affect our interpretation of the MMEN and make
comparisons to the Kepler observed catalog. In §2, we describe how we build
the MMEN, beginning with a summary of the SysSim population model developed
previously to reproduce the Kepler planet catalog (§2.1), and a review of how
planets are converted to solid surface densities using various prescriptions
for their feeding zones found in the literature (§2.2). We discuss the
standard method in which a power-law is fit to all the exoplanets in a given
catalog, and show differences in fitting to the simulated observed and
physical catalogs (§2.3). We then adopt a procedure modified from that of RC14
in which a power-law disk profile is fit to each multi-planet system to
directly construct the distribution of MMEN from the physical catalogs (§2.4).
In §3, this procedure is repeated for observed multi-transiting systems
(simulated and Kepler) to assess the effect of missing planets in altering the
inferred MMEN distribution. We discuss the implications for the distribution
of minimum disk masses in §4. Finally, we summarize and discuss our key
results in §5.

##  2 Constructing the Minimum Mass Extrasolar Nebula

###  2.1 Population model

Figure 1: Left: solid surface densities vs. semi-major axes for Kepler planet
candidates (black) and simulated-observed planets (blue) in multi-transiting
systems. The RC14 prescription is used here for all planets. The dashed lines
show the corresponding power-law fits over the observed planets. For
comparison, the solid green line shows an estimate of the minimum mass solar
nebula (MMSN; equation 2 of Chiang & Youdin 2010). For each power-law fit,
$\beta$ represents the slope while $\Sigma_{0}$ (given in g/cm2) represents
the normalization at $a_{0}=0.3$ AU (equation 7). The power-law fit parameters
using the other $\Delta{a}$ prescriptions are provided in Table 1. Right: same
as the left panel, but with simulated planets from a physical catalog, for
various prescriptions of $\Delta{a}$ (as labeled). Each set of colored points
denotes the same sample of simulated planets, with surface densities computed
using a given $\Delta{a}$ prescription. Likewise, each colored dashed line
denotes a power-law fit over all planets (observed or not) in the catalog. We
note that the RC14 prescription (purple) can only be computed for systems with
at least two planets.

The “maximum AMD model” (He et al. 2020; hereafter, the H20 model) was
developed to describe as many features of the Kepler planet catalog as
possible, with a focus on the correlated properties of planets in multi-
transiting systems, using a combination of statistical distributions and
conditions for dynamical stability (Laskar & Petit, 2017; Petit et al., 2017).
It provides a detailed parameterization of the underlying distribution of
planetary systems between orbital periods of $P=3-300$ days and planet sizes
of $R_{p}=0.5-10R_{\oplus}$ around a purified sample of FGK dwarf stars (see
Hsu et al. 2019 and H20 for a description of the stellar selection criteria).
The Kepler planetary systems in this sample were derived from the Kepler
Objects of Interest DR25 data set obtained from the NASA Exoplanet Archive
(NASA Exoplanet Archive, 2020)222Accessed on 2020-10-19 at 22:34. and include
2169 planet candidates of which 964 are in 388 multi-transiting systems. For
this planet sample, the model reproduces an astounding number of properties at
the population level (see H20; He et al. 2021; Millholland et al. 2021),
including but not limited to the following which are most relevant to this
study: (1) the overall number of planets per star and the observed
multiplicity distribution, (2) the period and period ratio distributions, and
(3) the intra-system size similarity patterns, of which the latter two are
commonly referred to as the “peas-in-a-pod” patterns (see Weiss et al. 2022
for a review).

_Mass-radius relation:_ a planetary mass-radius (M-R) relation is necessary to
estimate planet masses from their radii. The H20 model uses a probabilistic
M-R relation that consists of a lognormal distribution centered around the
“Earth-like rocky” model from Zeng et al. (2019) for small planets and a non-
parametric model defined by a series of Bernstein polynomials fit to a sample
of 127 Kepler planets with RV or TTV masses from Ning et al. (2018) for large
planets, with the transition radius from small to large chosen to be
$1.472R_{\oplus}$ such that both the mean prediction and scatter are
continuous across the entire range of radii. This M-R relation is thus
motivated by a combination of both physical (below the transition radius) and
empirical (above the transition radius) models and is more detailed than
previously adopted relations for computing the MMEN from the Kepler planets.

SysSim enables the generation of physical and observed catalog pairs from the
H20 model (the code is available at Ford et al. 2022; He & Ford 2022). A
physical catalog represents (one realization of) the true, underlying
distribution of planetary systems. Each system in this catalog consists of a
known Kepler target star (with stellar parameters from the Gaia-Kepler Stellar
Properties Catalog; Berger et al. 2020) and a set of planets with physical
radii, masses, orbital periods, and orbital elements drawn directly from the
H20 model. An observed catalog then represents one realization of the detected
transiting planets from the physical catalog under a Kepler-like (primary)
mission, simulated using a detailed model for the combined detection and
vetting efficiency of the Kepler pipeline that accounts for the window and
1-$\sigma$ depth functions of each individual target star (Burke & Catanzarite
2017; Christiansen et al. 2020; see Hsu et al. 2019 for details).

Table 1: Power-law fits (equation 7, with $a_{0}=0.3$ AU) for various
prescriptions of $\Delta{a}$ and ensembles of planets.

| CL13 | RC14 | 10Hill | S14
---|---|---|---|---
Fit all planets |  $\Sigma_{0}$ (g/cm2) | $\beta$ |  $\Sigma_{0}$ (g/cm2) | $\beta$ |  $\Sigma_{0}$ (g/cm2) | $\beta$ |  $\Sigma_{0}$ (g/cm2) | $\beta$
Physical catalogs | $65_{-8}^{+9}$ | $-2.08_{-0.03}^{+0.02}$ | $146_{-15}^{+16}$ | $-1.97_{-0.03}^{+0.05}$ | $568_{-45}^{+47}$ | $-2.05_{-0.02}^{+0.01}$ | $147_{-12}^{+13}$ | $-2.55_{-0.02}^{+0.01}$
Observed catalogs | $246_{-6}^{+7}$ | $-1.84_{-0.02}^{+0.04}$ | $369_{-27}^{+42}$ | $-2.03_{-0.05}^{+0.05}$ | $1396_{-33}^{+22}$ | $-1.89_{-0.02}^{+0.02}$ | $414_{-10}^{+12}$ | $-2.35_{-0.02}^{+0.02}$
Kepler catalog | $264_{-5}^{+5}$ | $-1.67_{-0.02}^{+0.02}$ | $413_{-13}^{+13}$ | $-1.90_{-0.03}^{+0.03}$ | $1466_{-17}^{+20}$ | $-1.77_{-0.01}^{+0.01}$ | $445_{-4}^{+4}$ | $-2.18_{-0.01}^{+0.01}$
Fit each system |  $\Sigma_{0}^{*}$ (g/cm2) |  $\beta$ |  $\Sigma_{0}^{*}$ (g/cm2) |  $\beta$ |  $\Sigma_{0}^{*}$ (g/cm2) |  $\beta$ |  $\Sigma_{0}^{*}$ (g/cm2) |  $\beta$
Physical catalog | $160_{-140}^{+313}$ | $-2.02_{-1.51}^{+1.50}$ | $336_{-291}^{+727}$ | $-1.98_{-1.52}^{+1.55}$ | $1038_{-771}^{+1125}$ | $-2.02_{-1.00}^{+1.00}$ | $263_{-196}^{+330}$ | $-2.53_{-0.99}^{+1.00}$
Observed catalog | $349_{-284}^{+597}$ | $-1.99_{-1.28}^{+1.36}$ | $519_{-420}^{+2100}$ | $-1.97_{-1.27}^{+1.54}$ | $1758_{-1214}^{+1766}$ | $-2.00_{-0.82}^{+0.96}$ | $496_{-328}^{+539}$ | $-2.39_{-0.79}^{+0.85}$
Kepler catalog | $388_{-308}^{+914}$ | $-1.84_{-1.25}^{+1.52}$ | $515_{-428}^{+2255}$ | $-1.80_{-1.28}^{+1.59}$ | $1776_{-1048}^{+2790}$ | $-1.82_{-0.72}^{+1.01}$ | $519_{-332}^{+723}$ | $-2.25_{-0.80}^{+0.91}$

**footnotetext: The fitted values of $\Sigma_{0}$ have been scaled up by a
factor $\alpha$ for each system, as described in §2.4.

Note. — The uncertainties for “Fit all planets” denote the 16th-84th
percentiles computed over many iterations of simulated catalogs (each with the
same number of target stars as the Kepler catalog), while the uncertainties
for “Fit each system” represent the 16th-84th percentiles of the distributions
over all the systems in a single catalog. Since there is only one Kepler
catalog, the uncertainties for “Fit all planets” in the Kepler catalog were
computed via re-samplings of the mass-radius relation for all of the planets.

###  2.2 Computing solid surface densities from planets

The minimum solid surface density ($\Sigma$) required to form a planet is
given by spreading the mass of the planet in solids ($M_{p}$) over an annulus
of width $\Delta{a}$ (which conceptually represents the feeding zone of the
planet) centered at a separation of $a$:

$\Sigma=\frac{M_{p}}{2\pi{a}\Delta{a}}.$ (1)

Since the solid mass of a planet is limited to its core mass and most planets
above $\sim 1.6R_{\oplus}$ have gaseous envelopes (Rogers, 2015), we also
limit $M_{p}$ to $10M_{\oplus}$. The choice of setting $10M_{\oplus}$ as the
maximum core mass is motivated by numerous studies on the critical core mass
for runaway gas accretion (Stevenson, 1982; Pollack et al., 1996; Rafikov,
2006; Lee et al., 2014; Piso et al., 2015), which find $M_{\rm core}\simeq
5-20M_{\oplus}$. Yet, it is possible for some planets to have solid masses
greater than 10 Earth masses (potentially including Jupiter and Saturn;
Guillot 1999; Wahl et al. 2017; Liu et al. 2019); thus, we also repeat our
analyses without setting any maximum core mass and find little change in our
results333For example, while the upper tail of our distribution in solid
surface density normalizations (later defined in §2.4) extends to modestly
higher values, the median value increases by only $\sim 2\%$. Thus, our
results are insensitive to our assumption of the exact core mass limit., due
to the infrequency of large/giant planets in the H20 model. For each planet,
we then compute $\Sigma$ using a number of previously adopted prescriptions
for the feeding zone width $\Delta{a}$:

  1. 1.

The simplest prescription is to adopt a width equal to the semi-major axis
(Chiang & Laughlin 2013; hereafter CL13),

$\Delta{a}=a.$ (2)

While convenient, this prescription tends to result in overlapping regions
between planets in the same system and likely overestimates $\Delta{a}$, thus
underestimating $\Sigma$.

  2. 2.

RC14 recommend using the geometric means of the semi-major axes for
neighboring planets as the dividing boundaries for their feeding zones,

$\displaystyle a_{{\rm sep},i}$ $\displaystyle=\sqrt{a_{i}a_{i+1}},\quad
i=1,\dotsc,N-1$ (3) $\displaystyle\Delta{a_{i}}$ $\displaystyle=a_{{\rm
sep},i}-a_{{\rm sep},i-1},$ (4)

where $a_{{\rm sep},i}$ is the boundary between the $i$ and $i+1$ planets,
$\Delta{a_{i}}$ is the width for the $i^{\rm th}$ planet, and $N$ is the
number of (physical or observed) planets in the system. For the inner (outer)
edge of the innermost (outermost) planet, we define it by enforcing the same
ratio in $a$ both interior and exterior to the planet. We note that this
prescription is only applicable to multi-planet systems and may overestimate
(or underestimate) $\Delta{a}$ when not all planets are detected.

  3. 3.

Alternatively, one can justify using a multiple of the planet’s Hill radius:

$\Delta{a}=kR_{\rm Hill}=k{a}\Big{(}\frac{M_{p}}{3M_{\star}}\Big{)}^{1/3},$
(5)

where $k$ is a constant factor and $M_{\star}$ is the mass of the host star.
Dai et al. (2020) chose $k=10$ as motivated by the spacings of Kepler observed
multi-planet systems (Weiss et al., 2018), which we denote hereafter as
10Hill. While this is appropriate for low-mass planets on low-eccentricity
orbits, it is expected to underestimate the width of the feeding zone for
systems with significant eccentricities.

  4. 4.

Finally, Schlichting (2014) (hereafter S14) suggests

$\Delta{a}=2^{3/2}a\sqrt{\frac{aM_{p}}{R_{p}M_{\star}}},$ (6)

motivated by considerations for the role of giant impacts in dictating a
planet’s effective feeding zone width.

In Figure 1 (left panel), we plot the solid surface density versus semi-major
axis for each planet in the Kepler catalog (black points) as well as in a
simulated observed catalog (blue points). To facilitate the rest of the paper,
the RC14 prescription is used for the purposes of this panel. Similarly, we
plot the solid surface densities for a sample of $10^{3}$ simulated planets
drawn from a simulated physical catalog (right panel). Each planet is repeated
as four points (once for each of the above prescriptions). While there is
significant scatter in $\Sigma$ due to both (1) the range of planet masses (up
to two orders of magnitude even after capping the core masses; $M_{p}\sim
0.1-10M_{\oplus}$) and (2) the varying prescriptions for $\Delta{a}$, there is
a clear linear trend in $\log{\Sigma}$ vs. $\log{a}$ indicative of a power-law
relation between $\Sigma$ and $a$ that is qualitatively consistent with
previous studies. In the next subsection, we describe our procedure for
fitting power-laws to both the total ensemble of planets as well as individual
multi-planet systems.

###  2.3 The canonical power-law model

A power-law model for the MMEN is typically fitted to the solid surface
densities computed from the planets as a function of the semi-major axis, of
the form:

$\Sigma(a)=\Sigma_{0}\bigg{(}\frac{a}{a_{0}}\bigg{)}^{\beta},$ (7)

where $\Sigma_{0}\equiv\Sigma(a_{0})$ is the normalization at separation
$a_{0}$ and $\beta$ is the slope. While $a_{0}$ is typically assumed to be 1
AU, we choose $a_{0}=0.3$ AU so that it is closer to the median $a$ of the
simulated and Kepler planets, and thus reduces the covariance of $\Sigma_{0}$
and $\beta$. Most previous studies have fitted equation 7 to an ensemble of
planets (from transit or RV surveys) to construct a single, “universal” MMEN.
To facilitate direct comparisons with these prior studies, we also fit a
power-law to all the planets in each of the catalogs (Kepler observed,
simulated observed, and simulated physical), as denoted by the various dashed
lines in Figure 1. The complete results for all prescriptions of $\Delta{a}$
are presented in Table 1 under “Fit all planets”.

First, we compare the results of the simulated observed catalogs to the
Kepler-observed catalog. While there is a qualitatively good agreement between
the two, the fits to the simulated planets give slightly lower normalizations
at 0.3 AU ($\sim 5-10\%$ smaller values of $\Sigma_{0}$) for all
prescriptions. Similarly, the values of $\beta$ are slightly steeper in the
simulated observed catalogs than in the actual Kepler observations. These
differences are likely due to having fewer simulated systems with positive
size ordering (“monotonicity”; see Figure 12 of H20) compared to the Kepler
data, i.e. there are slightly more larger planets at short periods relative to
the Kepler systems. However, while these differences are statistically
significant over many simulated catalogs, they are smaller than the
differences arising from the various prescriptions.

Next, we compare how the fits to the simulated planets change between the
physical and observed catalogs. This comparison quantifies how the detection
biases of the transit survey alter the inferred (mean) MMEN. We find that the
values of $\Sigma_{0}$ are a factor of $\sim 2.5-4$ lower for the physical
catalogs compared to the observed catalogs; this is readily explained by the
fact that transit detections are biased towards larger planets (which would
tend to be more massive) at all separations. Interestingly, $\beta$ is only
modestly affected, appearing to be slightly dampened by the detection biases;
surprisingly, the opposite may even be true for the RC14 prescription,
although both physical and observed fits are consistent with $\beta=-2$.

Finally, we focus on the physical catalogs and the differences between the
prescriptions for $\Delta{a}$. While the fits to the observed planets
illustrate how detection biases affect the results and serve as a more direct
comparison to previous studies, the fits to the physical planets should be
interpreted as the true MMEN if one knew of all the planets (within the range
probed by the simulations, 0.04 to 0.88 AU). We find that $\beta$ is steepest
for the S14 prescription, comparable between CL13 and 10Hill, and slightly
shallower for RC14. All of these MMEN slopes are steeper than the value of
$\beta=-1.5$ for the MMSN. The normalization ($\Sigma_{0}$) is highest for
10Hill, reflecting that it also tends to be the narrowest definition of
$\Delta{a}$ of the four prescriptions; intuitively, a smaller feeding zone
implies that a greater solid surface density is necessary to collect the same
amount of solid material for forming a planet. In contrast, the approximation
used by CL13 ($\Delta{a}=a$) results in the lowest value of $\Sigma_{0}$.
Broadly, these results indicate that the mean MMEN is more massive than the
MMSN for the innermost regions, although all of these (averaged) extrasolar
disk profiles must cross under the MMSN model at some point due to the steeper
$\beta$ (e.g., $\sim 0.3$ AU for CL13 and $\gtrsim 1$ AU for RC14).

###  2.4 A diversity of disk profiles: fitting power-laws to individual
systems

Figure 2: Top: power-law fits to the solid surface densities of planets in
each simulated multi-planet system from a physical catalog (a sample of
$10^{4}$ systems are shown). Here, each power-law has been scaled up to be
just above the $\Sigma_{i}$ computed for all planets in the same system; the
median scale factor is $\alpha=2.16$ for systems with at least three planets.
The RC14 prescription is used for all planets. Bottom: the corresponding
distribution of parameters of the power-law fits to each simulated system from
a physical catalog. $\Sigma_{0}^{*}$ is the solid surface density at 0.3 AU
(again, after scaling by $\alpha$) and $\beta$ is the power-law index on the
semi-major axis.

The single power-law model defined above averages over the global population
of exoplanets, but it does not describe the properties of any real or single
planet-forming disk. Furthermore, it fails to capture the diversity of
individual disks. A better approach is given by Raymond & Cossou (2014), who
fit the solid surface density profiles of each individual multi-planet system
to show that it is inconsistent to assume a universal disk profile. By
considering Kepler and RV systems with at least three planets smaller than
$5R_{\oplus}$ or $30M_{\oplus}$, RC14 showed that there is an extremely wide
range of slopes ranging from $\beta=-6.3$ to $5.8$.

Here, we adopt a very similar approach to RC14 by also fitting a power-law
(equation 7) to each system, starting with the multi-planet (2+) systems in
our physical catalog. One issue with this approach is that the resulting
power-law for a given system does not guarantee that there is enough solid
disk mass to form every planet in the system, since by design the power-law
fit will be above some points and below the other points in the $\Sigma$ vs.
$a$ space (the latter of which have under-predicted local solid surface
densities).444This is not a problem for the two-planet systems, since a power-
law can always be fit exactly through both points. To address this issue, we
then “scale up” each power-law such that all planets in a given system are at
or below the curve, by multiplying $\Sigma_{0}$ by a scale factor,
$\alpha={\rm max}\\{\Sigma_{i}/\Sigma(a_{i})\\}$, where $\Sigma_{i}$ is the
solid surface density for the $i^{\rm th}$ planet in the system and
$\Sigma(a_{i})$ is the solid surface density of the power-law fit evaluated at
that planet’s semi-major axis. Hereafter, we will use
$\Sigma_{0}^{*}\equiv\alpha\Sigma_{0}$ to denote the scaled-up solid surface
density normalization. This ensures that the resulting power-law profile
contains enough mass to form each planet in the system while being self-
consistent with the assumption that each planet accreted material from within
its feeding zone – a true “minimum mass” extrasolar nebula.

The resulting distribution of power-law fits to each multi-planet system in a
physical catalog is shown in Figure 2 (using the RC14 prescription). In the
top panel, we plot a sample of power-laws in $\Sigma$ vs. $a$. The bottom
panel shows the corresponding distribution of $\Sigma_{0}^{*}$ and $\beta$.
While the median MMEN slope ($\beta=-1.98$) is the same as the fit to the full
ensemble (Figure 1), there is a broad and symmetric distribution with the
16th-84th percentile ranging from $\beta=-3.5$ to $-0.43$. Similarly, there is
a wide distribution of $\Sigma_{0}^{*}$ (16th-84th percentile ranging from 45
to 1060 g/cm2); the median $\Sigma_{0}^{*}$ is also comparable to the fit over
all planets considering a scale factor was applied to each individual system
(median $\alpha=2.16$ for systems with 3+ planets) but not to the full
ensemble. We remind the reader that our values of $\Sigma_{0}^{*}$ are
normalized at 0.3 AU; projecting the power-laws to 1 AU gives $\Sigma(1\rm
AU)=31.5_{-29.2}^{+142.3}$ g/cm2. While the extreme ends of the distributions
are dominated by fits to systems with just two planets, restricting to 3+
systems still produces a substantially broad distribution:
$\beta=-1.96_{-1.17}^{+1.18}$ and $\Sigma_{0}^{*}=427_{-343}^{+642}$ g/cm2. We
find comparable results (a broad diversity of $\beta$ and $\Sigma_{0}^{*}$)
for the other prescriptions of $\Delta{a}$, as listed in Table 1 under “Fit
each system”. Remarkably, for all except the S14 prescription, the median
slope is highly reminiscent of the predictions from the peas-in-a-pod and
pair-wise energy-optimized configurations of planetary systems (Adams, 2019;
Weiss et al., 2022).

Figure 3: Examples of how missing planets can affect the fitting of MMEN to
individual systems. Each panel shows a comparison of a power-law fit to the
observed planets only (dashed blue line) versus to all the planets (solid blue
line) in a given simulated system. The alternating gray shaded regions denote
the feeding zone widths of each individual physical planet, using the RC14
prescription. Top: an illustration of how undetected planets can alter the
inferred disk surface density as demonstrated on a simulated system. Solid and
hollow circles that are vertically offset denote the same (observed) planet,
which has a biased surface density due to a neighboring missed planet that is
causing either an under- or over-estimated feeding zone width. In this case,
both $\beta$ and $\Sigma_{0}^{*}$ are biased due to the missed outer planets.
Bottom: four simulated systems showing different cases for how the fit to the
observed planets can diverge from the fit to the true system. In panel A,
$\beta$ is biased to near zero due to two missing outer planets, producing an
apparently flat disk profile. In panel B, only the two innermost planets are
observed, leading to an extremely positive value of $\beta$ (and
$\Sigma_{0}^{*}$ which increases arbitrarily towards longer separations) that
appears unphysical. In panel C, a large number of undetected planets can also
lead to an observed disk profile that is steeper than reality. Finally, panel
D shows an example where the inferred disk profile is unaffected despite a
missing innermost planet and a middle planet, both of which are among the most
massive in the system.

##  3 The effect of missing planets

Figure 4: Same as the bottom panel of Figure 2, but for a simulated observed
catalog (left panel) and the Kepler observed catalog (right panel). In other
words, the left panel shows the parameters of the power-law fits to each
simulated multi-planet system in an observed catalog (generated from a
physical catalog after a Kepler-like mission), while the right panel shows the
power-law parameters for Kepler multi-transiting systems. As before,
$\Sigma_{0}^{*}$ is the solid surface density at 0.3 AU and $\beta$ is the
power-law slope on the semi-major axis; the prescription from RC14 is used
here for all planets.

Figure 5: The extent to which the inferred MMEN slope and normalization can be
biased due to undetected planets. Top: Distribution of slope ratios,
$\beta_{\rm obs}/\beta_{\rm phys}$, for the observed systems compared to their
true underlying systems. Bottom: Same as the top, but for the normalization
ratios, $\Sigma_{0,\rm obs}^{*}/\Sigma_{0,\rm phys}^{*}$. In both panels, the
vertical dashed line denotes a ratio of one, i.e. unchanged $\beta$ or
$\Sigma_{0}^{*}$. The listed uncertainties denote the 16th-84th percentiles of
each distribution.

In §2.3, we showed how fitting a power-law to all planets in a physical
catalog leads to some notable differences in the inferred MMEN compared to
fitting to all planets in an observed catalog. While such a comparison
illustrates how the overall distribution of planet radii (and thus masses) and
semi-major axes is affected by transit detection biases, it does not
demonstrate how the biases of a Kepler-like survey affect individual systems,
in the form of missing planets due to their non-transiting geometries or
undetectably small sizes. On the other hand, the approach of constructing MMEN
from individual systems as described in §2.4 enables us to directly study this
effect, by comparing the power-law fits of the observed systems to those of
the physical systems. We note that for the RC14 prescription, the effect of
missing planets in a given system is two-fold: (1) the power-law fit (its
slope and/or normalization) may be significantly altered or biased, and (2)
the assumed feeding zone ($\Delta{a}$) of each planet may also be affected. We
illustrate these biases in Figure 3 (top panel) using a simulated system. In
this particular example, the $4.6M_{\oplus}$ observed planet has an under-
estimated surface density due to an over-estimated feeding zone width caused
by the unseen $1.5M_{\oplus}$ planet exterior to it, leading to a disk profile
that is biased steep. The power-law fit to the underlying system must also be
scaled up to ensure that it can form all the planets, including the
$6.1M_{\oplus}$ planet ($\alpha\simeq 3.8$ in this case). While this example
demonstrates a case in which $\beta$ can appear steeper than reality, a wide
range of other cases are also possible, from $\beta$ that are unchanged,
shallower, or even positive (bottom panels of Figure 3).

In this section, we aim to quantify how these biases affect the overall
observed distribution of MMEN disk profiles. We showed the wide diversity of
underlying disk profiles arising from a physical catalog in the previous
section. Here, we repeat the procedure of fitting a power-law to each
planetary system in the corresponding observed catalog (as well as the Kepler-
observed catalog, to test how well our procedure mimics the real Kepler
planetary systems). In Figure 4, we plot the distribution of power-law fit
parameters for each observed multi-transiting system in a simulated catalog
and the Kepler catalog, analogous to the bottom panel of Figure 2. The RC14
prescription is used for the figure; the results of the other prescriptions
are provided in the bottom half of Table 1 (under “Fit each system”).

###  3.1 Comparing fits to the observed versus physical planets in simulated
systems

Several conclusions can be drawn from comparing the distributions of
$\Sigma_{0}^{*}$ and $\beta$ fitted between the observed systems (left panel
of Figure 4) and the physical systems (bottom panel of Figure 2). First, it is
remarkable that the distributions of power-law slopes are virtually unchanged
between the physical and observed systems; the median $\beta$ and width of the
distribution remain similar for each prescription. This is in contrast to the
fits to all planets in a catalog, which tends to lead to slightly shallower
slopes for observed catalogs versus physical catalogs for all except the RC14
prescription. However, we find that the distributions of normalizations are
significantly skewed to larger values for all prescriptions. While
$\Sigma_{0}^{*}$ can be either under- or over-estimated for any given system,
the median $\Sigma_{0}^{*}$ is a factor of $\sim 1.5$ (RC14) to $\sim 2.2$
(CL13) higher for the observed systems than the physical systems. In addition,
the spread of the $\Sigma_{0}^{*}$ distribution is also substantially
increased. This is likely due to the fact that smaller planets tend to be
missed more often than larger ones, the latter of which tend to also be more
massive, thus biasing the disks inferred from only the observed planets
towards higher masses. We note that this is partly counter-acted by the
increased feeding zone widths due to missing planets in the RC14 prescription,
which explains why the median $\Sigma_{0}^{*}$ is least biased using this
prescription.

Although the overall distribution of $\beta$ is unaffected by detection
biases, both $\beta$ and $\Sigma_{0}^{*}$ can be significantly over- or under-
estimated for any given system due to missing planet(s), as already
illustrated in Figure 3. In Figure 5, we plot histograms of $\beta$ ratios
(top panel) and $\Sigma_{0}^{*}$ ratios (bottom panel) for the simulated
observed systems, where the ratio is the value of the fit to only the observed
planets compared to the fit to all planets in their true underlying systems.
The median $\beta$ ratio is close to unity, consistent with the previous
comparison showing how the median $\beta$ is insensitive to detection biases.
More than half of the observed systems have underestimated values of
$\Sigma_{0}^{*}$.

RC14 argued against a universal MMEN by showing that synthetic populations of
planetary systems generated from a single power-law model (i.e., a fixed value
of $\beta$) fail to produce the wide diversity of disk profiles of the Kepler-
observed systems even after simplistic simulations of transit detections.
While their reported distribution of surface density slopes resulted from fits
to the observed systems only, they also showed that a flat distribution of
$\beta\in[-2.5,0]$ or a Gaussian distribution centered around $\beta=-1.25$
(with standard deviation $\sigma=0.8$) for the underlying distribution appears
to roughly match the observed distribution. While we find a somewhat steeper
value for the median surface density slope ($\beta\simeq-2$) compared to RC14,
we have demonstrated that this is very insensitive to detection biases when
accounting for the role of missing planets using the H20 model, and the
diversity of slopes across systems remains. Additionally, our results show
that while detection biases do broaden (and shift) the observed distribution
of $\Sigma_{0}^{*}$, they cannot account for the full extent of the variance
across multi-planet systems. Thus, we further strengthen the conclusion of
RC14 that there is no universal MMEN, but rather a diversity of MMEN profiles
in the underlying population (Figure 2).

###  3.2 Comparing fits to the simulated versus Kepler systems

As seen in Figure 4, the distributions of MMEN power-laws are very similar
between the SysSim and Kepler observed systems. There is a systematic shift to
slightly shallower median values of $\beta$ for the Kepler systems regardless
of the prescription used (e.g. $\sim-2$ vs. $-1.8$ using RC14), similar to the
results of fitting a power-law to all planets as discussed in §2.3. The
distributions of $\Sigma_{0}^{*}$ are very similar, although the Kepler
catalog appears to have relatively more systems towards the high
$\Sigma_{0}^{*}$ tail as evidenced by the higher 84% percentile values for
some prescriptions. Nevertheless, the likeness of the power-law profiles
between the simulated systems and the Kepler systems illustrates the
robustness of the H20 model, and enables us to estimate how missing planets in
Kepler systems likely affect our inferences on the MMEN. As we have shown
here, the overall distribution of disk surface density slopes is not strongly
affected by detection biases and likely represents the true underlying
distribution. The true disk masses are often higher than what we would
conclude from the observed planets alone. In addition, the diversity of disk
masses is also reduced for the true systems, although there is still a wide
range (as summarized in §2.4). In the next section, we use the results from
the physical systems to infer the primordial distribution of total disk masses
within 1 AU.

##  4 Implications for minimum disk masses

A power-law model for the surface density profile can be integrated to give
the total mass in solids enclosed within a radius $r$ from a star:

$\displaystyle M_{r}$
$\displaystyle=\int_{0}^{2\pi}d\theta\int_{r_{0}}^{r}\Sigma(a)a\,{da}$ (8)
$\displaystyle=\begin{cases}\frac{2\pi\Sigma_{0}^{*}}{(2+\beta){a_{0}}^{\beta}}\big{(}r^{2+\beta}-{r_{0}}^{2+\beta}\big{)},&\beta\neq-2\\\
2\pi\Sigma_{0}^{*}{a_{0}}^{2}\ln(r/r_{0}),&\beta=-2,\end{cases}$ (9)

where $r_{0}$ denotes the inner edge of the disk where the density must
truncate. For $\beta\leq-2$, a non-zero value of $r_{0}$ is also necessary to
avoid an infinite amount of mass. We choose $r_{0}\simeq 0.04$ AU,
corresponding to a 3-day orbital period around a solar mass star (which is
also the minimum period of planets in the H20 model).

We compute $M_{r}$ for each fitted system in a physical catalog, for several
values of $r$, and plot (one minus) the cumulative distributions in Figure 6.
The y-axis should be interpreted as the fraction of planet-forming MMEN disks
with at least $M_{r}$ of solid mass enclosed within a given radius $r$. For
example, we find that most such disks have at least an Earth mass of solids
within even 0.1 AU (the dotted line). Over 40% (10%) of disks have over
$10M_{\oplus}$ ($100M_{\oplus}$) of solids within the same distance. Around a
third of disks have over $40M_{\oplus}$ within 0.5 AU (the dashed line), and
this fraction rises to a half within 1 AU (the solid line). The latter is very
similar to the median disk reported by RC14, which also contains $\sim
40M_{\oplus}$ within 1 AU (evaluating equation 9 using their median fit values
of $\Sigma(1\,{\rm AU})=116$ g/cm2 and $\beta=-1.45$). The observation that
the four curves ($r=0.1$ to 1 AU) closely approach each other past
$M_{r}\gtrsim 200M_{\oplus}$ suggests that for the most massive disks, most of
the material would be concentrated in the very inner regions if all the
planets formed in situ (i.e., these disks have high values of $\Sigma_{0}^{*}$
and steep negative values of $\beta$).

Figure 6: Cumulative distributions of the minimum total mass in solids
enclosed ($M_{r}$) at various separations from a star ($r$, as labeled). The
y-axis denotes the fraction of planet-forming disks with at least $M$ of solid
mass within radius $r$. For example, $\sim 70\%$ of disks have at least
$1M_{\oplus}$ of solids within 0.1 AU (dotted curve). About half of all disks
have more than $40M_{\oplus}$ of solids within 1 AU (solid curve).

##  5 Summary and Discussion

The MMEN remains an insightful framework for understanding the primordial
conditions of protoplanetary disks required to form the planets we see today
in their present locations. Yet despite numerous studies on inferring the
MMEN, two key limitations persist: (1) nearly all previous works rely on
rudimentary treatments of the transit detection biases, and (2) most studies
(with the exception of RC14) also attempt to fit all known exoplanets with a
single power-law model, which does not capture the diversity of planetary
systems and the disks from which they formed. The statistical models for
multi-planet system architectures developed via forward modeling of the Kepler
mission (e.g., He et al. 2019, 2020) provide a detailed and unprecedented way
of inferring the MMEN while overcoming both of the above issues.

In this article, we produce various constructions of MMEN using both simulated
physical and observed catalogs of exoplanetary systems in addition to the true
Kepler catalog, and using a variety of prescriptions for the feeding zone
width of each planet (§2.2). First, we follow the widely prevalent approach of
fitting a single power-law to each catalog to demonstrate the effect of
observational biases on the inferred “mean” MMEN and the differences between
various prescriptions (§2.3). We then fit a power-law to the minimum-mass
surface densities of the planets in each planetary system, producing a broad
distribution of MMEN (§2.4). By repeating this procedure for individual
physical and observed planetary systems, we show that undetected planets in
observed systems can significantly alter the inferred disk profile (§3).
Interestingly, while $\beta$ can be strongly affected for any given system,
the overall distribution is largely unaffected by missing planets. However,
the distribution of $\Sigma_{0}^{*}$ is biased and broader for the observed
systems compared to the physical systems. Altogether, these results
demonstrate that although detection biases do affect the inferred distribution
of solid disk profiles, they do not explain all of the variance in the
observed profiles. There is no universal MMEN if all planets formed in their
present locations.

While our approach is similar to that of RC14, one key difference is that we
scale up each power-law fit (i.e., $\Sigma_{0}^{*}$ $=\alpha\Sigma_{0}$ where
$\alpha\geq 1$) such that the surface density is no less than that local to
any individual planet in the system. This is necessary to ensure that each
resulting disk has enough solid mass to form each and every planet while being
self consistent with the assumption that each planet accreted solid material
from only within its local feeding zone. Our finding that the median
$\alpha\simeq 2$ (for systems with 3+ planets) suggests that this
consideration alone causes previous studies to typically underestimate the
minimum surface densities (and thus minimum disk masses) of planet-forming
disks by a factor of two.

RC14 used this wide diversity of MMEN to argue against the in situ formation
scenario. Viscous disk models generally predict surface density slopes of
$\beta\simeq 0$ to $-2$, depending on the temperature profile (Shakura &
Sunyaev 1973; Chiang & Goldreich 1997; Hure 1998; see RC14 for a concise
review). Observations of cold dust in disk structures also consistently yield
$\beta\simeq-0.4$ to $-1.1$ (Andrews et al., 2009, 2010). Thus, the wide swath
of slopes inferred from applying the MMEN framework to exoplanetary systems,
as we have shown, cannot be fully explained by either observations or theory.
Systems with slopes moderately steeper than $\beta\simeq-2$ may be indicative
of having witnessed a significant radial drift of solid materials prior to
accretion or the migration of planetary embryos or fully-formed planets toward
the inner regions of the disks. More sharply falling profiles in the inner
regions could be caused by a truncation of the inner planet-forming disk due
to the presence of an unknown, distant giant planet. Other extreme slope
profiles, including strongly positive values, that exist in the physical
systems can potentially be explained by more violent dynamical histories,
perhaps involving planetary collisions and/or ejections. Based on the number
of simulated physical systems with $\beta<-4$ or $\beta>0$, we estimate that
at least $\sim 23\%$ of planetary systems experienced a history of migration
and/or planet-planet interactions that prevent the final planet masses and
semi-major axes from conveying information about the initial disk profile.
Even among systems with a more typical disk profile, the great diversity of
disk masses necessary to enable in situ planet formation provides strong
evidence that a substantial fraction of these systems experienced a radial
drift of solids or substantial orbital migration. In any case, multi-stage
formation processes involving the migration, mergers, and scatterings of
planetary cores have also been proposed to explain the typical architectures
of compact planetary systems (e.g., Terquem & Papaloizou 2007; Cossou et al.
2014; Zawadzki et al. 2022). It appears that a combination of mechanisms
beyond the simple in situ accretion scenario is necessary to explain all
system outcomes.

Additionally, the broad range of inferred solid disk profiles implies that
there must have been some very massive or extreme disks, which would likely be
unstable. Previous studies have assessed the stability of the gaseous disk by
assuming a gas-to-dust ratio (typically $\Sigma_{\rm gas}/\Sigma\sim 200$;
e.g. Schlichting 2014). While outside the scope of this paper, future studies
may perform a detailed calculation of the stability of the disks inferred from
our H20 model using a range of gas-to-dust ratios.

As an alternative to invoking the large scale migration of planets post-
formation, it is also possible to consider the radial redistribution of solids
(in the form of dust or small pebbles) from which in situ planet formation
then occurs. Rather than starting from a smooth disk profile, the inward drift
of sub-meter sized “pebbles” may collect at pressure maxima creating a series
of gravitationally unstable rings, from which planets form “inside-out”
(Chatterjee & Tan, 2014; Tan et al., 2016). The pile-up of solids into narrow
annuli from fractions to a few AU which serve as the sites of planet formation
can produce steeper radial profiles than the initial disk profile (Drażkowska
et al., 2016), and has recently been applied to the inner solar system via the
silicate sublimation line at $\sim 1$ AU (Morbidelli et al., 2022).

In theory, one may also evaluate the feasibility of the MMEN model by
comparing the distribution of minimum disk masses from the exoplanet
population to the distribution of observed protoplanetary disks. Yet, this is
challenging for a number of reasons, including the biases inherent to either
population and the difficulty of measuring disk masses (see Drazkowska et al.
2022 and Miotello et al. 2022 for a review). Observations using the Atacama
Large Millimeter/sub-millimeter Array (ALMA) have revealed that Class 0
(young, $<0.5$ Myr) disks typically have $\sim 100M_{\oplus}$ of solids while
the older, Class I/II disks have less than $\sim 50M_{\oplus}$ (Tychoniec et
al., 2020; Anderson et al., 2022). However, these values are sensitive to
assumptions for the dust opacity, and extend out to $\sim 100$ AU. On the
other hand, ignoring the effects of photons scattering off dust grains can
lead to underestimated disk masses (Zhu et al., 2019). There may also be
differences in the stellar samples; ALMA mostly observed disks around K and M
dwarfs (e.g., Pascucci et al. 2016; Ansdell et al. 2016; Andrews et al. 2018),
while most exoplanets found by Kepler are around FGK dwarfs. Despite these
differences, it has also been suggested that the solid masses of planets and
of dust may be similar per star (Mulders et al., 2021). It remains that we
have a limited understanding of the total solid masses available in the
innermost regions of protoplanetary disks, and thus the efficiency of planet
formation.

Recent studies have also suggested that the MMEN is dependent on the stellar
mass, and to a smaller extent, stellar metallicity (Dai et al., 2020). We did
not find any correlation between the disk profiles (minimum mass or slope) and
stellar mass. However, it is possible that any underlying correlation is lost
due to the scrambling of fitted disk profiles caused by missing planets as we
have shown in §3. In any case, the inter-system variation in $\Sigma_{0}^{*}$
(which can change by a few orders of magnitude) is significantly greater than
the range of stellar masses (which at most change by a factor of $\sim 2$
across the FGK range). A more detailed analysis of how the MMEN varies with
host star properties is outside the scope of this paper and should be explored
in future work.

The distribution of primordial disk profiles serves as a fundamental input for
the initial conditions of planet formation simulations (e.g., Hansen & Murray
2013; Moriarty & Ballard 2016; Dawson et al. 2016; MacDonald et al. 2020). It
has been shown through simulations of the in situ assembly of planetesimals by
giant impacts that a diversity of solid disk normalizations can lead to an
ensemble of planetary systems resembling the Kepler exoplanetary systems
(MacDonald et al., 2020). However, these studies have typically fixed the
surface density slope (e.g., $\beta=-1.5$ or $-2.5$; Hansen & Murray 2013;
Moriarty & Ballard 2016; MacDonald et al. 2020) due to the large number of
tunable parameters. While it is unlikely that all exoplanetary systems formed
exclusively in their present locations, the in situ model may still explain a
wide variety of planetary system outcomes given an adequately flexible array
of initial conditions. The MMEN framework can continue to provide constraints
on these formation conditions.

The code for fitting MMEN to individual systems, as well as for reproducing
the figures and results of this manuscript, are available via the
“SysSimPyMMEN” package (He, 2022a). The simulated catalogs used in this study
and code for reproducing the H20 model are downloadable from the
“SysSimExClusters” package (He & Ford, 2022).

## Acknowledgements

We thank Darin Ragozzine, Danley Hsu, Robert Morehead, and Keir Ashby for
contributions to the broader SysSim project. We are grateful to Sarah
Millholland, Lauren Weiss, and Chao-Chin Yang for helpful discussions. We also
thank the anonymous referee for their constructive review and comments. M.Y.H.
acknowledges the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC), funding reference number PGSD3 - 516712 - 2018.
M.Y.H. and E.B.F. acknowledge support from the Penn State Eberly College of
Science and Department of Astronomy & Astrophysics, the Center for Exoplanets
and Habitable Worlds, and the Center for Astrostatistics. This research has
made use of the NASA Exoplanet Archive, which is operated by the California
Institute of Technology, under contract with the National Aeronautics and
Space Administration under the Exoplanet Exploration Program.

## References

  * Adams (2019) Adams, F. C. 2019, MNRAS, 488, 1446, doi: 10.1093/mnras/stz1832
  * Anderson et al. (2022) Anderson, A. R., Williams, J. P., van der Marel, N., et al. 2022, arXiv e-prints, arXiv:2204.08731.  https://arxiv.org/abs/2204.08731
  * Andrews et al. (2009) Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502, doi: 10.1088/0004-637X/700/2/1502
  * Andrews et al. (2010) —. 2010, ApJ, 723, 1241, doi: 10.1088/0004-637X/723/2/1241
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46, doi: 10.3847/0004-637X/828/1/46
  * Berger et al. (2020) Berger, T. A., Huber, D., van Saders, J. L., et al. 2020, AJ, 159, 280, doi: 10.3847/1538-3881/159/6/280
  * Burke & Catanzarite (2017) Burke, C. J., & Catanzarite, J. 2017, Planet Detection Metrics: Window and One-Sigma Depth Functions for Data Release 25, Kepler Science Document KSCI-19101-002, id. 14. Edited by Michael R. Haas and Natalie M. Batalha
  * Chatterjee & Tan (2014) Chatterjee, S., & Tan, J. C. 2014, ApJ, 780, 53, doi: 10.1088/0004-637X/780/1/53
  * Chiang & Laughlin (2013) Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444, doi: 10.1093/mnras/stt424
  * Chiang & Youdin (2010) Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493, doi: 10.1146/annurev-earth-040809-152513
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368, doi: 10.1086/304869
  * Christiansen et al. (2020) Christiansen, J. L., Clarke, B. D., Burke, C. J., et al. 2020, AJ, 160, 159, doi: 10.3847/1538-3881/abab0b
  * Cossou et al. (2014) Cossou, C., Raymond, S. N., & Pierens, A. 2014, in Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham, Vol. 299, 360–364, doi: 10.1017/S1743921313008958
  * Dai et al. (2020) Dai, F., Winn, J. N., Schlaufman, K., et al. 2020, AJ, 159, 247, doi: 10.3847/1538-3881/ab88b8
  * Dawson et al. (2016) Dawson, R. I., Lee, E. J., & Chiang, E. 2016, ApJ, 822, 54, doi: 10.3847/0004-637X/822/1/54
  * Drażkowska et al. (2016) Drażkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105, doi: 10.1051/0004-6361/201628983
  * Drazkowska et al. (2022) Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2022, arXiv e-prints, arXiv:2203.09759.  https://arxiv.org/abs/2203.09759
  * Ford (2019) Ford, E. 2019, ExoJulia/SysSimData: Initial release of Data Files for the Exoplanet System Simulator, v1.0.0, Zenodo, doi: 10.5281/zenodo.3255313
  * Ford et al. (2022) Ford, E., He, M. Y., Hsu, D., & TagBot, J. 2022, ExoJulia/ExoplanetsSysSim.jl: He et al. 2020, 2021ab, v1.1.1, Zenodo, doi: 10.5281/zenodo.5915004
  * Guillot (1999) Guillot, T. 1999, Planet. Space Sci., 47, 1183, doi: 10.1016/S0032-0633(99)00043-4
  * Hansen & Murray (2013) Hansen, B. M. S., & Murray, N. 2013, ApJ, 775, 53, doi: 10.1088/0004-637X/775/1/53
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hayashi et al. (1985) Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 1100–1153
  * He (2022a) He, M. Y. 2022a, hematthi/SysSimPyMMEN: Initial release, v1.0.0, Zenodo, doi: 10.5281/zenodo.7117309
  * He (2022b) —. 2022b, hematthi/SysSimPyPlots: v1.1.0: New functions for plotting galleries of systems, v1.1.0, Zenodo, doi: 10.5281/zenodo.7098044
  * He & Ford (2022) He, M. Y., & Ford, E. 2022, ExoJulia/SysSimExClusters: Code for He et al. (2020), v2.1.0, Zenodo, doi: 10.5281/zenodo.5963884
  * He et al. (2019) He, M. Y., Ford, E. B., & Ragozzine, D. 2019, MNRAS, 490, 4575, doi: 10.1093/mnras/stz2869
  * He et al. (2021) —. 2021, AJ, 161, 16, doi: 10.3847/1538-3881/abc68b
  * He et al. (2020) He, M. Y., Ford, E. B., Ragozzine, D., & Carrera, D. 2020, AJ, 160, 276, doi: 10.3847/1538-3881/abba18
  * Hsu et al. (2019) Hsu, D. C., Ford, E. B., Ragozzine, D., & Ashby, K. 2019, AJ, 158, 109, doi: 10.3847/1538-3881/ab31ab
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hure (1998) Hure, J.-M. 1998, A&A, 337, 625
  * Kuchner (2004) Kuchner, M. J. 2004, ApJ, 612, 1147, doi: 10.1086/422577
  * Laskar & Petit (2017) Laskar, J., & Petit, A. C. 2017, A&A, 605, A72, doi: 10.1051/0004-6361/201630022
  * Lee et al. (2014) Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95, doi: 10.1088/0004-637X/797/2/95
  * Liu et al. (2019) Liu, S.-F., Hori, Y., Müller, S., et al. 2019, Nature, 572, 355, doi: 10.1038/s41586-019-1470-2
  * MacDonald et al. (2020) MacDonald, M. G., Dawson, R. I., Morrison, S. J., Lee, E. J., & Khandelwal, A. 2020, ApJ, 891, 20, doi: 10.3847/1538-4357/ab6f04
  * Millholland et al. (2021) Millholland, S. C., He, M. Y., Ford, E. B., et al. 2021, AJ, 162, 166, doi: 10.3847/1538-3881/ac0f7a
  * Miotello et al. (2022) Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. I., & Kataoka, A. 2022, arXiv e-prints, arXiv:2203.09818.  https://arxiv.org/abs/2203.09818
  * Morbidelli et al. (2022) Morbidelli, A., Baillié, K., Batygin, K., et al. 2022, Nature Astronomy, 6, 72, doi: 10.1038/s41550-021-01517-7
  * Moriarty & Ballard (2016) Moriarty, J., & Ballard, S. 2016, ApJ, 832, 34, doi: 10.3847/0004-637X/832/1/34
  * Mulders et al. (2021) Mulders, G. D., Pascucci, I., Ciesla, F. J., & Fernandes, R. B. 2021, ApJ, 920, 66, doi: 10.3847/1538-4357/ac178e
  * NASA Exoplanet Archive (2020) NASA Exoplanet Archive. 2020, Kepler Objects of Interest DR25, Version: 2020-10-16 10:34, NExScI-Caltech/IPAC, doi: 10.26133/NEA5
  * Ning et al. (2018) Ning, B., Wolfgang, A., & Ghosh, S. 2018, ApJ, 869, 5, doi: 10.3847/1538-4357/aaeb31
  * Pascucci et al. (2016) Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125, doi: 10.3847/0004-637X/831/2/125
  * Petit et al. (2017) Petit, A. C., Laskar, J., & Boué, G. 2017, A&A, 607, A35, doi: 10.1051/0004-6361/201731196
  * Piso et al. (2015) Piso, A.-M. A., Youdin, A. N., & Murray-Clay, R. A. 2015, ApJ, 800, 82, doi: 10.1088/0004-637X/800/2/82
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
  * Rafikov (2006) Rafikov, R. R. 2006, ApJ, 648, 666, doi: 10.1086/505695
  * Raymond & Cossou (2014) Raymond, S. N., & Cossou, C. 2014, MNRAS, 440, L11, doi: 10.1093/mnrasl/slu011
  * Rogers (2015) Rogers, L. A. 2015, ApJ, 801, 41, doi: 10.1088/0004-637X/801/1/41
  * Schlichting (2014) Schlichting, H. E. 2014, ApJ, 795, L15, doi: 10.1088/2041-8205/795/1/L15
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Stevenson (1982) Stevenson, D. J. 1982, Planet. Space Sci., 30, 755, doi: 10.1016/0032-0633(82)90108-8
  * Tan et al. (2016) Tan, J. C., Chatterjee, S., Hu, X., Zhu, Z., & Mohanty, S. 2016, IAU Focus Meeting, 29A, 6, doi: 10.1017/S1743921316002313
  * Terquem & Papaloizou (2007) Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110, doi: 10.1086/509497
  * Tychoniec et al. (2020) Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19, doi: 10.1051/0004-6361/202037851
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649, doi: 10.1002/2017GL073160
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, Ap&SS, 51, 153, doi: 10.1007/BF00642464
  * Weiss et al. (2022) Weiss, L. M., Millholland, S. C., Petigura, E. A., et al. 2022, arXiv e-prints, arXiv:2203.10076.  https://arxiv.org/abs/2203.10076
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48, doi: 10.3847/1538-3881/aa9ff6
  * Zawadzki et al. (2022) Zawadzki, B., Carrera, D., & Ford, E. B. 2022, arXiv e-prints, arXiv:2202.05342.  https://arxiv.org/abs/2202.05342
  * Zeng et al. (2019) Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proceedings of the National Academy of Science, 116, 9723, doi: 10.1073/pnas.1812905116
  * Zhu et al. (2019) Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18, doi: 10.3847/2041-8213/ab1f8c