The Generation of the Distant Kuiper Belt by Planet Nine from an Initially Broad Perihelion Distribution

The continued unveiling of the outer solar system by astronomical surveys has revealed a puzzling orbital clustering among the distant Kuiper Belt objects (KBOs). Trujillo & Sheppard (2014) were the first to note a peculiar grouping in the argument of perihelion of long-period KBOs. Subsequently, Batygin & Brown (2016a) found that the orbits of KBOs with semimajor axes a ≳ 250 au and perihelion distances q > 30 au are in fact anomalously aligned in physical space, i.e., these objects exhibit statistically significant clustering in longitude of perihelion as well as ascending node. To explain this persistent orbital grouping, Batygin & Brown (2016a) proposed the existence of a Neptune-like perturber—"Planet Nine"—that is apsidally anti-aligned with the aforementioned KBOs and generates their orbital confinement through a complex interplay of resonant and secular effects (Batygin & Morbidelli 2017; Becker et al. 2017).

While this orbital clustering constitutes the most visually striking line of evidence for the existence of Planet Nine (Brown & Batygin 2016), it is by no means the only one. In particular, the detachment of high perihelion (q ≳ 36 au) KBOs from Neptune requires an extra gravitational influence, which a distant perturber can successfully provide. In addition, Planet Nine naturally explains the dynamics of highly inclined large semimajor axis centaurs (Gomes et al. 2015; Batygin & Brown 2016a), as well as retrograde KBOs with a < 100 au, such as Drac and Niku (Batygin & Brown 2016b). Lastly, the presence of an inclined distant planet in our solar system and the associated secular torque can explain the otherwise mysterious six-degree obliquity of the Sun (Bailey et al. 2016; Lai 2016; Gomes et al. 2017).

Given the inferred distant and eccentric orbit of Planet Nine, a natural question that arises concerns its formation. To this end, three main mechanisms have been considered in the literature. In the first scenario, Planet Nine is captured into its current orbit, either from another solar system, or as a free-floating planet (Li & Adams 2016; Mustill et al. 2016). In the second, the formation of Planet Nine is a two-step process. Particularly, originating in the inner solar system, Planet Nine is envisioned to have been scattered out to a wide orbit by the giant planets (possibly during the transient period of dynamical instability prescribed by the Nice model: Tsiganis et al. 2005; Nesvorný 2011; Batygin et al. 2012). As the hypothesized current orbit of Planet Nine has a high perihelion (q₉ > 250 au), gravitational influence of the solar system's birth cluster (Morbidelli & Levison 2004; Adams 2010) or dynamical friction from an extended massive planetesimal disk (Eriksson et al. 2018) is then invoked to perturb the orbit and raise Planet Nine's perihelion distance. The third mechanism is in situ formation, that is, the formation of Planet Nine in its current position (Bromley & Kenyon 2016).

Given how distant the orbit of Planet Nine is hypothesized to be, the third scenario would require efficient conglomeration of solids at a large stellocentric distance as well as a disk that extends to hundreds of au. The probability of the other two mechanisms was evaluated by Li & Adams (2016) and Eriksson et al. (2018); the authors show that in both cases, special conditions are required to reproduce and retain the orbit of Planet Nine. In light of these results, the formation mechanism of Planet Nine remains unclear and of considerable interest.

Unlike the in situ formation scenario, we note that both the capture and outward scattering mechanisms rely on the occurrence of stellar encounters. Importantly, the circumstances of these stellar interactions certainly affect the initial Kuiper Belt structure. In other words, the same perturbations that would have lifted Planet Nine's perihelion would have also emplaced a population of distant KBOs onto high-q orbits. As of yet, however, a systematic study of the evolution of the Kuiper Belt from different initial distributions has not been done.

