The Secular Dynamics of TNOs and Planet Nine Interactions

We consider an outer planet (Planet Nine) and a TNO orbiting the Sun, illustrating the configuration of the system in Figure 1. The mass of the TNO is much smaller than that of Planet Nine and our Sun, and thus the TNO can be treated as a test particle. We allow the orbit of Planet Nine to be eccentric and misaligned with the inner orbit, different from the circular restricted case.

Zoom In Zoom Out Reset image size

In the nonhierarchical configuration, when a/a₉ ≳ 0.1, the usual expansion in the semimajor axis ratio is not a good approximation. However, in the case when the outer perturber is much less massive than the center body, the perturbation from the outer companion is not strong enough to destabilize the system; i.e., the semimajor axes of the orbits are almost constant. In this case, one can consider the long-term secular evolution of the system by averaging out the fast varying orbital phase of the inner and outer orbits.

Specifically, the secular (averaged) Hamiltonian of the interaction energy can be expressed as

$$\begin{eqnarray}&&{H}_{\sec ,0}=-\displaystyle \frac{{{Gm}}_{9}}{4{\pi }^{2}}{\int }_{0}^{2\pi }{\int }_{0}^{2\pi }\left(\displaystyle \frac{1}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}| }-\displaystyle \frac{{\boldsymbol{r}}\cdot {{\boldsymbol{r}}}_{9}}{{r}_{9}^{3}}\right)\,{dl}\,{{dl}}_{9},\end{eqnarray} \tag{ 1 }$$

where r and r₉ are the distance to the test particle and the perturber from the central object, m₉ is the mass of the perturber as illustrated in Figure 1, and l and l₉ are the mean anomaly of the test particle and the perturber, respectively. Here ${\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}$ can be expressed as

$$\begin{eqnarray}&&| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{9}| =\sqrt{{r}^{2}+{r}_{9}^{2}-2{{rr}}_{9}\cos \phi },\end{eqnarray} \tag{ 2 }$$

and r, r₉, and ϕ can be expressed by orbital elements.

The evolution of the TNO's orbit can be well described using the averaged Hamiltonian, which converges in most configurations even for crossing orbits, as discussed in Gronchi & Milani (1999) and Gronchi (2002). For illustration, we compare the secular evolution with the N-body simulations in Figure 20 in Appendix A.

The TNOs in the outer solar system undergo perturbations from Planet Nine, as well as the known inner giant planets. The TNOs are quite far away from the inner giant planets; thus, the secular effects of the inner giant planets on the TNOs can be well approximated based on the Hamiltonian to the second order in the ratio of the TNO to giant planet semimajor axes. Combining the secular Hamiltonian for the interaction energy of the TNO with a perturber (Equation (1)), the secular Hamiltonian can be expressed as (e.g., Kaula 1964; Murray & Dermott 1999)

$$\begin{eqnarray}&&{H}_{\sec ,1}={H}_{\sec ,0}-\displaystyle \frac{1}{8}\displaystyle \frac{{GM}}{a}\displaystyle \frac{(3{\cos }^{2}i-1)}{{(1-{e}^{2})}^{3/2}}\displaystyle \sum _{i=1}^{4}\left(\displaystyle \frac{{m}_{{pi}}{a}_{{pi}}^{2}}{{{Ma}}^{2}}\right),\end{eqnarray} \tag{ 3 }$$

where M is the mass of the Sun; a, e, and i are the semimajor axis, eccentricity, and inclination of the TNO; and m_pi and a_pi are the masses and semimajor axes of the inner giant planets. Here we assume that the inner planets are coplanar with the orbit of Planet Nine.

The orbit of Planet Nine is also perturbed by the inner giant planets, which causes precession. The precession rate of Planet Nine's orbit in the low-inclination regime can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{d{\omega }_{9}}{{dt}}=\displaystyle \frac{3}{2}\displaystyle \frac{\sqrt{{GM}}}{{a}_{9}^{3/2}{(1-{e}_{9}^{2})}^{2}}\displaystyle \sum _{i=1}^{4}\left(\displaystyle \frac{{m}_{{pi}}{a}_{{pi}}^{2}}{{{Ma}}_{9}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Omega }}}_{9}}{{dt}}=-\displaystyle \frac{3}{4}\displaystyle \frac{\sqrt{{GM}}}{{a}_{9}^{3/2}{(1-{e}_{9}^{2})}^{2}}\displaystyle \sum _{i=1}^{4}\left(\displaystyle \frac{{m}_{{pi}}{a}_{{pi}}^{2}}{{{Ma}}_{9}^{2}}\right),\end{eqnarray} \tag{ 5 }$$

where m₉, a₉, e₉, ω₉, and Ω₉ are the mass, semimajor axis, eccentricity, argument of pericenter, and longitude of the ascending node of Planet Nine separately.

The Hamiltonian of the TNOs (Equation (3)) is implicitly expressed in terms of the canonical variables of the TNOs,

