Chaotic Dynamics of Trans-Neptunian Objects Perturbed by Planet Nine

As described in our introduction, one of the key outcomes of the orbital integrations undertaken by BB16a and BB16b, was that introducing a distant, eccentric Planet Nine into their simulations caused a clustering, at late times, of test-particle orbits into a configuration that was anti-aligned with the Planet Nine orbit. We provide an example of this from our own simulations.

We have conducted coplanar simulations in which Neptune is placed on circular orbit at a = 30 au and is accompanied by Planet Nine on a a = 500 au, e = 0.6 orbit. In our fiducial simulation, we take Planet Nine's mass to be $10\,{M}_{\oplus }.$ In Section 2.3, we compare our fiducial case to two additional simulations: one in which Neptune is replaced by a quadrupole approximation of its orbit-averaged contribution to the gravitational potential, and one in which Planet Nine's mass is reduced to $1\,{M}_{\oplus }$. We initialized 3000 test-particles in orbits with $a\in [150,550]\,\mathrm{au}$ and $q\in [33,50]\,\mathrm{au}$ and random mean longitudes. Particles' longitudes of perihelia, ϖ, are chosen to be either exactly aligned (${\rm{\Delta }}\varpi \equiv {\varpi }_{P9}-\varpi =0$) or anti-aligned (${\rm{\Delta }}\varpi =\pi$) with Planet Nine's periapse direction, ${\varpi }_{9}$. In other words, the test-particles are placed in a population similar to those used by BB16a and BB16b. We integrate our systems using the IAS15 integrator (Rein & Spiegel 2015) in the REBOUND code of Rein & Liu (2012).

Figure 1 shows the results of our fiducial simulation, where we plot the final 2.5 Gyr of data for particles that survive the full 5 Gyr of integration. Only ∼8% of the initial 3000 particles survive the full integration and are plotted. Two distinct dynamical populations are evident. The first population, consisting of particles interior to $a\sim 250\,\mathrm{au}$, evolves at approximately fixed semimajor axes with either apsidally aligned or circulating orbital orientations relative to Planet Nine. This interior aligned population exhibits clear structure in the a–${\rm{\Delta }}\varpi$ space that correlates with variations in perihelion distance. In Section 3, we demonstrate that this structure is well-explained by secular theory. The aligned population consists entirely of particles that are initialized with ${\rm{\Delta }}\varpi =0$. Some additional particles initialized with ${\rm{\Delta }}\varpi =0$ survive the full simulation after scattering to semimajor axes $a\gt 500\,\mathrm{au}$ and maintaining high pericenter distances. These particles appear as the vertical yellow lines seen at large a values in Figure 1. If TNOs existed on such orbits, their large heliocentric distances would render them essentially unobservable. We devote little time to these particles in the rest of this paper because they are unimportant for explaining current observations.

Zoom In Zoom Out Reset image size

The second population, by contrast, consists of particles on distant orbits ($a\gtrsim 250\,\mathrm{au}$) that maintain anti-alignment with the orbit of Planet Nine (${\rm{\Delta }}\varpi \sim \pi$) and low pericenter distances $q\lt 100\,\mathrm{au}$ while wandering in semimajor axis. In our simulations, this population consists entirely of particles that are initialized with ${\rm{\Delta }}\varpi =\pi$. Similar populations were observed in simulations conducted by BB16a and BB16b. This second population of anti-aligned particles is the fundamental means by which they explain the apparent orbital clustering of observed ESDOs. As noted by BB16a, the dynamical evolution of the anti-aligned population, specifically their stochastic wandering in semimajor axis, is indicative of chaotic diffusion through a web of overlapped resonances. In Section 4, we examine the structure of this chaotic web in detail, by illustrating how the mixture of resonant and chaotic trajectories are situated in the phase-space of the particles. We also address the mechanism by which the chaotically scattering particles maintain their persistent anti-alignment.

Malhotra et al. (2016) have noted that the periods of the four ESDOs, Sedna, 2010 GB174, 2004 VN112, and 2012 VP113, could place them in a series of $N/1$ and $N/2$ resonances with a putative Planet Nine at a semimajor axis ${a}_{9}\approx 665\,\mathrm{au}$. Millholland & Laughlin (2017) find that a similar Planet Nine semimajor axis leads to commensurabilities with a larger set of ESDOs. However, we note that, in our simulations, essentially all of the test particles initialized in the anti-aligned configuration experience significant semimajor diffusion for the full duration of the simulation (evident from the horizontal striations in the anti-aligned regions of Figure 1), and therefore are not found to occupy any long-term stable librating resonant configurations. Based on our numerical integrations, a significant population of long-term stable, resonant ESDOs seem unlikely to arise without dissipative forces that would aid permanent resonant capture. However, we often observe "resonance sticking" where particles temporarily occupy orbits in or near resonances in both the aligned and anti-aligned populations. Figure 2 illustrates resonance sticking for a test particle that, in the course of its chaotic evolution, is temporarily captured in various resonances. The figure plots the semimajor axis evolution of the particle, along with that of Planet Nine's mean anomaly, M₉, calculated when the test-particle comes to pericenter (i.e., the test-particle mean anomaly is M = 0). When there is an N:k resonance between the test-particle and Planet Nine, this angle librates about k discrete values when the particle is at pericenter. Otherwise, outside of resonance, the angle will uniformly fill the interval [0, 2π). The test-particle in Figure 2 is temporarily captured in the 5:9 interior MMR from 0.5 to 0.7 Gyr, the 71:5 exterior MMR from approximately 1–2 Gyr, and the exterior 5:1 MMR around 4.5 Gyr.

Zoom In Zoom Out Reset image size

