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

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 ≳ 30 Earth masses in today’s size distribution and 𝑇 mig ≳ 40Myr and (2) a top-heavy size distribution with ≳ 2000 Pluto-sized TNOs and 𝑇 mig ≳ 10Myr, where 𝑇 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.


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 wellcharacterized 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 ★ E-mail: arcelia@ucsc.eduSystem'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(Malhotra , 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 ≃ 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 ∼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. 1The 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 largescale 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 ∼ 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 ∼ 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.

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 ( 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.

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,  p how close they are to each other in semi-major axis, , and the encounter time, Δ.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, Δ, 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, Δ  , 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.≲ 102 yrs versus ∼ 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,  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  res , where C lib is a constant that depends on the resonance,  3 is the order of resonance,  p and  p are the semimajor axis and mass of the planet, respectively and  * is the mass of the Sun (for an expository derivation, see the textbook by Murray & Dermott 1999). 2 The libration width,  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 its Hill eccentricity is and its Hill velocity is where Ω p = (  * /3  ) 1/2 is the orbital angular velocity of the planet.The random velocity of a planetesimal, , with an eccentricity of  is where Ω is the planetesimal's orbital angular velocity and  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, M, which serves to parameterize the surface density of planetesimals in the disk, Σ  , such that when integrated over the surface, the total mass of the disk is M  p distributed among planetesimals of mass .The surface density is parameterized for one mass  because a single mass (typically the largest planetesimals) generally produce the most noise (MC2006).The mass contained in the disk is uniformly distributed between  d /2 to 3 d /2, where  d is the semi major axis at the middle of the disk, and so We keep   = 26au constant throughout this paper to keep a consistent translation between M and Σ  .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   ) since that will dictate the encounter rate.We use Equation ( 6) to intuitively relate this local Σ  to the amount of solid mass needed to form the planet.We define  ≡  −  p where R is this distance normalized by  H such that Figure 1.Numerical simulations (circles) validate analytical Equation (8a,8b) using C 1 = 6.27 and C 2 = 2.5 (lines) for the change in a planet's semimajor axis,  p , during a non-crossing encounter.We fit the lines to the simulations with initial /R  ≤ 25 since beyond that, the approximation that  =  −  p ≪ 1 begins to lose validity.We perform three sets of simulations, all with planetary mass  p = 10 −5  * and semi-major axis  p = 30AU.A planetesimal of mass  and eccentricity  starts with semimajor axis  p + .

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 || > .For these values of , the relative speed of the planetesimal is dominated by the Keplerian shear velocity ∼Ω p , such that the encounter time is approximately a dynamical time, Ω −1 p .The expected change in  p for a non-crossing orbit is if where C 1 and C 2 are order-unity numerical factors that were not determined during derivation in MC2006.
To find Δ 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 ( 8), as seen in Figure  10) and ( 12) (lines) for the change of a planet's semi-major axis during a crossing encounter when setting C 3 = 2.5.For both simulations, planetesimal semi-major axis, eccentricity, and mass are  = 25AU,  = 0.4, and  = 10 −10  * , respectively.The planet starts at semi-major axis  p = 28 au, with eccentricity  = 0 and masses  p = 10 −5  * for simulation set C(a) and  p = 5 × 10 −6  * for C(b).At impact parameters,  less than  min (Eq.11), Δ p is constant, as predicted by the maximum values given by Equation (10).Beyond that, Δ p falls off inversely proportionally to , as in Equation ( 12).8a) and (8b) with the linear least-squares method in log space.We perform the fit to simulation outputs with planetesimals initialized at /R  ≤ 25 since beyond that the approximation that  =  −  p ≪ 1 begins to break and the slope for Δ p begins to flatten out.The three simulations had different best-fit parameters of  1 = 6.27,  2 = 2.17, and  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 C 2 = 2.5.In a crossing orbit, || < , such that  > Ω so the relative velocity of the planetesimal is dominated by its epicylic velocity and its encounter time is less than Ω −1 .The impact parameter  can be rather different from ||, and  lies in the range where   is the impact parameter below which a collision occurs.This wide range can cause a wide spectrum of results for Δ p,c (where Δ p,c is the change of  p due to a crossing encounter) so MC2006 limit themselves to predict the maximum value of |Δ 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 generated when the impact parameter satisfies Once  does not satisfy Equation ( 11) we expect |Δ p,c | to decrease from the value in Equation ( 10), following the relationship As before, we find Δ 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   = 28au 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  =  p , while the planetesimal's orbit satisfies  = (1 −  2 )/(1 +  cos  ), so that an encounter occurs at angle  enc = cos −1 ((1 −  2 ) −   )/(  ) .We begin the simulation at time  = 2.5 years before the collision, when the mean anomalies of the objects are smaller by ΔM p = Ω p  for the planet and ΔM = Ω for the planetesimal.Since it is on a circular orbit, the planet's initial mean anomaly is M p =  enc − ΔM p .The planetesimal's initial mean anomaly is M = M enc − ΔM, where we calculate the mean anomaly at encounter, M enc , from  enc using Kepler's equation M enc =  −  sin  for the eccentric anomaly  = tan(  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, .
We run two sets of simulations with various impact parameters, , to verify Equations ( 10) and ( 12 10) but not its value at large impact parameters (Equation 12).Our simulations validate the equations for crossing orbits with C 3 = 2.5 (Figure 2 dashed lines).

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  p (i.e. its random walk in semi-major axis) is well modeled by Brownian motion.
MC2006 derive the expected random walk, Δ p,T , of a planet in a disk of planetesimals given by where  is the mean encounter rate,  is the migration duration, and Δ 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,  <   , the speed of the planetesimals relative to the planet is determined by Keplerian shear, and impact parameter  ∼  resulting in a mean encounter rate, Planetesimals in crossing orbits encounter the planet with super-Hill velocities ( >   ) at impact parameters  ≲    / 2 resulting in a mean encounter rate, where  ∼   Ω  .
We summarize Equation 13for the non-crossing and crossing cases as where C is a constant that encapsulates all factors of order unity that arose during derivations, and  is the eccentricity of the planetesimals.The non-crossing orbit scenario takes two forms (Equations 16a and 16b) since Δ p will vary depending on how |Δ| compares to the pre-encounter .Equation (16c) corresponds to crossing orbits.Note that Equation (16c) does not depend on 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 C (where C > 1 but not ≫ 1) and verify the total random walk of the planet, Δ 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,  H , since disks with  >  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  p = 10 −6  * , so that  H = 0.007.We explore the disk's low and high eccentricity regimes with simulation sets D(a) and D(b) where the initial  = 0.001 and  = 0.05, respectively.The planetesimals are all of equal mass  = 1.8 × 10 −10  * .Their semi-major axes  are uniformly distributed between ∼20 (21.5) au and ∼30 (28.5) au for D(a) (D(b)), inclinations between 0 and the value of  in radians, and remaining orbital angles are randomly generated between 0, 2.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.
For Equation ( 16) to be valid, the system must reach a pseudosteady 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.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 ∼ 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 Δ p,T given by Equation ( 16) as the standard   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.
deviation of the planet's semi-major axis across the 150 simulations,   p,T .
Because Equation ( 16) depends on several orbital parameters of the disk and properties of the planet, we must disentangle which subpopulation of simulated planetesimals causes the most randomness.To do that, we choose six regions of  −  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 R, , and M are shown in Table 1.
We calculate the surface density of both disks for each group centered at ± au away from the planet, where  = R   .For groups 1−3 and 6, the groups have radial width 2  and eccentricities within the limits expressed in Equations ( 16a) and (16b).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 (16c) does not depend on 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+||) au.The non-dimensional surface density parameter, M, is then determined by Equation ( 6) summing 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 Δ p,T : eccentricity of the planetesimals or proximity to the planet.
Inputting the values in Table 1 and Figure 3, we calculate Equation ( 16) for each group as a function of time and compare to the amount of migration from simulations, as shown in Figure 7.The simulationderived 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 C = 3.5.These groups have the largest impact on the planet's random walk (i.e. have largest values of Δ 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 C to 3.5 for Equation ( 16) provides a good fit to   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 Δ p ∝  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 ( 16) 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.

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 charac-  6 and Table 1 (solid lines), determined by Equation ( 16).The orbital parameters and local surface density from Groups 3 (dark purple) and 5 (light green) produce the largest Δ 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 Δ p,T to smallest and inline labels to clearly compare between low and high eccentricity groups.
teristic 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 M = M LG , which represents the surface density of the disk in planetesimals with radius equal to 780 km and  = 2g/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 M = 2, which corresponds to a total of ≳ 30 M ⊕ .Planetesimal disk masses between 20 to 35  ⊕ are typical for models of Neptune's migration (Luu & Jewitt 2002;Nesvorný 2018).Within this framework, M LG = 2, where M 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 M LG imply that only a fraction M LG /M of the total mass is in these large bodies.Thus, the value of M LG for a given time is directly related to the size distribution of the planetesimal disk.In addition to M  , we do this analysis for M  and M 2 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 M  , M  , and M 2 in Section 4.

