# Randomness and Retention: Using Weak Mean Motion Resonances to Constrain
Neptune’s Late-Stage Migration

Arcelia Hermosillo Ruiz,1 Harriet C.P. Lau,2 Ruth Murray-Clay1
1Department of Astronomy and Astrophysics, University of California Santa Cruz
2Department of Earth, Environmental, and Planetary Sciences, Brown University
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Planet-planetesimal interactions cause a planet to migrate, manifesting as a
random walk in semi-major axis. In models for Neptune’s migration involving a
gravitational upheaval, this planetesimal-driven migration is a side-effect of
the dynamical friction required to damp Neptune’s orbital eccentricity. This
migration is noisy, potentially causing Trans Neptunian Objects (TNOs) in mean
motion resonance to be lost. With Nbody simulations, we validate a previously-
derived analytic model for resonance retention and determine unknown
coefficients. We identify the impact of random-walk (noisy) migration on
resonance retention for resonances up to fourth order lying between 39 au and
75 au. Using a population estimate for the weak 7:3 resonance from the well-
characterized Outer Solar System Origins Survey (OSSOS), we rule out two
cases: (1) a planetesimal disk distributed between 13.3 and 39.9 au with
$\gtrsim$ 30 Earth masses in today’s size distribution and $T_{\rm
mig}\gtrsim$ 40Myr and (2) a top-heavy size distribution with $\gtrsim$2000
Pluto-sized TNOs and $T_{\rm mig}\gtrsim$10Myr, where $T_{\rm mig}$ is
Neptune’s migration timescale. We find that low-eccentricity TNOs in the
heavily populated 5:2 resonance are easily lost due to noisy migration.
Improved observations of the low-eccentricity region of the 5:2 resonance and
of weak mean motion resonances with Rubin Observatory’s Legacy Survey of Space
and Time (LSST) will provide better population estimates, allowing for
comparison with our model’s retention fractions and providing strong evidence
for or against Neptune’s random interactions with planetesimals.

###### keywords:

celestial mechanics—Kuiper belt—planets and satellites: dynamical evolution

††pubyear: 2023††pagerange: Randomness and Retention: Using Weak Mean Motion
Resonances to Constrain Neptune’s Late-Stage Migration–A

##  1 INTRODUCTION

We are approaching an era of unprecedented well-characterized data for the
outer solar system with the Vera C. Rubin Observatory’s Legacy Survey of Space
and Time (LSST). In preparation, we revisit a previously-derived analytic
model for the random walk of a planet, in semi-major axis, during
planetesimal-driven migration. This model has been used to constrain the Solar
System’s planetesimal disk and early planetary migration using the observed
population of trans-Neptunian objects (TNOs) in 3:2 first-order mean motion
resonance (MMR) with Neptune Murray-Clay & Chiang (2006) (hereafter MC2006).
We look ahead to LSST data which will observe tens of thousands of TNOs,
providing enough detections and well-characterized completeness to extend this
analysis across a broader range of resonances. In particular, the survey will
allow us to use the fact that higher-order (weaker) resonances lose objects
more easily to provide stronger constraints on the final, planetesimal-driven
stage of Neptune’s migration.

Over the last 20 years, many theoretical models have been developed to
reproduce the orbital structure and relative populations of TNOs in mean
motion resonances and in the classical and scattered regions of trans-
Neptunian space. While none have succeeded in producing all features of this
observed structure, they each have their virtues and approaches in solving the
puzzle of the early Solar System’s dynamical evolution (see reviews by
Nesvorný 2018; Malhotra 2019; Morbidelli & Nesvorný 2020; Gladman & Volk
2021). Early studies of TNO dynamics considered the outward migration of
Neptune due to angular momentum exchange with a residual disk of planetesimals
(Fernandez & Ip, 1984). This outward migration, later coined “smooth
migration", is efficient at capturing TNOs into mean motion resonances
(Malhotra, 1993, 1995; Hahn & Malhotra, 2005) but overpopulates these regions
in semi major axis space over the classical region between the 3:2 and 2:1
mean motion resonances with Neptune. In particular, TNOs in the 3:2 resonance
are observed to be $\simeq 2-4$ times less numerous than those in the hot
classical belt (classical objects with high eccentricity and inclination)
(Gladman et al., 2012), but Hahn & Malhotra (2005) “smooth migration" model
produces a 3:2 MMR population $\sim$1.5 times more numerous instead. Nesvorný
& Vokrouhlický (2016) finds that "graininess" introduced by placing 10-40% of
the planetesimal disk into Pluto-sized planetesimals overcomes this
difficulty, producing a ratio between the hot classicals and 3:2 resonant
population between 3-4 (see also Kaib & Sheppard (2016); Lawler et al.
(2019)). This graininess produces a random walk component in semi major axis
for Neptune, similar to the model discussed in this paper. Neptune’s random
walk causes resonances to lose previously captured TNOs, adding to the hot
classical population. While they solve this resonance overpopulation problem,
their Nbody simulation does not reproduce the orbital element distributions of
the hot classical population. In their best two migration simulations, the
original disk contains 1000 planetesimals, each with either the mass of Pluto
or twice the mass of Pluto.

Alternatively, evidence of a late heavy bombardment on Earth’s Moon inspired
models suggesting that the early Solar System may have had a giant planet
instability early in its lifetime (Tsiganis et al., 2005; Pike et al., 2017;
de Sousa et al., 2020, and references therein). In these instability models,
Neptune is launched onto an eccentric orbit, requiring a short period of
planetesimal-driven migration at the end to damp its eccentricity to its
current value (Levison et al., 2008; Dawson & Murray-Clay, 2012; Wolff et al.,
2012). Balaji et al. (2023) found that a planetary instability that scatters
planetesimals into semi major axis and eccentricity parameter space around the
3:2 MMR reproduces the semi major axis, eccentricity, and resonant angle
distributions of the observed 3:2 population, but not the libration amplitude
distribution. The mismatch with the libration amplitude distribution is modest
enough that it may result from the addition of transient sticking (Lykawka &
Mukai, 2007; Yu et al., 2018) or from a brief late-stage planetesimal-driven
migration, which was not rigorously modeled in their work.111The cumulative
distributions of the resonant angle did not differ greatly from the libration
amplitude distribution in their work; rather the statistical comparison was
less constraining for the resonant angle because it accounts for its periodic
nature. In addtion, the two distributions are biased by observations
differently, such that the block longitudes relative to Neptune place limits
on the libration amplitude magnitudes and range of phi that are observable
(Volk et al., 2016).  The model developed in this paper applies both to
scenarios in which dynamical interactions with planetesimals are the primary
driver of migration and to scenarios in which a large-scale dynamical upheaval
occurred among the giant planets, followed by a shorter epoch of planetesimal-
driven eccentricity damping and migration.

In this work, we quantify the efficiency of retaining (and thus, losing)
planetesimals in first- to fourth-order MMR with Neptune during an epoch of
noisy planetesimal-driven planetary migration. We evaluate how various
planetesimal size distributions affect migration to constrain the amount of
mass in large objects and typical eccentricity of planetesimals in the early
disk as well as the planets’ migration timescale. Simulating fully-realistic
planetesimal-driven migration and resonance retention requires billions of
massive planetesimals with a continuous mass distribution, since resonance
retention is a function of planetesimal size distribution. Computational
limitations hinder such simulations; for example, Nesvorný & Vokrouhlický
(2016) used $\sim 5000$ CPU days with 2000 particles and a size distribution
represented by a delta function. To overcome this challenge, MC2006
analytically model the planet-planetesimal encounters in analogy to Brownian
motion, where the migrating planet’s semi-major axis fluctuates about its mean
value. As demonstrated in MC2006, the random component of the migration is
dominated by planetesimal sizes with a maximum number density times mass
squared, which is generally at the high-mass end of the planetesimal size
distribution even though planetesimals with the largest sizes do not dominate
the total mass. They conclude that migration was likely smooth enough to
retain TNOs in first-order resonances for planetesimal disks with the bulk of
the mass in bodies having sizes of order $\sim 1$ km and maximum body size of
1000 km.

In this study, we revisit and extend the work of MC2006 with two aims in mind.
First, we wish to assess the validity of planetary migration as a form of
Brownian motion by comparing results from the analytical model of MC2006 with
numerical simulations. It is important that as the planet interacts with
planetesimals, it performs a random walk in semi-major axis. We also wish to
determine, by simulation, the unknown coefficients that arose during MC2006’s
order-of-magnitude derivations of the change in semi-major axes of the planet
and planetesimals perturbing each other during gravitational encounters.

Second, we wish to apply the theoretical expressions established in MC2006
(i.e. total random walk of a planet, retention fraction in MMR) to weaker
resonances and for a wider range of disk conditions and planetesimal size
distributions. While the relative population of different resonances depends
both on capture and retention processes, we can use retention alone as a
constraint because resonances observed to host large populations of TNOs must
have been able to retain these bodies. We investigate the conditions under
which weak resonances can retain TNOs. We find the 5:2 and 7:3 resonances most
promising for placing constraints on Neptune’s total migration timescale and
distance and the planetesimal disk’s size distribution since they are weak,
contain at least some observed objects, are not located in the cold classical
region, and are located near each other in semi-major axis. The 5:2 has a
large observed population, comparable to that in the 3:2 (Gladman et al.,
2012; Volk et al., 2016), while the 7:3 has an observed population that is
present but small (Gladman et al., 2012; Crompvoets et al., 2022). Neptune’s
migration is limited to an amount that would not produce a total loss to the
weakest 7:3 population and allows retention of many objects in the 5:2.
Comparing population estimates for resonances observed by the Outer Solar
System Origins Survey (OSSOS; Bannister et al. 2018) with retention fractions
produced by our model allows us to comment on the impact Neptune’s noisy
migration may have had on the resonances. If Neptune’s noisy migration had a
large impact, the mean motion resonances would preserve the retention fraction
pattern found in our work. We carry out this comparison for an assumed initial
surface denisty profile, assuming two commonly discussed methods for filling
the mean motion resonances.

In Section 2 we introduce the analytical theory describing the random walk
component of planetary migration. We then confirm that N-body integrations
produce a change in the planet’s semi-major axis that matches analytical
theory. In Section 3 we apply the analysis to all resonances (3:2, 11:7, 8:5,
5:3, 7:4, 2:1, 7:3, 5:2, 3:1, 4:1) and calculate the timescales for which
objects in resonance will be lost and their resulting retention fractions. In
Section 4 we use two end-member size distributions of the original disk during
Neptune’s migration to constrain the amount of largest planetsimals in the
disk. In Section 5, we use the relative retention probabilities to interpret
relative resonance populations constrained by OSSOS. In Section 6, we discuss
future work and additional considerations. Finally, in Section 7, we summarize
our results.

##  2 analytical model and numerical verification

Planetesimal-driven migration results from a superposition of discrete
scattering events due to multiple encounters between the planet and
planetesimals. This type of migration can be separated into average and random
components. The average component arises from the effects of dynamical
friction in a uniform sea of planetesimals, modified by asymmetries in mass,
eccentricity, and number of planetesimals interior and exterior to the
planet’s orbit. For a single planet, dynamical friction causes the planet to
lose orbital energy and move toward the star. In the solar system, global
transfer of planetesimals between the giant planets typically leads Jupiter to
migrate inward while Saturn, Uranus, and Neptune migrate outward (Fernandez &
Ip, 1984). In this work, we remain agnostic about the planet’s average
migration. The random component of migration is due to noise that arises from
variations in the orbital elements of planetesimals interacting with the
planet and variations in the number of encounters. Variations in the mass of
planetesimals also affect migration, such that tens of thousands of low-mass
planetesimals will produce less noise than thousands of high-mass
planetesimals. We treat this fact by applying the most massive planetesimal
mass to our analysis in this paper. This noise can cause resonances to lose
previously captured TNOs, which cannot adjust to abrupt changes in the
planet’s orbit that occur on timescales less than the resonant libration
period.

In Section 2.1 we introduce the theory of MC2006, which provides order-of-
magnitude expressions describing the random component of migration, which we
refer as a random walk for the remainder of the paper. We then turn to
numerical simulations to verify selected fundamental equations that predict
the change in the semi-major axis of a planet ($a_{\rm p}$) due to
interactions with single (Section 2.2) and multiple (Section 2.3)
planetesimals. Throughout this section, we choose planet and planetesimal
parameters to facilitate numerical simulation. These choices, which are
unrealistic for the young circumsolar disk, allow us to verify our analytic
expressions, which can then be applied confidently under realistic conditions,
such as those considered in Section 3. For all numerical simulations, we use
the N-body integrator _Mercury_ , version 6.2 (Chambers, 1999) with the
Bulirsch Stoer integrator.

###  2.1 Theoretical Introduction

As a planet and planetesimal approach each other, the planet will perturb the
planetesimal given its larger size and mass. The impulse of the planet on the
planetesimal will change its velocity, and the amount depends on the mass of
the planet, $M_{\rm p}$ how close they are to each other in semi-major axis,
$x$, and the encounter time, $\Delta t$. This change in velocity alters both
the angular momentum and the energy of the planetesimal. MC2006 make use of
the Jacobi integral to disentangle these changes and determine the
planetesimal’s new orbital energy. Given that energy is a function of semi-
major axis, we can then solve for the change in semi-major axis, $\Delta a$,
of the planetesimal. The perturbation on the planetesimal exerts a back
reaction onto the planet and due to conservation of energy of the system, we
can find the change in the planet’s semi-major axis, $\Delta a_{p}$, as a
result of the encounter.

