Feasibility of a Resonance-based Planet Nine Search

The primary line of evidence for the existence of "Planet Nine," a distant, massive planet in the solar system, stems from the physical confinement of KBO orbits with semimajor axis in excess of a ∼ 250 au and perihelion distance beyond Neptune. While the dominant mode of dynamical coupling between Planet Nine and KBO orbits remains a subject of active investigation, it has been suggested that mean-motion resonances (MMRs) may widely occur in the distant, eccentric confined population of the Kuiper Belt. In particular, Batygin & Brown (2016a) were the first to recognize that synthetic particles exhibiting anti-alignment with Planet Nine exhibit temporary capture into MMRs with Planet Nine persisting over hundred-Myr time intervals. Although Beust (2016), Malhotra et al. (2016) have suggested that even nonresonant bodies might survive in a detached, anti-aligned state for the lifetime of the solar system, Batygin & Morbidelli (2017) have found that, while secular effects are responsible for orbital clustering and perihelion detachment of KBOs in the confined population, MMRs are implicated in their long-term survival.

Explicit behavioral dependence on the period and phase of resonantly interacting bodies suggests that MMRs between Planet Nine and distant members of the scattered disk might encode information about the present-day semimajor axis and current position of Planet Nine, a possibility explored quantitatively by Millholland & Laughlin (2017). To this end, we note that the likelihood for a given object to occupy a specific resonance is fundamental to constraining the semimajor axis a₉ and subsequently the mean anomaly M₉ of Planet Nine. Therefore, one aim of this work is to produce a characterization of relative population sizes expected to occupy specific MMRs in the high-eccentricity case of Planet Nine.

While high-order mean-motion resonances¹ have negligible capture probability in the well-studied case of the circular restricted three-body problem, this tendency degrades in the high-eccentricity scenario of Planet Nine and, by extension, the population of KBOs sculpted by its gravity (Batygin & Brown 2016a). In other words, the typical ordering of terms found in the low-eccentricity expansion of the disturbing function (Murray & Dermott 1999) breaks down for the problem at hand, warranting a numerical evaluation of capture probabilities for specific resonances. Hence, in this work, we use numerical simulations to construct a prior distribution of MMRs among the anti-aligned KBOs. We note that this prior distribution is derived from a simplified planar model of the solar system, and therefore corresponds to a physical setup of the system where the probability of stable resonant capture is maximized (Batygin & Morbidelli 2017). Indeed, the detailed distribution of resonances depends on the specific semimajor axis and eccentricity of Planet Nine and the relative radial distance below which the dynamics become dominated by Neptune. Finally, in light of the obtained distribution, we explore the feasibility of updated resonant constraints on Planet Nine's present-day orbit and location.

The first potential resonance-based constraint is that on the semimajor axis of Planet Nine. Ideally, given the observed population of confined KBOs, it would be desirable to deduce likely period ratios with Planet Nine, thereby deriving the associated value of a₉. Using this brand of logic, Malhotra et al. (2016) reported an a₉ prediction of ∼665 au, produced by supposing that the six confined objects known at the time reside in P₉/P_KBO = N/1 and N/2 period ratios with Planet Nine. Subsequently, Millholland & Laughlin (2017) presented a more comprehensive analysis of the possibility of resonant constraints. Specifically, they constructed an a₉ distribution by assuming small integer ratio mean-motion resonances with 11 known KBOs, followed by a Monte Carlo test to confirm the significance of the highest peak in their distribution. Enticingly, the highest peak in their a₉ distribution occurs around 660 au, in apparent agreement with Malhotra et al. (2016). The second potential resonance-based constraint on the present-day, instantaneous position of the planet in its orbit requires knowledge of the resonant angles themselves. Crucially, such a constraint would be highly valuable in the observational search for Planet Nine.

