# Orbit injection of planet-crossing asteroids

Fathi Namouni

► To cite this version:

Fathi Namouni. Orbit injection of planet-crossing asteroids. Monthly Notices of the Royal Astronomical Society, inPress, 527 (3), pp.4889-4898. 10.1093/mnras/stad3570. hal-04290818

HAL Id: hal-04290818

https://hal.science/hal-04290818v1

Submitted on 8 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Orbit injection of planet-crossing asteroids

F. Namouni

Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, CS 34229, F-06304 Nice, France

Accepted 2023 November 15. Received 2023 November 3; in original form 2023 August 25

## ABSTRACT

Solar system Centaurs originate in trans-Neptunian space from where planet orbit crossing events inject their orbits inside the giant planets' domain. Here, we examine this injection process in the three-body problem by studying the orbital evolution of trans-Neptunian asteroids located at Neptune's collision singularity as a function of the Tisserand invariant, $T$. Two injection modes are found, one for $T > 0.1$, or equivalently prograde inclinations far from the planet, where unstable motion dominates injection, and another for $T \le 0.1$, or equivalently polar and retrograde inclinations far from the planet, where stable motion dominates injection. The injection modes are independent of the initial semimajor axis and the dynamical time at the collision singularity. The simulations uncovered a region in the polar corridor where the dynamical time exceeds the Solar system's age suggesting the possibility of long-lived primordial polar trans-Neptunian reservoirs that supply Centaurs to the giant planets' domain.

**Key words:** celestial mechanics – comets: general – Kuiper belt: general – minor planets, asteroids: general.

## 1 INTRODUCTION

Understanding precisely how the dynamical transfer of trans-Neptunian objects (TNOs) into the giant planets' domain takes place is key to identifying the origin of the small body populations in the outer Solar system (Fernández, Gallardo & Young 2016; Fraser et al. 2022; Kaib & Volk 2022).

In this regard, progress was recently spurred by the long-term statistical simulations of high-inclination Centaurs that demonstrated that their likely past orbits clustered perpendicular to the Solar system's invariable plane in a region that extends beyond Neptune's orbit termed 'the polar corridor' (Namouni & Morais 2018, 2020a,b).

The search for the polar corridor's dynamical origin and its role in the transfer of TNOs onto high-inclination Centaur orbits over Gyr time-scales was the object of our earlier work (Namouni 2022, hereafter Paper I). By studying analytically and numerically the long-term dynamics of a planet-crossing asteroid in the three-body problem, we demonstrated that despite the highly chaotic nature of the asteroid–planet interactions, the Tisserand invariant is conserved and that for moderate to high inclinations, this invariance makes the asteroid follow the Tisserand inclination pathway given as

$$I(a, T, a_p) = \arccos \left( \left[ T - \frac{a_p}{a} \right] \left[ 4 \left( 2 - \frac{a_p}{a} \right) \right]^{-\frac{1}{2}} \right), \quad (1)$$

where $a$ and $a_p$ are the semimajor axes of the asteroid's and the planet's orbits, $T$ is the Tisserand invariant,

$$T = \frac{a_p}{a} + 2 \left[ \frac{a(1 - e^2)}{a_p} \right]^{\frac{1}{2}} \cos I, \quad (2)$$

and $e$ and $I$ are the asteroid orbit's eccentricity and inclination. It was also demonstrated that the inclination pathway (equation 1) is independent of the initial perihelion and aphelion.

The motion of planet-crossing asteroids follows the Tisserand inclination pathway whether time flows forwards or backwards on the time-scale of Gyr. The pathway is followed strictly when the gravitational kicks imparted to the asteroid's orbit induce significant change to its semimajor axis. When the orbit's semimajor axis is not strongly perturbed, the asteroid's orbit wanders away from the inclination pathway at constant Tisserand invariant because of the build-up of secular perturbations, only to return to the Tisserand inclination pathway.

Paper I demonstrated orbit injection of incoming asteroids, from outside to inside the planet's orbit, onto polar or retrograde inclination orbits occurs only for $-1 < T < 2$. Translated in terms of the Tisserand pathway's inclination far from the planet,

$$I_{\infty}(T) = \arccos(T/\sqrt{8}), \quad (3)$$

the previous statement becomes: only TNOs with inclinations $45^{\circ} < I_{\infty} < 110^{\circ}$ may enter inside the planet's orbit and assume inclinations $I \ge 90^{\circ}$. TNOs with $I_{\infty} < 45^{\circ}$ may be injected inside Neptune's orbit but can only assume prograde orbits, whereas those with $I_{\infty} > 110^{\circ}$ may not be injected at all. The long-term statistical simulations of high-inclination Centaurs (Namouni & Morais 2020b, hereafter Paper II) utilized $2 \times 10^7$ clones to follow the Centaurs' evolution back to $-4.5$ Gyr from initial inclinations in the range $[62^{\circ}:173^{\circ}]$. They yielded average inclinations for the Centaurs' end states away from Neptune in the range $[60^{\circ}:90^{\circ}]$ with a standard deviation $\sim 10^{\circ}$ indicating that the polar corridor corresponds approximately to the analytically determined range of $-1 < T < 2$.

Paper I established a novel dynamical behaviour of planet-crossing asteroids. The decreasing semimajor axis of an incoming asteroid, whose planet encounters move it towards the planet, cannot be smaller than a minimum value that depends on its Tisserand invariant.

<footnote>* E-mail: namouni@oca.eu</footnote>
The asteroid's motion is reflected back at this minimum value, termed 'the reflection semimajor axis', $a_{\text{refl}}$, meaning that the asteroid's semimajor axis increases after approaching $a_{\text{refl}}$ that is given as

$$a_p a_{\text{refl}}^{-1}(T) = T - 2 + 2(3 - T)^{\frac{1}{2}}. \quad (4)$$

