© The Authors 2025

1 Introduction

Solar System Centaurs with inclinations greater than 60^° are conventionally known as high-inclination Centaurs. Like low-inclination Centaurs, they originate from planet-crossing orbits beyond Neptune, but unlike low-inclination Centaurs, their numbers cannot be explained by the protoplanetary disk in the early Solar System (Fernández et al. 2016; Fraser et al. 2022; Kaib & Volk 2022). High-resolution simulations of the past orbits of 19 real high-inclination Centaurs over 4.5 Gyr have shown that they tend to be polar with respect to the Solar System’s invariable plane, accumulating in the so-called polar corridor that extends beyond Neptune’s orbit (Namouni & Morais 2018, 2020, hereafter Papers I and II). The polar corridor’s existence is rooted in the conservation of the Tisserand parameter of Neptune-crossing trans-Neptunian objects (TNOs), which ensures that their dynamical evolution follows Tisserand inclination pathways over gigayear timescales regardless of whether time flows forward or backward (Namouni 2022, hereafter Paper III).

The semimajor axis spread of −4.5 Gyr high-inclination Centaur clone orbits in the polar corridor extends from the scattered disk’s inner edge to the Oort cloud, indicating a wide region of origin (Papers I and II). The Oort cloud, in particular, is a known possible source (Brasser et al. 2012b; Ito & Higuchi 2024). The presence of Oort cloud matter at −4.5 Gyr requires it to have been captured by the Solar System from the material in the stellar cluster of its birth (Fernández & Brunini 2000; Levison et al. 2010; Brasser et al. 2012b; Jílková et al. 2016; Hands et al. 2019; Kaib et al. 2019). Studies of the relaxation of the early flat proto-planetary disk through the orbital instability phase of the Solar System’s planets have demonstrated that the Oort cloud component that was formed through the scattering of disk matter later than −4.5 Gyr cannot explain the numbers of high-inclination Centaurs (Kaib et al. 2019; Nesvorný et al. 2019). Orbits of high-inclination Centaurs may be generated by hypothetical planets in the outer Solar System from an extended Kuiper belt (Batygin & Brown 2016; Lykawka & Ito 2023), but such models neglect the stellar environment of the Solar System known to induce instability over gigayear timescales (Kaib & Raymond 2025).

In this paper, we examine a new possible radial origin suggested by the semimajor axis spread of −4.5 Gyr high-inclination Centaur clone orbits in the polar corridor. The −4.5 Gyr semi-major axis probability distribution functions were found to peak outside Neptune’s orbit in the semimajor axis range [70:100] au (Paper II). These peaks potentially indicate a stable region where TNO reservoirs are currently present in the polar corridor and have been supplying high-inclination Centaurs to the giant planets’ domain. As a first step in understanding the dynamical origin of the semimajor axis distribution peaks, we studied the dynamics of Neptune-crossing TNOs in the Sun-Neptune-TNO three-body problem (Namouni 2024, hereafter Paper IV). This first step is important as the dynamics of Neptune-crossing TNOs is determined to a great extent by the conservation of the Tisserand parameter that gives rise to the concept of Tisserand inclination pathways. A Tisserand inclination pathway describes how orbital inclination evolves as Neptune is approached, and how the TNO is possibly injected as a Centaur before it eventually recedes from the planet to the outer Solar System if it does not experience collision or ejection. The Tisserand inclination pathway depends only on the TNO’s Tisserand parameter and the planet’s semimajor axis. In particular, it is independent of time. This makes the Tisserand parameter the perfect tool for probing original inclinations. The study of the TNO injection process in Paper IV, as a function of the Tisserand parameter, revealed the existence of a stable region in the polar corridor with a semi-major axis range of [40:160] au where the dynamical time of planet-crossing TNOs exceeds the Solar System’s age, suggesting the possible presence of long-lived Centaur-supplying TNO reservoirs. The fact that the stable region found in the three-body problem using time-forward integration is consistent with the −4.5 Gyr semimajor axis distribution peaks obtained in Paper II using time-backward integration in the presence of the four giant planets is encouraging. In this study, we numerically confirm the existence of this stable region in a comprehensive time-forward simulation of the Solar System; in doing so, we determined the likely location of the closest reservoirs to Neptune populated by TNOs that were captured in the polar corridor 4.5 Gyr ago and that produce high-inclination Centaurs.

In Sect. 2, we describe the initial conditions used in this work and their relationship to the stable region found in Paper IV. In particular, the parameter space is enlarged beyond the three-body stable region’s Tisserand parameter range and beyond the perihelion range of the past orbits of high-inclination Centaurs studied in Papers I and II. In Sect. 3, we describe the simulation method used to follow the time-forward evolution of TNO reservoirs subject to the perturbations of the four giant planets, the Galactic tide, and passing stars. The reservoirs’ state at the simulation’s end is discussed as a function of the initial semimajor axis and Tisserand parameter in Sect. 4. The reservoirs’ stability is confirmed, as is Neptune’s role in providing Tisserand inclination pathways for the reservoirs’ material to follow. It is also found that Saturn produces secondary structures in the polar corridor as it retains a fraction of TNO perihelia near its orbit. In Sect. 5, Centaur injection is described. In particular, the number of injected Centaurs is found to depend linearly on the Tisserand parameter. In Sect. 6, the features of Centaurs orbits are described. Additional linear relations with respect to the Tisserand parameter are found: for the Centaur minimum semimajor axis, for the Centaur perihelion at minimum semimajor axis, and for the Centaur inclination at minimum semimajor axis. The third explains why high-inclination Centaurs are unlikely to be produced from low-inclination reservoirs such as the early protoplanetary disk. Section 7 contains concluding remarks. In the remainder of this paper, the Tisserand parameter is specified with respect to Neptune unless stated otherwise.