Temporary resonance occupation has been seen, by various other authors, in simulations of TNOs in the presence of Planet Nine. As noted in the introduction, Batygin & Morbidelli (2017) observe resonance-hopping behavior when particles are given small orbital inclinations relative to Planet Nine. However, the test-particles in the simulations of Batygin & Morbidelli (2017), which do not fully account for Neptune's gravitational influence, appear to spend a larger fraction of their orbital evolution in resonances when compared with our fiducial simulation. Becker et al. (2017) observe resonance-hopping behavior in simulations of observed TNOs influenced by Neptune, Planet Nine (which, in their simulations, is inclined 30° to the ecliptic), and the quadrupole potential of the other outer solar system planets. In simulations from Becker et al. (2017), TNOs experience significant migration (changes in semimajor axis of ${\rm{\Delta }}a\gt 100\,\mathrm{au}$) roughly as often as they stay close their initial semimajor axis while residing in a single resonance or exhibiting small jumps between neighboring resonances. Millholland & Laughlin (2017) find persistent diffusion of semimajor axes for most/all particles, and hence only temporary residence near any given MMR in their numerical simulations that include Neptune's full gravitational potential and the quadrupole potential of the other giant planets (S. Millholland 2018, private communication).

In addition to our fiducial simulation, we have run one in which Neptune is replaced by its orbit-averaged quadrupole contribution to the gravitational potential and one with a $1\,{M}_{\oplus }$ Planet Nine. Hereafter, we will refer to these as the "quadrupole" and "low-mass" simulations, respectively. The results of these simulations are shown in Figures 3 and 4. Both show significant differences from the fiducial simulation presented in Section 2.1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Approximately ∼1/3 of the particles survive up to 5 Gyr in the quadrupole simulation (Figure 3), significantly more than the ∼8% survival fraction of the fiducial simulation. Significantly more anti-aligned particles remain at semimajor axes interior to Planet Nine in the quadrupole simulation than in the fiducial simulation. These particles maintain constant semimajor axes for most or all of the simulation, causing the vertical striations seen in Figure 3 near ${\rm{\Delta }}\varpi \approx \pi$ and $a\lesssim 500\,\mathrm{au}$. In contrast, anti-aligned population members undergo significant chaotic semimajor axis evolution in the fiducial simulations.

Batygin & Morbidelli (2017) investigate the dynamics of test-particles under the gravitational influence of Planet Nine and a quadrupole potential representing the solar system's giant planets. Thus, their simulations are qualitatively similar to our quadrupole simulation, and much of their analysis readily explains the results of our quadrupole simulation. In particular, they show that:

• 1.  
  Anti-aligned particles survive the duration of the simulation in MMRs with Planet Nine that protect them from close encounters.
• 2.  
  The secular evolution of these resonant particles maintains their anti-alignment with Planet Nine. Beust (2016) also demonstrated that resonant particles are expected to maintain anti-alignment over the course of their secular evolution in the presence of Planet Nine and the quadrupole potential of the outer solar system.

The quadrupole simulations adequately capture the features of the circulating/aligned population of test-particles seen in the fiducial simulation. This is because, as we show below in Section 3, their dynamics are well-described by a secular approximation. However, the dynamics of the anti-aligned population are significantly more chaotic in the fiducial simulation. Batygin & Morbidelli's (2017) explanation of the dynamical mechanism maintaining anti-alignment—secular evolution in resonance—appears to be incomplete: the inclusion of Neptune's full gravitational potential disrupts nearly all stable resonant librations in the fiducial simulation. Nonetheless, the fiducial simulation still exhibits a population of particles that maintain persistent anti-alignment with Planet Nine. In Section 4.5, we show that Neptune alone can induce significant chaotic behavior in particles with semimajor axes in our initial range $a\in [150,550]\,\mathrm{au}$, even for pericenter distances as large as $q\sim 50\,\mathrm{au}$. Therefore, it is unsurprising that the fiducial simulation shows significantly different evolution from the quadrupole simulation.

Figure 4 shows the results of our low-mass simulation. As in the fiducial simulation of Figure 1, the test-particles can roughly be separated into two populations: a "circulating/aligned" population at small semimajor axes and an "anti-aligned" population that experiences chaotic semimajor axis evolution. The most striking difference is that the "anti-aligned" population is no longer directly anti-aligned with Planet Nine's orbit, but rather appears to be centered around ${\rm{\Delta }}\varpi \sim 3\pi /2$. In Section 4.6, we argue that this concentration arises because the distribution of perihelion distances and longitudes does not relax to a steady state over the course of our 5 Gyr simulation.

We conducted some preliminary numerical simulations with Planet Nine's mass increased to $100\,{M}_{\oplus }$. These simulations produce orbital clustering similar to what is observed in the fiducial simulation, albeit with substantially fewer surviving test particles. We do not consider this high-mass simulation further, because observational constraints probably rule out such a large Planet Nine mass (Luhman 2014).

We have also run a simulation identical to the fiducial simulation presented in Section 2.1, but with each test particle given a random initial ϖ, rather than starting from perfectly aligned and anti-aligned configurations. We find that only 58 particles out of an initial population of 1500 (i.e., ∼4%) survive the full simulation when started from random initial apsidal alignments. Randomizing initial apsidal alignments reduces the particles' survival rate, because it places many more particles on secular trajectories that evolve through the crossing/non-crossing boundary (see Section 3 below). As we will see in Section 4, particles experience strong resonance overlap as they approach this boundary, and they tend to be ejected. We plan to explore the consequences of randomizing initial particle alignment more thoroughly in a forthcoming paper (G. Li et al. 2018, in preparation).

To summarize, the key features of our long-term simulations (which we seek to explain in Sections 3 and 4) are:

• 1.  
  Surviving particles inside $a\lesssim 250\,\mathrm{au}$ have apsidal lines that circulate or maintain alignment with Planet Nine's. These particles in this aligned population do not experience large excursions in semimajor axis.
