Long-term Dynamical Stability in the Outer Solar System I: The Regular and Chaotic Evolution of the 34 Largest Trans-Neptunian Objects

We begin our analysis by characterizing the overall orbital stability of the outer solar system. To this aim, we performed a frequency map analysis (FMA, Laskar, 1990; Laskar et al., 1992) of wide regions of the $a$-$e$ phase space plane (semimajor axis vs eccentricity). A similar analysis has been performed earlier by Robutel & Laskar (2001). Here we reproduce their analysis but we increase the resolution in the sampling of the phase-space values corresponding to the outer solar system, focusing in the regions occupied by the majority of the DPs used in this work.

To produce the diffusion maps of Figs. 1, 2, 3, and 4, we performed short-term numerical integrations (for the duration of $\sim 1.31\,$Myr, or approximately 8000 orbital periods of Neptune) of thousands of test particles distributed in different homogeneous grids that cover the heliocentric $a$-$e$ phase-space plane.

The test particles are subject to the gravitational potential generated by the Sun and the four giant planets, where the mass of all the interior planets (including the Moon and Ceres) are added to that of the Sun. The initial conditions for all the massive bodies, on the Julian day 2458176.5 corresponding to February 27, 2018, were retrieved from the JPL’s Horizons system²²2https://ssd.jpl.nasa.gov/horizons.cgi.

We create a diffusion map of the outer solar system by performing a frequency analysis of the evolution of each of our test particles; we focus on the quantity:

where $a$ and $\lambda$ are the semimajor axis and mean longitude of the particle, respectively, while the complex value of $z(t)$ represents a rough approximation of the position of the particle (the details regarding the relevance of $z(t)$ can be found elsewhere, e.g. Robutel & Laskar, 2001; Muñoz-Gutiérrez & Giuliatti Winter, 2017).

The algorithm used for the FMA was that of Šidlichovský & Nesvorný (1996), which provides an approximate decomposition of $z(t)$ as follows:

For particles that lie in keplerian orbits $\lambda(t)=nt$ (where $n$ is the mean motion of the particle); also $N=0$, $a=|\alpha_{0}|$, and $n=\nu_{0}$ will remain constant.

In general $N>0$, while $|\alpha_{0}|$ and $\nu_{0}$ need not remain constant; but, for moderately stable orbits, $|\alpha_{0}|$ and $\nu_{0}$ will remain nearly constant, and $a$ and $n$ will remain close to them (i.e. $\alpha_{0}\gg\alpha_{k}$); on the other hand, for unstable orbits these values will evolve quickly over time.

In particular $\nu_{0}$ (the mean frequency of the particle) will change for unstable orbits; a measure of the stability of the orbits can then be defined as the diffusion parameter:

where $\nu_{01}$ and $\nu_{02}$ are the main frequencies of a particle in each adjacent time interval of length $T$. Small values of $D$ indicate a stable trajectory, while larger values are indicative of unstable orbital evolution, resulting from strong variations of the main frequencies, which are characteristic of the orbital erratic nature. For a more detailed description of the FMA the interested reader is referred to the works by Laskar (1990); Laskar et al. (1992); Laskar (1993).

In order to compute the change of the main frequencies, we perform the frequency analysis for each particle in the intervals $0<T<0.655$ Myr and $0.655<T<1.31$ Myr.

In the maps of Figs. 1, 2, and 4 we color each ($a$,$e$) pair according to the logarithm of the diffusion parameter of its test particle such that redder colors indicate more unstable orbits, while bluer colors stand for the more stable orbits. Particles that are lost from the simulation before it finishes, mainly due to ejections or collisions with a planet, are colored white.

These diffusion maps let us identify regions where particles could survive in a long-term basis, while permitting us to easily identify the location and width of mean motion resonances with Neptune (MMRs). In the maps, the solid green lines indicate orbits with a perihelion $q=a_{\it Neptune}$ (i.e. it represents the crossing line of Neptune, with the implied risk of collision). The green dashed-lines indicate the distance from which the influence of Neptune starts to be dominant, as has been used in previous works (Gladman & Duncan, 1990; Muñoz-Gutiérrez et al., 2019); numerically the green dashed line can be expressed as: $q=a_{N}+2\sqrt{3}R_{\it HN}$, where $a_{N}$ and $R_{\it HN}$ are the semimajor axis and Hill radius of Neptune, respectively.

Figs. 1 and 2 show the diffusion maps produced from the integration of $34\,162$ particles distributed in two homogeneous grids with the following conditions: in Fig. 1 the heliocentric $a$ is sampled from 38 to 58 au, with a step-size $\Delta a=0.04\,$au, while $e$ is sampled from 0 to 0.3 with a step-size $\Delta e=0.01$. For Fig. 2 $a$ is sampled from 58 to 82 au while $e$ covers from 0.3 to 0.6 in steps of the same size as in Fig. 1.

