Dynamical Evolution Induced by Planet Nine

Over the course of the 170 years that followed the successful prediction of Neptune (Adams 1846; Le Verrier 1846), the presence of additional planets within the solar system has been contemplated by an extensive list of astronomers and celestial mechanicians (Hoyt 1980). Historically, the lines of evidence for the existence of distant, massive bodies that orbit the Sun have ranged from the (apparently) anomalous motion of Uranus (Pickering & Pickering 1909; Lowell 1915) and the unexpected orbital characteristics of long-period comets (Forbes 1880; Matese & Whitmire 1986), to the peculiar structure of the solar system's small body populations (Brunini & Melita 2002; Gladman & Chan 2006; Gomes et al. 2006; Lykawka & Mukai 2008; Trujillo & Sheppard 2014; Volk & Malhotra 2017). The predicted physical and orbital properties of the putative planets have been equally as varied, with inferred masses and semimajor axes spanning the Mars–Jupiter range and tens to thousands of astronomical units, respectively.

A recent addition to the aggregate of planetary predictions within the solar system is the Planet Nine hypothesis³ (Batygin & Brown 2016a). Within the framework of this model, the observed orbital clustering of $a\gtrsim 250\,\mathrm{au}$ Kuiper Belt objects (KBOs; Figure 1) is sculpted by an $m\sim 10\,{m}_{\oplus }$ planet residing on an appreciably eccentric ($e\sim 0.3\mbox{--}0.7$), large semimajor axis ($a\sim 300\mbox{--}700\,\mathrm{au}$) orbit, whose plane roughly coincides with the plane of the distant bodies, and is characterized by a perihelion direction that is anti-aligned with respect to the average apsidal orientation of the KBOs. In addition to (i) facilitating the orbital confinement of the aforementioned population of long-period KBOs and (ii) providing a physical mechanism for the perihelion detachment of Sedna-type orbits from Neptune, the presence of Planet Nine entails a series of additional consequences for the observed structure of the solar system (Brown & Batygin 2016). In particular, it has been shown that the dynamical effects of Planet Nine naturally explain (iii) the existence of highly inclined, large semimajor axis Centaurs (Gomes et al. 2015; Batygin & Brown 2016a), (iv) the six-degree obliquity of the Sun (Bailey et al. 2016; Lai 2016; Gomes et al. 2017), as well as (v) the origins of proximate ($a\lt 100\,\mathrm{au}$) retrograde KBOs (Batygin & Brown 2016b).

Zoom In Zoom Out Reset image size

A key characteristic that differentiates the various planetary proposals is the dynamical mechanism through which the envisaged planet generates its observational signatures. In this regard, the Planet Nine (P9) hypothesis remains incomplete. Although numerical simulations reveal that synthetic models of the solar system that include Planet Nine can provide a good match to the observational data (see, however, Nesvorny et al. 2017), the dynamical process through which the physical confinement of the distant orbits occurs remains poorly understood.

The original study of Batygin & Brown (2016a) suggested that mean-motion resonances (MMRs; including those of high order) are responsible for orbital clustering. Expanding on this idea, Malhotra et al. (2016) suggested that each of the distant KBOs are currently trapped in N:1 and N:2 MMRs with Planet Nine. More recently, Millholland & Laughlin (2017) and Becker et al. (2017) carried out a large-scale numerical exploration of the resonant hypothesis in relation to the existing data, further demonstrating its viability. Meanwhile, Beust (2016) showed that a purely secular treatment of the dynamics provides a good match to the simulation results of Batygin & Brown (2016a), suggesting that resonant dynamics may be irrelevant to the problem at hand. In light of these conflicting results, the dominant P9–KBO interaction mechanism remains elusive, and its identification is the primary purpose of this paper.

An analytical characterization of P9–KBO coupling is important for three reasons, the first being falsifiability. The specific perturbation mechanism (along with the data it aims to explain) draws the distinctions among the various planetary hypotheses and can be used to refute a specific model upon confrontation with the data. As an example, suppose that an observational survey discovers a planet at a radial separation of a few hundred astronomical units, but the gravitational effects of this planet do not facilitate a physical confinement among the distant KBO orbits through the envisioned dynamical process.⁴ Such a discovery would imply that the Planet Nine hypothesis, as formulated, is incorrect.

A second, more practical motivation for studying P9–KBO interactions is the link between existing observations and the predicted orbit of Planet Nine. That is, if all distant KBOs are presently trapped in MMRs with P9, their mean longitudes may contain information about the location of Planet Nine in its orbit (Malhotra et al. 2016; Millholland & Laughlin 2017). If, on the other hand, the interactions are purely secular, such a connection cannot be established.⁵ Therefore, understanding the physics of P9–KBO coupling can offer a more complete interpretation of the extant observational data set and may yield an avenue toward further constraining the astronomical search forPlanet Nine.

The third, and final, reason for this study is purely academic. While the exploration of the circular restricted three-body problem (as a paradigm for interactions between planets and small bodies in the solar system) dates back multiple centuries (Marquis de Laplace 1799; Poincaré 1902), its highly elliptical counterpart remains much more scarcely understood (see, e.g., Beaugé et al. 2006; Michtchenko et al. 2006; Pichierri et al. 2017). Indeed, a well-formulated theory for gravitational coupling between a planet and a test particle in the severe orbit-crossing regime does not currently exist. As a result, a perturbative study of P9 dynamics can yield interesting insights into the general mathematical structure of interactions among highly excited orbits. Such a framework would have considerable applications beyond the problem at hand, including the characterization of the remarkable dynamical states of highly eccentric, resonant exoplanets (Lee 2004; Tan et al. 2013).

Before delving into calculations, we delineate a list of specific questions we wish to answer in this work.

• 1.  
  What role (if any) do resonant interactions play within the dynamical evolution induced by Planet Nine? If resonances are prevalent, what order/multiplet harmonics dominate the dynamics, and what are their characteristic widths?
• 2.  
  What role (if any) do secular interactions play within the dynamical evolution induced by Planet Nine? If dominant, how are close encounters avoided on nearly coplanar, anti-aligned orbits? Moreover, if resonant interactions are relevant to the Planet Nine hypothesis, why does the purely secular phase-space portrait provide a good match to the results of numerical simulations?
• 3.  
  What parameters determine the critical semimajor axis corresponding to the transition between randomized and clustered longitudes of perihelion? What physical effect controls this transition?
• 4.  
  What is the qualitative behavior of inclination dynamics within the framework of P9-driven evolution? What dynamical process allows some of the objects to acquire exceptionally high inclinations in the distant Kuiper Belt?

We take these questions as an approximate guide to the logic of this paper, which is structured as follows. In Section 2, we revisit a purely secular description of P9-induced orbital evolution, carrying out the averaging procedure in closed form. In Section 3, we consider a pair of idealized N-body simulations and outline the key differences between numerical experiments and pure secular theory. In Section 4, we present a semi-analytical theory of resonant P9–KBO interactions and elucidate secular dynamics embedded within MMRs as the primary driver of apsidal confinement of distant KBOs. Subsequently, we discuss the onset of large-scale inclination oscillations of long-period bodies in Section 5. We then re-examine the existing observational data in light of the aforementioned theoretical developments in Section 6. We summarize and discuss the implications of our results in Section 7. A brief analysis of long-term stable, but apsidally unconfined, orbits is presented in Appendix B.

We begin our study of P9-induced dynamics within the framework of purely secular perturbation theory. In Batygin & Brown (2016a), our preliminary exploration of secular dynamics relied on an octupole-order expansion of the gravitational potential (Kaula 1962; Mardling 2013), which implicitly assumed that the orbits under consideration do not cross. Following Beust (2016), here we abandon the series-expansion approach to modeling P9–KBO coupling and carry out the phase-averaging procedure in closed form.

As in Batygin & Brown (2016a) and Beust (2016), we assume coplanarity and model the mean-field effects of the known giant planets, ignoring their eccentricities. Further, we choose to work in a slowly rotating coordinate frame that co-precesses with Planet Nine's perihelion (the corresponding contact transformation is spelled out in Batygin & Brown 2016a). The governing (doubly averaged) Hamiltonian of a test particle under planetary perturbations then has the form

$$\begin{eqnarray}{\overline{\overline{{ \mathcal H }}}}_{{\rm{s}}} & = & -\displaystyle \frac{1}{4}\displaystyle \frac{{ \mathcal G }\,{M}_{\odot }}{a}\displaystyle \frac{1}{{(1-{e}^{2})}^{3/2}}\displaystyle \sum _{j=5}^{8}\displaystyle \frac{{m}_{j}\,{a}_{j}^{2}}{{M}_{\odot }\,{a}^{2}}\\ & & +{\dot{\varpi }}_{9}\,\sqrt{{ \mathcal G }\,{M}_{\odot }a}(1-\sqrt{1-{e}^{2}})\\ & & -\displaystyle \frac{1}{4{\pi }^{2}}\oint \oint \displaystyle \frac{{ \mathcal G }\,{m}_{9}}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}| }\,d\lambda \,d{\lambda }_{9},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal G }$ is the gravitational constant, ${M}_{\odot }$ is the mass of the Sun, ${\boldsymbol{r}}$ is the position vector, λ is the mean longitude, a is semimajor axis, e is eccentricity, and ϖ is the longitude of perihelion. All quantities pertaining to the four canonical giant planets are labeled with indexes 5–8, while the values corresponding to Planet Nine are denoted with the subscript 9. The unlabeled variables correspond to the KBO. Finally, Planet Nine's orbit-averaged precession rate is given by the expression⁶

$$\begin{eqnarray}&&{\dot{\varpi }}_{9}=\displaystyle \frac{3}{4}\sqrt{\displaystyle \frac{{ \mathcal G }\,{M}_{\odot }}{{a}_{9}^{3}}}\displaystyle \frac{1}{{(1-{e}_{9}^{2})}^{2}}\displaystyle \sum _{j=5}^{8}\displaystyle \frac{{m}_{j}\,{a}_{j}^{2}}{{M}_{\odot }\,{a}_{9}^{2}}.\end{eqnarray} \tag{ 2 }$$

The three terms present in Equation (1) have simple physical interpretations. The first term governs the secular advance of the KBO's perihelion due to the phase-averaged gravitational potential of Jupiter, Saturn, Uranus, and Neptune (in direct parallel with Equation (2)). The second term accounts for the fact that the reference frame is slowly rotating. The third term governs secular P9–KBO interactions. Note that the indirect part of the disturbing potential is by default entirely averaged out within the framework of secular theory and need not be accounted for (Morbidelli 2002).

The only pair of dynamical variables present in Hamiltonian (1) is ($e,{\rm{\Delta }}\varpi$), meaning that the system is integrable. In other words, contours of the numerically averaged function (1) fully encapsulate the accompanying orbital evolution. Figure 2 shows the secular phase-space portraits of the system, projected onto $e\mbox{--}{\rm{\Delta }}\varpi$ space, for $a=50,150,\,\ldots ,\,550\,\mathrm{au}$, adopting P9 parameters from Batygin & Brown (2016a; specifically, ${a}_{9}=700\,\mathrm{au}$, e₉ = 0.6, and ${m}_{9}=10\,{m}_{\oplus }$). This figure can be readily compared with Figure 4 of Batygin & Brown (2016a) and confirms that the purely secular portrait provides a good match to the numerically computed portraits in the same dynamical regime (Beust 2016).

Zoom In Zoom Out Reset image size

In all panels denoted by $a\geqslant 250\,\mathrm{au}$, the $e\mbox{--}{\rm{\Delta }}\varpi$ space is characterized by two stable equilibrium points: one at ${\rm{\Delta }}\varpi =0$ and another at ${\rm{\Delta }}\varpi =180\,\deg$. In each diagram, the two libration regions surrounding these fixed points are separated by a solid curve, which corresponds to a tangential configuration of the KBO and P9 orbits. Note further that on these tangential contact curves, the derivatives of the Hamiltonian are discontinuous (i.e., ${\overline{\overline{{ \mathcal H }}}}_{{\rm{s}}}$ is locally class ${{ \mathcal C }}^{0}$), signaling a breakdown of the secular framework (Gronchi 2002).