• 2.  
  The majority of surviving particles outside $\sim 250\,\mathrm{au}$ in the fiducial simulation are anti-aligned with Planet Nine's orbit. These particles maintain their anti-alignment over the course of the simulation while diffusing in semimajor axis.
• 3.  
  Nearly all surviving anti-aligned particles were initialized with an anti-aligned configuration. Likewise, the aligned population is composed of particles initially in an aligned configuration.
• 4.  
  Particles often temporarily "stick" to MMRs, but no particles survive the full simulation librating stably in a single resonance. This is in contrast to the "quadrupole" simulation without an active Neptune, where many particles are observed to stably librate in various MMRs with Planet Nine for the duration of the simulation.
• 5.  
  Substantially more test particles are retained when Planet Nine's mass is decreased to $1{M}_{\oplus }$. The population of test particles can again be divided into two separate populations: one with circulating or aligned longitudes of perihelia that remain at roughly constant semimajor axis, and a second population that undergoes chaotic semimajor axis diffusion. The orbital alignments of the chaotically scattering population cluster around ${\rm{\Delta }}\varpi \sim 3\pi /2$, rather than ${\rm{\Delta }}\varpi \sim \pi$ as observed when Planet Nine has a mass of $10\,{M}_{\oplus }$.

Inside ${a}_{\mathrm{crit}}$, the secular evolution of most aligned particles from their initial configuration with perihelia in the range $q\in [35\,\mathrm{au},50\,\mathrm{au}]$ follow trajectories that librate about ${\rm{\Delta }}\varpi =0$ and avoid evolving onto crossing orbits. The disappearance of the aligned population beyond ${a}_{\mathrm{crit}}$ occurs where the secular trajectories on which aligned particles are initialized all evolve onto crossing orbits.³ This can be seen in Figure 5, where trajectories starting in the initial range terminate at the red line separating crossing from non-crossing orbits, rather than closing on themselves above this line. Orbits on the red crossing/non-crossing line meet Planet Nine's orbit tangentially. Orbits near this line are prone to close encounters with Planet Nine. The strength of these close encounters is greatly enhanced by the fact that they occur near the test particle's apoapse, when the particle is moving most slowly. These encounters tend to eject the test particles. In Section 4, we more thoroughly examine the onset of chaos induced by resonance overlap as the orbit-crossing line is approached.

Essentially all particles in the anti-aligned population experience chaotic semimajor axis evolution; hence, a purely secular treatment is not strictly applicable to their dynamics. Nonetheless, the morphology of anti-aligned libration islands appears to match the distribution of anti-aligned orbits. The anti-aligned particles are concentrated near the secular trajectories that correspond to their initial periapse distances. Few anti-aligned particles exist in the semimajor axis range ($145\,\mathrm{au}\lt a\lt 155\,\mathrm{au}$) shown in first panel of Figure 5. In this range, the initial secular trajectories evolve toward the boundary from crossing to non-crossing. Resonance overlap is strongest near this transition, as described previously and demonstrated in Section 4. We will return to the question of the anti-aligned population's apsidal confinement mechanism in Section 4.6.

Surviving particles in the aligned population maintain relatively constant semimajor axes. Consequently, secular averaging accurately approximates their long-term dynamics, as we saw in Section 2. However, a large fraction (88%) of initially aligned particles are ejected from the simulation, and many of the surviving particles experience limited chaotic variations in semimajor axis. Therefore, chaos caused by resonance overlap plays some role in shaping the distribution of the aligned population.

In Section 3, we saw that initially aligned particles do not survive on secular trajectories that lead to orbit-crossing with Planet Nine. We argued that this is because these particles suffer strong close encounters with Planet Nine that lead to their ejection. In fact, initially aligned particles that evolve secularly onto nearly crossing orbits experience significant chaotic variations in their semimajor axes and are often ejected. This is in agreement with Figure 6, where it can be seen that significant chaos occurs slightly before eccentricities reach orbit-crossing values.

In Figure 7, we examine an example of how secular evolution onto nearly crossing orbits brings particles into a chaotic region of phase space. Figure 7 shows a series of stability maps in test-particle semimajor axis a and mean longitude λ. The simulations are initialized with Planet Nine at its pericenter. Also shown in Figure 7, in the top-left panel, is a contour/density plot in the $(e\sin \varpi ,e\cos \varpi )$ plane reproduced from Figure 5, with a series of three points labeled "A" through "C" along a secular trajectory chosen from the range of our fiducial simulation's initial conditions. Each stability map is computed by initializing the test particles' eccentricities and longitudes of perihelia to the values at one of these points along the secular trajectory. From Figure 7, we see how particles' orbits slowly evolve into progressively more chaotic regions of phase space. Initially, orbits start near the top-most point on the secular trajectory shown in Figure 7. Here, orbits behave in an entirely regular fashion on short timescales. As secular evolution evolves orbits clockwise around the secular trajectory, they reach point "A." The stability map corresponding to point "A" indicates mostly regular trajectories (black), with faint hints of chaos beginning to appear (white). Subsequent panels "B" and "C" show a progressively more chaotic phase space punctuated by regular islands associated with libration in various MMRs. For example, there is a prominent stable black region associated with libration in the 3:1 MMR at $a\approx 240\,\mathrm{au}$. Upon reaching eccentricities and longitudes of perihelia near point "C," many of the simulation particles experience significant resonance overlap resulting in vigorous chaos and ejection. For particles that manage to survive, continued secular evolution brings their orbits back to regions of phase space where resonances are no longer overlapped. (The secular trajectory is symmetric about $e\sin {\rm{\Delta }}\varpi =0$, and stability maps computed with $e\sin {\rm{\Delta }}\varpi \lt 0$ mirror those with $e\sin {\rm{\Delta }}\varpi \gt 0$ but $\lambda \to -\lambda$, so that after reaching "C," particles effectively go through the same sequence in reverse order.)