In this paper, we investigate this question in detail. The paper is structured as follows. In Section 2, we introduce our simplified numerical model and outline the key differences in the evolution of the narrow and broad-q Kuiper Belt populations. In particular, we focus on the dynamics of the objects in the apsidally "aligned" and "anti-aligned" classes (with respect to Planet Nine's longitude of perihelion). We then contextualize our results within the framework of more realistic N-body simulations. We summarize and discuss the implications of our results in Section 3, in particular focusing on the extent to which our model can aid in interpreting current observations.

To study the dependence of the Planet Nine hypothesis on the initial conditions of the Kuiper Belt, we explore the evolution of two initially distinct Kuiper Belt distributions. Our choice of distributions is motivated by the proposed mechanisms for Planet Nine's formation. The first population, the "narrow" Kuiper Belt (q ∈ [30, 36]), corresponds to the conventional definition of the scattering disk (Gladman et al. 2008). The second "broad" population (q ∈ [30, 300]) is associated with mechanisms involving stellar encounters, which could lift the perihelion distance of both Planet Nine and the KBOs. In the scattering or capture scenarios, then, the Kuiper Belt q-distribution may range from Neptune's semimajor axis all the way to Planet Nine's perihelion distance and further.

With these considerations in mind, we begin our study of the evolution of these two populations by employing an idealized numerical model. We carry out two sets of N-body simulations that span 4 Gyr of integration time and utilize the mercury6 software package (Chambers 1999), employing the hybrid Wisdom–Holman/Bulirsch–Stoer algorithm.

As a first step, we consider a simplified planar model of the solar system, which includes Neptune and Planet Nine as active bodies, but replaces Jupiter, Saturn, and Uranus with a solar J₂ moment as in Batygin & Brown (2016a). That is, we account for the gravitational potential of the three planets by adding a nonzero quadrupolar field to the Sun and extending its radius to Uranus' semimajor axis. In this semi-averaged model, we introduce a 10 m_⊕ Planet Nine with semimajor axis a₉ = 700 au, eccentricity e₉ = 0.6, and inclination i₉ = 0°, and set the time step to 16 years (one-tenth of Neptune's orbital period).

While capture of Planet Nine by the solar system is unlikely to impulsively perturb the Kuiper Belt in any appreciable manner, the same is not necessarily true if Planet Nine attained its initial long-period orbit by scattering off of Jupiter or Saturn. Simulations of this process, however (Batygin et al. 2012), demonstrate that even in this violent scenario, the dynamical footprint of the scattering upon the Kuiper Belt is modest in comparison with that entailed by the Nice model instability (Levison et al. 2008).

The first set of numerical simulations contains a synthetic distant Kuiper Belt in which the perihelion distance of the test particles is uniform in the range q ∈ [30, 36] au (the "narrow" Kuiper Belt), consistent with the expectations of the Nice model in the absence of external perturbations. Conversely, in our second set of simulations, the perihelion distance of the test particles uniformly ranges between q ∈ [30, 300] au (the "broad" Kuiper Belt). The initial distributions of the other orbital parameters of the objects are the same in both experiments: the semimajor axes are drawn uniformly³ from a ∈ [150, 1000] au, inclinations follow a half-normal distribution centered at i = 0° with σ_i = 5°, and the argument of perihelion (ω), longitude of ascending node (Ω), and mean anomaly (M) are uniformly drawn from the full [0, 360) degree range. Each set of simulations contains about 10000 massless particles.

Objects surviving the full 4 Gyr integrations make up just a small fraction of the initially launched population, as can been seen in Figure 1. Rather than analyzing the q-distribution of the entire aggregate of test particles, here we focus on the behavior of specific dynamical classes of KBOs. Particularly, the two dynamical populations of interest consist of the objects that experience longitude of perihelion librations with respect to Planet Nine's longitude of perihelion, notated as Δϖ: the "anti-aligned" and the "aligned" populations. The objects anti-aligned to Planet Nine are identified as those which librate around Δϖ = 180°, while those that are aligned librate around Δϖ = 0°.

Zoom In Zoom Out Reset image size

In our analysis, we then track the evolution of the initially uniform perihelion distance distribution of the anti-aligned and aligned populations. Notably, we find that the final q-distribution of these objects in the narrow and broad Kuiper Belt populations is strikingly different. The top and bottom panels of Figure 1 show the final distribution of the narrow and broad populations, respectively. Evidently, very few aligned objects initialized within the narrow population survive for 4 Gyr (the surviving aligned population makes up only 0.14% of the original ∼10000 objects). On the other hand, a clear bimodal structure—with the higher q-peak corresponding to the aligned population (3.70% surviving), and the lower q-peak corresponding to the anti-aligned population (5.84% surviving)—emerges in the numerical experiments characterized by the initially broad q-distribution of KBOs. Of course, not all of these objects are observable; applying a d⁻⁴ detection bias, we compute that the median final perihelion distance of the detectable objects in the narrow distribution (top panel of Figure 1) is 38.13 au for the anti-aligned KBOs and 33.07 au for the aligned KBOs, and in the broad distribution (bottom panel), 44.16 au for the anti-aligned KBOs, and 100.29 au for the aligned KBOs.