2 Reservoirs’ initial extent

We built the reservoirs’ initial state using three results from the 4.5 Gyr time-backward numerical simulation of 19 real high-inclination Centaurs using ∼20 million clones (Papers I and II). The first result is the Centaurs’ past polar state. Starting with initial inclinations extending over [62^°:173^°], the 4.5 Gyr end state inclinations extended over [60^°:90^°] with a −standard deviation ∼10^° beyond Neptune’s orbit. An analytical analysis of the three-body problem (Paper III) demonstrated that this clustering is the manifestation of the Tisserand parameter’s conservation, which makes Centaurs and TNOs follow the Tisserand inclination pathway:

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

$$\text{ where }T=\frac{a_p}{a_0}+2{\left[\frac{a_0\left(1-e_0^2\right)}{a_p}\right]}^{\frac{1}{2}}\cos I_0\text{. }$$(2)

Here e₀ and I₀ are the asteroid’s initial eccentricity and inclination and a₀ ad a_p are the asteroid’s and the planet’s semimajor axes. The Tisserand inclination pathway is independent of time and perihelion and defines the polar corridor as the extent in aI-space between I_Tiss(a, 2, a_p) and I_Tiss(a, −1, a_p). This corresponds approximately to the inclination range of the −4.5 Gyr orbits of high-inclination Centaurs.

The second result is the extent of the 4.5 Gyr semimajor axis distributions peaks in the region −[70:100] au. This extent is smaller than the three-body stable region (Paper IV), [40:160] au, where TNO dynamics was followed at Neptune’s collision singularity. The discrepancy in the inner boundary has two reasons. First, the giant planets’ perturbations are stronger at the inner edge of the three-body stable region. Secular perturbations from the giant planets, in particular, increase instability as they can raise significantly the eccentricity of planet-crossing and non-crossing TNOs. Second, the initial conditions at Neptune’s singularity include, for a given Tisserand parameter, TNOs with small perihelia that would have crossed the un-modeled orbits of Jupiter, Saturn and Uranus and became unstable. In contrast, Papers I and II indicated that past high-inclination Centaur orbits in the polar corridor did not have arbitrary perihelia but they clustered mainly in the vicinity of Neptune’s orbit. The discrepancy in the outer boundary is unclear. Despite the proximity of the 100 au boundary to the location of the outer 1:6 mean motion resonance with Neptune, it was shown in Paper IV that resonances do not fully explain the stability of high-inclination Neptune-crossing TNOs.

The Tisserand parameter range of the three-body stable region where the dynamical time exceeds 4.5 Gyr is [−1:0] corresponding to an inclination far from Neptune 90^° ≤ I_∞ ≤ 110^° where I_∞ is given as (Paper III)

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

In Papers III and IV, we demonstrate that a TNO can be injected inside Neptune’s orbit only if its Tisserand parameter T ≥ − 1. Once injected, its Centaur inclination becomes retrograde if T ≤ 2 and remains prograde otherwise. The T = −1 limit is specific to the three-body problem and stems from the Tisserand inclination pathway. The pathway guides the TNO toward the planet up to a reflection semimajor axis that depends on T and Neptune’s semimajor axis a_N given as

$$a_{\text{refl}}=a_{\text{N}}{[T-2+2\sqrt{3-T}]}^{-1}.$$(4)

The TNO then reverses its radial motion and moves away from the planet. The limit T = −1 is when reflection occurs at Neptune’s semimajor axis. For T < − 1, the TNO is reflected outside the planet’s orbit and cannot become a Centaur. The addition of Jupiter, Saturn and Neptune to the simulation modifies the three-body injection process and a larger T-range is necessary. The Tisserand parameter of a TNO that can cross a planet’s orbit must satisfy $-\sqrt{8}\le T\le \sqrt{8}~2.82$ (Paper III). To probe the transfer of prograde orbits to high-inclination retrograde orbits, the reservoirs’ Tisserand parameter range is enlarged to [−2:2.5] with a 0.5-step from the three-body problem stability range of [−1:0]. The additional value T = 2.8 (I_∞ ∼ 8^°) was included to provide a means of comparison with low-inclination dynamics.

We examined two orbit types in the semimajor axis range [40:140] au with a 10 au-step: hot and cold TNOs. Hot TNOs are test particles with eccentric orbits and perihelia in the range [32:50] au with a 2 au-step¹. These reservoirs could be part of the material captured by the Solar System in its birth cluster. The perihelion range is based on the third result of the time-backward simulations that showed polar corridor orbits’ perihelia clustered at Neptune’s orbit. Enlarging the perihelion range beyond Neptune’s semimajor axis explores whether the giant planets’ combined perturbations produce high-inclination Centaurs from an initial population of non-planet-crossing TNOs.

Cold TNOs are test particles with circular orbits. These reservoirs have no plausible physical origin. They are used to assess how distant from Neptune can perihelia be before TNOs are no longer injected as Centaurs. Hot TNO ensembles with (a = 200 au and 32 ≤ q(au)≤ 50 and a 2 au q-step) and (a = 500 au, 50 ≤ q(au)≤ 90 and a 10 au q-step) along with cold ensembles at the same semimajor axes were also simulated to characterize the dynamics far from the reservoirs.

