# Uranus’s influence on Neptune’s exterior mean motion resonances

Severance Graham Severance Graham<EMAIL_ADDRESS>Lunar and Planetary
Laboratory, The University of Arizona, 1629 E University Blvd, Tucson, AZ
85721  Kathryn Volk Planetary Science Institute, 1700 East Fort Lowell, Suite
106, Tucson, AZ 85719, USA  Lunar and Planetary Laboratory, The University of
Arizona, 1629 E University Blvd, Tucson, AZ 85721

###### Abstract

Neptune’s external mean motion resonances play an important role in sculpting
the observed population of transneptunian objects (TNOs). The population of
scattering TNOs are known to ‘stick’ to Neptune’s resonances while evolving in
semimajor axis ($a$), though simulations show that resonance sticking is less
prevalent at $a\gtrsim 200-250$ au. Here we present an extensive numerical
exploration of the strengths of Neptune’s resonances for scattering TNOs with
perihelion distances $q=33$ au. We show that the drop-off in resonance
sticking for the large $a$ scattering TNOs is not a generic feature of
scattering dynamics, but can instead be attributed to the specific
configuration of Neptune and Uranus in our solar system. In simulations with
just Uranus removed from the giant planet system, Neptune’s resonances are
strong in the scattering population out to at least $\sim 300$ au. Uranus and
Neptune are near a 2:1 period ratio, and the variations in Neptune’s orbit
resulting from this near resonance are responsible for destabilizing Neptune’s
resonances for high-$e$ TNO orbits beyond the $\sim 20$:1 resonance at
$a\approx 220$. Direct interactions between Uranus and the scattering
population are responsible for slightly weakening Neptune’s closer-in
resonances. In simulations where Neptune and Uranus are placed in their mutual
2:1 resonance, we see almost no stable libration of scattering particles in
Neptune’s external resonances. Our results have important implications for how
the strengths of Neptune’s distant resonances varied during the epoch of
planet migration when the Neptune-Uranus period ratio was evolving. These
strength variations likely affected the distant resonant and detached TNO
populations.

Orbital resonances, Trans-Neptunian objects, Resonant Kuiper belt objects,
Scattered disk objects, celestial mechanics

††facilities: ADS††software: rebound

##  1 Introduction

Many known trans-neptunian objects (TNOs) are observed to orbit within
Neptune’s external mean-motion resonances (see, e.g., review by Gladman & Volk
2021; we refer to mean-motion resonances simply as resonances hereafter). In
particular, the observed set of TNOs on orbits with semimajor axes beyond the
classical belt region, $a\gtrsim 50$ au, imply the existence of large
intrinsic populations in many of Neptune’s distant resonances (see Crompvoets
et al. 2022 for recent population estimates). TNOs in Neptune’s resonances can
either be objects captured onto resonant orbits that are stable for a gigayear
or longer during the epoch of giant planet migration, or they can be objects
that have temporarily ‘stuck’ to resonances while evolving in semimajor axis
due to perturbations from the giant planets (see, e.g., review by Malhotra
2019). In either case, to understand the significance of the observed resonant
TNO populations, we need improved models of the dynamical extent of Neptune’s
resonances in the distant solar system. Volk & Malhotra (2022) recently showed
that, in the current solar system, Neptune’s resonances remain strong out to
surprisingly large semimajor axes of several hundred au in the high-perihelion
population of distant TNOs. Here we focus on the role of distant resonances in
the lower-perihelion scattering TNO population where the phenomenon of
resonance sticking is known to occur (e.g. Gallardo, 2006; Lykawka & Mukai,
2007; Yu et al., 2018).

Scattering TNOs are typically defined as objects whose semimajor axes
experience significant changes over short dynamical timescales due to
relatively close approaches to Neptune at perihelion; the widely used Gladman
et al. (2008) scheme defines them as objects that experience $\Delta a>1.5$ au
over a 10 Myr integration of their present-day orbit. Objects with unchanging,
non-resonant semimajor axes beyond the classical belt region (the classical
belt being typically defined as sunward of Neptune’s 2:1 resonance) belong to
the dynamically ‘detached’ TNO population. The perihelion distance ($q$)
boundary between detached and scattering TNOs is complex and semimajor axis
dependent (see, e.g., discussions in Saillenfest 2020 and Batygin et al.
2021), but non-resonant TNOs with $q\lesssim 37$ au are quite likely to
experience significant changes in orbital energy, and therefore changes in
$a$, as they evolve over time. This mobility in $a$ allows the particles to
encounter different resonances, and they often experience repeated bouts of
temporary libration in these resonances in between scattering events. Yu et
al. (2018) examined the detailed sticking behavior of scattering TNOs with
$a<100$ au and found that at any given time, $\sim 40\%$ of the scattering
population is temporarily librating in one of Neptune’s many resonances; these
temporary resonance sticks spanned anywhere from a few libration cycles ($\sim
10^{5}$ years) to a gigayear. We note that a reasonable lower boundary in $q$
for the population of scattering TNOs most likely to experience significant
resonance sticking is $q\approx 33$ au; scattering TNOs with lower $q$
experience very strong encounters with Neptune (and the other giant planets
for very low-$q$ TNOs) that shorten their dynamical lifetimes compared to the
$q\geq 33$ au TNOs (e.g., Tiscareno & Malhotra, 2003), reducing the time
available for resonance sticking.

Lykawka & Mukai (2007) performed numerical simulations of the scattering
population with the current giant planets and examined the frequency and
locations of resonance sticking events. They found that sticking was unlikely
for scattering particles with semimajor axes $a\gtrsim 250$ au. Presumably
this drop-off in resonance sticking is due to the resonances becoming weak at
these semimajor axes compared to the scattering perturbations. However,
analytical models (e.g. Pan & Sari, 2004; Batygin et al., 2021) of Neptune’s
resonances predict that the onset of resonance overlap, and thus the
destruction of stable libration zones, for perihelion distances similar to
those of scattering TNOs should occur at semimajor axes greater than the
$a\sim 250$ au drop-off seen in the simulations of resonance sticking. This
mismatch between the scattering simulation results and the analytical
resonance model expectations partly is due to the inherent limitations of
analytical approaches in terms of the number of resonances modeled; typically
analytical models only consider a subset of the lowest-order resonances, which
in the case of the scattering population are the N:1 and N:2 resonances.
Recent work by Hadden & Tremaine (2024) demonstrates how numerical mapping
approaches, which do not suffer from this limitation, can be used to model
resonant and scattering dynamics more accurately in the one-planet case.
However, we will show in this work that there is an additional reason for the
mismatch between models of the resonances and the findings of full scattering
simulations: dynamical interactions between Uranus and Neptune strongly
influence the strength of Neptune’s resonances and drive the $a\sim 250$ au
drop-off in resonance sticking in the present-day scattering population. In
the low-perihelion regime of the scattering TNOs, single-planet models
(whether analytical or numerical) fail to capture this dynamical influence
that weakens Neptune’s resonances. Resonance sticking is a very important
dynamical process in shaping the present-day scattering TNO population, but
variations in the strength of Neptune’s resonances during the epoch of planet
migration would also leave an imprint on the resulting resonant and detached
TNOs. It is thus important to have more complete models of Neptune’s external
resonances and a better understanding of what affects their strengths.

The rest of the paper is organized as follows. In Section 2 we describe how we
numerically model Neptune’s resonances, show how we identified Uranus’s
important role in affecting them, and describe how we quantify resonance
strengths from our simulations. In Section 3 we present our exploration of how
the Neptune-Uranus period ratio affects Neptune’s external resonances, and we
isolate the specific dynamical effects responsible for destabilizing the
distant resonances in today’s scattering population. Section 4 discusses the
implications of our findings for the scattering population and for using the
distant resonant and detached populations to constrain the migration history
of the giant planets. We summarize our results in Section 5.