In the following sections, we consider the difference between the final perihelion distance distributions of the narrow and broad populations in greater detail.

2.1. The Anti-aligned Population

As shown by Batygin & Morbidelli (2017), the observed dynamics of the anti-aligned population stems from a complex interplay between mean-motion resonances and secular interactions. That is, in face of orbit-crossing, the anti-aligned objects glean long-term stability by chaotically hopping (Becker et al. 2017) between phase-protected mean-motion commensurabilities (Malhotra et al. 2016; Millholland & Laughlin 2017). At the same time, the specific (corotation) nature of these resonances implies that the anti-alignment of the orbits itself is maintained via a phase-averaged flavor of dynamical coupling. Although this behavior is well understood, the effect of the initial q-distribution on this population is not immediately evident.

Interestingly, the distribution of the anti-aligned objects (light blue, Figure 1) is similar in both the narrow and broad populations. Yet, we see that these two distributions are not identical: the broad population is able to retain a greater number of anti-aligned objects, as many of these particles start out and remain detached from Neptune and thus experience minimal scattering interactions with the giant planets. Even so, the final q-distributions of both populations are dominated by objects with middling perihelion distances, with most of the objects possessing q_f < 100 au.

It is intriguing to note that anti-aligned objects in the broad Kuiper Belt population with initial q in excess of ∼100 au do not survive the full integration time. This is because such objects reside on Planet Nine-induced secular trajectories that experience high-amplitude perihelion distance oscillations. Once their maximum q is attained, the objects continue to approximately follow the level curves of the averaged Hamiltonian (computed as in Batygin & Morbidelli 2017) and are driven into Neptune-crossing orbits. At this point, these objects typically become centaurs, and as centaurs have short dynamical lifetimes (Tiscareno & Malhotra 2003; Horner et al. 2004), we allow them to be absorbed by the central body rather than resolving their interactions with the canonical giant planets.

We illustrate this effect in Figure 2: the light blue trajectories shown in the q–Δϖ phase plane attain high-q values as they librate in Δϖ. The initial (Δϖ₀, q₀) values of these objects are denoted with blue dots, and the red stars indicate the moment when these objects cross Uranus' orbit. These metastable objects survive less than a single libration in Δϖ, as they are driven to low q before returning to their initial locations in the phase plane (see the lower region of Figure 2). In fact, any real objects currently residing at any point on these light blue secular trajectories are expected to become unstable due to this effect.

Zoom In Zoom Out Reset image size

In contrast to this metastable evolution, the dark blue objects nested within the lighter blue loops are the trajectories of stable anti-aligned objects. A schematic of the geometry of the longitude of perihelion libration is shown in Figure 3. These particles experience only mild libration in both q and Δϖ, and are thus able to remain in a long-term stable configuration protected from Neptune-scattering. Accordingly, we find that after a 4 Gyr integration, all anti-aligned objects that attain high-q values in their trajectories are eventually ejected. This qualitatively explains why the perihelion distribution of stable anti-aligned objects is limited to q ≲ 100 au in both simulations of the broad and the narrow Kuiper Belt populations.

Zoom In Zoom Out Reset image size

2.2. The Aligned Population

In addition to the anti-aligned class of KBOs that executes resonant-secular dynamics, the Planet Nine hypothesis also entails a population of apsidally aligned objects. The biggest difference between the synthetic Kuiper Belts created in the two numerical experiments shown in Figure 1 lies in the behavior of these objects, whose dynamics are dominated by purely secular effects. Almost no aligned objects survive the full integration time in the narrow Kuiper Belt population; in fact, the few that do survive are found in high order resonances with Neptune. In contrast, the broad Kuiper Belt retains a large population of high-q aligned objects.