A simple examination of the contours of ${\overline{\overline{{ \mathcal H }}}}_{{\rm{s}}}$ depicted in Figure 2, reveals that KBOs initialized on nearly Neptune-crossing orbits (immediately below the horizontal lines labeled $q={a}_{8}$) will suffer drastically different secular evolutions depending on their starting value of ${\rm{\Delta }}\varpi$. Bodies initially close to ${\rm{\Delta }}\varpi \sim 0$ will be driven onto the tangential orbit-crossing curves through apsidal precession and will eventually be removed from the system by recurrent close encounters with Planet Nine. On the other hand, objects initialized at ${\rm{\Delta }}\varpi \sim \pi$ settle onto secular trajectories that encircle that anti-aligned equilibrium and never encounter the tangential collision curves. Thus, the dynamical lifetimes of apsidally anti-aligned objects can be naïvely envisioned to be longer than their aligned counterparts, and in time, a cluster of exclusively anti-aligned orbits should be carved out by Planet Nine.

We note further that there is an island of stable apsidal libration around ${\rm{\Delta }}\varpi =0$ that avoids crossing the tangential configuration curve and is thus protected from close encounters. However, the eccentricities along these protected orbits are moderate, and the corresponding perihelion distances are large. This means that even if KBOs somehow came to occupy these islands of stability, they would be difficult to detect observationally.

With this picture in mind, it is tempting to affirm that the agreement between theory and simulation is satisfactory, and proceed forward within the purely secular framework. This is, however, a misconception, facilitated by the apparent smoothness of the secular phase-space portraits shown in Figure 2. In particular, the fact that the contours of ${\overline{\overline{{ \mathcal H }}}}_{{\rm{s}}}$ do not show any kinks within the apsidally anti-aligned domain is simply a consequence of the integrability of the non-tangential singularity (Thomas & Morbidelli 1996; Gronchi & Milani 1998) and does not mean that the system can elude collisions. Instead, recalling that the physical setup of the orbits is planar, it is trivial to demonstrate that under the assumption of uncorrelated Keplerian motion, all objects entrained in the apsidally anti-aligned configuration with Planet Nine would suffer close encounters on timescales much shorter than the age of the solar system. Therefore, despite giving the illusion of agreement with N-body simulations, pure secular theory predicts that the entire distant Kuiper Belt is dynamically unstable and should have been cleared out on a timescale comparable to the orbital precession time.

Contrary to this view, published numerical experiments reveal that particles residing deep within the cores of the anti-aligned libration regions of Figure 2 remain stable over the multi-Gyr lifetime of the solar system (Batygin & Brown 2016a; Brown & Batygin 2016). This disparity suggests that secular theory alone is unlikely to represent a full dynamical description of P9-driven evolution and that an additional stabilizing mechanism is at play in the simulations. Let us now examine this point further.

To quantify the discrepancy between published simulations and pure secular theory, we carry out a pair of simplified numerical experiments that mirror those reported in Batygin & Brown (2016a) and Brown & Batygin (2016). In particular, we evolve an initially axisymmetric disk of 6000 eccentric test particles with $a\in (150,750)$ au and $q\in (30,36)$ au under the influence of the phase-averaged potential of the four inner giants as well as Planet Nine with $({a}_{9},{e}_{9})=(700\,\mathrm{au},0.6)$ and $({a}_{9},{e}_{9})=(600\,\mathrm{au},0.5)$. The test particles are initialized with a null vertical velocity dispersion and thus remain confined to Planet Nine's orbital plane throughout the simulation.

Unlike the simulation suite of Batygin & Brown (2016a), Brown & Batygin (2016), Millholland & Laughlin (2017), and Becker et al. (2017), where Neptune was modeled directly, or that of Batygin & Brown (2016b), where the Keplerian motion of all four inner giants was resolved, here we emulate the effects of Jupiter, Saturn, Uranus, and Neptune with an effective quadrupolar gravitational moment of the Sun,

$$\begin{eqnarray}&&{J}_{2}=\left(\displaystyle \frac{3{\cos }^{2}({i}_{9})-1}{4}\right)\displaystyle \sum _{j=5}^{8}\displaystyle \frac{{m}_{j}\,{a}_{j}^{2}}{{M}_{\odot }\,{{ \mathcal R }}^{2}},\end{eqnarray} \tag{ 3 }$$

setting the inner absorbing radius to ${ \mathcal R }=20\,\mathrm{au}$. This idealization is employed specifically to avoid contaminating P9-induced KBO evolution with chaotic dynamics that arise from scattering off Neptune, and to yield the closest point of comparison between numerical and semi-analytical results. Moreover, here we will assume that the inclination of Planet Nine relative to the Laplace plane of the canonical giant planets is sufficiently small to approximate $\cos ({i}_{9})\approx 1$, while keeping in mind that $\hat{z}\,$axis of our coordinate system coincides with the orbital plane of Planet Nine.

To carry out the simulations, we utilized the mercury6 gravitational dynamics software package (Chambers 1999). The integrations were performed using the hybrid Wisdom–Holman/Bulirsch–Stoer algorithm (Wisdom & Holman 1992; Press et al. 1992) with a time step of ${\rm{\Delta }}t$ = 3100 days (i.e., 1/10th of Uranus' orbital period), and spanned 4 Gyr. Any particle that attained a radial distance of $r\lt { \mathcal R }$ or r > 10,000 au was removed from the simulations.

Within the context of these idealized numerical experiments, all long-term stable particles retain nearly constant semimajor axes throughout the integration, and Figure 3 depicts an orbital histogram of the surviving bodies. Bins shown in black correspond to apsidally confined objects (defined by $| {\rm{\Delta }}\varpi -180\,\deg | \leqslant 90\,\deg$ throughout the simulation), while those shown in gray denote particles that are long-term stable but experience apsidal circulation. We note that here, the gray bins are stacked on top of the black bins, such that the relative size of the gray and black components of each column is a measure of the contamination of the apsidally confined population of objects (at a given semimajor axis) by those undergoing perihelion circulation.

Zoom In Zoom Out Reset image size

A strong propensity toward apsidal clustering for $P/{P}_{9}\gtrsim 1/3$ is clearly evident in both simulations. However, it is also important to note that the same orbital periods (e.g., those corresponding to the 3:1, 2:1, and 1:1 commensurabilities in Figure 3) can be simultaneously occupied by apsidally circulating and librating orbits. This means that even within the framework of highly idealized treatment of P9-induced dynamics, the clustering of the longitude of perihelion beyond a critical semimajor axis is not perfectly strict, and the existence of objects that do not conform to the general anti-aligned orbital pattern is an expected consequence of the model (Shankman et al. 2017).

In addition to the already established tendency toward orbital clustering with increasing orbital period, it can be clearly seen in Figure 3 that all long-lived objects that exhibit perihelion confinement have semimajor axes that correspond to mean-motion commensurabilities with Planet Nine. Particularly, the 1:1, 3:2, 2:1, and 5:2 resonances contain the largest populations of apsidally clustered KBOs in both simulations.⁷ This point undercuts a key difference between the results of N-body simulations and pure secular theory, and demonstrates that (at least within the framework of a planar physical setup) distant KBOs derive their long-term orbital stability from the phase-protection mechanism inherent to MMRs.

The typical dynamical evolutions of bodies trapped in the four aforementioned resonances are shown in Figure 4. Specifically, the left and middle panels depict the evolutions of orbital eccentricities and apsidal lines. On timescales of order ∼0.1–1 Gyr, orbital eccentricities experience considerable oscillations in concert with librations of ${\rm{\Delta }}\varpi$. This facilitates a periodic regression of the KBO perihelion distance, generating dynamically detached (Sedna-type) orbits (Gladman et al. 2002; Brown et al. 2004).

Zoom In Zoom Out Reset image size

Of particular importance is the behavior of the resonant angles,

$$\begin{eqnarray}{\psi }_{\mathrm{res}} & = & k\,{\lambda }_{9}-{\ell }\,\lambda -(k-{\ell })\varpi ,\\ {\phi }_{\mathrm{res}} & = & k\,{\lambda }_{9}-{\ell }\,\lambda -(k-{\ell }){\varpi }_{9},\end{eqnarray} \tag{ 4 }$$

where $k$ and ${\ell }$ are integers (note that for $| k-{\ell }| \gt 1$, additional resonant harmonics related to these angles through ${\rm{\Delta }}\varpi$ exist). In the right panels of Figure 4, ${\psi }_{\mathrm{res}}$ is plotted in green and ${\phi }_{\mathrm{res}}$ is plotted in black. The actual librations of ${\phi }_{\mathrm{res}}$ and ${\psi }_{\mathrm{res}}$ are short periodic (as portrayed by the inset within the top-right panel of Figure 4), and their amplitudes correspond to the thickness of the green and black "bands." Meanwhile, the long-periodic oscillations of these angles are mere reflections of the librations of ${\rm{\Delta }}\varpi$, which modulate the locations of the resonant equilibria associated with ${\phi }_{\mathrm{res}}$ and ${\psi }_{\mathrm{res}}$.

The fact that the amplitude of long-periodic oscillations of ${\phi }_{\mathrm{res}}$ is much smaller than that of ${\psi }_{\mathrm{res}}$ for orbits that exhibit libration in ${\rm{\Delta }}\varpi$ demonstrates that the resonant multiplets containing ${\phi }_{\mathrm{res}}$ represent a better approximation of the real dynamics than those containing ${\psi }_{\mathrm{res}}$. In particular, if the Hamiltonian depended only on ${\phi }_{\mathrm{res}}$, this angle would have no long-periodic oscillations while ${\psi }_{\mathrm{res}}$ (being equal to ${\phi }_{\mathrm{res}}+({\ell }-k){\rm{\Delta }}\varpi$) would oscillate with the amplitude and period of $| ({\ell }-k)| \,{\rm{\Delta }}\varpi$. The opposite would be true if the Hamiltonian depended only on ${\psi }_{\mathrm{res}}$. Thus, the relative amplitudes of the long-periodic oscillations of these angles are inversely correlated to the relative strengths of the corresponding terms. Accordingly, in our analytic approach in the next section, we will consider ${\phi }_{\mathrm{res}}$ as a reference angle, in contrast to Malhotra et al. (2016) and Beust (2016), who considered ${\psi }_{\mathrm{res}}$ instead.

Containing only the longitude of perihelion of Planet Nine and not that of the KBOs themselves, the (short-periodic) libration of ${\phi }_{\mathrm{res}}$ drives the oscillation of the particle's semimajor axis but does not affect the evolution of its eccentricity. Thus, the libration in ${\phi }_{\mathrm{res}}$, sometimes referred to as "corotation" resonance, already implies a certain disconnect among the degrees of freedom related to the particle's semimajor axis and its eccentricity. Moreover, the striking separation of timescales associated with the resonant (${\phi }_{\mathrm{res}}$ − a) dynamics and secular ($e\mbox{--}{\rm{\Delta }}\varpi$) dynamics motivates the construction of the semi-analytical model that will follow, based on the adiabatic approach.

The dominant dependence of the resonant dynamics on ${\phi }_{\mathrm{res}}$, as opposed to another harmonic that contains ϖ in its critical argument, is central to maintaining the apsidally anti-aligned libration of the orbits. In fact, long-term stable particles whose resonant dynamics are driven by any angle other than ${\phi }_{\mathrm{res}}$ exhibit circulation of their longitude of perihelion, and correspond to objects denoted by the gray bins in Figure 3. Although characterizing the dynamics in this transitionary semimajor axis range is not the primary purpose of this study, we present a brief analysis of this mode of orbital evolution in Appendix B.

As shown in the previous section, all surviving members of the synthetic scattered disk that exhibit persistent apsidal anti-alignment with Planet Nine's orbit are locked into MMRs with P9. The purpose of the following analysis is thus to explain this behavior on semi-analytic grounds. To achieve this goal, we consider the isolated resonant behavior first. As in Section 2, we will abandon traditional methods of celestial mechanics based on series expansions and cast our analysis of the governing Hamiltonian in closed form.

4.1. Mean Motion Resonances

In order to carry out the averaging process numerically and elucidate the ${\phi }_{\mathrm{res}}$–a resonant dynamics, we set up a rigorous Hamiltonian approach to the problem. As is typical for the planar restricted three-body problem, the canonical Poincaré action-angle variables for the test particle are (Murray & Dermott 1999)

$$\begin{eqnarray}{\rm{\Lambda }} & = & \sqrt{{ \mathcal G }{M}_{\odot }\,a}\qquad \lambda ={ \mathcal M }+\varpi \\ {\rm{\Gamma }} & = & \sqrt{{ \mathcal G }{M}_{\odot }\,a}(1-\sqrt{1-{e}^{2}})\qquad \gamma =-\varpi ,\end{eqnarray} \tag{ 5 }$$

where ${ \mathcal M }$ is the mean anomaly. In terms of these variables, the Hamiltonian of the problem is (Morbidelli 2002)

$$\begin{eqnarray}&&{{ \mathcal H }}_{{\rm{r}}}=-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{ \mathcal G }{M}_{\odot }}{{\rm{\Lambda }}}\right)}^{2}-{ \mathcal G }\,{m}_{9}\left(\displaystyle \frac{1}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}| }-\displaystyle \frac{{\boldsymbol{r}}\cdot {{\boldsymbol{r}}}_{9}}{| {{\boldsymbol{r}}}_{9}{| }^{3}}\right).\end{eqnarray} \tag{ 6 }$$

