Inclination pathways of planet-crossing asteroids

Close encounters of small bodies and planets are one of the major mechanisms that drive the orbital evolution of early disc planetesimals and current asteroid populations in the solar system such as comets, Centaurs and scattered trans-neptunian objects (TNOs). Planet-crossing asteroids in the outer solar system originated from two older reservoirs. One is the early planetesimal disc a small part of which was scattered beyond Neptune’s orbit to the Oort cloud region in the first Gyr or so after planet formation (Levison & Duncan, 1997; Tiscareno & Malhotra, 2003; Emel’yanenko et al., 2005; Di Sisto & Brunini, 2007; Brasser et al., 2012b; Nesvorný et al., 2017; Fernández et al., 2018; Vokrouhlický et al., 2019; Kaib et al., 2019). The other reservoir is the interstellar medium where the sun captured small bodies present in its birth cluster of stars (Fernández & Brunini, 2000; Levison et al., 2010; Brasser et al., 2012a; Jílková et al., 2016; Hands et al., 2019; Kaib et al., 2019; Namouni & Morais, 2018b, 2020c; Portegies Zwart et al., 2021). After a few Gyr residence in the outer solar system during which they experience the effects of the sun’s Galactic environment, asteroids and comets move back on planet-crossing orbits towards the inner solar system.

Centaurs, in particular, are known to have a chaotic dynamical evolution driven by resonance capture and close encounters with the planets (Tiscareno & Malhotra, 2003; Bailey & Malhotra, 2009; Volk & Malhotra, 2013; Morais & Namouni, 2013; Namouni & Morais, 2015; Morais & Namouni, 2016, 2017; Namouni & Morais, 2017, 2018a). Recently, long-term time-backward numerical simulations of high-inclination Centaurs have shown that their orbits follow distinct pathways in the inclination-semimajor axis plane on time-scales unrelated to their dynamical time (Namouni & Morais, 2018b, 2020c) (hereafter Papers I and II). It was found that 19 high-inclination Centaur orbits clustered around polar inclinations beyond Neptune’s orbit 4.5 Gyr in the past. The inclination clustering is the basis for the identification of the interstellar origin of high-inclination Centaurs because trans-neptunian space was empty of solar system material with polar inclinations with respect to the ecliptic $4.5$ Gyr in the past.

Clustering around polar inclinations in trans-neptunian space was first seen in the discovery article of TNO (471325) 2011 KT19 where Chen et al. (2016) simulated the time-forward evolution of $10^{3}$ TNO-clones over 1 Gyr. No attempt was made to investigate the inclination clustering’s dynamical origin probably because the object is currently on a nearly-polar orbit and the polar structure in the clone simulation could have simply reflected the stability of the TNO’s inclination. A first indication of the dynamical origin of the polar clustering was given in the million clone simulation of Jupiter’s co-orbital asteroid Ka‘epaoka‘āwela (Paper I) where the clone swarm was found to follow the asteroid’s Tisserand relation, a result later confirmed in (Köhne & Batygin, 2020). That region of polar motion in trans-neptunian space was termed the polar corridor and identified as a prime location for finding sun-bound asteroids of interstellar origin captured in the early solar system (Paper I). In (Namouni & Morais, 2020a) (hereafter Paper III), the Tisserand relation was further linked to the dynamics of high-inclination Centaurs 2016 YB13 and 1999 LE31 and the time-backward evolution of TNO (471325) 2011 KT19. The clustering of high-inclination Centaurs end-states in the polar corridor was also thought to be an artefact of simulating Centaur motion back in time because of a violation of the second law of thermodynamics (Morbidelli et al., 2020). In Paper III, it was explained using physical principles that no such violation occurred. However, the analytical foundation of the dynamical origin of the polar clustering, whether forward or backward in time, is still lacking.

The Tisserand relation is a time-independent function of the asteroid’s eccentricity and inclination, and the semimajor axis ratio of the asteroid and the planet. It was derived from the exact conservation of the Jacobi constant of the three-body problem (Jacobi, 1836) to identify returning comets whose orbital elements had changed following their encounters with the planets (Tisserand, 1896). The Tisserand relation was used to characterize small body encounter geometries with Jupiter (Kresák, 1980; Carusi et al., 1995), and to distinguish the orbital types of comets (ecliptic, Halley-type, Jupiter-family) (Levison & Duncan, 1994) and TNOs whose main perturber is Neptune (Elliot et al., 2005). If the Tisserand relation is somehow responsible for the inclination clustering of high-inclination Centaurs in a simulation with four giant planets and can indicate the inclination pathways of such planet-crossing asteroids over 4.5 Gyr, then this result must stem from the three-body problem where the consequences of the conservation of the Tisserand relation can be derived analytically, and the dynamical origin of the polar clustering can be established. This is the aim of this work.

Regardless of their origins and inclinations, we study analytically and numerically the long-term evolution of Centaurs whose main distinguishing feature we retain is their crossing of the planets’ orbits. In effect, whether they originate in the planetesimal disc or the interstellar medium, Centaurs were drawn back to the planets’ domain from trans-neptunian space through the process of close encounters with the planets.

In section 2, we recall the derivation of the Tisserand relation in the circular restricted three-body problem and apply it to planet-crossing asteroids. In section 3, we show analytically how the Tisserand relation defines a unique portrait of the planet’s influence on astroids crossing its orbit through specific inclination pathways that depend on the Tisserand parameter. The Tisserand relation of planet-crossing asteroids partitions parameter space into four dynamical regions. Asteroid injection inside the planet’s orbit onto retrograde and high-inclination prograde orbits occurs only in one region that is similar to and contains the polar corridor. In section 4, we determine analytically the secular dispersion of eccentricity and inclination across the inclination pathway caused by the planet’s secular perturbations through the Kozai-Lidov potential. This helps us confirm the precise agreement of the analytical theory with the numerical simulations and show how the Tisserand parameter still constrains the asteroid’s motion even after it temporarily leaves its Tisserand inclination pathway. In section 5, we use high-resolution statistical simulations on Gyr time-scales to confirm that the Tisserand inclination pathways describe accurately the long-term dynamical evolution of planet-crossing asteroids forward as well as backward in time because the Tisserand relation and the Kozai-Lidov potential are time-independent. We also show using asteroid orbits with maximal eccentricity that the partitioning of the inclination semimajor axis parameter space by the Tisserand relation is that which corresponds to encounters with the planet at perihelion or aphelion regardless of whether the asteroid’s initial perihelion or aphelion is different from the planet’s semimajor axis. In section 6, we illustrate asteroid motion in the Tisserand inclination pathways. In section 7, we discuss our findings and their implications for high-inclination Centaurs, comet delivery from the Oort cloud, extreme TNOs, and the presence of unknown planets in trans-neptunian space.