For each set (a, T, q), 100 orbits are generated by drawing the orbital angles, mean longitude, longitude of ascending node and argument of perihelion, from a uniform distribution. As the Tisserand inclination pathway is independent of any angles, the 100 orbits evolve on the same pathway in the three-body problem. Furthermore, regardless of the initial a and q, ensembles with the same T value evolve on the same Tisserand inclination pathway.

Figure 1 shows the reservoirs’ initial extent as a function of eccentricity, inclination and semimajor axis for T = −2, −1.5, −1, 0, 0.5, 2, 2.5, and 2.8. Hot and cold TNOs are shown in red and blue, respectively. Larger perihelia imply larger (smaller) inclinations than the Tisserand pathway’s for (retrograde) prograde motion. This deviation from the inclination pathway increases with the initial semimajor axis and has a limit value, I_circ, for circular orbits given as

$$I_{\text{circ }}\left(a,T,a_{\text{N}}\right)=\arccos \left(\left[T-\frac{a_{\text{N}}}{a}\right]{\left[\frac{4a}{a_{\text{N}}}\right]}^{-\frac{1}{2}}\right),$$(5)

and shown in the inclination panels of Fig. 1 as the blue curve along with the Tisserand inclination pathway (2) as the orange curve.

Fig. 1 Initial state of the TNO reservoirs for different Tisserand parameters, T, as indicated in the top-right corner of each panel. Hot TNOs are shown in red and cold TNOs in blue. The solid orange and blue curves in the inclination panels are, respectively, the Tisserand inclination pathway (2) and the initial inclination of circular orbits (5). In the eccentricity panels, the curves indicate perihelia that match the four giant planets’ semimajor axes.

3 Simulation method

Dynamical evolution is simulated 4.5 Gyr forward in time using the IAS15 N-body integrator of the REBOUND package (Everhart 1985; Rein & Spiegel 2015; Rein & Tamayo 2017). The simulations include the four giant planets, the TNO ensembles and passing stars. The terrestrial planets’ mass was added to the Sun’s. The three-dimensional Galactic tide (Heisler & Tremaine 1986) and relative inclination of the ecliptic and Galactic planes were included using the REBOUNDx package (Tamayo et al. 2020). The Oort constants, A = 15.3 km s⁻¹ kpc⁻¹, B = 11.9 km s⁻¹ kpc⁻¹ , and star density in the solar neighbor-hood, ρ₀ = 0.119M_⊙ pc⁻³, were obtained from Gaia’s first data release (Bovy 2017; Widmark & Monari 2019) and were used in Papers I and II.

The list of passing stars was generated using the algorithm of Heisler et al. (1987) and the stellar categories of Rickman et al. (2008) and their corresponding statistics (García-Sánchez et al. 2001). The 66 230 passing stars are introduced in the simulation according to their Solar System’s encounter epoch. Their evolution is integrated within 1 pc from the Solar System. Multiple stars can be present in the simulation at the same time and therefore interact with one another. The star residency time within 1 pc measured over 4.5 Gyr is 20 kyr ± 13 kyr with a maximum value of 415 kyr.

The orbits of 142 500 TNOs were integrated and their minimum semimajor axes monitored throughout the simulation. TNOs were removed from the simulation in the case of collision or ejection.

Fig. 2 TNO reservoirs at 4.5 Gyr for different Tisserand parameters, T. Hot TNOs are shown in red and cold TNOs in blue. Unstable TNOs are shown in black at the last sampling epoch before instability. The first and third columns display the orbital distributions in the (a, I) and (a, e) planes. The second and fourth columns show I, e, and a at minimum semimajor axis, where a vertical black line denotes Neptune’s semimajor axis. The solid orange and blue curves in the inclination panels are, respectively, the Tisserand inclination pathway ITiss(a, T, aN) (2) and the circular orbits’ initial inclination, Icirc(a, T, aN) (5), with respect to Neptune. The solid green curves for retrograde and prograde inclinations respectively are the Tisserand inclinations pathways with respect to Saturn ITiss(a, −2, aS) (135° inclination) and ITiss(a, 2, aS) (45° inclination). In the eccentricity panels, the curves indicate perihelia that match the four giant planets’ semimajor axes. The Tisserand inclination pathway with respect to Neptune is not shown for T = 2.8 because expression (2) is not valid for T > 2.7 (see Paper III for details).

4 Reservoirs’ state at 4.5 Gyr

The end states of the TNO reservoirs in the Ia- and ea-planes are shown in Fig. 2 for the Tisserand parameters, T = −2, −1.5, −1, 0, 0.5, 2, 2.5 and 2.8 (as in Fig. 1). Hot TNOs are shown in red and cold TNOs in blue. Unstable orbits are shown in black at the last sampled epoch before instability. Also shown are the TNOs’ distributions of minimum semimajor axes achieved over 4.5 Gyr and the corresponding eccentricities and inclinations for all orbital types stable and unstable.