The vector ${\boldsymbol{r}}$ can be written as a function of the canonical variables (5). Similarly, ${{\boldsymbol{r}}}_{9}$ can be expressed as a function of the parameters ${a}_{9},{e}_{9}$, which are assumed to be fixed in time, as well as the angles ${\lambda }_{9}$ and ${\varpi }_{9}$, which are assumed to advance with fixed frequencies, n₉ and ${\dot{\varpi }}_{9}$, respectively. Thus, the Hamiltonian (6) depends on time through these two planetary angles.

In order to remove the explicit time dependence from Equation (6), we extend the phase space by two degrees of freedom. That is, we consider ${\lambda }_{9}$ and ${\gamma }_{9}=-{\varpi }_{9}$ to be independent variables, with conjugated actions ${{\rm{\Lambda }}}_{9}$ and ${{\rm{\Gamma }}}_{9}$. Correspondingly, we add the term $({n}_{9}\,{{\rm{\Lambda }}}_{9}-{\dot{\varpi }}_{9}\,{{\rm{\Gamma }}}_{9})$ to expression (6), making the Hamiltonian autonomous. Of course, in doing so, we do not alter the dynamics in any way, since $d{\lambda }_{9}/{dt}=\partial {{ \mathcal H }}_{{\rm{r}}}/\partial {{\rm{\Lambda }}}_{9}={n}_{9}$ and $d{\varpi }_{9}/{dt}=-d{\gamma }_{9}/{dt}=-\partial {{ \mathcal H }}_{{\rm{r}}}/\partial {{\rm{\Gamma }}}_{9}={\dot{\varpi }}_{9}$.

We now denote the mean motion of the particle by $n={{ \mathcal G }}^{2}\,{M}_{\odot }^{2}/{{\rm{\Lambda }}}^{3}$ and assume that there is a resonance of the kind $k\,{n}_{9}-{\ell }\,n=0$ for some integers $k$ and ${\ell }$. It is then appropriate to make the following canonical transformation:⁸

$$\begin{eqnarray}{\rm{\Lambda }}^{\prime} ={\rm{\Lambda }}/{\ell } & {\phi }_{\mathrm{res}}^{{\prime} }=-{\phi }_{\mathrm{res}}\\ {{\rm{\Lambda }}}_{9}^{{\prime} }={{\rm{\Lambda }}}_{9}+(k/{\ell }) & {\rm{\Lambda }}{\lambda }_{9}^{{\prime} }={\lambda }_{9}\\ {\rm{\Gamma }}^{\prime} ={\rm{\Gamma }} & {\rm{\Delta }}\gamma =\gamma -{\gamma }_{9}=-{\rm{\Delta }}\varpi \\ {{\rm{\Gamma }}}_{9}^{{\prime} }={{\rm{\Gamma }}}_{9}+{\rm{\Gamma }}+((k-{\ell })/{\ell }){\rm{\Lambda }} & {\gamma }_{9}^{{\prime} }={\gamma }_{9}.\end{eqnarray} \tag{ 7 }$$

The unperturbed part of the Hamiltonian (the term independent of m₉) becomes

$$\begin{eqnarray}{({{ \mathcal H }}_{{\rm{r}}})}_{\mathrm{kep}} & = & -\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{ \mathcal G }{M}_{\odot }}{{\rm{\Lambda }}}\right)}^{2}+{n}_{9}({{\rm{\Lambda }}}_{9}^{{\prime} }-k{\rm{\Lambda }}^{\prime} )\\ & & -{\dot{\varpi }}_{9}({{\rm{\Gamma }}}_{9}^{{\prime} }-{\rm{\Gamma }}-(k-{\ell }){\rm{\Lambda }}^{\prime} ).\end{eqnarray} \tag{ 8 }$$

Because $\partial {({{ \mathcal H }}_{{\rm{r}}})}_{\mathrm{kep}}/\partial {\rm{\Lambda }}^{\prime} ={\ell }\,n-k\,{n}_{9}=0$, there is only one fast angle in the Hamiltonian, and it is ${\lambda }_{9}^{{\prime} }$. Thus, the perturbation can be averaged over this angle, leading to

$$\begin{eqnarray}&&{({\overline{{ \mathcal H }}}_{{\rm{r}}})}_{\mathrm{pert}}=-\displaystyle \frac{{ \mathcal G }\,{m}_{9}}{2\,\pi \,{\ell }}{\int }_{0}^{2\pi {\ell }}\left(\displaystyle \frac{1}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}| }-\displaystyle \frac{{\boldsymbol{r}}\cdot {{\boldsymbol{r}}}_{9}}{| {{\boldsymbol{r}}}_{9}{| }^{3}}\right)d{\lambda }_{9}.\end{eqnarray} \tag{ 9 }$$

It is easy to see from d'Alembert's rules that the averaged Hamiltonian can depend only on two angles, ${\phi }_{\mathrm{res}}^{{\prime} }$ and ${\rm{\Delta }}\gamma$. Thus, ${{\rm{\Lambda }}}_{9}^{{\prime} }$ and ${{\rm{\Gamma }}}_{9}^{{\prime} }$ are constants of motion and can be dropped. The Hamiltonian now comprises a two degree of freedom system and is not integrable. However, we note that the degree of freedom in (${\rm{\Gamma }}^{\prime} ,{\rm{\Delta }}\gamma$) is characterized by a frequency of order $\propto {m}_{9}/{M}_{\odot }$ and is slow relative to the other degree of freedom in (${\rm{\Lambda }}^{\prime} ,{\phi }_{\mathrm{res}}^{{\prime} }$), whose frequency is of the order of $\propto \sqrt{{m}_{9}/{M}_{\odot }}$ (Henrard & Caranicolas 1990).

Taking advantage of the aforementioned separation of timescales between the resonant and secular dynamics, we may evaluate the phase-space portrait associated with the resonant Hamiltonian in the adiabatic approximation by freezing the evolution of the KBO's eccentricity and apsidal line relative to the major axis of Planet Nine (Wisdom 1985). In particular, we compute the function

$$\begin{eqnarray}{\overline{{ \mathcal H }}}_{{\rm{r}}} & = & -\displaystyle \frac{{ \mathcal G }{M}_{\odot }}{2\,a}-{n}_{9}\left(\displaystyle \frac{k}{{\ell }}\right)\sqrt{{ \mathcal G }\,{M}_{\odot }\,a}\\ & & -{\left.\displaystyle \frac{{ \mathcal G }{m}_{9}}{2\pi {\ell }}{\displaystyle \int }_{0}^{2\pi {\ell }}\left(\displaystyle \frac{1}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}| }-\displaystyle \frac{{\boldsymbol{r}}\cdot {{\boldsymbol{r}}}_{9}}{| {{\boldsymbol{r}}}_{9}{| }^{3}}\right)\right|}_{{\phi }_{\mathrm{res}}}\,d{\lambda }_{9}\end{eqnarray} \tag{ 10 }$$

on a $({\phi }_{\mathrm{res}}-a)$ grid, setting the quantities Γ and ${\rm{\Delta }}\gamma$ to specific values. We note that unlike the doubly averaged Hamiltonian (1), here the averaging process is carried out only over ${\lambda }_{9}$, under the restriction of the resonant relationship dictated by Equation (4). Moreover, the indirect part of the disturbing function must be retained in this case.

Of particular interest to the problem at hand are the resonant phase-space portraits of KBOs in apsidally anti-aligned configurations with respect to Planet Nine, with perihelion distances characteristic of typical scattered disk objects (i.e., q ∼ a₈). Suitably, adopting ${\rm{\Delta }}\gamma =\pi$ and a value for Γ that corresponds to $q=35\,\mathrm{au}$ at $a={({\ell }/k)}^{2/3}\,{a}_{9}$ for the frozen degree of freedom, we computed the averaged Hamiltonian (10) by adopting parameters relevant to the 1:1, 3:2, 2:1, and 5:2 MMRs. The corresponding ${\phi }_{\mathrm{res}}$–a diagrams are presented in Figure 5, on which we also plot the N-body trajectories of particles shown in Figure 4 in red. It is clear that irrespective of the specific resonance argument, the topology of ${{ \mathcal H }}_{\mathrm{res}}$ is keenly reminiscent of a pendulum-like structure that has been cut into two separate domains by vertical lines. These lines depict collision curves—i.e., values of ${\phi }_{\mathrm{res}}$ for which the averaging process fails due to the inherent singularity. Importantly, this means that a trajectory residing within one domain cannot migrate to the other domain without compromising the phase-protection mechanism of MMRs.

Zoom In Zoom Out Reset image size

Within each of the phase-space portraits shown in Figure 5, one of the two domains contains an $\infty$-shaped separatrix characterized by a hyperbolic equilibrium at its center, while the other domain possesses an elliptic fixed point at its core. Numerical integrations reported in the previous section showed that resonant trapping typically occurs in the domain that does not host the separatrix. That is, even though objects that exhibit a stable libration of ${\phi }_{\mathrm{res}}$ = 180 degrees away from those shown in Figure 4 do exist, they are very rare and we will not consider their evolution in detail. We further note that strictly speaking, the phase-space evolution depicted in Figure 5 implies that the libration of ${\phi }_{\mathrm{res}}$ does not represent a formal resonance because it is not enclosed by a separatrix (Delisle et al. 2012). This point, however, is of little practical consequence, since the phase protection facilitated by this pseudo-resonance patently represents a stabilizing mechanism for the simulated KBOs.

4.2. Secular Dynamics Inside MMRs

Having characterized the evolution of the fast degree of freedom in the preceding subsection, we now consider secular dynamics facilitated by resonant interactions. Generally speaking, in order to compute an $e\mbox{--}{\rm{\Delta }}\varpi$ phase-space diagram, we must specify the state of a–${\phi }_{\mathrm{res}}$ dynamics everywhere on the domain. To do so, we once again rely on the principle of adiabatic invariance.

