Generation of Low-inclination, Neptune-crossing Trans-Neptunian Objects by Planet Nine

To construct the inter-Neptunian perihelion distribution, we followed previous studies (e.g., Becker et al. 2018; Hadden et al. 2018; Li et al. 2018; Batygin et al. 2019; Brown & Batygin 2021 and references therein) and examined the orbital footprints generated by particles satisfying the same orbital cuts as those adopted in Figure 1, within the final Gyr of the integration. These footprints were recorded at 1 Myr intervals, far exceeding the typical Lyapunov time of the particles. ⁵ Due to chaotic mixing, any two footprints, even if sequentially produced by the same particle, are effectively uncorrelated.

As a null hypothesis, we considered the P9-free cluster_2 simulation of Nesvorný et al. (2023). Despite having been conducted with a different integrator, this simulation includes all of the physical effects described in Section 2, with nearly identical implementation. Moreover, the Nesvorný et al. (2023) model, previously validated against the distribution of high-q TNOs, represents the current benchmark for the post-nebular evolution of the solar system. Given the larger particle count in this simulation compared to our P9 model, we sampled the final Gyr at 20 Myr intervals, yielding a similar number of footprints (between 2000 and 3000) satisfying a > 100 au, q < 30 au, i < 40°, and a sufficiently large sample to construct smooth histograms for both scenarios. We further verified that the histogram shapes remained consistent over time, indicating that the flux of Neptune-crossing objects had attained steady state by the final Gyr of the integration.

The left and right panels of Figure 3 compare the raw (unbiased) i < 40° semimajor axis-perihelion distributions from simulations with and without P9, respectively. While both models yield semimajor axis distributions that diminish with increasing a at long periods, the perihelion distributions are markedly different. The P9-free run shows a rapid decline in perihelion distribution with decreasing q, as Neptune's orbit forms a veritable dynamical barrier. In contrast, the simulation that includes P9 results in a relatively flat q distribution outside ∼16 au, with a notable dip at q ∼ 20 au (this feature can be attributed to strong gravitational interactions with Uranus, leading to a mild depletion in the perihelion distance distribution at this specific range).

Zoom In Zoom Out Reset image size

Qualitatively, these differing q distributions are tenable. In the P9-free scenario, objects with a ≲ 1000 au are too close to be significantly influenced by Galactic effects, leaving chaotic diffusion—which tends to preserve q ∼ a_N—as the primary driver of orbital evolution. On the other hand, P9-induced dynamics can continuously modulate the perihelion, even if the orbit dips well below Neptune's semimajor axis. Although the existing observational data indeed reveal a perihelion distribution that is relatively flat (Figure 1), we cannot compare the data to the modeled distributions without accounting for observational bias.

Well-known observational biases exist against detecting orbits with high inclinations, as well as objects at large heliocentric distances. The former effect arises because observational surveys—such as Pan-STARRS-1 and -2, which account for the discovery of a significant fraction of the observed objects—are often performed at low ecliptic latitudes, where high-inclination objects spend less time. Low-i objects are thus overrepresented in the catalogs. Fortunately, for the problem at hand, this bias is largely inconsequential because the numerical models do not show substantial difference in the perihelion distribution of q < 30 au objects as a function of inclination, in the i ∼ 0°–40° range. ⁶ Thus, by restricting our analysis to TNOs with inclinations lower than 40°, we mitigate the principal source of inclination bias. Addressing the heliocentric distance bias, on the other hand, requires a more nuanced treatment, which is developed below.

Bias correction informed by discovery distance. To account for the bias toward detecting objects with lower perihelia, we begin by examining the heliocentric distance at which each of the 17 known objects depicted in Figure 1 were discovered. To leading order, at the moment of discovery, each object serves as an unbiased probe of the entire perihelion distribution, up to its discovery distance. This notion effectively nullifies biases related to discovery distance or object size, as the increasing brightness of an object with diminishing heliocentric distance becomes irrelevant.

As a concrete example, consider a large sample of objects, all detected at 30 au. To first approximation, the perihelia of this collection of bodies constitute an unbiased probe of the q distribution inside of 30 au, regardless of the brightness of the object at the time of discovery ⁷ or the limiting magnitude of the survey. To advance beyond this estimate, one important correction must be made: objects with different q and a spend different fractions of their orbital period in the vicinity of 30 au. It is, however, straightforward to apply this geometric correction and weight the distribution accordingly (Brown 2001; Morbidelli et al. 2004).