The reservoirs’ semimajor axis extent is similar to the initial state’s. TNOs scattered by the combined effects of the giant planets, Galactic tide and passing stars are also present far outside the reservoirs’ domain. Secular perturbations cause a larger dispersion in eccentricity and inclination for hot TNOs. Cold TNOs remain mainly with low-eccentricity orbits as evidenced by their clustering along the inclination curve of circular orbits (5). As the Tisserand parameter increases toward the polar corridor boundary’s T = 2, instability is greater and the reservoirs are depleted. Examination of the minimum semimajor axis distributions shows that starting from T = −2 minimum perihelia are distributed between Saturn and Neptune but the perihelion range narrows around Neptune’s orbit outside the polar corridor as T increases and initial inclinations get smaller.

The Tisserand inclination pathways with respect to Neptune for the bulk of the reservoirs’ material are generally followed for all values of the Tisserand parameter. The pathways are shown in orange in the inclination panels of Fig. 2. A novel process related to Saturn adds secondary inclination pathways that are present for T ≤ 2 and prominent inside the polar corridor. Saturn appears to retain a small fraction of reservoir material with perihelia at Saturn’s orbit whereas their orbits are pushed back outside Neptune’s. This material follows two specific Tisserand inclination pathways with respect to Saturn regardless of the initial Tisserand parameter, T. These are I_Tiss(a, −2, a_S) and I_Tiss(a, 2, a_S) (2) where a_S is Saturn’s semimajor axis. The pathways correspond to inclinations far from Neptune I_∞ = 135^° and I_∞ = 45° (3) and appear as the green curves in the inclination panels of Fig. 2. These secondary populations are prominent for T = 0 and T = 0.5 in Fig. 2 and their density decreases outside the polar corridor as T tends to −2. Investigating Saturn’s secondary inclination pathways is beyond the scope of this paper.

Figure 3 shows the ratio of the number of surviving hot TNOs to their initial number, denoted f_TNO, for different initial semimajor axes as a function of the Tisserand parameter. For each semimajor axis, all the different perihelia are included in the statistic – the perihelion dependence is discussed in the next section. Four semimajor axis groups are visible. First, TNOs with initial a in the range [70:100] au follow a similar distribution to that of the dynamical time in the three-body problem found in Paper IV (Fig. 7 therein). The dynamical time can be considered a proxy for f_TNO. The three-body dynamical time versus Tisserand parameter, has three peaks: a prominent one near T ∼ −0.7, a moderate one near T ∼ 0.8 and a smaller one near T ∼ 2.2. The 6-body plus passing stars configuration of the present framework displaces the three peaks to T ∼ −1.5 with f_TNO ∼ 60% for a = 80 au, to T ∼ 0.5 with f_TNO ∼ 50% for a = 70 au and to T ∼ 2 with f_TNO ∼ 30% for a = 100 au. This group’s range of [70:100] au is reminiscent of the semimajor axis probability distributions peaks discussed in Sects. 1 and 2.

The ensembles with a = 50 au and a = 60 au form the second group whose retrograde peak has been all but erased in favor a prominent peak at T = 0.5 and f_TNO ∼ 60% for a = 60 au. The third group has a semimajor axis range [110:140] au whose f_TNO decrease steadily with increasing T. The fourth group is made up of the ensembles of a = 40 au, a = 200 au and 500 au that are mainly unstable confirming the existence of two stability boundaries, an inner one between 40 au and 50 au and an outer one between 140 au and 200 au.

The values of f_TNO in the range [50:140] au confirm that TNO reservoirs are long-lived in the presence of all giant planets, passing stars and the Galactic tide. In particular, the peaks at I_∞ = 122^° (T = −1.5) and I_∞ = 80^° (T = 0.5) single out inclinations that have the largest current populations. The average of f_TNO over all semimajor axes for hot and cold TNOs, denoted ⟨f_TNO⟩, is shown a function of the Tisserand parameter in Fig. 4. It indicates that cold TNOs on average have a 4.5 Gyr dynamical time as ⟨ f_TNO⟩ > 50% unsurprisingly since they mainly retain low-eccentricity orbits. Hot TNOs fractions decrease with T peaking at T ∼ −1.5 with ⟨ f_TNO⟩ ∼ 40%, T ∼ 0.5 with ⟨f_TNO⟩ ∼ 35% and T ∼ 2 with ⟨ f_TNO⟩ ∼ 22%.

Fig. 3 Fraction of hot TNOs that survived 4.5 Gyr as a function of T (and I∞) for semimajor axes 50 ≤ a(au) ≤ 100 (top panel), 110 ≤ a(au) ≤ 500, and a = 40 au (bottom panel). The number next to each curve is the initial semimajor axis.

Fig. 4 Average fraction of cold TNOs (in blue) and hot TNOs (in red) that survived 4.5 Gyr as a function of T (and I∞).

5 Centaur injection

The ratio of the number of TNOs whose semimajor axes became smaller than Neptune’s over 4.5 Gyr, to the initial TNO number is denoted f_cent(a, T). The initial TNO number of each semi-major axis includes the full perihelion range and the Centaur number includes stable and unstable orbits. The top panel of Fig. 5 shows f_cent(a, T) for hot TNOs and different semimajor axes as a function of the Tisserand parameter. The Centaur fraction increases with the Tisserand parameter at a larger rate as the semimajor axis is closer to Neptune. TNOs with low to moderate inclination are injected in greater numbers than high-inclination TNOs as their corresponding reservoirs are depleted because of instability. The hot TNO ensembles with a = 200 au and 500 au contribute smaller Centaur numbers also because of instability. On the other hand, instability for TNOs with a = 40 au translates into the largest injection rate because of its proximity to Neptune. For Tisserand parameters T < −1, Centaur fraction is finite indicating that the presence of the giant planets and passing stars breaks Neptune’s reflecting ability that was evidenced in the three-body problem with the existence of a_refl (4) (Papers III and IV) – further evidence is given in the next section.