In order to focus on the stability and shape of MMRs with Neptune, all of the remaining orbital elements: inclinations (from the plane of the ecliptic), $i$, arguments of pericenter, $\omega$, and mean anomalies, $M$, are set to zero. The longitudes of the ascending nodes, $\Omega$, are also nominally set to zero; while meaningless when $i=0^{\circ}$, they are a necessary input for the Mercury integrator.

The initial location of 30 out of the 34 DPs in our sample are indicated by black dots in both maps; the other 4 DPs all have $a>100$ au, and thus were left off of these figures.

We immediately notice that, for the majority of objects in our sample, a long-standing stability is expected based solely on their location on the phase-space. However these maps were made for an $i=0^{\circ}$; when studying orbits with $i\gtrsim 15^{\circ}$ the resonances move, slightly, outwards. We can see for example how two of the Plutinos (Pluto and Ixion) appear to be outside of the 3:2 MMR in Fig. 1, even more, they appear to be in a highly unstable region. However, using a map with a more appropriate inclination it is easy to see that this is not the case; in Fig. 3 we present a small diffusion map using $i=17.12^{\circ}$, the heliocentric inclination of Pluto; this map clearly shows that Pluto is firmly inside the 3:2 resonance. The shifting effect of the location of MMRs, due to the change in $i$, will occur for all resonances; for example, we can see the opposite case for $2015\,{\rm KH}_{162}$: in Fig. 2 it appears to be clearly located inside the 3:1 MMR, however, we were unable to identify a stable libration of any of the possible second-order resonant arguments for it; therefore, despite its suggestive location on the map, due to its high inclination, $2015\,{\rm KH}_{162}$ is found to be clearly non-resonant.

Based on the diffusion maps, it is possible to understand the low stability of some of the DPs in our simulations, for example, $2007\,{\rm OR}_{10}$ and $2015\,{\rm RR}_{245}$, are very clearly located in the broad unstable region (close to the dashed line: $a_{N}+2\sqrt{3}R_{\it HN}$), while $2002\,{\rm MS}_{4}$ and $2002\,{\rm UX}_{25}$ are located in the unstable regions near both sides of the 5:3 MMR. On the other hand $2005\,{\rm RN}_{43}$, Quaoar, Varuna, and $2010\,{\rm KZ}_{39}$ represent examples of objects with low $e$ located in very stable regions.

In order to better understand the effect that a non-zero inclination would have on the general aspect and on the stability picture drawn from the diffusion maps, we ran a set of simulations including $48\,235$ total test particles covering different patches of the $a$-$e$ phase space. In this case, the inclination of each particle was assigned at random, with values between 0^∘ and 50^∘, in order to represent the broad distribution of inclinations in our DP sample (where the minimum and maximum $i$ values are 7.552^∘ and 43.993^∘, respectively). Besides, the choosing of such values allow us to better identify the resonances for objects that orbit far away from the plane.

Since each of the 34 DPs also has its own value for the other angular parameters, $\omega$, $\Omega$, and $M$, we decided to assign random values between 0^∘ and 360^∘ for the corresponding initial conditions of each test particle in the map. The exact values used for these angles are probably not very important; but, where there is any difference, this approach has the advantage of showing us, locally, the range of stability available to objects with such $a$-$e$ values. We will study a more representative determination for each of the objects of our sample in section 3.4.

In Fig. 4 we show the resulting diffusion map for random inclinations and angular elements, as described above, for a larger region of the phase-space than those of Figs. 1 and 2 (though both these regions are covered in the simulations of this section). As in the previous maps, the solid green line delimits the collision region with Neptune, while the dashed green line delimits the region where Neptune’s gravity is dominant. The location of 30 out of 34 DPs is again marked by black dots.

It is interesting to note how the general form of the MMRs changes with respect to the zero inclination case. With random angles and inclinations, the definition of the resonances gets blurred; however, the stronger MMRs, whose ratios are labeled at the top of the figure, are still easily identifiable. It is evident that the locations of the resonances are perturbed by the value of the inclination, thus we can see how both Pluto and Ixion are indeed located within the limits of the 3:2 MMR.

This diffusion map constitutes a better representation of the stability of the Kuiper belt as a whole, not only for the DPs but for objects of all sizes in the various dynamical families present in the region covered by the map (namely the Classical belt, the Scattered disk and some Resonant populations). From the map, it is possible to estimate, statistically, the diffusion time of Kuiper belt objects as a function of their perihelion.