In practice, we do not have a large aggregate of objects that were all discovered at the same distance but rather a modest collection of objects, all discovered at different distances. Nevertheless, each of these detections amounts to an unbiased probe of the perihelion distribution interior to their discovery distance, as described above. We can thus compare each discovery made at a particular distance to our modeled perihelion distributions for all objects interior to this discovery distance, in the presence and in absence of P9.

Taking into account the geometric bias of the orbit, we use the simulation data to construct a PDF of the perihelion distribution for each discovery distance and then compute where each observed object falls within its individual cumulative distribution function (CDF). We label the resulting quantity ${\xi }_{j}={\mathrm{CDF}}_{{r}_{j}}({q}_{j})$. If the model perfectly matches the observations, ξ_j should be uniformly distributed between 0 and 1. Any deviation from this uniformity serves as an unbiased statistical measure of the congruence between the model and the observational data. The top panel of Figure 4 shows the distribution of ξ for the P9 and P9-free models. The Kolmogorov–Smirnov test for uniformity yields starkly distinct values, yielding p = 0.41 (P9) and p = 0.0034 (P9-free), thus significantly favoring the model that includes P9.

Zoom In Zoom Out Reset image size

To further quantify the statistical discrepancy, we define the statistic $\zeta ={{\rm{\Pi }}}_{j}^{{N}_{\mathrm{obj}}}\mathrm{log}({\xi }_{j})$. If the distribution of ξ is uniform, then in the limit of ${N}_{\mathrm{obj}}\to \inf$, the distribution of expected values of ζ approaches a Gaussian. Though our sample size (N_obj = 17) is not sufficiently large for this result to hold exactly, this approach still allows for a rigorous measure of uniformity of ξ. We begin by constructing an expected distribution of ζ by computing ${{\rm{\Pi }}}_{j}^{17}\mathrm{log}({\rm{U}}(0,1))$ one million times, thereby generating a smooth histogram. The resulting curve is shown on the bottom panel of Figure 4. We then compute the values of ζ for both simulations. Intriguingly, the P9 simulation's statistic (ζ = −7.9) aligns closely with the mean of this distribution (〈ζ〉 = −7.2), while the P9-free model's value (ζ = −16.5) deviates significantly, falling approximately 5 standard deviations (σ = 1.8) away from the peak. This comparison provides a quantitative assessment of the models, indicating a much higher likelihood of the P9 model given the current observational data.

As a check on the self-consistency of our statistical method, we conducted a validation exercise using simulated observations derived from a precisely defined synthetic distribution. This involved generating 100,000 orbits from on a synthetic distribution of TNOs characterized by a boxcar perihelion distribution within the 15–30 au interval, a Gaussian distribution of inclinations with a dispersion of 15°, and a semimajor axis distribution proportional ⁸ to ${dN}/{da}\propto {a}^{-3/2}$. Subsequently, we simulated the observation of 20 objects using the OSSOS survey simulator (Lawler et al. 2018) and subjected these observations to our statistical analysis.

The outcomes of this test confirmed the method's robustness: applying our metric to the synthetic population resulted in a KS p-value of 0.72 for the uniformity of ξ, and the value of ζ (−9.8) was only 0.6σ away from the mean value (−8.7) of the expected distribution. For completeness, we also acknowledge the potential for ecliptic longitude-dependent bias in the perihelion distribution due to the survey footprint of individual surveys like Pan-STARRS. Nevertheless, our simulations do not show any substantial dependence of the perihelion distribution on the longitude, meaning that this form of survey bias is almost certainly secondary and is unlikely to meaningfully impact our conclusions.

Magnitude-limited bias correction. A distinct method to address observational bias in simulation data is the magnitude-limited correction approach. This technique simulates the detectability of objects within a generated orbital distribution for a survey with a specific limiting magnitude (see, e.g., Morbidelli et al. 2004). While our current data set comprises objects discovered across various surveys ⁹ —which precludes a straightforward application of this method—understanding its implications remains beneficial. Specifically, this approach serves as a predictive tool for future uniform surveys, such as those planned for the Vera Rubin Observatory (VRO).