Realistically, a planet doesn’t experience single encounters solely, for it is
embedded in a sea of planetesimals. The planet will undergo a large number of
encounters over a given encounter time, and the cumulative impact of each
individual shift will produce a random step in the planet’s semi-major axis.
Derivations in MC2006 assume that the number of planetesimals encountering the
planet over a given dynamical time (approximately an encounter time) is random
and follows a Poisson distribution, with the average determined by the disk’s
surface density. Planet-planetesimal encounters are random and occur on
timescales shorter than the libration period of a resonant object (i.e.
$\lesssim 10^{2}$ yrs versus $\sim 10^{4}-10^{5}$ yrs).

Objects initially in resonance have the potential to be lost since the
semimajor axis of a resonance will take a random step along with the planet,
and the object might not have enough time to adjust to remain in resonance. To
avoid loss, the total cumulative change in the planet’s semi-major axis due to
the random walk cannot exceed half of a resonant libration width, $\delta
a_{\rm p,lib}$. We comment on the relationship between each individual random
walk step, the cumulative distance of the random walk (which is superimposed
on the average migration), and the total distance migrated in Sections 2.3 and
3.2. For a resonant planetesimal with an eccentricity of $e_{\rm res}$,

$\delta a_{\rm p,lib}\approx 2\mathcal{C}_{\rm lib}a_{\rm
p}\left(\frac{M_{\rm{p}}e_{\rm res}^{J_{3}}}{M_{*}}\right)^{1/2}\,\,,$ (1)

where $\mathcal{C}_{\rm lib}$ is a constant that depends on the resonance,
$J_{3}$ is the order of resonance, $a_{\rm p}$ and $M_{\rm p}$ are the
semimajor axis and mass of the planet, respectively and $M_{*}$ is the mass of
the Sun (for an expository derivation, see the textbook by Murray & Dermott,
1999).222Comparing the analytical libration width for the 3:2 MMR with the
numerical width of the resonance in panel 1 of Figure 2 in Balaji et al.
(2023) produced comparable results. The libration width, $\delta a_{\rm
p,lib}$, depends on the resonance, the eccentricity of the planetesimal, and
the order of the resonance, such that loss is not uniform across the
resonances. Rather, it is preferential: different mean motion resonances will
lose more or fewer planetesimals during planetary migration, making these
population relative-loss ratios a great tracer for planetesimal-driven
migration.

Before proceeding to consider the distance of a planet’s random walk, first we
introduce a few equations and define variables. A planet’s Hill radius is

$R_{\rm H}=a_{\rm p}\left(\frac{M_{\rm p}}{3M_{*}}\right)^{1/3}\,\,,$ (2)

its Hill eccentricity is

$e_{\rm H}\equiv R_{\rm H}/a_{\rm p}\,\,,$ (3)

and its Hill velocity is

$v_{\rm H}\equiv R_{\rm H}\Omega_{\rm p}\,\,,$ (4)

where $\Omega_{\rm p}=(GM_{*}/a_{p}^{3})^{1/2}$ is the orbital angular
velocity of the planet. The random velocity of a planetesimal, $u$, with an
eccentricity of $e$ is

$u\sim e\Omega a\,\,,$ (5)

where $\Omega$ is the planetesimal’s orbital angular velocity and $a$ is its
semimajor axis. Equation (5) approximates the difference in velocity between
an eccentric and a circular orbit at the same semi-major axis, often called
the epicyclic velocity. The eccentricity and random velocity are often used
interchangeably in dynamics.

We introduce a dimensionless number of order unity, $\mathcal{M}$, which
serves to parameterize the surface density of planetesimals in the disk,
$\Sigma_{m}$, such that when integrated over the surface, the total mass of
the disk is $\mathcal{M}M_{\rm p}$ distributed among planetesimals of mass
$m$. The surface density is parameterized for one mass $m$ because a single
mass (typically the largest planetesimals) generally produce the most noise
(MC2006). The mass contained in the disk is uniformly distributed between
$a_{\rm d}/2$ to $3a_{\rm d}/2$, where $a_{\rm d}$ is the semi major axis at
the middle of the disk, and so

$\Sigma_{m}=\frac{\mathcal{M}M_{\rm p}}{2\pi a_{\rm d}^{2}}\,\,.$ (6)

We keep $a_{d}=26$au constant throughout this paper to keep a consistent
translation between $\mathcal{M}$ and $\Sigma_{m}$. We note that the real
planetesimal disk need not have this simplified global structure. The random
walk of the planet depends primarily on the surface density local to the
planet (i.e. within several $R_{H}$) since that will dictate the encounter
rate. We use Equation (6) to intuitively relate this local $\Sigma_{m}$ to the
amount of solid mass needed to form the planet. We define $x\equiv a-a_{\rm
p}$ where $\mathcal{R}$ is this distance normalized by $R_{\rm H}$ such that

$\mathcal{R}\equiv\frac{|x|}{R_{\rm H}}\,\,.$ (7)

###  2.2 Single encounters

We separate encounters between the planet and a planetesimal into non-crossing
and crossing depending on whether the orbits of the two bodies cross one
another. In a non-crossing orbit, the planet and planetesimal are far enough
from each other, such that $|x|>ae$. For these values of $x$, the relative
speed of the planetesimal is dominated by the Keplerian shear velocity
$\sim$$\Omega_{\rm p}x$, such that the encounter time is approximately a
dynamical time, $\Omega_{\rm p}^{-1}$. The expected change in $a_{\rm p}$ for
a non-crossing orbit is

where $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are order-unity numerical
factors that were not determined during derivation in MC2006.

Figure 1:  Numerical simulations (circles) validate analytical Equation
(LABEL:NC1,LABEL:NC2) using $\mathcal{C}_{1}=6.27$ and $\mathcal{C}_{2}=2.5$
(lines) for the change in a planet’s semi-major axis, $a_{\rm p}$, during a
non-crossing encounter. We fit the lines to the simulations with initial
$x/\mathcal{R}_{H}\leq 25$ since beyond that, the approximation that
$x=a-a_{\rm p}\ll 1$ begins to lose validity. We perform three sets of
simulations, all with planetary mass $M_{\rm p}=10^{-5}M_{*}$ and semi-major
axis $a_{\rm p}=30$AU. A planetesimal of mass $m$ and eccentricity $e$ starts
with semi-major axis $a_{\rm p}+x$. NC(a) (top panel) verifies Equation
(LABEL:NC1) with $e=0$ and $m=5\times 10^{-10}M_{*}$. NC(b) and NC(c) (bottom
panel) verify Equation (LABEL:NC2) with $e=0.01$ for masses $m=5\times
10^{-10}M_{*}$ (b) and $m=10^{-12}M_{*}$ (c). Figure 2: Simulations (circles)
match analytical Equations (8j) and (8l) (lines) for the change of a planet’s
semi-major axis during a crossing encounter when setting
$\mathcal{C}_{3}=2.5$. For both simulations, planetesimal semi-major axis,
eccentricity, and mass are $a$ = 25AU, $e=0.4$, and $m=10^{-10}M_{*}$,
respectively. The planet starts at semi-major axis $a_{\rm p}=28$ au, with
eccentricity $e=0$ and masses $M_{\rm p}=10^{-5}M_{*}$ for simulation set C(a)
and $M_{\rm p}=5\times 10^{-6}M_{*}$ for C(b). At impact parameters, $b$ less
than $b_{\rm min}$ (Eq. 8k), $\Delta a_{\rm p}$ is constant, as predicted by
the maximum values given by Equation (8j). Beyond that, $\Delta a_{\rm p}$
falls off inversely proportionally to $b$, as in Equation (8l).

To find $\Delta a_{\rm p,n}$ numerically, we run simulations with the Sun, a
planet, and a planetesimal. For the simulation to begin with the minimum
amount of influence experienced by each body on the other, the planet and
planetesimal begin at opposite sides of the Sun. The simulation runs for a
synodic period so that a close encounter will occur and enough time after the
encounter is left for the change in the semi-major axis to settle to its final
value. We test three non-crossing orbits, verifying the two regimes in
Equation (2.2), as seen in Figure 1. All simulations have $M_{\rm
p}=10^{-5}M_{*}$, $a_{\rm p}=30$AU, and planetary eccentricity $e_{p}=0$.
Simulation set NC(a) has $e=0$, appropriate for Equation (LABEL:NC1), while
NC(b) and (c) have $e=0.01$, appropriate for Equation (LABEL:NC2).
Planetesimals are given mass $m=5\times 10^{-10}M_{*}$ for NC(a) and (b) or
$m=10^{-12}M_{*}$ for (c) and have initial semi-major axes $a_{\rm p}+x$ for a
range of values of $x$ shown in Figure 1. We determine the numerical factors
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ by fitting the simulation output to
Equations (LABEL:NC1) and (LABEL:NC2) with the linear least-squares method in
log space. We perform the fit to simulation outputs with planetesimals
initialized at $x/\mathcal{R}_{H}\leq 25$ since beyond that the approximation
that $x=a-a_{\rm p}\ll 1$ begins to break and the slope for $\Delta a_{\rm p}$
begins to flatten out. The three simulations had different best-fit parameters
of $C_{1}=6.27$, $C_{2}=2.17$, and $C_{2}=2.78$, respectively. For simulation
sets (b) and (c), we expect the same coefficient value, so we adopt the
average of the latter and $\mathcal{C}_{2}=2.5$.

In a crossing orbit, $|x|<ae$, such that $u>\Omega x$ so the relative velocity
of the planetesimal is dominated by its epicylic velocity and its encounter
time is less than $\Omega^{-1}$. The impact parameter $b$ can be rather
different from $|x|$, and $b$ lies in the range

$b_{c}<b\lesssim ae\,\,,$ (8i)

where $b_{c}$ is the impact parameter below which a collision occurs. This
wide range can cause a wide spectrum of results for $\Delta a_{\rm p,c}$
(where $\Delta a_{\rm p,c}$ is the change of $a_{\rm p}$ due to a crossing
encounter) so MC2006 limit themselves to predict the maximum value of $|\Delta
a_{\rm p,c}|$ possible. Weaker encounters contribute to the random walk of the
planet at a comparable rate, generating a Coulomb logarithm, which we fold
into the coefficient calculated in Section 2.3. The maximum $|\Delta a_{\rm
p,c}|$ over an encounter is

$\max\left|\Delta a_{\rm p,c}\right|\sim\mathcal{C}_{3}\frac{m}{M_{\rm
p}}a_{\rm p}e\,\,,$ (8j)

generated when the impact parameter satisfies

$b\lesssim GM_{\rm p}/u^{2}\,\,,$ (8k)

Once $b$ does not satisfy Equation (8k) we expect $|\Delta a_{\rm p,c}|$ to
decrease from the value in Equation (8j), following the relationship

$\left|\Delta a_{\rm
p,c}\right|\sim\mathcal{C}_{3}\frac{Gm}{bae\Omega^{2}}\propto\frac{1}{b}\,\,.$
(8l)

As before, we find $\Delta a_{\rm p,c}$ numerically with simulations
containing the Sun, a planet, and a planetesimal. To ensure the planet and
planetesimal are in a crossing orbit, we begin with the planet on a circular
orbit with $a_{p}=28$au and place an eccentric planetesimal in the same plane.
We set the planet’s and the planetesimal’s longitude of ascending node and
argument of pericenter to be zero. Under these conditions, a collision occurs
if the true anomalies of both objects coincide simultaneously at a place where
the two orbits cross. Since the planet’s orbit is circular, this crossing must
occur at a distance from the sun $r=a_{\rm p}$, while the planetesimal’s orbit
satisfies $r=a(1-e^{2})/(1+e\cos{f})$, so that an encounter occurs at angle
$f_{\rm enc}=\cos^{-1}\left[(a(1-e^{2})-a_{p})/(a_{p}e)\right]$. We begin the
simulation at time $\delta t$ = 2.5 years before the collision, when the mean
anomalies of the objects are smaller by $\Delta\mathbb{M}_{\rm p}=\Omega_{\rm
p}\delta t$ for the planet and $\Delta\mathbb{M}=\Omega\delta t$ for the
planetesimal. Since it is on a circular orbit, the planet’s initial mean
anomaly is $\mathbb{M}_{\rm p}=f_{\rm enc}-\Delta\mathbb{M}_{\rm p}$. The
planetesimal’s initial mean anomaly is $\mathbb{M}=\mathbb{M}_{\rm
enc}-\Delta\mathbb{M}$, where we calculate the mean anomaly at encounter,
$\mathbb{M}_{\rm enc}$, from $f_{\rm enc}$ using Kepler’s equation
$\mathbb{M}_{\rm enc}=E-e\sin{E}$ for the eccentric anomaly $E=\tan(f_{\rm
enc}/2)$. Because we want close encounters rather than perfect collisions, we
add a small random value of order $10^{-4}$ radians to the planetesimal’s
initial mean anomaly and, for each simulation, we record the distance of
closest approach, $b$.

We run two sets of simulations with various impact parameters, $b$, to verify
Equations (8j) and (8l) (Figure 2). The planetesimals in both simulations have
$a=25$ au, $e=0.4$, and $m=10^{-10}M_{*}$. Simulation set C(a) has
$M_{p}=10^{-5}M_{*}$ whereas C(b) has $M_{\rm p}=5\times 10^{-6}M_{*}$, thus
highlighting the fact that the planet’s mass determines the maximum value of
$|\Delta a_{\rm p,c}|$ (Equation 8j) but not its value at large impact
parameters (Equation 8l). Our simulations validate the equations for crossing
orbits with $\mathcal{C}_{3}=2.5$ (Figure 2 dashed lines).

###  2.3 Multiple encounters