$$\begin{eqnarray}\begin{array}{rcl}q & = & (M,\omega ,{\rm{\Omega }})\\ p & = & (\sqrt{{GMa}},\sqrt{{GMa}(1-{e}^{2})},\sqrt{{GMa}(1-{e}^{2})}\cos (i)),\end{array}\end{eqnarray} \tag{ 6 }$$

where M is the mean anomaly of the TNO, $J\,=\sqrt{{GMa}(1-{e}^{2})}$ and ${J}_{z}=\sqrt{{GMa}(1-{e}^{2})}\cos (i)$ are the angular momentum and z-component of angular momentum of the TNO, and $L=\sqrt{{GMa}}$ is a constant in the secular regime.

To obtain the alignment of the orbit of the TNO relative to that of Planet Nine, we transfer the coordinates of the TNOs to the difference in the argument of pericenter and the longitude of the ascending node between the TNO and Planet Nine. Using the type-III generating function, G₃ = −(Δω+ω₉)J −(ΔΩ+Ω₉)J_z, we transform to new canonical angles Δω = ω − ω₉ and ΔΩ = Ω−Ω₉. Then, the new Hamiltonian becomes

$$\begin{eqnarray}&&{H}_{\sec }={H}_{\sec ,1}-{\dot{\omega }}_{9}J-{\dot{{\rm{\Omega }}}}_{9}{J}_{z}\,.\end{eqnarray} \tag{ 7 }$$

The secular result following Equation (7) is a good approximation until close encounters occur between the TNO and Planet Nine or the inner giant planets. For illustration, Figure 2 compares the secular and N-body results, which show the evolution of a TNO starting with a = 368.75 au, e =0.867, Δϖ = 180°, and i = 10° relative to the ecliptic and with Planet Nine and the inner four giant planets all coplanar in the ecliptic plane. The semimajor axis and eccentricity of Planet Nine are a = 500 au and e = 0.6. The blue crosses are the secular results, and the red dashed lines are the N-body results. To isolate the effects of close encounters with Planet Nine, we substitute the inner giant planets (Jupiter, Saturn, Uranus, and Neptune) with an equivalent J₂ term in the N-body simulation. We use the package Mercury for the N-body simulation, with the "hybrid" Wisdom–Holman/Bulirsch–Stoer integrator (Press et al. 1992; Wisdom & Holman 1992; Chambers 1999) and a time step of dt = 3000 days, which is roughly 5% of Neptune's orbital period.

Zoom In Zoom Out Reset image size

Figure 2 shows that the secular result is a good approximation up to ∼50 Myr, at which point the TNO has a close encounter with Planet Nine, causing changes in the semimajor axis of the TNO. However, the secular effects are not immediately suppressed after the close encounter. In particular, the libration of Δϖ around 180° shown in the secular results can still be observed in the N-body results after the close encounter until the end of simulation. Neither Δω nor ΔΩ librates for both the secular and N-body results. The secular results agree qualitatively with the N-body results after the change in the semimajor axis. This is due to the weak dependence of the secular interactions on the semimajor axis for TNOs with large semimajor axes (≳300 au), as shown in Figure 4 and Hadden et al. (2018).

Note that the precession of the orbits increases the chance of close encounters between TNO and Planet Nine, which causes the secular approximation to deviate from the N-body results. Thus, the secular integration is a better approximation for the three-body interactions without J₂ precession, as illustrated in Appendix A. For instance, the J₂ precession timescale due to the inner giants in Ω for Planet Nine is 13 Gyr, and for the TNO, it is 276 Myr.

2.2.1. Alignment in Pericenter Orientation

It has previously been found (e.g., Hadden et al. 2018 and our Section 3) that particles starting with anti-aligned pericenters are more likely to survive and play an important role in sculpting the overall orbital architecture of the TNOs. Thus, we start the analysis by focusing on the secular evolution of an (initially) anti-aligned population. We consider a set of 500 test particles initialized with 150 au < a < 550 au, 30 au < r_p < 50 au, i = 10°, ϖ = 180°, and ω uniformly distributed between 0° and 360°. As in the illustrative example in Figure 2, Planet Nine lies in the ecliptic plane, together with the inner four giant planets, and we substitute the inner giant planets with the equivalent J₂ potential. The semimajor axis and eccentricity of Planet Nine are set to be a₉ = 500 au and e₉ = 0.6, and the longitude of pericenter of Planet Nine is ϖ₉ = 0°. Thus, the pericenters of the test particles are initially anti-aligned with that of Planet Nine. We follow the particles' secular evolution by integrating equations of motion generated from the secular Hamiltonian, Equation (7). Details of our numerical method are given in Appendix A.