The smallest reflection semimajor axis occurs for $T = 2$ at $a_{\text{refl}} = a_p/2$, whereas for $T < -1$, $a_{\text{refl}} > a_p$ and reflection occurs outside the planet's orbit thus precluding orbit injection. The validity of the Tisserand inclination pathway and the reflection semimajor axis holds for $T \lesssim 2.7$ or $I_{\infty} > 17^{\circ}$ even though the Tisserand invariant is always conserved regardless of its value. This limitation is related to the effect of strong mean motion resonances at low inclinations (see Paper I).

This work is the continuation of the analysis in Paper I. We examine precisely how Neptune-crossing asteroids are injected from outside to inside its orbit. Studying the dynamical pathways numerically in Paper I utilized the cloning of an asteroid's orbit with a perihelion (aphelion) equal to or smaller (larger) than Neptune's semimajor axis. Although this approach is adequate to establish the conservation of the Tisserand invariant and study the dynamics of the Tisserand inclination pathways, it does not provide a comprehensive analysis of the injection process over a large set of diverse initial conditions. In this work, we follow a different approach by studying the dynamics at Neptune's collision singularity.

In the three-body problem, the collision singularity is defined by the orbital intersection of the Neptune-crossing asteroid with the planet. For an initial Tisserand invariant and an asteroid semimajor axis, a diverse range of eccentricities, inclinations, and orbital angles can be generated to cover the collision singularity, and followed over time spans larger than the Solar system's age. As the Tisserand invariant is conserved, the diverse asteroid orbits that cover the collision singularity should in principle follow the same Tisserand inclination pathway and allow us to examine the orbit injection process as a function of the Tisserand invariant.

In Section 2, the initial conditions set-up at the collision singularity is explained. In Section 3, we examine the dynamics at the collision singularity for stable and unstable orbits and find that the Tisserand invariant conservation at the collision singularity ensures that the statistical ensemble's initially diverse orbits follow the same Tisserand inclination pathway. The dynamical time at the collision singularity is measured and found to exceed the Solar system's age for certain negative Tisserand invariant values or equivalently retrograde inclinations at infinity. The evolution of the minimum semimajor axis at the collision singularity is examined and reveals a semimajor axis-independent two-mode injection process, one at retrograde inclinations at infinity where stable motion dominates injection and another at prograde inclinations at infinity where unstable motion dominates injection. We also examine the role of mean motion resonance capture in the collision singularity and find that it cannot explain overall stability and the two-mode injection process. In Section 4, we measure the extent of the collision singularity's stable region in the polar corridor and identify the possible location of Centaur supplying small body reservoirs beyond Neptune's orbit. In Section 5, we sum up our findings and point out how the presence of long-lived primordial polar TNO reservoirs was already indicated in the earlier numerical simulations with all giant planets in Paper II.

## 2 INITIAL CONDITIONS AT THE COLLISION SINGULARITY

The collision singularity is defined by the intersection of the three-dimensional orbit of the planet-crossing asteroid with the planet's circular orbit. Intersection may occur at the ascending or descending nodes and is expressed by equating the asteroid's radial distance at its nodes with the planet's semimajor axis as follows:

$$r = \frac{a(1 - e^2)}{1 \pm \cos \omega} = a_p, \quad (5)$$

where $a$ and $a_p$ are the asteroid's and the planet's semimajor axes, and $e$ and $\omega$ are the asteroid's eccentricity and argument of perihelion. The $\pm$ signs correspond to the ascending or descending nodes, respectively, with $0^{\circ} \le \omega \le 180^{\circ}$.

For an initial Tisserand invariant $T$ and an initial semimajor axis $a$, we generate random initial conditions that place the asteroid on an intersecting trajectory as follows. The asteroid's mean longitude $\lambda$, longitude of ascending node $\Omega$, and $\omega$ are drawn from uniform distributions in the intervals $[-180^{\circ}:180^{\circ}]$ for the first two and $[0^{\circ}:180^{\circ}]$ for the third angle. Eccentricity is obtained from the nodal distance equation (5) to ensure orbit intersection as

$$e_{\text{coll}} = \frac{\left[ 4a(a - a_p) + a_p^2 \cos^2 \omega \right]^{\frac{1}{2}} \pm a_p \cos \omega}{2a}, \quad (6)$$

where the $-$ (+) sign corresponds to a crossing at the ascending (descending) node. The orbit's inclination is obtained from the Tisserand invariant (equation 2) using the initial semimajor axis and the collision eccentricity (equation 6) as follows:

$$I = \arccos \left( (T - a_p/a) \left[ 4a(1 - e_{\text{coll}}^2)/a_p \right]^{-\frac{1}{2}} \right). \quad (7)$$

The minimum eccentricity of the collision singularity ensemble generated by this method is that of an encounter at perihelion $e_{\min} = 1 - a_p/a$. The maximum eccentricity is determined by the Tisserand invariant as the amplitude of the argument in the brackets of equation (7) should be $\le 1$. This leads to the maximum eccentricity encountered in Paper I and given by

$$e_{\max} = \left( 1 - \frac{a_p}{4a} \left[ T - \frac{a_p}{a} \right]^2 \right)^{\frac{1}{2}}. \quad (8)$$

The existence of a maximum eccentricity reduces the argument of perihelion range depending on the value of $T$. Fig. 1 shows the collision eccentricity for encounters at the ascending and descending nodes for an initial semimajor axis 200 au that intersects Neptune's orbit with $a_p = 30.1$ au. Maximal eccentricity is shown for $T = 1.5, 2.5$, and $2.15$. The latter value is the limit below which there is no gap in the argument of perihelion's range. It is given generally as $T = 2 + a_p a^{-1}$. The $\omega$ gap is shown in the figure as the shaded region for $T = 2.5$.

The inclination range is related to the eccentricity range by the Tisserand invariant. At minimum eccentricity, the inclination is that of the Tisserand inclination pathway (equation 1), whereas at maximum eccentricity, the inclination is $0^{\circ}$ for prograde orbits and $180^{\circ}$ for retrograde orbits.

## 3 COLLISION SINGULARITY DYNAMICS