Because the period of oscillation of ${\rm{\Delta }}\varpi$ greatly exceeds that of ${\phi }_{\mathrm{res}}$, we can define the adiabatic invariant (Neishtadt 1984)

$$\begin{eqnarray}&&{ \mathcal J }=\oint {\rm{\Lambda }}^{\prime} \,d{\phi }_{\mathrm{res}}^{{\prime} }=-\displaystyle \frac{1}{{\ell }}\oint {\rm{\Lambda }}\,d{\phi }_{\mathrm{res}}\end{eqnarray} \tag{ 11 }$$

associated with motion in the a–${\phi }_{\mathrm{res}}$ plane. Physically, this quasi-integral corresponds to the phase-space area occupied by the trajectory and is conserved to an excellent approximation, as long as the system does not encounter any criticality (i.e., hyperbolic fixed points, collision curves, etc.; Henrard 1993). As can be seen from Figure 4, bodies entrained in MMRs with P9 can have substantial libration amplitudes that are in essence determined by the state of the system at t = 0. For definitiveness, here we ignore this complication and instead assume that the libration amplitude is null, meaning ${ \mathcal J }=0$. From a computational point of view, the ${ \mathcal J }=0$ assumption is simplifying, since rather than finding the correct phase-space trajectory in the a–${\phi }_{\mathrm{res}}$ plane (specified by a given non-zero value of ${ \mathcal J }$) at every combination of e and ${\rm{\Delta }}\varpi$, we instead suppose that the system is adiabatically confined to the resonant fixed point and carry out the averaging process under the resonant equilibrium condition (Morbidelli & Moons 1993).

Because the resonant angle of interest ${\phi }_{\mathrm{res}}$ only contains the longitude of perihelion of Planet Nine and not that of the KBO (Equation (4)), the resonant equilibrium in the a–${\phi }_{\mathrm{res}}$ plane always resides on the ${a}_{\mathrm{eq}}={({\ell }/k)}^{2/3}\,{a}_{9}$ line. On the contrary, the equilibrium value of ${\phi }_{\mathrm{res}}$ itself shifts away from ${\phi }_{\mathrm{res}}=0\,\deg$ or ${\phi }_{\mathrm{res}}=180\,\deg$, in concert with oscillations of ${\rm{\Delta }}\varpi$. Correspondingly, in order to compute the equilibrium value of ${\phi }_{\mathrm{res}}$ as a function of e and ${\rm{\Delta }}\varpi$, we evaluated the function ${{ \mathcal H }}_{\mathrm{res}}$ (Equation (10)) along the $a={a}_{\mathrm{eq}}$ axis at every grid point on the ($e,{\rm{\Delta }}\varpi$) plane and found its local maximum.

Within certain regions of the $e\mbox{--}{\rm{\Delta }}\varpi$ diagram, the collision curves, shown as black vertical lines on Figure 5, can approach one another, shrinking the domain within which the resonant trajectory resides. This process is demonstrated in Figure 6, which shows a series of resonant a–${\phi }_{\mathrm{res}}$ diagrams corresponding to the 3:2 MMR at various values of ${\rm{\Delta }}\varpi$. Moreover, for specific values of e and ${\rm{\Delta }}\varpi$, the two collision curves can cross, engulfing the pseudo-resonant equilibrium point, which we assume the system occupies. In other words, the phase-space portrait in a–${\phi }_{\mathrm{res}}$ dictates a locus on the $e\mbox{--}{\rm{\Delta }}\varpi$ plane that cannot be crossed without compromising the phase-protection mechanism inherent to mean-motion commensurabilities.⁹ Consequently, even prior to computing the secular dynamics explicitly, we can identify a restricted domain on the $e\mbox{--}{\rm{\Delta }}\varpi$ plane, bounded by the collision locus, that can be stably explored by an orbit trapped at the center of an MMR.

Zoom In Zoom Out Reset image size

Generally speaking, the admissible domain shrinks as the integer (k − ℓ) increases. Accordingly, the full $e\mbox{--}{\rm{\Delta }}\varpi$ plane is stable for the 1:1 MMR, but the stability region is tightly confined about the $e\sim 0.85$, ${\rm{\Delta }}\varpi \sim \pi$ point for the 5:2 resonance. For clarity, here we only compute the secular diagram within this stability domain.

Expressing λ as a function of ${\lambda }_{9}$ through the pre-computed stationary value of ${\phi }_{\mathrm{res}}$, we calculate the averaged Hamiltonian,

$$\begin{eqnarray}{\overline{{ \mathcal H }}}_{\mathrm{rs}} & = & -\displaystyle \frac{1}{4}\displaystyle \frac{{ \mathcal G }\,{M}_{\odot }}{a}\displaystyle \frac{1}{{(1-{e}^{2})}^{3/2}}\displaystyle \sum _{j=5}^{8}\displaystyle \frac{{m}_{j}\,{a}_{j}^{2}}{{M}_{\odot }\,{a}^{2}}\\ & & +{\dot{\varpi }}_{9}\,\sqrt{{ \mathcal G }\,{M}_{\odot }a}(1-\sqrt{1-{e}^{2}})\\ & & -{\left.\displaystyle \frac{{ \mathcal G }{m}_{9}}{2\pi {\ell }}{\displaystyle \int }_{0}^{2\pi {\ell }}\left(\displaystyle \frac{1}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}| }-\displaystyle \frac{{\boldsymbol{r}}\cdot {{\boldsymbol{r}}}_{9}}{| {{\boldsymbol{r}}}_{9}{| }^{3}}\right)\right|}_{{\phi }_{\mathrm{res}}}d{\lambda }_{9},\end{eqnarray} \tag{ 12 }$$

on the admissible domain. Figure 7 shows the contours of the numerically averaged function (12) for the 1:1, 3:2, 2:1, and 5:2 MMRs (recall that the corresponding resonant phase-space portraits are shown in Figure 5). In addition to the semi-analytic level curves of ${\overline{{ \mathcal H }}}_{\mathrm{rs}}$, depicted with black lines, and the background color scale of the figure, we also overplotted the $e\mbox{--}{\rm{\Delta }}\varpi$ evolutions of the four particles emphasized in Figure 4 as orange curves. Clearly, the agreement between the semi-analytic results and N-body simulations is satisfactory, and any mild quantitative discrepancies can likely be attributed to our assumption of null resonant libration amplitude in the perturbative calculation.

Zoom In Zoom Out Reset image size

Taken together, these results yield the following insight into the dynamical evolution of distant KBOs induced by Planet Nine. First and foremost, long-term dynamical stability is facilitated by capture into MMRs. Objects that are not (fortuitously) scattered into mean-motion commensurabilities with Planet Nine initially are removed from the system by way of close encounters. Meanwhile, due to the specific ("corotation") nature of the resonant multiplets that guide the resonant motion, the evolution of distant KBO eccentricities and longitudes of perihelion are dominated by secular dynamics that ensure inside resonances. In turn, this separation between the degrees of freedom qualitatively explains why the purely secular phase-space portraits shown in Figure 2 approximately match the results of large-scale numerical simulations.

4.3. Critical Semimajor Axis

In light of the analysis presented above, it is evident that even though the purely secular treatment of dynamics outlined in Section 2 does not formally account for the full dynamical evolution observed within N-body calculations, it does provide a satisfactory approximation for $e\mbox{--}{\rm{\Delta }}\varpi$ evolution (Beust 2016). Accordingly, we now take advantage of this simplified framework to explore the dependence of the critical semimajor axis, ${a}_{\mathrm{crit}}$ (beyond which apsidal confinement ensues), on the parameters of Planet Nine.

As a proxy for ${a}_{\mathrm{crit}}$, we adopt the minimum value of a at which the anti-aligned secular equilibrium exists on the secular $e\mbox{--}{\rm{\Delta }}\varpi$ diagram.¹⁰ Practically, we calculate this quantity by computing ${\overline{\overline{{ \mathcal H }}}}_{{\rm{s}}}$ as a function of e along the ${\rm{\Delta }}\varpi =\pi$ line, gradually increasing a from 50 au, and noting the first instance where a local maximum appears between e = 0 and e = 1. Correspondingly, this calculation is carried out on an (${a}_{9},{e}_{9}$) grid for a given value of m₉.

We computed ${a}_{\mathrm{crit}}$ as a function of a₉ and e₉ for ${m}_{9}=5\,{m}_{\oplus }$, ${m}_{9}=10\,{m}_{\oplus }$, and ${m}_{9}=20\,{m}_{\oplus }$. Figure 8 depicts curves corresponding to ${a}_{\mathrm{crit}}=150\,\mathrm{au}$, ${a}_{\mathrm{crit}}=200\,\mathrm{au}$, and ${a}_{\mathrm{crit}}\,=250\,\mathrm{au}$ for these choices of m₉. As can be clearly seen in this figure, ${a}_{\mathrm{crit}}$ exhibits a rather mild dependence of m₉, and follows a shallow relationship between e₉ and a₉.

Zoom In Zoom Out Reset image size

It is important to note that observationally, the specific value of ${a}_{\mathrm{crit}}$ is not unequivocally determined. It is indeed possible to demonstrate that the clustering of orbits with $a\gtrsim 250\,\mathrm{au}$ in their respective longitudes of perihelion (ϖ) is statistically significant (Brown 2017). However, given that the confinement in the argument of perihelion (ω) may persist down to $a\sim 150\,\mathrm{au}$ (Trujillo & Sheppard 2014), it is possible that anti-aligned dynamics in fact ensues beyond $a\gt 150\,\mathrm{au}$, and the observed ϖ distribution in the $150\lt a\lt 250\,\mathrm{au}$ range is contaminated by metastable objects residing within the apsidally aligned regions of the $e\mbox{--}{\rm{\Delta }}\varpi$ diagram. Although perfectly plausible, this scenario is not as strongly supported in a raw statistical sense by the current data set (Brown 2017; Shankman et al. 2017), leaving some ambiguity as to where in the $150\lt a\lt 250\,\mathrm{au}$ interval the true value of ${a}_{\mathrm{crit}}$ resides.

Whatever the exact value of ${a}_{\mathrm{crit}}$ may be, it is worth noting that P9 parameters within N-body simulations that were found to match the observational data with a comparatively high probability by Brown & Batygin (2016) all lie between the ${a}_{\mathrm{crit}}=150\,\mathrm{au}$ and ${a}_{\mathrm{crit}}=250\,\mathrm{au}$ contours in Figure 8. In other words, these contours delineate a region of orbital element space that yields simulation results that compare well with the real solar system. At the same time, the semi-analytic calculations presented herein do not suffer from limitations in resolution (i.e., particle count) inherent to numerical experiments, and therefore likely inform a broader range of acceptable P9 parameters than was reported in Brown & Batygin (2016). Exploring this parameter range with an expanded suite of high-resolution $(N\gg 5000)$ N-body simulations constitutes an appealing avenue for future development of the Planet Nine hypothesis.

Up until this point, we have considered the orbital evolution induced upon the Kuiper Belt by Planet Nine under the strict assumption of coplanarity. As already discussed above, this assumption leads to somewhat idealized behavior and fails to capture three important aspects of the dynamics that emerge within the framework of full-fledged N-body simulations. First, rather than exhibiting stable resonance capture and remaining locked to a particular commensurability for the duration of the solar system's lifetime, real KBOs experience chaotic semimajor axis evolution and therefore explore a wide range of orbital period ratios with Planet Nine (this is evident, for example, in Figure 11 of Millholland & Laughlin 2017). Second, in addition to confinement in the longitude of perihelion, real long-period KBOs exhibit a clustering in the longitude of ascending node as well, which jointly leads a clustering in the argument of perihelion, $\omega =\varpi -{\rm{\Omega }}$. Third, the N-body simulations of Batygin & Brown (2016a) show that KBOs initially residing close to P9's orbital plane can occasionally undergo high-amplitude oscillations of the inclination, leading to the generation of retrograde KBOs (Batygin & Brown 2016b). Accordingly, the goal of this section is to characterize these features qualitatively.