This work is organized as follows. Section 2 describes our suite of semi-averaged n-body simulations. Section 3 describes the manner in which these simplified simulations capture the essence of the outer Kuiper Belt's interactions with Planet Nine, and delineates the resulting numerically derived prior distribution of MMRs. The feasibility of resonance-based determination of Planet Nine's current location is further discussed in Section 4. Concluding remarks are provided in Section 5.

By now, a considerable number of studies, employing variable levels of approximation aimed at simulating the dynamical evolution induced by Planet Nine, have been published in the literature (de la Fuente Marcos et al. 2016; Bailey et al. 2016; Batygin & Brown 2016a, 2016b; Brown & Batygin 2016; Shankman et al. 2016; Batygin & Morbidelli 2017; Becker et al. 2017; Hadden et al. 2017; Lawler et al. 2017; Millholland & Laughlin 2017). It has been found, in n-body simulations accounting for the observed inclinations of distant KBOs and the ∼20–30 degree inclination of Planet Nine, that objects tend to chaotically skip among commensurabilities. The primary aim of this work is to characterize a prior distribution of mean-motion resonances in order to ascertain the feasibility of resonance-based constraints on Planet Nine, and such transitional behavior obfuscates the classification of specific MMRs. Thus, we employ a simplified, two-dimensional model of the solar system to understand the degree of resonance-based constraints that can be made. If significant resonance-based Planet Nine constraints can be obtained from present observational data, this capability should be best reflected in this highly idealized two-dimensional model.

Within the framework of this two-dimensional model, we confine all objects to the plane and average over the Keplerian motion of the known giant planets. Accordingly, the solar system interior to 30 au is treated as a central mass with a J₂ gravitational moment having magnitude equivalent to the mean-field contribution of the canonical giant planets to the secular evolution of exterior bodies (Burns 1976; Batygin & Brown 2016a). Hence, Planet Nine is the only massive perturber in these simulations. This model omits various realistic details. Notably, modulations in the eccentricity and inclination of KBOs due to close-range interactions with Neptune, as well as dynamics induced by the mutual inclination of the KBOs with Planet Nine, are absent from our calculations. Crucially, however, due to the lack of repeated transitions between resonances induced among surviving objects, these simplified simulations allow rigorous identification of the resonances in which objects reside, and their capacity to reveal Planet Nine's parameters.

We implemented direct n-body simulations using the mercury6 integration package (Chambers 1999), employing the built-in Hybrid symplectic/Bulirsch-Stoer integrator (Press et al. 1992; Wisdom & Holman 1992), with time step chosen to be 1/8 the orbital period of Neptune. In the simulations, we evolve an initially axisymmetric disk of eccentric test particles having uniformly random angular distribution and perihelion distance and semimajor axis randomly drawn from the q ∈ [30, 50] au and a ∈ [50, 1000] au range, respectively. While the initial distribution of test particles does not reflect the complete evolution of KBOs into resonance with Planet Nine, it serves as a probe of relative strengths of resonances. For each of the eight values of e₉ tested, 6000 such test particles were randomly initialized and simulated. Particles attaining radial distances r < 10,000 au or r < 30 au were removed.

In principle, the relative strengths of resonances are not expected to vary significantly with a₉, as the relative strengths of individual terms associated with specific resonances in the usual expansion of the disturbing function only depend on the semimajor axis ratio of the interacting bodies (Murray & Dermott 1999). Therefore, Planet Nine was assigned a single, nominal semimajor axis, a₉ = 600 au. However, we note that different semimajor axes of a₉ would, in reality, subject the innermost resonances to variable levels of secular coupling with the canonical giant planets, altering the resonant widths slightly. Still, a₉ = 600 au is roughly in keeping with the semimajor axis predictions of Millholland & Laughlin (2017), Malhotra et al. (2016); thus, we choose this value of a₉. Moreover, eccentricities e₉ were tested ranging from 0 to 0.7 in increments of 0.1. The simulations in this work span 4 Gyr in approximate accordance with the solar system's lifetime.