The equations of motion of the circular restricted three-body problem in the inertial centre-of-mass frame are written as:

where $M_{\odot}$ and $M_{p}$ are the masses of the star and the planet ($M_{p}\ll M_{\odot}$), and ${x}$, ${x}_{\odot}$ and ${x}_{p}$ are the positions of the asteroid, the sun and the planet respectively. As the asteroid is assumed to be massless, the orbits of the sun and planet are circular. Their rotation vector perpendicular to their orbital plane is denoted by ${n}$. The dynamical system has a constant of motion, the Jacobi constant (Jacobi, 1836), which is written as (Murray & Dermott, 2000):

After the planet-crossing asteroid moves far away from the planet, its motion around the Sun is purely Keplerian. The first two terms in Eq. (1) then equal the opposite of the asteroid’s Keplerian orbital energy $GM_{\odot}/2a$, where $a$ is the orbital semimajor axis. The third term vanishes and the fourth is the product of the asteroid’s vertical component of its specific angular momentum, $[GM_{\odot}a(1-e^{2})]^{1/2}\cos I$, and the planet’s mean motion, $[GM_{\odot}a_{p}^{-3}]^{1/2}$ where $e$ and $I$ are the asteroid’s orbital eccentricity and inclination respectively and $a_{p}$ the planet’s semimajor axis. Rearranging the constants $GM_{p}$ and $a_{p}$ in $C_{J}$ leads to the Tisserand relation:

The asteroid’s encounter with the planet is brief implying that the new post-encounter values of $a$, $e$ and $I$ away from the planet conserve the Tisserand relation (Tisserand, 1896). Close encounters with a planet tend to conserve the perihelion, $q=a(1-e)$, or the aphelion $Q=a(1+e)$ depending on the encounter’s geometry (Duncan et al., 1987). The perihelion and aphelion of a planet-crossing asteroid satisfy $0\leq q\leq a_{p}$ and $a_{p}\leq Q\leq 2a_{p}$. The strongest perturbations occur when $q\sim a_{p}$ or $Q\sim a_{p}$ and lead to the fastest evolution of the asteroid’s orbit.

In terms of the constant $q$ and $Q$, the Tisserand relation may be written as:

Three features of the Tisserand relation should be emphasized. First, the Tisserand relation is not an exactly conserved quantity along the asteroid’s orbit as it is based on taking the limit of the Jacobi constant away from the planet. Secondly, the Tisserand relation of a planet-crossing asteroid given by (3,4) is based on the conservation of perihelion (aphelion) distance of a planet-crossing asteroid. Processes such as the secular perturbations of the planet and its mean motion resonances can modify $q$ and $Q$, so that the asteroid’s orbit may no longer satisfy (3,4). We show in section 4 that the effect of secular perturbations is to cause inclination and eccentricity dispersions that can be accounted for analytically as the asteroid temporarily moves away from the orbit-crossing conditions. When mean motion resonances affect the orbit-crossing conditions, they lead to a radial drift away from the Tisserand inclination pathway (3,4). In section 5 we show that this occurs only for small inclinations and large Tisserand parameters ($T>2.7$) whereas secular perturbations dominate at moderate to high inclinations. Thirdly, the Tisserand relation defines the inclination pathway of a planet-crossing asteroid but does not indicate its direction of motion along the pathway. Examples of asteroid motion in the inclination pathways forward and backward in time are given in section 6.

A planet-crossing asteroid with initial orbital elements, $a_{0}$, $q$ and $I_{0}$, has a Tisserand parameter:

To derive the pathway the asteroid follows in the inclination-semimajor axis plane, the Tisserand relation (3) is inverted to obtain:

For asteroids inside the planet’s orbit with a constant aphelion, the inclination pathway’s expression is identical to (6) after substituting $Q$ for $q$. The inclination pathway $I(a,T,q$ or $Q,a_{p})$ is shown in Fig. (1) as a function of $a/a_{p}$ and various values of the Tisserand parameter for three values of $q$ and $Q$. In the left-hand panel, the asteroid’s initial orbit may be outside the planet’s orbit with $q=a_{p}$ ($a\geq a_{p}$) or inside it with $Q=a_{p}$ ($a\leq a_{p}$). In the middle panel, the asteroid is initially outside the planet’s orbit with $q=0.6\,a_{p}$. In the right-hand panel, the asteroid is initially inside the planet’s orbit with $Q=1.5\,a_{p}$.

To avoid cumbersome notation in describing the mathematical properties of (6), we define the normalized perihelion or aphelion distance as $\bar{q}=q/a_{p}$ or $Q/a_{p}$ because the expression of $I(a,T,q,a_{p})$ is identical to $I(a,T,Q,a_{p})$. From the perihelion’s and aphelion’s definitions, $\bar{q}=q/a_{p}$ in the range $0\leq\bar{q}\leq 1$ and $\bar{q}=Q/a_{p}$ in the range $1\leq\bar{q}\leq 2$.

The inclination pathway $I(a,T,\bar{q},a_{p})$ of a planet-crossing asteroid is defined only for $a/a_{p}\geq\bar{q}/2$ reflecting the fact that eccentricity is $\leq 1$.¹¹1$\bar{q}$ here stands for $Q/a_{p}$. $I(a,T,\bar{q},a_{p})$ has an asymptote far from the planet given by:

for $-\sqrt{8\bar{q}}\leq T\leq\sqrt{8\bar{q}}$. The inclination pathway reaches a maximum given by:

for the following values of the Tisserand parameter: $2/\bar{q}\leq T\leq\bar{q}^{2}+2/\bar{q}$ for $0\leq\bar{q}\leq 2^{1/3}$ and $2/\bar{q}\leq T\leq 4/\bar{q}$ for $2^{1/3}\leq\bar{q}\leq 2$. The inclination pathway leads the asteroid’s orbit to coplanarity before it is reflected back on the same pathway. The smallest reflection semimajor axis is given as:

defined in the ranges $-\sqrt{8\bar{q}}\leq T\leq\bar{q}^{2}+2/\bar{q}$ for $0\leq\bar{q}\leq 2^{1/3}$ and $-\sqrt{8\bar{q}}\leq T\leq\sqrt{8\bar{q}}$ for $2^{1/3}\leq\bar{q}\leq 2$. The second reflection semimajor axis is given by:

and exists only in the range $\sqrt{8\bar{q}}\leq T\leq\bar{q}^{2}+2/\bar{q}$ and $0\leq\bar{q}\leq 2^{1/3}$.

In this work, we seek the analytical foundation of the long-term simulation results of high-inclination Centaurs. Such objects feel the strongest perturbations from the planets and evolve onto Neptune-crossing orbits with $q\sim a_{p}$ and $Q\sim a_{p}$ (see Figure 2 of Paper II). We therefore specialize in planet-crossing asteroids with $q=a_{p}$ and $Q=a_{p}$. In section 4, we examine what happens to the inclination pathways when the conditions $q=a_{p}$ and $Q=a_{p}$ are no longer satisfied. In section 5, we show using asteroids with maximal eccentricity orbits that the physically-relevant partitioning of the inclination-semimajor axis parameter space by the Tisserand relation is that with $q=a_{p}$ and $Q=a_{p}$ regardless of whether the initial perihelion and aphelion distances are different from $a_{p}$.

For asteroids scattered at perihelion or aphelion, The Tisserand inclination pathway is written as:

The Tisserand relation partitions the inclination-semimajor axis plane into four regions (left-hand panel of Figure 1). Region I is that which includes only retrograde motion (region above the solid blue curve). The corresponding Tisserand parameter $-\sqrt{8}\leq T\leq-1$. In this region, the orbit’s inclination $I(a,T,a_{p})$ has an asymptote:

A far-away asteroid inclined by $I_{\infty}$ moving towards the planet cannot enter inside its orbit as $a>a_{p}$. The smallest possible semimajor axis that can be reached is the reflection semimajor axis $a_{180^{\circ}}$:

In the vicinity $a_{180^{\circ}}$, the incoming asteroid nears retrograde coplanarity ($I=180^{\circ}$) before it is reflected back on the same Tisserand inclination pathway from where it came. At the lower boundary of Region I given by $T=-1$, $a_{180^{\circ}}=a_{p}$ and $I_{\infty,1}\simeq 110^{\circ}$.

Region II is defined by $-1\leq T\leq 2$ (region between the solid and dashed blue curves). In this region, asteroids have moderate prograde to polar inclinations far from the planet and may enter inside its orbit. If an asteroid follows the inclination pathway towards its reflection semimajor axis $a_{180^{\circ}}$, its motion is forced to become retrograde. In the vicinity of $a_{180^{\circ}}$, the incoming asteroid is reflected back on the same Tisserand inclination pathway and may move away from the planet. No other region allows such evolution. Asteroids starting with aphelia at the planet may move away from its orbit towards the inclination $I_{\infty}$. For the lower boundary of Region II, $T=2$ (dashed blue curve) yields $a_{180^{\circ}}=a_{p}/2$ and $I_{\infty,2}=45^{\circ}$. Region II is a high-inclination corridor that is reminiscent of the polar corridor of high-inclination Centaurs of Papers I–III that it contains.

Region III is defined by $2\leq T\leq\sqrt{8}$ and motion is prograde everywhere (region between the dashed blue and dotted red curves). Starting from $I_{\infty}$ away from the planet, inclination increases along the pathway as the asteroid gets closer to the planet reaching a maximum at (8). If the asteroid moves closer to the planet, its orbit tends towards the sun-planet orbital plane $I=0^{\circ}$ at the reflection semimajor axis $a_{0^{\circ},1}=a_{180^{\circ}}$ (13) where it is sent back on the same Tisserand inclination pathway. The upper boundary of Region III has $T=2$ (dashed blue curve) which yields $a_{m}=0.5\,a_{p}$, and $I_{m}=90^{\circ}$. For the lower boundary of Region III, $T=\sqrt{8}$ (dotted red curve) yields $I_{\infty}=0^{\circ}$, $a_{m}\simeq 0.85\,a_{p}$, and $I_{m}\simeq 24^{\circ}$.

Region IV is defined by $\sqrt{8}\leq T\leq 3$ (below the dotted red curve). The Tisserand relation indicates that motion may be radially confined to a region defined by $a_{0^{\circ},1}\leq a\leq a_{0^{\circ},2}=a_{0^{\circ}}$ (10) where inclination reaches its maximum at (8). If the asteroid approaches one of the two edges of Region IV, its orbit becomes coplanar with the planet’s.

The constant perihelion (aphelion) distance may not always be satisfied. If the semimajor axis’ change after successive encounters of the asteroid and planet is small, the planet’s pull will slowly modify the eccentricity and inclination of the asteroid and may temporarily invalidate the orbit-crossing conditions upon which the inclination pathway is based. This aspect is examined analytically in the next section and numerically in section 5. The orbit-crossing conditions can also be modified by mean motion resonances with the planet. The strongest mean motion resonances are close to the planet and will affect motion in Region IV where the semimajor confinement at small inclination and large Tisserand parameters will not occur (section 5).

When the initial orbit-crossing condition differs from $\bar{q}=1$, the regions defined by the Tisserand relation are modified except one (see Figure 1). The no-injection Region I exists for all values of $\bar{q}$. The lower inclination boundary of Region I is given by the equality $a_{\rm 180^{\circ}}=a_{p}$ yielding $T=1-2\sqrt{\bar{q}(2-\bar{q})}$ and :