The Tisserand invariant conservation was established in Paper I by generating clones of a number of planet-crossing asteroid orbits. For one such asteroid with a Tisserand invariant $T_0$, the cloned orbits differ slightly from the nominal orbit and the deviations of their
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Collision eccentricity (equation 6) as a function of argument of perihelion for ascending (solid blue) and descending (solid dark blue) nodes with an asteroid at *a* = 200 au. Maximal eccentricity (equation 8) is shown as a horizontal dashed line for three values of the Tisserand invariant. The shaded region is the *ω* gap below *T* = 2.5.0</figcaption></figure>

<table><caption>**Table 1.** Initial eccentricity and inclination ranges [$e_{min}:e_{max}$] and [$I_{min}:I_{max}$]. $e_{min}$ is the eccentricity for perihelion crossing, $e_{max}$ is given by equation (8), $I_{min}$ and $I_{max}$ are given by equation (1) and either of 0° (prograde motion) or 180° (retrograde motion). $I_{∞}$ is given by equation (3).</caption><thead><tr><th>$T_0$</th><th>$e_{max}$</th><th>$I_{min}$ (°)</th><th>$I_{max}$ (°)</th><th>$I_{\infty}$ (°)</th></tr><tr><td></td><td colspan="4">$a_0$ = 100 au, $e_{min}$ = 0.699</td></tr></thead><tbody><tr><td>-0.7</td><td>0.961</td><td>113</td><td>180</td><td>104</td></tr><tr><td>-0.1</td><td>0.993</td><td>99</td><td>180</td><td>92</td></tr><tr><td>1.2</td><td>0.969</td><td>0</td><td>70</td><td>65</td></tr><tr><td>2.7</td><td>0.753</td><td>0</td><td>23</td><td>17</td></tr><tr><td></td><td colspan="4">$a_0$ = 200 au, $e_{min}$ = 0.849</td></tr><tr><td>-0.8</td><td>0.983</td><td>104</td><td>180</td><td>106</td></tr><tr><td>-0.3</td><td>0.996</td><td>99</td><td>180</td><td>96</td></tr><tr><td>1.0</td><td>0.986</td><td>0</td><td>72</td><td>70</td></tr><tr><td>2.5</td><td>0.980</td><td>0</td><td>30</td><td>28</td></tr><tr><td></td><td colspan="4">$a_0$ = 73 au, $e_{min}$ = 0.588</td></tr><tr><td>-0.5</td><td>0.956</td><td>111</td><td>180</td><td>100</td></tr><tr><td>0</td><td>0.991</td><td>99</td><td>180</td><td>90</td></tr><tr><td>1.4</td><td>0.948</td><td>0</td><td>67</td><td>60</td></tr><tr><td>2.8</td><td>0.642</td><td>0</td><td>19</td><td>8</td></tr></tbody></table>

Tisserand invariants were measured with respect to the nominal orbit. In this work, the initial conditions for an initial $T_0$ are different as there is no nominal orbit. By letting the initial angle variables sample their maximum possible range subject to the limits set by the Tisserand invariant $T_0$, the generated statistical ensemble with an initial semimajor axis $a_0$ is made up of different asteroids with diverse eccentricities and inclinations and not clones of the same asteroid.

The equations of motion of the circular restricted three-body problem with Neptune at $a_p = 30.1$ au and collision singularity asteroids at $a_0 = 73, 100, \text{and } 200$ au were integrated forward in time using the Bulirsch–Stoer algorithm with an error tolerance of $10^{-11}$. In this simulation, the Tisserand invariant takes values in the range $-1.2 \le T_0 \le 2.8$ with a 0.1 step. The choice of this range is motivated by two particular values of the Tisserand invariant. One is $T_0 = -1$ the upper inclination boundary of the polar corridor below which asteroids are reflected back outside the planet's orbit $a_{refl} > a_p$ and injection is not possible. The other particular value is $T_0 = 2.7$ above

which mean motion resonances induce a significant departure of the inclination pathway from $I_∞$ even though the Tisserand invariant is conserved (see Paper I).

For each $T_0$, 1000 orbits were generated. The eccentricity and inclination ranges for each semimajor axis $a_0$ produced by the method in Section 2 are given in Table 1. The integration time span of 5.4 Gyr is larger than the Solar system's age because the simulation unveiled a region in parameter space where the dynamical time exceeds 4.5 Gyr. As in Paper I, orbital evolution stops when the asteroid is ejected, collides with the Sun or Neptune, or reaches an outer semimajor axis boundary at $10^4$ au. For each asteroid orbit, the minimum semimajor axis is monitored.

## 3.1 Tisserand inclination pathways

Fig. 2 shows snapshots of the orbit distribution of the collision singularity ensembles with $a_0 = 100$ au for four values of the Tisserand invariant $T_0 = 2.7, 1.2, -0.1, \text{and } -0.7$. These examples illustrate the dynamics corresponding to $I_∞$ from small and moderate prograde to polar and retrograde motion. Stable orbits are shown in green and unstable orbits in red.<fnref n="1" /> Also shown are the Tisserand inclination pathway (equation 1), the maximum eccentricity curve (equation 8), and the reflection semimajor axis (equation 4).

The top row panels show the eccentricity and inclination as a function of semimajor axis at 4.5 Gyr for stable orbits and at the last sampling epoch before instability for unstable orbits. The conservation of the Tisserand invariant makes the diverse ensemble of collision singularity asteroids follow the unique Tisserand inclination pathway given by their initial $T_0$ as predicted by the analytical theory of Paper I. Their eccentricity and inclination dispersions are explained by the von Zeipel–Lidov–Kozai secular potential (von Zeipel 1910; Kozai 1962; Lidov 1962) as evidenced by maximum eccentricity curve of each $T_0$ (see Paper I for an analytical discussion of this effect). Orbits are predictably reflected at $a_{refl}(T_0)$. These conclusions are valid for stable and unstable orbits as the latter are simply former stable orbits for which an ultimate gravitational kick knocked them out of their stable motion by ejecting them from the system or forcing them to collide with the Sun or Neptune. Instability is stronger as $T_0$ tends to 3.