Figure 3 shows the secular evolution of the TNOs. Dividing the particles into semimajor axis regions, we notice that Δϖ circulates when TNO semimajor axes are low: a ≲ 200 au (upper left panel), where we illustrate the trajectories of a few representative particles. This is the region where the J₂ precession dominates. For instance, the precession caused by the J₂ term is roughly $\dot{\varpi }=4200^\circ /\mathrm{Gyr}$, calculated following Equation (4). On the other hand, the change rate of ϖ caused by Planet Nine is $\dot{\varpi }=-360^\circ /\mathrm{Gyr}$, following the averaged Hamiltonian (Equation (1)) without the J₂ term in a coplanar configuration. Therefore, the J₂ precession leads to the circulation in ϖ in the positive direction.

Zoom In Zoom Out Reset image size

When a increases, the dynamical influence due to the inner giant planets becomes weaker and the perturbation due to the outer planet becomes stronger. Thus, the J₂ term no longer dominates the evolution. As shown in the upper middle and left panels of Figure 3, when a increases to ∼210–300 au, Δϖ stops circulating and starts to librate around 180°. In the upper panels, where a < 250 au, the inclinations of the TNOs stay low for the entire 4 Gyr simulation.

As a further increases (a ≳ 250 au), the initial eccentricities approach closer to unity because we chose a fixed range of initial pericenter distances (30 au < r_p < 50 au). This is shown in the lower panels in Figure 3: as the eccentricity increases, the J₂ precession increases, while the precession due to Planet Nine decreases, leading to the circulation in Δϖ at high eccentricity.

Thus far, the secular dynamics of anti-aligned particles we have described are qualitatively the same as those seen in the strictly coplanar case previously studied by Beust (2016), Batygin & Morbidelli (2017), and Hadden et al. (2018). In particular, an island of librating trajectories appears at a critical semimajor axis where the apsidal precession induced by Planet Nine is able to balance that from the solar system giants.

Figure 3 shows that Planet Nine can also excite the inclinations of some TNOs, a dynamical effect that, by construction, is absent from previous coplanar studies. When TNOs reach very high eccentricity orbits, Planet Nine is able to excite extreme inclinations, often leading to orbital flips. We will study the dynamics of these orbital flips in more detail in Section 4.

During the flips, the pericenter orientation of the TNO in the ecliptic plane (θ_e, the geometrical longitude of pericenter³ ) librates (as illustrated in Section 4). This is similar to the coplanar flips shown in the hierarchical octopole limit (Li et al. 2014b). In contrast to the nearly coplanar cluster, for high-inclination orbits (i > 140°), Δϖ circulates, as can be seen in the lower panels of Figure 3. This differs from Figure 10 in Batygin & Morbidelli (2017), where Δϖ remains confined as the inclination becomes larger. This is because the objects selected in Batygin & Morbidelli (2017) only stay briefly in the high-inclination region, before circulation is completed. We can see both types of objects in our secular and N-body simulations (e.g., some trajectories in Figures 3 and 26); some of them stay only briefly in the retrograde stage, and some of them stay much longer. The circulation of Δϖ when the TNO stays at high inclinations can also be seen in Figure 16. During the circulation of Δϖ, the eccentricity of the TNO becomes high, and the pericenter distance is reduced to the point where TNOs may be scattered and ejected by the inner giant planets, an effect that is not captured by our secular calculations. However, we do find instances in our our N-body simulations where TNOs survive such phases of high eccentricity and inclination, during which Δϖ rapidly circulates. Figure 26 illustrates two such examples, near t ∼1.5 Gyr for the TNO in red and t ∼ 3.5 Gyr for the TNO in yellow.

In addition, we note that in the upper right panel of Figure 3, where 250 au < a < 280 au, there are two teal objects displaying an interesting evolution, where their eccentricities vary with large amplitudes. There are only two teal objects in this regime; one of them has ϖ centered in and oscillating around 0°.

To illustrate the libration of the geometrical longitude of pericenter in the ecliptic plane for some of the flipped orbits, we show in Figure 4 the evolution of TNO pericenter distances versus Δθ_e in the same semimajor axis ranges as those in the lower panels of Figure 3, where the inclination of the TNOs can be excited. At high inclinations (e.g., i ≳ 150°), Δθ_e librates, while Δϖ circulates. We note that θ_e ∼ ϖ when inclinations are low (i ∼ 0°), and θ_e ∼ 2Ω−ϖ when TNOs counterorbit with regard to Planet Nine (i ∼ 180°). The libration of θ_e is consistent with the libration of 2Ω−ϖ noted by Batygin & Morbidelli (2017) for high-inclination orbits.

Zoom In Zoom Out Reset image size

Out of the 500 particles, 105 have their inclinations excited to retrograde configurations during the 4 Gyr simulation. It is more likely for the inclination to flip if the particles start with large semimajor axes and small pericenter distances, but the detailed dependencies are more complicated, as we demonstrate in Section 4.