Because this work addresses the distribution of closely spaced, high-order mean-motion resonances, which have finite width in semimajor axis, the period ratio alone is insufficient to confirm a specific resonance. Instead, we confirmed specific mean-motion resonances with Planet Nine among the surviving objects by searching for a librating resonant argument (Figure 1). The general form of such a resonant argument can be stated as φ = j₁λ + j₂λ₉ + j₃ϖ + j₄ϖ₉, where the d'Alembert relation, following from rotational symmetry, restricts the integer coefficients j to satisfy ${\sum }_{i=1}^{4}{j}_{i}=0$ (Murray & Dermott 1999). (For the two-dimensional case, we adopt the standard convention that longitude of ascending node Ω = 0, thus ϖ = ω.) Based on the values found for coefficients j₁ and j₂, identification of a critical argument informs the individual resonance in which a particle resides. Specifically, for an object in p/q resonance with Planet Nine, the resonant argument takes the form φ = qλ − pλ₉ + j₃ϖ + j₄ϖ₉.

Zoom In Zoom Out Reset image size

For the range e₉ ∈ [0, 0.7] of Planet Nine eccentricities examined, the range of longitude of perihelion explored by test particles surviving for the full length of the 4 Gyr simulations is shown in Figure 2. For all cases of nonzero perturber eccentricity, the surviving bodies occupy several specific regimes. Below a critical semimajor axis a_crit, the vast majority of surviving bodies exist in a circulating regime, in which Δϖ covers all angular values. Notably, the value of a_crit depends the semimajor axis, eccentricity, and mass of Planet Nine, and does not simply correspond to the perihelion distance of Planet Nine.

Zoom In Zoom Out Reset image size