We defined ⟨ f_cent⟩ as the semimajor axis averaged Centaur fraction in the stable region. The TNO ensembles of a = 40 au, 200 au and 500 au are excluded from the hot TNO statistics because they lie outside the stability region. This will be applied in the remainder of the paper when averaging over semimajor axes. The bottom panel of Fig. 5 shows that ⟨f_cent⟩ of hot and cold TNOs are remarkably similar. The hot TNO ⟨f_cent⟩ follows a simple linear expression that depends on the Tisserand parameter as follows ⟨f_cent⟩(%) = 17 + 4.3T.

More information on the injection process is found by examining the dependence of the Centaur number on perihelion for hot TNOs and on semimajor axis for cold TNOs. Figure 6 shows non-averaged f_cent(a, T, q) along with the corresponding dynamical time, τ_d(a, T, q), for the ensembles with an initial a = 80 au. As stated in Sect. 2, each (a, T, q) ensemble has 100 TNOs and f_cent(a, T, q) is the fraction of those among them that became Centaurs. An ensemble’s dynamical time is its median lifetime. It is meant as a stability indicator and not a dynamical time in the classical sense that is associated with a single orbit and derived from its clones’ evolution. As such, the TNOs’ dynamical time, τ_d, shown here is an average of the classical dynamical times of the orbits that make up the (a, T, q) ensemble. We find, unsurprisingly, that Centaur injection tends to increase with smaller perihelia for low to moderate inclinations (T = 2.8). However, the perihelion dependence become shallower for reservoirs inside the polar corridor (−1 ≤ T ≤ 2) suggesting that it is secular perturbations and not the initial perihelion that control Centaur injection – see the panel of T = 1. The dynamical time is an increasing function of perihelion and reaches 4.5 Gyr in the polar corridor and for T > −1. Its dependence on the Tisserand parameter inside the polar corridor is unclear but could be depend on the vicinity of T to the boundaries T = −1 and T = 2. Different initial semimajor axes of hot TNOs in the region [50:140] au have similar f_cent(a, T, q) and dynamical times as a function of perihelion. The median and maximum dynamical times for the different perihelia respectively are 3 Gyr and 4.5 Gyr, with the latter achieved in the polar corridor. The median and maximum dynamical times are 2.4 Gyr and 3.5 Gyr for a = 200 au and 500 au, whereas for a = 40 au, they are 0.48 Gyr and 1.5 Gyr.

The Centaur fraction, f_cent(a, T, q), and dynamical time of cold TNOs, τ_d(a, T, q), are shown in Fig. 7. The number of injected Centaurs decreases exponentially with semimajor axis so that most injected Centaurs come for the smallest semimajor axes. This is visible on the dynamical times as they reach 4.5 Gyr as soon as a ≥ 60 au regardless of the Tisserand parameter. It is interesting that regardless of which semimajor axis and perihelion contributes most to Centaur injection, the average Centaur fractions, ⟨ f_cent⟩, of the hot and cold reservoirs are similar (Fig. 5 bottom panel). Understanding this similarity’s origin may reveal the common dynamical process underlying Centaur injection.

Fig. 5 Fraction of hot TNOs that became Centaurs as a function of T (and I∞). In the top panel, fcent, corresponds to the initial semimajor axis given next to each curve. In the bottom panel, the cold and hot TNO fractions averaged over initial semi-major axes, ⟨fcent⟩, are shown in blue and red, respectively. The black line is ⟨fcent⟩ (%) = 17 + 4.3T.

Fig. 6 Dynamical time (green) and Centaur fraction (purple) for initially hot TNOs with an initial a = 80 au as a function of perihelion for different Tisserand parameters. The dashed blue line is τd = 4.5 Gyr.

6 Injected Centaur orbits

The analysis of orbit injection in the three-body problem indicates that a Centaur orbit acquires its largest inclination at its minimum semimajor axis as it follows its Tisserand inclination pathway inside the planet’s orbit (Paper III). This makes the inclination at minimum semimajor axis a useful diagnostic tool for accessing pre-injection original inclinations far from the planet (3) through the Tisserand parameter.

Figure 8 shows the Centaurs’ inclination at minimum semi-major axis as a function of the Tisserand parameter for hot TNO reservoirs and the different initial semimajor axes. The statistics for an initial semimajor axis include all initial perihelia. The top panel displays the mean inclinations and the bottom panel the standard deviations. The different initial semimajor axes are not indicated because the curves are close.

The mean inclinations remarkably share the same functional dependence on the Tisserand parameter regardless of semima-jor axis. There is a Tisserand value around T ∼ −0.5, for each semimajor axis, below which the mean inclination is constant and beyond which inclination decreases linearly with the Tisserand parameter. The constant inclination value for T < −0.5 decreases form ∼150◦ for 40 au to ∼120◦ for 140 au.

Figure 9 shows the semimajor axis-averaged inclinations at minimum semimajor axis of Centaurs originating in the stable region as a function of the Tisserand parameter for hot and cold TNO reservoirs along with their standard deviations. The initially hot TNO Centaur-inclination is approximated by ⟨I⟩ = 110◦ − 30◦T for T ≥ −0.67 and ⟨I⟩ = 135◦ otherwise. The cold TNO Centaur-inclination has a similar dependence. The inclination standard deviation increases steadily from σ_I ∼ 25◦ at T = 2.8 to σ_I ∼ 45◦ for T = −2.

