Orbital Precession in the Distant Solar System: Further Constraining the Planet Nine Hypothesis with Numerical Simulations

For our numerical simulations, we use a version of the SWIFT N-body package (Levison & Duncan 1994) that is modified to include the effects of the galactic tide and stellar encounters (Kaib et al. 2019). The galactic tide is modeled with a radial component based on the Oort constants (Oort 1927) and a vertical term derived from the local density of matter in the Milky Way's disk (fixed at 0.1 M_⊙ pc⁻³ in our simulations; Levison et al. 2001). Stellar encounters are generated at random assuming the present-day mass function of Reid et al. (2002), the local velocity dispersions of Rickman et al. (2008), and a stellar density of 0.034 M_⊙ pc⁻³ (Kaib et al. 2019). Each stellar encounter is integrated directly for r < 1 pc. While these parameters have not remained constant in the solar neighborhood over the past 4 Gyr (e.g., Brasser et al. 2006; Kaib et al. 2011), we select this simplified model for a first-order approximation of the dynamic environment in the distant solar system.

Each of our dynamical simulations includes the four giant planets on their modern orbits, the hypothetical Planet Nine, and 1000 massless test particles. For completeness, we test two different proposed orbits (Batygin & Brown 2016a; Batygin et al. 2019) for Planet Nine (Table 1). Angular orbital elements for our test particles (specifically, argument of perihelion, longitude of ascending node, and mean anomaly, M) are drawn randomly from uniform distributions, while inclinations are selected from the following distribution (in order to resemble the modern distribution of KBOs; e.g., Brown 2001; Kaib & Sheppard 2016):

$$\begin{eqnarray}&&f\left(i\right)=\sin \left(i\right)\exp \left(-\displaystyle \frac{i}{2{\sigma }_{i}^{2}}\right).\end{eqnarray} \tag{ 1 }$$

Our simulations are designed to study the evolution of orbital clustering of ETNOs, IOCOs (Table 2), and clones (Table 3) of specific significant (e.g., Sheppard et al. 2019) TNOs under the influence of the external perturber (Batygin & Brown 2016a). For our simulations aiming to investigate the bulk populations of ETNOs and IOCOs, we select semimajor axes and perihelia randomly from the ranges given in Table 2. We also perform one set of control simulations that do not include Planet Nine and make use of the same initial conditions as in our ETNO, low-i batch. For our object clone runs, each test particle is randomly assigned a value within 1 au of the object's actual semimajor axis, 0.5 au of its pericenter, and 05 of its inclination (Table 3). We then perform 10 4 Gyr simulations of each initial condition set.

We follow a debiasing scheme similar to that of Brown (2017) and Brown & Batygin (2019) for determining the cumulative bias with respect to ϖ and Ω for survey detections in the distant solar system. This method makes the fundamental assumption that the probability of detecting a TNO at a given ϖ or Ω can be inferred from the complete MPC catalog of KBOs with q ≥ 30 au (for this work, queried on 2019 December 18). We first calculate the ecliptic latitude, longitude, and brightness for each KBO at the time it was discovered. This can be considered a sample of survey times, limiting magnitudes, and positions on the sky that might have been capable of detecting one of our distant TNOs had it possessed a different set of orbital angles. To calculate the longitude of perihelion bias for a given object of interest (see Section 2.3), we iterate through each possible ϖ/M combination (assuming a uniform distribution of angles) for each KBO detection in our aforementioned "sample of surveys." Through this process, we determine the range of ϖ and M for which each extreme object would be bright enough and at the correct position (within 1°) to be "detected" by any of the surveys. This allows us to build a probability distribution of detection with respect to ϖ for each distant TNO (an example for Sedna is plotted in Figure 1). As discussed in more depth in Brown (2017), we also attempt to account for several of our fundamental assumption's fallacies (though it should be noted that Brown & Batygin 2019 determined that the final debiased distributions do not depend strongly on these assumptions). Specifically, these are related to the overabundance of resonant plutinos detected near pericenter, the dependency of detection latitude on the KBO inclination distribution, and the fact that all detection surveys are not sensitive to extremely faint objects. In practice, we account for these inherent biases by limiting the underlying population of KBOs to those with a > 40 au (thus removing the plutinos), weighting each detection by the expected density of KBOs at that latitude (by converting Equation (1) to a latitudinal distribution, assuming σ_i = 14.9; Brown 2001), and only considering objects "detectable" at heliocentric distances <90 au. Our method for determining the bias in Ω follows in a similar manner. In that case, simulated detections have fixed values of ϖ, and Ω is the angle that is iterated over. Our results should be interpreted with an understanding that our debiasing scheme is inherently flawed but still a reasonable approach to capture the cumulative bias that is entangled in decades of different surveys utilizing various methodologies and reporting differing parameters.