To illustrate the clustering of high-inclination TNO orbits, Figure 5 presents trajectories in the plane of (Δθ_e, i), color-coded in the pericenter distance. There is a moderate clustering of trajectories around Δθ_e ∼ 90° and ∼270° when the inclinations are high, i ≳ 60°–120°, and the pericenter distances are low, ≲80 au. Around i ≳ 140°, Δθ_e is clustered around 150°–210°.

Zoom In Zoom Out Reset image size

Combining all the particles in the secular simulations with different semimajor axes, we show the scatter in semimajor axis versus geometrical longitude of pericenter in Figure 6. The top panel is color-coded by the pericenter distance, and the bottom panel is color-coded by the inclination of the TNOs. To approximate detectable objects, we select only particles with t > 3 Gyr. In addition, we select only particles with 30 au <r_p < 80 au to focus on the closer-in objects, which are more likely to be detectable. All 500 particles are started with anti-aligned pericenters (Δϖ = 180°) and inclinations of 10° from the ecliptic. Each point corresponds to a test particle, with snapshots taken at 1 Myr time steps. There is a clear alignment in Δθ_e that begins when the semimajor axes of the TNOs are around ≳200–300 au, corresponding to the libration region in Figure 3. In addition, when we consider only the lower-inclination objects (≲40°), the alignment near Δθ_e ∼Δϖ ∼ 180° is stronger, particularly when a ≳ 300 au. This shows that the observed clustering can be produced by pure secular effects alone when starting with TNOs in the near-coplanar configuration. We note that the clustering in Δθ_e is stronger than that in Δϖ because Δϖ no longer librates for the high-inclination population. This is consistent with what we find in Figures 3–5. The clustering in both Δω and ΔΩ is very weak and cannot be seen in the scatter plots (as shown in the histogram plots in Figure 8). Thus, we do not include them here.

Zoom In Zoom Out Reset image size

2.2.2. Alignment in Argument of Pericenter ω and Longitude of Node Ω

We now consider the alignment of the orbital plane in argument of pericenter and longitude of ascending node. We focus on TNOs with a ≳ 300 au, which allows the inclination of the TNOs to be excited to large values. Figure 7 presents the evolution of the TNOs in the plane of (Δω, i) and (ΔΩ, i), color-coded in the pericenter distance r_p. The figure illustrates "looping" trajectories followed by TNOs that reach high inclination, resulting in the clustering of ω near Δω ∼ 0° and 180° and Ω near Δω ∼ 90° and 270° among the high-inclination TNOs when the pericenter distance is low, ≲80 au.

Zoom In Zoom Out Reset image size

In the upper row of Figure 7, the dynamical regions where the trajectories librate around Δω ∼ 90° and ∼270° for 40° <i < 90° and 90° < i < 140° are analogous to the quadrupole-order Kozai–Lidov resonant regions, which are seen for both interior and exterior test particles (e.g., Kozai 1962; Lidov 1962; Naoz et al. 2017; Vinson & Chiang 2018). These resonant regions are also shown in the secular study by Saillenfest et al. (2017) on TNOs perturbed by Planet Nine using surface of sections. However, due to the close separation between the TNOs and Planet Nine, some of the hierarchical approximation breaks down, as shown in Batygin & Brown (2016a). In the lower row of Figure 7, in the (ΔΩ, i) plane, the regions where the trajectories librate around ΔΩ ∼ 180° are analogous to the octopole Kozai–Lidov resonances, where i can be excited from near zero inclination. The cat-eye-shaped regions centered around inclination i ∼ 90° and ΔΩ ∼ 0° and ∼ 180° are also analogous to the octopole Kozai–Lidov resonances when the interaction energies are higher (e.g., Figure 4 of Li et al. 2014a, panel H = −1 and H = −0.5).

We summarize the orbital alignment in Figure 8, which shows histograms of Δθ_e, Δω, and ΔΩ. We choose only particles with t > 3 Gyr and 30 au < r_p < 80 au, since these objects are more likely to be detectable. We break up the particles into three inclination bins in Figure 8: the upper panel selects particles with a > 250 au and i < 40°, the middle panel focuses on the high-inclination (60° < i < 120°) particles, and the lower panel focuses on the retrograde particles with high inclination (i > 140°). For the low-inclination population (i < 40°), there is a strong clustering around Δθ_e ∼ 180°, and the clusterings in Δω and ΔΩ are weak. The high-inclination population, shown in the middle panel, exhibits a clear deficit of particles around Δθ_e ∼ 180° and peaks near Δθ_e = 90° and 270°. We note that the clustering in Δϖ ∼ 180° for the high-inclination TNOs (middle panel of Figure 8) is missing due to the detectability cut, but clustering in Δϖ ∼ 180° can be seen for high-pericenter objects in Figure 5. In addition, there is an excess of objects with ΔΩ ∼ 90° and 270° and Δω ∼ 0° and 180°, which illustrate clustering of high-inclination orbits for the low-pericenter TNOs. The counterorbiting particles (i > 140°) show a double-peaked clustering near Δθ_e ∼ 150° and 210°, with a slight deficit near Δθ_e ∼ 180°, while Δω and ΔΩ are roughly uniform.