##  2 Simulating Neptune’s resonances

A conceptually simple way to explore the strength and boundaries of Neptune’s
resonances is to produce Poincarè maps, also called surfaces of section. In a
simplified system with just the Sun and Neptune present, these maps of any
particular resonance can be produced by simulating test particles with fixed
perihelion distances over a small range of semimajor axis, $a$, near the exact
resonant value and recording the evolution of these particles’ orbits and the
angle $\psi$ between the particle and Neptune at every perihelion passage
(following, e.g., Wang & Malhotra 2017; Malhotra et al. 2018; Lan & Malhotra
2019); $\psi$ is illustrated in the top panel of Figure 1. For particles with
perihelion distances consistent with scattering TNOs ($q\approx 33-37$ au),
those on resonant orbits remain relatively stable in semimajor axis, tracing
out distinct resonant islands in $a$-$\psi$ space, while those on scattering
orbits wander randomly through that space filling in a sea of background
points.

Figure 1 shows example $a$-$\psi$ evolution for $q=33$ au particles in and
near the 2:1, 7:3, and 5:2 resonances in the circular restricted three-body
problem with the Sun and Neptune as perturbers. The details of our methods and
how these maps are generated are described in Section 2.1, but we briefly
highlight here a few key features of Neptune’s external resonances. In
general, a particle in a $p$:$q$ external resonance will visit $q$ distinct
resonant islands over a resonant cycle; so an N:2 resonant particle (such as
the 5:2 particles in Figure 1) trace out two islands, N:3 particles trace out
three islands, and so on. Particles in external N:1 resonances follow this
pattern in that individual particles trace out a single island, but the
structure of these resonances is distinct in that there are three possible
libration islands to choose from. These are the two so-called asymmetric
islands centered on $\psi\sim 90^{\circ}$ and $\psi\sim 270^{\circ}$ (note in
Figure 1 that these two islands are traced by distinct, differently colored
particles) and the larger symmetric island centered on $\psi=180^{\circ}$. We
refer the reader to Volk & Malhotra (2022) for additional details on the
structure of these resonances. In our later maps of Neptune’s resonances in
various simulations, we examine how all of these resonant islands become
stronger, weaker, or disappear entirely.

In Section 2.1, we describe how the simplified three-body model shows that
Neptune’s resonances should remain strong out to semimajor axes of several
hundred au in the scattering population. We then expand our simulations to map
Neptune’s resonances under the influence of the other giant planets (following
the methods described in Volk & Malhotra 2022) to show that the distant
resonances become much weaker for scattering TNOs. In Section 2.2 we show that
Uranus is primarily responsible for this weakening. In Section 2.3 we describe
how we quantify resonance strengths in our simulations and present the
present-day strengths of Neptune’s resonances in the scattering population.

---
Figure 1: Top: Illustration of the angle $\psi$, which describes a TNO’s
position relative to Neptune when the TNO is at perihelion. Bottom: Example
evolution of resonant (colored) and non-resonant (gray) particles with $q=33$
au in $a$-$\psi$ space for the 2:1, 7:3, and 5:2 resonances in the circular
restricted three-body problem (Sun + Neptune + massless test particles).
Within each panel, individual resonant particles are given distinct colors;
note the clear islands traced out by the resonant particles.

###  2.1 Exploring resonances in the scattering TNO population

We start by comparing the $a$-$\psi$ evolution of particles in Neptune’s
resonances in the simple Sun-Neptune-TNO problem to those in simulations with
all four giant planets. To do this, we created an illustrative $a$-$\psi$
mosaic from 165 au to 285 au in the circular restricted three-body problem, as
shown in the left panel of Figure 2. We included the Sun and Neptune as the
only perturbers in our simulations, with Neptune on a circular orbit with its
current semimajor axis. We then initialized massless test particles with
perihelia of 33 au and inclinations of 0∘ along lines of constant semimajor
axis spanning from $\psi=0^{\circ}$ to $\psi=360^{\circ}$ at the expected
location of all the N:1, N:2, and N:3 resonances in this semimajor axis range.
We ran the simulation forward through 1000 resonant cycles for each target
resonance in each simulation using the rebound ias15 integrator (Rein & Liu,
2012; Rein & Spiegel, 2015), recording the orbit of each particle every time
it passed through perihelion. The left panel in Figure 2 shows the position of
the particles at each perihelion passage in $a$-$\psi$ space. We plotted any
non-resonant particles, loosely defined as any that deviated 5 au or more in
$a$ from their starting points, in light gray and the particles that remain
resonant in darker colors; this highlights the resonant particles against the
background chaos of non-resonant particles.

We then generated a similar $a$-$\psi$ map from 165 au to 285 au for Neptune’s
resonances in the present-day solar system with the Sun and all four giant
planets as perturbers for the test particles. To do this, we queried JPL
Horizons for the positions, masses, and velocities of the Sun, Jupiter,
Saturn, Uranus, and Neptune. We then initialized massless test particles with
perihelia of 33 au and inclinations of $0^{\circ}$ along four lines of
constant $\psi$ spanning the entire semimajor axis range: $\psi=70^{\circ}$
and $\psi 290^{\circ}$ (corresponding to the centers of the N:1 resonances),
$\psi=90^{\circ}$ (corresponding to the center of the N:2 resonances), and
$\psi=180^{\circ}$ (corresponding to the center of the N:3 resonances).
Initializing particles at constant $\psi$ rather than along lines at the
expected resonant semimajor axes is more efficient in the non-idealized
problem; we found that with the four choices of $\psi$, a resolution of 0.06
au in $a$ was sufficient to resolve all the resonance of interest. We note
that because the planets are not co-planar, test particles initialized with
$i=0$ in the full simulations results in a small inclination dispersion
relative to Neptune. As above, we integrated these particles with the ias15
integrator and recorded their orbits at every perihelion passage. The
resulting $a$-$\psi$ map is shown in the right panel of Figure 2 with the same
color scheme as above.

Neptune only | all four giant planets
---|---
|
Figure 2: Semimajor axis (a) vs. $\psi$ for particles in Neptune’s distant
resonances, with perihelia of 33 au. On the left, only the Sun and Neptune are
included as perturbers. On the right, the Sun and all four giant planets are
included as perturbers. Note how the N:1 resonances in the 4-planet problem
get weaker as the semimajor axis increases and begin to drop out entirely at
the 20:1 position, with a brief resurgence around the 26:1 position. This
contrasts the more gradual shrinking of distant resonances in the three-body
problem

This comparison reveals that the overall structure of Neptune’s resonances is
significantly affected by the presence of the other giant planets. Most
notably, the strengths of the N:1 resonances drop off sharply by $a\sim 220$
au, which is where the 20:1 would be located, before briefly reappearing at
$a\sim 260$ au. Hand-measured widths of the stable libration zones of the N:1
resonances in the full simulations show that they are significantly weaker in
this region when compared to the same measurements from the simplified
simulations (as seen when comparing the left and right panels of Figure 2).
Additionally, in the full problem the N:2 resonances disappear slightly past
the position of the 20:1 at $\sim 225$ au (41:2), and the N:3 resonances drop
out at $\sim 180$ au (44:3). In the simplified model, the N:2 resonances
continue over our entire semimajor axis range, and the N:3 resonances persist
to $a\sim 235$ au.

###  2.2 Identifying what is destabilizing resonances in the distant
scattering population