Zoom In Zoom Out Reset image size

At this point, it is worth discussing the selection criteria used by authors analyzing orbital clustering in the distant solar system. Recent work by Brown & Batygin (2019) and Batygin et al. (2019) use delimiting values of q > 30 au, a > 250 au, and i < 40° (14 objects). However, this excludes the high-inclination (i = 54°, q = 35 au, a = 863 au) object 2015 BP519, discovered by the Dark Energy Survey (Dark Energy Survey Collaboration et al. 2016), the extreme orbit of which seems to be explained by interactions with Planet Nine (Becker et al. 2018). Additionally, the object uo5m93, reported in the full OSSOS data release (Bannister et al. 2018), would also meet these criteria, though the number of observations was insufficient to adequately constrain its orbit in accordance with the standards established by OSSOS. Finally, two additional objects' orbits have been determined and added to the MPC database since the work of Brown & Batygin (2019) and Batygin et al. (2019): 2013 RA109 (q = 46 au, a = 315 au, i = 65) and 2013 SL102 (q = 35 au, a = 863 au, i = 54°). While each of these new objects' longitude of perihelion loosely opposes the proposed value for Planet Nine, they are both at angles near the extreme of the apparent observed clustering

Sheppard et al. (2019) studied 13 objects that meet a more conservative minimum perihelia value (q = 40 au) and more relaxed semimajor axis limit (a > 150 au). As a result, that manuscript considers 2013 UT15 (q = 44 au, a = 196 au), 2013 GP136 (q = 41 au, a = 155 au), and 2000 CR105 (q = 44 au, a = 227 au) but does not include the objects with q ≃ 35–38 au (2014 FE72, 2007 TG422, 2013 RF98, and 2015 GT50) addressed by Batygin et al. (2019). While these imposed limits are indeed somewhat arbitrary, it is also readily apparent that bodies with perihelia near ∼30 au can interact strongly with Neptune. The degree to which Neptune perturbs these low-q objects decreases with increasing semimajor axis (Morbidelli et al. 2008). Thus, as an object's dynamical separation from Neptune is of utmost importance when considering orbital shepherding by an external perturber, we focus our dynamical simulations on detached objects satisfying the more conservative minimum perihelion distance of 40 au (Table 2). We take this one step further in our simulations designed to study object clones (Table 3) and select only objects with q > 45 au (thus, the objects in common to both the listings of Batygin et al. 2019 nand Sheppard et al. 2019). For consistency, when we compare our simulation-generated systems of clustered orbits to the observed objects using the test of Brown & Batygin (2019), we employ a similar object selection criterion: q > 30 au and a > 250 au. Therefore, with the inclusion of 2015 BP519, 2013 RA109, and 2013 SL102, our study considers the observed orbital clustering of 17 objects.