The Centaur inclination at minimum semimajor axis shows that high-inclination Centaurs may not be produced by reservoirs with arbitrary initial inclinations. For instance, a low-inclination reservoir with T = 2.8 has ⟨I⟩ ∼ 30^°. Retrograde Centaurs with inclinations 100^° and 155^° respectively lie at 3σ_I and 5σ_I from the mean value. This explains the inability of relaxation models of the early planar protoplanetary disk to account for the numbers of high-inclination Centaurs. On the other hand, reservoirs in the polar corridor can produce polar and retrograde inclination Centaurs with statistical significance.

The minimum semimajor axes (mean and standard deviation) of hot TNOs that achieved a Centaur orbit are shown in Fig. 10 as a function of the Tisserand parameter for different initial semimajor axes. The minimum semimajor axis is averaged over all initial perihelia. Centaur injection occurs across all values of the Tisserand parameter in contrast to the dynamics of the three-body problem. The minimum semimajor axes decrease with increasing T in a similar fashion. The large variations of the TNO ensembles with an initial a = 500 au comes for the small statistics of injected Centaurs. The minimum semi-major axes averaged over the different initial semimajor axes of Centaurs originating in the stable region are shown in Fig. 11 for hot and cold TNOs. The two curves are similar and can be approximated by the linear expression ⟨a_min⟩ (au)= 18.5 − T with a mostly constant standard deviation σ_a ∼ 8au. The three-body reflection semimajor axis of Neptune (4) shown in Fig. 11 is of the same order as ⟨a_min⟩ only for T > 0.4 as it does not account for the presence of the giant planets.

The perihelia at minimum semimajor axis (mean and standard deviation) of hot TNOs that achieved a Centaur orbit are shown in Fig. 12 as a function of the Tisserand parameter for different initial semimajor axes. Mean perihelia and their standard deviations increase with increasing Tisserand parameter (decreasing initial inclination). It is also found that perihelia decrease with increasing initial semimajor axis. In particular, hot TNOs with a = 500 au are injected as Centaurs with perihelia around Jupiter’s orbit. The semimajor axis averaged Centaur perihelia originating in the stable region are shown in Fig. 13. The hot TNO’s mean perihelion as a function of the Tisserand parameter follows the linear relation ⟨q⟩(au) = 6.9 + 0.4T. The standard deviation increases from 3 au at T = −2 to 4 au at T = 2.8.

Fig. 7 Dynamical time (green) and Centaur fraction (purple) for initially cold TNOs as a function of the semimajor axis for different Tisserand parameters. The dashed blue line is τd = 4.5 Gyr. The empty semimajor axis ranges indicate the absence of Centaur injection.

Fig. 8 Centaur inclination at minimum semimajor axis of initially hot TNOs as a function of T (and I∞) for semimajor axes 50 ≤ a (au) ≤ 140 (solid curves) inside the stable region, and a = 40, 200, 500 au (dashed curves) outside the stable region. The color codes are those of Figs. 3 and 5. The average inclination, ⟨I⟩, and standard deviation, σI, are shown in the top and bottom panels, respectively.

Fig. 9 Semimajor axis-averaged Centaur inclination at minimum semi-major axis of cold TNOs (blue) and hot TNOs (red) as a function of T (and I∞). Dashed blue and red lines are the corresponding standard deviations. The solid black lines are ⟨I⟩ = 110° − 30°T for T ≥−0.67 and ⟨I⟩ = 135° otherwise.

Fig. 10 Centaur minimum semimajor axis of initially hot TNOs as a function of T (and I∞) for semimajor axes 50 ≤ a (au) ≤ 140 (solid curves) inside the stable region, and a = 40, 200, 500 au (dashed curves) outside the stable region. The color codes are those of Figs. 3 and 5. The average semimajor axis, ⟨amin⟩, and standard deviation, σa, are shown in the top and bottom panels, respectively.

Fig. 11 Average Centaur minimum semimajor axis of cold TNOs (blue) and hot TNOs (red) as a function of T (and I∞). Dashed blue and red curves are the corresponding standard deviations. The solid black line is ⟨amin⟩ (au) = 18.5 − T. The solid orange curve is the three-body reflection semimajor axis of Neptune (4). The planets’ positions are indicated with dashed black lines.

Fig. 12 Centaur perihelion at minimum semimajor axis of initially hot TNOs as a function of T (and I∞) for semimajor axes 50 ≤ a (au) ≤ 140 (solid curves) inside the stable region, and a = 40, 200, 500 au (dashed curves) outside the stable region. The color codes are those of Figs. 3 and 5. The average semimajor axis, ⟨q⟩, and standard deviation, σq, are shown in the top and bottom panels, respectively.

Fig. 13 Average Centaur perihelion at minimum semimajor axis of cold TNOs (blue) and hot TNOs (red) as a function of T (and I∞). Dashed blue and red curves are the corresponding standard deviations. The solid black line is ⟨q⟩ (au) = 6.9 + 0.4T. The planets’ positions are indicated with dashed black lines.

7 Conclusion