To investigate the mechanism responsible for destabilizing Neptune’s distant
resonances, we consider a case study centered on Neptune’s 20:1 resonance, the
lowest-$a$ N:1 resonance that lacks a clear stable libration zone in the full
simulations (right panel of Figure 2). We first ran a set of simulations
targeting just the 20:1 resonance over a range of perihelion distances
$q=33-45$ au in both the simplified model and the full, four giant planet
model. For each value of $q$ (varied in 1 au increments), we initialized 240
massless test particles evenly spaced in $\psi$ from $0^{\circ}$ to
$360^{\circ}$ at the center of the 20:1 resonance ($a\approx 221.6$ au). The
inclination of all test particles were set to zero; in the simplified model
this makes them truly co-planar with Neptune while in the full model they have
very small inclinations relative to Neptune. We again used the ias15
integrator to integrate forward over 1000 resonances cycles (20,000 orbits of
Neptune for the 20:1 resonances), recording the test particle’s orbits every
perihelion passage to produce $a$-$\psi$ maps. The boundaries of the 20:1
resonance’s stable libration zones in these $a$-$\psi$ maps of both the
simplified and full models were hand measured and are plotted in $a$-$e$ space
in Figure 3 to show the width of the 20:1 resonance as a function of $e$; this
figure also shows an extended eccentricity range to higher-$q$ for the
simplified model (which is much less computationally intensive to run). While
the maximum width of the 20:1 resonance is similar in both models, its
location is shifted to higher $q$ in the full model compared to the Neptune
and Sun only case. There is a significant decrease in the width and strength
of the 20:1 resonance at low $q$ when all giant planets are included in the
simulation.

Figure 3: Eccentricity vs. semimajor axis of the 20:1 resonance with Neptune.
In gray are resonance widths with only Neptune and the Sun acting as
perturbers. In blue are resonance widths with all four giant planets included
as perturbers. Notice how resonance width decreases with higher eccentricity
in the four-planet case relative to teh simplified case.

This led us to investigate whether one of the three giant planets other than
Neptune is destabilizing the 20:1 resonance at higher eccentricities. Using
the same simulation set up as described above in Section 2.1 for generating
the four-planet model $a$-$\psi$ maps, we ran six simulations of Neptune’s
resonance in the semimajor axis region surrounding the 20:1 with the following
objects included:

  * •

Sun + Saturn + Uranus + Neptune + TNOs

  * •

Sun + Jupiter + Uranus + Neptune + TNOs

  * •

Sun + Jupiter + Saturn + Neptune + TNOs

  * •

Sun + Jupiter + Neptune + TNOs

  * •

Sun + Saturn + Neptune + TNOs

  * •

Sun + Uranus + Neptune + TNOs

where in each case the Sun and planets are massive perturbers and the TNOs are
massless test particles.

In doing this, we were able to isolate the effects of each planet on Neptune’s
resonances. We found that when just Jupiter and/or Saturn were removed from
the simulation, the 20:1 resonance remained weak at $q=33$ au. However,
whenever Uranus was removed from the simulation, the 20:1 resonance was nearly
as strong as predicted by the simple three-body problem and did not drop out
relative to the neighboring resonances. We also ran a simulation with the Sun,
Jupiter, Saturn, and Neptune as massive perturbers but added a J2 term to
represent the averaged secular effect of Uranus on the other planets and the
TNO test particles (we used Malhotra & Ito (2022)’s J2 value and formulation
for the simulation modifications). The 20:1 resonance in this simulation
behaved very similarly to the simulation without Uranus, though some of the
higher-order resonances nearby appear less distinct. Figure 4 shows a
comparison between a few of these test simulations.

|
---|---
|

Figure 4: Barycentric semimajor axis ($a$) vs. $\psi$ for particles near
Neptune’s 20:1 resonance ($a\approx 221$ au) with perihelion distance $q=33$
au. In all panels, we plot particles that remain within 10 au of their initial
semimajor axis in colored points, with each particle plotted in a different
color (these are the particles that tend to interact with the resonances);
less stable particles that quickly scatter away from their starting $a$ are
plotted as gray background points. In the top panel, only the Sun and Neptune
act as perturbers on the test particles. In the middle left panel, the Sun and
all major planets except for Uranus act as perturbers. In the middle right
panel, the Sun and all 4 major planets and the Sun act as perturbers. In the
bottom left panel, all giant planets except Uranus act as perturbers, with the
addition of a J2 parameter to simulate Uranus’ secular effects. In the bottom
right panel, the Sun and all four giant planets are included as perturbers in
the integration of the massive bodies, but the test particles are only
perturbed by the Sun, Jupiter, Saturn, and Neptune. Note the similarity
between the simplified (Sun + Neptune) model at the top and the middle left
panel with Uranus excluded from the simulation. The destabilization of the
20:1 resonance only occurs when Uranus is included in the simulation; the
destabilization remains when Uranus perturbs Neptune’s orbit but is not
included as a perturber for the test particles themselves.

From these tests, we can conclude that the weakening of the 20:1 resonance at
$q=33$ au is related to the presence of Uranus, as the resonance is weak
whenever Uranus is present in the simulation, but unaffected by the presence
of Jupiter or Saturn. We also conclude that the weakening is not purely or
dominantly due to the secular effects of Uranus. It is clear that
gravitational perturbations by Uranus can have a significant impact on
Neptune’s distant resonances for high eccentricity TNOs–a regime typically
assumed to be dominated only by the gravitational effects of Neptune. We ran
one more test simulation near the 20:1 resonance to examine whether the
weakening of the resonances when Uranus is present is due to Uranus’s direct
influence on the test particles (via, e.g., resonant interactions with Uranus)
or is due to Uranus’s influence on the other giant planets’ orbits; as we will
discuss in Section 3, Neptune and Uranus are near a mutual resonance, so
Neptune’s orbital evolution is influenced significantly by non-secular effects
from Uranus. In this test simulation, we modified rebound’s whfast integrator
(Rein & Tamayo, 2015) such that Uranus was included as a perturber for the Sun
and the other giant planets, but was not included in the force calculations
for the test particles111We chose to use whfast for this modification because
it involved just a very simple alteration to the ‘kick’ subroutine for the
test particles within the drift-kick-drift scheme for the whfast algorithm. We
ran a test to ensure that the whfast results for a full Sun plus four giant
planets plus test particles simulation are comparable to the results of our
ias15 simulations. The switch in integrator for this test does not affect the
results.. The results of this test simulation are also shown in Figure 4, and
the resonance structures are the most similar to the full Sun plus all four
giant planets simulations. This hints that, near the 20:1 resonance at least,
the weakening of Neptune’s resonances is caused by variations in Neptune’s
orbit due to its near resonance with Uranus.

We explore in more detail how Uranus is affecting a wider range of Neptune’s
resonances in Section 3. But, first, we must establish a more quantitative
method to measure resonance strengths in our simulations.

###  2.3 Quantifying resonance strengths

Previously, we have used resonance width (in semimajor axis) to characterize
the strength and stability of Neptune’s external resonances. This is a useful
metric in many circumstances, such as in the Sun-Neptune-TNO three-body
problem. In the three-body case, Neptune’s resonances are have generally clear
boundaries in $a$-$\psi$ phase space. These well-defined boundaries break down
and become more chaotic when all giant planets are included in our
simulations. This makes identifying the boundaries of resonances in the four-
planet model more subjective, as the dimensions change depending on the
parameters used to plot the test particles (such as limiting test particles by
change in semimajor axis over a specified simulation time). This is especially
difficult when studying scattering (low-$q$) populations. Because of this we
sought a different, more concretely defined, metric for measuring the strength
and stability of resonances in the full six-body problem (Sun + four giant
planets + TNO).