The location $a_0 = 100$ au is close to the 1:6 outer resonance ($a_{res} = 99$ au) where some asteroids are captured mostly at large inclinations ($T_0 = 1.2, -0.1, \text{and } -0.7$) as mean motion resonance capture is particularly efficient at large inclinations (Namouni & Morais 2015, 2017). Some collision singularity asteroids with retrograde inclinations become unstable while at resonance. Resonance capture is examined in Section 3.4 as a function of the Tisserand invariant.

The panels in the second row from the top show the eccentricity and inclination at minimum semimajor axis, $a_{min}$. The asteroids follow their Tisserand inclination pathway in ($a_{min}, I$) space just as they do in ($a, I$) space and are bounded by $e_{max}$ in ($a_{min}, e$) space. Stable and unstable asteroids alike fill the possible area in ($a_{min}, I$) and ($a_{min}, e$) showing that not all asteroids reach $a_{refl}(T_0)$ or are injected inside the planet's orbit.

The third row panels show the epoch of minimum semimajor axis, $T_{a_{min}}$. This epoch is bounded from below by a power-like function $\log(\tau_{min}/\tau_0) = (a_0 - a)^\alpha$, where $a$ and $a_0$ are expressed in au, $\alpha$

<footnote n="1">¹Stable orbits live the full duration of the simulation (≥4.5 Gyr). Unstable orbits either collide with the Sun or the planet or are ejected from the system. Their lifetimes range from a few 1000 yr to <4.5 Gyr.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Distribution of collision singularity ensemble orbits at 4.5 Gyr for different initial Tisserand invariants $T_0$ and the initial semimajor axis $a_0 = 100$ au. First row from top: distributions in the $(a, I)$ and $(a, e)$ planes of stable (green) and unstable (red) orbits. The latter are shown at the last sampling epoch before instability. The inclination pathway (equation 1) is the solid curve in the inclination panels. In the eccentricity panels, the top curve is the maximum eccentricity dispersion (equation 8) and the bottom ones are the perihelion and aphelion collision conditions $a(1 \pm e) = a_p$. The solid and dashed vertical lines indicate the reflection semimajor axis $a_{\text{refl}}$ (equation 4) and the planet's position, respectively. Second row: distributions in the $(a_{\text{min}}, I)$ and $(a_{\text{min}}, e)$ planes. Third row: epoch at minimum semimajor axis as a function of $a_{\text{min}}$. Bottom row: initial conditions in the mean longitude argument of perihelion plane. The blue squares indicate injection was achieved.</figcaption></figure>

~ 0.5 and $\tau_0 \sim 10^3-10^5$ yr with decreasing $T_0$. The traveltime to the minimum semimajor axis is longer for retrograde orbits because planet encounters generate weaker perturbations as they occur at larger speeds in briefer times than prograde orbits. Asteroids captured in resonance have some of the largest $T_{a\text{min}}$. Resonant asteroid statistics are given in Section 3.4.

The bottom row panels show the initial conditions' distribution in the mean longitude, argument of perihelion plane for stable and unstable orbits. Asteroids that achieved orbit injection are shown with blue squares. An $\omega$ gap is present for $T_0 = 2.7$. There are no preferred regions or initial condition clustering with respect to stability or orbit injection.

In Section 3.3, we examine the dynamical time in the collision singularity and show that there is a Tisserand invariant range where the dynamical time exceeds the Solar system's age for $a_0 = 100$ au. We therefore set up two further high-resolution simulations, one with $a_0 = 200$ au where the dynamical time nowhere exceeds 4.5 Gyr and another with $a_0 = 73$ au so that the initial location is not close to a major mean motion resonance like 100 au is to the 1:6 resonance. The scope of these two simulations is to find out if the dynamics at different semimajor axes differ qualitatively from that of $a_0 = 100$ au because of stability and proximity to resonance.

The simulation results with $a_0 = 200$ au are shown in Fig. 3 for $T_0 = 2.5, 1, 0.3$, and $-0.8$. They confirm that the colli-
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Same as Fig. 2 but with $a_0 = 200$ au.</figcaption></figure>

sion singularity asteroids follow the Tisserand inclination pathways in the $(a, I)$ and $(a_{\min}, I)$ planes for stable and unstable orbits. Their eccentricities are bounded by the maximal eccentricity (equation 8). The location $a_0 = 200$ au is close to the 1:17 mean motion resonance ($a_{res} = 199$ au) where similar to the previous simulation some asteroids are captured mostly at large inclinations ($T_0 = 1$ and $-0.3$) whereas some collision singularity asteroids become unstable while at resonance. The dynamics at $a_0 = 200$ au is qualitatively similar to that of $a_0 = 100$ au.

The value $a_0 = 73$ au was chosen instead of the rounded off 70 au to avoid the 2:7 resonance. Instead $a_0 = 73$ au lies within the 9:34 resonance of order 25. The simulation results are shown in Fig. 4 for $T_0 = 2.8, 1.4, 0$, and $-0.5$. $T_0 = 2.8$ is below the validity limit of the Tisserand inclination pathways $T = 2.7$ (see Paper I). It is shown to illustrate as in Paper I that reflection still occurs even if $T_0 > 2.7$ but the actual reflection semimajor axis becomes slightly smaller

than equation (4) (for details see Paper I). The results in Fig. 4 are qualitatively similar to those of the previous two semimajor axes. Inclination pathways and maximum eccentricity curves are followed precisely for $T_0 \le 2.7$ and there is no preferred orbit geometry that maximizes orbit injection (bottom panels in Figs 2–4). Some asteroids are captured in the 9:34 resonance (see Section 3.4).