Zoom In Zoom Out Reset image size

We note that we make the detectability cut at r_p < 80 au to facilitate comparison with observational results, in order to clarify the role of secular interactions in producing the observed clustering of the TNO orbits. Many dynamical effects, e.g., scatterings with Neptune, Planet Nine, and MMR, are neglected in the secular approach. Thus, we do not intend to reproduce the full dynamical interactions between the TNOs and Planet Nine using the secular methods. In addition, we note that when including randomly initialized TNOs with ω and Ω uniformly distributed, the secular dynamics is similar to that of the exact coplanar case, where there are two clusterings of ϖ around 0° and 180°. The 0° clustering is unstable if one considers the full dynamics, including close encounters with Neptune and MMR with Planet Nine (e.g., Hadden et al. 2018).

For simplicity, we start with the secular approximation. To characterize the range of parameter space where the TNOs could flip, we investigate the flips in the plane of TNO semimajor axes and initial eccentricity using the secular integration in Figure 13. Similar to the secular simulations discussed above, we set a₉ = 500 au, e₉ = 0.6, and i₉ = ω₉ = Ω₉ = 0. We set the initial longitude of pericenter of the TNO to be ϖ₀ = 180° (ω₀ = π, Ω₀ = 0) here, since the surviving TNOs mostly have Δϖ₀ ∼ 180° (as shown in Figure 10). The inclinations of the TNOs are all set to be 5°, slightly misaligned from the ecliptic plane. The parameter space that allows the TNOs to flip for the anti-aligned TNOs does not depend sensitively on the initial mutual inclination, as long as it is nonzero and near-coplanar.

Zoom In Zoom Out Reset image size

To understand the flips better, we break down the problem into pieces and consider the three-body interactions first, without the J₂ precession. Including only the central star, a TNO, and the perturbing Planet Nine, the orbit of the TNO can be easily flipped both when the TNO is inside the orbit of Planet Nine and when it is farther. For illustration, we represent the runs that could not flip using plus signs and then use crosses to represent the runs that can flip when we ignore J₂ precession.

Next, we marked with circles the runs that can flip in the presence of the J₂ precession. Including J₂ precession, the flip of the orbits can be suppressed. This is determined by the libration and J₂ precession timescales of the TNOs. To compare the timescales, we color-coded the crosses using the ratio of the libration timescale of Δθ_e and the J₂ precession timescale. To obtain the timescales, we numerically calculated the libration timescale of the geometrical longitude of pericenter following the secular integration of the three-body interaction, and we calculated the J₂ precession due to the inner four giant planets analytically based on Equation (4).

When the libration timescale exceeds 4 Gyr, we color the crosses red. As expected, the runs can still flip when the libration timescale is shorter than the J₂ precession timescales in general. Most of the TNOs can still flip when they are farther from the inner giant planets with a ≳ 300 au and r_p ≳ 100 au. Note that we only focus on the anti-aligned configurations here, since these TNOs are more likely to survive. The TNO orbits are less likely to flip if their pericenters are aligned with that of Planet Nine, analogous to the hierarchical limit (Li et al. 2014b). The general flip condition in the nonhierarchical limit is quite complicated and beyond the scope of our paper.

As illustrated in Figure 13, the overall dynamics can be divided in the following three regions.

• (1)  
  Small semimajor axis region (a ≲ 150 au), where the precession timescale is much shorter than that of the Δθ_e libration timescale. Both the flip and the libration in Δθ_e can be suppressed.
• (2)  
  Large semimajor axis region (a ≳ 150 au) within the two black lines inside 30 au ≲ r_p < ≲ 100 au in Figure 13, where the two timescales are comparable. The flips are suppressed while the libration in Δθ_e still persists.
• (3)  
  Large semimajor axis region (a ≳ 150 au), where the libration timescale is much shorter. Both the libration and the flip remain.

To illustrate the flips in more detail, we show examples of arbitrarily selected trajectories in these three regions below.

Figure 14 illustrates the evolution in region 1, where the semimajor axis of the TNO is small, and thus the J₂ term from the inner four giant planets dominates over the perturbation from the outer Planet Nine. The left panel shows the case with J₂ precession, and the right panel shows the case without J₂ precession. Without J₂ precession, the flip is analogous to that in the hierarchical limit (Li et al. 2014b). The precession timescale is much shorter compared with the Δθ_e libration timescale (as shown in Figure 13), and thus the libration of Δθ_e is suppressed when the J₂ term is included. The flip of the orbit is also suppressed in this case.

Zoom In Zoom Out Reset image size