The lower inclination boundary of Region I is shown in Figure (2). $I_{\infty}^{\rm noinject}\geq 90^{\circ}$ for most perihelion distances and reaches a maximum of $112.5^{\circ}$ at $\bar{q}=0.689$. Therefore regardless of perihelion distance and the planet’s position, planet-crossing asteroids cannot be injected inside the planet’s orbit if their inclination far from the planet $I>I_{\infty}^{\rm noinject}$. Instead, they are reflected near the semimajor axis $a_{\rm 180^{\circ}}$ (9) and sent back to the outer solar system.

The change to the asteroid’s orbital elements following a planet encounter is random and may be large or small. When the semimajor axis change is small over several encounters, the secular perturbations from the planet will build up giving rise to secular periodic changes in the orbital eccentricity and inclination. This implies that in the secular oscillation cycle the eccentricity will move away from the value that satisfies the orbit-crossing conditions, $q=a_{p}$ or $Q=a_{p}$, only to return to it at the cycle’s end. This in turn implies that the orbit will temporarily leave the inclination pathway (11) as the latter is based on the orbit-crossing conditions $q=a_{p}$ or $Q=a_{p}$. The asteroid’s orbit eventually leaves the secular cycles it experiences at constant semimajor axis, and returns to the inclination pathway (11) when the amplitude of the semimajor axis change after a planet encounter is large. The secular oscillations at constant semimajor axis will appear as an inclination dispersion across the inclination pathway but also as a dispersion in the eccentricity-semimajor axis plane with respect to the orbit-crossing conditions.

It is possible to find the amplitude of the eccentricity and inclination dispersions analytically. In the circular restricted three-body problem, the asteroid’s orbit secular evolution is driven by the Kozai-Lidov potential (Kozai, 1962; Lidov, 1962; Quinn et al., 1990):

where $K$ is the complete elliptic integral of the first kind, $r=a(1-e^{2})/(1+e\cos f)$ is the radial position of the asteroid and $f$ its true anomaly. The Kozai-Lidov potential conserves the projection of the asteroid’s orbital angular momentum on the rotation vector ${n}$ of the sun-planet system, $L_{z}=\sqrt{a(1-e^{2})}\cos I$ (section 2). Figure 3 shows the level curves of the Kozai-Lidov potential. In each row, the inclination-argument of pericentre $(I,\omega)$ level curves are shown for a prograde inclination $I_{0}$ and its polar symmetric retrograde inclination $180^{\circ}-I_{0}$. As the Kozai-Lidov potential has reflection symmetry with respect to exactly polar motion, the two inclination panels are symmetrical with respect to $I=90^{\circ}$ and the corresponding $(e,\omega)$ level curves are identical. The dashed curves are the three-dimensional orbit-crossing conditions $a_{p}=a(1-e^{2})/(1\pm e\cos\omega)$. The Kozai-Lidov level curves show that for an initially prograde inclination $I_{0}$ and eccentricity $e_{0}$ on the orbit-crossing curve at perihelion ($\omega=0^{\circ}$ and $a_{p}=a(1-e)$), the inclination travels on a secular cycle that at maximum amplitude makes the asteroid coplanar with the planet. Similarly, a retrograde orbit will travel towards retrograde coplanarity that is a maximum inclination of $180^{\circ}$ before returning to the initial inclination and eccentricity values.

An asteroid in the inclination pathway with a Tisserand parameter $T$ whose orbital semimajor axis is no longer strongly perturbed by the planet will leave the pathway with an initial inclination $I(a,T,a_{p})$ (11) and a vertical angular momentum component given by (2) as:

The inclination dispersion by the Kozai-Lidov potential will therefore occur in the range $[0^{\circ},I(a,T,a_{p})]$ for prograde orbits ($I\leq 90^{\circ}$) and in the range $[I(a,T,a_{p}),180^{\circ}]$ for retrograde orbits ($I\geq 90^{\circ}$).

Coplanarity (i.e. $I=0^{\circ}$ or $I=180^{\circ}$) is achieved for the largest possible eccentricity value. This means that the maximal eccentricity that can be achieved is given by the coplanarity relation $\pm\sqrt{a(1-e_{\rm max}^{2})}=L_{z}$ because $L_{z}$ is conserved as soon as the asteroid moves away from the inclination pathway onto the Kozai-Lidov secular cycle. The plus and minus signs indicate prograde and retrograde motion respectively. This relationship and (16) yield the eccentricity dispersion amplitude:

Like the vertical angular momentum component, the eccentricity dispersion does not depend explicitly on the perihelion or aphelion distance but only implicitly through the value of the Tisserand parameter.

To the maximal eccentricity corresponds a minimal perihelion distance $q_{\rm min}=a(1-e_{\rm max})$. As the asteroid moves away from the planet $a\rightarrow\infty$ and $e_{\rm max}\rightarrow 1$ but the minimal perihelion distance is finite and tends to:

The eccentricity dispersion should be viewed alongside the orbit-crossing conditions $q=a_{p}$ and $Q=a_{p}$ as these bound the eccentricity from below. Examples are shown in Fig. 4. In the no-injection Region I, eccentricity dispersion is small and confined to the vicinity of the orbit-crossing condition (panel $T=-2$). There is a minimum eccentricity $e_{l}$ below which the asteroid’s orbit cannot go given by $e_{\rm max}(a_{180^{\circ}},T,a_{p})$ at the reflection semimajor axis $a_{180^{\circ}}$. For $T=-2$, $e_{l}=0.53$. Dispersion increases as Region II is approached ($T=-1$). Inside Regions II and III the maximal eccentricity may reach 1 for small semimajor axes (panels $T=1$ and $2$). The limit eccentricity $e_{l}$ steadily decreases to zero as $T$ reaches 1 then increases steadily until $T=2$ only to decease afterwards. Near the lower boundary of Region III, eccentricity dispersion starts to decrease with increasing Tisserand parameter (panels $T=2$ and $2.85$) to ultimately be confined around the planet’s semimajor axis.

It should be noted that since the Kozai-Lidov secular potential (15) and the Tisserand inclination pathway (11) are both time-independent, time-forward asteroid motion along and across the inclination pathway will be similar to its time-backward motion. This is confirmed using statistical simulations in the next section.