In Figure 5 we present the fraction of stable orbits as a function of their perihelion, for diffusion times varying between 0.1 to 5 Gyr. The diffusion time, $T_{\mathrm{Diff}}$, is obtained from the diffusion parameter as $T_{\mathrm{Diff}}=1/(DP)$, where P is the period of the orbit (see for example Gaslac Gallardo et al., 2020; Roberts & Muñoz-Gutiérrez, 2021, for other recent applications of the diffusion time on different systems). To determine the stable fraction, we counted the number of particles, within our sample, for each perihelion band between 28 and 50 au, whose orbits have diffusion times above 6 different time limits: 0.1, 0.2, 0.5, 1, 2, and 5 Gyr. These limits are significant in the sense that any object with a diffusion time below 0.1 Gyr is expected to be lost very quickly, and thus those orbits are not expected to have a meaningful contribution to the present day population; on the other hand a diffusion time above 5 Gyr represent longer than the lifetime of the solar system.

We note that our sample is incomplete; besides having a hard limit at a semimajor axis of 82 au, there are additionally a few gaps below 30 au and above 40.6 au. Nonetheless, the sample is numerous enough as to provide a good representation of the expected stable fraction of Kuiper belt objects as a function of perihelia. For example, we expect that, from an initial population with perihelia of 38 au: almost all objects would be stable for more than 100 Myr, 10% of this initial population would be unstable in time scales of 1 Gyr, and nearly 30% of the population would be lost over the age of the solar system.

In order to obtain a better estimate, regarding the specific stability of each object in our sample of large TNOs (i.e. beyond the global stability suggested by the general configuration of the diffusion map of Fig. 4), we performed additional stability simulations of 100 mass-less clones of each DP. The clones were generated by randomly assigning orbital elements within the intervals $a_{\it DP}\pm 0.04$ au and $e_{\it DP}\pm 0.01$, while the angular elements, $i_{\it DP}$, $\omega_{\it DP}$, $\Omega_{\it DP}$, and $M_{\it DP}$, were randomly generated within intervals of $\pm 1^{\circ}$ around the initial values of each DP’s orbit.

We calculated the value of the diffusion parameter for each DP, as the average of the parameters of its 100 clones. The results are shown in Table 2. In the last column of Table 2, we give a prediction on the stability of each object based on the value of its diffusion time.

Our definition of $T_{\mathrm{Diff}}$ represents an estimate of the time scale for a change of unity in the period of an orbit; from Kepler’s Third Law it is clear that changes in period (or mean motion) originate from changes in the semimajor axis of the orbit. In this work we assume that a drastic change in the orbital motion would be a quarter of this, and we will set our criteria by comparing our stability time scales with $T_{\mathrm{Diff}}/4$.

Based on the previous argument, we call “Very Stable” those orbits with $T_{\mathrm{Diff}}>20$ Gyr, this is, those for which a significant change in their mean motion would not occur over the remaining life of the solar system. We call “Stable” those orbits for which $4<T_{\mathrm{Diff}}<20$ Gyr, this is, an appreciable change in their orbits can occur over the age of the solar system, but requiring more than 1 Gyr. Finally, “Unstable” are those orbits with $T_{\mathrm{Diff}}<4$ Gyr, i.e. orbits which will experience drastic changes in less than 1 Gyr of evolution.

The 1 Gyr limit may seem arbitrary, however, many dynamical works which consider a time scale for the study of continuous and stationary phenomena in the solar system, a so called steady state, use this same scale of 1 Gyr (e.g. Levison & Duncan, 1997; Tiscareno & Malhotra, 2009; Nesvorný et al., 2017; Yu et al., 2018; Nesvorný et al., 2019). Indeed, the present study arises from 1 Gyr integrations, which were performed for the characterization of the contribution of DPs to the delivery of short-period comets, from the Kuiper belt to the inner solar system (Muñoz-Gutiérrez et al., 2019). We will present the results for the DPs of those simulations, as well as several complementary 1 Gyr integrations, in the next sections.

From Table 2 we can see that $2007\,{\rm OR}_{10}$ is the DP with fastest diffusion rate, with important changes in a time scale of 150 Myrs. Also, Haumea, $2002\,{\rm MS}_{4}$, $2007\,{\rm UK}_{126}$, and $2010\,{\rm EK}_{139}$ have significant changes in less than 1 Gyr.

On the other extreme we find that the diffusion time scale for $2012\,{\rm VP}_{113}$ is greater than 5000 Gyr (5 Tyr!). In fact 21 DPs (out of 34) have stabilities beyond 5Gyr (beyond the remaining lifetime of the solar system).