Zoom In Zoom Out Reset image size

The critical semimajor axis, ${a}_{\mathrm{crit}}$, beyond which initially Planet Nine-aligned orbits do not survive, is the semimajor axis at which secular trajectories bring all particles sufficiently deep into the chaotic zone near orbit-crossing that they are all ejected. The mass of Planet Nine influences where this boundary occurs via two effects: first, a larger Planet Nine mass changes the contours of the secular Hamiltonian, and therefore, the initial pericenter distances of orbits that evolve to become nearly Planet Nine-crossing. Second, a larger Planet Nine mass causes the onset of large-scale chaos to occur further from the orbit-crossing boundary.

In contrast with the aligned population, nearly all surviving members of the anti-aligned population are initialized on Planet Nine-crossing orbits. From Figure 6, we see that this places the anti-aligned particles in a region of phase-space suffused by an intricate mixture of regular and chaotic trajectories generated by the overlap of Planet Nine MMRs. Figure 8 illustrates some slices through this chaotic web in greater detail, showing two slices of phase space for particles with constant perihelion distances and semimajor axes both interior and exterior to Planet 9. In both panels, the initial conditions are such that they are anti-aligned with Planet Nine's orbit (i.e., $\varpi ={\varpi }_{9}+\pi$). In the top panel, corresponding to semimajor axes interior to Planet Nine's orbit, Planet Nine is initially at pericenter and the test particle's initial mean anomaly is plotted on the vertical axes, which we have labeled ϕ to match notation introduced in Section 4.3. In the bottom panel, showing orbits exterior to Planet Nine, the roles are reversed: the test particle is initially at pericenter and Planet Nine's initial mean anomaly is plotted on the vertical axes and labeled ϕ.

Zoom In Zoom Out Reset image size

We can understand much of the structure observed in Figure 8 by considering initial conditions that lead to a close encounter with Planet Nine. Let ${M}_{9,+}^{* }$ and ${M}_{9,-}^{* }$ be the mean anomalies of Planet Nine at the two points where the test particle's orbit intersects Planet Nine's. Also, let ${M}_{+}^{* }$ and ${M}_{-}^{* }$ be the mean anomalies of the test particle at these same intersection points. The test particle will experience a close encounter if it reaches ${M}_{\pm }^{* }$ at the same time Planet Nine reaches ${M}_{9,\pm }^{* }$ after completing an integer number of orbits. In other words, there is a time t that satisfies

$$\begin{eqnarray}&&{nt}={M}_{\pm }^{* }+2\pi j\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{n}_{9}t={M}_{9,\pm }^{* }+2\pi {j}_{9}\end{eqnarray} \tag{ 2 }$$

for some pair of integers j and j₉, where n and n₉ are the mean motions of the particle and Planet Nine. This condition is equivalent to

$$\begin{eqnarray}&&{M}_{9,\pm }^{* \pm }-\displaystyle \frac{{n}_{9}}{n}{M}_{\pm }^{* }=2\pi \displaystyle \frac{{n}_{9}}{n}j\mathrm{mod}2\pi \,.\end{eqnarray} \tag{ 3 }$$

If the orbital frequency ratio $n/{n}_{9}$ is irrational, then the test particle will come arbitrarily close to Planet Nine after a sufficient amount of time has passed. Close encounters can only be avoided if frequency ratio is of the form $n/{n}_{9}=k/N$ for integer N and k, i.e., the particle is in resonance with Planet Nine. In Figure 8, we have traced the initial test particle orbital phases, as a function of semimajor axis, that lead to collisions with Planet Nine in less than two orbits going forward or backward in time. These lines trace clear features in the chaotic web of white pixels shown in both panels. By tracing initial conditions that lead to collisions after more and more orbits, we would fill in the "backbone" of the chaotic web. Continuing this tracing process indefinitely would densely fill the $a\mbox{--}\phi$ plane in Figure 8, except at exactly resonant period ratios. Clearly, the entire phase-space is not densely filled by chaotic trajectories. This is because, sufficiently close to a resonant period ratio, gravitational interactions cause the test particle's instantaneous orbital period to oscillate about the resonant value. Accordingly, this provides resonances with finite widths in semimajor axis space. A multitude of resonances are evident in Figure 8 as stable "islands" of black pixels in the white chaotic "sea." The MMRs in Figure 8 exhibit obvious structure in the $a\mbox{--}\phi$ plane, which we describe in the next section.

Before discussing chaos driven by overlap of MMRs, we first describe the dynamics of a single isolated MMR. For concreteness, we will restrict our discussion to test particles in exterior resonances with Planet Nine, though interior resonances can easily be considered by analogous means. We will refer to MMRs occurring at period ratios $P/{P}_{9}=N/k$, where N and k are relatively prime integers, as kth order resonances. We caution the reader that this is not conventional terminology when referring to MMRs but, as Pan & Sari (2004) note, it is an appropriate redefinition when considering highly eccentric orbits as we are here. While the dynamics of MMRs are often examined by considering the libration or circulation of the resonant angles

$$\begin{eqnarray}&&\phi =\displaystyle \frac{1}{k}(N\lambda -k{\lambda }_{9}-(N-k)\varpi )\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\psi =\displaystyle \frac{1}{k}(N\lambda -k{\lambda }_{9}-(N-k){\varpi }_{9})\end{eqnarray} \tag{ 5 }$$

and their linear combinations, it is helpful to relate these angles to a more physical picture. To this end, we will consider stroboscopic snapshots of the system at each pericenter passage of the test particle. When the test particle is at pericenter, $\lambda =\varpi$, and hence the resonant angle reduces to $\phi \,=(\varpi -{\lambda }_{9})$, the angular separation between the test particle's longitude of pericenter and Planet Nine's mean longitude.