We may conclude that regardless of the semimajor axis and proximity to mean motion resonance, the collision singularity asteroids are found to follow the Tisserand inclination pathways<fnref n="2" /> and the maximal eccentricity curves associated with their Tisserand invariant. Their orbits fill the available parameter space in the $(a_{\min}$,

<footnote n="2">²The Tisserand inclination pathway is followed in the sense explained in the Introduction either strictly if semimajor axis changes are large or as an inclination boundary if secular perturbations build up. See examples in Fig. 10.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same as Fig. 2 but with $a_0 = 73$ au.</figcaption></figure>

l) and (a<sub>min</sub>, e) planes. There are no preferred orbital angles that maximize orbit injection.

## 3.2 Tisserand invariant conservation

Deviation from the initial value of the Tisserand invariant is measured at 4.5 Gyr and shown in Fig. 5. The solid lines show the deviation of the mean Tisserand invariant $\langle T \rangle$ from its initial value $T_0$ for stable and unstable orbits as a function of $T_0$ for $a_0 = 73, 100,$ and $200$ au. The deviation $|\langle T \rangle - T_0|$ is of order a few $10^{-3}$ and is smaller for unstable orbits than stable orbits. The reason is twofold. First, unstable orbits are more numerous as can be seen in Fig. 6 where the ratio of the number of stable orbits $N_s$ to that of unstable orbits $N_u$ is shown as a function of $T_0$ (and $I_\infty$). Second, the accumulation of numerical error is larger for stable orbits that survived up to 4.5 Gyr, whereas unstable orbits live over smaller time spans from $\sim 10^3$ yr to

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Conservation of the Tisserand invariant at 4.5 Gyr. Solid lines indicate the mean value of *T* with respect to *T*₀, and dashed lines indicate the Tisserand invariant’s standard deviation σ_T. Green (red) indicates stable (unstable) orbits.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Ratio of the number of stable orbits present at 4.5 Gyr, $N_s$, to the number of all unstable orbits from the simulation's start up to 4.5 Gyr, $N_u$, as a function of the initial Tisserand invariant (and $I_\infty$). The solid curves correspond, respectively, to $a_0 = 200$ au (red), 100 au (blue), and 73 au (green).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Dynamical time $\tau_d$ as a function of the initial Tisserand invariant (and $I_\infty$) for $a_0 = 200$ au (red), 100 au (blue), and 73 au (green).</figcaption></figure>

less than 4.5 Gyr. This explains why $|\langle T \rangle - T_0|$ increases as $T_0$ tends to -1 where stable orbits become more numerous.

The standard deviation $\sigma_T$ of $T$ at 4.5 Gyr is shown in Fig. 5 with dashed lines for stable and unstable motion. It is of order a few $10^{-3}$ across all values of $T_0$. This deviation amplitude and that of $|\langle T \rangle - T_0|$ do not affect quantitatively the dynamics of the Tisserand inclination pathways as can be seen in Figs 2–4.

## 3.3 Dynamical times

For an initial semimajor axis $a_0$ and Tisserand invariant $T_0$, we define the singularity collision dynamical time, $\tau_d(a_0, T_0)$, as the median lifetime of the statistical ensemble generated by the method in Section 2. This definition is inspired by the long-term evolution of the ensemble's diverse orbits that cover the available parameter space defined by $a_0$ and $T_0$ (Figs 2–4) much like the clones of single asteroid orbits in Paper I. Fig. 7 shows the dynamical time as a function of $T_0$ (and $I_\infty$) for $a_0 = 73, 100$, and 200 au.

Regardless of the semimajor axis, the dynamical time has a distinct shape with three peaks: (i) a prominent retrograde-inclination peak with $-0.85 \le T_0 \le -0.6$ equivalent to an inclination at infinity $107^\circ \ge I_\infty \ge 102^\circ$; (ii) an intermediate prograde-inclination peak with 0.7

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Ratio of the resonant stable orbits, $N_{res}$, to the total number of stable orbits, $N_s$, at epoch $t$ as a function of the Tisserand invariant. The epochs are $t = 0.1$ Gyr (dashed black), 1 Gyr (solid green), 3 Gyr (solid blue), and 4.5 Gyr (solid red). Relevant resonances are shown next to the initial semimajor axis $a_0$.</figcaption></figure>

$\le T_0 \le 1$ ($76^\circ \ge I_\infty \ge 69^\circ$); and (iii) a shallow prograde-inclination peak with $2.1 \le T_0 \le 2.3$ ($42^\circ \ge I_\infty \ge 36^\circ$).

The dynamical time is of order 0.1–0.4 Gyr near the shallow peak. It increases to 1.5–2 Gyr at the intermediate peak. The retrograde-inclination peak occurs at $\tau_d = 4.11$ Gyr and $T_0 = -0.9$ for $a_0 = 200$ au. It exceeds the age of the Solar system for $a_0 = 73$ au and $a_0 = 100$ au with, respectively, $\tau_d = 5.37$ Gyr at $T_0 = -0.5$ and $\tau_d = 5.20$ Gyr at $T_0 = -0.8$. Regardless of the semimajor axis, the dynamical time drops sharply from a few Gyr to ~1 Gyr for polar orbits ($T_0 = 0$).

The region where the dynamical time exceeds 4.5 Gyr is $-1.1 \le T_0 \le -0.3$ for $a_0 = 73$ au and $-1 \le T_0 \le -0.4$ for $a_0 = 100$ au. This raises the possibility of long-lived nearly polar TNOs that cross Neptune's orbit. Injection is examined in Section 3.5 and the stable region's semimajor axis extent is estimated in Section 4.

## 3.4 Resonance capture

The simulations show that resonance capture may occur for collision singularity asteroids. It is possible that the large dynamical times found at retrograde inclinations that exceed the Solar system's age are the result of such temporary resonance capture as asteroid orbits are protected from the planet's destabilizing perturbations. We examined quantitatively the simulations of the three semimajor axes and found that capture occurs mainly in the mean motion resonance nearest to the initial semimajor axis. These are 9:34 for $a_0 = 73$ au, 1:6 for $a_0 = 100$ au, and 1:17 for $a_0 = 200$ au. Fig. 8 shows for each initial semimajor axis, the ratio of the number of stable asteroid orbits in mean motion resonance with Neptune, $N_{res}$, to the number of stable asteroid orbits, $N_s$, as a function of the initial Tisserand invariant $T_0$ at epochs 0.1, 1, 3, and 4.5 Gyr.