Our dynamical model is composed, in all cases, by a central star with a mass equal to $(1+\epsilon)$M_⊙, where $\epsilon$ is a small quantity that represents the mass of all the terrestrial planets plus the Moon and Ceres; also, we consider the gravitational influence of the 4 giant planets and the 34 massive DPs.

We ran a total of 200, 1 Gyr long simulations, using the hybrid symplectic integrator from the Mercury package (Chambers, 1999). The initial conditions were the same as in section 3.1.

Most of the simulations we are using (177 out of 200) were designed for the study of the evolution of cometary nuclei and each one included a different set of some hundreds of test particles; the results for the evolution of test particles on many of these simulations were reported in Muñoz-Gutiérrez et al. (2019).

The initial conditions of the massive objects in all 177 simulations were identical, but each simulation differs from the others by the specific set of test particles which were included; while the nominal time step for all integrations was 400 days, the presence of close interactions forced the integrator to modify the time steps when required, this implies that, each specific set of test particles resulted in a different specific set of time steps for each simulation, which in turn resulted in different evolutions for the massive objects in each case. The differences in evolution are barely noticeable for the giant planets, but sometimes resulted in wildly different positions for the DPs; since the error tolerances for the integrations were stringent (with a tolerance accuracy parameter of $10^{-10}$), these differences represent the evolution of very close orbits in a chaotic region of the phase-space, and indeed all these results are consistent with the very same initial conditions within very small error bars.

To the previous sample we added 23 extra simulations, to end up with an even 200. Each of the new simulations had the same initial conditions for the massive objects, and included a token test particle at 100 au; instead of depending on interactions on test particles, which slow down the integrations, each new simulation has a slightly different initial time-step (between 389-412 days), which again resulted in different evolutions for each of our DPs.

The statistical studies we present here are possible due to the different results obtained in different simulations, despite the identical initial conditions used for the massive objects. Those differences could naïvely put into question the validity of our results, or be attributed simply to the finite precision of the integrator; however, our usage of a highly reliable and widely accepted integrator (the Mercury package of Chambers, 1999, with over 1000 citations), could easily quell at least the later argument.

On the other hand, it is understood that numerical errors are unavoidable in integrations and even more, while studying a very unstable system, any small imprecision may lead to important divergences over time, i.e. dynamical chaos (e.g. Murray & Dermott, 1999). This doesn’t invalidate the meaning of such integrations; in fact, the numerical errors in the short-term are smaller that the uncertainties in the observational data for the initial conditions of the objects. On the longer term, however, one could argue, that only the very stable orbits will exist in the real system, since any unstable orbit would be expected to have been expelled from the system long ago; leading to the false conclusion that the divergences found in our numerical experiments cannot have an observational equivalent nowadays. On the face of this argument, we must remember that the solar system has thousands of objects that we are not including in our simulations, which can chaotically interact at different levels, and a small interaction with another close or even not so close component may lead to a divergence more significant that those produced by numerical errors.

To these we must add the presence of extra solar perturbations which are not easily predictable neither accounted for, such as: moderately close stellar flybys or interstellar visitors (’Oumuamua, 2I/Borisov); any of which also can destroy the stable orbits expected in a perfectly steady-state system, or the ones produced by a perfect hypothetical integrator.

In a sense, the instability of orbits in the trans-Neptunian region is hardly unexpected, since it has been long known that some apparently stable orbits, over time spans of Gyr, can suddenly and quickly become chaotic on time spans of Myr (Torbett, 1989; Holman & Wisdom, 1993; Duncan et al., 1995); this phenomenon supplies comets to the inner solar system. At first glance the instability of larger objects, such as a DP, may seem surprising, however, there is nothing that prevents, in principle, such massive objects from suffering the same fate as cometary nuclei, at least in a statistical sense. To put it simply, the constant presence of new comets shows that the solar system, and in particular the trans-Neptunian region, is not in a perfectly smooth and stable steady-state.

Of course, we do not expect all DPs to be unstable (see Table 2, as well as Section 5); some of the most stable orbits (such as Pluto’s) will not be destabilized by inaccuracies resulting form numerical errors. Still, for objects surrounded by a not-so stable phase space, the presence of slight (but constant) perturbations can disrupt their stability.

We now present the statistical results from our long-term simulations; this is, the results from 200 different realizations for each of the 34 DPs (for clarity, in the rest of this work, we will call ”orbit” each of the 200 realizations for the evolution of each DP; i.e. we will be studying 6800 different orbits). As we did in section 3.4, we are interested in providing a general classification for each object (stable or unstable), as well as their possible evolutions, since, as we will show, some of the DPs are likely to leave the trans-Neptunian region during the remaining lifetime of the Sun.