Figure 9 shows representative phase-space portraits for N:1 and N:2 exterior resonances with Planet Nine in the $\phi \mbox{--}{\rm{\Delta }}a/a$ plane. The phase-space portraits show constant energy contours of resonant Hamiltonians constructed via an averaging procedure for particle orbits anti-aligned with Planet Nine's (see Appendix A for details). The portraits are accompanied by schematic diagrams showing the location of Planet Nine in its orbit when the test particle is at pericenter for the various critical points of the phase space. On timescales sufficiently shorter than the apsidal precession timescale, resonant test-particle trajectories are confined to the constant energy contours.⁵ The dashed red curves show the separatrix dividing librating and circulating trajectories.

Zoom In Zoom Out Reset image size

First, we will consider a test particle in an exterior N:1 resonance with Planet Nine. Ignoring any orbital perturbations, after one complete test-particle orbit, Planet Nine will have returned to its initial orbital phase if the particle's orbital period is exactly commensurate with Planet Nine's. In actuality, Planet Nine's gravitational effects will cause changes in the test particle's semimajor axis and apsidal precession. Ignoring apsidal precession for the moment, the induced change in semimajor axis dictates the dynamics of the resonance on short timescales by altering the particle's orbital period. An increase (decrease) in the test particle's orbital period causes a corresponding decrease (increase) in ϕ at the test particle's next pericenter passage, as Planet Nine will have advanced through slightly more (less) than N complete orbits. A fixed point of the resonance occurs when there is no net change to the test particle's semimajor axis over the course of its orbit, so ϕ is unchanged when the test particle returns to pericenter. For perfectly anti-aligned test particles orbits, $\phi =0$ and π are always fixed points by symmetry. The point $\phi =\pi$, i.e., when Planet Nine is at apoapse during the particle's pericenter passage, is a stable fixed point: initial conditions that are slightly displaced from $(a,\phi )=({N}^{2/3}{a}_{9},\pi )$, where a and a₉ are the semimajor axes of the particle and Planet Nine, respectively, will oscillate about this point. If these oscillations are too large, however, the test particle will collide with Planet Nine at the point where the orbits cross. The loci of points corresponding to collisions with Planet Nine are shown in Figure 9 with red lines. These collision points set the maximum amplitude of stable oscillations in ϕ about the fixed point at π. There are two additional stable "asymmetric libration" fixed points for N:1 resonances.⁶ The ϕ location of the asymmetric libration points depends on the particle's perihelion distance.

The fixed points of an N:k resonance with k > 1, correspond to k-periodic orbits when viewed in terms of stroboscopic snapshots at pericenter: the test particle sees Planet Nine at a repeating series of k distinct points at sequential pericenter passages. For example, in an N:2 MMR like the one shown at the bottom of Figure 9, there is a stable fixed point associated with the test particle's pericenter passages occuring alternatively when Planet Nine is at ${M}_{9}=\pi /2$ then ${M}_{9}=3\pi /2$. Additionally, there is an unstable fixed point corresponding to Planet Nine's orbital phase occuring successively at $\phi =0,\pi$ as well as two pairs of ϕ values associated with asymmetric librations.

For interior resonances where the test particle's orbital period is less than Planet Nine's, essentially every aspect of the dynamics described above are the same when one considers stroboscopic snapshots capturing Planet Nine's pericenter passage rather than the test particle's.

On a longer timescale, apsidal precession slowly modulates the phase-space portrait, changing the location of the resonance's fixed points and the shape of the resonant islands. Anti-aligned test particles in our fiducial numerical simulations do not precess far from perfect anti-alignment, so these effects are modest.

Having examined the dynamics of individual resonances, we now turn to a more global picture of the dynamics of the anti-aligned particles. In order to do so, we will approximate the dynamics of these particles by means of an area-preserving map. This mapping captures the effects of the infinite number of MMRs with Planet Nine in any given semimajor axis interval and the chaos driven by their overlap. The full derivation of our map, which we summarize here, is presented in Appendix B. The mapping approximates stroboscopic snapshots of a test particle's orbital energy E and Planet Nine's mean anomaly ϕ at each perihelion passage of the particle. Similar maps have been derived and used previously in the study of highly eccentric particles subject to a planetary perturber by a number of authors (e.g., Petrosky & Broucke 1988; Malyshkin & Tremaine 1999; Pan & Sari 2004; Shevchenko 2011). Between each perihelion passage, the gravitational influence of Planet Nine imparts a change in energy to the test particle. Because the gravitational influence of Planet Nine is minimal when the particle is near aphelion and far from Planet Nine, we approximate the particle's orbit as a parabolic orbit with the same perihelion distance as the true particle orbit. The change to the particle's orbital energy changes its orbital period, which in turn affects Planet Nine's orbital phase at the time of the particle's next perihelion passage.

The mapping is given by

$$\begin{eqnarray}E^{\prime} & = & E+{\mu }_{9}\delta E(\phi ;q,{\rm{\Delta }}\varpi )\\ \phi ^{\prime} & = & \phi +2\pi {\left(\displaystyle \frac{{E}_{0}}{E^{\prime} }\right)}^{3/2}\mathrm{mod}2\pi \end{eqnarray} \tag{ 6 }$$

where primed and un-primed variables represent values before and after a mapping step, respectively, $\delta E$ is a function that gives the increment in E as a function of Planet Nine's orbital phase, as well as the test-particle perihelion distance q and apsidal alignment ${\rm{\Delta }}\varpi =\varpi -{\varpi }_{9}$, and E₀ is the orbital energy of a test particle with the same semimajor axis as Planet Nine. We describe our method for determining $\delta E$ via numerical integration in Appendix A. Because particles' eccentricities and orbital alignments evolve slowly, compared to their semimajor axes, we treat q and ${\rm{\Delta }}\varpi$ as fixed parameters when applying the mapping.