5.1. The Anti-aligned Population

To carry out the exploration of non-planar dynamics within the framework of our simplified model, we repeated the numerical experiment with ${a}_{9}=700\,\mathrm{au}$ described in Section 3, allowing for a small inclination dispersion among the particles. In particular, initial particle inclinations were drawn from a half-normal distribution with a standard deviation of ${\sigma }_{i}\,=5\,\deg$, while longitudes of ascending node Ω were assumed to be uniformly distributed between 0 and $360\,\deg$. As before, the simulations are performed in a frame coplanar with the orbit of Planet Nine (meaning that ${i}_{9}=0$), and the Keplerian motion of the four inner giants is averaged out. Meanwhile, the ecliptic plane is envisioned to be inclined with respect to the orbit of P9 by a small amount (e.g., $10\mbox{--}20\,\deg$), although the corresponding reduction in the effective magnitude of the J₂ moment (Equation (3)) is ignored. Remarkably, this simple modification to the physical setup of the problem is sufficient to approximately recover the full breadth of the dynamical behavior observed within more detailed numerical experiments.

Figure 9 shows the orbital evolution of six (uniquely colored) representative objects, whose orbits remain roughly confined to the orbital plane of Planet Nine. The top panel depicts the longitude of perihelion (relative to the apsidal line of Planet Nine) as a function of semimajor axis. Clearly, allowing the KBOs to possess even a small inclination breaks the immutability of MMRs seen in the strictly coplanar model and renders the evolution of the semimajor axes chaotic. At the same time, the clustering of the orbits in the longitude of perihelion persists in spite of irregular semimajor axis evolution, and indeed, the secular portrait of orbital eccentricity depicted in the middle panel of the figure agrees well with the resonant-secular $e\mbox{--}{\rm{\Delta }}\varpi$ diagrams shown in Figure 7.

Zoom In Zoom Out Reset image size

We note that although the introduction of a finite inclination dispersion is sufficient to drive chaotic mixing in the semimajor axis, its efficiency is underestimated in our simplified model, since stochastic evolution is driven solely by P9 resonances. In reality, scattering facilitated by Neptune significantly enhances the rate of semimajor axis diffusion (particularly during the low perihelion phases of the secular cycle) and modifies the observed behavior on a quantitative level. As a consequence, the results obtained in our semi-averaged simulations depict a somewhat idealized realization of distant Kuiper Belt evolution.

The apsidally clustered resonance-hopping behavior observed in Figure 9 can be qualitatively understood within the framework of the semi-analytical theory described above. The introduction of a new ($i-{\rm{\Omega }}$) degree of freedom into the dynamics implies that rather than being forced by a single secular harmonic, the modulation of the a–${\phi }_{\mathrm{res}}$ phase-space portrait is now driven by two separate angles. This modulation transforms the collision loci depicted in Figure 7 into fuzzy chaotic bands because any given point ($e,{\rm{\Delta }}\varpi$) near the original collision locus may or may not lead to a collision, depending on the values of i and Ω. To this end, notice that in Figure 7, the curves showing the secular evolution are essentially tangent to the collision loci, which means that in a strictly planar physical setup, secular evolution almost never leads to close encounters between the particles and Planet Nine (i.e., most of the plotted trajectories do not intersect the collision loci). On the other hand, if the collision locus becomes a chaotic band, many more secular trajectories can infiltrate this irregular region, and with appropriate values of i and Ω, close encounters can ensue. When this happens, the particle enters a stochastic dynamical regime and hops in the semimajor axis until it locks into a new resonance, and the secular dynamics drives it away from the collision locus of the new resonance. Importantly, this type of behavior is observed in numerical integrations of real objects under the influence of Planet Nine (Becker et al. 2017; Millholland & Laughlin 2017).