We set out to confirm the possible existence of current TNO reservoirs beyond Neptune’s orbit that supply high-inclination Centaurs to the inner Solar System. We did so by simulating the gravitational perturbations of TNO reservoirs by the four giant planets, passing stars, and the Galactic tide. We based the initial extent of the TNO reservoirs on two previous results. The first is the existence of a stable region in the semimajor range axis [40:160] au in the Sun-Neptune-TNO three-body problem (Paper IV). The second is the proximity of past Centaur perihelia to Neptune’s orbit found in the 4.5 Gyr time-backward simulations of 19 real high-inclination Centaurs that involved 20 million clones (Papers I and II). These extents were enlarged to better understand Centaur injection: the perihelion range was enlarged to [32:50] au, and a cold population with circular orbits was added to examine if Centaur injection is strongly dependent on perihelion. The Tisserand parameter range was enlarged to [−2:2.8], or equivalently to TNO inclinations far from Nep-tune in the range [8^°:135^°], to examine Centaur injection beyond the polar corridor’s retrograde inclination boundary in the Sun-Neptune-TNO three-body problem, I_∞ ∼ 110^° (T = −1), and also at low inclinations of relevance to the relaxation of the early protoplanetary disk, I_∞ ∼ 8^° (T = 2.8).

The numerical simulations confirmed the existence of current TNO reservoirs that provide high-inclination Centaurs to the Solar System with ∼40% of the initial hot TNOs in the polar corridor surviving 4.5 Gyr (Fig. 4) in the semimajor axis range [50:140] au. It was also found that two specific Tisserand parameter values, T = 0.5 and T = 1.5 (initial inclinations at infinity of 80^° and 122^°), have the −largest current populations (Fig. 3). A finer analysis unveiled regions in perihelion space where the dynamical time is 4.5 Gyr and, hence, surviving hot TNOs constitute more than half the original reservoir (Fig. 6).

The role of the Tisserand inclination pathways with respect to Neptune in the reservoirs’ evolution is confirmed (Fig. 2). The presence of all the giant planets modifies the polar corridor’s structure as Saturn becomes the main perturber of a fraction of TNO material, adding secondary Tisserand inclination pathways to the polar corridor (Fig. 2). These pathways are defined only by Saturn regardless of the TNOs’ initial Tisserand parameters. Neptune’s reflection ability of incoming TNOs is also modified by the giant planets because Centaur injection occurs outside the polar corridor’s high-inclination boundary of T = −1. Although the presence of more planets causes some significant departures from three-body predictions, the Tisserand inclination pathway concept is robust as TNO dynamics remain based on it.

The existence of linear relations with respect to the Tisserand parameter that describe the number of injected Centaurs (Fig. 5), the Centaur inclination at minimum semimajor axis (Fig. 9), and the Centaur minimum semimajor axis (Fig. 11) affirms the role of the Tisserand parameter as the main control parameter in TNO dynamics. These relations were obtained by averaging the orbital elements over initial perihelia and semimajor axes for hot and cold TNOs. That the two populations follow similar Tisserand parameter-based linear relations indicates that the Centaur injection process is largely independent of the initial perihelia. The dependence of the mean Centaur perihelia on the initial semimajor axis shows that the larger the TNO initial semimajor axis, the smaller the perihelion (Fig. 12). Beyond this characterization of the injected Centaur orbits of high-inclination-TNO reservoirs, the linear relations will help us discriminate between different high-inclination Centaur origins, such as the Oort cloud component formed by the protoplanetary disk and the excitation caused by additional planets beyond Neptune – provided that Centaur injection is analyzed in terms of the Tisserand parameter using orbital elements at minimum semimajor axis.

The most compelling result regarding Centaur injection is perhaps how the Centaur inclination at minimum semimajor axis depends on the Tisserand parameter and how it is independent of the initial semimajor axis (Fig. 8). This explains why high-inclination Centaurs are not produced with statistical significance from a low-inclination reservoir such as the early protoplanetary disk, and why matter in the polar corridor can yield Centaurs with retrograde orbits.

Matter in the stable region at −4.5 Gyr is unlikely to originate in the early planetesimal disk because it was flat, and to explain Neptune’s current location, it must not have extended beyond ∼30 au (Gomes et al. 2004). The natural source of that matter is the stellar environment of the early Solar System. According to simulations of interstellar matter capture in the Sun’s birth cluster, the inner Oort cloud formed with an inner edge near 100 au (Brasser et al. 2012a), which falls in the middle of the stable region found in this work. Simulations of star encounters in the presence of protoplanetary disks show that it is possible to capture planetesimals at high inclinations while preserving a star’s flat planetesimal disk (Hands et al. 2019). After matter exchange in the Solar System’s early environment has settled, the stable region of this work is where the TNO reservoirs can currently be found nearest to the giant planets. Like the Trojans asteroids found at Jupiter’s stable Lagrange points that are relics of early planet formation, the high-inclination reservoirs are relics of the Solar System’s interactions with the stars of its birth cluster. The Legacy Survey of Space and Time (LSST) of the Vera Rubin Observatory will be able to constrain their extent and population size through its coverage of high ecliptic latitudes necessary to detect eccentric polar TNOs that spend more time above the ecliptic plane (Kurlander et al. 2025).

Acknowledgements

The author thanks the reviewer for their comments. The simulations in this work were performed using HPC resources of the National Facility GENCI-CINES (Grants 2024-AD010415952 and 2025-AD010416591) and on the SIGAMM Cluster hosted at the Observatoire de la Côte d’Azur.