To ascertain the validity of the analytical developments in the previous sections, we perform long-term high-resolution statistical simulations of the full equations of motion of the three-body problem with Neptune as the perturbing planet located at $a_{p}=30.1$ au on a circular orbit. The asteroid has initial inclination $I_{0}$ and semimajor axis $a_{0}$. section 5.1 is devoted to asteroids initially located far outside the planet’s orbit with a perihelion distance $q_{0}=a_{p}$. In section 5.2 we examine asteroids inside the planet’s orbit with an initial aphelion distance $Q_{0}=a_{p}$. In section 5.3, the initial perihelion and aphelion are chosen using the maximal eccentricity curve (17). The initial mean longitude relative to Neptune is set to $\lambda=180^{\circ}$ and the initial argument of perihelion and longitude of ascending node are set to $\omega=0^{\circ}$ and $\Omega=0^{\circ}$. For each initial inclination and semimajor axis, the asteroid’s initial orbit is cloned $10^{4}$ times using a multivariate normal distribution with standard deviations in its equinoctial orbital elements of eccentricity, inclination and mean longitude of $10^{-6}$. For the semimajor axis, a relative standard deviation $\sigma_{a}=6\times 10^{-7}a_{0}$ (au) was chosen so that the simulation outcomes of the different initial semimajor axes $a_{0}$ may be compared to one another particularly in terms of the number of surviving clones. The clone swarm’s orbital elements’ deviations from the nominal orbit generate a dispersion in the initial Tisserand parameter, $T_{0}$, of $\sim 10^{-6}$ to $10^{-5}$ (Table I). The orbital elements’ standard deviations are similar to the observational error bars of high-inclination Centaurs’ orbits (see Paper II).

The equations of motion were integrated for 1 Gyr forward and 1  Gyr backward in time for asteroids that encounter the planet at perihelion or aphelion, and for 2 Gyr forward and 2  Gyr backward in time for asteroids with maximal eccentricities on account of their slower dynamics. The integration uses the Bulirsch and Stoer algorithm with an error tolerance of $10^{-11}$. Clone orbital evolution stops when one of the following events occurs: collision with the Sun, collision with Neptune, ejection from the system, reaching the inner 1 au semimajor axis boundary and reaching the outer boundary at $10^{4}$ au. The choice of the outer boundary in the inner Oort cloud region is motivated by showing the inclination pathway at a large distance from the planet. The simulation does not describe the evolution of asteroids in that region as the gravitational effect of the Galactic environment was not modelled.

Each set of initial conditions may be characterized by the initial nominal Tisserand parameter, $T_{0}$ (5), its standard deviation $\sigma_{T_{0}}$, the inclination at infinity, $I_{\infty}$ (12) and the reflection semimajor axis $a_{180^{\circ}}$ (13). These are given in Table I. The simulation outcome may be characterized by the mean Tisserand parameter at the simulation’s end, $\langle T\rangle$, its standard deviation $\sigma_{T}$ and the inclination standard deviation from the analytical Tisserand pathway inclination (11). The latter is given for two semimajor axis’ ranges: the first is [1:$10^{4}$ au] which, as will become apparent below, is equivalent to [$a_{180^{\circ}}$:$10^{4}$ au] and gives the global standard deviation. We also use the range [$3\,a_{180^{\circ}}$:$10^{4}$ au] that indicates the inclination’s standard deviation from the analytical pathway inclination away from the planet. The factor $3$ is related to the fact that regardless of inclination, the smallest reflection semimajor axis $a_{180^{\circ}}=a_{p}/2$ (section 3). This puts $3\,a_{180^{\circ}}$, but not $2\,a_{180^{\circ}}$, outside the planet’s neighbourhood and beyond its strongest perturbations. Another reason for using the second semimajor axis range is that the high-inclination Centaur simulations (Papers I and II) showed that these objects move away from the planet on polar orbits. The inclination’s dispersion observed in those simulations concern semimajor axes larger than 50 au which fall in the asymptotic part of the inclination pathway (11). The two inclination standard deviations are denoted by $\sigma_{I}$ and $\sigma_{I,a\gg a_{p}}$. The simulation outcomes are given in Table I. Unlike our previous works on the origin of high-inclination Centaurs (Papers I and II), we do not quote the collision and ejection fractions of unstable clones because of the presence of the simulation’s outer boundary at $10^{4}$ au and the absence of three giant planets. We only quote the number of stable clones to indicate the orbit’s stability in this work’s setting.

The exact orbit of a planet-crossing asteroid cannot be known on Gyr time-scales. The asteroid’s long-term evolution forward and backward in time can be ascertained only statistically. It is still instructive to examine examples of clone motion in order to see the different dynamical processes present as well as the orbits’ diversity in the inclination pathway. The examples will also illustrate how the conserved Tisserand parameter conditions the long-term time-forward and time-backward dynamical evolution of planet-crossing asteroids whether they are in the inclination pathway or leave it temporarily because of secular perturbations.

Figure 9 contains orbit examples for different initial conditions and different time directions. In order to give an idea about the time residency in parameter space, the orbits are shown in purple for the first 0.1 Gyr and in green for the subsequent 0.9 Gyr for which a smaller data sampling is used so as not to overcrowd the Figure.

The first three examples in Panels (a), (b) and (c) illustrate motion in the polar corridor (Papers I–III). Panel (a) shows an asteroid clone with initial nominal parameters $I_{0}=170^{\circ}$, $a_{0}=20$ au and an aphelion at the planet’s orbit moving backward in time to a polar orbit of semimajor axis $452$ au. The clone initially moves along the inclination pathway (11) and across it through secular oscillations seen as the vertical segments indicating semimajor axis’ conservation. These encounters take the clone to the vicinity of the planet and consequently to nearly circular motion. The next planet encounters send the asteroid back on its tracks before expelling it with a single kick away from the planet to 200 au where it travels back and forth on the inclination pathway. This occurs in the first 0.1 Gyr. For the subsequent 0.9 Gyr, the gravitational kicks are small enough so that the clone moves smoothly along the inclination pathway to 452 au.