An ideally stable object would present 200 regular orbits. This is, each orbit would have only small variations on its main orbital parameters during our entire simulation; this also implies that the orbits on all of the 200 simulations would be equivalent.

In order to define what a regular (or irregular) orbit is, we need to determine the limits in the deviations of the orbital parameters that we will tolerate and still call an orbit “regular”. We first calculate the standard deviations of the orbital parameters $a$, $e$, and $i$ in the first 25 Myr of each integration; this is, for each object we calculate the $\sigma_{a}$, $\sigma_{e}$, and $\sigma_{i}$ using the data of the first 25 Myr for all 200 of its orbits. In Fig. 6 we plot the values of $\sigma_{a}/a_{0}$, $\sigma_{e}$, and $\sigma_{i}$ versus the initial orbital elements of each DP, this is, against $a_{0}$, $e_{0}$, and $i_{0}$. From Fig. 6 we can see that most of the DPs have similar values of their corresponding dispersion. Using the values of Fig. 6 we now define a characteristic deviation for each orbital element simply as the average of the above standard deviations, thus we have $\left<\sigma_{a}/a_{0}\right>=0.00372$, $\left<\sigma_{e}\right>=0.00836$, and $\left<\sigma_{i}\right>=1.148^{\circ}$. With these characteristic values, we are now able to determine, for each DP, how many of the orbits are regular in the corresponding parameter; we do this by defining an individual orbit as regular if its orbital parameters ($a$, $e$, and $i$) remain within 5 characteristic deviations from their original value for the entire simulation, i.e. they remain within the intervals $a_{0}\left(1\pm 5\left<\sigma_{a}/a_{0}\right>\right)$, $e_{0}\pm 5\left<\sigma_{e}\right>$, and $i_{0}\pm 5\left<\sigma_{i}\right>$ for 1 Gyr. Note that we do not use any restrictions on the evolution of $\omega$, $\Omega$, or $M$.

In Table 3 we show the number of regular orbits in each orbital parameter (out of 200), according to the previous procedure. We also include the number of orbits which were regular in all 3 parameters simultaneously. We will consider a DP to be stable if it fulfills the following two conditions: if at least 80% of its orbits are stable, and if none of its orbits ended up being lost in the simulation. We consider 3 kinds of losses: objects whose $a$ increases beyond $10\,000$ au (escape from the solar system), objects whose distances to the sun decrease under 1 au (big comet), or collisions with a planet.

Of the 15 objects classified as unstable in Table 3, 13 were classified as such for having between 0 and 95 stable orbits (clearly less than 80%); the other two unstable objects, $2010\,{\rm RF}_{43}$ and $2002\,{\rm UX}_{25}$, were classified as such since, despite having 170 and 173 stable orbits, some of their orbits were lost (5 and 1, respectively).

An additional result of the analysis of unstable and semi-stable orbits is the presence of resonant orbits; we have added an additional category to allow for the specifics of the behavior of such orbits. There are 6 objects in Table 3 that have been classified as resonant; 3 of those (Pluto, $2010\,{\rm JO}_{179}$, and $2003\,{\rm AZ}_{84}$) could easily be considered as stable, while the other 3 (Ixion, $2002\,{\rm TC}_{302}$, and $2010\,{\rm EK}_{139}$) would be considered unstable according to our previous selection criteria. It is worth to note that, although the last 3 objects could be classified as unstable, they are not really fully unstable, as they remain locked in their resonance 94% - 100% of the time (as shown by the number of orbits regular in $a$). The reduced number of regular orbits for such objects (between 0 and 2) is actually due to their resonant nature, which leads to large $e$ variations (and sometimes of the inclination) above the 5-$\sigma$ level defined before. This new category highlights their resonant behaviour which dominates for most of their orbits and is the most defining characteristic of their evolution; it is for this reason that we have also included the three more-stable resonant objects in this category.

In figure 7 we present the initial and final perihelion vs eccentricity for all 200 orbits for 32 DPs (Sedna and $2012\,{\rm VP}_{113}$, have perihelion of 76.2 and 80.5 au, outside the range of this figure). We can see that for some DPs (the stable ones; blue) the final values for each orbit (dots) lie all very close to one another, and also lie very close to the initial condition (circle); for the resonant DPs (green) we see that the initial and final values lie in long lines corresponding to constant $a$, showing the evolution in $e$ common in resonant orbits; for both stable and resonant orbits there are frequently a few stragglers far from the initial values. On the other hand there are DPs (unstable DPs; red) where the final points change both in $e$ and $a$, sometimes spreading over large areas of this plane. This highlights the differences of the 3 classifications defined for Table 3.