With few exceptions, consistent with Batygin et al. (2019), we find our P9b configurations (5.0 M_⊕ planet with a = 500 au, e = 0.25, and i = 200) to be significantly more successful than the corresponding P9a sets (10.0 M_⊕ planet with a = 700 au, e = 0.60, and i = 300). Thus, we focus the majority of our discussion in the subsequent sections on P9b simulations. We begin our analysis by repeating the Monte Carlo orbital clustering significance test of Brown & Batygin (2019). We first calculate the bias distribution function independently for each TNO as described in Section 2.2 and then perform a transformation to the orthogonal basis of canonical Poincaré variables (e.g., Morbidelli 2002). Beginning from the modified Delaunay variables,

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{j}={\mu }_{j}\sqrt{G({m}_{0}+{m}_{j}){a}_{j}},\quad {\lambda }_{j}={\omega }_{j}+{{\rm{\Omega }}}_{j}+{M}_{j},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{j}={{\rm{\Lambda }}}_{j}\left(1-\sqrt{1-{e}_{j}^{2}}\right),\quad {\gamma }_{j}=-{\omega }_{j}-{{\rm{\Omega }}}_{j},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{Z}_{j}={{\rm{\Lambda }}}_{j}\sqrt{1-{e}_{j}^{2}}\left(1-\cos {i}_{j}\right),\quad {z}_{j}=-{{\rm{\Omega }}}_{j}.\end{eqnarray} \tag{ 4 }$$

The Poincaré variables are defined as follows:

$$\begin{eqnarray}&&{x}_{j}=\sqrt{2{{\rm{\Gamma }}}_{j}}\cos -{\gamma }_{j},\quad {y}_{j}=\sqrt{2{{\rm{\Gamma }}}_{j}}\sin -{\gamma }_{j},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{p}_{j}=\sqrt{2{Z}_{j}}\cos -{z}_{j},\quad {q}_{j}=\sqrt{2{Z}_{j}}\sin -{z}_{j}.\end{eqnarray} \tag{ 6 }$$

When scaled by the quantity of Λ, a vector in the x/y plane points in the direction of the longitude of perihelion with magnitude proportional to the eccentricity. Similarly, vectors in p/q space point toward the orbital node and scale directly with inclination and indirectly with eccentricity. The top left panel of Figure 2 reproduces the test of Brown & Batygin (2019) in x/y space. By performing 100,000 iterations of randomly selecting a longitude of perihelion for each of the 17 distant (this includes three additional objects not considered in that work; see Section 2.3), clustered objects from each body's respective bias distribution, we confirm the result of Brown & Batygin (2019). In x/y space, the observed clustering (red point and circle) falls outside of ∼98.3% of the randomly selected sets of orbits. In p/q space, the result is similar (∼97.7%). It is important to point out that the significance of the observed clustering is higher with the addition of 2015 BP519, 2013 RA109, and 2013 SL102. When we remove these objects and thus consider the same 14 orbits as Brown & Batygin (2019), we replicate the results of that work (∼96.0% of iterations outside the red circle in x/y space and ∼96.5% in p/q space). In this manner, Brown & Batygin (2019) concluded that (by combining the x/y and p/q clustering) the probability that the 14 observed orbits are drawn from a uniform distribution of orientations is just ∼0.2%. However, it should be noted here that the sample of OSSOS-detected (Bannister et al. 2016) objects has been shown to be consistent with originating from a uniform distribution in independent studies (for example, Brown & Batygin 2019; Kavelaars et al. 2019).

Zoom In Zoom Out Reset image size

These results should be taken in the appropriate context given that, of the 17 objects discussed in the previous paragraph, only six have dynamical lifetimes significantly longer than the age of the solar system (Batygin et al. 2019). Thus, it is likely that several of the 17 detached objects are "interlopers" on orbits that coincidentally cluster with those of the other TNOs. The middle panel of Figure 2 plots the same test of Brown & Batygin (2019) for these six objects. In this case, the observed clustering (red point and circle) falls outside of ∼97.6% and ∼90.8% of the randomly selected sets of orbits in x/y and p/q space, respectively. While these six objects are all highly clustered in ϖ (see, for example, Figure 6 of Batygin et al. 2019), the overall statistical significance is lessened because of the small number of objects.

Next, we ask the same question in reverse. Given the cumulative biases, to what degree would the orbits of 17 objects detected from our underlying distribution of simulated test particles be clustered in p/q and x/y space? This is essentially the test used by authors investigating the most probable possible orbit for Planet Nine (e.g., Batygin et al. 2019). As our simulations include the effects of the galactic tide and stellar encounters, we argue that this analysis and comparison is an important step in assessing the strength of the Planet Nine hypothesis and the success of the most probable orbits. Furthermore, analyzing the differences between the expected clustering of the different object classes (ETNOs and IOCOs) might provide inferences into the significance of new detections.