Ranking of loss times
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 (Δ p,T <  p,lib /2).We call this timescale, the resonance loss timescale,  loss , which can also be interpreted as the migration duration for which a planet will random walk a total of more than  p,lib /2, thus causing loss from resonances for given parameters of the disk.Using Equation (1) and Equation (16a), appropriate for  <  H /R 2 (other regimes produce analogous results), the resonance loss timescale is Figure 8 presents  loss for our collection of resonances, evaluated using Equation ( 17) with R = 1,  = 0.01,   =   = 26.6 au,  = 0.3 PL ,  p =  N , and  * =  ⊙ , where  PL ,  N and  ⊙ are the masses of Pluto, Neptune and the Sun, respectively, so that  H ∼ 0.69 au and  H ∼ 0.026.We set C = 3.5, the value determined in Section 2.3, and plug in typical eccentricities of TNOs in MMR for  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 2016).We determine  res for the resonances missing observed eccentricities by assuming the resonances possibly have a similar pericenter distribution,  free , motivated by the observed eccentricities increasing as a function of semi-major axis.Similarly stable libration zone studies in (, ) 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  res = 1 −  free / res produced a best-fit value  free =∼ 34.5 au, resulting in appropriate  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  loss as a function of M = M 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  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 ( 17) are the same across the disk except those relating to the resonance.The strength of the resonance therefore depends highly on  lib ,  res and  3 .Since  loss ∝ ( p,lib ) 2 ∝   3 res , a small change in  res for higher order resonances ( 3 > 1) will significantly affect the ranking of the resonances.
The observed eccentricities  res used to calculate  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,  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  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  3 where the order of the resonance,  3 , is 3, for a given M LG ).Though the resonance width decreases sharply at low eccentricities, a loweccentricity 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),  res = 1−/ res for a resonance at semi-major axis  res hosting an object with pericenter distance .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  2), as a function of M = M LG (see Equation 6).For this calculation, R = 1,  = 0.01,   =   = 26.6 au,  = 0.3 PL ,  p =  N , and times change by a constant factor given a fixed set of disk conditions (this factor is simply a combination of R and  H /).
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  res during migration and transiently stuck populations (Yu et al. 2018) should also be considered.

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  of migration is where   = () 2 /(2 2 p,lib ).The diffusion coefficient describing the cumulative random-walk change in semimajor axis, Δ p,T due to many encounters over a total time,  is As in all diffusive processes, the diffusion coefficient relates this large-scale behavior to the individual random steps in semi-major axis Δ p resulting from single encounters spaced in time by an average Δ ≡ 1/ , where  is the mean encounter rate given by Equations ( 14) and ( 15) for non-crossing and crossing orbits.The diffusion coefficient for the small-scale case is then  ∼ C 4 (Δ p ) 2 /Δ.
For our results to be self-consistent, we expect the coefficient C in Equation ( 16) 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 Δ p,T = Δ p √︁ C 4 (/Δ) equal to Equation ( 16) for each regime and solving for C 4 .We find that for cases (a), (b), and (c), respectively.Since we found that C 1 is 2.5 times larger than C 2 and C 3 (see Section 2.2), these results verify that 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  keep across the resonances we study for the case  ∼  H , corresponding to a diffusion coefficient calculated with Equation (16b).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  res = 0.25 (bottom panel) produces different strengths.In particular, the 7:4 MMR becomes stronger than the other weakest resonances for lower  res .The 7:3 MMR and in particular the 5:2 MMR produce significantly lower retention fractions for smaller  res , highlighting the correlation with  res .If Neptune traveled in a disk with  ≤  H , we expect the same retention fraction behavior since Equations ( 16a) and (16b) are the same when  =  H and R = 1, thus producing the same amount of noise.For the case where Neptune migrates across a hot disk (i.e. >   ), the timescale to lose the same fraction of objects as the first two cases increases by a factor of  2 / 2  .Increasing the number of large objects in the disk (e.g., increasing M 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 M  .

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  2 where  is the number of planetesimals and  is the mass of an individual planetesimal (c.f.Equation ( 16), noting that  ∝ M/).It turns out that larger-radius planetesimals, , produce the highest  2 .We summarize why that's the case here.The Kuiper Belt size distribution is described by a power law, / ∝  − , and integrating that gives  () ∝  1− .Assuming spherical, uniform-density TNOs, then  () 2 ∝  7− so that for  < 7 (typical of Kuiper Belt size distributions)  2 is larger for larger .Given that few TNOs have direct size estimates, their measured absolute () magnitude distribution is typically discussed instead, where / ∝ 10  .The size distribution exponent  is related to  by  = 5 + 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 ≲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 endmember 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 Increasing the number of large planetesimals in the disk produces more noise, causing a smaller retention fraction (each resonance gradient gets lighter for higher M 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.
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.

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  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  distribution for the dynamically hot population from Lawler et al. ( 2018) has a steep slope ( = 0.9) for large objects (low-) and a break magnitude   = 8.3, corresponding to a radius of 75 km (assuming an albedo,  = 0.04).Most of the mass resides in planetesimals with sizes near the break.We calculate the parameterized mass in large objects, M LG , for this size distribution with a simple approximation: where M = 2 (appropriate for ∼ 30 M ⊕ in the disk),   is the radius at the break of the broken power law (75 km) and   is the maximum radius covered in the distribution (780 km), and  is 5.5 for sizes ranging from the break to the maximum.As a result, M 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 ∼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 ( 20) and ( 21) for a planetesimal mass  =  pluto and   = 1100 km.Considering the same size distribution again, now M 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.

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, M  /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  = 780 km and  = 2g/cm 2 , then M 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 − ∼ 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 ∼ 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, 780km 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  ,1 to  ,2 will be (M 2 /M 1 )/( 1 / 2 ) where the masses  are proportional to  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  = 60 Myr, 3000 Plutos with  = 10 Myr, and 200 twice-Plutos with  = 10 Myr.Requiring 25-30 % retention for the 7:3 MMR produces stronger constraints: 2000 780-km object with  = 60 Myr, 1333 Plutos for  = 10 Myr, 80 twice-Plutos for  = 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  keep (   ) is proportional to M.

COMPARISON WITH OBSERVATIONS
De-biased absolute population estimates of Neptune's mean motion resonances from Gladman et al. (2012), Volk et al. (2016), andCrompvoets 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 (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 -dependence of   (), so we consider two distributions (  () ∝  2  lib /, top panel) and (  () ∝  3.5  lib /, 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   () ∝  which is normalized to the 3:2 estimated population (middle panel).
().The number distribution   () is determined by the two resonance capture models and respective normalization.For example, if   () ∝  and we normalize to the 3:2 estimated population,  est,3:2 , such that the predicted population matches the estimated population, then  est,3:2 =  ,3:2  keep,3:2 .Calculating the remaining MMR populations is then a simple relation, which in this case is  ,5:3 = ( ,3:2 / 3:2 ) 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,   ( res ) 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   () ∝  7/2 which they confirm through an empirical fit to simulations from Kaib et al. (2011). 3On 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  −1/2 , therefore   () ∝  −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' -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  is therefore   () ∝    p,lib /, where the number of objects at semi-major axis  in the scattered population scales as   .We find that normalizing to the 3:2 MMR population, a power, , from -1/2 to 2, and applying  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 ( = 2), and  = 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,   () ∝ 3.5  p,lib / provides a more reasonable fit compared to the same  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,   () picked up into a resonance that ends this migration with final semi-major axis   can be estimated by  () = 2ΔΣ() where the distance swept by the resonance Δ is given by the distance Δ  migrated by Neptune (the same for all resonances) multiplied by /  , 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 Σ() ∝ 1/, we find  () ∝ . 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 res-onances.With these new data, it will be possible to use patterns of resonance occupation to determine whether stochasticity didor equally interestingly whether it didn't-play a significant role in shaping the dynamical structure of the trans-Neptunian region.

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.

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). 4 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.

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 ∼ 500 au on distant  : 1 resonances between ∼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  : 1 resonances that are significantly perturbed by a planet at 500 au might offer an interesting way to constrain the Planet Nine hypothesis.

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  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.

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)

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 M and by the fiducial semi-major axis of the disk,   , which we hold constant (see Equation 6).We note that the semimajor 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, Σ  ∝  − , with  > 0. From Equation ( 18),  keep declines exponentially when    ∝  1/2−  increases, where  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  > 1/2.In other words, for typically-assumed surface density power law distributions with 1 ≲  ≲ 3/2,  1/2−  is largest at small  as long as  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 ).Thus, the beginning of the epoch of smooth migration (when chaotic planetary motion ceases and  becomes long) is likely to dominate loss of TNOs from resonance due to noise.The details of this time evolution merit further work.