The last column in Table 3 gives us an idea of the quality of the determinations presented in Table 2. Of the 28 non-resonant objects, 21 remain in the same category (considering “Very Stable” objects to be in the “Stable” category). Of the other 7 objects: $2010\,{\rm RF}_{43}$ changes, but it is near the limit of both classifications, while the other 6 change category completely (Orcus, $2013\,{\rm FY}_{27}$, $2015\,{\rm RR}_{245}$, $2014\,{\rm UZ}_{224}$, $2002\,{\rm UX}_{25}$, $2005\,{\rm QU}_{182}$); in all 7 cases the new determination is unstable while the old one was either “Stable” (2 cases) or “Very Stable” (5 cases). This difference occurs mostly because of the limitations of the FMA, which extrapolates the amount that the orbit will change based on a short-term integration in which the orbit does not deviate substantially from its initial conditions, i.e. while the perturbations are small their behaviour can be modeled by the FMA, but when they grow large enough to reach the direct influence of Neptune (within 3 hill radii of Neptune) the stability decreases dramatically and the short term estimates are no longer predictive.

Regarding the resonant DPs: we find that Pluto, $2010\,{\rm JO}_{179}$, and $2003\,{\rm AZ}_{84}$ were stable in the FMA and remain stable in the long-term integration. On the other hand Ixion and $2002\,{\rm TC}_{302}$ have regular evolution for $a$ and $i$ but their $e$ changes drastically (as is expected for many resonant objects). Finally $2010\,{\rm EK}_{139}$ was unstable in the FMA and shows erratic behavior, in both $e$ and $i$, in the long-term integrations while remaining in the resonance (with constant $a$).

The 4 Plutinos represent an interesting subset of DPs, showing a wide variety of behaviors. While Pluto is completely (100%) stable, $2003\,{\rm AZ}_{84}$ would be considered stable (94%), and Ixion would be considered unstable (0% stable); yet, out of the 600 orbits corresponding to these 3 DPs only 3 orbits leave the resonance. Orcus, on the other hand, is truly unstable; not only 39 orbits leave the resonance (which, in itself negates the argument that it is locked in the 3:2 resonance); but, out of those 39 orbits, 26 leave the simulation entirely (18 leaving the solar system and 8 falling to the inner solar system).

A final note should be made about the definition of stability: one could measure stability by the consistency in $e$ and $i$, or by the consistency of the orbital period (or $a$), or by the ability to remain in the Kuiper belt; the least stable object would be different in each case. Ixion is the DP with the smallest amount of regular orbits (0%) while remaining in resonance, $2005\,{\rm QU}_{182}$ is the DP that is less likely to remain close to its original $a$ (0%) while Haumea is the most unstable DP considering the number of orbits that remain in the solar system (37.5%).

Table 4 shows the list of objects ejected and the end state of their ejections. Overall 378 (out of 6800; 5.56%) of the orbits were ejected from the solar system; this implies that the solar system is still evolving and 1 Gyr into the future will have lost $\sim$2 of these objects. Most of the ejections (305 out of 378; 80.7%) were into outer space, with an important fraction (69 out of 378; 18.2%) falling into the inner solar system. Finally a few of the orbits (4 out of 378; 1.1%) end up colliding with a giant planet (2 with Jupiter and 2 with Neptune).

Encounters between DPs and giant planets are to be expected over the lifetime of the solar system. In fact, collisions of even larger objects are frequently used to explain some of the striking characteristics of the giant planets such as: the capture of Triton by Neptune (e.g. Agnor & Hamilton, 2006), the tipping of Uranus’s rotational axis(e.g. Slattery et al., 1992), Jupiter’s core composition (e.g. Liu et al., 2019) and even a possible origin for Saturn’s rings(e.g. Dubinski, 2019).

While the 4 collisions presented in table 4 lead to interesting possibilities, it is potentially much more interesting (and important) the presence of 69 orbits that fall into the inner solar system. In the case of such an event, the infalling DP would become the equivalent of a huge short period comet, with a longer physical life span that would permit tens of thousands of orbits and could potentially endanger Earth as we know it. Unfortunately, we cannot follow the evolution of these 69 orbits, since the setup of our simulations does not allow us to follow the evolution of objects with orbits inside 5 au (inside Jupiter’s orbit).

These 69 orbits represent an average infalling rate of 1.01% per DP per Gyr, and an overall 34.5% probability when considering all 34 DPs we are tracking. This probability is slightly smaller than the rate of cometary infall from the Kuiper belt found by Muñoz-Gutiérrez et al. (2019), that amounts to 4.6% per cometary nuclei per Gyr; this difference is not entirely unexpected since lighter objects are easier to move around. However, the presence of a DP in the inner solar system, would certainly be much more spectacular.