Interestingly, the mechanism that depletes the lower-q aligned population in both distributions is a combination of secular dynamics and scattering interactions with Planet Nine and Neptune. To envision this process in physical space, consider the geometry of a typical orbital configuration occupied by a stable particle belonging to the dominantly secular, apsidally aligned population of long-period KBOs (Figure 3). The pink (outer) ellipse represents Planet Nine's orbit, while the blue (inner) ellipse is an example orbit of an apsidally aligned test particle (note that we are operating under the assumption that Planet Nine and the test particle lie in the same plane). Clearly, as the relative longitude of perihelion of the object librates, the inner blue orbit explores some portion of the region enveloped within Planet Nine's pink orbit. Once the orbit of the object crosses the orbit of Planet Nine, however, the object is in danger of experiencing nearly tangential close encounters with Planet Nine.

For a fixed a < a₉ of the test particle, there exists a protected region in the e–Δϖ phase plane, in which the object can experience safe librations within the orbit of Planet Nine. If a > a₉, there exists a similar protected region in which the orbits of the aligned objects engulf the orbit of Planet Nine without crossing it. In configurations other than these apsidally aligned geometries, however, close encounters must occur.

Due to secular interactions with Planet Nine, these aligned objects execute libration cycles in the e–Δϖ phase plane. Depending on the initial q, some of these loops are fully contained within the protected region of phase space, that is, the colinear orbits of these objects remain entirely within or outside of Planet Nine's orbit, and in this case the objects survive the full integration time. An example of such a trajectory is shown in purple in the left panel of Figure 4.

Zoom In Zoom Out Reset image size

For objects initialized on Neptune-hugging orbits, however, secular evolution drives the particle toward an orbit-crossing configuration with Planet Nine. Near this boundary, secular theory formally breaks down, as close encounters ensue. Intriguingly, our simulations show that despite repeated close encounters, nearly tangential objects are not removed from the system by Planet Nine. Instead, these interactions chaotically shift their positions in the phase plane. As a result, the object crosses the tangential collision curve and wanders into the orbit-crossing regime with Planet Nine. Once the object moves away from the tangential collision curve, secular theory regains its accuracy (Gronchi 2002): the object finds itself on a relatively smooth secular trajectory that drives it into the inner solar system via a rapid decrease in q. An example of such dynamics is shown in pink in the left panel of Figure 4. This transition from an initially stable secular cycle to an unstable one is the primary mechanism which depletes the lower-q population of aligned objects, and qualitatively explains the lack of long-term stable aligned objects in the top panel of Figure 1.

The tangential collision curve (blue boundary in the left panel of Figure 4) and the entire e–Δϖ phase plane presented in this plot are mapped for a fixed value of a (in this case, the semimajor axis of the objects shown is a = 345 au). To compute the collision curve for a particular a, we find the value of the eccentricity of the test particle for which its orbit first touches Planet Nine's orbit, given a fixed Δϖ (marked by stars in Figure 3). We then iterate this procedure for varying values of Δϖ. The location of the tangential collision curve is thus determined by the semimajor axis of the object. In fact, the area of the non-orbit-crossing region and the width of the allowed longitude of perihelion librations is dependent on the a/a₉ ratio. As the semimajor axis of the test particles approaches that of Planet Nine, the region in which the orbits of these objects can remain without crossing Planet Nine's orbit becomes more and more limited, until it vanishes in the limit a → a₉.

Indeed, this effect is evident in our simulation results. In the right panel of Figure 4, we show the Δϖ–a phase plane for all objects that survive the full 4 Gyr integrations in the planar simulations. The dark blue objects, which librate about Δϖ = 180°, are the anti-aligned objects. The aligned objects, which librate about Δϖ = 0°, are denoted in pink. The blue curve overlaid on the simulation results represents the tangential collision region with Planet Nine. To compute this boundary, we construct the appropriate tangential collision curve for a fixed a in the e–Δϖ phase plane, as shown in Figure 4, locate the maximum libration width (denoted with a dark blue dot), and map it onto the Δϖ–a phase plane.