The following describes the basic behavior of some typical low-$q$ test
particles simulated in the full six-body case. Test particles initialized near
the edge of more stable, closer-in resonances will often move in and out of
resonance over the span of a simulation due to chaotic gravitation
interactions with the giant planets. Therefore, these test particles cannot be
defined as explicitly in or out of resonance, but rather are transient. This
is one difficulty when it comes to measuring widths of resonances in the six-
body problem: do we count transient particles as resonant or not? In the case
of further out resonances (such as the 20:1 at $q=33$ au and beyond), there
may be no simulated test particles that are persistently resonant over the
timescale we are considering. But that does not mean there is no resonant
interaction, and it would be inaccurate to say that these resonances have a
width or strength of zero. So somehow, we must account for test particles that
transiently pass in and out of, or stick, to each resonance; example typical
sticking behavior is demonstrated in Figure 5. Additionally, even test
particles that spend most of their time librating in a resonance will often
‘slip’ in $\psi$ at chaotic boundaries, so purely considering confinement of
$\psi$ doesn’t provide a full measure of the stability and strength of a
resonance Overall, using the width to measure these resonances based only on
long-term stable libration does not accurately capture the strength of these
sticking interactions that are very important in the scattering TNO
population.

Figure 5: Example evolution of a test particle that scatters and sticks to
Neptune’s resonances. We show the evolution in $a$ vs time (upper left),
$\psi$ vs time (upper right), $q$ vs time (lower left) and in $a$-$\psi$ space
(lower right). We highlight sticking events in the 35:2 ( $\psi$ alternates
between $\sim 90^{\circ}$ and $\sim 270^{\circ}$ each resonant cycle) and the
17:1 ($\psi$ first librates in the leading asymmetric island then the trailing
island).

In these situations, it is far more useful to look at the total resonant
lifespan of test particles in a given simulation. That is, how long do test
particles spend interacting with a given resonance as a fraction of the total
simulation time, even if these interactions are not continuous. We quantify
this fractional resonance occupation time averaged over all the test particles
in a given region around a resonance as our resonance strength parameter
(exact details of its calculation are provided in Appendix A). Figure 6 shows
the relative strengths of Neptune’s distant resonances at $q=33$ au using this
metric. This method is also utilized in the next section to quantify resonance
strength.

Figure 6: Relative strength of the N:1 (black), N:2 (pink), and N:3 (yellow)
resonances with Neptune for low-inclination praticles with $q=33$ au in a
model of the current solar system (with the Sun and all four giant planets as
perturbers). Our strength metric is normalized to the largest value across all
resonances and is a measure of how confined test particles are to each
resonance (see Appendix A).

The dramatic drop-off in our resonance strength metric is consistent with
scattering population simulations that show a dramatic reduction in resonance
sticking at large semimajor axes ($a\gtrsim 250$ au; Lykawka & Mukai 2007). We
have established that Uranus is directly linked to this drop-off in resonance
strength (Figure 4). Neptune’s current orbital period is 1.96 times Uranus’s
orbital period, which is close to their mutual 2:1 resonance. In the next
section, we explore how this near-resonance might contribute to weakening
Neptune’s distant external resonances.

##  3 Experimenting with the Neptune-Uranus period ratio

As discussed in Section 2.2, it is clear that Uranus has a significant
influence on Neptune’s distant resonances, and there are hints that a lot of
that influence is due to how Uranus affect’s Neptune’s orbital evolution.
Uranus and Neptune are not far from a mutual 2:1 resonance, so here we we
explore how the proximity to that mutual resonance can reduce stability in
Neptune’s distant resonances. To test this we ran a series of simulations with
Neptune’s semimajor axis shifted slightly from its present day value to change
its mutual orbital period ratio with Uranus. We chose to simulate the phase
space between Neptune’s 2:1 and 28:1 resonances at q = 33 au, for Neptune to
Uranus (N:U) period ratios of 1.96 (present-day value), 2.00, 2.02, and 2.36.
A period ratio of 2.02 shows us how resonance stability behaves slightly
beyond this mutual 2:1, and 2.36 shows us resonance stability when Neptune and
Uranus are far from their mutual 2:1.

For each N:U period ratio simulation we included all four giant planets with
their orbital parameters being queried from JPL Horizons in barycentric
coordinates. We then changed only Neptune’s semimajor axis to correspond with
our desired N:U period ratio. All of Neptune’s other orbital parameters were
left at their present day values, and we recomputed the new barycenter of the
simulated system. Before adding test particles to this modified system, we
made sure that the giant planet orbits remained stable when changing Neptune’s
orbital parameters. We found that, for each of our four chosen period ratio
values, Neptune’s semimajor axis and eccentricity remained stable for the
maximum duration of our 1000 resonant cycle simulations, which is 4.6 Myr for
the most distant resonance we tested. There was no significant destabilization
of orbital parameters for any of the other giant planets on this timescale.

For each of these simulated N:U period ratios, we then ran a set of
simulations with test particles spanning the expected locations of Neptune’s
N:1 resonances from the 2:1 to the 28:1, with a $\pm 0.5$ period ratio range
around each N:1 to include the nearby higher-order resonances. For example,
the simulation centered on the 10:1 resonance includes particles spanning from
from Neptune’s 19:2 through 21:2 resonances. We introduced 1000 massless test
particles to each simulation. For the N:U=1.96 simulation (Neptune’s current
day position), we assigned all particles a perihelion distance of 33 au to
match our previous simulations of scattering bodies. When moving Neptune’s
semimajor axis to change the N:U period ratio, we scaled the test particle
perihelia: $q=33$ au with Neptune at $a=30.06$ au (N:U$=1.96$), $q=33.3$ au
with Neptune at $a=30.37$ au (N:U$=2.00$), $q=33.65$ au with Neptune at
$a=30.65$ au (N:U$=2.02$), and $q=37.3$ au with Neptune at $a=34$ au
(N:U$=2.36$). For each simulation, the 1000 particles were allocated uniformly
across the N$\pm 0.5$ period ratio range along lines of constant $\psi$ values
as follows: 500 particles targeted the centers of the asymmetric N:1 resonant
islands (250 at $\psi=70^{\circ}$ for the leading island, and 250 at
$\psi=70^{\circ}$ for the trailing island), 250 targeted the centers of the
N:2 resonances at $\psi=90^{\circ}$, and 250 targeted the centers of N:3
resonances at $\psi=180^{\circ}$. The particles were initialized at perihelion
(mean anomaly of zero), so their longitude of perihelia were determined by the
desired $\psi$ value. Their semimajor axis and eccentricity were set by the
period ratio with Neptune and fixed perihelion distance. We set inclinations
and longitudes of ascending node to zero; as noted in Section 2.1, because the
planets are not co-planar in our simulations, this amounts to allowing the
test particles to have a small inclination dispersion relative to Neptune (a
few degrees at most over the timescale of our simulations). We found that
these initial parameters for test particles were adequate to resolve the
general structure of Neptune’s resonances when plotted in phase space and
calculate resonance strengths.

As before, we ran these simulations for a duration of 1000 resonant cycles (of
the target N:1 resonance in each simulation) using ias15 for each desired N:U
period ratio. We analyzed and plotted resonance strengths using the methods
described in Section 2.3 and Appendix A. The resulting mosaics of individual
particle residence times near the investigated resonances are shown Figure 7
for these four N:U period ratios, and the calculated strengths of each N:1
resonance for the same N:U period ratios are shown in Figure 8.