CONCLUSION
A planetesimal disk with fixed surface mass density Σ 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 latestage 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  2 / 2  for planetesimal eccentricities larger than   .In addition, changing the eccentricity of the resonance from  res,1 to  res,2 changes the loss time by a factor of  res,2 / res,1  3 , where  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 M = 2.We considered size distribution end-member cases accounting for the mass of large planetesimals with maximum radii   = 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  ≤   : Figure1.Numerical simulations (circles) validate analytical Equation (8a,8b) using C 1 = 6.27 and C 2 = 2.5 (lines) for the change in a planet's semimajor axis,  p , during a non-crossing encounter.We fit the lines to the simulations with initial /R  ≤ 25 since beyond that, the approximation that  =  −  p ≪ 1 begins to lose validity.We perform three sets of simulations, all with planetary mass  p = 10 −5  * and semi-major axis  p = 30AU.A planetesimal of mass  and eccentricity  starts with semimajor axis  p + .NC(a) (top panel) verifies Equation (8a) with  = 0 and  = 5 × 10 −10  * .NC(b) and NC(c) (bottom panel) verify Equation (8b) with  = 0.01 for masses  = 5 × 10 −10  * (b) and  = 10 −12  * (c).

Figure 3 .
Figure 3.Typical  −  snapshot after 10 6 years for one of the 150 loweccentricity (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×10 −10  * )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  p = 10 −6  * 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 .Figure 5 .
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.