The bottom panels of Figure 2 depict this test for our control simulations. We find that our control simulations provide a good match to the expected range of clustering that might be observed in x/y space assuming the sample of detected objects possesses orbital angles drawn from a uniform distribution (top panel). Thus, perturbations from the galactic tide and stellar encounters do not appear to induce a significant bias in our simulations. Moreover, this result implies that any bias observed in our simulations that includes Planet Nine is a result of perturbations from the hypothetical distant planet. On the other hand, the tight clustering of our control orbits in p/q space (bottom right panel of Figure 2) is a result of the region's unconstrained inclination distribution's effect on orbital pole clustering. Since the magnitude of vectors in p/q space scales with inclination, we expect objects from our low-i control runs to possess p/q values near the origin.

Figures 3 and 4 plot this same test for our simulations designed to study bulk populations of ETNOs and IOCOs, respectively. To generate these plots, we randomly draw 17 orbits from our debiased population of simulated test particles (here we define an object as "detectable" if it has q < 100 au after 4 Gyr, e.g., Sheppard et al. 2019; see further discussion in Section 3.2). While our P9b simulations consistently outperform their P9a counterparts, as expected, both sets of initial conditions are successful at generating clustered systems of distant orbits (Batygin & Brown 2016a; Batygin et al. 2019). However, it is important to point out the significant range of possible clusterings that can be generated when selecting just 17 orbits (thus, the noticeable overlap of Figures 2–4). This speaks to the problem of small number statistics that has become an intense topic of debate in the literature. Thus, while the observed clustering is inconsistent with having originated from a uniform distribution of orbits at the ≲1% level (Figure 2; Batygin & Brown 2016a; Brown & Batygin 2019), it is not possible to detect the difference between a uniform distribution and a Planet Nine–generated distribution at greater than the 1σ level. To collapse these uncertainties and achieve a 2σ distinction between our model-generated orbital clusterings and a random sample of uniformly distributed objects requires ≳100 detections. In other words, a small number of detected objects carry a greater inherent bias, which manifests in Figures 2–4 as a large range of deviations from the true underlying clustering (green point). Kaib et al. (2019) reached a similar conclusion when studying the genesis of scattered TNOs as constrained by the OSSOS data set. In Figure 6 of that work, the authors estimated that ≳60 detections of high-q objects would be required to distinguish a Planet Nine model from the uniform case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The precise degree of deviation from the origin for our clustered ETNO and IOCO orbits plotted in Figures 3 and 4 is also a function of the underlying distribution of eccentricities and inclinations. Because the magnitudes of vectors in x/y and p/q space scale with eccentricity and inclination, respectively, the range of possible detected clusterings (black points) scales likewise. This manifests most significantly in p/q space for our simulations, given the fact that we test different inclination distributions. This has the obvious effect of increasing the degree of scatter in systems that possess more high-inclination objects (right panels of Figures 3 and 4). Our plotted "detections" are bias-informed such that low-inclination objects that spend more time at the ecliptic are more likely to be detected. In spite of this correction, the effect of greater scatter with increasing σ_i persists because high-i particles still cross the ecliptic and can thus still be "detected," albeit at a lower frequency.

In Table 4, we provide relevant statistics for our various P9b simulation sets. It is important to note that these extremely distant objects' orbits are highly perturbed by the galactic tide and the giant planets over Gyr timescales. External forces cause the objects' perihelia to fluctuate by tens of au, often driving them onto orbits where they interact weakly with the giant planets and can be ejected (e.g., Sheppard et al. 2019). The extremely low fractions (∼10%–25%) of test particles and object clones that remain detectable after 4 Gyr speak to how difficult it is to study this region of the solar system with such a small sample of observed objects.