Regardless of the semimajor axis and resonance order (5, 16, and 25), the resonant asteroid ratio has three major peaks. For $a_0 = 73$ au, the resonant ratio peaks ($T_0, N_{res}/N_s$) are (0, 16 per cent), (0.7, 19 per cent), and (1.4, 6 per cent). For $a_0 = 100$ au, the peaks are (−0.1, 30 per cent), (0.7, 34 per cent), and (1.4, 36 per cent) and for $a_0 = 200$ au, they are (−0.3, 20 per cent), (0.5, 28 per cent), and (1.5, 6 per cent). The similar peak locations with respect to the Tisserand invariant indicate the existence of preferred inclinations for resonance capture within the polar corridor.

Another property of the resonant ratio is how it conserves the same profile from epoch 0.1 to 4.5 Gyr for $a_0 = 100$ and 200 au and from 1 to 4.5 Gyr for $a_0 = 73$ au. The profile amplitudes increase with time confirming that resonances increase the stability of trapped collision singularity asteroids. The reason why at $a_0 = 73$ au the resonant ratio profile of epoch 0.1 Gyr is different at later epochs is
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Ratio of the number of injected stable orbits present at 4.5 Gyr, $N_s^i$, to the number of all injected unstable orbits from the simulation's start up to 4.5 Gyr, $N_u^i$, as a function of the initial Tisserand invariant (and $I_\infty$). The solid curves correspond, respectively, to $a_0 = 200$ au (red), 100 au (blue), and 73 au (green). The dashed line is $\log N_s^i/N_u^i = -T_0/0.61$.</figcaption></figure>

that the initial semimajor axis lies inside the 9:34 resonance region [72.94,73.08] au. Dynamical relaxation over ~1 Gyr shifts the main retrograde-inclination peak from $T_0 \sim -0.3$ to 0.

For low inclinations, it is the 1:6 resonance that is more efficient at resonance capture as seen by the sharp increase in the resonant ratio. Of the 1000 asteroids at $a_0 = 100$ au and $T_0 = 2.8$, four asteroids survived at 4.5 Gyr one of which in the 1:6 resonance. The 1:17 and 9:34 resonant ratios essentially vanish beyond $T_0 \sim 2$.

The resonant ratios show that resonance capture is not the main reason for the longevity of retrograde inclination collision singularity asteroids. In the ≥4.5 Gyr $T_0$ range of $a_0 = 73$ au, [-1.1: -0.3], the resonant ratio $N_{\text{res}}/N_s \sim 5$ per cent, whereas for the equivalent range of $a_0 = 100$ au, [-1: -0.4], $N_{\text{res}}/N_s \sim 20$ per cent. These estimates do not explain the longevity of the bulk of orbits with $\tau_d \ge 4.5$ Gyr.

## 3.5 Orbit injection

Paper I showed that orbit injection occurs only if $T > -1$. Injected TNOs from the polar corridor ($-1 \le T \le 2$ or equivalently $45^\circ \le I_\infty \le 110^\circ$) may assume polar and retrograde orbits, whereas those with $T > 2$ (or equivalently $I_\infty \le 45^\circ$) may be injected only onto prograde orbits.

Fig. 9 shows the ratio, $N_s^i/N_u^i$, of the number of injected stable orbits present at epoch 4.5 Gyr to the number of all injected unstable orbits from the simulation's start to 4.5 Gyr for the three semimajor axes $a_0 = 73, 100, \text{ and } 200$ au. The Tisserand invariant range starts from $T_0 = -0.9$, because injection is not possible for $T_0 \le -1$ for stable and unstable orbits alike ($N_s^i = N_u^i = 0$) as $a_{\text{refl}} > a_p$ (see Paper I).

Regardless of the semimajor axis, there are two injection modes depending on the Tisserand invariant. For $T_0 > 0.1$ ($I_\infty < 88^\circ$), most injected asteroids have unstable orbits. For $T_0 \le 0.1$ ($I_\infty \ge 88^\circ$), most injected orbits are stable. This property is independent of the dynamical time being larger than the Solar system's age. For instance, whereas for $a_0 = 200$ au, $\tau_d = 3.67$ Gyr, and $N_s/N_u = 0.8$ at $T_0 = -0.6$, the injection ratio $N_s^i/N_u^i = 8$. Regardless of the semimajor axis, the injection ratio may be approximated by $\log N_s^i/N_u^i = -T_0/0.61$ over a wide range of Tisserand invariants. It should also be noted that the two injection modes are present within the polar corridor and that the unstable injection mode for prograde

motion is unchanged as the polar corridor's $T_0 = 2$ boundary is crossed.

Combined with the >4.5 Gyr dynamical times for negative $T_0$ at certain semimajor axes, the stable injection mode implies the possible existence of long-lived primordial nearly polar TNO reservoirs that supply high-inclination Centaurs to the giant planet's domain. The semimajor extent of this stable region is examined in the next section. We further illustrate the stable injection mode in Fig. 10 with two examples of 4.5 Gyr stable collision singularity TNOs with an initial semimajor axis $a_0 = 100$ au and an initial Tisserand invariant $T_0 = -0.5$ that are transferred onto retrograde Centaur orbits. In the left (right) two panels, the TNO's initial parameters are $e_0 = 0.7112$, $I_0 = 108^\circ$, and $\omega_0 = 154^\circ$ ($e_0 = 0.8313$, $I_0 = 113^\circ$, and $\omega_0 = 92^\circ$). The first TNO (panels a and b) is injected early on during its interaction with the planet in two Centaur phases. The first occurs in the interval [0.14:0.2] Gyr where its orbital parameters $e \sim 0.1-0.2$ and $I = 150^\circ$. The second phase occurs in [0.26:0.39] Gyr with $e \sim 0.05-0.65$ and $I = 140^\circ-178^\circ$. At 4.5 Gyr the TNO's semimajor axis $a = 375$ au, the location of the 1:44 resonance with Neptune that the asteroid entered at epoch 4.12 Gyr. The conservation of the Tisserand invariant ensures that the TNO follows the inclination pathway (equation 1) and its eccentricity is limited by the maximal eccentricity curve (equation 8).