References

All Figures

Fig. 1 Initial state of the TNO reservoirs for different Tisserand parameters, T, as indicated in the top-right corner of each panel. Hot TNOs are shown in red and cold TNOs in blue. The solid orange and blue curves in the inclination panels are, respectively, the Tisserand inclination pathway (2) and the initial inclination of circular orbits (5). In the eccentricity panels, the curves indicate perihelia that match the four giant planets’ semimajor axes.

In the text

Fig. 2 TNO reservoirs at 4.5 Gyr for different Tisserand parameters, T. Hot TNOs are shown in red and cold TNOs in blue. Unstable TNOs are shown in black at the last sampling epoch before instability. The first and third columns display the orbital distributions in the (a, I) and (a, e) planes. The second and fourth columns show I, e, and a at minimum semimajor axis, where a vertical black line denotes Neptune’s semimajor axis. The solid orange and blue curves in the inclination panels are, respectively, the Tisserand inclination pathway ITiss(a, T, aN) (2) and the circular orbits’ initial inclination, Icirc(a, T, aN) (5), with respect to Neptune. The solid green curves for retrograde and prograde inclinations respectively are the Tisserand inclinations pathways with respect to Saturn ITiss(a, −2, aS) (135° inclination) and ITiss(a, 2, aS) (45° inclination). In the eccentricity panels, the curves indicate perihelia that match the four giant planets’ semimajor axes. The Tisserand inclination pathway with respect to Neptune is not shown for T = 2.8 because expression (2) is not valid for T > 2.7 (see Paper III for details).

In the text

	Fig. 3 Fraction of hot TNOs that survived 4.5 Gyr as a function of T (and I∞) for semimajor axes 50 ≤ a(au) ≤ 100 (top panel), 110 ≤ a(au) ≤ 500, and a = 40 au (bottom panel). The number next to each curve is the initial semimajor axis.
In the text

	Fig. 4 Average fraction of cold TNOs (in blue) and hot TNOs (in red) that survived 4.5 Gyr as a function of T (and I∞).
In the text

	Fig. 5 Fraction of hot TNOs that became Centaurs as a function of T (and I∞). In the top panel, fcent, corresponds to the initial semimajor axis given next to each curve. In the bottom panel, the cold and hot TNO fractions averaged over initial semi-major axes, ⟨fcent⟩, are shown in blue and red, respectively. The black line is ⟨fcent⟩ (%) = 17 + 4.3T.
In the text

	Fig. 6 Dynamical time (green) and Centaur fraction (purple) for initially hot TNOs with an initial a = 80 au as a function of perihelion for different Tisserand parameters. The dashed blue line is τd = 4.5 Gyr.
In the text

	Fig. 7 Dynamical time (green) and Centaur fraction (purple) for initially cold TNOs as a function of the semimajor axis for different Tisserand parameters. The dashed blue line is τd = 4.5 Gyr. The empty semimajor axis ranges indicate the absence of Centaur injection.
In the text

Fig. 8 Centaur inclination at minimum semimajor axis of initially hot TNOs as a function of T (and I∞) for semimajor axes 50 ≤ a (au) ≤ 140 (solid curves) inside the stable region, and a = 40, 200, 500 au (dashed curves) outside the stable region. The color codes are those of Figs. 3 and 5. The average inclination, ⟨I⟩, and standard deviation, σI, are shown in the top and bottom panels, respectively.

In the text

	Fig. 9 Semimajor axis-averaged Centaur inclination at minimum semi-major axis of cold TNOs (blue) and hot TNOs (red) as a function of T (and I∞). Dashed blue and red lines are the corresponding standard deviations. The solid black lines are ⟨I⟩ = 110° − 30°T for T ≥−0.67 and ⟨I⟩ = 135° otherwise.
In the text

	Fig. 10 Centaur minimum semimajor axis of initially hot TNOs as a function of T (and I∞) for semimajor axes 50 ≤ a (au) ≤ 140 (solid curves) inside the stable region, and a = 40, 200, 500 au (dashed curves) outside the stable region. The color codes are those of Figs. 3 and 5. The average semimajor axis, ⟨amin⟩, and standard deviation, σa, are shown in the top and bottom panels, respectively.
In the text

	Fig. 11 Average Centaur minimum semimajor axis of cold TNOs (blue) and hot TNOs (red) as a function of T (and I∞). Dashed blue and red curves are the corresponding standard deviations. The solid black line is ⟨amin⟩ (au) = 18.5 − T. The solid orange curve is the three-body reflection semimajor axis of Neptune (4). The planets’ positions are indicated with dashed black lines.
In the text

Fig. 12 Centaur perihelion at minimum semimajor axis of initially hot TNOs as a function of T (and I∞) for semimajor axes 50 ≤ a (au) ≤ 140 (solid curves) inside the stable region, and a = 40, 200, 500 au (dashed curves) outside the stable region. The color codes are those of Figs. 3 and 5. The average semimajor axis, ⟨q⟩, and standard deviation, σq, are shown in the top and bottom panels, respectively.

In the text

	Fig. 13 Average Centaur perihelion at minimum semimajor axis of cold TNOs (blue) and hot TNOs (red) as a function of T (and I∞). Dashed blue and red curves are the corresponding standard deviations. The solid black line is ⟨q⟩ (au) = 6.9 + 0.4T. The planets’ positions are indicated with dashed black lines.
In the text