The regions outside of this collision curve allow for orbital configurations in which orbit crossing does not occur. Within the region bounded by the curve, however, all orbital configurations are marked by orbit crossing with Planet Nine. As we discussed earlier, it is important to note that close encounters with Planet Nine do not directly destabilize these objects. Instead, these interactions simply act to stochastically modulate the KBOs within the Planet Nine orbit-crossing regime. Once in this region, however, objects are carried into the inner solar system via secular trajectories and are ejected or absorbed by the Sun shortly after. Thus, although the tangential collision curve itself does not signal instability, this boundary represents a useful proxy for dynamical lifetime, as objects that cross this curve are destabilized in tens of millions of years, a timescale much shorter than the lifetime of the solar system.

As a consequence of this evolution, the only objects that find long-term stability within the orbit-crossing region are ones trapped in mean-motion resonances with Planet Nine. Interestingly, all of the aligned objects (pink) whose longitude of perihelion librations protrude into the orbit-crossing region in our simulations are in fact in such resonances. The most visible ones are those objects at a = 700 au (1:1 resonance) and at a = 920 au (2:3 resonance). Such resonances are nevertheless rare, meaning that Planet Nine carves out a high-q aligned population due to primarily secular effects.

We construct similar tangential collision curves for different Planet Nine configurations in Figure 5. Approximately fixing q₉, it is easy to see that lowering the semimajor axis and eccentricity of Planet Nine results in a shrinking orbital crossing region. Thus, closer in and more circular Planet Nine's would carve out the aligned and anti-aligned populations for a smaller range of test particle semimajor axes.

Zoom In Zoom Out Reset image size

While our procedure yields well-defined vertical borders on the a < a₉ and a > a₉ sides of Figure 5, it is important to note that the left and right boundaries of these orbit-crossing regions may not be physically meaningful for the aligned population. In the right panel of Figure 4, we still find aligned objects at semimajor axis values of a < 300 au and a > 1100 au. Outside of this region, however, we are now in principle able to generate a population that circulates in Δϖ that does not require resonance for stability. In the a < 300 au region, such objects have low eccentricities and aphelion distances less than the perihelion distance of Planet Nine. Similarly, in the a > 1100 au region, this nonresonant circulating population is composed of objects with perihelion distances greater than Planet Nine's aphelion. In contrast, circulating objects in the orbit-crossing regime must reside in resonances to attain stability.

2.3. The Inclined Model

The results described thus far were found in our simplified, semi-averaged simulations. In order to verify that the computed dynamics is robust, we considered a more realistic numerical experiment. In this set of simulations, we ran test particles distributed in the same manner as above, but with a uniform inclination distribution i ∈ [0°, 40°], and in the presence of an inclined Planet Nine, with a₉ = 700 au, e₉ = 0.6, i₉ = 20°, and ω₉ = 110°. In addition, we treated the four giant planets as active bodies, foregoing the J₂ approximation made in our planar simulations, and appropriately reducing the time step to 1 year. Both the narrow and broad Kuiper Belt simulations contained about 10000 particles, as before.

Figure 6 shows the inclined simulation results on the Δϖ–a plane and can be readily compared to the planar case shown in the right panel of Figure 4. Clearly, the simulated particles exhibit a greater degree of stochasticity in their evolution; however, the general dynamical behavior is preserved. Interestingly, the libration width of the anti-aligned objects (blue) in Figure 6 appears to be smaller than that in the corresponding right panel of Figure 4. A possible explanation for this is that due to modulation from Planet Nine's inclination, the edges of the anti-aligned orbital cluster breakup into a rapidly chaotic zone, reducing the width of the stable region (Batygin & Morbidelli 2017). A more detailed analysis of this phenomenon would present an interesting avenue for future work.

Zoom In Zoom Out Reset image size

Note that the dashed blue curve is the same planar collision curve as that shown in Figure 4. The meaning of this boundary is now more complicated. Because Planet Nine and the test particles could have substantial mutual inclinations in these more realistic simulations, there are orbital configurations within the region bounded by this curve in which the orbits of the object and Planet Nine do not cross. Even so, due to the oscillations in inclination and longitude of ascending node, objects within the dashed curve predominantly reside in a domain in which close encounters are more likely. As a result, despite possessing prolonged dynamical lifetimes as compared to their counterparts in the strictly planar simulations, nonresonant objects found in this region will ultimately experience instability as well, and will be removed from the system.