In choosing a limiting magnitude, we set ${V}_{\mathrm{lim}}=24$, coinciding with the anticipated capabilities of VRO. The size distribution assumes a power-law exponent of η = 2/3 drawing on the results of Fraser et al. (2014) for objects with absolute magnitude range of H ∼ 6–9 though we find that the results are only weakly dependent on this choice. By accounting for the fraction of the orbit visible and the time spent traversing it as above, we generate biased distributions of a and q for both models and show them as histograms on Figure 3.

The biased semimajor axis distributions for both P9 and P9-free scenarios show a similar decay with increasing a, rendering them more akin to each other than in their unprocessed form. However, the perihelion distributions reveal a persistent, notable disparity: even after accounting for observational biases, the perihelion distribution in the P9-inclusive model retains a relatively flat distribution beyond approximately 16 au. As already discussed above, this characteristic bears considerable resemblance to the distribution of the actual data presented in Figure 1. In stark contrast, the P9-free model continues to exhibit a pronounced peak around 30 au. This analysis indicates that the perihelion distribution of low-inclination TNOs is very likely to remain an important indicator for the existence of P9 in forthcoming data sets.

As important as the comparison with existing observations, the results presented herein offer a set of readily falsifiable predictions, with near-term prospects for resolution. To this end, we note that any comparison with the current data, even when biases are accounted for, is inherently imperfect, and a more uniformly acquired set of observations would provide a superior basis for testing our model. Fortunately, with the expected commencement of operations by the VRO, the orbital distribution of the class of objects considered here (Figure 1) will come into much sharper focus, and the a − q orbital distribution depicted in Figure 3 will be tested directly.

Another observational handle is provided by examining the absolute ratio of Neptune-crossing objects to those with q > 30 au. Though the number of particles shown in the two panels of Figure 3 is approximately equal, the number of total orbital footprints that were used to generate the P9-free panel significantly exceeds that of the P9-inclusive panel. This discrepancy is a result of the more efficient injection of objects into the q < 30 au region in the presence of P9. More quantitatively, the ratio of Neptune-crossing objects with inclination i < 40° and semimajor axis between 100 and 1000 au to those with q > 30 au is ∼3% in the P9 scenario, compared to only ∼0.5% in the P9-free case. While the current observational census does not provide a rigorous means to quantify this value, the advent of a comprehensive survey conducted by VRO will offer a more definitive opportunity to evaluate this prediction.

Finally, the predicted inclination distribution provides an avenue of inquiry. Although we have not delved into inclination biases in detail here, it is noteworthy that the inclination distributions produced by P9 and P9-free simulations are strikingly different. Specifically, the distribution in the presence of P9 shows a steep rise with i, below 30°. Meanwhile, the P9-free model exhibits a considerably flatter dispersion (Figure 5). This prediction will be put on solid statistical footing with forthcoming results from VRO.

Zoom In Zoom Out Reset image size

As concluding points, we note that while P9 explains the anomalous structure of the outer solar system in a unified framework, several alternative theories have been proposed to account for individual aspects of the P9 hypothesis. We briefly review these theories here and discuss how our work fits into this broader context. First, as already mentioned above, Nesvorný et al. (2023) have shown how cluster-induced dynamics can generate the perihelion-detached population of TNOs, indicating that the q-broadening process may have been primordial and that P9 is not strictly necessary to explain this observation. ¹⁰ In a parallel vein, Huang et al. (2022) have proposed the possibility that a few-Earth-mass rogue trans-Neptunian planet could have influenced the outer solar system's structure for hundreds of Myr, before being removed by some process. ¹¹

With respect to orbital clustering, Shankman et al. (2017), Bernardinelli et al. (2020), and Napier et al. (2021) have shown that individual surveys, which have examined limited areas of the sky, generally struggle to overcome their inherent observational biases sufficiently to rigorously determine the presence or absence of orbital alignment. This limitation has led some authors to interpret the observed orbital alignment as being illusory. Despite these challenges, it is important to recognize that, even with the strong biases of the Dark Energy Survey, Bernardinelli et al. (2020) reported a ∼2σ level of significance in the clustering of the longitude of ascending nodes. Additionally, a comprehensive observability analysis of all available data indicates that, after accounting for observational biases, distant Kuiper Belt Objects are jointly clustered in Runge–Lenz (eccentricity) and angular momentum vectors with a significance level of approximately 99.6% (Brown 2017; Brown & Batygin 2019, 2021). Finally, the observed anticorrelation between the rate of orbital diffusion and clustering within the data (as discussed in Batygin et al. 2019) is unlikely to be attributable to observational bias alone. In a separate development, Huang & Gladman (2024) have recently proposed that the alignment of three specific TNOs—Sedna, 2012 VP₁₁₃, and Leleakuhonua—is real but is not created by P9. Instead, in their picture, these objects' orbits could have been shaped by an early event in the solar system's history and have since precessed just enough to realign in the present epoch.