As far as life on Earth is concerned, any object larger than about 20 km could, potentially, become a game changer. There are expected to be approximately $10\,000$ times more objects between 20 and 200 km than those on our sample (Fraser & Kavelaars, 2009; Dones et al., 2015), this implies an infall rate of a few thousand objects larger than 20 km each Gyr (or equivalently, a few every Myr).

Specifically, the most likely object to reach the inner solar system within the next Gyr is Haumea, while $2002\,{\rm MS}_{4}$, $2007\,{\rm OR}_{10}$, and Orcus, have a modest probability and $2015\,{\rm RR}_{245}$, $2013\,{\rm FY}_{27}$, $2007\,{\rm UK}_{126}$, $2003\,{\rm AZ}_{84}$, and $2010\,{\rm RF}_{43}$ have a minor probability of reaching the inner solar system (see table 4). It must be also noted that if Haumea (and to a slightly lesser extent $2007\,{\rm OR}_{10}$) don’t reach the inner solar system within the next Gyr they become less likely to do so, as they are also very likely to be ejected from the solar system. In general a critical step for infalling orbits includes a close encounter with Neptune, but such an encounter is approximately four times more likely to result in an expulsion of the solar system than in an infalling trajectory.

The most likely fate for objects ejected from the Kuiper belt is the expulsion of the solar system. Haumea and $2007\,{\rm OR}_{10}$ are the objects more likely to do so, with a half life of about 1 Gyr; $2002\,{\rm MS}_{4}$, $2015\,{\rm RR}_{245}$, $2013\,{\rm FY}_{27}$, and Orcus are also quite likely to be ejected with a half life of 2 - 4 Gyr, and $2007\,{\rm UK}_{126}$, $2010\,{\rm RF}_{43}$, $2002\,{\rm TC}_{302}$, $2003\,{\rm AZ}_{84}$, $2002\,{\rm UX}_{25}$, and $2010\,{\rm EK}_{139}$, having a non negligible probability of ejection, but with a half life longer than that of the solar system. In summary we expect about 1.5 DPs to be expelled from the solar system in the next Gyr.

When trying to paint a picture of the Kuiper belt 1 Gyr into the future, the most striking feature is the overall loss of approximately 2 objects (1.89 objects lost).

Another characteristic is that there is no object for which its 200 orbits converge in the same point after 1 Gyr (even for Sedna, the slowest and most stable of the DPs, the behavior of the different runs shows we can not keep track of the mean anomaly beyond 10 Myrs). If we ignore the mean anomaly, we can only have some confidence on the argument of pericenter and longitude of the ascending node for Sedna and $2012\,{\rm VP}_{113}$, for all other objects we lose coherence for both quantities before the 50 Myr mark; on the other hand, our simulations do not consider the interactions with external potentials (such as stellar flybys, or the effect of Planet 9); any such encounter would affect more severely precisely those objects that are particularly well behaved in our simulations.

This does not mean that we will be completely unable to recognize some of the orbital characteristics of the trans-Neptunian DPs in the next Gyr:

Focusing only on the 6 DPs we have classified as resonant, we see that they are extremely likely to remain in the Solar System (5.96 out of 6, 99.3%); in fact, they are likely to remain inside their current MMRs (5.915 out of 6, 98.6%), but only slightly less than half will remain with similar values of their main orbital parameters (2.89 out of 6, 48.2%).

Although the 28 non-resonant objects are less likely to remain in the solar system (26.15 out of 28, 93.4%), they are more likely to maintain similar values of their main orbital parameters (20.135 out of 28, 71.9%).

Overall, by focusing on the main orbital parameters, we would be able to recognize about two thirds of the DPs (23.025 out of 34, 67.7%), and about three fourths (26.05 out of 34, 76.6%) if we allow for the resonant objects that remain inside their MMR.

When considering together this work (which deals with DPs and TNOs greater than about 400km) and our previous work (Muñoz-Gutiérrez et al., 2019, which deals with small TNOs in the 2-10 km range), there is a gap for a population of intermediate sized objects (20-400 km). These objects are able to also have a gravitational effect on other objects, but overall we expect their disruptive effect to be much smaller than the one produced by our sample; we would expect that including 100 additional objects (the most massive within this population) would increase the disruption by less than 10% but with an increase of a factor of $\sim$ 5 on the computational effort, while a massive simulation that includes the 1000 most massive objects of this population, would include nigh all of the disruptive effect due to the rapid decay in the mass distribution, and increase the disruption by less than 20% but with an egregious increase on the computer expense.