In this work, we considered two synthetic Kuiper Belts with initially different perihelion distance distributions: a "narrow" (q₀ ∈ [30, 36] au) and a "broad" (q₀ ∈ [30, 300] au) population. Our numerical simulations show that after a 4 Gyr integration, the final distribution of perihelion distance of these two populations remains distinct. Although both Kuiper Belt incarnations retain a low-q anti-aligned population, only the broad distribution develops a high-q aligned population.

Qualitatively speaking, this discrepancy can be explained by Planet Nine-induced secular dynamics. The very high-q anti-aligned objects do not survive the full integration time because they develop large librations and their trajectories are carried into the inner solar system. Analogously, the low-q aligned objects are ejected once driven to q < 30 au. Thus, we identify two permanently stable populations in the outer Kuiper Belt: the low q anti-aligned objects and the high-q aligned objects.

While our initial simulations contained only a low inclination Kuiper Belt and a noninclined Planet Nine, we ran a more realistic set of simulations, in which we considered an inclined Planet Nine and the effect of the giant planets as active bodies. Due to the good agreement between the results in the planar and inclined models, we are able to conclude that the planar approximation provides a reasonable representation of the dynamical behavior induced by Planet Nine. We can thus consider the implications of the results of our planar simulations as being representative.

As different parameters of Planet Nine "carve out" differing stable aligned populations, as shown in Figure 5, studying the distribution of Δϖ as a function of semimajor axis for distant KBOs as more are found presents an interesting constraint on the semimajor axis of Planet Nine. In the census of available observations, however, the sample of librating aligned objects is contaminated by objects that have a Δϖ ∼ 0, but are in fact circulating. The population of observable circulating objects dominates the observable aligned population, as the aligned objects are only stable at very high (and thus hard to observe) q values. In light of this, a closer examination of the distribution of circulating objects may allow for stricter conclusions on the semimajor axis of Planet Nine and presents an important avenue for future work.

More importantly, we find that the two incarnations of the distant Kuiper Belt are indeed distinct. As all low-q objects are depleted due to secular interactions with Planet Nine, the main difference between the narrow and broad Kuiper Belt distributions is in the existence of high-q aligned objects. Thus, the presence of high-q aligned objects in our solar system would not only lend support for the existence of Planet Nine but would indicate that our current Kuiper Belt stems from an initially wide q-distribution. This, in turn, would signal that Planet Nine's formation likely involved stellar encounters, as this is the most likely mechanism that would create such a spread out Kuiper Belt.

We note that the above analysis entails a caveat. In our model, we only consider the behavior of the outer Kuiper Belt; once we initialize the ∼10000 particles, those are the only objects whose evolution we consider. That is, we do not take into account objects that may have originated in the inner Kuiper Belt or Scattered Disk, and later were scattered out to sufficiently high a, thus falling into our population of interest. These objects could replenish certain regions of parameter space, such as the low-q aligned or the high-q anti-aligned regions, adding an influx of objects to the existing populations.

The objects that may be scattered into the outer solar system, however, will not become members of the two permanently stable populations identified in this work (low-q anti-aligned and high-q aligned). Instead, these objects will replenish the orbital configurations that are eventually driven to centaur-like orbits and are then rapidly ejected from the solar system. As a result, these objects are only metastable, as their lifetimes are generally shorter than one libration in Δϖ, a process which lasts <1 Gyr. Even so, carefully modeling the flux of these scattered objects in the outer solar system under the influence of Planet Nine would be a possible direction of future research.

Under the assumption that our current and future observational census is dominated by the permanently stable objects, we hope that upon discovery of additional distant KBOs, results identified in this study can be used to place meaningful constraints upon the initial distribution of the outer Kuiper Belt. In turn, the structure of this primordial population will allow us to draw conclusions regarding the formation pathway of Planet Nine, and by extension, the Sun's birth environment.

We thank Elizabeth Bailey, Christopher Spalding, and Juliette Becker for useful conversations, and Fred Adams for a careful review of the draft. We thank the anonymous referee for valuable comments that led to the improvement of the manuscript. T.K. is grateful to the Johnson and Johnson WAVE Fellowship for funding as well as to the Caltech Student-Faculty Program for their support of this work. K.B. acknowledges the generous support of the David and Lucile Packard Foundation.