There is a clear difference between our ETNO and IOCO sets in terms of the stability of test particles and clustering in ϖ with respect to inclination. In general, longitude of perihelion clustering is most efficient in our lowest-inclination sets (also visible in Figures 3 and 4). This is a result of high-inclination dynamics randomizing objects' ϖ and Ω as the secular shepherding weakens. We direct the reader to Equation (8) of Batygin et al. (2019) for a full derivation and explanation.

In general, objects become decoupled from Planet Nine when their perihelion is driven into the planetary regime. This can be understood in terms of orbital precession induced by the Sun's artificial J₂ component (for a full derivation, see, for example, Murray & Dermott 1999):

$$\begin{eqnarray}&&{J}_{2}=\displaystyle \frac{1}{2}\sum _{i=5}^{8}\displaystyle \frac{{m}_{i}{a}_{i}^{2}}{{M}_{\odot }{R}_{\odot }^{2}}.\end{eqnarray} \tag{ 7 }$$

The longitude of perihelion and nodal precessions driven by J₂ are (e.g., Brasser et al. 2006)

$$\begin{eqnarray}&&\dot{\varpi }=\displaystyle \frac{3n}{4}{\left(\displaystyle \frac{{R}_{\odot }}{a}\right)}^{2}{J}_{2}\left(\displaystyle \frac{5{\cos }^{2}i-1}{{\left(1-{e}^{2}\right)}^{2}}\right)\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\dot{{\rm{\Omega }}}=-\displaystyle \frac{3n}{2}{\left(\displaystyle \frac{{R}_{\odot }}{a}\right)}^{2}{J}_{2}\left(\displaystyle \frac{\cos i}{{\left(1-{e}^{2}\right)}^{2}}\right).\end{eqnarray} \tag{ 9 }$$

Once the perihelion of an object is lowered into the planetary regime, stronger J₂ perturbations can result in its value of ϖ becoming decoupled from Planet Nine. This effect leads to correspondingly lower fractions of N_ϖ/N_detect for higher-inclination objects in Table 4. This process is more efficient for ETNOs (56% of the low-i batch clustered in ϖ compared with 38% of the high-i) than for IOCOs, since they are already significantly detached from the giant planet system (56% of the low-i batch clustered in ϖ compared with 51% of the high-i). In the same manner, only 12% of our ETNO high-i set and 6% of the low-i set still have q < 100 au after 4 Gyr of evolution. Consistent with previous authors, we conclude that, as a broad trend, the strength of longitude of perihelia shepherding weakens with increasing inclination and decreasing perihelion distance.

In general, we find similar fractions of objects detectable (∼10%–25%) and clustered in ϖ (∼50%–65%) and Ω (∼35%–50%) after 4 Gyr from within our various sets of clone simulations. However, we note extremely high levels of clustering for our surviving 2013 SY99 clones (78% clustered in ϖ) and relatively poor levels for the 2014 SR349 and 2004 VN112 clones. This is a result of Planet Nine's orbital shepherding having a higher efficiency for the most detached objects or those that spend the greatest fractions of time at large distances away from the giant planets. Thus it follows that, of the five ETNO clones we test (each with similar perihelia), SY99 has the largest aphelion distance (>1400 au) and also anti-aligns most efficiently with Planet Nine. Contrarily, SR349 and VN112 have aphelia close to 600 au. We demonstrate this trend in Figure 5 by plotting the aphelia of objects (from all of our simulations that include Planet Nine) that are detectable and opposed in ϖ with Planet Nine after 4 Gyr against the remainder of detectable particles. This is the reason that the most probable orbits (e.g., Batygin et al. 2019) for the hypothetical planet place its longitude of perihelion almost exactly 180° in opposition of the detached TNOs with the most extreme aphelia (namely SY99 and TG387).

Zoom In Zoom Out Reset image size

Absent an external perturber, the orbits of ETNOs and IOCOs precess exceedingly slow as the result of perturbations from the giant planets. Indeed, we would expect that, only under the influence of J₂, the longitudes of perihelia for most of the known objects in the region would traverse a full 360° only a few times over the life of the solar system (Figure 6). The most extreme objects (for example, Sedna; Brasser et al. 2006) would only be expected to make ∼one precession in 4 Gyr.