Panel (b) shows an asteroid clone with the same initial nominal parameters in Panel (a) moving backward in time to a polar orbit of semimajor axis 633 au. Unlike the previous example, the clone sticks around its initial orbit for the first 0.1 Gyr then exhibits high-eccentricity secular oscillations bounded by $e_{\rm max}$ (17). This results in large inclination oscillations caused by the Kozai-Lidov perturbation and reminiscent of the cloud around the planet’s orbit in Fig. 7 ($I_{0}=170^{\circ}$). During these cycles, the clone moves steadily outside the planet’s neighbourhood until it is removed at high eccentricity by a single kick onto a polar orbit where it makes its way to its final position.

Panel (c) shows an incoming polar asteroid clone located at 300 au with a perihelion at the planet’s orbit. The clone makes it way forward in time to an orbit inclined by $170^{\circ}$ of semimajor axis 20 au similar to the initial orbit of the previous two examples. The clone follows the inclination pathway and the orbit-crossing conditions all the way to its final position passing through a nearly-circular orbit in the vicinity of the planet’s location.

Panel (d) shows an incoming asteroid clone located initially at 300 au with a perihelion at the planet’s orbit moving forward in time and being reflected by the planet in the no-injection Region I. Its initial nominal inclination is $130^{\circ}$ and its final semimajor axis and inclination are respectively 2390 au and $127^{\circ}$. In the first 0.1 Gyr, the clone is transferred near a semimajor axis of 106 au where it exhibits Kozai-Lidov oscillations. After reaching a minimum semimajor axis of 59 au, the clone is sent back on its way far from the planet. The eccentricity evolution is strictly confined between the perihelion-encounter condition and $e_{\rm max}$ (17).

Panel (e) shows an initially polar asteroid clone with a nominal semimajor axis of 20 au and an aphelion at the planet’s orbit moving backward in time to a semimajor axis of 913 au and a prograde inclination of $60^{\circ}$. The clone travels along the inclination pathway and the orbit-crossing conditions. In the first 0.1 Gyr, it reaches a maximum semimajor axis 368 au. During the next 0.9 Gyr, it achieves a minimum and maximum semimajor axes of 49 au and 1691 au respectively before arriving at its final semimajor axis.

Panel (f) shows an incoming asteroid clone located at 300 au with an initial nominal inclination of $60^{\circ}$ and a perihelion at the planet’s orbit moving forward in time to reach a polar orbit at 21 au at 1 Gyr similar to the initial orbit of Panel (e) albeit with a larger eccentricity. In the first 0.1 Gyr, its orbit jumps to 50 au and explores motion at small eccentricity in the planet’s neighbourhood. In the next 0.9 Gyr, the clone travels on the inclination pathway to a maximum semimajor axis of 182 au, turns around and moves near the orbit-crossing condition towards its final position.

The last two panels show examples of motion at maximal eccentricity forward and backward in time and further illustrate the conservation of the Tisserand parameter. Panel (g) shows an asteroid clone moving forward in time from an initial nominal inclination of $60^{\circ}$, semimajor axis of 300 au and a perihelion at the planet’s orbit. The clone follows the orbit-crossing conditions and the inclination pathway in the first 0.1 Gyr to nearly-circular motion. In the planet’s coorbital region, the eccentricity is raised almost to $e_{\rm max}$ (17). As the clone’s motion is reflected by the planet’s gravitational pull and its semi-axis increases, the inclination is lowered to $4^{\circ}$ at 1 Gyr corresponding to a semimajor axis of 145 au whereas the eccentricity follows precisely the eccentricity dispersion curve defined by the asteroid’s Tisserand parameter. In time the clone will return to its original Tisserand inclination pathway.

Panel (h) shows a returning orbit moving backward in time from a low inclination and a Tisserand parameter equal to that in Panel (g). The asteroid clone starts out with initial nominal inclination $0.46^{\circ}$, semimajor axis $a_{0}=100$ au and perihelion distance $q=5.612$ au on the maximal eccentricity curve (17). In the first 0.1 Gyr it travels back and forth with respect to its initial location on the maximal eccentricity curve. As the clone reaches the planet’s vicinity its eccentricity is lowered by the planet’s gravitational kicks as well as the secular perturbations (vertical segments) until the asteroid reaches the inclination pathway of its Tisserand parameter (11) that corresponds to the constant perihelion $q=30.1$ au as explained in section 5.3. The subsequent motion occurs on this pathway where the clone travels to a maximum distance of 332 au and ends up at 221 au with an inclination $\sim 60^{\circ}$ at 1 Gyr near the initial orbit of the example of Panel (g). The clone’s motion is unrelated to the inclination pathway of its initial perihelion shown in Panel (h) as the dotted orange curve.

The examples in this section showed how the conserved Tisserand parameter determines the long-term evolution of planet-crossing asteroids forward and backward in time. The examples had Tisserand parameters $T<2.7$ where the effect of mean motion resonances is negligible. They illustrated how the Tisserand inclination pathway is followed, and how maximal eccentricity low inclination planet-crossing orbits may evolve into large inclination orbits. In all our simulations of low inclination asteroid orbits with $T>2.7$ near the lower boundary of Region III and in Region IV, no clone achieved large prograde or retrograde inclinations.

In this work, we set out to explain analytically the dynamical origin of the polar clustering found in the time-backward statistical simulations of high inclination Centaurs. To do so, we reduced the dynamical system where high-inclination Centaurs evolve to its minimum: a planet-crossing asteroid in the three-body problem. In this simplified setting, the planet randomly imparts gravitational kicks to the asteroid’s orbit. Their amplitudes may be large or small, and may affect the different orbital elements differently. The orbit may also be influenced by mean motion resonances and secular perturbations. All of this makes the long-term orbital evolution of a planet-crossing asteroid chaotic. In this work, we demonstrated analytically and numerically that for moderate- to high-inclination asteroids, there is some order within that chaos. Order is ensured by the Tisserand relation that defines the inclination pathway the planet-crossing asteroid follows (11), and the Kozai-Lidov secular potential that explains the inclination and eccentricity dispersions across the inclination pathway (15,17). The Tisserand parameter, in particular, is inscribed on the planet-crossing asteroid’s orbit and governs its long-term evolution on Gyr-time-scales whether time flows forwards or backwards.

