Neptune’s resonances in the Scattered Disk
The Scattered Disk Objects (SDOs) are thought to be a small fraction of the ancient population of leftover planetesimals in the outer solar system that were gravitationally scattered by the giant planets and have managed to survive primarily by capture and sticking in Neptune’s exterior mean motion resonances (MMRs). In order to advance understanding of the role of MMRs in the dynamics of the SDOs, we investigate the phase space structure of a large number of Neptune’s MMRs in the semi-major axis range 33–140 au by use of Poincaré sections of the circular planar restricted three body model for the full range of particle eccentricity pertinent to SDOs. We find that, for eccentricities corresponding to perihelion distances near Neptune’s orbit, distant MMRs have stable widths that are surprisingly large and of similar size to those of the closer-in MMRs. We identify a phase-shifted second resonance zone that exists in the phase space at planet-crossing eccentricities but not at lower eccentricities; this second resonance zone plays an important role in the dynamics of SDOs in lengthening their dynamical lifetimes. Our non-perturbative measurements of the sizes of the stable resonance zones provide an explanation for the prominence of the :1 sequence of MMRs over the :2, :3 sequences and other MMRs in the population statistics of SDOs, yield a theoretical understanding of the outer boundary of SDOs at perihelion distance near 36 au, and also provide a tool to more easily identify resonant objects.
The population of near- and beyond-Neptune minor planets overed over the past two-and-a-half decades provides important diagnostics about the origin and evolution of our solar system (Luu & Jewitt, 2002; Morbidelli & Brown, 2004). A number of these objects are recognized to be part of an extensive “scattered” and “scattering” population that has incurred gravitational scattering encounters with Neptune over its past history; this population inhabits a region of orbital parameter space defined by somewhat fuzzy boundaries of semimajor axis in the range and perihelion distance in the range 37 au (Gladman et al., 2008). Theoretical models suggest that the present-day Scattered Disk Objects (SDOs) are just a small, , fraction of the ancient population of leftover planetesimals formed in the outer solar system that were gravitationally scattered by the giant planets (Levison & Duncan, 1997; Duncan & Levison, 1997; Lykawka & Mukai, 2007b; Gomes et al., 2008). The SDOs are postulated to be a major source of the transient populations of Centaurs and the Jupiter-family comets and possibly also contribute to the Halley-type comets and the Oort Cloud (Duncan & Levison, 1997; Fernández et al., 2004; Emel’yanenko et al., 2005; Levison et al., 2006; Volk & Malhotra, 2008).
The median dynamical lifetime of a minor planet subsequent to a strong gravitational scattering encounter with Neptune is about ten million years (Levison & Duncan, 1997; Tiscareno & Malhotra, 2003; Di Sisto & Brunini, 2007). Typically, such encounters result in further strong scattering by the other giant planets, leading to either ejection of the small body from the solar system or to collisions with the planets or the Sun. The survival and persistence of a remnant of the ancient planetesimal disk in the form of the present-day Scattered Disk population over the age of the solar system is currently understood to be owed to a phenomenon called “resonance sticking” which ensures that when these objects approach Neptune’s orbit near their perihelion their true longitude is usually well-separated from Neptune’s (Gomes et al., 2008). This is not dissimilar to the well-understood mechanism for the long term stability of Pluto’s Neptune-crossing orbit (Milani et al., 1989; Malhotra & Williams, 1997). However, unlike the stable libration of Pluto’s perihelion well away from Neptune’s true longitude, the SDOs’ orbits evolve near the chaotic boundaries of MMRs, such that they are not strictly prohibited from close encounters with Neptune. Numerical models find that the SDOs have only occasional gravitational scattering encounters with Neptune but otherwise spend a large fraction of their dynamical lifetime in the vicinity of mean motion resonances with that planet (Duncan & Levison, 1997; Lykawka & Mukai, 2007a; Volk et al., 2018; Yu et al., 2018). In this “resonance sticking” behavior, the SDOs visit many different MMRs for varying (random) lengths of time. The simulations indicate that, in order of decreasing importance, the :1, :2 and :3 resonances play the most important roles in the resonance sticking phenomenon, in terms of both capture probabilities and the residence time in the MMRs. This pattern for the SDOs’ resonance sticking is consistent with a prior theoretical study by Pan & Sari (2004) who showed that, in the simple model of the circular planar restricted three body problem, the :1 exterior resonances have the largest resonance widths when the Jacobi constant is close to 3; the :2, :3, etc. resonances have decreasing size of resonance widths when the particle eccentricity is in the scattering regime. However, it is in contrast with the classical perturbation theoretic analyses of mean motion resonances at lower eccentricities which anticipates that resonance width decreases exponentially with the order of a resonance (Murray & Dermott, 1999). In the high eccentricity (scattering) regime, Pan & Sari (2004) suggested a redefinition of the “order” of resonance as “”-th order for a resonance.
The “resonance sticking” mechanism is complex and we still have only limited understanding of it. A number of previous studies have examined the dynamics of Neptune’s MMRs with semi-analytical and numerical models in various approximations (e.g. Beauge, 1994; Morbidelli et al., 1995; Malhotra, 1996; Morbidelli, 1997; Nesvorný & Roig, 2000, 2001; Robutel & Laskar, 2001; Chiang & Jordan, 2002; Kotoulas & Voyatzis, 2004a, b; Kotoulas, 2005; Voyatzis & Kotoulas, 2005; Murray-Clay & Chiang, 2005; Tiscareno & Malhotra, 2009; Yu et al., 2018; Gallardo, 2018, 2019). The present paper provides another view on Neptune’s exterior resonances pertinent to understanding their role in the resonance sticking dynamics of the SDOs. We begin by recognizing that, in the high eccentricity regime of the SDOs, the chaotic separatrices at the boundaries of the stable resonance zones are undoubtedly important in the “resonance sticking” dynamics but their role in defining the widths of resonances is not accurately computed in analyses that isolate a single resonance and/or average over short period perturbations. The nature of the phase space is revealed more clearly with non-perturbative numerical analyses of the behavior of planet-crossing orbits. We therefore investigate the phase space structure of a wide range of Neptune’s MMRs by examination of Poincaré surfaces of section in the circular planar restricted three body model of the Sun, Neptune and a test particle (the latter representing an SDO). The advantage of this simple model is that its two dimensional Poincaré sections serve to accurately compute and visualize the chaotic boundaries of resonances and provide a non-perturbative measure of the widths of the stable libration zones near resonances for the full range of particle eccentricity. (A disadvantage is that non-co-planarity, the effects of a non-circular Neptune orbit, and the perturbations of other planets are not included; still, the map of stable libration zones of moderate-to-high-eccentricity MMRs in the simplified model yields new insights, as we will show.) In the context of the outer solar system, this non-perturbative method has been previously employed by Pan & Sari (2004) for their study of exterior resonant orbits of Jacobi constant close to 3 (pericenter distance close to the planet’s orbit); Wang & Malhotra (2017) employed it for a general study of the 2:1 and 3:2 interior resonances in the high eccentricity regime. Malhotra (1996) and Malhotra et al. (2018) made use of this method for a small number of Neptune’s exterior MMRs. Here we use this method for a large number of Neptune’s exterior MMRs pertinent to the dynamics of SDOs, in the semi-major axis range au (particle-to-Neptune period ratio from 6:5 to 10:1 and test particle eccentricity from to ). We carry out comparisons of different resonances to identify patterns that can help towards a better understanding of the resonance sticking phenomenon. Additionally, because our investigation yields fairly accurate boundaries of the libration zones of MMRs in the semi-major axis–eccentricity plane, a by-product of our investigation is a new tool to identify resonant minor planets beyond Neptune; this tool is complementary to the current method of arduous examination of the behavior of critical resonant angles in megayear long numerical integrations (e.g., Gladman et al., 2008; Volk et al., 2016).
2 Poincaré sections for Neptune’s exterior resonances
To study the phase space structure and dynamics near Neptune’s mean motion resonances, we generated Poincaré sections of the circular planar restricted three body model of the Sun, Neptune and a massless test particle (the latter representing a Kuiper belt object). Because Neptune’s orbit is nearly circular (its eccentricity does not exceed over long timescales (Murray & Dermott, 1999)), and because the masses of the Kuiper belt objects are very small ( of Neptune’s), this simplified model is sufficient to identify the basic phase structure and dynamics near Neptune’s exterior mean motion resonances. This model neglects the influence of Neptune’s eccentricity and the solar system planets interior to Neptune, as well as non-coplanar motion, but affords the advantage of visualization of the phase space structure in two-dimensional sections. These neglected effects are not critical for the present purpose. In a study of Neptune’s 5:2 MMR, Malhotra et al. (2018) showed that perturbations of the other giant planets (Jupiter, Saturn and Uranus) and non-co-planarity exert only a mild influence on the boundaries of the resonance zone in the plane.
Natural units are adopted for this model: we set the unit of mass to be (where and denote the masses of Sun and Neptune, respectively), the unit of length to be the constant orbital separation of the two primaries (i.e., Neptune’s semi-major axis, ), and the unit of time to be their orbital period divided by . With these units, Newton’s constant of gravitation and the orbital angular velocity of the two primaries about their barycenter are unity. We employed the equations of motion for the test particle in the rotating reference frame of constant unit angular velocity with origin at the barycenter of and ; in this frame, the two primaries are fixed at locations, and , respectively, where is the fractional mass of Neptune. For the numerical integrations of the equations of motion we used the adaptive step size seventh-order Runge-Kutta method (Fehlberg, 1968), with relative and absolute error tolerances of .
We adopted a systematic strategy for the exploration of near-resonant particle initial conditions as in Malhotra et al. (2018). Briefly, all the test particle trajectories represented in a Poincarè section have the same value of the Jacobi integral and the same initial values of semi-major axis and eccentricity, but different values of the initial longitude of perihelion. The particle’s initial osculating semi-major axis is always set equal to the nominal resonant value, , where and are mutually prime integers, and is the period ratio of test particle and Neptune. The initial location of each particle in its orbit is always at its perihelion but different particles have different initial pericenter longitudes spanning the full range ). We integrate the particle motion for several thousand orbital periods, recording its state vector () at every perihelion passage. We then transform the state vector into osculating orbital elements and generate the Poincaré section as a plot of (), where is the osculating barycentric semi-major axis, and is the true longitude separation between Neptune and the test particle when the particle is at perihelion. The phase angle is related to the usual critical resonant angle for an : resonance as follows:
Figures 2 and 9 show examples of Poincaré sections for Neptune’s 2:1 and 3:2 exterior MMRs. In these plots, we see stable resonant orbits that make closed smooth curves and we also see chaotic orbits that wander over a wider area and do not stay confined to smooth curves. The region containing a sequence of closed smooth curves is a stable resonance island. The center of a stable island indicates the exact resonant orbit (corresponding to a periodic orbit), whereas the smooth curves surrounding it are orbits librating about the exact resonance (corresponding to quasi-periodic orbits). The centers of the stable resonant islands are located at a few rete values of , each with a semi-major axis value close to . Each Poincaré section is labeled with the value of the initial eccentricity common to all the trajectories represented in that section; this value is very close to the eccentricity of the periodic orbit at the center of the stable islands in that section. For each initial eccentricity value, we measure the width of the stable island by its range in semi-major axis, , where and are the upper and lower semi-major axis boundaries of the stable island. For most MMRs, the chaotic orbits are not visible at low eccentricities, but become visible at moderate and high eccentricity, as scattered points filling an area bounding the stable islands. As illustrated in Figures 2 and 9, the phase space structure changes with eccentricity, with at least one particularly significant transition when new stable islands appear when eccentricity exceeds the Neptune-crossing value. We show in Sections 4 and 5 that these phase space transitions are related to the shape of the resonant orbit in the rotating frame and its relationship to Neptune’s orbit (see also Wang & Malhotra (2017) and Malhotra et al. (2018)).
After generating many Poincaré sections and measuring the stable resonance boundaries, and , for the full range of eccentricities for many MMRs, we gather the results and plot the resonance boundaries in the parameter plane. Figures 1, 6 and 12 show the boundaries of resonances with the particle-to-Neptune period ratios from 6:5 to 10:1. Each MMR has two distinct libration zones, indicated as the “first resonance zone” (black curves) and the “second resonance zone” (red curves). These two zones have stable islands that are centered at different values of the angular separation, , of the particle’s perihelion longitude relative to Neptune. The stable island centers of the first resonance zone avoid , whereas the second resonance zone has a stable island centered at . Physically this means that the periodic orbits at the center of each of these zones correspond to traces in the rotating frame that have the same shape but different geometrical orientations relative to the Sun-Neptune line. In terms of the critical resonant angle, (Eq. 1), the first resonance zone has librations of centered at whereas the second resonance zone has librations of centered at 0.
For all the MMRs investigated here, the first resonance zone exists over almost the entire range of eccentricity; for some resonances, it disappears for a small range of high eccentricities, then reappears again at an even higher eccentricity. The second resonance zone emerges at an eccentricity slightly exceeding the critical value, , for Neptune-crossing,
For the :1 resonances, the stable islands of the second resonance zone expand monotonically with increasing eccentricity, but this is not the case for the :2, the :3 and higher resonances; for latter, the second resonance zone has several transitions with increasing eccentricity, first expanding then shrinking, even disappearing (e.g., the 3:2, 5:3 and 8:5 resonances).
3 N:1 resonances
We computed the Poincaré sections for the :1 sequence of Neptune’s exterior MMRs, from 2:1 to 10:1, for the full range of particle eccentricities. The 2:1 MMR is the closest-in of this sequence, and its phase space structure is exemplary of all the :1 sequence of MMRs. Below we describe this MMR in some detail.
3.1 Neptune’s 2:1 Exterior MMR
Figure 2 plots several representative surfaces of section in the neighborhood of Neptune’s 2:1 exterior MMR (semi-major axes around ). We observe that at the lower value of eccentricity, near , only smooth closed curves are present, and there is no visible chaotic region. The libration zone has semi-major axis centered at , and it contains three families of orbits, each of which has different centers for : a pair of families confined to small libration amplitudes has librating about a center near and ( is eccentricity-dependent and varies from near to near – see Fig. 5), and the third family has larger amplitude librations of centered at . The boundary, or separatrix, between the latter and the former is a curve shaped like the outline of a butterfly. The latter family is called the “symmetric” librators, while the former pair are called the “asymmetric” librators (Beauge, 1994; Malhotra, 1996; Nesvorný & Roig, 2001; Chiang & Jordan, 2002; Pan & Sari, 2004). This nomenclature is also descriptive of the center of librations of the critical resonant angle, . Physically, it is descriptive of the orientation of the particle’s perihelion relative to Neptune (see Fig. 3a): for the symmetric libration center, the particle’s perihelion longitude is away from Neptune and the trace of its orbit in the rotating frame is symmetric relative to the line joining the fixed locations of the Sun and Neptune; for the asymmetric libration centers, the angular separation of the particle’s perihelion from Neptune is approximately (for particle eccentricity ) and the trace of its orbit in the rotating frame is not symmetric relative to the line joining the fixed locations of the Sun and Neptune. As evident in the surface of section, Figure 2a, the two families of asymmetric librators are confined to small libration amplitudes; the symmetric librators surround the two zones of the asymmetric librators and are accordingly confined to large libration amplitudes. We report the width in semi-major axis, , of the first resonance zone as measured by the boundaries of the symmetric libration zone; it includes within it the asymmetric libration zone.
At larger eccentricities, the test particle’s perihelion is closer to Neptune’s orbit. The Poincaré sections show that the stable island is larger in semi-major axis width, . The center of the asymmetric librations drifts further away from , and their maximum libration amplitude also increases, squeezing the libration amplitudes of the symmetric librators to a smaller range. As seen in Figure 2(b), for , the orbits near the outer boundary of the libration zone are visibly chaotic. When the eccentricity exceeds the critical Neptune-crossing value, , the particle’s perihelion is interior to Neptune’s orbit. At eccentricity , we find that the width, , of the libration zone reaches a maximum. At a somewhat higher eccentricity, , a new stable island, which is centered at , becomes visible in the Poincaré section; this can be seen in Figure 2(c) for . We call this the second resonance zone. For increasing eccentricity above , the stable island of the second resonance zone expands at the expense of the size of the islands of the first resonance zone (Figure 2(d,e,f)). The emergence of a new stable island signals the existence of a pair of new periodic orbits in the 2:1 exterior resonance. The trace of the periodic orbits in the rotating frame has the same shape for the new and old periodic orbits, but their orientation relative to the fixed locations of the Sun and Neptune in the rotating fame are actually different. This is illustrated in Figure 3 which shows the paths of the particle’s 2:1 resonant orbit in the rotating frame, at different particle eccentricities. We observe that, for , the resonant orbit intersects itself and creates a “perihelion lobe”. The circle of unit radius about the Sun (which is the circular orbit of Neptune) cuts the perihelion lobe of the particle’s trajectory. For the second resonance zone, the fixed location of Neptune in the rotating frame is inside the perihelion lobe, whereas for the resonant orbits of the first resonance zone, Neptune’s location is outside the perihelion lobe. We observe that the length of arc of the unit circle which is interior to the perihelion lobe correlates with the range of librations of in the second resonance zone: the longer the length of this arc, the wider the stable islands of the second resonance zone. At an eccentricity of , the widths, in semi-major axis, of the new and old stable islands are nearly equal. For eccentricities above , the semi-major axis width of the new stable island exceeds that of the old one, and the sizes of the libration islands of the first resonance zone shrink rapidly. The two asymmetric islands of the first resonance zone merge into a single, symmetric island when the eccentricity exceeds .
As noted above, the centers of the two asymmetric islands drift away from as the eccentricity increases from small values to , and then drift back toward for eccentricity increasing from to . This behavior is shown in Figure 5. Nesvorný & Roig (2001) have previously reported the locations of the asymmetric centers of the 2:1 resonance for eccentricity values up to 0.7; our results reported here are similar for that case but extend to eccentricity values up to unity. We additionally report the locations of the asymmetric centers of the 3:1 and 4:1 resonances. For increasing of the :1 resonances, it is interesting to note that the asymmetric islands first become visible at a higher threshold eccentricity.
A dynamical quantity of interest is the libration period of resonant orbits. A librating orbit is not at exact resonance so it is not a periodic orbit; it can be thought of as tracing almost the same shape as the exact resonant orbit, but slowly drifting and librating about the exact periodic orbit. We measured the libration periods by measuring the time for successive points in the Poincaré section to complete one circuit around the center of a resonant island. For Neptune’s 2:1 exterior resonance, we note that each smooth closed curve in a Poincaré section is uniquely identified by the minimum value, , of the critical resonant angle on a librating orbit. These libration periods are plotted as a function of libration amplitude of the resonant angle, , in Figure 4. In the first resonance zone, the libration period is not a monotonic function of eccentricity, nor of the libration amplitude. As a function of eccentricity, the libration period decreases as eccentricity increases in the range 0.1 to 0.7, but the trend reverses for higher eccentricities. The shortest libration periods occur for librating orbits with close to ; the libration period increases very slowly with increasing when this amplitude is small; indeed, no variation is measurable for . At higher , the libration period rises sharply and presumably diverges at the boundary between the symmetric and asymmetric families of librating orbits. In the zone of the symmetric librators, the libration period drops sharply as we move away from its inner boundary and towards larger amplitude librations. In the region of the sharp spike of the libration period, the numerical solutions show that the test particle orbit librates very slowly in the narrow boundary region which connects the two wings of the asymmetric librators, leading to the quite long libration periods. We conjecture that this is one of the likely origins of “resonance sticking” behavior of SDOs.
In the second resonance zone, the libration periods are shorter than in the first resonance zone (for similar values of eccentricity), and they decrease slightly with increasing libration amplitude.
3.2 Patterns in the :1 sequence of exterior MMRs
The phase space structure of all the :1 exterior resonances is qualitatively similar to that of the 2:1 resonance. The width, , of the first resonance zone (measured by the outer boundaries of symmetric librators) increases with eccentricity from , reaches a maximum near (Eq. 2, corresponding to perihelion distance equal to Neptune’s orbital radius, au), then slowly shrinks with increasing particle eccentricity beyond . For larger , the asymmetric libration zone first appears at larger values of eccentricity. The centers of the asymmetric librations as a function of particle eccentricity are shown in Figure 5 for the 2:1, 3:1 and 4:1 MMRs. Similar to the 2:1 MMR, these libration centers drift away swiftly from near towards an extremum near , then drift back towards as the eccentricity approaches unity. The second resonance zone first emerges at eccentricity (corresponding to perihelion distance au), and expands in width as the eccentricity approaches unity.
The boundaries in the plane of the first and second resonance zones are plotted in Figure 6 for the sequence of :1 resonances, ; in this figure we also plot the boundaries of the asymmetric libration zone (which is a subset of the first resonance zone). It is apparent that the resonance widths in semi-major axis, , depend on both and eccentricity. This dependence is illustrated further in Figure 7. Perhaps the most notable result is that the maximum width of the first resonance zone (which occurs near eccentricity ) is rather similar for all ; it changes only very slowly with increasing period ratio (see Figure 7a). The 3:1 resonance has the largest width, even larger than 2:1 resonance, but all the :1 resonances we examined (up to ) have quite large maximum widths. In physical units, these maximum widths are all near au. For eccentricity near , with increasing the size of the asymmetric libration zone is an increasing fraction of the size of the first resonance zone.
For the 3:1 resonance, the first and second resonance zones have similar maximum widths, but for the larger particle-to-Neptune period ratios, the maximum width of the second resonance zone significantly exceeds that of the first resonance zone (Figure 7a). In contrast with the behavior near the planet-crossing values of eccentricities, the first resonance zone’s widths at lower eccentricities decrease rapidly with the increasing period ratio, as shown in Figure 7b. For period ratios larger than 6:1 and eccentricities up to , the width of the first resonance zone does not exceed 0.1 au.
We also observe in Figure 7a that the maximum width of the second resonance zone increases rapidly with increasing period ratio. For the 2:1 resonance, the first and second resonance zones have similar maximum widths, but for the larger particle-to-Neptune period ratios, the maximum width of the second resonance zone significantly exceeds that of the first resonance zone. The best-fit power law, , provides a very good empirical approximation for the width of the second resonance zone for eccentricities near the critical Neptune-crossing value.
The behavior of the resonance libration periods with libration amplitude and with eccentricity is qualitatively similar to that of the 2:1 MMR (Figure 4). The small amplitude libration periods of the asymmetric librations generally increase with and decrease with increasing eccentricity, as shown in Figure 8. Near the eccentricity , the small amplitude libration period can be approximated with a double power law: , for :1 resonances of ; for larger , we find . At other values of eccentricity, the small amplitude libration period increases exponentially with (Figure 8b).
4 N:2 resonances
We computed the Poincaré sections for the :2 sequence of Neptune’s exterior MMRs, from 3:2 to 19:2, for the full range of particle eccentricities. The 3:2 MMR is the closest-in of this sequence, and its phase space structure is exemplary of all the :2 sequence of MMRs. We first describe this resonance in some detail, and then discuss the patterns of behavior of the whole sequence.
4.1 The 3:2 MMR
Examining the Poincaré sections near the 3:2 MMR, we observe some similarities and some differences with the 2:1 MMR. As seen in Figure 9(a), at low eccentricity, there is no visible chaotic region in the Poincaré section in () which has a pair of stable islands composed of smooth closed curves. This is the first resonance zone. The pair of islands is centered at the two values of and , but share a common central value of the semi-major axis, . Unlike in the case of the 2:1 MMR, these libration centers of do not drift with the change of eccentricity. Also, unlike the case of the 2:1 MMR, these two islands are not separate families but form a chain of two islands, such that a single trajectory visits alternately each island in the Poincaré section. That is to say, the exact 3:2 resonant orbit is a “period two” orbit. In the inertial frame, this orbit completes two circuits around the Sun in the time it closes one complete trace in the rotating frame. A pair of unstable points exists at the separatrix of this chain of two islands. These two unstable points are located at and , but at slightly different values of the semi-major axis; the slight difference in the semi-major axis is due to the different positions of Neptune at the two successive perihelion passages of the test particle in this configuration. This difference gives rise to a slight difference in perihelion distance and perihelion velocity, hence a slight difference in the osculating orbital parameters at alternate perihelion passages.
The widths of the stable islands of the 3:2 MMR increase with eccentricity, reaching a maximum at eccentricity , where the perihelion distance is au. At an eccentricity , where the perihelion distance approaches 27 au, a new pair of stable islands appear in the surface of section (see Figure 9(b)). These are centered at and , in-between the old pair of stable islands observed at lower eccentricity. We call this pair the second resonance zone. (As with the old pair, this pair is also a period-two chain of islands.) With increasing eccentricity, the semi-major axis widths, , of these new stable islands expand at the expense of the old islands, reaching a maximum at eccentricity . In this regime, the width, , of the first resonance zone decreases with increasing eccentricity and vanishes at (Figure 9(c,d)). At an eccentricity , which is slightly larger than , the stable islands of the first resonance reappear, and their sizes increase with increasing eccentricity. Meanwhile, the width of the second resonance zone shrinks slightly (Figure 9(e)). Thus, when the eccentricity approaches unity, the width of the second resonance zone is only slightly larger than that of the first resonance zone (Figure 9(f)).
The transitions with increasing eccentricity in the Poincaré sections of Neptune’s 3:2 MMR described above are correlated with the changes in the geometry of the trace of the resonant orbit in the rotating frame. In Figure 10, we show the trajectories of the particle’s 3:2 resonant orbit in the rotating frame at different particle eccentricities. We observe that the first appearance of the second resonance zone occurs when the eccentricity slightly exceeds the Neptune-crossing eccentricity, and the small initial size of this zone is due to the small length of the arc of the Sun-centered unit circle which is cut by the particle’s perihelion lobes. With increasing eccentricity, the length of arc enclosed by the perihelion lobes increases, and the size of the second resonance zone increases correspondingly, at the expense of the first resonance zone. The vanishing of the first resonance zone occurs when the two perihelion lobes touch each other, when the eccentricity approaches . At even higher eccentricity, the two perihelion lobes intersect each other, creating new arcs for the reappearance of the first resonance zone.
The behavior of the libration period of the resonant angle, , is rather complicated, as shown in Figure 11. In the first resonance zone, for the eccentricity range 0.1–0.7, the libration period decreases with increasing eccentricity and also decreases with increasing libration amplitude, ; however, the librating orbits of have longer libration periods than those of . In the second resonance zone, the libration period decreases with increasing libration amplitude, (similar to the behavior in the first resonance zone), but increases with increasing eccentricity (unlike in the first resonance zone).
4.2 Patterns in the :2 sequence of exterior MMRs
By examining the Poincaré sections for the full range of eccentricities for the :2 sequence of Neptune’s exterior MMRs, 3:2 to 19:2, we ascertained the boundaries of the stable libration zones. These are shown in Figure 12 in the () plane. Similar to the 3:2 resonance, every :2 resonance has a first resonance zone and a second resonance zone, each with a pair of stable islands. With increasing values of the test particle eccentricity, the resonance structure in the Poincaré sections show five transitions at eccentricity values which are similar to those described above for the 3:2 MMR. The first resonance zone consists of a pair of stable islands centered at and . Starting at a low eccentricity, we observe that its width, , increases with increasing value of eccentricity, its largest extent occurs at an eccentricity, (corresponding to perihelion distance au). The stable islands then shrink and finally vanish at an eccentricity . Meanwhile, at an eccentricity, (corresponding to perihelion distance au), the second resonance zone emerges, consisting of a new pair of stable islands centered at and (located in-between the stable islands of the first resonance). The second resonance zone expands with increasing eccentricity from , reaches a maximum size at eccentricity , and then shrinks with increasing eccentricity. At an eccentricity, , the pair of stable islands of the first resonance zone (centered at their original locations, and ) emerges again. These expand in width, while the islands of the second resonance zone shrink as the eccentricity approaches unity. A detail that we note in Figure 12 is that the boundaries of the second resonance zones are not symmetric about . This is due to the different positions of Neptune at the two successive perihelion passages of the test particle, as noted above (Section 4.1).
The variation of the resonance width with particle-to-Neptune period ratio is illustrated in Figure 13. We plot the maximum widths (which occur at , perihelion distance au) of the first and second resonance zones in Figure 13a. Notably, the maximum width of the first resonance zone is almost the same for all :2 resonances; it decreases very slowly with increasing particle-to-Neptune period ratio. Even more surprisingly, the maximum width of the second resonance zone increases markedly with increasing particle-to-Neptune period ratio. For the 9:2 resonance, the first and the second resonance zones have almost the same maximum widths. For the higher period ratios, the maximum width of the second resonance zone even exceeds that of the first resonance zone.
In Figure 13b, we plot the :2 resonance widths at lower and moderate eccentricities ( and ). In this eccentricity regime, the resonance widths drop quickly with increasing particle-to-Neptune period ratio; only those :2 resonances with period ratio below about 6 have widths exceeding 0.1 au.
The variation of the libration period with libration amplitude of the :2 resonances is qualitatively similar to that of the 3:2 MMR (Figure 11), but with an overall scale factor that can be deduced from the small amplitude libration periods (); the latter are plotted in Figure 14 as a function of the particle-to-Neptune period ratio, for several values of eccentricity. We observe that, in the first resonance zone, the small amplitude libration period, , at eccentricity (where the width of the first resonance zone reaches a maximum) approximately follows a power-law relation with particle-to-Neptune period ratio: ( denotes period ratio, ). However, for eccentricity away from the small amplitude libration period has an exponential relation with particle-to-Neptune period ratio (Figure 14b); the best-fit functions are found to be for , and for . For the second resonance zone, the small amplitude libration period at (where the width of the second resonance reaches a maximum) also is approximately a power-law relation with particle-to-Neptune period ratio: .
5 Summary and Discussion
We investigated the phase space structure of many of Neptune’s exterior MMRs over nearly the entire range of eccentricities, using non-perturbative numerical analysis and visualizations with Poincaré sections of the circular planar restricted three body model of the Sun, Neptune and test particle. We presented the Poincaré sections in the plane, where is the longitude separation of the planet from the test particle and is the test particle’s osculating semi-major axis at the particle’s perihelion passage; this allows easy identification of the resonance libration zones and measurement of their widths in semi-major axis. Our investigation spanned all the :1 and :2 sequence of resonances in the semi-major axis range 33–140 au (particle-Neptune period ratios up to 10), and a few :3, :4 and :5 resonances. We find that, in general, there are two distinct resonance zones in phase space. The first resonance zone exists over almost the entire range of eccentricities but the second resonance zone exists only in the high eccentricity regime for eccentricities exceeding the planet-crossing value. The first and second resonance zones differ from each other in that their libration centers are shifted in phase. In the first resonance zone, conjunctions with the planet are stable near the aphelion of the particle, whereas in the second resonance zone the conjunctions are stable near the particle’s perihelion. In terms of the usual resonant angle, Eq. 1, the first resonance zone librations are centered at , whereas the second resonance zone librations are centered at . (For the :1 resonances, the first resonance zone is further split into symmetric and asymmetric librators; see below.) In the high eccentricity planet-crossing regime, a particle in the first resonance zone avoids close encounters with the planet by reaching perihelion at true longitudes well separated from the planet’s true longitude. This is the well known “phase protection” mechanism for many resonant populations in the Kuiper belt, such as the Plutinos and Twotinos (Malhotra, 1995). In the second resonance zone, the particle avoids encounters with the planet differently: by “looping around the planet” during the part of its orbit that is interior to the planet’s; this physical explanation is illustrated by the trace of the high eccentricity resonant orbit in the rotating frame (for examples, see Figures 3 and 10).
For the first resonance zone, the width in semi-major axis, , is quite small at low eccentricity, rises sharply with increasing eccentricity, reaches a maximum near (or slightly below) the critical planet-crossing eccentricity (Eq. 2), then decreases again; for closer-in resonances, it vanishes and then reappears again at higher eccentricities whereas for the more distant resonances its width monotonically decreases as eccentricity approaches unity (Figures 1, 6, 12). Previous studies based on perturbative approaches with a truncated or numerically averaged disturbing potential to isolate individual resonances (e.g. Murray & Dermott, 1999; Morbidelli et al., 1995) have also found similar trends of resonance width with eccentricity, but with two significant differences: (i) the resonance width was found to diverge at the critical eccentricity, , and (ii) for , the non-monotonic variation of resonance width of the closer-in resonances was not observed. Our non-perturbative method finds finite, non divergent, widths near because it naturally accounts for the interaction of neighboring resonances; these interactions quench the singularity at which occurs in the numerically averaged disturbing potential. The non-monotonic variation of the resonance width is owed to the existence of the second resonance zone; this zone has not been discussed in previous studies.
The second resonance zone exists at eccentricities exceeding the critical planet-crossing value. For distant resonances, its width increases monotonically with increasing eccentricity, but for closer-in resonances, its width reaches a maximum then decreases and even vanishes and then reappears at higher eccentricities. The second resonance zone co-exists with but competes for phase space volume with the first resonance zone. We observe that the various transitions in the phase space structure, including the changes in the resonance widths with eccentricity, are related to the geometrical properties of the trace of the resonant orbit in the rotating frame. These results – particularly the characteristics of the second resonance zone and its influence on the first resonance zone – shed new light on Neptune’s exterior resonances, beyond what has been previously studied with either perturbative analytical theory in the low eccentricity regime (e.g. Murray & Dermott, 1999) or with numerical methods to estimate resonance strengths and widths (e.g. Morbidelli et al., 1995; Robutel & Laskar, 2001; Gallardo, 2018).
Below we summarize and discuss our results for the :1 and the :2 sequence of MMRs and then discuss implications for the dynamics of scattered disk objects.
5.1 The :1 sequence of MMRs
The 2:1 resonance is the closest-in of the :1 sequence of Neptune’s exterior resonances, and its phase space structure is also exemplary of these MMRs. With increasing eccentricity, we find that the structure in the Poincaré sections of in the neighborhood of the 2:1 exterior MMR undergoes two transitions, with the appearance and disappearance of different libration centers (Figure 2). The first resonance zone exists in almost the entire range of eccentricity from to 1. It contains three families of orbits. A pair of families of asymmetric librators with relatively small libration amplitudes has librating about two different centers displaced symmetrically from . These two libration centers are eccentricity-dependent and their locations vary from near at a lower eccentricity to near as the eccentricity approaches the planet-crossing minimum, . At even higher eccentricity, these two centers shift towards each other and merge again when the eccentricity exceeds , as shown in Figure 5. The boundary – or separatrix – of the asymmetric librators looks like the outline of a butterfly. The third family consists of symmetric librators surrounding the separatrix; it has a single libration center, at , for all values of particle eccentricity. The width in semi-major axis of the first resonance zone is measured by the outer boundary of the symmetric librators. The largest extent of this width occurs close to the planet-crossing limit, . The second resonance zone exists in the eccentricity range ; it has a single center of libration, at , and it co-exists with the first resonance zone (Figure 2). The width of the second resonance zone increases as eccentricity approaches unity.
All the :1 exterior resonances have similar phase space structure as the 2:1 resonance; their boundaries in the plane are summarized in Figure 6. The first resonance zone expands with increasing eccentricity, reaching a maximum size at the planet-crossing minimum, , where the perihelion distance is . For , the first resonance zone shrinks in width as the eccentricity approaches unity. The asymmetric libration zone first appears at a threshold value of the eccentricity which is near 0.04 for the 2:1 resonance; this threshold eccentricity is larger for larger (Figures 5 and 6). The second resonance zone first appears at eccentricity , where the perihelion distance is approximately . It expands in width as the eccentricity approaches unity (Figure 6).
With increasing value of the particle-to-Neptune period ratio, , the maximum width of the first resonance zone increases slightly from the 2:1 to the 3:1 MMR, then declines slowly for larger . In contrast, the maximum width of the second resonance zone rises steeply with ; its largest extent is smaller than that of the first resonance zone only at the period ratio of 2:1; for higher period ratios, it exceeds that of the first resonance zone (Figure 7a). At lower eccentricities, the resonance width decreases rapidly with (Figure 7b).
The dependence of the libration period on libration amplitude and eccentricity is rather complex (Figure 4). Notably, the libration period decreases with libration amplitude, in contrast with the behavior of the common pendulum model for resonance (see, e.g., Murray & Dermott (1999)) and also in contrast with the behavior of the “second fundamental model for resonance” (Henrard & Lemaitre, 1983); in both of these simplified models for resonance, the libration period monotonically increases with increasing libration amplitude, indicating that these models are oversimplified for the cases studied here. The small amplitude libration period in the first resonance zone increases monotonically with ; it has power-law dependence on for the planet-crossing minimum eccentricity , but has exponential dependence on for eccentricities away from (Figure 8). The small amplitude libration periods are on the order of times Neptune’s orbital period, shorter in the second resonance zone than in the first resonance zone.
For the special case of eccentricities close to , we can compare our numerical results with the analytical estimates derived by Pan & Sari (2004). The centers of libration seen in our Poincaré sections correspond to the “generalized Lagrangian points” discussed by these authors, who also used the terminology “generalized tadpole” and “generalized horseshoe” orbits for the asymmetric and symmetric librators of the first resonance zone, respectively. Pan & Sari (2004) estimated the resonance widths of tadpoles (asymmetric librators) and of horseshoes (symmetric librators) as and , respectively, and they estimated the minimum tadpole libration period as which is equivalent to . (We have multiplied by a factor of two the resonance half-widths stated in Pan & Sari (2004).) In our numerical studies with the fixed value of , the best-fit power-law relation that we found for the resonance width of the asymmetric librators (genearalized tadpoles) is ; this power-law index, 0.44, is lower than the estimate of 0.5 in Pan & Sari (2004), but overall their analytic estimate predicts widths similar to those we measured for the tadpoles (compare the curves labeled “Asymmetric librators” and “Tadpoles(Pan-2004) in Figure 7a). For the symmetric librators (generalized horseshoes), we found the widths to be nearly independent of , contrary to Pan & Sari (2004)’s estimate which predicts larger widths than we measured (compare the curves labeled “First resonance zone” and “Horseshoes(Pan-2004) in Figure 7a). Our numerical results on the small amplitude libration period (in the first resonance zones at particle eccentricity ) are generally similar to but slightly lower than Pan & Sari (2004)’s estimates (Figure 8(a)).
5.2 The :2 sequence of MMRs
The 3:2 resonance is the closest-in of the :2 sequence of Neptune’s exterior resonances, and its phase space structure is also exemplary of these MMRs. The Poincaré sections of in the neighborhood of the 3:2 resonance exhibit five transitions with increasing eccentricity. Its first resonance zone consists of a pair of stable libration islands centered at and and its second resonance zone consists of a pair of stable libration islands centered at and .(Figure 9). Starting from low eccentricity, the width of the first resonance zone increases, reaches a maximum at eccentricity , then shrinks with further increase of eccentricity; it vanishes at eccentricity but reappears again in the range of eccentricity exceeding . The second resonance zone exists in the range of eccentricity exceeding . The maximum width in semi-major axis of this resonance zone occurs at an eccentricity, .
All the :2 exterior resonances have similar structure to the 3:2 resonance (Figure 12), and similarly have five transitions in their phase space structure with increasing eccentricity. The first resonance zone expands with increasing eccentricity, reaching a maximum at , where the perihelion distance is approximately . The second resonance zone first appears at the second transition at , where the perihelion distance is approximately ; it expands in width with increasing eccentricity. The first resonance zone shrinks and vanishes at the third transition at . At a slightly higher eccentricity, , the widths of the stable islands of the second resonance zone reach a maximum at the fourth transition. At an even higher eccentricity, , the fifth transition occurs when the stable islands of the first resonance zone reappear and expand in width as the eccentricity approaches unity, while the islands of the second resonance zone shrink slightly. These last three transitions occur when the eccentricity reaches values that cause the two perihelion lobes in the rotating frame to intersect each other (for illustrations, see the lower middle panel in Figure 10 and Figure 2f in Malhotra et al. (2018)). The three eccentricity values, differ from each other only slightly; they correspond to perihelion distances of approximately and , respectively.
For the :2 resonances, the maximum width of the first resonance zone (which occurs at eccentricity , equivalent to perihelion distance ) decreases slowly with the increasing value of the particle-to-Neptune period ratio, but that of the second resonance zone increases steeply (Figure 13). The maximum widths of the first and second resonance zones achieve parity near the 9:2 resonance, and at higher particle-to-Neptune period ratios, the second resonance zone exceeds the maximum width of the first resonance zone.
We find that the dependence of the libration period on libration amplitude and eccentricity is rather complex, and departs from expectations based on the pendulum model. For fixed eccentricity , the libration period decreases with increasing libration amplitude, up to the boundary of the stable libration zone (Figure 11). Also in common with the behavior of :1 resonances, we find that the :2 resonances have small amplitude libration periods on the order of times Neptune’s orbital period, and these increase monotonically with ; the libration period has power-law dependence on for eccentricity values where the resonance zone has maximum width, but has exponential dependence on for other eccentricities (Figure 14).
5.3 Implications for SDOs
5.3.1 The second resonance zone seriously affects the dynamics of SDOs
At planet-crossing eccentricities, the existence of the second resonance zone affects the dynamics of the SDOs in two ways. Firstly, it squeezes the angular extent of the first resonance zone such that (as well as the resonant angle ) is limited to values well away from 0 (see Figures 2 and 9); this has the effect of pushing the resonance boundary far away from the location of the planet. From this result, we learn that the main reason that resonance sticking increases the dynamical lifetimes of SDOs is that the angular extent of the first resonance zone is small enough that resonant SDOs are effectively excluded from close approaches with Neptune. (Only the action of long term chaotic diffusion away from the resonance boundaries, possibly with significant influence of the collective effects of all the planets, allows eventual destabilizing close approaches to Neptune.) Secondly, it offers a new mechanism for planet-crossing orbits to avoid destabilizing close encounters with Neptune by “looping around the planet” during the portion of the orbit interior to the planet’s orbit; this is distinct from the previously understood “phase protection” mechanism in the first resonance zone which maintains the perihelion of the particle away from the longitude of the planet. This creates additional sticking opportunities for SDOs, in addition to those associated with the previously recognized first resonance zone.
5.3.2 An explanation for the outer perihelion distance boundary of SDOs
The largest extent of the first resonance zone of the :1 sequence and the :2 sequence of resonances occurs at eccentricity and , corresponding to perihelion distances 30.1 au and 33.1 au, respectively. We also find that the largest extent of the :3 resonances occurs at eccentricity , where the perihelion distance is about 34.6 au. The trend is that the largest extent of the first resonance zone of the :() resonances (integer , prime with respect to ) occurs at a perihelion distance slightly larger than that of : resonances (integer , prime with respect to ), and that the difference between two adjacent critical perihelion distances shrinks with the increasing value of . We observe that the critical perihelion distances of 30.1, 33.1 and 34.6 are in a geometric progression which converges to perihelion distance au. Because every :1 resonance is surrounded by an infinity of resonances, the implication is that the interaction of all these neighboring resonances produces a chaotic zone around each resonance which extends in perihelion distance from at least au up to au. An SDO which has perihelion distance near 30 au early in its scattering history, can thereby diffuse chaotically near the boundaries of Neptune’s resonances and lift its perihelion distance to au. This accounts for the most active region of resonance sticking of SDOs of perihelion distance –36) au, as first noted by Duncan & Levison (1997) in their simulations.
In a more realistic solar system model, this perihelion distance boundary is fuzzy and, on timescales of the age of the solar system, it is slightly larger than what is predicted by our simple restricted three body model; for example, Gladman et al. (2002) suggest that the boundary is at –38 au, a value commonly adopted in the literature.
|MMR||area (au)||cumulative area (au)||area (au)||cumulative area (au)|
|(up to au)||( au au)|
5.3.3 A large fraction of the SDOs’ phase space is stable resonant zones
Because “resonance sticking” in Neptune’s exterior MMRs appears to be the main mechanism for lengthening the dynamical lifetimes of scattered disk objects (Duncan & Levison, 1997; Lykawka & Mukai, 2007a; Yu et al., 2018), it is useful to measure the sizes of the resonance zone boundaries. As a rough indication of the sizes of resonance zone boundaries, we measured the area of the stable resonance zones in the plane for all the MMRs that we investigated. In measuring the resonant areas, we adopted a cutoff perihelion distance au because Neptune-resonant orbits of lower perihelion distance are subject to destabilizing perturbations from the other giant planets (Uranus, Saturn, etc.) on timescales of less than 10 megayears (e.g. Malhotra et al., 2018). Within this perihelion cutoff, most MMRs have only a small segment of the second resonance zone therefore we focussed on measuring the area of the first resonance zone. The measured resonant areas are listed in Table 1 and also plotted in Figure 15. We also measured and tabulated the resonant areas in the sub-region bounded by 26 au au; the upper boundary of this perihelion distance range is chosen to generously cover the fuzzy boundary of SDOs noted above.
The measured cumulative resonant area (Figure 15b) is found to be approximately 22% of the entire area in the plane bounded by 26 au au and 33 au au. This fraction should be considered a lower limit because we have not measured many resonances with particle-to-Neptune period ratios in-between the resonances that we have measured in our study. Furthermore, considering that the sizes of the first resonance zones decrease only very slowly with increasing particle-to-Neptune period ratio (Figures 7 and 13), we infer that more distant resonances (at au) also have stable zones of significant sizes. This implies that in the orbital parameter regime occupied by the SDOs, the long term stable resonant zones occupy a large fraction of the phase space, even at very large semi-major axes, explaining the prominence of the resonance sticking behavior of SDOs in numerical simulations.
Examining the tabulated results, we observe that, individually, the 2:1 and 3:2 resonances have very similar area in the plane, the :1 sequence of MMRs have the largest areas, and the resonant areas decrease rather slowly with increasing . The areas of the :2 resonances are comparable to but lower than those of the nearby :1 resonances; these also decrease rather slowly with increasing period ratio. These trends explain the general result from numerical simulations that the :1 resonances are the strongest/stickiest in the scattered disk, followed by the :2 resonances (Lykawka & Mukai, 2007a; Yu et al., 2018).
5.3.4 A tool to identify candidate resonant objects
The current method of identifying and classifying observed resonant objects (including resonance sticking SDOs) involves arduous visual examination of the libration behavior of the critical resonant angles in -megayear-long numerical integrations with a solar system model that includes the perturbations of the giant planets (e.g. Gladman et al., 2008; Volk et al., 2016). The boundaries of Neptune’s resonances in the plane that we have computed here can be used to identify candidate SDOs from knowledge of just their observed orbital elements, specifically, their barycentric semi-major axis and eccentricity, thus greatly reducing the computational and human effort. (The data for the resonance boundaries in the plane are available in electronic form upon request.) In the future, this tool can be made more useful by computing the resonance boundaries of a more complete set of Neptune’s resonances than we have attempted here.
: We thank Kathryn Volk for helpful discussions. LL acknowledges funding from National Natural Science Foundation of China (11572166) and China Scholarship Council. RM acknowledges funding from NASA (grant NNX14AG93G) and NSF (grant AST-1824869).
- Beauge (1994) Beauge, C. 1994, Celestial Mechanics and Dynamical Astronomy, 60, 225
- Chiang & Jordan (2002) Chiang, E. I., & Jordan, A. B. 2002, AJ, 124, 3430
- Di Sisto & Brunini (2007) Di Sisto, R. P., & Brunini, A. 2007, Icarus, 190, 224
- Duncan & Levison (1997) Duncan, M. J., & Levison, H. F. 1997, Science, 276, 1670
- Emel’yanenko et al. (2005) Emel’yanenko, V. V., Asher, D. J., & Bailey, M. E. 2005, MNRAS, 361, 1345
- Fehlberg (1968) Fehlberg, E. 1968, NASA TR, R-287, 1
- Fernández et al. (2004) Fernández, J. A., Gallardo, T., & Brunini, A. 2004, Icarus, 172, 372
- Gallardo (2018) Gallardo, T. 2018, Planet. Space Sci., 157, 96
- Gallardo (2019) —. 2019, Icarus, 317, 121
- Gladman et al. (2002) Gladman, B., Holman, M., Grav, T., et al. 2002, Icarus, 157, 269
- Gladman et al. (2008) Gladman, B., Marsden, B. G., & Vanlaerhoven, C. 2008, Nomenclature in the Outer Solar System (in The Solar System Beyond Neptune, Barucci, M. A. and Boehnhardt, H. and Cruikshank, D. P. and Morbidelli, A. and Dotson, R. (eds.), University of Arizona Press, Tucson), 43–57
- Gomes et al. (2008) Gomes, R. S., Fernández, J. A., Gallardo, T., & Brunini, A. 2008, The Scattered Disk: Origins, Dynamics, and End States, ed. M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, A. Morbidelli, & R. Dotson, 259–273
- Henrard & Lemaitre (1983) Henrard, J., & Lemaitre, A. 1983, Celestial Mechanics, 30, 197
- Kotoulas & Voyatzis (2004a) Kotoulas, T., & Voyatzis, G. 2004a, Celestial Mechanics and Dynamical Astronomy, 88, 343
- Kotoulas & Voyatzis (2004b) —. 2004b, Celestial Mechanics and Dynamical Astronomy, 88, 343
- Kotoulas (2005) Kotoulas, T. A. 2005, A&A, 429, 1107
- Levison & Duncan (1997) Levison, H. F., & Duncan, M. J. 1997, Icarus, 127, 13
- Levison et al. (2006) Levison, H. F., Duncan, M. J., Dones, L., & Gladman, B. J. 2006, Icarus, 184, 619
- Luu & Jewitt (2002) Luu, J. X., & Jewitt, D. C. 2002, ARA&A, 40, 63
- Lykawka & Mukai (2007a) Lykawka, P. S., & Mukai, T. 2007a, Icarus, 189, 213
- Lykawka & Mukai (2007b) —. 2007b, Icarus, 192, 238
- Malhotra (1995) Malhotra, R. 1995, AJ, 110, 420
- Malhotra (1996) —. 1996, AJ, 111, 504
- Malhotra et al. (2018) Malhotra, R., Lan, L., Volk, K., & Wang, X. 2018, ArXiv e-prints, arXiv:1804.01209
- Malhotra & Williams (1997) Malhotra, R., & Williams, J. G. 1997, Pluto’s Heliocentric Orbit (in Pluto and Charon, Stern, S. A. and Tholen, D. J. (eds.), University of Arizona Press, Tucson), 127
- Milani et al. (1989) Milani, A., Nobili, A. M., & Carpino, M. 1989, Icarus, 82, 200
- Morbidelli (1997) Morbidelli, A. 1997, Icarus, 127, 1
- Morbidelli & Brown (2004) Morbidelli, A., & Brown, M. E. 2004, The kuiper belt and the primordial evolution of the solar system (Comets II, M. C. Festou, H. U. Keller, and H. A. Weaver (eds.), University of Arizona Press, Tucson), 175–191
- Morbidelli et al. (1995) Morbidelli, A., Thomas, F., & Moons, M. 1995, Icarus, 118, 322
- Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics, 1st edn. (New York, New York: Cambridge University Press)
- Murray-Clay & Chiang (2005) Murray-Clay, R. A., & Chiang, E. I. 2005, ApJ, 619, 623
- Nesvorný & Roig (2000) Nesvorný, D., & Roig, F. 2000, Icarus, 148, 282
- Nesvorný & Roig (2001) —. 2001, Icarus, 150, 104
- Pan & Sari (2004) Pan, M., & Sari, R. 2004, AJ, 128, 1418
- Robutel & Laskar (2001) Robutel, P., & Laskar, J. 2001, Icarus, 152, 4
- Tiscareno & Malhotra (2003) Tiscareno, M. S., & Malhotra, R. 2003, AJ, 126, 3122
- Tiscareno & Malhotra (2009) —. 2009, AJ, 138, 827
- Volk & Malhotra (2008) Volk, K., & Malhotra, R. 2008, ApJ, 687, 714
- Volk et al. (2016) Volk, K., Murray-Clay, R., Gladman, B., et al. 2016, AJ, 152, 23
- Volk et al. (2018) Volk, K., Murray-Clay, R. A., Gladman, B. J., et al. 2018, AJ, 155, 260
- Voyatzis & Kotoulas (2005) Voyatzis, G., & Kotoulas, T. 2005, Planet. Space Sci., 53, 1189
- Wang & Malhotra (2017) Wang, X., & Malhotra, R. 2017, AJ, 154, 20
- Yu et al. (2018) Yu, T. Y. M., Murray-Clay, R., & Volk, K. 2018, AJ, 156, 33