The broad-scale level of noise for a planet in a sea of planetesimals relies
on the single encounters discussed above, as the cumulative effects are
essentially a randomized distribution of those events. When multiple
planetesimals encounter the planet simultaneously, the planet and
planetesimals exchange angular momentum, and the planet migrates due to the
back reaction it experiences. We now consider this scenario, where the change
of $a_{\rm p}$ (i.e. its random walk in semi-major axis) is well modeled by
Brownian motion.

MC2006 derive the expected random walk, $\Delta a_{\rm p,T}$, of a planet in a
disk of planetesimals given by

$\Delta a_{\rm p,T}=(\dot{\overline{N}}T)^{1/2}\Delta a_{\rm p}$ (8m)

where $\dot{\overline{N}}$ is the mean encounter rate, $T$ is the migration
duration, and $\Delta a_{\rm p}$ is the change in semi-major axis of the
planet due to a single encounter as shown in the previous two subsections. For
a non-crossing orbit, $u<v_{H}$, the speed of the planetesimals relative to
the planet is determined by Keplerian shear, and impact parameter $b\sim x$
resulting in a mean encounter rate,

$\dot{\overline{N}}\sim\frac{\Sigma_{m}}{m}\Omega_{p}x^{2}$ (8n)

Planetesimals in crossing orbits encounter the planet with super-Hill
velocities ($u>v_{H}$) at impact parameters $b\lesssim GM_{p}/u^{2}$ resulting
in a mean encounter rate,

$\dot{\overline{N}}\sim\frac{\Sigma_{m}\Omega_{p}}{mu}\Big{(}\frac{GM_{p}}{u^{2}}\Big{)}^{2}u$
(8o)

where $u\sim a_{p}e\Omega_{p}$.

We summarize Equation 8m for the non-crossing and crossing cases as

$\displaystyle\Delta a_{\rm p,T}\sim$
$\displaystyle\mathcal{CR}^{-4}\left(\frac{\mathcal{M}m}{M_{\rm
p}}\right)^{1/2}\frac{a_{\rm p}}{a_{\rm d}}e_{\rm H}v_{\rm
H}\left(\frac{T}{\Omega_{\rm p}}\right)^{1/2},$ (8pa) $\displaystyle\Delta
a_{\rm p,T}\sim$ $\displaystyle\Delta a_{\rm p,T}\sim$ if $e<e_{\rm
H}/\mathcal{R}^{2}$;  $\displaystyle\Delta a_{\rm p,T}\sim$
$\displaystyle\mathcal{CR}^{-2}\left(\frac{\mathcal{M}m}{M_{\rm
p}}\right)^{1/2}\frac{a_{\rm p}}{a_{\rm d}}ev_{\rm
H}\left(\frac{T}{\Omega_{\rm p}}\right)^{1/2},$ (8pb) $\displaystyle\Delta
a_{\rm p,T}\sim$ $\displaystyle\Delta a_{\rm p,T}\sim$ if $e_{\rm
H}/\mathcal{R}^{2}<e<\mathcal{R}e_{\rm H}$;  $\displaystyle\Delta a_{\rm
p,T}\sim$ $\displaystyle\mathcal{C}\left(\frac{\mathcal{M}m}{M_{\rm
p}}\right)^{1/2}\frac{a_{\rm p}}{a_{\rm d}}\frac{e_{\rm H}^{2}}{e}v_{\rm
H}\left(\frac{T}{\Omega_{\rm p}}\right)^{1/2},$ (8pc) $\displaystyle\Delta
a_{\rm p,T}\sim$ $\displaystyle\Delta a_{\rm p,T}\sim$ if $e>\mathcal{R}e_{\rm
H}$,

where $\mathcal{C}$ is a constant that encapsulates all factors of order unity
that arose during derivations, and $e$ is the eccentricity of the
planetesimals. The non-crossing orbit scenario takes two forms (Equations 8pa
and 8pb) since $\Delta a_{\rm p}$ will vary depending on how $|\Delta e|$
compares to the pre-encounter $e$. Equation (8pc) corresponds to crossing
orbits. Note that Equation (8pc) does not depend on $\mathcal{R}$, so all
planetesimals with crossing orbits and the same eccentricities will affect the
random walk of the planet equally.

With numerical simulations, we determine the unknown coefficient $\mathcal{C}$
(where $\mathcal{C}>1$ but not $\gg 1$) and verify the total random walk of
the planet, $\Delta a_{\rm p,T}$, due to continuous interactions with
planetesimals.We run two sets of 150 simulations, each lasting one million
years with the Sun, a planet, and 9000 planetesimals. We do not know the
typical eccentricity of the disk of planetesimals at the time of Neptune’s
migration, so we initialize one set of simulations with low eccentricity
planetesimals and the other with high eccentricity planetesimals, which we
call D(a) and D(b), respectively. Still, these properties are set by the
disk’s global dynamics, and both regimes likely occur. Here, the
characterization of "low" and "high" is relative to the Hill eccentricity,
$e_{\rm H}$, since disks with $e>e_{\rm H}$ host planetesimals on crossing
orbits even if the region within a Hill radius of the planet is dynamically
cleared.

For all simulations the planet begins at 25 au with $M_{\rm p}=10^{-6}M_{*}$,
so that $e_{\rm H}=0.007$. We explore the disk’s low and high eccentricity
regimes with simulation sets D(a) and D(b) where the initial $e=0.001$ and
$e=0.05$, respectively. The planetesimals are all of equal mass $m=1.8\times
10^{-10}M_{*}$. Their semi-major axes $a$ are uniformly distributed between
$\sim$20 (21.5) au and $\sim$30 (28.5) au for D(a) (D(b)), inclinations
between 0 and the value of $e$ in radians, and remaining orbital angles are
randomly generated between 0, 2$\pi$. Gravitational interactions between the
planet and the planetesimals are simulated, but interactions between
planetesimals themselves are not. The planet and planetesimal disk have a
comparable mass, so one might wonder why neglecting planetesimal-planetesimal
interactions here is appropriate. We discuss the reason for this assumption in
Section 6.

Figure 3: Typical $a-e$ snapshot after $10^{6}$ years for one of the 150 low-
eccentricity (top) and high-eccentricty (bottom) simulations discussed in
Section 2.3. For these simulations, we place a planet within a disk of 9000
planetesimals with low eccentricities of 0.001 (top) and high eccentricities
of 0.05 (bottom). The planetesimals (grey dots) have equal mass ($1.8\times
10^{-10}M_{*}$) and are uniformly distributed in semi-major axis space between
20 au and 30 au. Their inclinations are between 0 and the value of their
respective eccentricity in radians. The planet (blue dot) begins the
simulation at 25 au with $M_{\rm p}=10^{-6}M_{*}$ and will migrate through the
disk due to multiple encounters with the planetesimals. We confirm that the
disk width is sufficiently wide to test for migration and that the systems
have reached pseudo steady-state. Excitation at some resonances occurs
(vertical lines). Figure 4: Typical random walks in semi-major axis of the
planet from 5 simualtions from the low eccentricity (light brown) and high
eccentricity (dark brown) ensembles. The planet starts off at 25 au and random
walks for the duration of the simulation. The trajectories resulting from
migration within a disk of planetesimals follow a random walk.

For Equation (2.3) to be valid, the system must reach a pseudo-steady state by
the end of the simulation, meaning over the timescale we are simulating, the
distribution of particles and surface density of the disk is not changing
substantially. We confirm that the disk is wide enough such that the
planetesimals towards the edge of the disk are not heavily excited by the
planet and the planet does not migrate outside the disk. After several
thousands of years, the system converges to a stable configuration. The
snapshot at $10^{5}$ years looks very similar to the snapshot at $10^{6}$
years. The one-million-year snapshot in Figure 3 shows the extent of influence
of the planet along with the excitation of planetesimals at positions of mean
motion resonances. In addition, we illustrate in Figure 4 that the typical
paths taken by the planet across several simulations are indeed random.

Figure 5: Histograms showing the final semi-major axis of the planet, $a_{\rm
p}$, after $10^{6}$ years for the 150 runs performed for the low eccentricity
simulation D(a) (dark blue) and high eccentricity simulation D(b) (light blue)
described in Figure 3. The planet’s initial $a_{\rm p}=25$ au for all runs.
The low eccentricity regime produces both a larger shift in the mean and a
larger spread of final positions. The mean moves inward due to dynamical
friction because our disk contains only one planet, eliminating global
asymmetries which can result in outward migration. The distribution of final
positions is well-modeled by Gaussians (dark pink lines), which for the low
eccentricity (high eccentricity) case have mean $\mu$ = 24.9585 au (24.9847
au) and standard deviation $\sigma$ = 0.0253 au (0.0103 au). This demonstrates
that the planet’s migration is well modeled as a diffusion process.

We average over 150 simulations for each regime to determine the randomness of
this migration process. At $10^{6}$ years, the planet’s final semi-major axis
for the low eccentricity (high eccentricity) case is dispersed among an
ensemble of simulations with mean value $\sim 24.9585$ au (24.9847 au) and
standard deviation 0.0253 au (0.0103 au). The mean semi-major axis moves
towards the Sun after $10^{6}$ years (Figure 5), as expected for a simulation
containing a single planet. In the Solar System, planetesimals scattered by
the outermost giant planets are preferentially removed from Neptune’s
dynamical environment by interactions with Jupiter and Saturn. As a result,
Neptune migrates outward (Fernandez & Ip, 1984). Because the random walk
arises due to Poisson noise on top of the mean number of planetesimal
encounters on one or both sides of the planet (MC2006), the magnitude of the
random walk is not expected to depend on the direction of motion. In Appendix
A, we illustrate that the average migration seen in our simulations matches
that expected for a single planet in a planetesimal disk.

For the remainder of this work, we consider only the random walk component. We
interpret the typical order-of-magnitude random walk distance $\Delta a_{\rm
p,T}$ given by Equation (2.3) as the standard deviation of the planet’s semi-
major axis across the 150 simulations, $\sigma_{a_{\rm p,\rm T}}$.

Because Equation (2.3) depends on several orbital parameters of the disk and
properties of the planet, we must disentangle which sub-population of
simulated planetesimals causes the most randomness. To do that, we choose six
regions of $a-e$ space, as shown in Figure 6, for the low eccentricity (D(a),
top) and high eccentricity (D(b), bottom) simulation sets. All parameters that
are the same across the six groups are described in Figure 3, whereas
$\mathcal{R}$, $e$, and $\mathcal{M}$ are shown in Table 1.

We calculate the surface density of both disks for each group centered at $\pm
x$ au away from the planet, where $x=\mathcal{R}R_{H}$. For groups $1-3$ and
6, the groups have radial width $2R_{H}$ and eccentricities within the limits
expressed in Equations (8pa) and (8pb). For groups 4 and 5, we include all
planetesimals in semi-major axis space appropriate for the crossing orbits,
(denoted by the dark pink dashed line in Figure 6), since Equation (8pc) does
not depend on $\mathcal{R}$. The two groups are cut off at eccentricity 0.03,
and the width used to calculate the surface densities for the groups is
$(25+|x|)e$ au. The non-dimensional surface density parameter, $\mathcal{M}$,
is then determined by Equation (6) summing $\mathcal{M}$ pertaining to the
inner and outer disk for each group. With this exercise, we wish to determine
which factor plays the larger role in producing a higher $\Delta a_{\rm p,T}$:
eccentricity of the planetesimals or proximity to the planet.

| group no. | $e$ | $\mathcal{R}$ | $\mathcal{M}$ | equation
---|---|---|---|---|---
D(a) | 1 | 0.001 ($e_{i}$) | 5 | 14.72 | 8pb
2 | 0.0069 ($e_{\rm H}$) | 3 | 2.56 | 8pb
3 | 0.001 ($e_{i}$) | 1 | 6.75 | 8pa
| 4 | 0.05 ($e_{i}$) | 1 | 6.38 | 8pc
D(b) | 5 | 0.01 ($e_{i}$) | 1 | 4.16 | 8pc
| 6 | 0.05 ($e_{i}$) | 10 | 8.56 | 8pb
Table 1: Summary of the parameters that differ across groups 1-6 (illustrated
by the colored groups in Figure 6). Groups 1-3 probe the low eccentricity
simulation ensemble, whereas 4-6 probe the high eccentricity simulation
ensemble. The planetesimals in each group occupy non-crossing and crossing
orbits, and will therefore interact with the planet differently and produce
varying magnitudes of random walk distance, determined by the equation denoted
in the last column.

Inputting the values in Table 1 and Figure 3, we calculate Equation (2.3) for
each group as a function of time and compare to the amount of migration from
simulations, as shown in Figure 7. The simulation-derived migration for low
eccentricity (dashed pink line) and high eccentricity (dashed purple line)
ensembles are best matched by the analytical curves with parameters from D(a)
group 3 and D(b) group 5 and $\mathcal{C}=3.5$. These groups have the largest
impact on the planet’s random walk (i.e. have largest values of $\Delta a_{\rm
p,T}$). In low eccentricity case (a) groups 1 and 2 underestimate the
migration by factors greater than 10, compared to group 3, so the
planetesimals in group 3 clearly dominate the migration. This is likely
because they are closely packed at one Hill radius from the planet at low
eccentricities. For the high eccentricity case (b), groups 4 and 6 are lower
by factors less than 10, compared to group 5, so we can say that planetesimals
across the disk still affect migration but not as greatly as those closer to
the planet with lower eccentricities.