Figure 15 illustrates the evolution in region 2, where the semimajor axis of the TNO is larger, and thus the J₂ precession from the inner four giant planets becomes comparable to the perturbation from the outer Planet Nine. Similar to Figure 14, the left panel shows the case with the J₂ term, and the right panel shows the case without the J₂ term. The J₂ precession timescale is similar to the Δθ_e libration timescale (as shown in Figure 13). The libration of Δθ_e is not suppressed when the J₂ term is included; however, the flip of the orbit is suppressed in this case. Orbits lying in this region are detectable with r_p ≲ 100 au, and they contribute to the alignment of the orbits around Δθ_e ∼ 180°.

Zoom In Zoom Out Reset image size

Finally, Figure 16 illustrates the evolution in region 3, where the eccentricity of the TNO is lower and the perturbation from the outer Planet Nine is more dominant. Again, the left panel shows the case with the J₂ term, and the right panel shows the case without the J₂ term. The J₂ precession timescale is lower than the Δθ_e libration timescale (as shown in Figure 13), and neither the libration of Δθ_e nor the flip of the orbit are suppressed when the J₂ term is included. Orbits lying in this region contribute to the flipped orbits.

Zoom In Zoom Out Reset image size

In addition, the TNO could stay in the high-inclination regime for many cycles of low-amplitude 2Ω−ω libration, and ϖ circulates during this high-inclination phase. This is similar to the N-body results presented in Figure 23, but this is different from the particles selected in Figure 10 of Batygin & Morbidelli (2017), where the particles fall back to the low-inclination regime quickly and ϖ is still confined. We note that there are also particles that only spend one libration cycle at high inclination in our simulation, similar to Figure 10 of Batygin & Morbidelli (2017).

Interestingly, we notice that orbits within 30 au ≲ r_p ≲100 au and a ≳ 200 au inside the solid and dashed black lines in Figure 13 do not flip, and the geometrical longitude of pericenter Δθ_e still librates in the presence of the J₂ term, since the two timescales are comparable to each other. These TNOs lie within our selection criteria based on observational limits, and they lead to the alignment of the low-inclination TNO orbits with Δθ_e ∼ 180°, as shown in the N-body and secular simulations discussed in the previous sections.

The parameter space of clustered low-inclination TNOs corresponds to the central libration region in the plane of (Δϖ, r_p) (e.g., Figure 3), near Δϖ = 180°. Thus, characterizing the dependence of the libration region on the properties of Planet Nine can help constrain the possible outer planet based on the detected clustering of low-inclination orbits.

In Figure 17, we show the corresponding pericenter distances of the fixed points as a function of semimajor axis for different Planet Nine orbital parameters. The fixed point of the libration region is located at Δϖ = 180°, and we numerically calculate the pericenter distances of the fixed point at different semimajor axes by searching for eccentricity corresponding to the maximum energy when Δϖ = 180°.

Zoom In Zoom Out Reset image size

The libration regions exist around the fixed points. The secular resonances around Δϖ = 180° disappear for low TNO semimajor axes (e.g., upper left panel of Figure 3). We marked the minimum semimajor axes when the libration regions appear with black circles. There are no fixed points when the semimajor axis is small, in particular for wider and more circular Planet Nine's orbit. In other words, when the orbit of Planet Nine is more circular, the anti-aligned libration appears only for TNOs with larger semimajor axes. In addition, when Planet Nine is more distant, it also requires a more eccentric Planet Nine orbit to produce the alignment for the closer-in TNOs. This is consistent with the N-body results shown in Brown & Batygin (2016), where Planet Nine favors a more eccentric orbit around a ∼ 500 au to produce the clustering.

As shown in the previous section, the flip of the orbits depends on the existence of the libration region of Δθ_e ∼ 180°, and this itself depends on the properties of Planet Nine. Thus, the retrograde TNO orbits and the signatures of the inclination distribution can provide valuable constraints on the properties of any possible outer planet. In this section, we illustrate the dependence of the inclination distribution using N-body simulations. First, we present the inclination distribution of the TNOs under perturbations from a Planet Nine with m₉ = 10 M_⊕, a₉ = 500 au, e₉ = 0.6, and i₉ = 3°, similar to that included in Section 3. We find that some high-inclination objects have small pericenter distances (r_p < 30 au). Thus, to obtain a more accurate TNO inclination distribution, we rerun the N-body simulation including Jupiter, Saturn, and Uranus as massive point particles, instead of their equivalent J₂ potential.

In Figure 18, we show the scatter in inclination as a function of semimajor axis from our N-body simulations. The colors represent the density of the TNOs (log of the number of TNOs per au per degree). Particles are included with t > 3 Gyr, plotted at snapshots taken every 1 Myr. To focus on the high-inclination population, we selected TNOs with i > 30°, and to select those more likely to be observable, we require r_p < 80 au. We find that the TNOs become retrograde (i > 90°) when a ≳ 30 au. This is closer than in the secular approximation, since the semimajor axes of the particles can drift due to close encounters with Neptune and overlap of MMR, both of which are not captured in the secular approximation.

Zoom In Zoom Out Reset image size