The Tisserand inclination pathway depends only on the Tisserand parameter regardless of the asteroid’s initial perihelion or aphelion distance. This was illustrated for asteroids with maximal eccentricity and hence the smallest perihelion and largest aphelion possible (section 5.3). It was found that the long-term evolution of such asteroids ignores the inclination pathways that correspond to their initial perihelia and aphelia and follows the Tisserand pathways of $q=a_{p}$ and $Q=a_{p}$ (11). It can be checked numerically that this result holds for intermediate perihelion and aphelion distances between the orbit crossing conditions and the maximal eccentricity curve.

The partitioning of the inclination-semimajor axis parameter space by the Tisserand relation is reminiscent of the partitioning of motion space by the Jacobi constant and its zero velocity curves (Murray & Dermott, 2000). The latter applies to any kind of motion and gives the regions of permissible and forbidden motion depending on the value of the Jacobi constant. The former applies only to planet-crossing asteroids and gives their inclination pathways. The partitioning of parameter space (left-hand panel of Figure 1) unveiled several novel aspects of the dynamics of planet-crossing asteroids. One is orbital reflection at the semimajor axis $a_{180^{\circ}}$ (13). Like the Tisserand inclination pathway (11), the reflection semimajor axis does not depend on the initial perihelion or aphelion.

Although its evolution is chaotic and the asteroid may travel towards or away from the planet (Figure 9), in the long-term it will be reflected by the planet at or outside of $a_{180^{\circ}}$. Even for small inclinations and $T>2.7$, leakage below the reflection radius was minimal and confined to the immediate vicinity of the analytical value $a_{180^{\circ}}$ indicating that modelling the effect of mean motion resonances will only displace $a_{180^{\circ}}$ slightly.

Another novel result is the existence of the no-injection Region I. For Tisserand parameters $T\leq-1$ (or equivalently inclinations at infinity in the range $110^{\circ}\leq I_{\infty}\leq 180^{\circ}$), $a_{180^{\circ}}$ is larger the planet’s semimajor axis implying that planet-crossing asteroids will be reflected outside the planet’s orbit and cannot become Centaurs or short-period comets. This finding will have consequences for the delivery of comets from the spherically shaped Oort cloud as the existence of the no-injection semimajor axis will discriminate between the processes that make comets cross the planets’ orbits and the order in which they occur. For instance, Neptune reflects back all comets and asteroids with Tisserand parameters $T_{\rm Neptune}\leq-1$ (or inclinations $I_{\infty}>110^{\circ}$). Jupiter the innermost giant planet will do the same and reflect comets with $T_{\rm Jupiter}\leq-1$ but those objects with semimajor axes between Jupiter and Neptune will still enter the giant planets’ domain and may become Centaurs and short-period comets. This gives a first limit on the Tisserand parameter of comets driven by close encounters with Jupiter. Setting $a_{180^{\circ}}(a_{\rm Jupiter})=a_{\rm Neptune}$ where $a_{\rm Jupiter}$ and $a_{\rm Neptune}$ are the planets’ semimajor axes, leads to the Tisserand parameter cutoff $T_{\rm Jupiter}\leq-2.53$ (or $I_{\infty}\geq 153^{\circ}$) below which no comet from the Oort cloud may enter the giant planets’ domain. As of 2021, September 3, the JPL Small-Body Database Search Engine³³3https://ssd.jpl.nasa.gov indicates that short-period comets (with periods $<200$  yr) have Tisserand parameters $T_{\rm Jupiter}>-1.65$ in accordance with the previous cutoff. More interestingly, of the listed 609 long period comets with eccentricities $<1$ (and periods $\geq 200$ yr), only 5 have Tisserand parameters $T_{\rm Jupiter}<-2.53$. This suggests that the large majority of the known long period comets will be or have been able to enter the giant planets’ domain in that their semimajor axes could have been or will become smaller than Neptune’s.

The reflection semimajor axis $a_{\rm 180^{\circ}}$ and the no-injection Region I may be used to constrain the existence of unknown planets in the outer solar system (Trujillo & Sheppard, 2014; Batygin & Brown, 2016; Malhotra et al., 2016; Volk & Malhotra, 2017; Silsbee & Tremaine, 2018; Oldroyd & Trujillo, 2021). The hypothetical planets are usually assumed to have eccentric and inclined orbits so that they may explain possible orbital alignments of the dynamically-detached TNOs. Inclined planets beyond Neptune would be able to inject Oort cloud comets only with inclinations $0\leq I_{\infty}\leq 110^{\circ}$ relative to their orbital plane which, depending on their inclination would in turn affect the inclination distribution of comets that reach inside Neptune’s orbit. Furthermore, the hypothetical planets are assumed to have semimajor axes $100\alpha$ au from the Sun with $\alpha$ of order a few. Therefore TNOs that cross an outer planet’s orbit cannot travel further inwards than the semimajor axis $50\alpha$ (see section 3). The presence of such a TNO population with an inner semimajor axis boundary could be indicative of the planet’s position. In particular, it would be useful to examine if the known extreme TNOs form a component of such a population whose dynamics is described by maximal eccentricity motion (section 5.3). The analysis developed in this paper for a planet on a circular orbit will need to be extended to eccentric planets.

Going back to the original scope of this article, the polar clustering of high inclination Centaurs can be explained by the analytical and numerical results in this work.