Setting $\mathcal{C}$ to 3.5 for Equation (2.3) provides a good fit to
$\sigma_{a_{\rm P,T}}$ for cases (a) and (b). Note that, particularly for the
low eccentricity case, the slope of the numerical curve matches our
theoretical expectation of $\Delta a_{\rm p}\propto T^{1/2}$, both near the
beginning of the simulation and near the end, but the normalization changes at
intermediate times. At early times, the disk has not yet reached a pseudo
steady-state, in which some regions near the disk have been cleared by
scattering (see Figure 6), leading to enhanced but transient dispersion.

For the high eccentricity case, the late time slope is less convincingly
consistent with our theoretical expectation, but again, the early time slope
matches. The transition to a true pseudo steady state may simply be taking
longer in this case—longer integrations would clarify this picture. Our
results provide good verification of Equation (2.3) and convince us that
planetary migration within a disk of planetesimals can be modelled by Brownian
motion and that no cumulative effects build up, i.e each encounter is truly
random.

Figure 6:  The planetesimals’ relative distances to the planet,
eccentricities, local disk surface density, and more, greatly affect the
magnitude of the planet’s random walk. We show the same simulation snapshots
as Figure 3 to visualize the orbital differences between six regions in the
disk, to probe which region produces the higher $\Delta a_{\rm p,T}$. The six
groups–denoted by the colors (numbers) dark blue (1), light blue (2), dark
purple (3), dark green (4), light green (5), and dark pink (6)–are shown in
Table 1 and Figure 7. The planetesimals in group 3 have eccentricities
$\lesssim e_{H}/\mathcal{R}^{2}$ (blue solid line) so Equation (8pa)
determines their impact on the planet’s noisiness. Planetesimals in groups 1,
2, and 6 have eccentricities $\gtrsim e_{H}/R^{2}$ but
$\lesssim\mathcal{R}e_{H}$ (dark pink dashed line) so Equation (8pb)
determines their impact. Planetesimals in groups 4 and 5 occupy crossing
orbits with eccentricities $\gtrsim\mathcal{R}e_{H}$, and their impact on the
planet’s random walk is determined by Equation (8pc). The Hill eccentricity of
the planet is shown for reference (blue dashed line).  Figure 7: The
numerical, expected planet’s random walk, determined by the standard
deviation, $\sigma$, of the distribution of semi major axes over time for
simulation sets D(a) (light pink dashed line) and D(b) (light purple dashed
line) compared to the analytical random walk for each group shown in Figure 6
and Table 1 (solid lines), determined by Equation (2.3). The orbital
parameters and local surface density from Groups 3 (dark purple) and 5 (light
green) produce the largest $\Delta a_{\rm p,T}$ for the low and high
eccentricity case, respectively, highlighting the role of proximity of the
planetesimals to the planet and their eccentricity. The group labels are shown
in order from largest $\Delta a_{\rm p,T}$ to smallest and inline labels to
clearly compare between low and high eccentricity groups.

##  3 analysis applied to all resonances

With a confirmed random walk analytical model and normalized coefficients, we
apply the analysis from Section 2 to the Solar Systems’ early planetesimal
disk present during Neptune’s migration. The random walk of a planet is
dominated by the largest planetesimals driving migration (MC2006), so we
proceed to consider an average size for the largest planetesimals observed
today in the Kuiper Belt (Brown, 2008). We use the analytics from Section 2 to
find the characteristic times for loss from various resonances with Neptune
(Section 3.1) and then calculate the retention fractions for these resonances
(Section 3.2). We analyze the following characteristic loss times and
retention fractions as a function of $\mathcal{M}=\mathcal{M}_{\rm LG}$, which
represents the surface density of the disk in planetesimals with radius equal
to 780 km and $\rho=2\rm g/cm^{2}$, resulting in a mass that is roughly one-
third of Pluto’s mass. We choose this radius because it represents the average
of the largest 8 observed TNOs in our solar system from Brown (2008). The
characteristic loss times and retention fractions can give us an upper limit
migration timescale constraint before objects are lost from resonance.

For our analysis, we assume an early planetesimal disk with total mass given
by $\mathcal{M}=2$, which corresponds to a total of $\gtrsim 30$ M⊕.
Planetesimal disk masses between 20 to 35 $M_{\oplus}$ are typical for models
of Neptune’s migration (Luu & Jewitt, 2002; Nesvorný, 2018). Within this
framework, $\mathcal{M}_{\rm{LG}}=2$, where $\mathcal{M}_{\rm{LG}}$ is the
mass in large planetesimals, would imply that 100% of the mass in the disk is
contained in large objects, while smaller values of $\mathcal{M}_{\rm LG}$
imply that only a fraction $\mathcal{M}_{\rm LG}/\mathcal{M}$ of the total
mass is in these large bodies. Thus, the value of $\mathcal{M}_{\rm LG}$ for a
given time is directly related to the size distribution of the planetesimal
disk. In addition to $\mathcal{M}_{LG}$, we do this analysis for
$\mathcal{M}_{PL}$ and $\mathcal{M}_{2PL}$ corresponding to Pluto size and
twice the size of Pluto objects. We discuss constraints for the size
distribution in relation to retention fractions of MMR during the epoch of
Neptune’s migration for values of $\mathcal{M}_{LG}$, $\mathcal{M}_{PL}$, and
$\mathcal{M}_{2PL}$ in Section 4.

###  3.1 Ranking of loss times

$p:q$ |  $a_{\rm res}$ (au) | $\mathcal{C}_{\rm lib}$ | $e_{\rm res}$ |  $e_{\rm res}$ reference
---|---|---|---|---
3:2 | 39.4 | 3.64 | 0.175 | Volk et al. (2016)
11:7 | 40.7 | 23.66 | 0.18 |  $q$-fit
8:5 | 41.2 | 11.16 | 0.19 |  $q$-fit
5:3 | 42.3 | 5.51 | 0.16 | Gladman et al. (2012)
7:4 | 43.7 | 9.03 | 0.12 | Gladman et al. (2012)
9:5 | 44.5 | 15.39 | 0.22 |  $q$-fit
2:1 | 47.8 | 3.00 | 0.275 | Chen et al. (2019)
7:3 | 53.0 | 8.87 | 0.30 | Gladman et al. (2012)
5:2 | 54.4 | 5.26 | 0.40 | Volk et al. (2016)
3:1 | 62.5 | 3.25 | 0.48 | Alexandersen et al. (2016)
4:1 | 75.7 | 3.64 | 0.57 | Alexandersen et al. (2016)
Table 2: Mean motion resonances up to forth order lying between the 3:2 and
4:1 resonances (denoted by $p:q$), listed in order of ascending semi major
axis, $a_{\rm res}$. Values of $\mathcal{C}_{\rm lib}$ are provided for
reference (see Murray & Dermott (1999)). For resonances with a citation listed
in column five, we quote the average value of eccentricities of observed
resonant TNOs, $e_{\rm res}$, which may be affected by selection effects. The
remaining resonances’ eccentricities were determined by fitting the data to
$e_{\rm res}=1-q_{\rm free}/a_{\rm res}$, $q_{\rm free}$ is a free parameter
we fit. The best-fit pericenter distance, $q_{\rm free}$, value was then used
to set a typical $e_{\rm res}$.

We calculate the maximum migration timescale Neptune can migrate before losing
objects in resonance given that a planet must random walk a distance less than
half the libration width of the resonance in order to maintain planetesimals
in a mean motion resonance ($\Delta a_{\rm p,T}<\delta a_{\rm p,lib}/2$). We
call this timescale, the resonance loss timescale, $T_{\rm loss}$, which can
also be interpreted as the migration duration for which a planet will random
walk a total of more than $\delta a_{\rm p,lib}/2$, thus causing loss from
resonances for given parameters of the disk. Using Equation (1) and Equation
(8pa), appropriate for $e<e_{\rm H}/\mathcal{R}^{2}$ (other regimes produce
analogous results), the resonance loss timescale is

$T_{\rm loss}\sim\mathcal{R}^{8}\left(\frac{\mathcal{C}_{\rm
lib}}{\mathcal{C}}\right)^{2}\frac{M_{\rm
p}^{2}}{\mathcal{M}M_{*}m}\frac{e_{\rm
res}^{J_{3}}}{e_{H}^{2}}\left(\frac{a_{\rm d}}{v_{\rm
H}}\right)^{2}\Omega_{\rm p}\,\,.$ (8aj)

Figure 8 presents $T_{\rm loss}$ for our collection of resonances, evaluated
using Equation (8aj) with $\mathcal{R}=1$, $e=0.01$, $a_{p}=a_{d}=26.6$ au,
$m=0.3M_{\rm PL}$, $M_{\rm p}=M_{\rm N}$, and $M_{*}=M_{\odot}$, where $M_{\rm
PL}$, $M_{\rm N}$ and $M_{\odot}$ are the masses of Pluto, Neptune and the
Sun, respectively, so that $R_{\rm H}\sim 0.69$ au and $e_{\rm H}\sim 0.026$.
We set $\mathcal{C}=3.5$, the value determined in Section 2.3, and plug in
typical eccentricities of TNOs in MMR for $e_{\rm res}$ (see Table 2). We do
this study for resonances up to 4th order that lie between the 3:2 and 4:1.
For the resonances with enough observations to produce a typical eccentricity
range for the TNOs, we quote their Gaussian centered or average value in Table
2 based on surveys done by Gladman et al. (2012); Volk et al. (2016);
Crompvoets et al. (2022); Bannister et al. (2018); Alexandersen et al. (2016);
Chen et al. (2019); Volk et al. (2016). We determine $e_{\rm res}$ for the
resonances missing observed eccentricities by assuming the resonances possibly
have a similar pericenter distribution, $q_{\rm free}$, motivated by the
observed eccentricities increasing as a function of semi-major axis. Similarly
stable libration zone studies in $(a,e)$ space for Neptune’s exterior MMRs
show that the widest part of the libration width (i.e. where more TNOs are
expected to reside) increases along a line of constant pericenter (Lan &
Malhotra, 2019). Fitting the data in Table 2 to $e_{\rm res}=1-q_{\rm
free}/a_{\rm res}$ produced a best-fit value $q_{\rm free}=\sim 34.5$ au,
resulting in appropriate $e_{\rm res}$ values equal to 0.18, 0.19, and 0.22
for the 11:7, 8:5, and 9:5 MMRs, respectively.

Because the noise due to the planet’s random walk is dominated by the largest
planetesimals driving migration, we show $T_{\rm loss}$ as a function of
$\mathcal{M}=\mathcal{M}_{\rm LG}$. The MMRs requiring the longest time to
lose objects (2:1, 3:1, 4:1, 3:2; Figure 8) can be interpreted as being the
strongest resonances, whereas the 7:3, 11:7, 9:5, and 7:4 are the weakest. The
order in Figure 8 matches the order of $\delta a_{\rm p,lib}$ for the values
in Table 2, which is in agreement with the width of a resonance being
positively correlated with the strength of the resonance (Gallardo, 2019). The
five strongest resonances in Figure 8 correspond to the resonances with
highest width in semi-major axis in Lan & Malhotra (2019), whereas the rest
differ in ranking.

All parameters in Equation (8aj) are the same across the disk except those
relating to the resonance. The strength of the resonance therefore depends
highly on $C_{\rm lib}$, $e_{\rm res}$ and $J_{3}$. Since $T_{\rm
loss}\propto(\delta a_{\rm p,lib})^{2}\propto e_{\rm res}^{J_{3}}$, a small
change in $e_{\rm res}$ for higher order resonances ($J_{3}>1$) will
significantly affect the ranking of the resonances.

The observed eccentricities $e_{\rm res}$ used to calculate $T_{\rm loss}$ are
averages or Gaussian-centered values from observations (see Table 2). TNOs are
harder to detect the further their semi-major axes and pericenters are from
the Sun and so these observational biases make discovering low-eccentricity
distant TNOs difficult. In addition, $e_{\rm res}$ for the resonances without
constrained values might deviate from the pericenter fit we calculated. As a
result, the ranking of resonance strength in Figure 8 is an estimate, and
varying $e_{\rm res}$ would change the ordering of the resonances.

We highlight that because we pick the eccentricities of objects observed to
occupy each resonance today, we know that those objects have not been lost,
and we can use this verified resonance retention to place limits on the
planetesimal population driving migration. However, it may well be that
resonances—particularly higher-order resonances—preferentially lost particles
at low eccentricities, so that particles are not observed at those
eccentricities today. This behavior generates another pattern that can be
explored via comparison with future observational data better able to
constrain the low-eccentriciy populations of distant resonances. For example,
the third-order 5:2 resonance has a large observed population composed of
bodies with large typical eccentricities of 0.4 (Gladman et al., 2012; Volk et
al., 2016). As illustrated in Volk et al. (2016), the low-eccentricity
population of the 5:2 MMR is not well-constrained given current data. At
observed eccentricity of 0.2, the 5:2 MMR is significantly weaker than at 0.4
(i.e., the migration timescale decreases by $2^{J_{3}}$ where the order of the
resonance, $J_{3}$, is 3, for a given $\mathcal{M_{\rm LG}}$). Though the
resonance width decreases sharply at low eccentricities, a low-eccentricity
population in the 5:2 resonance is expected given its similar stable libration
zone to the 3:2 and 2:1 MMRs (Malhotra et al., 2018). If future observations
indicate that the 5:2 MMR does not host a low-eccentricity population, then
either (1) it never had one because, for example, the population originated
from a planetary gravitational upheaval that scattered planetesimals along
trajectories with pericenters in the scattering region and the pericenter
distances of low-eccentricity 5:2 resonant objects are too large, or (2)
stochasticity during migration caused preferential loss of low-eccentricity
objects in the 5:2 MMR. Conveniently, these two scenarios predict different
patterns in the current eccentricities of objects in distant resonaces. In
scenario (1), $e_{\rm res}=1-q/a_{\rm res}$ for a resonance at semi-major axis
$a_{\rm res}$ hosting an object with pericenter distance $q$. While resonant
libration and chaotic diffusion smear observed pericenter distances over time,
observed eccentricities in resonance are expected to depend primarily on the
semi-major axis of the resonance. In scenario (2), loss of low-eccentriciy
objects due to stochasticity is stronger for higher-order resonances,
generating a pattern that depends strongly not only on a resonance’s location,
but also on its order. Future characterization of distant resonances at low
eccentricities will enable this analysis.