In addition to the limiting semimajor axis for which TNO orbits can be flipped, the inclination distribution shows other interesting features. As shown in the density map in Figure 18, within a ≲ 100 au, there is an overdensity of inclined orbits around 150°. This is due to the nature of the flips when the semimajor axis is small, crosses over 90° fast, and then spends more time near ∼150° before flipping back. As a continues to increase (a ≳ 300 au), the distribution of the inclination is more uniform.

The inclination distribution depends on the properties of Planet Nine. For instance, if Planet Nine is less massive (∼M_⊕), the perturbations on the TNOs are weaker, and the J₂ precession due to the inner planets dominates. Thus, it is difficult to flip the TNO orbits, and none of the TNO orbits are flipped over 90°.

Since the flips are associated with the libration of the geometrical longitude of pericenter (Δθ_e) around 180°, and since this only occurs for Planet Nine orbits that have high eccentricity and are wide (a ∼ 400–800 au; see Figure 17), we expect that the TNO orbits can flip only for such Planet Nine orbits. To test this, we performed more N-body simulations to calculate the fraction of TNOs that can have their inclinations excited using different Planet Nine properties. We included a 10 M_⊕ Planet Nine and varied the Planet Nine orbital parameters in these simulations.

Figure 19 shows the fraction of TNOs that have inclinations over 60° and 90° at 4 Gyr. Specifically, it is more likely to flip the TNO orbits for an eccentric Planet Nine (e ≳ 0.4) with a large semimajor axis of ∼400–800 au. This agrees with our expectation based on the appearance of the libration region around Δθ_e ∼ 180°. The recent discovery and dynamical analysis of high-inclination objects, e.g., 2015 BP₅₁₉ (Becker et al. 2018), supports the existence of possible outer planets. Future surveys of the inclination distribution of Centaurs can help further constrain the parameters of possible outer planets (e.g., Petit et al. 2017; Lawler et al. 2018).

Zoom In Zoom Out Reset image size

In this section, we present the detailed TNO orbital evolution starting in the near-coplanar configuration as mentioned in Section 3 based on the N-body results. First, we illustrate the presence of high-order MMRs during the evolution of TNOs. As mentioned in the main text, Figure 9 shows that the TNOs could survive outside of the lower-order MMRs. To illustrate any presence of higher-order MMRs, we plot the mean anomaly of the TNOs when that of Planet Nine is zero in Figure 21.

Zoom In Zoom Out Reset image size

The evenly spaced blank regions in Figure 21 illustrate the presence of MMRs due to commensurability between the TNO and Planet Nine. Two out of the three randomly selected TNOs spend only a negligible fraction of time in MMR, while the TNO in the middle panel spends around ∼20%–30% of the time in MMR. All of them stay out of MMR for a large fraction of time, indicating the importance of pure secular interactions in their orbital evolution. Performing a similar analysis for 100 surviving TNOs, we find that more than ∼80% of them remain outside all MMRs for more than ∼80% of the time.

In addition, we adopt a more systematic approach, where we cut the time series of M_TNO at M₉ = 0 to 1 and 10 Myr segments and perform Kolmogorov–Smirnov (ks) tests in each segment. If the ks test shows that the M_TNO agrees with a uniform distribution, we mark the TNO to be out of MMR for the time segment. Using time segments of 1 Myr, 95% of the TNOs spend more than 80% of the time outside of MMR, and using time segments of 10 Myr, 85% of the TNOs spend more than 80% of the time outside of MMR. We note that this approach only gives a rough estimate and is very sensitive on the segment widths.

Next, we illustrate the clusterings of the TNO orbits in more detail and compare with the secular results. The alignment of the test particles in pericenter (Δϖ and Δθ_e) as a function of particle semimajor axis is shown in Figure 22. We record the longitude of pericenter of the test particles every 1 Myr for t > 3 Gyr. To illustrate the dependence of the alignment on different parameters of the system, we include three different color codes. The upper left panel of Figure 22 color-codes the TNOs in the pericenter distance, and it shows the orbital distribution of the TNOs in the plane of (a, Δθ_e). It illustrates that there is a clear alignment of the test particles around the geometrical longitude of pericenter, Δθ_e ∼ 180°, for small-pericenter objects (30 au < r_p < 80 au), consistent with Brown & Batygin (2016) and the observational results of Trujillo & Sheppard (2014). Particles with smaller pericenter distances, r_p ≲ 30 au, exhibit faster precession in the longitude of pericenter, suppressing any clustering, and the clustering is weak for longer pericenter distance TNOs, ≳80. In addition, the anti-aligned particles exhibit drifts in their semimajor axes at approximately constant pericenters. This is due to close encounters with Neptune when the particles are close to their pericenter, as shown in the coplanar case investigation by Hadden et al. (2018).

Zoom In Zoom Out Reset image size