We ran additional simulations from various other N:U period ratios to get a
better sense of the overall trends in Neptune’s distant resonances. For N:U
period ratios between 1.96 and 2.00, as well as between 2.00 and 2.02,
resonance stability drops off steeply. For example, a residence time plot at
N:U$=1.98$ looks very similar to N:U$=2.02$. Between N:U period ratios of 2.02
and 2.36, resonance strength returns steadily. In order to visualize these
trends in overall resonance strength with N:U period ratios, we averaged
individual resonance strength values across a range of N:1 resonances for each
set of N:U period ratio simulations. This “average strength” value represents
the general weakening of all simulated resonances for various N:U period
ratios. This is shown in Figure 9, where we average the strengths of Neptune’s
closer-in N:1 resonances (2:1 through 14:1; left panel) and distant N:1
resonances (15:1 through 28:1; right panel). To provide a baseline for
comparison, at each modified value of Neptune’s semimajor axis we ran a set of
simulations with Uranus removed from the simulation and calculated the
strengths of Neptune’s 2:1 through 28:1 in the absence of Uranus’s
destabilizing effects.

Figure 9 shows that the reduction in resonance strength due to Uranus is much
more dramatic for the distant resonances than the closer-in resonances,
regardless of the specific N:U period ratio. The weak trend in the baseline
simulation’s average distant resonance strength (right panel) shows that the
other giant planets do influence the strength of Neptune’s distant resonances
but to a much lesser extent than Uranus. For both the closer-in and distant
resonances, the present day N:U period ratio is actually more favorable for
Neptune’s resonance strengths than N:U period ratios smaller but farther from
their mutual 2:1; the general proximity of Neptune to Uranus, which increases
Uranus’s influence on Neptune’s orbital evolution, is a significant factor in
the destabilization of Neptune’s resonances. It is also clear that if Neptune
and Uranus were any closer to their mutual 2:1 resonance than they are today,
all of Neptune’s external resonances would be significantly weakened. We
discuss the implications of this further in Section 4.

|  |  |  |
---|---|---|---|---
Figure 7: Test particle period ratio with Neptune vs. $\psi$ in the phase
space from Neptune’s 2:1 to 28:1 resonances ($\sim 50-280$ au) for four
different Neptune : Uranus period ratios. In each simulation, particles are
initialized with the equivalent of $q=33$ au for present-day Neptune’s orbit.
Dark blue test particles indicate a longer normalized residence time, which is
defined as a function of its change in semimajor axis ($\Delta$a) across
multiple small time intervals (see Appendix A). Lighter yellow test particles
indicate a shorter normalized residence time. Particles assigned a residence
time of 0 (i.e., particles that scatter for the entire simulation) are plotted
in light gray for clarity to separate them from the resonant-interacting
particles. The Neptune : Uranus period ratio for each set of simulations are,
from left to right: 1.96 (their present-day period ratio), 1.96, 2, 2.02, and
2.36. Uranus is included as a full perturber in all but the left-most plot; in
the limited Uranus effects simulation, Uranus interacts with the other massive
bodies but not directly with the test particles. Notice the dramatic
destabilization of Neptune’s external resonances when Neptune is set to be in
the 2:1 resonance with Uranus. We discuss in Section 3 how the stability of
Neptune’s external resonances depends on the the mutual period ratio between
Uranus and Neptune. In the leftmost figure, note that the destabilization
remains when Uranus perturbs Neptune, but is not included in the force
calculations for the test particles. Figure 8: Relative sticking strengths of
Neptune’s N:1 resonances for different Neptune : Uranus period ratios. The
black triangles and dashed line show the relative resonance strengths in the
current solar system (Sun plus all four giant planets on their present-day
orbits); the gray circles and dashed line shows the strengths with Uranus
removed. The four sets of colored points and lines show strengths from
simulations with Jupiter, Saturn, and Uranus on their current orbits, and
Neptune’s semimajor axis adjusted to the labeled period ratio with Uranus
(these are the same simulations shown in Figure 7). The purple colored points
and lines show strengths from simulations where Uranus is included as a
perturber for the other massive bodies but not for the test particles. In all
simulations, the strengths have been normalized such that Neptune’s external
2:1 resonance has a strength of 100. There is a very rapid drop-off in N:1
resonance strengths with increasing N at the exact N:U 2:1 resonance. The
resonance strengths plateau past the $\sim 5:1$ in the simulation with Uranus
removed. In the simulations with Uranus, the destabilization of Neptune’s
distant N:1 resonances is more pronounced when the N:U period ratio is closer
to 2:1. The simulations with limited gravitational effects from Uranus on the
test particles show similar destabilization. Figure 9: Resonance strengths as
a function of the Neptune : Uranus period ratio averaged across Neptune’s 2:1
through 14:1 resonances (left panel) and 15:1 through 28:1 resonances (right
panel). The blue points and dashed lines show the averaged strengths in
simulations with the Sun and all four giant planets, with Neptune’s semimajor
axis increased or decreased to change its period ratio with Uranus. The black
points and dashed lines show the averaged strengths in simulations with the
same values of Neptune’s semimajor axes but with Uranus removed (i.e.,
simulations with the Sun, Jupiter, Saturn, and Neptune as perturbers). The red
stars highlight the present day N:U period ratio (1.96) in the simulations
with Uranus removed; the average resonance strengths measured each simulation
are normalized relative to these two values (both red stars are at y=100). The
purple triangles represent the averaged strength from a set of simulations
where Uranus perturbs the orbits of Neptune and the other giant planets, but
is not included in the force calculations on the test particles. Note that
these values are similar to the those of the 1.96 N:U period ratio simulations
with the full effects of Uranus included. Neptune’s more distant resonances
(right panel) are much more dramatically weakened by the presence of Uranus
across all tested N:U period ratios. All of Neptune’s resonances show the most
dramatic weakening at period ratios close to the N:U 2:1, and less weakening
at higher N:U period ratios.

We performed a few additional simulations to better isolate which aspect of
Uranus’s influence on Neptune is the most important, and whether any of the
weakening is due to more direct effects of Uranus on the test particles
themselves. One obvious influence Uranus has on Neptune is that it causes
Neptune’s eccentricity to vary over time, and the amplitude of these
variations increases as the N:U period ratio approaches 2:1; we tested whether
increasing Neptune’s eccentricity variations in the absence of Uranus itself
weakened Neptune’s resonances. We ran two sets of simulations with Neptune at
its present day semimajor axis and only the Sun, Jupiter, and Saturn as
additional perturbers. In these simulations, we imposed sinusoidal variations
in Neptune’s eccentricity and precession of its longitude of perihelion to
mimic its present day evolution as well as amplified variations to mimic
proximity to the N:U 2:1 without actually including Uranus in the simulation;
in the latter case, we based the variation on the N:U$=1.98$ simulations in
which Neptune’s eccentricity maximum was nearly double its present day value.
To do this, we used fictitious forces within rebound similar to those
described in Wolff et al. (2012) for the Mercury integrator (the
implementation of these integrator modifications within rebound is described
in Hermosillo Ruiz et al. in prep). In both scenarios, we modeled test
particles in Neptune’s 2:1 through 28:1 resonances, calculated the associated
resonance strengths, and visually inspected maps like those shown in Figure 7.
There were no significant differences in Neptune’s resonances as a result of
increasing Neptune’s eccentricity; both of these test simulations behaved
essentially the same as our baseline simulations with Uranus removed.