Changing the eccentricity of the planetesimals near the planet driving
migration— instead of the resonant particle—does not affect the ranking of the
resonances. Instead, the absolute values of the loss times change by a
constant factor given a fixed set of disk conditions (this factor is simply a
combination of $\mathcal{R}$ and $e_{\rm H}/e$).

The large observed population in the 3rd order 5:2 resonance and the weakness
and small population of the 7:3 resonance allow us to place new constraints on
planetesimal-driven migration. The timescale and amount of Neptune’s migration
are confined to a scenario where the weak 7:3 MMR population would not be
totally lost and where many objects in the 5:2 MMR would still be retained.
Figure 8 shows that the 7:4 MMR is actually the weakest, but we do not
consider this resonance in our study due to its proximity to the cold
classical region and the ease with which it can pick up its neighbors for
short timescales (Lykawka & Mukai, 2005, 2007). As a result, the 7:4 MMR might
not be entirely or even dominantly populated by an era of planetesimal-driven
migration. We also do not consider the 11:7 and 9:5 MMRs since they do not
have constrained population estimates. We discuss constraints due to the 7:3
and 5:2 resonances in Section 4. This discussion merits renewed attention as
we get more data from OSSOS (Bannister et al., 2018), DES (Bernardinelli et
al., 2022) and the highly anticipated LSST (LSST Science Collaboration et al.,
2009), at which point the evolution of $e_{\rm res}$ during migration and
transiently stuck populations (Yu et al., 2018) should also be considered.

Figure 8: Characteristic resonance loss time, $T_{\rm loss}$, as defined by
Equation (8aj) for the exterior mean motion resonances of Neptune that we
study (see Table 2), as a function of $\mathcal{M}=\mathcal{M}_{\rm LG}$ (see
Equation 6). For this calculation, $\mathcal{R}=1$, $e=0.01$,
$a_{p}=a_{d}=26.6$ au, $m=0.3M_{\rm PL}$, $M_{\rm p}=M_{\rm N}$, and
$M_{*}=M_{\odot}$, where $M_{\rm PL}$, $M_{\rm N}$ and $M_{\odot}$ are the
masses of Pluto, Neptune and the Sun, respectively, so that $R_{\rm H}\sim
0.69$ au and $e_{\rm H}\sim 0.026$. The eccentricity of particles in
resonance, $e_{\rm res}$, is set to the typical value observed today for each
resonance or a value determined by fitting the data to $e_{\rm res}=1-q_{\rm
res}/a_{\rm res}$ for those without measurements (see Table 2 and text).

###  3.2 Retention Fraction

We now translate our results from Section 3.1 into the fraction of initially
captured objects that are retained in each resonance as a function of disk
properties and the total migration time. MC2006 provide analytical solutions
for the retention fraction produced by a disk with uniform planetesimal
masses. The probability that a captured TNO is retained in resonance over the
duration $T$ of migration is

$P_{\rm
keep}=\sum_{n=1}^{\infty}\frac{4}{n\pi}\sin{{}^{3}\left(\frac{n\pi}{2}\right)}e^{-\lambda_{n}T}\,\,,$
(8ak)

where $\lambda_{n}=(n\pi)^{2}D/(2\delta a_{\rm p,lib}^{2})$. The diffusion
coefficient describing the cumulative random-walk change in semimajor axis,
$\Delta a_{\rm p,T}$ due to many encounters over a total time, $T$ is

$D=\frac{\Delta a_{p,T}^{2}}{T}\,\,,$ (8al)

As in all diffusive processes, the diffusion coefficient relates this large-
scale behavior to the individual random steps in semi-major axis $\Delta
a_{\rm p}$ resulting from single encounters spaced in time by an average
$\Delta t\equiv 1/\dot{\overline{N}}$, where $\dot{\overline{N}}$ is the mean
encounter rate given by Equations (8n) and (8o) for non-crossing and crossing
orbits. The diffusion coefficient for the small-scale case is then
$D\sim\mathcal{C}_{4}(\Delta a_{\rm p})^{2}/\Delta t$. For our results to be
self-consistent, we expect the coefficient $\mathcal{C}$ in Equation (2.3) to
be the same for all three regimes since this choice generates smooth
transitions at regime boundaries. Equating the large-scale and small-scale
diffusion coefficients requires setting $\Delta a_{\rm p,T}=\Delta a_{\rm
p}\sqrt{\mathcal{C}_{4}(T/\Delta t)}$ equal to Equation (2.3) for each regime
and solving for $\mathcal{C}_{4}$. We find that
$\mathcal{C}_{4}=\Big{(}\frac{\mathcal{C}}{\mathcal{C}_{1}}\Big{)}^{2}\frac{2\pi}{81}$,
$\Big{(}\frac{\mathcal{C}}{\mathcal{C}_{2}}\Big{)}^{2}\frac{2\pi}{9}$, and
$\Big{(}\frac{\mathcal{C}}{\mathcal{C}_{3}}\Big{)}^{2}\frac{2\pi}{9}$ for
cases (a), (b), and (c), respectively. Since we found that $\mathcal{C}_{1}$
is 2.5 times larger than $\mathcal{C}_{2}$ and $\mathcal{C}_{3}$ (see Section
2.2), these results verify that $\mathcal{C}$ is approximately the same for
all three regimes. This assures us that the normalization for the micro-scale
diffusion coefficient is correct and this process is in fact well modeled by a
diffusion process.

In Figure 9, we illustrate the variation in $P_{\rm keep}$ across the
resonances we study for the case $e\sim e_{\rm H}$, corresponding to a
diffusion coefficient calculated with Equation (8pb). Inputting values that
correspond to 2000 Plutos in the primordial planetesimal disk produces varying
amounts of loss across the MMRs, with the 7:4 MMR completely obliterated (top
panel). Inputting the same values except $e_{\rm res}=0.25$ (bottom panel)
produces different strengths. In particular, the 7:4 MMR becomes stronger than
the other weakest resonances for lower $e_{\rm res}$. The 7:3 MMR and in
particular the 5:2 MMR produce significantly lower retention fractions for
smaller $e_{\rm res}$, highlighting the correlation with $e_{\rm res}$. If
Neptune traveled in a disk with $e\leq e_{\rm H}$, we expect the same
retention fraction behavior since Equations (8pa) and (8pb) are the same when
$e=e_{\rm H}$ and $\mathcal{R}=1$, thus producing the same amount of noise.
For the case where Neptune migrates across a hot disk (i.e. $e>e_{H}$), the
timescale to lose the same fraction of objects as the first two cases
increases by a factor of $e^{2}/e_{H}^{2}$.

Increasing the number of large objects in the disk (e.g., increasing
$\mathcal{M_{\rm LG}}$ by assuming a more top-heavy size distribution)
produces more noise, and fewer objects are retained in resonance (Figure 10).
As the timescale for migration increases from 10 Myr (top panel) to 50 Myr
(bottom panel), the total accumulated noise increases, and significantly more
loss is seen for a given $\mathcal{M}_{LG}$.

Figure 9: Retention probabilities predicted by Equation 8ak for Neptune’s
exterior mean motion resonances (between panels) for the recorded $e_{\rm
res}$ (top panel) and $e_{\rm res}=0.25$ (bottom panel). For this calculation,
$T=10$ Myr $\mathcal{R}=1$, $e=e_{H}=0.026$, $a_{p}=a_{d}=26.6$ au, $m=M_{\rm
PL}$, $M_{\rm p}=M_{\rm N}$, and $M_{*}=M_{\odot}$, and $\mathcal{M}=0.22$
which corresponds to 2000 Plutos in the planetesimal disk. Varying $e_{\rm
res}$ changes the strength of the resonances compared to each other, thus
highlighting the high dependence on the eccentricity of the objects in MMR in
understanding the expected loss from each resonance. Note that if all
resonances have a retention fraction of 1, this does not mean that all
resonances have equal populations; just that the same proportion is retained.
Figure 10: Retention fractions spanning 0 (white) to 1 (dark green) for
Neptune’s exterior mean motion resonances (between panels) as a function of
$\mathcal{M}_{LG}$ for migration timescales 10 million years (top panel) and
50 million years (bottom panel). $P_{\rm keep}$ is calculated with the
parameters described in Figure 9 (top panel) for a range of surface densities
$\mathcal{M}_{LG}$ representing a planetesimal disk with the largest objects
having $m=0.3m_{\rm PL}$. Increasing the number of large planetesimals in the
disk produces more noise, causing a smaller retention fraction (each resonance
gradient gets lighter for higher $\mathcal{M}_{\rm LG}$). A similar effect is
seen when increasing the migration timescale; more time allows for noise to
accumulate, causing more objects to be lost from the resonances.

##  4 Influence of Planetesimal Size Distribution

For a given size distribution of planetesimals, noisiness in the random walk
of a planet is dominated by those with maximum $Nm^{2}$ where $N$ is the
number of planetesimals and $m$ is the mass of an individual planetesimal
(c.f. Equation (2.3), noting that $N\propto\mathcal{M}/m$). It turns out that
larger-radius planetesimals, $s$, produce the highest $Nm^{2}$. We summarize
why that’s the case here. The Kuiper Belt size distribution is described by a
power law, $dN/ds\propto s^{-q}$, and integrating that gives $N(s)\propto
s^{1-q}$. Assuming spherical, uniform-density TNOs, then $N(s)m^{2}\propto
s^{7-q}$ so that for $q<7$ (typical of Kuiper Belt size distributions)
$Nm^{2}$ is larger for larger $s$. Given that few TNOs have direct size
estimates, their measured absolute ($H$) magnitude distribution is typically
discussed instead, where $dN/dH\propto 10^{\alpha H}$. The size distribution
exponent $q$ is related to $\alpha$ by $q=5\alpha+1$.

The number of large planetesimals in Neptune’s vicinity at the time of its
migration depends both on the total mass in planetesimals at that time and on
their size distribution. The observed mass in Kuiper belt objects (e.g.,
Chiang et al., 2007) is $\lesssim$1% of the likely initial mass in solids in
the Neptunian region. Given uncertainties in the evolution of the population
of large TNOs, we consider two end-member possibilities: (1) the surface
density of TNOs was depleted in a size-independent way, perhaps through
dynamical ejection during the epoch of planetary migration, so that the
current observed size distribution matches the distribution at the time of
migration; and (2) smaller TNOs were cleared from the outer solar system more
easily, perhaps through grinding down to dust (e.g., Pan & Sari, 2005), and
the large TNOs observed today are the only large objects that ever occupied
the migration region. These possibilities provide and upper limit and lower
limit to the expected stochasticity resulting from planetesimal interactions,
respectively.

We note that planetesimal-driven migration and damping of the planet’s
eccentricity result from interactions with all planetesimals and are not
dominated by large bodies that house a minority of the mass. Thus, the
stochastic component of the migration cannot be directly inferred from the
distance of planetesimal-driven migration or amount of eccentricity damping
without knowledge of the planetesimal size distribution.

###  4.1 Current size distribution

Measurements for the absolute magnitude distribution of the Kuiper Belt are
often separated into dynamically "cold" and "hot" populations or MMR
populations (Fraser et al., 2014; Adams et al., 2014; Lawler et al., 2018;
Kavelaars et al., 2021; Petit et al., 2023), and recent work demonstrates that
a full accounting of each population may require several power law slopes with
several breaks (Morbidelli & Nesvorný, 2020). We choose the broken power law
$H$ distribution for the dynamically hot population provided by Lawler et al.
(2018), to comment on constraints for Neptune’s migration, motivated by the
idea that the hot population is most likely to have scattered from the
planet’s vicinity during migration.

The $H$ distribution for the dynamically hot population from Lawler et al.
(2018) has a steep slope ($\alpha=0.9$) for large objects (low-$H$) and a
break magnitude $H_{b}=8.3$, corresponding to a radius of 75 km (assuming an
albedo, $a=0.04$). Most of the mass resides in planetesimals with sizes near
the break. We calculate the parameterized mass in large objects,
$\mathcal{M_{\rm LG}}$, for this size distribution with a simple
approximation:

$\displaystyle\mathcal{M}M_{\rm disk}=M_{\rm total}=\mathcal{C}s_{b}^{4-q}$
(8am) $\displaystyle\mathcal{M_{\rm large}}M_{\rm disk}=M_{\rm
large}=\mathcal{C}s_{m}^{4-q}$

such that

$\displaystyle\mathcal{M_{\rm large}}=\mathcal{M}s_{m}^{4-q}/s_{b}^{4-q}$
(8an)

where $\mathcal{M}=2$ (appropriate for $\sim$ 30 M⊕ in the disk), $s_{b}$ is
the radius at the break of the broken power law (75 km) and $s_{m}$ is the
maximum radius covered in the distribution (780 km), and $q$ is 5.5 for sizes
ranging from the break to the maximum. As a result, $\mathcal{M}_{\rm LG}$ for
this size distribution is 0.060, which transaltes to 3% of the mass in the
disk in large (780 km) planetesimals (or 1528 objects)