A different possibility would be the presence of additional large TNOs, not included in our sample (larger than 400 km, or even 1000 km). For example Schwamb et al. (2009) and Shankman et al. (2017) estimate between 40 and 80 Sedna-like objects in the outer solar system; these objects are expected to be similar in size, eccentricity, and perihelia to Sedna. Such population would have very little effect on our 34 object sample; we can see that Sedna is very stable, this implies that no other body from our sample interacts with it, equivalently Sedna does not interact with any other object; also with a perihelion of 76 au it never gets close enough to perturb objects that could potentially get close to Neptune to be strongly disrupted. Also, most of the effect between DPs is secular in nature, and any object with such orbital characteristics would spend too little time close to other objects for the secular interactions to be significant. Interactions between pairs of Sedna-like objects would also be unimportant due to the large volume at such large distances from the Sun.

A potentially more disruptive population would be objects with semimajor axis similar to Sedna, but with perihelia able to reach the Classical Kuiper Belt region; any such object would be much less efficient than the DPs from our sample since they spend only a very small fraction of their time within the Kuiper belt range. Since only about 1% of their time would be spent in the Classical Kuiper Belt region, approximately 100 such objects would be needed to contribute as much as any single object from our sample, and a few hundred would be needed to change our overall conclusions. While objects with this characteristics can not be ruled out, such a large population is not expected to be hidden from us.

While the existence of additional large and distant DPs would decrease the stability of our sample, and our results thus represent a lower limit to the rate at which DPs will be lost, we do not expect there are enough additional interactions to change our conclusions significantly.

This study explores the future evolution of the trans-Neptunian region. It would be interesting to know how this region has evolved to reach its present configuration. It is obvious the number of DPs is falling with time, this also means that the evolution speed is slowing down with time.

According to Muñoz-Gutiérrez et al. (2017) the evaporation rate is proportional to the square of the mass of the disk, in systems with an interior giant planet; however such models ignore the possibility of additional perturbations, such as flybys or giant planet migration. However, while it is possible to do a statistical study of the evolution tempo, it is not possible to do a detailed study of the number of objects that were lost in the previous Gyr, since the detailed efficiency depends on the detailed configuration of a past forgotten time (be it 1 Gyr into the past, or at the beginning of the solar system).

While we have made use of complex simulations to explore the future evolution of the outer solar system, it is arguably worth considering if a simpler (and faster) model is able to provide us with the conclusions about the future stability of the system reached so far. This can be accomplished by treating our DP sample as mass-less test particles, recognizing that the orbital dynamics are dominated by the giant planets.

We ran 200 additional simulations with identical parameters and initial conditions as those of Section 4.1, but considering only the gravitational perturbations of the four giant planets over the 34 DPs, which were treated as test particles. Each simulation differs from the others by the initial time step of the integrator, varied from 350 to 450 days on 0.5 day intervals. The simplified simulations were approximately 4 times faster that our full N-body simulations.

An analysis over the 200 simplified simulations shows results that are in general equivalent to the full N-body runs. There were a total 2221 irregular orbits and 372 ejections in the simplified runs, as opposed to 2195 irregular orbits and 378 ejections in the full runs. Although, globally, these results are interchangeable, when looking at the behaviour of individual objects there are some irreconcilable differences.

First, out of the 22 objects that have at least 1 irregular orbit in table 3, 10 present statistically different results within the simplified simulations. When comparing the full runs to the simplified runs we find that 3 DPs are more stable, 6 DPs are less stable, and 1 DP has more irregular orbits but less ejections. Two objects in particular stand out, namely Haumea and Orcus. For Haumea, the simplified runs produce 146 ejections versus 125 in the full runs; for Orcus, not a single ejection is obtained in the simplified runs, additionally it remains well locked inside the 3:2 MMR in all cases; in comparison the full runs produce 26 ejections and only 161 regular orbits in $a$. The latter result highlights the importance of self interactions between massive DPs, particularly inside MMRs, where mutual distances between the objects can be significantly reduced and secular perturbations accumulate more efficiently. This shows that the effect of DPs mutual perturbations on the secular evolution of the Kuiper belt cannot be in general ignored, otherwise some of the long-term conclusions will be flawed.

Finally, besides the statistical analysis of the whole sample, we consider that the evident chaotic nature of some of the DPs (e.g. Haumea, $2007\,{\rm OR}_{10}$, or $2015\,{\rm RR}_{245}$) calls for an individualized analysis of their expected evolution; we intend to pursue such analysis in a future paper (Muñoz-Gutiérrez et al. in preparation).