After moving away from Neptune as far as 230 au, the second TNO in panels (c) and (d) is injected onto a Centaur orbit in the interval [3.97:4.11] Gyr where its mean orbital elements are $e \sim 0.4$ and $I \sim 160^\circ$. At 4.5 Gyr the TNO's semimajor axis $a = 54$ au. The TNO's evolution differs significantly from that of the previous one except in one aspect, the conservation of the Tisserand invariant makes them follow the same Tisserand inclination pathway (equation 1) and their eccentricity is always limited by equation (8).

## 4 STABLE COLLISION SINGULARITY REGION

We examine the possible extent of the long-lived region with negative Tisserand invariants in the Sun–Neptune–asteroid three-body problem by following the evolution of collision singularity ensembles within the grid $-1 \le T_0 \le 0$ with a 0.1 step and $40 \le a_0 \le 200$ au with a 10 au step. Each grid point ensemble includes 100 asteroids<fnref n="3" /> generated by the method of Section 2 and integrated forward in time to 5 Gyr because of the computational cost in the grid's stable regions.

The dynamical time portrait as a function of the initial semimajor axis $a_0$ and the Tisserand invariant $T_0$ is shown in Fig. 11. Colour codes indicate the dynamical time $\tau_d$. The simulation time span 5 Gyr is also the maximal value of $\tau_d$. Grid points where $\tau_d = 5$ Gyr may have larger dynamical times. For our purposes of searching for possible 4.5-Gyr-old TNO reservoirs that supply high-inclination Centaurs to the giant planets' domain, the 5 Gyr time span is sufficient.

The stable region ($\tau_d \ge 4.5$ Gyr) can be seen to extend up to $a_0 = 170$ au. Beyond this limit, possible material present in the early Solar system would have been depleted as $\tau_d < 4.5$ Gyr. For each semimajor axis inside the stable region, there are islands of Tisserand invariants where the dynamical time is much larger than the Solar system's age. The region enclosed by $-0.9 \le T_0 \le -0.5$ ($\le I_\infty \le$) and $80 \le a \le 140$ au is a continuous domain where $\tau_d \ge 5$ Gyr. In this domain, the initial semimajor axes are close to the mean motion

<footnote n="3">³Except for $a_0 = 100$ and 200 au whose ensembles include 1000 asteroids.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Examples of collision singularity TNOs in the stable region of $a_0 = 100$ au. The two TNOs have initial parameters $e_0 = 0.7112$, $I_0 = 108^\circ$, $\omega_0 = 154^\circ$ (panels a and b), and $e_0 = 0.8313$, $I_0 = 113^\circ$, $\omega_0 = 92^\circ$ (panels c and d). Time evolution over 4.5 Gyr is shown in the left-hand panels of each asteroid (panels a and c) sampled every 1 Myr. Neptune's position is indicated by the dashed black line in the semimajor axis panels. The TNOs' pathways in (a, I) and (a, e) planes are shown in the right-hand panels (panels b and d). The red and green full circles are the initial and final positions, respectively. Also shown in panels (b) and (d) are the inclination pathway (equation 1), the perihelion and aphelion crossing conditions, the maximal eccentricity (equation 8), the planet's position (vertical dashed line), and the reflection semimajor axis (equation 4) (orange vertical line).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Dynamical time portrait as a function of the initial semimajor axis $a_0$ and Tisserand invariant $T_0$. Colour codes indicate the value of the dynamical time $\tau_d$ in Gyr. The simulation time span is 5 Gyr implying that $\tau_d$ at the green positions can exceed 5 Gyr as shown for instance in Fig. 7.</figcaption></figure>

<caption>**Figure 12.** Ratio of the number of injected stable orbits at 4.5 Gyr, $N_s^1$, to the initial asteroid number as a function of the initial semimajor axis $a_0$ and Tisserand invariant $T_0$ for those locations with dynamical times $\tau_d \ge 4.5$ Gyr and $N_s^1 > 0$.</caption>

resonances 3:13 at 80 au and 1:n with $5 \le n \le 10$ for $a_0 = 90–140$ au. In Section 3.4, it was shown that mean motion resonances do not explain the overall stability of collision singularity asteroids for $-1 \le T_0 \le 0$ as only a fraction of the long-lived asteroids are captured in resonance.

In order to identify the possible primordial reservoirs that inject TNOs onto Centaur orbits, we show in Fig. 12 the ratio of the number of injected stable orbits at 4.5 Gyr, $N_s^1$, to the initial asteroid number, 100 (or 1000), only for locations where $\tau_d \ge 4.5$ Gyr and the injected number of stable orbits is finite ($N_s^1 > 0$).

Two regions are of particular interest. The first is nearest to the planet $40 \le a_0 \le 80$ au (blue and green squares) where injection ratios are largest, 15–30 per cent. The second region is the 5 Gyr continuous stable domain. Here the injection ratio peaks at 10 per cent for $a_0 = 90$ and 100 au and $-0.9 \le T_0 \le -0.5$. The finite injection ratios shown in Fig. 12 indicate that all such locations are prime contenders for hosting long-lived primordial nearly polar TNO reservoirs that supply Centaurs to the Solar system. However, as we explain in the next section, the stable region nearest the planet is likely to lose its stability and large injection ability when all giant planets are added to the Sun–Neptune–TNO three-body problem.

# 5 SUMMARY

In this work, we set out to study the injection dynamics of trans-Neptunian asteroids that cross Neptune’s orbit in the three-body problem. Injection in this framework means the process of lowering the value of the asteroid’s semimajor axis so that it becomes smaller than the planet’s. In order to model this process comprehensively, we examined the dynamics of Neptune’s collision singularity at different trans-Neptunian semimajor axes by using the conservation of the Tisserand invariant for planet-crossing asteroids established in Paper I to generate statistical ensembles at the collision singularity and follow their evolution.