When an object gets fortuitously trapped in some resonance (by entering it through the chaotic layer that surrounds the collision curves), its orbital evolution can be temporarily stabilized by the secular evolution in e and ${\rm{\Delta }}\varpi$. In other words, test particles have the tendency to exhibit prolonged periods of resonance locking (on timescales similar to the secular libration period) before breaking out and jumping to another commensurability. In fact, looking again at Figure 7 and imagining a chaotic band near each depicted collision locus, it is evident that the ($e,{\rm{\Delta }}\varpi$) evolution can drive a body away from the band at the peaks of the eccentricity cycle of the 3:2, 5:3, and 2:1 resonances (and also at the bottom of the eccentricity cycle in the 2:1 resonance). During this secular phase, close encounters between the particles and Planet Nine are no longer possible. However, as the eccentricity–perihelion cycle unfolds, the particle must eventually plunge back into the chaotic band. The object can thus experience new close encounters with Planet Nine and hop to another resonance where this process repeats. In reality, this sequence of events is further complicated by the fact that at the peak of the eccentricity cycle, objects are brought to an orbital state where q ∼ a₈ and suffer enhanced semimajor axis diffusion due to scattering off Neptune (or, equivalently, overlap with Neptune's exterior MMRs; Gomes et al. 2008). As a result, it is reasonable to assert that only dynamically "detached" objects that are not actively scattering off Neptune are presently entrained in MMRs with Planet Nine.

The dynamical evolution of the orbital inclination observed in the bottom panel of Figure 9 is a consequence of the chaotic rotation of the angular momentum vectors around Planet Nine's orbit normal. That is, their longitudes of ascending node relative to the orbital plane of Planet Nine are in circulation. However, if the inclination of P9 is sufficiently large relative to the ecliptic, then particles that are less inclined, with respect to Planet Nine's plane, than the inclination of P9 itself, will appear to have a librating node relative to the ecliptic. In other words, viewed from a coordinate system that coincides with the ecliptic plane, the orbits of these particles execute a libration around the forced $i-{\rm{\Delta }}{\rm{\Omega }}$ equilibrium forced by P9's inclination.

Taken together, our calculations suggest that the clustering of the longitudes of ascending node first noted in Batygin & Brown (2016a) is nothing but a trivial consequence of the bending of the Laplace plane away from the solar system's mean plane by Planet Nine. Moreover, the apparent libration of the ascending node Ω, together with the true libration of longitude of perihelion ϖ, produces the apparent libration of the argument of perihelion ω, as observed in the real data (Trujillo & Sheppard 2014). This implies that the orbital inclination of Planet Nine must simultaneously be sufficiently large for apparent nodal clustering to ensue (e.g., ${i}_{9}\gtrsim 10\mbox{--}20\,\deg$), but not be so large as to disrupt the stable confinement of the longitudes of perihelion (${i}_{9}\lesssim 40\,\deg ;$ Brown & Batygin 2016; Saillenfest et al. 2017).

5.2. The Highly Inclined Population

Perhaps the most remarkable consequence of P9-driven dynamics is exemplified by the induction of large-amplitude oscillations in the orbital inclinations of distant KBOs. Not only is this mode of orbital evolution a unique prediction of the Planet Nine hypothesis (Batygin & Brown 2016a), real objects presently entrained in this pattern of perturbation now comprise a firmly established part of the observational data set (Gomes et al. 2015). Accordingly, this regime of P9-induced dynamics constitutes one of the strongest lines of evidence for the existence of Planet Nine as no other dynamical model can reasonably account for the origin of the observed highly inclined trans-Neptunian object population. Let us now examine these extreme orbital excursions within the framework of our simplified numerical model.

The top panel of Figure 10 shows the inclination time series of six simulated particles with initial semimajor axes between 500 au and 600 au that do not remain bound to Planet Nine's orbital plane for the entire duration of the integration. As can be readily seen in this figure, during the latter half of the solar system's lifetime, each object abruptly enters a phase of extreme orbital variation, and upon experiencing a single large-scale oscillation of the inclination, rejoins the low-i population of apsidally anti-aligned bodies. Although representative, we note that orbital excursions of this sort are not always limited to a single cycle—some objects within the simulation suite experience a multitude of sequential oscillations. Moreover, the onset of high-i excursions is not limited to a small subset of particles—in our simulation, 38% of all stable objects experience at least one such excursion.

Zoom In Zoom Out Reset image size

An intriguing feature of the depicted evolution is that the pattern of $e\mbox{--}{\rm{\Delta }}\varpi$ dynamics changes drastically when a particle enters the highly inclined regime. As shown in the middle panel of Figure 10, rather than encircling an elliptic equilibrium located at ${\rm{\Delta }}\varpi =180\,\deg$ as in Figure 9, the $e\mbox{--}{\rm{\Delta }}\varpi$ projection of the phase-space portrait acquires a three-lobed shape with eccentricity maxima located $\sim \pm 70\,\deg$ away from perfect apsidal anti-alignment with Planet Nine. The corresponding instants where the eccentricities are maximized (and perihelion distance is minimized) are shown with points on the top panel of Figure 10 and lie almost exactly at $i\approx 90\,\deg$.

Taken together, the top and middle panels of Figure 10 show that the model predicts the highly inclined population to be most readily observable in a state that is approximately perpendicular to the ecliptic and is slightly sub-orthogonal in apsidal orientation with respect to the anti-aligned cluster of distant orbits. Note, however, that during phases of lower eccentricities (and higher perihelion distance), this population remains relatively well localized in the longitude of perihelion¹¹ around ${\rm{\Delta }}\varpi =180\,\deg$. This means that KBOs that lie beyond the current observational frontier exhibit an even more complex dynamical structure than those comprising the known long-period data set.

Although it is tempting to attribute the mode of orbital evolution shown in Figure 10 to the oft-cited Kozai–Lidov (KL) resonance (Kozai 1962; Lidov 1962), it is crucial to understand that the flavor of secular dynamics executed by the simulated particles is keenly distinct. Specifically, in contrast with conventional KL evolution (where orbital inclination is traded for eccentricity such that the particle's eccentricity is minimized when orbits become orthogonal), the simulated particles reach their peak eccentricities near $i\approx 90\,\deg$. Moreover, instead of being constrained by the conservation of the $\hat{z}\,$ component of the specific angular momentum vector $h=\sqrt{1-{e}^{2}}\cos (i)$ (Kinoshita & Nakai 1999), the orbital excursions observed in our simulation are accompanied by large-amplitude variations of this quantity (ranging from ∼1 to −1). This further implies that the KL resonance does not represent the primary driver of the depicted evolution.

Because the orbit of Planet Nine is assumed to have ${i}_{9}=0$ in our calculations, the only angles that appear in the secular Hamiltonian are ${\rm{\Delta }}\varpi$, ω, and their various linear combinations. We already argued above that the KL mechanism is not responsible for the observed evolution, so the libration of ω alone cannot facilitate the observed dynamics.¹² Simultaneously, ${\rm{\Delta }}\varpi$ cannot force oscillations in inclination, meaning that the secular resonance at play must contain both ${\rm{\Delta }}\varpi$ and ω. Correspondingly, we propose that the large-scale orbital variations depicted in Figure 10 are driven by the libration of the secular angle,

$$\begin{eqnarray}&&\theta ={\rm{\Delta }}\varpi -2\omega =2{\rm{\Omega }}-\varpi -{\varpi }_{9}.\end{eqnarray} \tag{ 13 }$$

Incidentally, this angle arises at the octupole order of expansion when the Hamiltonian is expressed as a series in semimajor axis ratios (Mardling 2010).

Adopting ${\rm{\Delta }}\gamma =-{\rm{\Delta }}\varpi$ and θ as the secular angles of the test particles, we identify the quantity

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{\sqrt{1-{e}^{2}}}{2}(1-\cos i)\end{eqnarray} \tag{ 14 }$$

as the action conjugate to θ. In the bottom panel of Figure 10, we show the evolution of canonical Cartesian coordinates related to the $(\theta ,{\rm{\Theta }})$ action-angle variables. From this figure, it is evident that during the high-inclination phase of orbital evolution, θ executes a bounded oscillation and the secular trajectory traces out the shape of a typical resonant separatrix (Henrard & Lamaitre 1983). Thus, the dynamics shown in Figure 10 is characterized by the simultaneous libration of the relative longitude of perihelion ${\rm{\Delta }}\varpi$, as well as the angle θ (which jointly leads to the libration of the longitude of ascending node) representing an exceptionally strong form of secular coupling.

Due to the lack of separation of timescales on which the two secular degrees of freedom operate, we cannot study the depicted large-amplitude variations of the inclination using the same flavor of semi-analytic perturbation theory as that outlined in Section 4. Nevertheless, we note that the onset of these large-scale oscillations almost always corresponds to the eccentricity minimum in the $e\mbox{--}{\rm{\Delta }}\varpi$ cycle that precedes the oscillation. Qualitatively, this implies that resonant-secular $e\mbox{--}{\rm{\Delta }}\varpi$ dynamics modulates the proximity parameter of the secular $\theta -{\rm{\Theta }}$ resonance, such that when e approaches a critical value, a dynamical gateway toward the capture of low-inclination objects into the secular resonance characterized by the libration of θ temporarily opens. Accordingly, the coercion of low-i objects onto trajectories that experience large-amplitude oscillations of inclination is a fundamentally stochastic process. In turn, this means that some of the currently observed members of the apsidally clustered population can in principle join the highly inclined population in the future, and vice versa.

Given that our spatial model reproduces the key features of more detailed numerical simulations, it warrants a rudimentary comparison with the observational data in light of the semi-analytical insight into the governing dynamics developed above. Because we are resolving neither the Keplerian motion of the canonical giant planets nor their inclination with respect to Planet Nine, here we will not consider the clustering of the longitudes of ascending node at low i. Instead, we will focus exclusively on the confinement of the longitude of perihelion as well as the behavior of the highly inclined long-period objects.

As shown in Figure 1, there are currently 10 known objects with $a\gt 250\,\mathrm{au}$ (shown in purple) that comprise the primary ϖ cluster. The observational data set also shows the existence of two objects that are diametrically opposed to the mean orientation of this cluster (shown in green), as well as a single outlier, 2015 GT₅₀ (shown in gray), that does not fall within either the apsidally aligned or anti-aligned subpopulations of objects. In order to meaningfully compare the expectations of the model with the data, we extended our integrations such that the initial semimajor axis distribution of the particle disk stretches out to 850 au.

The orbital footprint of all simulated long-term stable and metastable particles is shown in Figure 11, where the top panel depicts the longitude of perihelion as a function of the semimajor axis (as in the top panel of Figure 9). Meanwhile, the middle panel shows the orbital inclination as a function of the argument of perihelion and the bottom panel elucidates the action Θ as a function of its conjugate angle, θ. As a crude proxy for observability, we adopt simple cuts of the numerical output at $q\leqslant 100\,\mathrm{au}$ and $i\leqslant 40\,\deg$. Points corresponding to bodies with dynamical lifetimes in excess of 4 Gyr that do not simultaneously satisfy these criteria are shown in gray, while those that do are shown in blue or red, depending on their inclination evolution. In particular, objects that remain confined to the plane of Planet Nine throughout the integration are shown in blue. On the other hand, bodies that experience large-amplitude inclination oscillations at any point in their evolution are depicted in red.

Zoom In Zoom Out Reset image size

In agreement with the results of Batygin & Brown (2016a), the top panel of Figure 11 shows the emergence of a well-defined cluster of apsidally anti-aligned orbits that are contaminated by trajectories that circulate in the longitude of perihelion. One striking example of such a circulating trajectory is shown as a vertical blue line with a semimajor axis of $a\approx 335\,\mathrm{au}$. As discussed in Sections 3 and 4, the dynamics of apsidally confined trajectories are driven by the resonant harmonic ${\phi }_{\mathrm{res}}$ (see Equation (4)), while the comparatively less frequent apsidally circulating trajectories tend to reside within resonant multiplets characterized by the low-amplitude libration of other resonant angles that contain the particle's longitude of perihelion, ϖ. The aforementioned observational data are overplotted on this panel and are color-coded in the same way as in Figure 1.

As expected, all observed KBOs that comprise the primary cluster (shown in purple) are seamlessly explained by the simulation results. More remarkably, however, the outlier within the data (shown in gray) is also naturally reproduced by the model as an object that belongs to the class of observable stable particles that exhibit apsidal circulation, i.e., those entrained in resonances characterized by angles that contain ϖ (see Appendix B for an analysis). In the present example, this agreement stems from the fact the observed object is very close to 3:1 resonance with Planet Nine. While this correspondence may be purely accidental, it emphasizes that the mere existence of a small number of apsidally unconfined objects that do follow the overall pattern exhibited by the data does not constitute strong evidence against the Planet Nine hypothesis.

Although a similar narrative could in principle be invoked for the two objects that are the apsidally aligned with Planet Nine (shown in green), it can be more reasonably speculated that these bodies are in fact subject to purely secular interactions with P9. Recall from Section 2 (Figure 2) that apsidally aligned objects residing in the secular domain are protected from close encounters by the geometric collinearity of the orbits. Correspondingly, objects that never attain low perihelion distances at the top of their eccentricity cycle and are therefore close to the secular equilibrium point exhibit long-term stable apsidal libration about ${\rm{\Delta }}\varpi =0$. Conversely apsidally aligned objects with low perihelion distances are metastable, as they only begin to scatter off Planet Nine once they precess onto their respective tangential collision curves, yielding dynamical lifetimes that are on the order of (a fraction of) the precession timescale—i.e., ∼few × 100 Myr.

Because in our numerical experiments we restricted the initial perihelion range of all particles to $q\in (30,36)$ au, our simulations do not produce any long-term stable apsidally aligned orbits. Instead, secular dynamics observed within our N-body simulations are dominated by the metastable particles that originate at high eccentricity and eventually precess toward the tangential configuration, where they are scattered out. To demonstrate the apsidal behavior of this subpopulation of orbits, objects with dynamical lifetimes between 100 and 500 Myr are shown as orange dots in Figure 11. These particles clearly cluster around ${\rm{\Delta }}\varpi =0$ and provide an excellent match to the aligned (green) data points shown in the figure.

Naturally, this explanation would not be sensible if the entire distant Kuiper Belt had been generated ∼4 Gyr ago and never replenished since then. This, however. is not the case in the real solar system: just as the highly inclined long-period KBOs are routinely scattered inwards to create the population of retrograde bodies with $a\lt 100\,\mathrm{au}$ (Batygin & Brown 2016b), scattered disk objects with $a\lt 250\,\mathrm{au}$ are continuously scattered outward by Neptune, resupplying the distant trans-Neptunian region with metastable KBOs (Gomes et al. 2008). Thus, our model points to the possibility that 2013 FT₂₈ and 2015 KG₁₆₃ are relative newcomers to the distant Kuiper Belt and will eventually be destabilized by short-periodic interactions with Planet Nine.

There exists yet another interpretation of the apsidally aligned data points as well, which is not well represented by Figure 11 due to our choice of low-q initial conditions. Particularly, rather than belonging to the aforementioned metastable subpopulation of bodies that experience pure secular evolution, 2013 FT₂₈ and 2015 KG₁₆₃ could be entrained in the stable secular libration island around ${\rm{\Delta }}\varpi =0$ and are presently observed near the peak of their respective eccentricity cycles. In order to unequivocally distinguish between these two interpretations, we would need to know the exact orbital parameters of Planet Nine. However, we simultaneously note that within the context of the long-term stable interpretation, some additional mechanism other than scattering off Neptune (such as, say, interactions with the birth cluster; Morbidelli & Levison 2004; Adams 2010) would likely be required to initially raise the perihelion distances of these objects and lock them into the apsidally aligned secular libration island. This is because bodies scattered to distant elliptic orbits by the traditional giant planets would necessarily have low perihelia (like the initial conditions of our simulation) and therefore could not reside within the stable libration island.

The yellow data points shown in the middle and bottom panels of Figure 11 represent the population of distant ($a\gt 250\,\mathrm{au}$), highly inclined ($i\gt 40\,\deg$) objects with $q\lt 30\,\mathrm{au}$ discussed in Section 5.2. Although these objects conform to the dynamical streamlines traced out by the simulated particles relatively well, it is important to keep in mind that (by virtue of having $q\lt 30\,\mathrm{au}$) these bodies have eccentricities that are much closer to unity than any of the particles in our simulations. This means that the observed objects are drawn from the extreme end of the broader high-i population and have likely had their orbits somewhat perturbed by the canonical giant planets. As a result, we expect that the agreement between theory and observations will be even better for (yet undiscovered) highly inclined long-period KBOs with $q\gt 30\,\mathrm{au}$. Certainly, continued observational monitoring of the distant Kuiper Belt outside of the ecliptic plane constitutes a viable avenue toward further characterization of the long-term dynamical evolution induced by Planet Nine.

As a final point, it is instructive to remark on the extent to which the calculations described above are in agreement with more detailed simulations that fully resolve the orbital motion of the inner giants. In particular, the top panel of Figure 11 can be readily compared with Figures 5 and 8 of Batygin & Brown (2016a), which depict simulations with initial conditions very similar to those considered herein. As expected, upon comparison of these numerical experiments, we find that the general survival rate of particles is somewhat lower in simulations that model Neptune directly. Specifically, in our semi-averaged calculations, 19% of particles initialized between 250 au and 550 au remain stable over 4 Gyr, while only 7% of objects in this initial semimajor axis and perihelion range survive the full integration span in a simulation where Neptune is modeled directly. This number increases slightly to 10% when Planet Nine is endowed with a ${i}_{9}=30\,\deg$ inclination with respect to the inner solar system.

Another notable difference between the direct and semi-averaged calculations lies in the fact that simulated objects that circulate in perihelion and match our observability criteria are somewhat less prevalent in the simulations that include Neptune's short-periodic perturbations. This is likely because apsidally circulating particles tend to experience diminished eccentricity variations (see Appendix B for details) and thus retain low perihelion distances, where they are more likely to be removed by Neptune. Simultaneously, we note that for the exact same reason, such objects are more readily discoverable by astronomical surveys and are thus bound to be overrepresented within the observational census of long-period KBOs. Accordingly, further characterization of the P9-sculpted orbital distribution, fully accounting for the overlying observational biases using high-resolution numerical experiments, constitutes an important step toward the continued evaluation of the Planet Nine hypothesis within the framework of the emergent data set.

Within the current observational census of trans-Neptunian objects, the longest-period orbits exhibit an unexpected collective structure that is most readily attributed to gravitational perturbations induced by a yet unseen, massive planet. Although numerical simulations aimed at reproducing the Kuiper Belt's orbital makeup through gravitational interactions with Planet Nine are now plentiful in the literature (Batygin & Brown 2016a; Brown & Batygin 2016; Becker et al. 2017; Millholland & Laughlin 2017), the physics of the dynamical processes responsible for shaping the distant Kuiper Belt remains largely unclear (Beust 2016). In this work, we sought to resolve this problem and characterize the dynamical evolution induced by Planet Nine in long-period KBOs on semi-analytic grounds.

The specific aim of this work has been to qualitatively understand the three primary lines of evidence for the existence of Planet Nine. They are the (i) orbital clustering of long-period KBOs, (ii) dynamical detachment of KBO orbits from Neptune, and (iii) generation of highly inclined/retrograde bodies within the solar system. We note that none of these effects are new and have already been pointed out in the work of Batygin & Brown (2016a). Accordingly, the primary purpose of this study has been to create an analytical guide for the interpretation of existing and future numerical results, rather than to generate new ones. In doing so, we were guided by a series of queries that we outlined in the introduction. Let us now recall and provide the answers to these questions.

• 1.  
  What role (if any) do resonant interactions play within the dynamical evolution induced by Planet Nine? If resonances are prevalent, what order/multiplet harmonics dominate the dynamics, and what are their characteristic widths?

In the most idealized case of a strictly planar physical setup, all high-eccentricity particles that occupy stable orbits within Planet Nine's gravitational domain of influence derive their prolonged dynamical lifetimes from the phase-protection mechanism inherent to MMRs.¹³ The critical angles associated with these resonances typically have the form ${\phi }_{\mathrm{res}}\,=k\,{\lambda }_{9}-{\ell }\,\lambda -(k-{\ell }){\varpi }_{9}$, and due to the absence of ϖ (i.e., the longitude of perihelion of the KBO) from the expression, these harmonics only modulate the mean anomalies and semimajor axes of the KBOs.¹⁴ The most prevalent interior resonances that arise in our calculations correspond to the 1:1, 3:2, 2:1, and 5:2 period ratios, with secondary contributions from the 5:4, 4:3, 5:3, and 3:1 commensurabilities. Provided nominal Planet Nine parameters, the widths of these resonances lie on the order of $\delta a\sim 2\mbox{--}15\,\mathrm{au}$, and scale as $\delta a\propto {m}_{9}^{1/2}$ (Henrard & Lamaitre 1983).

If the strict coplanarity restriction is lifted and the particles are endowed with a small inclination dispersion relative to the plane of Planet Nine's orbit (as in the real Kuiper Belt), a large portion of resonant dynamics become chaotic and facilitate an essentially diffusive semimajor axis evolution, with particles hopping from one resonance to another. Moreover, for high-eccentricity orbits that reach $q\lesssim 36\,\mathrm{au}$, stochastic semimajor axis transport is further enhanced by gravitational scattering off Neptune (Gomes et al. 2008). These complications imply that only dynamically detached objects with $q\gtrsim 40\,\mathrm{au}$ can reasonably be speculated to currently reside in MMRs with Planet Nine, and any attempt to calculate the present-day semimajor axis of Planet Nine exclusively from resonant relationships with the observed KBOs (Malhotra et al. 2016) may be spoiled by the chaotic nature of the underlying dynamics.

• 2.  
  What role (if any) do secular interactions play within the dynamical evolution induced by Planet Nine? If dominant, how are close encounters avoided on coplanar, anti-aligned orbits? Moreover, if resonant interactions are relevant to the Planet Nine hypothesis, why does the purely secular phase-space portrait provide a good match to the results of numerical simulations?

While MMRs stabilize the orbits against close encounters, they do little to modulate the eccentricities and longitudes of perihelia of the affected bodies. As a result, the $e\mbox{--}{\rm{\Delta }}\varpi$ dynamics induced in distant KBOs by Planet Nine is largely secular. This explains why the doubly averaged treatment of the dynamics outlined in Section 2 provides a good approximation to the results of the numerical simulations discussed in Section 3 (Beust 2016). At the same time, it is important to keep in mind that the true resonant-secular evolution facilitated by Planet Nine's mean gravitational potential is subtly different from the purely secular limit, since the orbital averaging process itself is subject to the resonant relationships among the orbital phases of the bodies. The delicate differences between purely secular and resonant-secular dynamics induced by Planet Nine can be noted by comparing Figures 2 and 7.

The eccentricity–perihelion projection of the resonant-secular phase-space portraits outlined in Figure 7 further shows that the apsidal clustering of distant KBOs and the detachment of perihelion from Neptune's orbit exemplified by objects such as Sedna and $2012\,$ VP₁₁₃ (Brown et al. 2004; Trujillo & Sheppard 2014) are actually the same physical effect. That is, as the orbits of KBOs evolve along the level curves of the resonant-secular Hamiltonian, they are forced to encircle a stable equilibrium that resides at ${\rm{\Delta }}\varpi =180\,\deg$. Thus, the libration of a KBO's longitude of perihelion around the apsidally anti-aligned configuration with respect to Planet Nine is accompanied by conjugate oscillations of the eccentricity that periodically detach (and reattach) the perihelion from (to) Neptune.

Although the gravitational influence of P9 alone provides a perfectly adequate mechanism for perihelion detachment of long-period objects, we cannot exclude the possibility that KBOs were additionally affected by other dynamical processes (such as interactions with the birth cluster; Morbidelli & Levison 2004) during the solar system's infancy. If so, some fraction of the distant clustered population may occupy secular cycles that will never bring their perihelia sufficiently close to the orbit of Neptune for scattering to ensue. In either case, our calculations reveal that the maximum width of the perihelion cluster is limited to $| {\rm{\Delta }}\varpi -180\,\deg | \lesssim 90\,\deg$—a restriction facilitated by both the character of the dynamics itself as well as the collision locus that limits the domain where close encounters can be avoided on the $e\mbox{--}{\rm{\Delta }}\varpi$ plane. Within this framework, Planet Nine's mass merely regulates the size of the chaotic layer and the timescale on which the perihelion cluster get sculpted.

• 3.  
  What parameters determine the critical semimajor axis corresponding to the transition between the randomized and clustered longitudes of perihelion? What physical effect controls this transition?

Given that the long-term dynamics induced in KBOs by Planet Nine are essentially secular in nature, the critical semimajor axis beyond which orbital clustering ensues corresponds to a point where quadrupolar torques induced by Planet Nine begin to dominate over those arising from the canonical giant planets. Computed in this framework, curves corresponding to critical semimajor axes of ${a}_{\mathrm{crit}}=150,200$ and 250 au are delineated in ${e}_{9}-{a}_{9}$ space for ${m}_{9}=5,10,$ and $20\,{m}_{\oplus }$ in Figure 8. While these loci provide an approximate measure of ${a}_{\mathrm{crit}}$, we note that the transition between apsidally confined and randomized orbits in our N-body simulations is somewhat gradual (as shown in Figure 3).

With decreasing semimajor axis, the relative number of resonant bodies that are apsidally clustered decreases with respect to those that are not.¹⁵ This means that although ${a}_{\mathrm{crit}}$ provides a characteristic semimajor axis that corresponds to the onset of orbital clustering, the real Kuiper Belt at $a\lesssim {a}_{\mathrm{crit}}$ will show orbital clustering that is increasingly contaminated by non-anti-aligned bodies—an effect seen in the real data (Trujillo & Sheppard 2014; Shankman et al. 2017). Meanwhile, our calculations suggest that the clustering of the orbital planes of distant KBOs and the corresponding nodal confinement is a simple consequence of the tilting of the Laplace plane away from the ecliptic by Planet Nine.

• 4.  
  What is the qualitative behavior of the inclination dynamics within the framework of P9-driven evolution? What dynamical process allows some of the objects to acquire exceptionally high inclinations in the distant Kuiper Belt?

Not all objects affected by Planet Nine remain confined to the Laplace plane on multi-Gyr timescales. Instead, a subset of long-period KBOs execute large-scale oscillations in the orbital inclination as well as eccentricity. This mode of P9-induced evolution is distinct from the Kozai–Lidov mechanism and is driven by a high-order secular resonance that is characterized by the simultaneous libration of the critical angle $\theta ={\rm{\Delta }}\varpi -2\omega$ as well as ${\rm{\Delta }}\varpi$. This doubly resonant form of secular coupling forces a particular strenuous exchange of angular momentum and generally leads to the acute orbit-flipping behavior of distant KBOs (Li et al. 2014).

The onset of these large-amplitude orbital excursion is fundamentally chaotic and is facilitated by variations in the angular momentum deficit that is driven by resonant-secular $e\mbox{--}{\rm{\Delta }}\varpi$ oscillations. That is, oscillations in the eccentricity modulate the system's proximity to the $\theta$ resonance, periodically allowing low-inclination orbits to enter the highly inclined dynamical regime. As the subsequent orbital evolution unfolds, the eccentricity reaches a peak of its cycle when the inclination is approximately $i\approx 90\,\deg$, meaning that bodies belonging to the highly inclined population are most readily observable in a state roughly perpendicular to the plane of the solar system. This qualitatively explains the observed dynamical state of large semimajor axis centaurs, which constitutes the third major line of evidence for the existence of Planet Nine.

There exist two other secondary lines of evidence for the existence of Planet Nine, which we did not discuss in this paper. The first is the obliquity of the Sun (Bailey et al. 2016; Lai 2016; Gomes et al. 2017). One reason we chose to not discuss this dynamical effect is because it is intrinsically trivial. It is well-known that secular coupling between two mutually inclined orbits results in a regression of the node of both orbits, meaning that the twisting of the giant planets' orbital plane out of alignment with the solar spin axis is an inescapable consequence of Planet Nine's existence. As a result, the genuinely remarkable aspect of this calculation is not that a spin–orbit misalignment can be excited, but the fact that a Planet Nine configuration close to the one deduced from distant Kuiper Belt constraints can adequately reproduce the solar obliquity when its gravitational influence is exerted over the entire lifetime of the solar system.

At the same time, we note that strictly speaking, Planet Nine is not required to explain the obliquity of the Sun, and other theoretical models exist. For example, Batygin (2012) argued that a primordial binary companion to the solar system could have excited the observed spin–orbit misalignment through the same exact mechanism, while Lai et al. (2011) proposed that magnetic interactions between the Sun and the inner regions of the protosolar nebula could have accomplished the same task. We note, however, that while a multitude of processes could have contributed to the observed obliquity of the Sun, the gravitational influence of Planet Nine provides the only dynamical mechanism that is directly testable.

Another population of objects that Planet Nine naturally generates is the highly inclined component of the proximate ($a\lt 100\,\mathrm{au}$) Kuiper Belt. This subset of bodies includes all objects with $i\gtrsim 35\,\deg$ that do not naturally emerge from simulations of the solar system's primordial evolution (Levison et al. 2008; Batygin et al. 2011; Nesvorný 2015) and include the retrograde objects "Drac" (Gladman et al. 2009) and "Niku" (Chen et al. 2016). Although these objects appear observationally distinct from the nearly perpendicular high semimajor axis centaurs discussed above, full-fledged N-body simulations reported in Batygin & Brown (2016b) show that the highly inclined $a\lt 100\,\mathrm{au}$ objects are simply the large-a objects that have been scattered inwards by Neptune. As a consequence, they do not require a separate dynamical explanation from the objects undergoing the large-scale orbital excursions discussed in Section 5.2 of this paper.

Cumulatively, the aforementioned lines of evidence constitute a compelling case for the existence of Planet Nine, as none of the objects within the current observational data set exert any significant tension upon the model. In light of this broad-ranging agreement (short of a direct detection of Planet Nine), the most direct avenue toward further reinforcement or falsification of our theory is the continued detection of $a\gtrsim 250\,\mathrm{au}$ KBOs with the aim to better establish the statistical significance of the clustering of longitudes of perihelion and ascending node (Brown 2017; Shankman et al. 2017). To this end, we reiterate that even though the contamination of the clustered orbital pattern by unconfined particles is an expected result of P9's gravitational influence, a notable grouping of long-period trajectories in physical space remains a key feature of the dynamical model.

Although the evidence for P9 remains strong, as already discussed in the introduction, a simple proposition of an extant planet beyond Neptune does not amount to a meaningful theoretical prediction. Instead, the Planet Nine hypothesis is uniquely defined by the combination of observational signatures it explains and the specific dynamical mechanisms through which these astronomical patterns arise. Thus, the final evaluation of the Planet Nine hypothesis will not simply correspond to a detection of a planet beyond Neptune, but the confrontation of the outlined theory with the dynamical evolution induced by this planet. Thankfully, the observational prospects for the direct detection of Planet Nine either through ongoing or future surveys are quite promising (Brown & Batygin 2016; Fortney et al. 2016; Holman & Payne 2016a, 2016b; Millholland & Laughlin 2017), and it is likely that the concluding assessment of the theoretical model outlined in this paper will occur on a timescale considerably shorter than a decade.

We are thankful to Mike Brown, Greg Laughlin, Chris Spalding, Matt Holman, Gongjie Li, Tali Khain, and Elizabeth Bailey for illuminating discussions. Additionally, we are grateful to Sarah Millholland for providing a thorough and insightful referee report, as well as to Charles Fairchild, whose generous support facilitated this collaboration.

Throughout this paper, we relied on various closed-form computations of the averaged Hamiltonian (specifically, expressions (1), (10) and (12)), without explicitly stating how the calculations were carried out. Let us now comment on the practical details inherent to these evaluations. In the plane, the Cartesian coordinates of the particle's position vector, ${\boldsymbol{r}}$, are given by the well-known relations (Morbidelli 2002)

$$\begin{eqnarray}x & = & a(\cos E-e)\cos \varpi -a\sqrt{1-{e}^{2}}\sin E\sin \varpi \\ y & = & a(\cos E-e)\sin \varpi -a\sqrt{1-{e}^{2}}\sin E\cos \varpi ,\end{eqnarray} \tag{ 15 }$$

where E is the eccentric anomaly. Identical expressions (with subscript 9) apply to the position vector of Planet Nine.

While the aforementioned Cartesian coordinates are most naturally expressed in terms of the eccentric anomaly, canonical averaging of the Hamiltonian is carried out with respect to the mean longitude, λ. The two quantities are related via Kepler's equation. Therefore, rather than solving Kepler's equation at every computational step to express x and y in terms of λ, it is more convenient to integrate directly with respect to E by introducing the Jacobian (Pichierri et al. 2017)

$$\begin{eqnarray}&&d\lambda =(1-e\,\cos E){dE}.\end{eqnarray} \tag{ 16 }$$

For the case of the doubly averaged secular Hamiltonian (1), it is appropriate to carry out the integral assuming that λ and ${\lambda }_{9}$ are not correlated. This is, however, not the case for the singly averaged resonant and resonant-secular Hamiltonians (10) and (12). Particularly, in these instances, the mean longitudes of the particle and Planet Nine are linked to one another through the resonant relationship

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{\phi }_{\mathrm{res}}-(k-{\ell }){\varpi }_{9}+k\,{\lambda }_{9}}{{\ell }}.\end{eqnarray} \tag{ 17 }$$

Accordingly, in this case, Kepler's equation must be solved at every iteration to obtain the eccentric anomaly of the particle as a function of ${\lambda }_{9}$.

In order to construct resonant-secular phase-space portraits in the ${ \mathcal J }\to 0$ limit, it is necessary to first map out the equilibrium value of the resonant angle ${\phi }_{\mathrm{res}}$ on the $(e,{\rm{\Delta }}\varpi )$ domain. For the specific problem at hand, this can be done by sampling the resonant Hamiltonian (10) as a function of ${\phi }_{\mathrm{res}}$ along the $a={({\ell }/k)}^{2/3}\,{a}_{9}$ line and numerically finding the relevant local maximum. Incidentally, the equilibrium value of ${\phi }_{\mathrm{res}}$ typically lies halfway between the collision curves, which are easily obtained by solving the simultaneous equations $x={x}_{9},y={y}_{9}$, and substituting the resulting values of E and E₉ into Kepler's equation to yield the values of ${\phi }_{\mathrm{res}}$ that correspond to collisional trajectories.

In Section 4 of the main text, we constructed a semi-analytic model for apsidally confined resonant orbits residing in interior 5:2, 2:1, 3:2, and 1:1 MMRs with P9. As demonstrated in Figure 3, however, particles residing at somewhat lower values of semimajor axes (particularly those entrained in the 3:1 and 4:1 commensurabilities) predominantly exhibit apsidal circulation. Let us briefly consider the dynamical behavior of such trajectories in greater detail.

Representative trajectories of particles locked in 3:1 and 4:1 resonances, drawn from the $a=700\,\mathrm{au}$ simulation described in Section 3, are shown in Figure 12. The panels in this figure are analogous to those depicted in Figure 4 of the main text, with the exception that the librating resonant angle that drives the dynamics has the form

$$\begin{eqnarray}&&{\varphi }_{\mathrm{res}}=k\,{\lambda }_{9}-{\ell }\,\lambda -\varpi -(k-{\ell }-1){\varpi }_{9}={\phi }_{\mathrm{res}}-{\rm{\Delta }}\varpi .\end{eqnarray} \tag{ 18 }$$

Because these orbits are characterized by the circulation of ${\rm{\Delta }}\varpi$, no resonant multiplets other than ${\varphi }_{\mathrm{res}}$ are in libration. Note further that eccentricity modulation associated with the evolution of the apsidal angle ${\rm{\Delta }}\varpi$ is in this case very mild, especially compared to that shown in Figure 4.

Zoom In Zoom Out Reset image size

In order to more closely examine the phase-space evolution of the apsidally circulating orbits, we follow the same semi-analytic procedure as that outlined in Section 4. To avoid redundancy, we focus exclusively on the 3:1 MMR, since the 4:1 MMR exhibits a very similar behavior. Similarly to Equation (7), we define the canonical coordinates that identify ${\varphi }_{\mathrm{res}}$ as a reference angle:

$$\begin{eqnarray}{\rm{\Lambda }}^{\prime\prime} ={\rm{\Lambda }}/{\ell } & {\varphi }_{\mathrm{res}}^{{\prime\prime} }=-{\varphi }_{\mathrm{res}}\\ {{\rm{\Lambda }}}_{9}^{{\prime\prime} }={{\rm{\Lambda }}}_{9}+(k/{\ell }){\rm{\Lambda }} & {\lambda }_{9}^{{\prime\prime} }={\lambda }_{9}\\ {\rm{\Gamma }}^{\prime\prime} ={\rm{\Gamma }}+{\rm{\Lambda }}/{\ell } & {\rm{\Delta }}\gamma =\gamma -{\gamma }_{9}=-{\rm{\Delta }}\varpi \\ {{\rm{\Gamma }}}_{9}^{{\prime\prime} }={{\rm{\Gamma }}}_{9}+{\rm{\Gamma }}+((k-{\ell })/{\ell }){\rm{\Lambda }} & {\gamma }_{9}^{{\prime\prime} }={\gamma }_{9}.\end{eqnarray} \tag{ 19 }$$

As before, upon averaging the Hamiltonian with respect to the fast angle ${\lambda }_{9}^{{\prime\prime} }$, we are left behind with an adiabatic system that is characterized by a resonant degree of freedom in (${\rm{\Lambda }}^{\prime\prime} -{\varphi }_{\mathrm{res}}^{{\prime\prime} }$) as well as a secular degree of freedom in (${\rm{\Gamma }}^{\prime\prime} -{\rm{\Delta }}\gamma$).

Freezing ${\rm{\Gamma }}^{\prime\prime}$ at a value that corresponds to $q=35\,\mathrm{au}$ at $a={({\ell }/k)}^{2/3}\,{a}_{9}$ and setting ${\rm{\Delta }}\varpi =-{\rm{\Delta }}\gamma =\pi$, we computed the resonant phase-space diagram akin to those presented in Figure 5 for the interior 3:1 resonance. However, in this case, the Hamiltonian (10) was computed under the constraint of the resonant relationship (18) instead of Equation (4). The resulting (a−${\varphi }_{\mathrm{res}}$) diagram is shown in the left panel of Figure 13 with the result of our N-body simulation overplotted in red.

Zoom In Zoom Out Reset image size

Unlike the long-term stable orbits that exhibit steady perihelion clustering, apsidally circulating trajectories shown in Figure 12 reside in the subdomain of the (a−${\varphi }_{\mathrm{res}}$) diagram that is occupied by the $\infty$-shaped separatrix and encircle this homoclinic curve from the outside. In contrast with the neighborhood of the elliptic equilibrium point discussed in Section 4.1, this phase-space domain is never fully swept by the collision curves as ${\rm{\Delta }}\varpi$ swings from 0 to $2\pi$. As a result, trajectories that surround this resonance can safely rotate through all possible values of ${\rm{\Delta }}\varpi$ without compromising the resonant phase-protection mechanism.

This means that in the case of resonances characterized by the libration of ${\varphi }_{\mathrm{res}}$, there is no equivalent of the $(e\mbox{--}{\rm{\Delta }}\varpi )$ collision locus that arises within the framework of MMRs characterized by the libration of ${\phi }_{\mathrm{res}}$ (shown as the bounding curves in Figure 7). Nevertheless, the resonance that resides at the core of the apsidally circulating trajectory (as shown in Figure 13) undergoes topological changes driven by the evolution of ${\rm{\Delta }}\varpi$. In particular, as ${\rm{\Delta }}\varpi$ shifts away from π toward 0, the $\infty$-shaped separatrix deforms asymmetrically and eventually vanishes, as two of the three associated fixed points disappear. Finally, as ${\rm{\Delta }}\varpi$ tends closer to 0, the remaining (elliptic) fixed point returns to the origin, such that the trajectory once again encircles ${\varphi }_{\mathrm{res}}=0$ at ${\rm{\Delta }}\varpi =0$.¹⁶ In other words, as ${\rm{\Delta }}\varpi$ circulates from 0 to $2\pi$, the resonant center of ${\varphi }_{\mathrm{res}}$ experiences a concurrent oscillation about ${\varphi }_{\mathrm{res}}=0$ with an amplitude of order $\sim \pi /2$ (as shown in the right panel of Figure 12).

In order to approximately elucidate the resonant-secular dynamics that unfolds in the ${\rm{\Delta }}\varpi$ circulating regime, we computed the averaged Hamiltonian (12) under the constraint of the resonant relationship (18). Following the discussion outlined in Section 4.2, we assume that the adiabatic invariant ${ \mathcal J }=0$, and thereby confine ourselves to the a–${\varphi }_{\mathrm{res}}$ equilibrium point that remains extant (and stable) in the ${\rm{\Delta }}\varpi \in (0,\pi )$ range, and simply reflect the computed portrait onto the ${\rm{\Delta }}\varpi \in (\pi ,2\pi )$ domain. The corresponding $(e\mbox{--}{\rm{\Delta }}\varpi )$ diagram is shown in the right panel of Figure 13, with the N-body trajectory overplotted in orange.

Although the semi-analytical secular portrait matches the numerically computed evolution well, we note that for these apsidally circulating trajectories, the ${ \mathcal J }=0$ assumption is a relatively crude one, since the orbit itself resides outside the separatrix and does not directly encircle the equilibrium point. To this end, the kinks in the semi-analytical curves shown in Figure 13 at ${\rm{\Delta }}\varpi =\pi$ are an unphysical consequence of this assumption and would disappear if the portrait were more carefully computed by carrying out the averaging process along the contour of an (a−${\varphi }_{\mathrm{res}}$) trajectory that encircles a phase-space area greater than that occupied by the separatrix. In fact, the reason why the ${ \mathcal J }=0$ approximation works relatively well in our calculation is that the area engulfed by the $\infty$-shaped curve is never big. Therefore, a more rigorous treatment of adiabatic theory can be carried out by assuming a non-zero, but nevertheless small, value of ${ \mathcal J }$.

Irrespective of the details associated with computation of the Hamiltonian (12) for ${\varphi }_{\mathrm{res}}$-type commensurabilities, the results of the N-body simulation shown in Figure 12 as well as the semi-analytic calculations presented in Figure 13, point to the fact that resonantly protected trajectories that circulate in the longitude of perihelion naturally arise within the framework of the Planet Nine hypothesis. Correspondingly, given sufficient data, characterization of the specific semimajor axes that such trajectories occupy may provide key constraints on the present-day orbital state of Planet Nine.