Figure 6 .
Figure6.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 Figure3to visualize the orbital differences between six regions in the disk, to probe which region produces the higher Δ 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 Table1 and Figure 7.The planetesimals in group 3 have eccentricities ≲   /R 2 (blue solid line) so Equation (16a) determines their impact on the planet's noisiness.Planetesimals in groups 1, 2, and 6 have eccentricities ≳   / 2 but ≲ R  (dark pink dashed line) so Equation (16b) determines their impact.Planetesimals in groups 4 and 5 occupy crossing orbits with eccentricities ≳ R  , and their impact on the planet's random walk is determined by Equation (16c).The Hill eccentricity of the planet is shown for reference (blue dashed line).

Figure 7 .
Figure 7.The numerical, expected planet's random walk, determined by the standard deviation, , 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 Figure6and Table1(solid lines), determined by Equation (16).The orbital parameters and local surface density from Groups 3 (dark purple) and 5 (light green) produce the largest Δ 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 Δ p,T to smallest and inline labels to clearly compare between low and high eccentricity groups.
Figure8presents  loss for our collection of resonances, evaluated using Equation (17) with R = 1,  = 0.01,   =   = 26.6 au,  = 0.3 PL ,  p =  N , and  * =  ⊙ , where  PL ,  N and  ⊙ are the masses of Pluto, Neptune and the Sun, respectively, so that  H ∼ 0.69 au and  H ∼ 0.026.We set C = 3.5, the value determined in Section 2.3, and plug in typical eccentricities of TNOs in MMR for  res (see Table2).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 Table2 basedon 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  res for the resonances missing observed eccentricities by assuming the resonances possibly have a similar pericenter distribution,  free , motivated by the observed eccentricities increasing as a function of semi-major axis.Similarly stable libration zone studies in (, ) 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 Table2to  res = 1 −  free / res produced a best-fit value  free =∼ 34.5 au, resulting in appropriate  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  loss as a function of M = M LG .The MMRs requiring the longest time to lose objects (2:1, 3:1, 4:1, 3:2; Figure8) 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 Figure8matches the order of  p,lib for the values in Table2, which is in agreement with the width