If the primordial planetesimal disk had the same size distribution as today’s
hot TNOs, Neptune’s noisiness during 80 Myr migration would cause the 7
strongest resonances (as shown in Figure 8) to retain 70-80% of their resonant
objects, while 7:3, 5:3, and 9:5 retain between 30-45% of objects. The 7:4 MMR
loses objects extremely fast compared to other resonances, so we do not
include this resonance when discussing "the weakest resonances" further (see
Figure 10). A longer migration timescale produces more loss across the
resonances. For the 7:3 MMR to retain at least $\sim$30% of its objects, this
size distribution results in a maximum migration timescale of 80 Myr.

The retention fraction is highly dependent on the mass of the planetesimal, so
we do a similar analysis as described in Equations (8am) and (8an) for a
planetesimal mass $m=m_{\rm pluto}$ and $s_{m}=1100$ km. Considering the same
size distribution again, now $\mathcal{M_{\rm PL}}=0.036$, which translates to
1.5% of the total disk’s mass in Pluto-sized objects (i.e. 325 objects). For
this scenario, the same retention behavior as discussed above occurs, but for
a migration timescale of 40 Myrs.

###  4.2 Number of Large TNOs

We switch gears to considering a scenario in which the large TNOs observed
today are the only large objects that existed in the disk during Neptune’s
migration. The fraction of large objects, $\mathcal{M}_{LG}/\mathcal{M}$, is
simply the fraction of the total mass in large objects over the total mass of
the disk. Constraints for the timescale and number of large objects in the
disk are confined to scenarios where the weakest 7:3 MMR population would not
be totally lost given its estimated current population (see Figure 11 and
Section 5 for further discussion). If the disk of planetesimals had the same 8
largest objects from Brown (2008), corresponding to an average radius $s=780$
km and $\rho=2\rm g/cm^{2}$, then $\mathcal{M}_{\rm LG}=0.0003$. This surface
density does not cause any loss from the resonances for reasonable migration
timescales. To drop a significant fraction of objects in the 7:3 during an
epoch of Neptune’s planetesimal-driven migration, the disk would need several
hundred times more planetesimals with a radius of 780 km, than what exist
today.

According to the minimum mass solar nebula model, $10^{2}-\sim 10^{3}$ times
more mass existed during early Solar System formation than now (Chiang et al.,
2007). If the expulsion of material out of the outer Solar System was due to a
mass-independent chaotic, gravitational upheaval, then one can expect up to
$\sim 10^{3}$ times more large TNOs than those that exist today. Even in the
absence of a dynamical upheaval event, long-term dynamical evolution leads to
some depletion of TNOs in MMR over time (see studies for 3:2 and 2:1, e.g.,
Tiscareno & Malhotra 2009; Balaji et al. 2023). We comment on the relative
retention fractions for long term Nbody effects and Neptune’s noisy migration
in Section 6.

In addition, there may exist more than 8 large TNOS that simply remain
undiscovered. We find the minimum number of large, 780-km objects needed to
produce significant noise-induced loss from the MMRs during planetesimal-
driven migration. Though only a couple of Pluto-sized TNOs currently are
known, we perform this calculation for planetesimals with a Pluto size or
twice the size of Pluto (hereafter twice-Pluto) to provide a larger range of
constraints. A few thousand 780 km objects will produce the similar noise, and
thus loss from resonance, as hundreds of Plutos and tens of twice-Plutos.
Specifically, the timescale to lose the same fraction of objects for a disk
extended from a size $s_{m,1}$ to $s_{m,2}$ will be
$(\mathcal{M}_{2}/\mathcal{M}_{1})/(m_{1}/m_{2})$ where the masses $m$ are
proportional to $s^{3}$.

The migration timescale and number of large objects that would deplete the 7:3
MMR up to 99%, in strong disagreement with observations, is as follows: 4500
780-km planetesimals with $T=60$ Myr, 3000 Plutos with $T=10$ Myr, and 200
twice-Plutos with $T=10$ Myr. Requiring 25-30 % retention for the 7:3 MMR
produces stronger constraints: 2000 780-km object with $T=60$ Myr, 1333 Plutos
for $T=10$ Myr, 80 twice-Plutos for $T=10$ Myr. Note that increasing the
number of large objects by 2 and decreasing the migration timescale by the
same factor produces the same results since the argument in the exponential in
$P_{\rm keep}$ ($\lambda_{n}T$ ) is proportional to $\mathcal{M}T$.

Figure 11:  The debiased absolute population estimates of resonant TNOs
brighter than $H_{g}=9.16$ and $H_{r}=8.66$ (light blue lines with 95%
confidence error bars) provided by Gladman et al. (2012) (5:3, 7:3, 7:4), Volk
et al. (2016) (3:2, 2:1, 5:2), and Crompvoets et al. (2022) (3:1, 4:1)
compared to noise-reduced populations due to 2000 Plutos around Neptune as it
migrates for 10 Myr (circles). We consider two scenarios for filling the mean
motion resonances with an initial number of objects $N_{i}(a)$: (1) a
scattering scenario that fills $a-e$ phase space with planetesimals over a
region corresponding to orbital pericenters in the scattering region (top and
bottom panel) and (2) a substantial amount of “smooth" planetesimal-driven
migration moved Neptune outward, picking up objects into the various
resonances (middle panel). Scenario (1) requires more work to understand the
$a$-dependence of $N_{i}(a)$, so we consider two distributions
($N_{i}(a)\propto a^{2}\delta a_{\rm lib}/a$, top panel) and ($N_{i}(a)\propto
a^{3.5}\delta a_{\rm lib}/a$, bottom panel), where the first is normalized to
the 3:2 estimated population and the latter is normalized to the 5:2 estimated
population. See text for further explanation. Scenario (2) results in
$N_{i}(a)\propto a$ which is normalized to the 3:2 estimated population
(middle panel).

##  5 Comparison with Observations

De-biased absolute population estimates of Neptune’s mean motion resonances
from Gladman et al. (2012),Volk et al. (2016), and Crompvoets et al. (2022)
provide a good opportunity to directly compare with populations affected by
the stochastic model discussed in this paper. These comparisons allow us to
comment on the possibility of these resonances’ current populations being set
during a planetesimal-driven migration era. We consider two methods for
filling the mean motion resonances with an initial number of objects
$N_{i}(a)$. The number distribution $N_{i}(a)$ is determined by the two
resonance capture models and respective normalization. For example, if
$N_{i}(a)\propto a$ and we normalize to the 3:2 estimated population,
$N_{\rm{est},3:2}$, such that the predicted population matches the estimated
population, then $N_{\rm est,3:2}=N_{i,3:2}P_{\rm{keep},3:2}$. Calculating the
remaining MMR populations is then a simple relation, which in this case is
$N_{i,5:3}=(N_{i,3:2}/a_{3:2})a_{5:3}$. After calculating the initial
population for each resonance, we apply the retention fraction calculated in
Figure 9 to get a population affected by noise during planetesimal-driven
migration, $N_{i}(a_{\rm res})P_{\rm keep,res}$

Scenario (1) appeals to a planetary upheaval that fills the phase space beyond
Neptune for a range of pericenters in the scattering region as described in
Balaji et al. (2023), motivated by the simulation results in Levison et al.
(2008). This pattern results from an early random walk in energy and angular
momentum space due to kicks at every perihelion passage with Neptune (or, in
principle, with another giant planet, depending on the details of the
scattering scenario). The time it takes for Neptune to change an object’s
angular momentum by an order of itself is proportional to the number of
objects at a given semi-major axis for a steady state. Yu et al. (2018)
reports that the number of scattering objects at a given semi-major axis from
30 to 100 au is $N_{i}(a)\propto a^{7/2}$ which they confirm through an
empirical fit to simulations from Kaib et al. (2011). 333Note that Yu et al.
(2018) give an incorrect explanation for the a-dependence in their Section 3.2
since the pericenter velocity at high eccentricity and fixed pericenter is not
proportional to $a_{p}$. On the other hand, Levison et al. (2006) approaches
the same problem by looking at planetesimals’ random walk in energy due to
perihelion passages with Neptune. They find the time it takes to change an
object’s energy by order of itself is proportional to $a^{-1/2}$, therefore
$N_{i}(a)\propto a^{-1/2}$. More work is needed to understand the numerical
distribution found by Kaib et al. (2011). A clue may come from the fact that
Yu et al. (2018) find that today’s scattering disk of TNOs and today’s
population of objects transiently-stuck in MMRs comprise a single population,
40% of which are “pausing" their random walk in a resonance at any given time.

In this scenario, we model the relative resonance populations as having
populations that follow the background scattering populations’
$a$-distribution. The scattering occurs before Neptune reaches its final
orbital location so that when this final location is reached, subsets of the
scattered objects happen to be in resonant orbits. The initial number of
objects in a resonance at semi-major axis $a$ is therefore $N_{i}(a)\propto
a^{k}\delta a_{\rm p,lib}/a$, where the number of objects at semi-major axis
$a$ in the scattered population scales as $a^{k}$. We find that normalizing to
the 3:2 MMR population, a power, $k$, from -1/2 to 2, and applying $P_{\rm
keep}$, provides a reasonable match to the observed populations, with the 2:1
MMR prediction the most discrepant. Figure 11 (top panel) illustrates the
upper end of this range ($k=2$), and $k=3.5$ provides a poor match. If we
instead normalize to the 5:2 MMR and do not worry about matching the 3:2 MMR
population due to complications resulting from secular resonances (Volk &
Malhotra, 2019), Yu et al. (2018)’s empirical fit, $N_{i}(a)\propto
a^{3.5}\delta a_{\rm p,lib}/a$ provides a more reasonable fit compared to the
same $k$ normalized to the 3:2 estimated population (green circles, bottom
panel). The 7:4 fraction is close to 0, but this resonance tends to capture
nearby objects easily for short periods of time (Lykawka & Mukai, 2005, 2007).

Scenario (2) represents efficient, preferential capture of objects into MMRs
due to the “smooth" average component of outward, planetesimal-driven
migration in a disk with a radial gradient of planetesimal number density. The
number of objects, $N_{i}(a)$ picked up into a resonance that ends this
migration with final semi-major axis $a_{f}$ can be estimated by $N(a)=2\pi
a\Delta a\Sigma(a)$ where the distance swept by the resonance $\Delta a$ is
given by the distance $\Delta a_{N}$ migrated by Neptune (the same for all
resonances) multiplied by $a/a_{N}$, which stays constant in time for each
resonance through Kepler’s third law and hence may be evaluated at the final
semi-major axes of both bodies. For a fiducial surface density of
planetesimals in the disk prior to resoance capture $\Sigma(a)\propto 1/a$, we
find $N(a)\propto a$. Figure 11 (middle panel) provides results for a smooth
migration as discussed in Nesvorný & Vokrouhlický (2016) with "grainiess" due
to 2000 Pluto sized planetesimals migrating over 10 million years. Other than
the 7:4 (which is likely to be repopulated by neighbors; Lykawka & Mukai 2007)
and 2:1 MMRs, the retention fractions reproduce observed population estimates.

The results in Figure 11 are not conclusive, but the fits are promising. Data
from the Rubin telescope’s LSST survey is expected to constrain resonance
populations for a much larger number of resonances. With these new data, it
will be possible to use patterns of resonance occupation to determine whether
stochasticity did—or equally interestingly whether it didn’t—play a
significant role in shaping the dynamical structure of the trans-Neptunian
region.

##  6 Discussion

The results presented in Section 5 are promising as we await unprecedented
data from LSST. Better observations will lead to improved constraints for the
relative populations and eccentricities of the objects in resonance. Such data
will merit revisiting this analysis and improving on our assumptions. We
discuss future work and additional considerations below.

###  6.1 Comparing the Imprint from Neptune’s Noisy Migration and Long Term
Evolution

We have calculated relative retention fractions for ten of Neptune’s exterior
mean motion resonances for the case of Neptune’s noisy planetesimal-driven
migration lasting 10 or 50 Myr. Between the end of such an epoch of migration
and now, additional small-scale perturbations imparted onto TNOs by the
planets cause TNOs in MMR to be lost from resonance. To determine whether the
pattern of resonance population we aim to match in this work (Figure 11)
indeed arises from noisy early migration, it will be necessary to understand
the pattern of retention across resonances produced by this long-term
dynamical sculpting. Tiscareno & Malhotra (2009) and Balaji et al. (2023)
studied the effect of small scale perturbations over billions of years between
the giant planets and objects in 3:2 and 2:1 MMRs. Tiscareno & Malhotra (2009)
quotes retention of 37% (27%) and 24% (15%) after 1 Gyr (4 Gyr) for the 3:2
and 2:1, respectively. For the 3:2 MMR long term evolution simulaitons from
Balaji et al. (2023), 29% (19%) were retained after 1 Gyr (4.5 Gyr).
444Tiscareno & Malhotra (2009)’s 3:2 initial population are tighter in
resonance than Balaji et al. (2023)’s, so a more fair comparison is possible
by finding the retention fraction relative to the 10 Myr population in Balaji
et al. (2023)’s simulations. If we count the resonant objects at 10 Myr as the
initial population, thus more tight in resonance, then 39% (25%) remain in
resonance after 1 Gyr (4.5 Gyr), providing a closer match to Tiscareno &
Malhotra (2009). The retention percentages in the text were not quoted in
Balaji et al. (2023)—we calculated them based on simulation data from Balaji
et al. (2023).  Both studies get comparable magnitudes for the retention
fractions, but relative retention fraction across a larger collection of MMRs
has not been systematically explored. Comparing the relative population
estimates across the resonances, rather than the percentages directly, will
provide insight into which method left a larger impact: long term evolution or
noisy migration. Future population estimates and studies similar to this work
and Tiscareno & Malhotra (2009); Balaji et al. (2023) extended to a larger
range of MMRs will help constrain which method had a larger role in shaping
the structure in the outer Solar System.