To get a better understanding of the local phase space structure in a given semimajor axis range, we expand the mapping, Equation (6), about a semimajor axis, ${a}_{\mathrm{res}}$, corresponding to an N:1 MMR with Planet Nine. We define a new variable $x=2N(a\mbox{--}{a}_{\mathrm{res}})/3{a}_{\mathrm{res}}$ and obtain the new mapping

$$\begin{eqnarray}x^{\prime} & = & x+\epsilon \delta E(\phi ;q,{\rm{\Delta }}\varpi )\\ \phi ^{\prime} & = & \phi -2\pi x^{\prime} \mathrm{mod}2\pi \end{eqnarray} \tag{ 7 }$$

where $\epsilon =-(3/2){\mu }_{9}{\left(\tfrac{{a}_{\mathrm{res}}}{{a}_{9}}\right)}^{5/2}$. Particles become more loosely bound as the semimajor axis, ${a}_{\mathrm{res}}$, increases, such that a smaller ${\mu }_{9}$ can have the same dynamical effect at large ${a}_{\mathrm{res}}$ as a greater ${\mu }_{9}$ does at a smaller ${a}_{\mathrm{res}}$. This is reflected in 's scaling with ${\mu }_{9}$ and ${a}_{\mathrm{res}}$. In the new mapping, Equation (7), first-order resonances occur at integer values of x where ϕ advances by an integer multiple of 2π. The mapping in Equation (7) is unchanged if we add an integer to x so we can obtain a complete picture of the phase-space structure by taking x modulo 1. Therefore, phase-space regions between successive N:1 MMRs are identical to the extent that the linearization approximation we used to derive Equation (7) remains valid. Repeating, nearly identical phase-space structure between successive N:1 MMRs is evident in the stability map shown in the bottom panel of Figure 8. Equation (7) provides a good approximation of a test particle's dynamics so long as $| x| \ll N$ because we assumed that $| a\,-\,{a}_{\mathrm{res}}| \ll {a}_{\mathrm{res}}$ to derive it from Equation (6).

Some sample mapping trajectories of Equation (7) are plotted in Figure 10. These trajectories reveal a number of notable features of the dynamics of the anti-aligned population:

• 1.  
  Chaotic trajectories fill much of the space for all parameter values. There are no KAM curves spanning the full phase space from $\phi =0$ to $\phi =2\pi$. Therefore, there should be no barriers in phase space bounding the chaotic diffusion of test particles' semimajor axes. The lack of bounding KAM curves reflects the fact that, for any initial period ratio, there are initial values of ϕ that result immediately in a close encounter with Planet Nine.
• 2.  
  Increasing increases the size of low-order resonant islands while subsuming many of the higher-order resonant islands in the chaotic sea.
• 3.  
  The middle panels of Figure 10 demonstrate the effects of a modest deviation from exact anti-alignment: the resonance centers become shifted in ϕ and the resonant islands become slightly distorted in shape. Otherwise, much of the resonant structure remains unchanged.
• 4.  
  Comparing the left- and right-hand columns of Figure 10 illustrates the influence of the perihelion distance, q, on the test-particle dynamics. Here, q determines the extent of resonant islands in the ϕ coordinate by setting the ϕ value at which collisions occur. This same effect was seen to determine the extent of individual resonances in Section 4.3.

The lack of bounding KAM curves for any set of parameters mapping Equation (7) is atypical, compared to the behavior of traditional perturbed twist mappings such as the standard map. The map lacks bounding KAM curves because the perturbation, $\delta E$, is discontinuous at ϕ values corresponding to collision, violating smoothness conditions necessary for such curves to exist (e.g., Lichtenberg & Lieberman 1983).⁷ A mutual inclination between Planet Nine and would smooth out the discontinuities $\delta E$, possibly allowing for KAM curves to appear. We return to this point below in our discussion in Section 5.

Zoom In Zoom Out Reset image size

We saw in Section 2.1 that including an active Neptune in numerical simulations leads to significantly more chaotic test-particle behavior. When Neptune's influence is approximated only by its quadrupole contribution to the gravitational potential, most of the long-lived test-particles librate stably in MMRs with Planet Nine. Batygin & Morbidelli (2017) find similar stable, librating behavior in their simulations that approximate the effect of all of the giant planets by a quadrupole potential. By contrast, when the full effect of Neptune is taken into account, most of the test particles' semimajor axes diffuse chaotically.

We can apply the same mapping derived above in Equation (7) in order to assess Neptune's contribution to particles' chaotic semimajor axis evolution. Figure 11 show some example chaotic trajectories of particles perturbed solely by Neptune near three different semimajor axis values. The trajectories were computed using Equation (7), but with parameters appropriate to Neptune.⁸ Each panel shows trajectories of points started near the unstable fixed points of resonances up to fifth order. We compute the "chaotic fraction" of the mapping's phase space by dividing the phase space into a 100 × 100 grid and counting the fraction of grid cells that contain at least one point of any trajectory. Given the fractal nature of the chaotic phase space, the measured fraction is sensitive to the size of grid cell used. Nonetheless, the computed chaotic fractions illustrate that the full phase space of the mapping becomes progressively more chaotic with increasing semimajor axis, corresponding to larger values of . This point is further emphasized in Figure 12, where we plot chaos fractions computed for mapping Equation (7) as a function of semimajor axis for different perihelion distances. Clearly, a significant portion of the initial conditions in our fiducial simulation ($q\in [33,50]\,\mathrm{au}$, $a\in [150,550]\,\mathrm{au}$) occupy regions of phase space where the motion is expected to be chaotic under the influence of Neptune alone.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Gravitational interactions with Neptune and Planet Nine affect test particles' angular momentum, as well as their energy. Therefore, test particles' eccentricities and longitudes of pericenter should evolve at the same time as they experience chaotic semimajor axis evolution. In the numerical integrations presented in Section 2, we saw that the chaotically scattering population maintains orbital anti-alignment with Planet Nine to within roughly $\sim \pm 60^\circ$. In Section 3, we showed that, in the secular approximation, a stable anti-aligned libration island exists where torques from Planet Nine and Neptune's averaged gravitational potentials balance to maintain anti-aligned particle orbits. However, these particles' chaotic semimajor evolution invalidates the usual assumptions made when employing a secular approach. Nonetheless, we argue that the evolution of test particles' eccentricities and longitudes of perihelion should, on average, follow the predictions of secular theory. In part, this is because, at high eccentricity, the torques causing changes in angular momentum and longitude of perihelion vary little with semimajor axis. Consequently, the particles' chaotic semimajor axis diffusion has little effect on their angular momentum dynamics. The torques are largely independent of semimajor axis because the gravitational effects of Neptune and Planet Nine are concentrated near pericenter when interactions are strongest. In particular, the average precession rate induced by Neptune reduces to