Figure 8 .
Figure 8. Characteristic resonance loss time,  loss , as defined by Equation (17) for the exterior mean motion resonances of Neptune that we study (see Table2), as a function of M = M LG (see Equation6).For this calculation, R = 1,  = 0.01,   =   = 26.6 au,  = 0.3 PL ,  p =  N , and

Figure 9 .
Figure 9. Retention probabilities predicted by Equation 18 for Neptune's exterior mean motion resonances (between panels) for the recorded  res (top panel) and  res = 0.25 (bottom panel).For this calculation,  = 10 Myr R = 1,  =   = 0.026,   =   = 26.6 au,  =  PL ,  p =  N , and  * =  ⊙ , and M = 0.22 which corresponds to 2000 Plutos in the planetesimal disk.Varying  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 .
Figure10.Retention fractions spanning 0 (white) to 1 (dark green) for Neptune's exterior mean motion resonances (between panels) as a function of M  for migration timescales 10 million years (top panel) and 50 million years (bottom panel). keep is calculated with the parameters described in Figure9(top panel) for a range of surface densities M  representing a planetesimal disk with the largest objects having  = 0.3 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 M 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
Tiscareno & Malhotra (2009)'s 3:2 initial population are tighter in resonance thanBalaji et al. (2023)'s, so a more fair comparison is possible by finding the retention fraction relative to the 10 Myr population inBalaji 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 toTiscareno & Malhotra (2009).The retention percentages in the text were not quoted in Balaji et al. (2023)-we calculated them based on simulation data fromBalaji et al. (2023).

Table
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  = 24.9585au (24.9847 au) and standard deviation  = 0.0253 au (0.0103 au).This demonstrates that the planet's migration is well modeled as a diffusion process.

Table 2 .
Murray & Dermott (1999)up to forth order lying between the 3:2 and 4:1 resonances (denoted by  : ), listed in order of ascending semi major axis,  res .Values of C lib are provided for reference (seeMurray & Dermott (1999)).For resonances with a citation listed in column five, we quote the average value of eccentricities of observed resonant TNOs,  res , which may be affected by selection effects.The remaining resonances' eccentricities were determined by fitting the data to  res = 1 −  free / res ,  free is a free parameter we fit.The best-fit pericenter distance,  free , value was then used to set a typical  res .
where  PL ,  N and  ⊙ are the masses of Pluto, Neptune and the Sun, respectively, so that  H ∼ 0.69 au and  H ∼ 0.026.The eccentricity of particles in resonance,  res , is set to the typical value observed today for each resonance or a value determined by fitting the data to  res = 1 −  res / res for those without measurements (see Table2 and text).