The other notable aspect of Neptune’s evolution that can be attributed to
Uranus is that Neptune’s barycentric semimajor axis experiences sinusoidal
variations with a period of $\approx 4300$ years; this is the period
associated with the near 2:1 N:U resonance, and the amplitude of these
variations are larger when the planets are closer to the resonance. This
period is equal to a bit more than 26 times Neptune’s orbital period and thus
equal to the orbital periods of TNOs near our most distant set of N:1
resonances where the resonance strengths drop off (though we see the drop-off
being closer in than the 26:1). We note that there is also a longer-period
modulation of the amplitude of Neptune’s barycentric semimajor axis variations
that has a period of $\sim 1.15$ Myr, though that is significantly longer than
the orbital or resonant libration timescales for the N:1 resonances we
consider. To see whether Neptune’s resonances regained any strength at larger
distances where the orbital periods are better separated from this short-term
period in Neptune’s $a$ variations, we extended our present-day solar system
configuration simulations out to Neptune’s 50:1 resonance; there was no
apparent return of resonant stability for particles with $q=33$.

We ran one last set of simulations to determine how much of the observed
weakening of Neptune’s present-day resonances in the scattering population is
attributable just to the Uranus-induced variations in Neptune’s orbit. As
described in Section 2.2 for simulations near Neptune’s 20:1 resonance, we ran
a set of simulations for the 2:1 through 28:1 resonances in which the Sun and
all four giant planets were fully interacting on their present-day orbits, but
the test particles did not feel any direct gravitational effects from Uranus.
The residence time mosaic for this set of simulations is shown in Figure 7 (we
labeled these as ‘limited Uranus’ in the figures) immediately alongside the
full simulation. There is a remarkable similarity between the two scenarios,
especially at and beyond Neptune’s 20:1 resonance. The intermediate N:1
resonances (between the 10:1 and 20:1) do appear stronger when Uranus does not
directly perturb the test particles, and this is reflected in the resonance
strength measurements shown in Figure 8; we also include the averaged
resonance strengths from the limited Uranus simulations in Figure 9. When we
compare the resonance strengths in Figure 8 from the simulations with Neptune
on its current orbit without Uranus, with limited Uranus effects, and with the
full influence of Uranus, we note only small differences in the strengths of
Neptune’s resonances out to the $\sim 10$:1, at which point the full
simulation resonances strengths begin to decrease significantly; the limited
Uranus simulation resonance strengths maintain similarity to the no Uranus
simulations a bit further out before dramatically decreasing at the $\sim
15$:1. It is clear that the variations in Neptune’s orbit due to Uranus are
responsible for most of the weakening in Neptune’s distant resonances, though
more direct perturbations from Uranus have some influence on the intermediate
resonances.

To briefly summarize what we conclude from our numerical simulations:

  * •

In the absence of Uranus, Neptune’s exterior mean motion resonances in the
scattering population would be significantly stronger and extend out to much
larger period ratios.

  * •

At the present day Neptune to Uranus period ratio (N:U = 1.96), nearly all of
the weakening in Neptune’s most distant resonances (beyond the $\sim 20$:1) in
the scattering population can be attributed to the variations in Neptune’s
orbit induced by the near-resonance with Uranus (particularly the $a$
variations) rather than to Uranus’s direct gravitational influence on the
scattering population

  * •