$$\begin{eqnarray}&&\dot{\varpi }=\displaystyle \frac{3}{4}n{\mu }_{N}{\left(\displaystyle \frac{{a}_{N}}{a(1-{e}^{2})}\right)}^{2}\approx \displaystyle \frac{3}{4}n{\mu }_{N}{\left(\displaystyle \frac{{a}_{N}}{2q}\right)}^{2}\end{eqnarray} \tag{ 8 }$$

for large a at fixed q. According to Equation (8), the average perihelion precession induced by Neptune per orbit depends only on q. Similarly, the per-orbit precession induced by Planet Nine depends only on q for a sufficiently large particle a. Figure 13 compares the increment in ϖ induced by Planet Nine after a single orbit, as measured in N-body simulations of particles with semimajor axes spanning the range $a\in [1000,2000]\,\mathrm{au}$. An analytic prediction that approximates the test-particle orbit as parabolic is also shown (see Appendix B for details). The analytic prediction depends only on the perihelion distance, q, and provides an excellent fit to the N-body simulation results.

Zoom In Zoom Out Reset image size

The near-independence of torques on semimajor axis alone does not fully justify the secular approximation as an explanation for why anti-alignment is maintained. In the usual secular averaging procedure, one performs a double integral over the orbital phases of the test particle and Planet Nine. The validity of this double averaging requires that the orbital phases of the particle and Planet Nine are uncorrelated, e.g., there are no MMRs between Planet Nine and the test particle. The mapping shows that the uncorrelated assumption is valid for the chaotically diffusing particles as well. It is straightforward to show, by making a histogram of the mapping variable ϕ for some chaotic trajectories of the mapping Equation (6), that the distribution of ϕ for these trajectories is nearly uniform, demonstrating that Planet Nine's orbital phase is not correlated with the test particles' phases. Therefore, the double-integral over orbital phase used in the usual secular approximation is equally valid in the case of the chaotically diffusing particles.

In Figure 14, we plot perihelion distance, q, versus orbital alignment ${\rm{\Delta }}\varpi$ for particles with $a\gt 1000\,\mathrm{au}$ from our fiducial simulation. We also plot contours of secular Hamiltonians computed for a = 1500 and 3000 au. These contours are produced in essentially the same way as the secular contours plotted in Figure 3 (see also Beust 2016; Batygin & Morbidelli 2017), except now we plot q versus ${\rm{\Delta }}\varpi$ instead of $e\cos {\rm{\Delta }}\varpi$ and $e\sin {\rm{\Delta }}\varpi$. As argued above, the shape of these contours is insensitive to the choice of a when a is large, so these contours should be representative of the secular trajectories followed by all of the particles plotted in Figure 14. Indeed, we see that the tracks followed by particles in the $q\mbox{--}{\rm{\Delta }}\varpi$ plane closely follow the predicted secular contours. This demonstrates that secular theory can be used to successfully predict the apsidal evolution of these particles despite their significant diffusion in semimajor axis.

Zoom In Zoom Out Reset image size

The low-q points in Figure 14 are not uniformly spaced on their secular trajectories, but rather cluster at q values above the equilibrium point located near $(q,{\rm{\Delta }}\varpi )=(\sim 35\,\mathrm{au},\pi )$. In our simulations, we observe that when particles' perihelion distances dip below $q=30\,\mathrm{au}$, they are quickly ejected by Neptune. Eventually, most of the low-q particles shown in Figure 14 should evolve onto Neptune-crossing orbits and be removed. However, the timescale to do so is longer than our simulation duration and we are left only with particles whose secular evolution has stalled at perihelion distances $q\gt 30\,\mathrm{au}$, causing the observed clustering.

In Section 2.3, we saw that when Planet Nine's mass was reduced to $1{M}_{\oplus }$, the chaotically scattering particles clustered near ${\rm{\Delta }}\varpi =3\pi /2$ (rather than ${\rm{\Delta }}\varpi =\pi$). We argue that this clustering reflects the fact that the particles have not had sufficient time to relax from their initial conditions. Figure 15 plots q versus ${\rm{\Delta }}\varpi$ for particles both interior and exterior to Planet Nine in the $1{M}_{\oplus }$ Planet Nine simulation. Representative secular Hamiltonian contours are plotted as well. In order to demonstrate that secular evolution has not had sufficient time to randomize particle alignments, we integrate the equations of motion of the secularly averaged Hamiltonian at some representative values of particle semimajor axes. At a given semimajor axis, we initialize a set of secular trajectories on contours contained within the anti-aligned libration region. These trajectories are initialized at their minimum q values, which occur at ${\rm{\Delta }}\varpi =\pi$. The trajectories' evolution proceeds counter-clockwise about the secular contours, and we plot the final q and ${\rm{\Delta }}\varpi$ after 5 Gyr as red lines in Figure 15. In the left panel of that figure, we see that the secular trajectories of particles interior to Planet Nine complete approximately ∼1/4 of a precession cycle over the course of 5 Gyr, leaving them with ${\rm{\Delta }}\varpi \sim 3\pi /2$.⁹ Because at least one secular precession cycle must occur for initial ${\rm{\Delta }}\varpi$ values to be effectively randomized, at least $\gtrsim 20\,\mathrm{Gyr}$ of evolution would be necessary for particles to relax to an equilibrium distribution. Trajectories exterior to Planet Nine shown in the right panel experience an even smaller fraction of their secular precession cycle. While the apsidal evolution experienced by the N-body particles plotted in Figure 15 is more complicated than these simple secular trajectories because of semimajor axis diffusion, the N-body particles nonetheless fall within a region of the $(q,{\rm{\Delta }}\varpi )$ plane similar to the secular trajectories.