In Paper II, a search for the past stable orbits of 19 high inclination Centaurs was performed by integrating the equations of motion with all giant planets and the Galactic tide 4.5 Gyr in the past using $\sim 2\times 10^{7}$ clones. The logic of that search was that since Centaurs were believed to have originated from the early planetesimal disc, their orbits are 4.5 Gyr old. As the gravitational equations of motion are time-reversible, it possible to search statistically for the region of parameter space such objects came from. If no clear dynamical pathways were followed by those Centaurs, then the statistical search would be inconclusive and past orbits would be scattered all over parameter space. Instead, high inclination Centaur orbits clustered mainly in the polar corridor with small inclination dispersions. In their time-backward evolution high-inclination Centaurs crossed the orbits of the giant planets as they moved outward away from Jupiter’s reflection semimajor axis. In the last Gyrs of evolution their clones crossed Neptune’s orbit (Paper II, Figure 2). From an initial inclination range of [62^∘:173^∘], the 19 high-inclination Centaurs evolved to an average inclination range of $\sim$[60^∘:90^∘] with a dispersion $\sim 10^{\circ}$ in the scattered disc region at $-4.5$ Gyr (Paper II Figure 4). Interpreting the latter as a range of inclinations at infinity (12) gives the Tisserand parameter range with respect to Neptune as $0\leq T\leq 1.5$ which is enlarged to $-0.5\leq T\leq 1.9$ when the inclination dispersion is taken into account. These ranges correspond to Neptune’s Region II (Figure I, left-hand panel). The small inclination dispersion comes from the polar symmetry of the Kozai-Lidov potential not only of Neptune but of all giant planets because the relative inclinations of the giant planets’ orbital planes are negligible.

A recurring question about the simulations of high-inclination Centaurs is whether the inclination pathways they follow in their time-backward evolution are related to those they follow in a time-forward simulation that starts following their motion far away from the planets as they arrive from the Oort cloud early in solar system history. The simulations in Figures 5 and 7 answer this question. Consider the planet-crossing asteroid in Figure 7 (top row, second panel from left) located at $a_{0}=20$ au and inclined by $I_{0}=170^{\circ}$ with a Tisserand parameter $T=0.12$. The asteroid’s evolution is propagated backward in time. After 1 Gyr the clone swarm disperses along the Tisserand inclination pathway away from the planet. This asteroid is representative of Centaurs with retrograde motion such as those examined in Paper II. The swarm’s inclination clusters around $88^{\circ}\pm 3^{\circ}$ (Table I). In order to determine if such dispersed clones away from the planet can actually follow the same inclination pathway when their motion is reversed in time, we may simulate the time-forward evolution of a clone swarm with an inclination of $90^{\circ}$ with a semimajor axis far from the planet. This is done in Figure 5 (bottom row, left-hand panel) and Figure 6 (top row) for $a_{0}=300$ au. The clone swarm disperses as predicted on the same Tisserand inclination pathway with inclination dispersions similar to those of the inner retrograde asteroid. Examples of individual clones of these two simulations are given in Figure 9, panels (a), (b) and (c). Further corresponding examples on either side of the planet with similar Tisserand parameters are given in Figures 7 and 5: ($I_{0}=90^{\circ}$, $a_{0}=20$ au, negative time) $\leftrightarrow$ ($I_{0}=60^{\circ}$, $a_{0}=300$ au, positive time), and ($I_{0}=60^{\circ}$, $a_{0}=20$ au, negative time) $\leftrightarrow$ ($I_{0}=40^{\circ}$, $a_{0}=300$ au, positive time).

To understand more precisely how the past orbits of high-inclination Centaurs evolve in time, it is necessary to examine their motion in the presence of all the giant planets. This is beyond the scope of this paper but it is possible to identify how the current results derived in the context of the three-body problem will be modified. First, the Centaur’s Tisserand parameter with any planet will not be conserved as it is a purely three-body invariant. The evolution of the Tisserand parameters with the giant planets on the dynamical time-scale occurs as the Centaur’s clone swarm disperses and moves towards Neptune. The Centaur’s inclination pathway will depend on the evolution of its perihelion as it too is no longer conserved but will secularly move towards Neptune’s semimajor axis. When the clone swarm’s mean perihelion has reached Neptune, the Tisserand parameter becomes constant in the scattered disc region and the Tisserand inclination pathways apply. In the inner Oort could region, however, the Galactic environment modifies the Centaur’s eccentricity and inclination and consequently the inclination pathway. The second way this work’s results are modified concerns the collision and ejection rates. Centaurs are more stable in the three-body problem than they are in the solar system. The asteroid in Figure 7 with an initial inclination of $170^{\circ}$ has a median lifetime $>1$ Gyr as it still had 7611 clones at $+1$ Gyr out of the initial $10^{4}$. In Paper II the two Centaurs with a similar inclination (330759) 2008 SO218 and (434620) 2005 VD have median lifetimes $\sim 2$ Myr. The presence of all the planets therefore does not change the inclination pathways in the scattered disc but by decreasing the Centaurs’ stability it modifies the determination of the original population’s size with respect to a three-body treatment.

The analysis of the agreement between the analytical and numerical results in this work was possible because of the high-resolution statistical simulation that concentrates a large number of clones in a minuscule neighbourhood of the asteroid’s orbit. This approach was initiated with the million clone simulation of Jupiter’s co-orbital asteroid (Paper I). Its importance as illustrated in this work goes beyond the chaotic dynamics of high-inclination Centaurs to the evolution of the early planetesimal disc. In the early solar system, planets are thought to have undergone an unstable dynamical phase where Neptune is pushed into the early planetesimal disc whose planetesimals it scatters from the immediate vicinity of its orbit to the Oort cloud (Pfalzner et al., 2015). The example in Figure 6 (bottom row) illustrates what happens to a low-inclination asteroid scattered by Neptune. After an evolution over 1 Gyr, the $10^{4}$-clone swarm initially inclined by $10^{\circ}$ and confined to a range of $10^{-4}$ au produces a disc that extends from 20  au to $10^{4}$ au (our arbitrary outer simulation boundary in the absence of Galactic effects). It is interesting to note the distribution similarity of our Figure 6 with Figure 6 Panels (a) and (b) in (Nesvorný et al., 2017) where the authors study the dynamics of short period comets from the effect of the planets’ instability on a 10 au-wide disc of a million planetesimals. Despite the inclination amplitude difference and the presence of four planets in the short-period comet simulation, the inclination similarity likely indicates the presence of inclination pathways the disc planetesimals follow when excited by Neptune’s perturbations. Analytical modelling of the effect of mean motion resonances on the Tisserand relation for small inclinations may help uncover the pathways followed by disc planetesimals in the early solar system.

I thank the reviewer for constructive comments that helped improve the clarity of the paper. The numerical simulations were done at the Mésocentre SIGAMM hosted at the Observatoire de la Côte d’Azur.

The data underlying this article will be shared on reasonable request to the author.