Objects with a > a_crit exhibit orbital confinement, with the longitude of perihelion of the object relative to that of Planet Nine, Δϖ, librating about 180 degrees over the course of the simulation (i.e., in an anti-aligned configuration with respect to Planet Nine's orbit). This population is taken to be analogous to the clustered population of the observed distant Kuiper Belt objects (Batygin & Brown 2016a), and will be the focus of the following sections of this paper.

While the majority of surviving bodies found in our two-dimensional simulations occupy the two aforementioned regimes, there are some exceptions. In particular, circulating objects are also found at distances significantly beyond a_crit–specifically, within a neighborhood of the 1/2, 1/1, 3/2, and 2/1 resonances with Planet Nine. Another class of objects observed is those that migrate stochastically through a variety of semimajor axes and values of Δϖ while avoiding ejection. Although fewer than 1 in 600 bodies exhibit this behavior in our two-dimensional simulations, such migratory evolution bears resemblance to the stochastic transport between commensurabilities observed among objects in full three-dimensional simulations.

The dynamical regimes occupied by test particles show broad-ranging consistency across all nonzero Planet Nine eccentricities examined. For example, Figure 3 shows the maximum width in a and widest apsidal excursion $| \omega -180^\circ |$ exhibited among objects in identified resonances. Notably, the width in a is similar across the e₉ tested, and there appears not to be a strict relation between e₉ and confinement in Δϖ. Instead, the primary indicators of e₉ in these simulations are the location of a_crit and, in the two-dimensional case, the prevalence of circulating objects at the 1/2, 1/1, 3/2, and 2/1 resonances for lower-e₉ cases (Figure 2).

Zoom In Zoom Out Reset image size

3.1. Numerically Derived Period Ratio Distribution

All confined objects in the simulation are in confirmed mean-motion resonance with Planet Nine. Moreover, the identified resonant angles persist over the full 4 Gyr simulation. This stable behavior can be attributed to the lack of interference from close encounters with Neptune occurring in unaveraged three-dimensional simulations. Figure 4 shows the number of bodies found in individual resonances across the seven simulations with e₉ ranging from 0.1 to 0.7. In all simulations with e₉ > 0, the individual resonances with the greatest number of occupied objects were the 1/1 and exterior 2/1 resonances. However, as the next section will discuss, most objects do not occupy these particular commensurabilities.

Zoom In Zoom Out Reset image size

As Planet Nine's semimajor axis defines its period ratio with resonant KBOs, it is natural to attempt to predict a₉ through commensurabilities. Such predictions have been produced by Malhotra et al. (2016), Millholland & Laughlin (2017), in rough agreement at a₉ ∼ 665 and a₉ ∼ 654 au respectively. In particular, Malhotra et al. (2016) postulated that a₉ ∼ 665 au, placing Sedna at interior 3/2 resonance and five other objects at N/1 and N/2 interior resonances with Planet Nine. However, the prior distribution of MMRs derived in this work appears to suggest that the majority of objects do not have period ratios with Planet Nine of the form N/1 or N/2 (Figure 4). Furthermore, the probability that any six independently chosen objects all have N/1 or N/2 period ratios is less than 0.05 (Figure 5). We therefore conclude the assumption that a predominance of distant KBOs are likely to reside in N/1 and N/2 resonances with Planet Nine is not supported by these numerical simulations.

Zoom In Zoom Out Reset image size

We now turn to our numerically derived resonance distribution from our two-dimensional simulations and consider what information can be gathered about Planet Nine's whereabouts, following the statistical approach developed by Millholland & Laughlin (2017). Each iteration of their method proceeds as follows: First, they draw a sample m_i from a truncated Gaussian mass distribution for Planet Nine, with mean 10 M_⊕ and bounds 5–20 M_⊕. Then, they randomly select one of the known distant KBO semimajor axes, together with one orbital period ratio from a chosen period ratio distribution. The implied Planet Nine semimajor axis is then calculated. If the implied semimajor axis lies in the range [200 au + 30m_i/M_⊕, 600 au + 20m_i/M_⊕] (derived based on Brown & Batygin 2016), a Gaussian centered at the implied semimajor axis, with σ equal to the approximate resonance width, is added to the a₉ distribution being constructed. These steps are repeated until the resulting a₉ distribution is converged upon. In this way, Millholland & Laughlin (2017) have provided a rigorous means of estimating a₉ by inferring mean-motion resonances with the observed Kuiper Belt objects. However, the period ratio distribution from which they drew their sample was simplified; they considered period ratios in the Farey sequence F_N, which includes all fractions with denominators up to N (for example, ${F}_{4}=\{\tfrac{1}{2},\tfrac{1}{3},\tfrac{1}{4},\tfrac{2}{3},\tfrac{3}{4}\}$). Critically, Millholland & Laughlin (2017) made the assumption that period ratios in F_N are equally occupied.

In contrast, Figure 6 illustrates the difference between the a₉ distribution reported by Millholland & Laughlin (2017) using the Farey sequence F₅ distribution, versus the distribution obtained using the identical methodology, as well as the identical set of observed objects with the same values as those used by Millholland & Laughlin (2017), but employing the period ratio distribution obtained in this work. We note that, because Planet Nine's eccentricity is poorly constrained, and because the period ratio distributions for varied e₉ resemble each other across the board, we have considered the sum distribution of period ratios found for all cases e₉ > 0. Due to the close spacing of many occupied high-order resonances, the prominent peak of the distribution of Millholland & Laughlin (2017) is replaced in our a₉ distribution by a broad plateau. We note that reduction of peak prominence also occurs when the resonance distribution is taken to be an equally weighted Farey sequence of higher order, although such treatment neglects the relative population occupying each resonance.

Zoom In Zoom Out Reset image size

Indeed, our resulting plateau-shaped distribution for a₉ demonstrates that when the prevalence of high-order MMRs is considered, hopes of a resonance-based constraint on a₉ all but vanish, at least given current observations. Furthermore, compared to the simulations of a simplified Planet Nine system described in this work, additional behavior arises in fully inclined simulations which include the canonical giant planets (especially Neptune). Because the behavior exhibited by eccentric test particles in these more comprehensive simulations includes repeated transitions between resonances, fully rigorous determination of the distribution of objects among mean-motion resonances is challenged. Moreover, in contrast to the planar simulations in which all surviving anti-aligned objects are in resonance, anti-aligned objects in full three-dimensional simulations are not all in identifiable resonances with Planet Nine. Therefore, the prevalence of low-order resonances found in the two-dimensional case is greater than what can be realistically expected. Moreover, a running average of period ratios with Planet Nine among objects in full simulations shows no obvious prevalence of objects occupying particular resonances (Figure 7). Indeed, this lack of obvious resonant populations suggests that N/1 and N/2 period ratios such as those endorsed by Malhotra et al. (2016) are substantially less common in the more realistic three-dimensional case than in planar simulations (Figure 5). Without constraints on the period ratios between observed objects and Planet Nine, constraints on Planet Nine's present-day location along its orbit remain elusive.

Zoom In Zoom Out Reset image size

Using two- and three-dimensional direct n-body simulations, we have highlighted several features inherent to the resonant dynamics of distant KBOs. Overall, we have illustrated that the resonance-based search for Planet Nine is frustrated by the increased strength of high-order resonances that follows from the high eccentricity of Planet Nine and the test particles. Existing work aiming to predict a₉ by considering possible period ratios between Planet Nine and observed objects (Malhotra et al. 2016; Millholland & Laughlin 2017) unfortunately neglects the prevalence of high-order resonances and chaotic transport–traits of the system, which, in practice, challenge characterization of the specific resonances occupied.

A constraint on a₉, to a degree allowing inference of the specific resonances in which observed objects reside, appears crucial for constraining the true anomaly of Planet Nine through resonant means. The degeneracy between mass and semimajor axis in simulations further complicates such efforts. Moreover, because knowledge of the specific period ratios between objects and Planet Nine are necessary to rule out Planet Nine's instantaneous position from specific regions of its orbit, resonance-based constraints on the current orbital position of Planet Nine do not appear feasible at this time.

The prevalence of high-order resonances in our simulations, brought about by the high eccentricity of the system, serves as a reminder that care should be taken before assuming the dynamics of circular orbits will hold for an eccentric case, such as Planet Nine's interaction with the confined population. For example, the lack of orbital angle clustering found by Lawler et al. (2017) in n-body simulations can be explained in part as resulting from their condition that e₉ = 0. It has already been established (Batygin & Brown 2016a; Brown & Batygin 2016) that Planet Nine's orbit must be eccentric to produce the confinement of distant, eccentric Kuiper Belt objects. Although it follows from symmetry, the results of our direct n-body simulations provide additional illustration of this fact (Figure 2). Because confinement occurs at the lowest nonzero Planet Nine eccentricity tested, e₉ = 0.1, the lowest eccentricity necessary to produce confinement in the planar eccentric three-body system is a possible subject of interest for future work.

In summary, the expected phase-protected mean-motion resonances between observed distant KBOs and Planet Nine offer a tantalizing connection to the semimajor axis and current position of Planet Nine. Due to the high eccentricities involved, the underlying resonant dynamics are fundamentally different from the circular case, and high-order resonances with Planet Nine appear to dominate in the anti-aligned population. Considering this and taking the updated MMR distribution into account, the obtainable constraints on a₉ appear far less useful than suggested by Malhotra et al. (2016), Millholland & Laughlin (2017). A resonance-based constraint on the mean anomaly M₉ appears to require a constraint on a₉, so that correct resonant angles of observed objects can be deduced and specific regions of the orbit excluded. Thus, in addition to its unprecedented nature among planets in the solar system, the large eccentricity of Planet Nine currently postpones precise prediction of its position.

We wish to thank Fred Adams and Sarah Millholland for useful discussions. In addition, we thank Caltech High Performance Computing for managing the Beowulf cluster Fram and providing technical assistance, as well as the anonymous reviewer, for a thorough and insightful review that has led to a substantial improvement of the manuscript.