###  6.2 Retention Fractions in the Presence of a Ninth Planet

We do not include the presence of a hypothetical Planet Nine in our analysis
to calculate the retention fractions of objects in MMR between 39 and 75 au.
Clement & Sheppard (2021) considered the impact of an undetected planet around
$\sim 500$ au on distant $n:1$ resonances between $\sim$60-180 au (3:1-14:1).
In their study, the farthest resonances we study (3:1 and 4:1), are similarly
eroded in simulations with and without Planet Nine. In addition, TNOs with
perihelia less than around 40 au are expected to remain coupled with Neptune
(Clement & Sheppard, 2021). Therefore, we expect the relative retention rates
from this work to not be affected from an undetected distant planet. Future
work similar to Clement & Sheppard (2021) and that done here for even more
distant $n:1$ resonances that are significantly perturbed by a planet at 500
au might offer an interesting way to constrain the Planet Nine hypothesis.

###  6.3 Equivalent Analysis with Inclination-Type Resonances

The noisy migration framework in this paper was originally derived for the
first order 3:2 MMR with Neptune, where the resonance width (Equation (1)) and
$C_{\rm lib}$ are derived for low eccentricities to first order. In this
paper, we have shown that higher order, weaker resonances, provide strong
constraints for Neptune’s migration given their expected smaller population
and relative retention fraction. To second and third order, additional
inclination and mixed eccentricity/inclination weak resonances exist that we
have not considered here (See Chapters 6 and 8 in Murray & Dermott (1999)). An
analysis of these resonances analogous to that presented here would be
interesting. We leave this for future work, once LSST provides a larger amount
of data to constrain eccentricities, inclinations, and population estimates.

###  6.4 Neglecting Planetesimal-Planetesimal Interactions

We do not include gravitational interactions between planetsimals in the
simulations described in this work in an effort to reduce computation time. We
justify this choice here. Including the self-gravity of the disk affects the
dynamics of a planetary system in a number of ways: (1) a planet’s orbit will
precess due to the disk’s effect on the planet, (2) the planetesimals’
relative eccentricities will be both excited and damped through mutual
gravitational interactions, thus changing the dynamics of the system, and (3)
the total migration of a planet will be affected. For (1), it is not obvious
that the planet precessing would cause a significant change to the random walk
of the planet and thus loss of TNOs from resonance since the encounter times
between the planet and planetsimals driving the migration are shorter than the
precession time. For (2), to reproduce an accurate model for the excitement
and damping of planetesimals, we would require collision and gas damping
physics which mostly impacts the smallest planetesimals of which there are a
larger amount. To accurately describe the relative eccentricities we would
need hundreds of thousands of particles with a mass distribution, which is
currently computationally infeasible. We leave the relative eccentricities of
the planetesimals as a free parameter, shown by the low-eccentricity and high-
eccentricity simulation sets in 2.3 and do our analysis for the largest
planetsimals since they yield the most noise. For (3), we are only interested
in simulating the noisy part of the migration, not the overall behavior. The
timescale for the noisy migration component is shorter than the total
migration, so it is sufficient to consider these separately. Fan & Batygin
(2017) found that for a Nice model type of dynamical evolution, a self-
gravitating disk did not lead to significantly different final planetary
orbital behavior compared to a non-self-gravitatating disk. In addition, the
planetesimal evolution for both cases did not appear to be significantly
distinguishable after an instability was triggered. On the other hand, Quarles
& Kaib (2019) found that the self stirring of the disk has a greater impact on
the giant planets’ dynamics than previously thought. The effect of the self-
gravity of a disk on the dynamics of planetary systems is an active area of
research where improved computational capabilities such as GPUs are providing
an avenue to study this further (Grimm & Stadel, 2014; Grimm et al., 2022;
Kaib et al., 2024)

###  6.5 Surface Density Profile Effect on Resonance Retention Fraction

In calculating the relative retention fractions for each mean motion resonance
in this work, we assume a constant disk surface density parameterized by
$\mathcal{M}$ and by the fiducial semi-major axis of the disk, $a_{d}$, which
we hold constant (see Equation 6). We note that the semi-major axis of the
planet does not change substantially over the course of our simulations (see
Figure 4), so we are simulating a snapshot in time during the planet’s
migration. In reality, the total resonance loss will be produced by noise
properties integrated over the full migration.

Consider a disk surface density distribution that follows a negative power
law, $\Sigma_{m}\propto a^{-\beta}$, with $\beta>0$. From Equation (8ak),
$P_{\rm keep}$ declines exponentially when $\lambda_{n}T\propto
a^{1/2-\beta}T$ increases, where $T$ is the total migration timescale. Within
a given disk, an equivalent amount of time spent at small semi-major axis
produces more noise than at large semi-major axis as long as $\beta>1/2$. In
other words, for typically-assumed surface density power law distributions
with $1\lesssim\beta\lesssim 3/2$, $a^{1/2-\beta}T$ is largest at small $a$ as
long as $T$ remains roughly constant. As the overall disk surface density
decreases with time due to scattering by the planets, the surface density
encountered by Neptune as a function of semi-major axis likely declines more
rapidly than the initial disk surface density distribution (effectively
producing a larger $\beta$). Thus, the beginning of the epoch of smooth
migration (when chaotic planetary motion ceases and $T$ becomes long) is
likely to dominate loss of TNOs from resonance due to noise. The details of
this time evolution merit further work.

##  7 Conclusion

A planetesimal disk with fixed surface mass density $\Sigma$ generates more
stochasticity when composed of larger, fewer planetesimals. An epoch of
planetesimal-driven planetary migration, capturing objects into MMR and losing
a fraction due to noise, will leave clues behind regarding its dynamic past.
The analysis shown in this paper appeals to migration produced by
planetesimals over timescales of up to 100 Myr or to scenarios after a
gravitational upheaval, where late-stage planetesimal-driven migration is seen
in numerical simulations (e.g., Tsiganis et al., 2005) on timescales of up to
a few tens of Myr. The weak MMRs lose objects more easily, so we use them to
comment on the size distribution of the planetesimal disk such that the
weakest resonances with currently observed occupants do not lose their
populations due to the stochasticity of planetesimal-driven migration.

In this paper, we tested the validity of modeling the randomness in planetary
migration as Brownian motion by numerically confirming results from MC2006 and
deriving coefficients that were left out during order-of-magnitude
derivations. We analytically determine the relative retention fraction of
populations of TNOs captured by Neptune into a range of MMRs during migration.
Increasing the mass in large planetesimals contained in the disk or increasing
the migration timescale produces more loss across all resonances. However, the
proportion of this reduction is not equal across the resonances, with
particularly large losses occurring in the 5:3, 7:3, 7:4, and 9:5 resonances.
The smallest losses occur in the 2:1, 3:1, 4:1, and 3:2 resonances. The
‘strength’ of resonances can be ranked in terms of their loss times and
retention fractions (see Figures 8 and 9). We note that increasing the
eccentricity of the planetesimal disk increases the timescale needed to
produce the same amount of loss by a factor of $e^{2}/e_{H}^{2}$ for
planetesimal eccentricities larger than $e_{H}$. In addition, changing the
eccentricity of the resonance from $e_{\rm res,1}$ to $e_{\rm res,2}$ changes
the loss time by a factor of $\big{(}e_{\rm res,2}/e_{\rm
res,1}\big{)}^{J_{3}}$, where $J_{3}$ is the order of the resonance.

We apply our analysis to the case of the Solar System with a planetesimal disk
twice the mass of Neptune, corresponding to $\mathcal{M}=2$. We considered
size distribution end-member cases accounting for the mass of large
planetesimals with maximum radii $s_{m}=780$, 1100, and 2200 km. If part of
the 7:3 MMR’s population was obtained before or during the epoch of
planetesimal-driven migration, it cannot have a retention fraction of zero.
Given this condition, we find four conditions that must have been met for a
given maximum size of planetesimal to retain at least 30 % of TNOs in 7:3 MMR
for a disk with $e\leq e_{H}$:

  * •

the size distribution with $q=5.5$ extended to a maximum radius of 780 km
requires $\lesssim$80 Myr migration

  * •

the size distribution with $q=5.5$ extended to a maximum radius of 1100 km
requires $\lesssim$40 Myr migration

  * •

1333 Plutos require a 10 Myr migration

  * •

80 twice Pluto-size objects require $\lesssim$10 Myr migration

Additional scenarios may be extrapolated using these results using the facts
that for constant $\mathcal{M}T$ and constant
$(\mathcal{M}_{2}/\mathcal{M}_{1})/(m_{1}/m_{2})$, the same level of resonance
retention is produced. A more massive disk ($\mathcal{M}>2$) would produce
constraints for shorter timescales. Observational constraints on the
population and size distribution of bodies with radii $\gtrsim$780km are
limited by small number statistics, leaving uncertainty in the expected early
population and sizes of large bodies. In the limiting case that the large
bodies currently observed in the TNO region are the only large planetesimals
that ever existed in the outer disk, planetesimal-driven migration is
effectively smooth, with stochasticity producing no meaningful resonace loss.
The unprecedented data for TNOs from the Vera C. Rubin Observatory will
improve these parameters, and a similar analysis will cover the populations of
a significantly larger collection of distant resonances, producing stronger
constraints.

The values stated above that decrease the 7:3 population by 70% would decrease
the 5:2 by 25 % only. Observations show a large population in the 5:2
resonance—comparable to the number in the 3:2 MMR, with high eccentricities
ranging from approximately 0.2 to 0.6 (Gladman et al., 2012; Volk et al.,
2016; Malhotra et al., 2018). Presumably, if the 5:2 resonance was populated
before or during the epoch of planetesimal-driven migration, a large
proportion of the captured population must have been retained. This large
proportion of high eccentricity bodies could be due to detection biases since
beyond 46AU, resonant bodies with low eccentricities are difficult to detect.
As noted in Section 4.1, the 5:2 resonance is relatively weak for $e_{\rm
res}=0.25$ but substantially stronger for $e_{res}=0.4$, the typical value for
observed 5:2 objects measured (Figure 9). If future surveys sensitive to
fainter objects find many low-eccentricity objects in the 5:2, the low
eccentricity 5:2 will generate constraints similar to or better than those
from the 7:3. The Minor Planet Center555https://www.minorplanetcenter.net/data
shows that a few objects at low eccentricities have been found, but not
through a survey that allows a debiased population estimate.

Future investigations of migration scenarios, resonance capture efficiency,
and unprecedented data from the Vera C. Rubin Observatory will allow for
direct comparison between resonance retention fractions and improved
population estimates. We expect decent population estimates for mean motion
resonances up to 200 au; with this data, we will be able to investigate the
trends found in Figure 11. In particular, patterns in resonance population
that depend on the order of resonances as calculated in this work will be a
strong indicator that stochasticity played a significant role in sculpting the
trans-Neptunean, region. Conversely, the lack of such patterns would indicate
that this random process did not shape the populations we observe today.

## Acknowledgements

We thank Kat Volk for their comments on an earlier draft of this paper and the
anonymous reviewer for their comments which led to an improvement in this
paper and addition of Section 6. AHR thanks the LSSTC Data Science Fellowship
Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the
Brinson Foundation, and the Moore Foundation; her participation in the program
has benefited this work. This material is based upon work supported by the
National Science Foundation Graduate Research Fellowship Program under Grant
No. 1842400. Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the author(s) and do not necessarily
reflect the views of the National Science Foundation. AHR and RMC acknowledge
support from NASA grant number 80NSSC21K0376. We acknowledge use of the lux
supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315.

## Data Availability

The code for MERCURY that were used to produce the simulation data and plots
can be found at https://github.com/ahermosillo/stochastic-randommigration.