For the high-inclination population, Kaib et al. (2019) have showed that the flux of retrograde Centaurs, as inferred from the OSSOS survey, is too large to be accounted for by a solar system model that excludes P9. Their calculations further showed that the presence of P9 could reconcile this discrepancy although the adopted P9 parameters led to a median inclination of simulated detections that is a few degrees higher than the observed value. Still, as an alternative explanation, Kaib et al. (2019) have also proposed that past migration of the Sun through the galaxy could have enriched the Oort cloud, thereby enhancing the flux of retrograde Centaurs.

In contrast to all of the above, the Neptune-crossing objects we have focused on in this work are distinct in a crucial way: due to their low inclinations and perihelia, these objects experience rapid orbital chaos and have short dynamical lifetimes (a simple rebound simulation illustrates that, in the absence of P9, objects shown in Figure 1 have a dynamical lifetime on the order of ∼100 Myr; Rein & Tamayo 2015). This implies that the dynamical process responsible for their current orbits is ongoing, not a relic of the distant past. Accordingly, any alternative to P9 that aims to explain this population must invoke active perturbations beyond those accounted for in calculations of Nesvorný et al. (2023). One such possibility is modified Newtonian dynamics (MOND), and recent proposals by Brown & Mathur (2023) and Migaszewski (2023) have suggested that MOND might explain some phenomena attributed to P9. However, this hypothesis faces significant challenges, as Vokrouhlický et al. (2024) have shown that MOND significantly disrupts the observed specific energy distribution of long-period comets and have further illustrated that certain variants of MOND fail to accurately account for the dynamics of detached objects. Additionally, Banik et al. (2024) have used wide binary data from the Gaia DR3 data set to demonstrate that Newtonian gravity is very strongly favored over MOND on outer solar system scales. Equivalent conclusions were reached by Fienga et al. (2016) and Fienga & Minazzoli (2024) from the vantage point of planetary ephemerides.

Another class of alternative theories involves the self-gravitational dynamics. Madigan & McCourt (2016) were the first to propose the inclination instability as a mechanism for shaping the present-day architecture of the outer solar system. Follow-up studies have refined this picture, arguing that this instability could naturally manifest within a disk of initially planar but highly eccentric minor bodies, totaling around 10 Earth masses (see Zderic & Madigan 2023 and the references therein). Yet, as shown in the GPU-accelerated simulations of Das & Batygin (2023), even with such a massive disk, the inclination instability is fully suppressed if Neptune scattering is modeled self-consistently. In an unrelated effort, Sefilian & Touma (2019) explored the idea of a similarly massive, mildly lopsided disk of planetesimals extending to about 700 au as a potential driver of P9-like dynamics. Nevertheless, a suitable explanation for the origin of such a shepherding disk—and more importantly—its capacity to remain coherent on multi-Gyr timescales remains elusive.

In a study closely related to our work, Nesvorný et al. (2023) showed that the same interplay between giant-planet scattering and cluster-driven evolution that would have raised the perihelia of distant TNOs could also trap as much as ∼3 Earth masses of material within the inner Oort Cloud—potentially facilitating nontrivial orbital evolution. In a recent paper (Batygin & Nesvorný 2024), we analyzed the emergent phase-averaged dynamics of this scenario and found that the physical picture is qualitatively identical to that of vZLK cycles. Our calculations also indicate that unless the mass of the inner Oort Cloud is taken to be unreasonably large (i.e., tens to hundreds of Earth masses), the characteristic timescale of these cycles would far exceed the age of the Sun. Therefore, at present, P9 remains the only plausible explanation for the observed distribution of long-period Neptune crossers.

In summary, this work has introduced a new line of evidence supporting the P9 hypothesis. Excitingly, the dynamics described here, along with all other lines of evidence for P9, will soon face a rigorous test with the operational commencement of the VRO. This upcoming phase of exploration promises to provide critical insights into the mysteries of our solar system's outer reaches.