Neptune’s intermediate resonances in the present-day solar system (the $\sim
10$:1 through the $\sim 19$:1 are slightly weakened by Uranus’s direct
gravitational influence.

  * •

All of Neptune’s external resonances would be significantly weaker if Neptune
and Uranus were closer to their mutual 2:1 resonance, essentially vanishing at
N:U=2.

We discuss the implications of these points in Section 4.

##  4 Discussion

Our exploration of Neptune’s resonance strengths in the current giant planet
configuration and alternative configurations has several important
implications for the current and past evolution of the TNO populations. First,
we now have an explanation for the drop-off in resonance sticking at large $a$
in the scattering TNO population (seen by, e.g., Lykawka & Mukai 2007 in
numerical models of the scattering population). The drop-off occurs because
variations in Neptune’s orbit (particularly its semimajor axis) induced by the
near 2:1 resonance with Uranus essentially destabilizes Neptune’s resonances
beyond $a\sim 200-250$ au entirely for nearly Neptune-crossing eccentricities
($q=33$ au). If Neptune had migrated further outward by several au in
semimajor axis in the early solar system and was thus better dynamically
separated from Uranus, the present-day scattering TNO population would
experience more resonance sticking at large semimajor axes. The variations in
Neptune’s orbit also slightly weaken Neptune’s closer-in resonances, though
there is additional weakening from Uranus’s direct gravitational influence on
the scattering population.

The simplest explanation for this additional destabilization of Neptune’s
present-day resonances by Uranus is overlap between Neptune’s external
resonances and Uranus’s external resonances. This likely explains the rapid
deterioration of stable libration in Neptune’s resonances for the simulations
in which we placed Neptune and Uranus at their exact 2:1 period ratio. When
Neptune and Uranus are in their mutual 2:1 resonance, each of Neptune’s
exterior N:1 and N:2 resonances will be overlapped by an exterior N:1
resonance of Uranus. This overlap between multiple strong resonances causes
the strengths of Neptune’s N:1 resonances to drop off very quickly in Figure
8; in Figure 7 there is no apparent stable libration in any resonance with
Neptune beyond the $\sim 5:1$ for the simulation with N:U = 2. When Neptune
and Uranus are not at their exact 2:1 period ratio, there will still be places
where their exterior resonances happen to overlap or otherwise interact.
Gallardo (2020) discusses how Uranus’s resonances are faintly visible in
dynamical maps of Neptune’s resonances (they highlight Uranus’s 4:1 resonance
just exterior to Neptune’s 2:1 in their Figure 23). They note that Uranus’s
resonances are stronger at higher eccentricities, such as those in the low-$q$
scattering population we consider here (see also Figure 4 in Gallardo 2018).
This is likely why Neptune’s distant resonances regain strength at large $a$
for the lower-$e$, higher-$q$ orbits in the detached population (see Figure 3
and Volk & Malhotra 2022).

Overlap between Neptune’s external resonances and three-body resonances
involving Neptune, Uranus, and TNOs also potentially contribute to the
weakening of Neptune’s resonances. Many main belt asteroids are known to exist
in three-body resonances, mostly involving Jupiter and Saturn (e.g. Gallardo,
2018). Given the short orbital periods of main belt asteroids, their orbital
parameters are typically much more tightly constrained than TNOs; it is thus
more rare for TNOs to have orbit-fits accurate enough for secure
identification in resonances as narrow as three-body resonances tend to be.
Gallardo (2018) discuss 3 TNOs likely in the the 2 +1U-3N three-body resonance
with Neptune and Uranus (a resonance identified by Smirnov & Markov 2017);
they conclude that the influence of three-body resonances on Neptune’s two-
body resonances leads to chaotic diffusion but is not a problem that has been
sufficiently well studied. Overlap between Neptune’s external two-body
resonances, Uranus’s external two-body resonances, and potential three-body
resonances with both planets likely plays a role in setting some of the
smaller scale, noisier structures seen in the patterns of resonance strengths
in Figure 8.

Beyond resonance sticking in the present-day scattering population, the strong
dependence of Neptune’s distant resonance strengths on the Neptune to Uranus
period ratio (and on Neptune’s orbital variations due to interactions with
Uranus) has very important implications for using Neptune’s current resonant
populations to constrain the migration history of the giant planets. As the
orbits of both Neptune and Uranus evolved during the epoch of giant planet
migration, their period ratio and amount of dynamical coupling would vary over
time. Thus the strengths of Neptune’s distant external resonances could vary
significantly depending on the exact interactions between Uranus and Neptune;
this could dramatically change the likelihood of capturing TNOs into Neptune’s
resonances at different times during migration. In particular, overall
predictions for distant resonance populations from planet migration
simulations could be very sensitive to how close Neptune and Uranus approach
their mutual 2:1 resonance, with distant resonance capture probabilities
likely dropping to near zero if the two planets spend significant time too
close to that resonance.

The observed semimajor axis distribution of TNOs that are stably in Neptune’s
distant resonances and those on stable orbits that are near but not in
resonance provide powerful constraints on the migration history of the
planets. The stable distant resonant TNOs are likely scattering TNOs that
stuck to the resonances while Neptune was still migrating and were pushed into
more stable phase space. The strength of the resonances during migration, and
thus the probability of sticking to them and being captured, varies depending
on the N:U period ratio (Figure 9). Thus the relative populations of captured
TNOs across a wide range of Neptune’s distant resonances will be a sensitive
probe of the period ratio history of the two planets, though the observed
populations are a combination of stable and transiently resonant objects (e.g.
Crompvoets et al., 2022) that must be disentangled. The non-resonant TNOs on
stable (i.e., not scattering orbits) sunward of these distant resonances are
often interpreted as resonant drop-outs. These could be TNOs that stuck to the
resonances and evolved through secular interactions to low-enough
eccentricities to be left behind as Neptune migrated (e.g. Gomes et al.,
2005). They also could have been dropped from resonance by discrete jumps in
Neptune’s semimajor axis during migration as Neptune interacted with the most
massive planetestimals in the pre-migration disk (e.g. Nesvorný et al., 2016;
Kaib & Sheppard, 2016; Lawler et al., 2019). The sensitivity of Neptune’s
resonance strengths to the N:U period ratio hints at an additional
possibility: resonance drop-outs could also occur if the two planets approach
too closely to their mutual 2:1 resonance during migration. While our
investigation only focused on the low-$q$ population, and Neptune’s distant
resonances at the present-day N:U period ratio remain strong for the
higher-$q$ orbits associated with the resonance drop-outs, the destabilization
of Neptune’s resonances is very dramatic as the planets approach their mutual
resonance; it is likely that Neptune’s resonances would be significantly
weaker even for those lower-$e$ orbits.

##  5 Summary and Conclusions

We have explored the phase space of Neptune’s mean motion resonances from
$a\approx 50-285$ au for eccentricities typical of scattering TNOs, finding:

  * •

In the three-body problem (Sun, Neptune, and test particles) as well as in a
solar system with Uranus removed (Sun, Jupiter, Saturn, Neptune, and test
particles) Neptune’s external mean motion resonances would remain strong out
to several hundred au in the scattering TNO population.

  * •

In the current solar system (Sun, Jupiter, Saturn, Uranus, Neptune, and test
particles), Neptune’s resonances become unstable beyond $a\approx 200-250$ au
for low-$q$ orbits, explaining the lack of resonance sticking seen in
simulations of the high-$a$ scattering TNO population.

  * •

The destabilization of Neptune’s most distant resonances in the scattering
population is due to variations in Neptune’s orbit that are caused by the
proximity of Uranus and Neptune to their mutual 2:1 resonance.

  * •

Uranus has a more direct influence on Neptune’s closer-in resonances,
weakening them sightly compared to simulations in which Uranus is not present;
this is likely the result of overlap between Neptune and Uranus’s two-body
resonances as well as potentially weaker three-body resonances.

  * •

The strengths of all of Neptune’s external resonances are quite sensitive to
Neptune’s period ratio with Uranus; if Neptune and Uranus are placed at their
2:1 resonance, there would be almost no stable libration in Neptune’s
resonances.

Our results imply that the drop-off in resonance sticking in the high-$a$
scattering TNO population (e.g., Lykawka & Mukai, 2007) is not a generic
outcome of small body scattering dynamics, but is instead a result of Uranus’s
perturbations on Neptune. The exact pattern of destabilization for Neptune’s
resonances as a result of Neptune-Uranus interactions is complex and has
important implications for the behavior of Neptune’s resonances during the era
of planet migration when interactions between Neptune and Uranus were
evolving. The overall number of stable distant resonant TNOs captured during
migration, the relative stable populations between the distant resonances, as
well as the number and distribution of resonant dropouts could be very
sensitive to the history of Neptune and Uranus’s period ratio. Future, large-
scale surveys such as the Vera Rubin Observatory’s Legacy Survey of Space and
Time (LSST) will dramatically improve observational constraints on the distant
resonant, scattering, and detached TNO populations (e.g. Ivezić & et al.,
2019). These observations should reveal complex patterns in the number of
temporarily resonant scattering TNOs that reflect the relative strengths of
the resonances. The revealed distributions of stably resonant and detached
TNOs will provide important constraints on the relative orbital histories of
Uranus and Neptune.

Acknowledgements: This work was supported by NASA (grant 80NSSC19K0785). KV
acknowledges additional support from NASA (grant 80NSSC21K0376) and NSF
(AST-1824869). This work used High Performance Computing (HPC) resources
supported by the University of Arizona TRIF, UITS, and Research, Innovation,
and Impact (RII) and maintained by the UArizona Research Technologies
department.

##  Appendix A Resonance strength calculation

As was explained in Section 3.3, we wanted to better quantify the strength of
resonances that accounts for test particles that ‘stick’ transiently in
addition to stably librating particles. We do this by considering total
resonance occupation timescales as a fraction of simulation time. For example,
a test particle may start as non-resonant, then move into Neptune’s 20:1
resonance, librate for 200 20:1 resonant cycles, then be scattered away. If
the simulation length is 1000 resonant cycles, then that test particle can be
said to be resonant 20% of the simulation time (with respect to the 20:1). For
this particular resonance, this one particle can be assigned a strength value
of 0.2. As discussed in Section 2.3, the nature of the sticking events and
resonance behavior we are interested in makes it difficult to rely solely on
the confinement of $\psi$ to determine periods of resonant interactions; the
symmetric libration zones of Neptune’s N:1 resonances in particular span
nearly the entire $\psi=0-360^{\circ}$ range, and even particles in the the
more confined resonant islands will experience slips in $\psi$ through the
full range. So instead we rely on relative confinement of semimajor axis to
determine periods of resonant interactions. We note that this method of
determining resonant interactions only works for the low-$q$ scattering
particles we are considering. High-$q$ particles can remain at constant $a$
over significant timescales without any resonant interactions, whereas low-$q$
particles are very mobile in $a$ if they are not interacting with resonances.
Even if low-$q$ particles stay at constant $a$ without true resonant
libration, they typically only do so because they are interacting with a
resonance boundary, so it is reasonable to include them in the strength
calculation.

We extend our analysis of individual test particle resonance interaction
timespans to assign an overall strength value to a given resonance based on
the entire set of test particles using the following steps:

  * •

Define a window in period ratio with Neptune around a target resonance. We
define the boundaries of N:1 resonances with a window size of $\pm 0.15$ in
period ration with Neptune on either side of the exact resonant ratio; for N:2
resonances we use a period ratio window of $\pm 0.1$, and $\pm 0.075$ for N:3
resonances. An example of this for the 20:1 resonance would be a period ratio
with Neptune from 19.85 to 20.15, corresponding to $a=220.36-222.6$. These
values adequately cover the expected width of each resonance type, without
overlapping onto neighboring strong resonances. We only limit our target
windows by period ratio, not $\psi$.

  * •

Identify every test particle that passes through this window throughout the
duration of the simulation and assign that particle a ‘residence time ratio’,
$s_{i}$, that describes the resonant interaction with a target resonance. This
is done as follows:

    * –

Define a sliding timescale window ($\Delta t$). We set this to length of time
that is 10% of the overall simulation time ($T=1000$ resonant cycles), so 100
resonant cycles ($\Delta t=100$). For each particle that passes though our
period ratio region of interest, we start the sliding time window at $t=0$ and
find the maximum and minimum period ratio recorded for that particle over
$\Delta t$. If these values are wholly within the period ratio bounds, then we
consider this particle ‘bound’ to that phase space window in that given
timespan of the simulation.

    * –

The start of the sliding time window is then set to a $t=1$ resonant cycle,
then $t=2$, then $t=3$, and so on until (for a simulation of 1000 total
resonant cycles) the window start reaches $t=900$, with the above
determination made for each window.

    * –

Once have recorded the number $N_{res}$ of those 900 time windows in which the
test particle was bound to the period ratio range, we calculate the fraction
of the total simulation time test particle $i$ was resonant:
$s_{i}=N_{res}/(T-\Delta t)$.

  * •

A unitless strength value for a particular resonance can then be determined by
the sum of residence time ratios ($s_{i}$) for all the particles that enter
that region of phase space throughout the duration of the simulation. This is
then normalized to the total number of test particles that are initialized in
that region at the beginning of the simulation time ($n_{tp,0}$). The strength
value $S$ is thus:

$S=\frac{\sum_{0}^{N}{s_{i}}}{n_{tp,0}}.$ (A1)

We note that some of the particles whose $s_{i}$ values are included in $S$
were initialized outside of the period ratio range from which $n_{tp,0}$ is
determined. But because we are considering a low-$q$ scattering population,
and many of the particle initialized near a resonance will scatter fairly
quickly, the value of $S$ in our simulations is always between 0 and 1, and
$S$ is approximately the average of all residence time ratios of for particles
in the target period ratio window. This normalization scheme is tailored to
our scheme of initializing particles along lines of constant $\psi$ and
equally spaced in the semimajor axis; other normalizations would need to be
considered for alternative test particle distributions.

As noted above, our strength parameter $S$ is best suited for the low-$q$
scattering population we focus on. It is also best used as a relative rather
than absolute strength measurement; this is why we typically normalize
resonance strengths to the value of $S$ for the strongest resonance in a
simulation (e.g., Figure 6).

When plotting $a$-$\psi$ phase diagrams in Figure 7, we use a modified version
of a particle residence time to better highlight general particle stability.
We used the same sliding window method as described above for examining
particle behavior for the duration of the simulation, but looked only at
semimajor axis stability in each window rather than confinement to a specific
period ratio range. If the semimajor axis of a particle changed by no more
than 2 au (roughly twice the width of typical resonances in the range we are
interested in) we considered it stable for that time window. The fraction of
all examined time windows for which a particle was stable is then the
‘normalized residence time’ parameter used to color-code particles in the
$a$-$\psi$ maps in Figure 7. We found this scheme the most effective for
visually illustrating the changes in resonant behavior in our different
simulations.

## References

  * Batygin et al. (2021) Batygin, K., Mardling, R. A., & Nesvorný, D. 2021, ApJ, 920, 148, doi: 10.3847/1538-4357/ac19a4
  * Crompvoets et al. (2022) Crompvoets, B. L., Lawler, S. M., Volk, K., et al. 2022, , 3, 113, doi: 10.3847/PSJ/ac67e0
  * Gallardo (2006) Gallardo, T. 2006, Icarus, 181, 205, doi: 10.1016/j.icarus.2005.11.011
  * Gallardo (2018) —. 2018, Planet. Space Sci., 157, 96, doi: 10.1016/j.pss.2018.03.007
  * Gallardo (2020) —. 2020, Celestial Mechanics and Dynamical Astronomy, 132, 9, doi: 10.1007/s10569-019-9948-7
  * Gladman et al. (2008) Gladman, B., Marsden, B. G., & Vanlaerhoven, C. 2008, in The Solar System Beyond Neptune, ed. M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, A. Morbidelli, & R. Dotson (Univ. Arizona Press), 43–57
  * Gladman & Volk (2021) Gladman, B., & Volk, K. 2021, ARA&A, 59, 203, doi: 10.1146/annurev-astro-120920-010005
  * Gomes et al. (2005) Gomes, R. S., Gallardo, T., Fernández, J. A., & Brunini, A. 2005, Celestial Mechanics and Dynamical Astronomy, 91, 109, doi: 10.1007/s10569-004-4623-y
  * Hadden & Tremaine (2024) Hadden, S., & Tremaine, S. 2024, MNRAS, 527, 3054, doi: 10.1093/mnras/stad3478
  * Hermosillo Ruiz et al. (in prep) Hermosillo Ruiz, A., Volk, K., Murray-Clay, R., & Pike, R. in prep
  * Ivezić & et al. (2019) Ivezić, Ž., & et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Kaib & Sheppard (2016) Kaib, N. A., & Sheppard, S. S. 2016, AJ, 152, 133, doi: 10.3847/0004-6256/152/5/133
  * Lan & Malhotra (2019) Lan, L., & Malhotra, R. 2019, Celestial Mechanics and Dynamical Astronomy, 131, 39, doi: 10.1007/s10569-019-9917-1
  * Lawler et al. (2019) Lawler, S. M., Pike, R. E., Kaib, N., et al. 2019, AJ, 157, 253, doi: 10.3847/1538-3881/ab1c4c
  * Lykawka & Mukai (2007) Lykawka, P. S., & Mukai, T. 2007, Icarus, 192, 238, doi: 10.1016/j.icarus.2007.06.007
  * Malhotra (2019) Malhotra, R. 2019, Geoscience Letters, 6, 12, doi: 10.1186/s40562-019-0142-2
  * Malhotra & Ito (2022) Malhotra, R., & Ito, T. 2022, Proceedings of the National Academy of Science, 119, 2118692119, doi: 10.1073/pnas.2118692119
  * Malhotra et al. (2018) Malhotra, R., Lan, L., Volk, K., & Wang, X. 2018, AJ, 156, 55, doi: 10.3847/1538-3881/aac9c3
  * Nesvorný et al. (2016) Nesvorný, D., Vokrouhlický, D., & Roig, F. 2016, ApJ, 827, L35, doi: 10.3847/2041-8205/827/2/L35
  * Pan & Sari (2004) Pan, M., & Sari, R. 2004, AJ, 128, 1418, doi: 10.1086/423214
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rein & Spiegel (2015) Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424, doi: 10.1093/mnras/stu2164
  * Rein & Tamayo (2015) Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376, doi: 10.1093/mnras/stv1257
  * Saillenfest (2020) Saillenfest, M. 2020, Celestial Mechanics and Dynamical Astronomy, 132, 12, doi: 10.1007/s10569-020-9954-9
  * Smirnov & Markov (2017) Smirnov, E. A., & Markov, A. B. 2017, MNRAS, 469, 2024, doi: 10.1093/mnras/stx999
  * Tiscareno & Malhotra (2003) Tiscareno, M. S., & Malhotra, R. 2003, AJ, 126, 3122, doi: 10.1086/379554
  * Volk & Malhotra (2022) Volk, K., & Malhotra, R. 2022, ApJ, 937, 119, doi: 10.3847/1538-4357/ac866b
  * Wang & Malhotra (2017) Wang, X., & Malhotra, R. 2017, AJ, 154, 20, doi: 10.3847/1538-3881/aa762b
  * Wolff et al. (2012) Wolff, S., Dawson, R. I., & Murray-Clay, R. A. 2012, ApJ, 746, 171, doi: 10.1088/0004-637X/746/2/171
  * Yu et al. (2018) Yu, T. Y. M., Murray-Clay, R., & Volk, K. 2018, AJ, 156, 33, doi: 10.3847/1538-3881/aac6cd