## References

  * Adams et al. (2014) Adams E. R., Gulbis A. A. S., Elliot J. L., Benecchi S. D., Buie M. W., Trilling D. E., Wasserman L. H., 2014, AJ, 148, 55
  * Alexandersen et al. (2016) Alexandersen M., Gladman B., Kavelaars J. J., Petit J.-M., Gwyn S. D. J., Shankman C. J., Pike R. E., 2016, AJ, 152, 111
  * Balaji et al. (2023) Balaji S., et al., 2023, MNRAS, 524, 3039
  * Bannister et al. (2018) Bannister M. T., et al., 2018, ApJS, 236, 18
  * Bernardinelli et al. (2022) Bernardinelli P. H., et al., 2022, ApJS, 258, 41
  * Brown (2008) Brown M. E., 2008, in Barucci M. A., Boehnhardt H., Cruikshank D. P., Morbidelli A., Dotson R., eds, , The Solar System Beyond Neptune. The University of Arizona Press, p. 335
  * Chambers (1999) Chambers J. E., 1999, MNRAS, 304, 793
  * Chen et al. (2019) Chen Y.-T., et al., 2019, AJ, 158, 214
  * Chiang et al. (2007) Chiang E., Lithwick Y., Murray-Clay R., Buie M., Grundy W., Holman M., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. p. 895 (arXiv:astro-ph/0601654)
  * Clement & Sheppard (2021) Clement M. S., Sheppard S. S., 2021, AJ, 162, 27
  * Crompvoets et al. (2022) Crompvoets B. L., et al., 2022, , 3, 113
  * Dawson & Murray-Clay (2012) Dawson R. I., Murray-Clay R., 2012, ApJ, 750, 43
  * Fan & Batygin (2017) Fan S., Batygin K., 2017, ApJ, 851, L37
  * Fernandez & Ip (1984) Fernandez J. A., Ip W. H., 1984, Icarus, 58, 109
  * Fraser et al. (2014) Fraser W. C., Brown M. E., Morbidelli A., Parker A., Batygin K., 2014, ApJ, 782, 100
  * Gallardo (2019) Gallardo T., 2019, Icarus, 317, 121
  * Gladman & Volk (2021) Gladman B., Volk K., 2021, ARA&A, 59, 203
  * Gladman et al. (2012) Gladman B., et al., 2012, AJ, 144, 23
  * Grimm & Stadel (2014) Grimm S. L., Stadel J. G., 2014, ApJ, 796, 23
  * Grimm et al. (2022) Grimm S. L., Stadel J. G., Brasser R., Meier M. M. M., Mordasini C., 2022, ApJ, 932, 124
  * Hahn & Malhotra (2005) Hahn J. M., Malhotra R., 2005, AJ, 130, 2392
  * Kaib & Sheppard (2016) Kaib N. A., Sheppard S. S., 2016, AJ, 152, 133
  * Kaib et al. (2011) Kaib N. A., Roškar R., Quinn T., 2011, Icarus, 215, 491
  * Kaib et al. (2024) Kaib N. A., Parsells A., Grimm S., Quarles B., Clement M. S., 2024, Icarus, 415, 116057
  * Kavelaars et al. (2021) Kavelaars J. J., Petit J.-M., Gladman B., Bannister M. T., Alexandersen M., Chen Y.-T., Gwyn S. D. J., Volk K., 2021, ApJ, 920, L28
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lan & Malhotra (2019) Lan L., Malhotra R., 2019, Celestial Mechanics and Dynamical Astronomy, 131, 39
  * Lawler et al. (2018) Lawler S. M., et al., 2018, AJ, 155, 197
  * Lawler et al. (2019) Lawler S. M., et al., 2019, AJ, 157, 253
  * Levison et al. (2006) Levison H. F., Duncan M. J., Dones L., Gladman B. J., 2006, Icarus, 184, 619
  * Levison et al. (2008) Levison H. F., Morbidelli A., Van Laerhoven C., Gomes R., Tsiganis K., 2008, Icarus, 196, 258
  * Luu & Jewitt (2002) Luu J. X., Jewitt D. C., 2002, ARA&A, 40, 63
  * Lykawka & Mukai (2005) Lykawka P. S., Mukai T., 2005, Planet. Space Sci., 53, 1175
  * Lykawka & Mukai (2007) Lykawka P. S., Mukai T., 2007, Icarus, 192, 238
  * Malhotra (1993) Malhotra R., 1993, Nature, 365, 819
  * Malhotra (1995) Malhotra R., 1995, AJ, 110, 420
  * Malhotra (2019) Malhotra R., 2019, Geoscience Letters, 6, 12
  * Malhotra et al. (2018) Malhotra R., Lan L., Volk K., Wang X., 2018, AJ, 156, 55
  * Morbidelli & Nesvorný (2020) Morbidelli A., Nesvorný D., 2020, in Prialnik D., Barucci M. A., Young L. A., eds, , The Trans-Neptunian Solar System. , pp 25–59
  * Murray & Dermott (1999) Murray C. D., Dermott S. F., 1999, Solar system dynamics. Cambridge University Press
  * Murray-Clay & Chiang (2006) Murray-Clay R. A., Chiang E. I., 2006, ApJ, 651, 1194
  * Nesvorný (2018) Nesvorný D., 2018, ARA&A, 56, 137
  * Nesvorný & Vokrouhlický (2016) Nesvorný D., Vokrouhlický D., 2016, ApJ, 825, 94
  * Pan & Sari (2005) Pan M., Sari R., 2005, Icarus, 173, 342
  * Petit et al. (2023) Petit J.-M., Gladman B., Kavelaars J. J., Bannister M. T., Alexandersen M., Volk K., Chen Y.-T., 2023, ApJ, 947, L4
  * Pike et al. (2017) Pike R. E., Lawler S., Brasser R., Shankman C. J., Alexandersen M., Kavelaars J. J., 2017, AJ, 153, 127
  * Quarles & Kaib (2019) Quarles B., Kaib N., 2019, AJ, 157, 67
  * Tiscareno & Malhotra (2009) Tiscareno M. S., Malhotra R., 2009, AJ, 138, 827
  * Tsiganis et al. (2005) Tsiganis K., Gomes R., Morbidelli A., Levison H. F., 2005, Nature, 435, 459
  * Volk & Malhotra (2019) Volk K., Malhotra R., 2019, AJ, 158, 64
  * Volk et al. (2016) Volk K., et al., 2016, AJ, 152, 23
  * Wolff et al. (2012) Wolff S., Dawson R. I., Murray-Clay R. A., 2012, ApJ, 746, 171
  * Yu et al. (2018) Yu T. Y. M., Murray-Clay R., Volk K., 2018, AJ, 156, 33
  * de Sousa et al. (2020) de Sousa R. R., Morbidelli A., Raymond S. N., Izidoro A., Gomes R., Vieira Neto E., 2020, Icarus, 339, 113605

##  Appendix A Average Migration

The average component of planetesimal-driven migration in the outer solar
system is determined by global asymmetries set up as planetesimals are
transeferred between the gravitational zones of influece of the various giant
planets (Fernandez & Ip, 1984). Since the model presented here only includes
one planet, this migration is not reproduced. We show here that the average
migration seen in our simulations may be derived using the same analytic
framework that we use to model stochasticity, verifying the self-consistency
of our results.

The mean component of planetesimal driven migration may be thought of as
resulting from a form of dynamical friction. Regardless of planetesimal
eccentricity, a population of planetesimals with semi-major axis deviation
$x=a-a_{p}$ shears past the planet with bulk relative velocity $v_{\rm
rel}\approx-(3/2)\Omega_{p}x$. Thus, the dynamical friction felt by the planet
from planetesimals with $x=x_{0}$ and those with $x=-x_{0}$ roughly cancels.
Migration results from the small difference felt between the torque from
planetesimals at positive and negative $x$. The difference in relative
velocity on either side of the planet ($x=\pm x_{0}$) due to disk curvature is
of order
$a_{p}\left(\frac{d^{2}\Omega}{da^{2}}\right)2x_{0}\sim(15/2)\Omega_{p}x(x/a_{p})$.
Power-law variations in other quantities on either side of the planet
similarly lead to differences of order the single-sided torque multiplied by
$x/a_{p}$. This second-order behavior is why the average distance migrated in
Figure 5 is not much larger than the standard deviation.

Considering planetesimals with a single value of $x$, the average (one-sided,
Sun-centered) torque experienced by the planet is $\tau_{1}\approx
a_{p}M_{p}\Delta v_{\parallel,p}\dot{\overline{N}}$, where
$\dot{\overline{N}}\approx(\Sigma_{m}/m)\Omega_{p}x^{2}$ is the average rate
at which the planetesimals shear past the planet at distance $x$. The change
in orbital velocity of the planet per encounter in the aximuthal direction is
$\Delta v_{\parallel,p}\sim-(m/M_{p})\Delta v_{\parallel}$ where $\Delta
v_{\parallel}$ and $\Delta v_{\perp}$ are the velocity changes experienced by
the planetesimal in the azimuthal and radial directions, respectively,
resulting from a single encounter. From the impulse approximation, $\Delta
v_{\perp}\sim GM_{p}/(xv_{\rm rel})$ and in the limit that tidal gravity from
the star does not change the energy of the planetesimal substantially over the
encounter,

$\Delta v_{\parallel}\approx\begin{cases}\Delta v_{\perp}^{2}/v_{\rm
rel}&\text{, $e<\Delta v_{\perp}/(\Omega_{p}a_{p})$,}\\\ v_{\rm ep}\Delta
v_{\perp}/v_{\rm rel}&\text{, $e>\Delta
v_{\perp}/(\Omega_{p}a_{p})$}\end{cases}$ (8ao)

where $v_{\rm ep}\sim e\Omega_{p}a_{p}$ is the particle’s epicyclic velocity
so that

$\Delta v_{\parallel}\sim\begin{cases}G^{2}M_{p}^{2}/(x^{2}v_{\rm
rel}^{3})&\text{, $e<e_{H}^{3}(a/x)^{2}$,}\\\ GM_{p}v_{\rm ep}/(xv_{\rm
rel}^{2})&\text{, $e>e_{H}^{3}(a/x)^{2}$.}\end{cases}$ (8ap)

Note that the criteria separating the two cases matches the criteria
separating the two cases in Equation (2.2).

Combining these expressions yields

$\tau_{1}\approx-\Sigma_{m}\Omega_{p}^{2}a_{p}^{4}\times\begin{cases}(M_{p}/M_{*})^{2}(a_{p}/x)^{3}&\text{,
$e<e_{H}^{3}(a/x)^{2}$,}\\\ e(M_{p}/M_{*})(a_{p}/x)&\text{,
$e>e_{H}^{3}(a/x)^{2}$.}\end{cases}$ (8aq)

The small-eccentricity case of Equation (8aq) is equivalent to 2D dynamical
friction resulting from a population that does not extend arbitrarily close to
the planet. This one-sided torque yields a rate of change of the planet’s
semi-major axis given by the approximation $\dot{a}_{p}\sim(\dot{L}/L)a_{p}$,
resulting in

$\frac{da_{p}}{dt}\sim\frac{\tau_{1}}{M_{p}a_{p}\Omega_{p}}\sim\dot{\overline{N}}\Delta
a_{\rm p,n}$ (8ar)

where $\Delta a_{\rm p,n}$ is given by Equation (2.2).

The net torque resulting from the difference between torques from particles at
$\pm x$ is $\tau\sim f\tau_{1}(x/a_{p})$, where $f$ is an order-unity
coefficient resulting from the combined powerlaw exponents of all quantities
changing in the disk on either side of the planet. For the first case in
equation 8aq that gives

$\displaystyle\tau\sim
f\Sigma_{m}\Omega_{p}^{2}a_{p}^{4}\times\begin{cases}(M_{p}/M_{*})^{2}(a_{p}/x)^{2}&\text{,
$e<e_{H}^{3}(a/x)^{2}$,}\\\ e(M_{p}/M_{*})&\text{,
$e>e_{H}^{3}(a/x)^{2}$.}\end{cases}$ (8as)

The torque from the lower-eccentricity case of Equation (8as) is familiar from
Type I migration, where $x$ is replaced by $H$ in a gas disk since gas
perturbations are limited to move at the sound speed, meaning that the closest
coherent Lindbland resonances pile up at $x\approx H$.

We can apply our normalization coefficents to this calculation in the
following way. The accumulated encounters influence the average distance
migrated, such that $\Delta a_{\rm avg}=f(x/a_{p})\overline{N}\Delta a_{p}$
where $\Delta a_{p}$ represents the change in semi-major axis due to a single
encounter and $\overline{N}=\dot{\overline{N}}T$ is the average number of
encounters. For all eccentricities, this average migration comes from the
average bulk relative velocity, due to shear, so that
$\dot{\overline{N}}=(\Sigma_{m}/m)\Omega_{p}x^{2}$ (folding any coefficient
into the overall coefficient we are calculating). The random walk distance
$\Delta a_{\rm p,T}=(\overline{N})^{1/2}\Delta a_{p}$ so that $\Delta a_{\rm
avg}=f(x/a_{p})\Delta a_{\rm p,T}^{2}/\Delta
a_{p}=\mathcal{C}_{5}f(x/a_{p})\dot{\overline{N}}T\Delta a_{\rm p,n}$ for non-
crossing cases 1 and 2, respectively. Plugging in Equations (8pa) and (8pb)
reveals the value of $\mathcal{C}_{5}$. For non-crossing case 1, the
coefficient is $\mathcal{C}_{5}=2\pi\mathcal{C}^{2}/(81\mathcal{C}_{1}^{2})$
and $\mathcal{C}_{5}=2\pi\mathcal{C}^{2}/(9\mathcal{C}_{2}^{2})$ for non
crossing case 2 (this coefficient applies for all high-eccentricity orbits,
including crossing orbits). Plugging in, we get

$\displaystyle\Delta a_{\rm avg}$ $\displaystyle=$ $\displaystyle
f\frac{x}{a_{p}}\frac{2\pi}{9}\mathcal{C}^{2}\frac{\Sigma_{m}}{m}\Omega_{p}x^{2}T$
(8at) $\displaystyle\times$
$\displaystyle\begin{cases}\frac{1}{9\mathcal{C}_{1}}\frac{mM_{\rm
p}}{M_{*}^{2}}\frac{a_{\rm p}^{6}}{x^{5}}&\text{, $e<e_{H}^{3}(a/x)^{2}$,}\\\
\frac{1}{\mathcal{C}_{2}}\frac{m}{M_{*}}\frac{a_{p}^{4}}{x^{3}}e&\text{,
$e>e_{H}^{3}(a/x)^{2}$.}\end{cases}$

Calculating $\Delta a_{\rm avg}$ at 1 million years for the simulations
plotted in Figure 5, using parameters provided for groups 3 and 5 in Table 1
gives $\Delta a_{\rm avg}=0.0039f$ au for the low-eccentricity simulaiton and
$\Delta a_{\rm avg}=0.0026f$ au for the high-eccentricity simulation. These
values match the simulation results $\Delta a_{\rm avg}=0.041$ au and 0.015
au, respectively, for $f=$10 and 6, respectively. These values are reasonable
(comparable to the 7.5 estimated at the beginning of this appendix) and
similar to one-another. We consider this reasonable verification that our
interpretation of our numerical results is self-consistent.