We examined the dynamics of stable and unstable motion at the collision singularity and uncovered a two-mode injection process regardless of the initial semimajor axis and the dynamical time. For polar and retrograde inclinations, stable motion dominates injection, whereas for prograde motion unstable motion dominates. There are regions outside Neptune’s orbit where the median lifetime of retrograde collision singularity asteroids exceeds the Solar system’s age. We also found that the collision singularity’s dynamical time has a unique dependence on the Tisserand invariant regardless of the
semimajor axis. Similarly, resonance capture of collision singularity asteroids is more efficient at preferred inclinations in the polar corridor regardless of semimajor axis.

These results' validity holds in the three-body problem. The Solar system's giant planets will affect the injection process and the collision singularity's stability. This influence can be constructive or destructive. Pluto's Neptune-crossing orbit provides an example of Jupiter's constructive influence on perihelion precession that helps confine its latitudinal motion (Malhotra & Ito 2022). In the case of the polar TNO reservoirs uncovered in this work, Paper II indicates that the presence of all giant planets modifies the dynamics at high inclination both constructively and destructively.

To understand this, we recall that Paper II aimed at searching for past stable orbits of 19 high-inclination Centaurs by integrating the equations of motion back in time to $-4.5$ Gyr with all giant planets and the Galactic tide using $\sim 2 \times 10^7$ clones. The idea behind that search was that since Centaurs were believed to have originated from the early planetesimal disc, their orbits are 4.5 Gyr old. Searching statistically for the region in parameter space such objects came from is possible in principle because of the time reversibility of the gravitational equations of motion. In fact, the conservation of the Tisserand invariant makes the inclination pathway taken by a planet-crossing asteroid unchanged whether its motion is followed forwards or backwards in time on Gyr time-scales as shown later in Paper I. The statistical search for stable orbits did not indicate that high-inclination Centaur orbits returned to the planetesimal disc, nor was it inconclusive such as when past orbits scattered all over parameter space. Instead, high-inclination Centaur clustered mainly in the polar corridor with small inclination dispersions on Neptune-crossing orbits (figs 2 and 4 of Paper II).

One of the search's unexpected results was that the Centaur semimajor axis probability density at $-4.5$ Gyr has a prominent peak between 70 and 100 au (fig. 3 in Paper II) indicating the possible existence of a TNO reservoir where the three-body analysis in this paper predicts. However, despite the simulation's high resolution of $2 \times 10^7$ clones, there was no accumulation in the region nearest to Neptune ($40 \le a \le 70$ au) where the three-body analysis of this paper found the largest injection rates. This indicates that the presence of all giant planets likely destroys the stability of that region.

The large inclination of the possible 4.5-Gyr-old TNO reservoirs poses the problem of their origin as at that early epoch, trans-Neptunian space was devoid of Solar system-born material orbiting the Sun perpendicular to the planetesimal disc. Evidence from planet instability-based disc relaxation simulations over 4.5 Gyr indicate that high-inclination Centaurs cannot be produced in sufficient numbers from the early planetesimal disc, and that other sources are necessary to explain their presence (Kaib et al. 2019; Nesvorný et al. 2019). The natural source is the material captured by the early Solar system from the Sun's birth cluster (Fernández & Brunini 2000; Levison et al. 2010; Brasser et al. 2012; Jílková et al. 2016;

Kaib et al. 2019). This work demonstrates, in the context of the three-body problem, the possible existence of pockets of dynamical stability where some primordial interstellar material is protected in the Solar system and interacts with its planets as polar and retrograde Centaurs. Understanding further collision singularity dynamics with the giant planets and how they influence the stable region, found in this paper using the three-body problem, will help narrow the extent of the stability region and ascertain the population size of the possible Centaur-producing polar TNO reservoirs in the Solar system.

## ACKNOWLEDGEMENTS

I thank the reviewer for their useful comments. The numerical simulations were done at the Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur.

## DATA AVAILABILITY

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

## REFERENCES

Brasser R., Duncan M. J., Levison H. F., Schwamb M. E., Brown M. E., 2012, Icarus, 217, 1

Fernández J. A., Brunini A., 2000, Icarus, 145, 580

Fernández J. A., Gallardo T., Young J. D., 2016, MNRAS, 461, 3075

Fraser W. C., Dones L., Volk K., Womack M., Nesvorný D., 2022, in K. Meech, M. Combi, eds, Comets III. University of Arizona Press

Hands T. O., Dehnen W., Gration A., Stadel J., Moore B., 2019, MNRAS, 419, 1064

Jílková L., Hamers A. S., Hammer M., Portegies Zwart S., 2016, MNRAS, 457, 4218

Kaib N. A., Volk K., 2022, in K. Meech, M. Combi, eds, Comets III. University of Arizona Press

Kaib N. A. et al., 2019, AJ, 158, 43

Kozai Y., 1962, AJ, 67, 591

Levison H. F., Duncan M. J., Brasser R., Kaufmann D. E., 2010, Science, 329, 187

Lidov M. L., 1962, Planet. Space Sci., 9, 719

Malhotra R., Ito T., 2022, Proc. Natl. Acad. Sci., 119, 2118692119

Namouni F., 2022, MNRAS, 510, 276 (Paper I)

Namouni F., Morais M. H. M., 2015, MNRAS, 446, 1998

Namouni F., Morais M. H. M., 2017, MNRAS, 467, 2673

Namouni F., Morais M. H. M., 2018, MNRAS, 477, L117

Namouni F., Morais M. H. M., 2020a, preprint (arXiv:2009.09773)

Namouni F., Morais M. H. M., 2020b, MNRAS, 494, 2191 (Paper II)

Nesvorný D. et al., 2019, AJ, 158, 132

von Zeipel H., 1910, Astron. Nachr., 183, 345

This paper has been typeset from a TeX/Latex file prepared by the author.