Zoom In Zoom Out Reset image size

The clustering is not just an artifact of initializing test particles in exactly (anti-)aligned configurations. We observe a similar clustering in simulations initialized from randomized initial alignments, ϖ. In these randomized simulations, particles precess relatively quickly through the bottom portion of the secular trajectories plotted in Figure 15, where their precession is driven by mainly by Neptune. Their evolution slows dramatically upon reaching ${\rm{\Delta }}\varpi \sim 3\pi /2$, where their qs begin to increase and the rate of Neptune-induced precession consequently decreases (Equation (8)). Thus, even though the particles are initially randomized in ϖ, they have not relaxed from their initial narrow range of qs, which place them near the bottom of the secular trajectories plotted in Figure 15.

Our simulations show that a $1{M}_{\oplus }$ Planet Nine is capable of inducing clustering among ESDOs via a dynamical mechanism very different from the one operating in our simulations with a $10{M}_{\oplus }$ Planet Nine. Taken at face value, this poses a serious impediment to deriving indirect dynamical constraints on Planet Nine's mass and orbit because there is no obvious way to observationally distinguish which mechanism is responsible for the ESDOs' apsidal clustering. However, it remains to be seen whether the $1{M}_{\oplus }$ Planet Nine clustering mechanism operates effectively in more realistic simulations allowing for mutual inclinations. The results of BB16b, who disfavor such a small Planet Nine mass based on an extensive grid of N-body simulations, suggest that inclinations probably render the $1\,{M}_{\oplus }$ Planet Nine clustering mechanism less effective.

To construct the secular Hamiltonian used in Section 3, we also include the quadrupole contribution of Neptune,

$$\begin{eqnarray}&&{H}_{{\rm{N}},\sec }=-\displaystyle \frac{{\mu }_{N}{a}_{N}^{2}}{4{{\rm{\Lambda }}}^{3}{G}^{3}},\end{eqnarray} \tag{ 11 }$$

so that the full Hamiltonian is $H({\rm{\Lambda }},\lambda ,L,l,G,g)\,=\,{H}_{\mathrm{three} \mbox{-} \mathrm{body}}\,+{H}_{{\rm{N}},\sec }$.

The secular Hamiltonian ${\bar{H}}_{\sec }$ is computed by averaging over the orbital phases of both Planet Nine and the test particle and is given by the double integral

$$\begin{eqnarray}&&{\bar{H}}_{\sec }(G,g)=\displaystyle \frac{1}{4{\pi }^{2}}{\int }_{0}^{2\pi }{\int }_{0}^{2\pi }H({\rm{\Lambda }},\lambda ,L,l,G,g){dld}\lambda .\end{eqnarray} \tag{ 12 }$$

After dropping constant terms, which do not enter the equation of motion, Equation (12) can be rewritten as

$$\begin{eqnarray}{\bar{H}}_{\sec }(G,g) & = & -\displaystyle \frac{{\mu }_{N}{a}_{N}^{2}}{4{{aG}}^{3}}-\displaystyle \frac{{\mu }_{9}}{4{\pi }^{2}}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{2\pi }\\ & & \times \displaystyle \frac{(1-e\cos u)(1-{e}_{9}\cos {u}_{9})}{| {r}_{9}-r| }{{dudu}}_{9}.\end{eqnarray} \tag{ 13 }$$

where u and u₉ are the eccentric anomalies of the particle and Planet Nine, respectively.¹¹ We numerically evaluate the double integral in Equation (13).

In order to study the dynamics of an N:k resonance between a test particle and Planet Nine, we introduce new canonical coordinate-momentum pairs

$$\begin{eqnarray}&&\sigma =l-\lambda ;\quad {\rm{\Sigma }}=\displaystyle \frac{{NL}+k{\rm{\Lambda }}}{N-k}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\phi =\displaystyle \frac{1}{k}(N\lambda -{kl}-(N-k)g);\quad {\rm{\Phi }}=\displaystyle \frac{{\rm{\Lambda }}+L}{N-k}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\gamma =g;\quad {\rm{\Gamma }}=G+{\rm{\Lambda }}+L.\end{eqnarray} \tag{ 16 }$$

We can construct a model with two degrees of freedom for the resonant dynamics by averaging over σ, yielding the new Hamiltonian

$$\begin{eqnarray}&&{\bar{H}}_{\mathrm{res}.}({\rm{\Phi }},\phi ,{\rm{\Gamma }},\gamma )=-\displaystyle \frac{1}{2}{({\rm{\Sigma }}-N{\rm{\Phi }})}^{-2}+{{kn}}_{9}{\rm{\Phi }}-{\mu }_{9}\bar{R},\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&\bar{R}=\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{Rd}\sigma .\end{eqnarray} \tag{ 18 }$$

The Hamiltonian of Equation (17) represents a two-degree-of-freedom system that depends parametrically on the conserved quantity Σ. In Section 4.3, we have plotted contours of constant ${\bar{H}}_{\mathrm{res}}$ while keeping Γ and $\gamma =\pi$ fixed.