Zoom In Zoom Out Reset image size

Our simulations (Section 3.1) lead us to conclude that the possible levels of orbital clustering that might be observed from a Planet Nine–generated distribution of TNOs and that of a uniform sample overlap significantly (because of the small number of known extreme objects). Thus, while the observed degree of clustering is highly significant in spite of small number statistics and observational biases (Brown & Batygin 2019), its specific connection with a 5.0 M_⊕ planet orbiting around 500 au is not. However, there is an obvious absence of equally robust explanations in the literature (Batygin et al. 2019). While a chance or coincidental alignment might seem compelling given the minute degree of J₂-induced precession in the distant solar system (Figure 6), the probability of such an event and the corresponding dispersion timescale for a system of spontaneously clustered orbits has yet to be characterized. As such, we conclude this manuscript by commenting on the slow dispersion timescale for clustered orbits in the distant solar system.

For this phase of our study, we select the Mercury6 Hybrid integrator (Chambers 1999) and utilize a 180 day time step (∼4% of Jupiter's orbital period). In each simulation, we include the solar system's four giant planets and 50 Pluto-mass ETNOs (here defined as 300 au < a < 1000 au and 45 au < q < 80 au, loosely based on the sample of detected ETNOs and IOCOs; e.g., Table 3). The ETNOs begin with identical values of Ω and ϖ and randomly selected semimajor axes, perihelia, and mean anomalies. All simulation particles interact gravitationally with one another. In 25 simulations, we analyze ETNOs with coplanar inclinations (i = 0°). In an additional set of 25 integrations, we assign inclinations in accordance with Equation (1) (σ_i = 100, akin to our mod-i simulations of Section 3.1). Note also that these simulations do not include algorithms that account for the galactic tide or stellar encounters, as we seek to minimize unnecessary particle loss.

Figure 7 plots the average dispersal of longitudes of perihelia in our 50 simulations. At 1 Myr intervals, we calculate the angle (ϖ_test) that maximizes the ratio of objects with ϖ within ±90° of ϖ_test to those anti-aligned. Consistent with Figure 6, our higher-inclination population of TNOs remain "aligned" for longer (∼250 Myr), as they are less affected by the giant planets. It is interesting that, in both cases, the orbital alignment persists for over 1 Gyr, requiring ∼2 Gyr to reach an equilibrium value of ∼60% aligned objects. Furthermore, the deviation from a 50% equilibrium value is intriguing in and of itself. This is a direct representation of the effects of small number statistics on a binary division of angles where ϖ (or Ω for clustering in the orbital pole) is taken as a free parameter (e.g., circular statistics). This speaks to the fundamental assumption of attempts to debias the population of ETNO and IOCO detections: the precondition that, in the absence of an external perturber, the objects' orbits would be drawn from a uniform distribution. Or rather, as in our case, supposing that the detection of all KBOs with respect to ϖ or Ω is representative of the probability of detecting any extreme object with a given value of ϖ or Ω. As shown in Figure 7, when the underlying population of objects is finite, and ϖ is taken as a free parameter, the underlying distribution of angles is inherently biased and nonuniform (this idea is analogous to the observed scatter in our iterative selection of 17 random orbits from within our control set in Figure 2). This is easily exploited by selecting a subset of just 14 objects (gray line). With only 14 randomly drawn angles in the data set, it is easy to find an angle that a large fraction (equilibrium value close to 80%) appears to cluster around. This argument would obviously not hold for the classical Kuiper Belt as a whole, as it is fairly well characterized observationally. However, the underlying population from which the observed population of detached ETNOs and IOCOs originates is not the entire ETNO and IOCO constituency or bodies. Rather, it is the subset of objects (each with orbital periods of order ∼10⁴ yr) that are close to perihelion now. Such a collection of objects is finite in number and thus likely to be concentrated around some orbital elemental angles when "alignment angle" is taken as a free parameter.

Zoom In Zoom Out Reset image size