Figure 22 shows that there is no clustering around Δθ_e ∼ 0°, which differs from the secular case initialized with a uniform distribution of ϖ. This is because the TNOs with Δϖ₀ ∼ 0 are ejected due to instability caused by the overlap of MMRs, as discussed in Hadden et al. (2018). This is illustrated in the upper right panel in Figure 22, which is color-coded in Δϖ₀. It shows that most of the surviving TNOs that are clustered around Δθ_e ∼ 180° started with Δϖ₀ ∼ 180°. Meanwhile, there is a population of particles that exhibit drifts in pericenter distances at constant semimajor axes (most clearly seen as the multicolored vertical paths in the top left panel of Figure 22), which arise due to secular effects.

The lower panels of Figure 22 color-code the TNOs in inclination. The lower left one shows the alignment in Δϖ. Here Δϖ clusters around ∼180°, and the clustering is strongest when the test particle inclinations are low, ≲40°. The majority (∼90%) of the anti-aligned (120° < ϖ < 240°), small-pericenter (30 au < r_p < 80 au) particles maintain a low inclination throughout their 4 Gyr evolution. On the other hand, the lower right panel shows the alignment in the geometrical longitude of pericenter (Δθ_e), color-coded in inclination. In contrast to the result for Δϖ, Δθ_e shows tight clustering, even for high-inclination TNOs. This is because Δθ_e librates during the flips of the orbits, similar to the flips of the inner orbit in the hierarchical three-body interactions (Li et al. 2014b; see more details in Section 4).

As shown in the secular results of Section 2, the high-inclination population also shows interesting clustering in the orbital orientation. To illustrate this, we plot in Figure 23 the alignment in different orbital orientations as a function of semimajor axis. All of the long-lived (t > 3 Gyr), high-inclination (60° < i < 120°) particles are selected: there are 303 long-lived TNOs that reached above 60°, and 255 of them reached retrograde configurations. The clustering in longitude of node, ΔΩ, is shown in the upper left panel; argument of pericenter, Δω, in the upper right panel; longitude of pericenter, Δϖ, in the lower left panel; and geometrical longitude of pericenter, Δθ_e, in the lower right panel.

Zoom In Zoom Out Reset image size

Consistent with the secular results, there is a strong deficit of particles with Δϖ ∼ 180° and Δθ_e ∼ 180°. In addition, there are clusters around Δϖ ∼ 90° and 270°, and the clusters around the geometrical longitude of pericenter are slightly tighter (Δθ_e ∼ 90° and 270°) for the low-pericenter TNOs. This is very different from the low-inclination population. The clusterings in Δω ∼ 0° and 180° and ΔΩ ∼ 90° and 270° are very similar to the secular results. The existence of such clustering in future TNO detection can help constrain the orbital parameters of any outer planet.

To illustrate the evolution of the particles and Planet Nine, we plot Planet Nine and selected representative TNO orbital elements as a function of time in Figures 24–26. Figure 24 shows the evolution of Planet Nine, Figure 25 shows the evolution of some examples of low-inclination objects, and Figure 26 shows the evolution of some high-inclination objects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The semimajor axis and eccentricity of Planet Nine show little variation, as illustrated in Figure 24. The inclination of Planet Nine has low-amplitude, long-timescale variations due to interactions with Neptune. The red dashed lines in the lower three panels represent the precession due to the J₂ component by the inner giant planets. The J₂ precession agrees well with the changes of ω, Ω, and ϖ as a function of time.

Figure 25 shows the trajectories of low-inclination TNOs. The semimajor axes of the particles are perturbed by Planet Nine and Neptune, which cause large period-ratio variations. The eccentricity of the particles also varies throughout the 4 Gyr of integration, which allows the pericenter distance to oscillate. The inclinations stay low (≲30°) within the 4 Gyr integration. Here Δϖ librates for the particles with large semimajor axes, as illustrated by the particle represented by the yellow line. At smaller semimajor axes, Δϖ starts to circulate, and Δω and ΔΩ generally circulate within 4 Gyr.

Figure 26 shows the trajectories of some representative high-inclination TNOs. The semimajor axes of the particles are perturbed by Planet Nine and Neptune, which cause large period-ratio variations, similar to the low-inclination case. The eccentricity of the particles also varies throughout the 4 Gyr of integration, while the inclination can be excited to high values. Only ∼20% of the trajectories flip back to low inclination, ≲20°, within the 4 Gyr integration. Different from the low-inclination particles, the Δϖ of the particles does not always librate around 180°, even when the semimajor axes are large. Decreases in Δϖ can be seen when Δϖ does not librate, which explains the overdensity for Δϖ ∼ 90° compared with Δϖ ∼ 270°, since the long-lived particles starting near Δϖ ∼ 180° tend to reach ∼90° earlier and slightly more often than ∼270°. There are libration regions in Δω and ΔΩ as their inclinations increase that lead to clusterings in Δω and ΔΩ for high-inclination objects (60° < i < 120°).