© M. Saillenfest et al. 2019

1 Introduction

Beyond Neptune, the orbits of distant small bodies around the barycentre of the solar system are subject to two kinds of perturbations: an internal perturbation from the planets (mainly the giant ones), and an external perturbation from the galactic tides, passing stars, and molecular clouds. Historically, a distinction is made between the trans-Neptunian or Kuiper belt objects (with all their subclasses) and the Oort cloud. This distinction was made because the trans-Neptunian population is indeed observed on orbits lying beyond or close to Neptune, whereas the long-period comets coming from the Oort cloud are only observed when they are injected into the inner solar system, making them observable from Earth. These different classes of objects are thought to have been initially populated through distinct mechanisms (see e.g. the recent review by Morbidelli & Nesvorny 2019). However, there is no dynamical boundary between the trans-Neptunian and the Oort cloud populations, and numerical simulations show a continuous transfer of objects in both directions (Fouchard et al. 2017; Kaib et al. 2019). This means that objects that were initially dominantly perturbed by the planets are driven into a region where the galactic tides dominate, and vice versa.

However, the external perturbations are often neglected in simulations of trans-Neptunian objects, even when they feature very distant orbits (Gallardo et al. 2012; Saillenfest et al. 2017a,b; Batygin & Brown 2016; Becker et al. 2017), whereas the internal perturbations are usually neglected in simulations of the Oort cloud, at least beyond a distance threshold (see e.g. Higuchi et al. 2007; Fouchard et al. 2018). These simplifications are not necessarily wrong, but a clear understanding of where the transition occurs is still missing, as well as the behaviour of small bodies when they cross the limit.

In reality, there necessarily exists an intermediate region where perturbations from the planets and from the galactic tides have the same order of magnitude. This region, which is itself a continuous transition rather than a clear boundary, can be thought of as the dynamical frontier between the trans-Neptunian and the Oort cloud populations. Since both types of perturbations are expected to be small in this region, we call it the “inert Oort cloud” throughout the article. Authors generally consider that nothing has happened in this region since the formation of the solar system, excluding a very unlikely star passage going completely through, or an even more unlikely close encounter with a giant molecular cloud. Strong orbital perturbations could only have occurred there in the very early evolutionary stages of the solar system, when it was still in a dense stellar cluster. For this reason, the inert Oort cloud is sometimes called “fossilised”, or “detached”, meaning that the objects it contains could have been placed very early on their current orbits through the interaction with neighbour solar siblings (see e.g. Brasser et al. 2012; Jílková et al. 2015). For such a frozen configuration to be achieved, the objects within this region should have a perihelion far away from the giant planets, and a semi-major axis small enough for the Galactic tides not to be able to significantly change the perihelion distance over long timescales.

A rough idea of the location of the inert Oort cloud can be obtained from previous works. Gladman et al. (2002) showed that the scattering effect by Neptune is significant over long timescales only for perihelion distances below about 45 astronomical units (au). The precise limit actually increases with the semi-major axis (Gallardo et al. 2012), because energy kicks result in larger variations of semi-major axis if the semi-major axis is large. In fact, some observed objects with perihelion beyond 45 au are known to experience scattering (Bannister et al. 2017). In any case, we look here for a rough limit only. The scattering process mostly affects the semi-major axis of small bodies, which diffuses chaotically, while the perihelion distance does not vary much. Later on, Gomes et al. (2005), Gallardo et al. (2012), and Saillenfest et al. (2016, 2017a), showed that the Lidov–Kozai mechanism raised by the giantplanets inside a mean-motion resonance with Neptune is able to raise the perihelion of small bodies beyond 60 au in a few thousands of million years. Contrary to scattering effects, this mechanism induces a variation of perihelion distance and inclination, while the semi-major axis remains at the resonance location. This mechanism, however, is only efficient for semi-major axes smaller than about 500 au. From these studies, one can deduce that the action of the planets is limited to orbits with perihelion distances smaller than about 80 au, and that for perihelion beyond 45 au, the semi-major axis should be smaller than 500 au for the planets to possibly have a substantial effect through mean-motion resonances. As regards the effects of the galactic tides, Fouchard et al. (2017) showed that an object with perihelion in the Jupiter-Saturn region, that is, below 15 au from the Sun, should have a semi-major axis larger than 1600 au for the tides to be able to raise its perihelion beyond 45 au in less than the age of the solar system. In other words, the tides can move its perihelion out of reach of any significant planetary scattering.

Consequently, the inert Oort cloud can be considered as the region where the semi-major axis is smaller than 1600 au and the perihelion distance is larger than 45 au, but the semi-major axis should be larger than 500 au if the perihelion distance is smaller than 80 au. The resulting zone is schematised in Fig. 1.

At the early stages of the solar system history, most of the small bodies had nearly circular and coplanar orbits that were close to, or even intersecting, the trajectories of the planets (see e.g. Tsiganis et al. 2005). By scattering, their semi-major axes then spanned a large range of values, populating the bottom part of Fig. 1. It is therefore common in simulations to see small bodies wander around the inert Oort cloud, roughly following the straight lines of Fig. 1 (see Dones et al. 2004 or Gomes et al. 2015). However, a few bodies have been discovered within this region: (90377) Sedna (Brown et al. 2004), 2012 VP₁₁₃ (Trujillo & Sheppard 2014), and 2015 TG₃₈₇ (Sheppard et al. 2019). After the discovery of 2012 VP₁₁₃, Trujillo & Sheppard (2014) conjectured the existence of a massive stable population lying in this region. However, one can already notice that these bodies all have very eccentric orbits, implying that either dramatic events occurred during the early evolutionary stages of the solar system (e.g. as a result of a dense stellar environment, see Brasser et al. 2012), or this inert Oort cloud may not be as inert as it appears at first glance. The Planet 9 hypothesis could point in this direction (Trujillo & Sheppard 2014; Batygin & Brown 2016), but strong external perturbations would still be required to emplace the planet itself on its distant orbit. Anyway, it should not distract us from studying the complex interplay between planets and galactic tides in this transitional regime. Even though planets and galactic tides have almost negligible effects on the inert Oort cloud over long timescales, their combined effects could pile up and still induce substantial orbital changes over a 4.5 Gyr evolution.

The aim of the present paper is to characterise and explore the long-term dynamics of the inert Oort cloud, driven by the perturbations from both the galactic tides and the giant planets. We will investigate the dynamical mechanisms at play in this region (resonances, chaos) and draw a quantitative picture of the relevant timescales.

Section 2 is devoted to the dynamical model used and its underlying simplifications. General considerations about the long-term dynamics are exposed in Sect. 3. They are followed in Sect. 4 by a detailed exploration of the trajectories allowed through Poincaré surfaces of section. In Sect. 5, we discuss the implicationsof this mixed-type dynamics for real objects, and we map the inert region in the space of orbital elements. We finally conclude in Sect. 6.

Fig. 1 Naiveview of the inert Oort cloud. It is defined as the region where neither the planets nor the galactic tides have substantial effects on the orbit of small bodies. In this schematic view, the planet scattering makes small bodies move horizontally, whereas the galactic tides and the isolated mean-motion resonances with planets (labelled “planet resonances” on the graph) make them move vertically. The orbital inclination of small bodies is not considered in this picture, though it is known to play a role as well (Saillenfest et al. 2017a).

2 Unified model of planets and galactic tides

We considera small body of negligible mass with respect to the giant planets of the solar system. The Hamiltonian function governing its orbital motion can be decomposed into the Sun-body Keplerian part¹, a perturbation due to the planets, and a perturbation due to the galactic tides: $$\mathcal{H}={\mathcal{H}}_0+{\epsilon }_{\text{P}}{\mathcal{H}}_{\text{P}}+{\epsilon }_{\text{G}}{\mathcal{H}}_{\text{G}}.$$(1)

Expressed using Keplerian elements, the two-body part is $${\mathcal{H}}_0=-\frac{\mu }{2a},$$(2)

where a is the semi-major axis of the small body and μ is the gravitational parameter of the Sun.

We assume that the small body never goes inside the orbits of the planets. The Hamiltonian ${\epsilon }_{\text{P}}{\mathcal{H}}_{\text{P}}$ can therefore be expanded in Legendre polynomials. As explained in the Introduction, mean-motion resonances are inefficient in the inert Oort cloud, such that we are allowed to use the averaged perturbation from the planets (whose orbital periods are much smaller thanthe one of the small body). At this level of approximation, effects coming from the small eccentricities and mutual inclinations of the planets are perfectly negligible for such distant small bodies. Consequently, using circular and coplanar orbits for the planets, we obtain $${\epsilon }_{\text{P}}{\mathcal{H}}_{\text{P}}={\epsilon }_{{\text{P}}_0}{\mathcal{H}}_{{\text{P}}_0}+{\epsilon }_{{\text{P}}_2}{\mathcal{H}}_{{\text{P}}_2}+{\epsilon }_{{\text{P}}_4}{\mathcal{H}}_{{\text{P}}_4}+\mathcal{O}({\epsilon }_{{\text{P}}_6}).$$(3)

These terms correspond to the monopole (index 0), quadrupole (index 2), and hexadecapole (index 4), respectively. Their expressions can be taken from Saillenfest et al. (2016), or Laskar & Boué (2010) in a general context: $$\{{\epsilon }_{{\text{P}}_0}{\mathcal{H}}_{{\text{P}}_0}=-\frac{1}{r}{\sum }_{i=1}^N{\mu }_i,{\epsilon }_{{\text{P}}_2}{\mathcal{H}}_{{\text{P}}_2}=\frac{1}{2}\left(\frac{1}{r^3}-\frac{3}{2}\frac{x^2+y^2}{r^5}\right){\sum }_{i=1}^N{\mu }_i{a}_i^2,{\epsilon }_{{\text{P}}_4}{\mathcal{H}}_{{\text{P}}_4}=-\frac{3}{64}\left(\frac{8}{r^5}-40\frac{x^2+y^2}{r^7}+35\frac{{(x^2+y^2)}^2}{r^9}\right){\sum }_{i=1}^N{\mu }_i{a}_i^4.$$(4)

In these expressions, (x, y, z) are the coordinates of the small body in a reference frame centred on the Sun, where the (x, y) is the orbital plane of the planets, and $r\equiv \sqrt{x^2+y^2+z^2}$. We call this frame the “ecliptic” reference frame. The quantities μ_i and a_i are the gravitational parameter and the semi-major axis of the planet i, for a total of N planets. Because of their small semi-major axes, the planets are supposed to be unaffected by the galactic tides, such that this reference frame is inertial (we consider no precession of the ecliptic pole around the galactic pole).

We now consider the coordinates (X, Y, Z) of the small body in a fixed reference frame centred on the Sun, where the (X, Y) plane is the galactic plane. We call it the “galactic” reference frame. We note (X′, Y′, Z′) the coordinates of the small body in an analogous reference frame, but for which at any time the X′ axis points towards the galactic centre. Because of the motion of the Sun in the Galaxy, the latter reference frame is rotating. At lowest-order of approximation, the Sun describes a circular orbit with constant velocity lying in the galactic plane (e.g. Fouchard 2004). We have in this case the relation $$\left(\begin{array}{c}{X}^{\prime }\\ Y^{\prime }\\ Z^{\prime }\end{array}\right)=\left(\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ -\sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right),$$(5)

where the time derivative of θ is a constant and corresponds to the angular velocity of the galactic centre seen from the Sun. In the following, we write it ν_G. In the quadrupolar approximation, the Hamiltonian function describing the orbital perturbations of the small body from the galactic tides can be written $${\epsilon }_{\text{G}}{\mathcal{H}}_{\text{G}}={\nu }_{\text{G}}{P}_{\theta }+{\mathcal{G}}_1\frac{{X^{\prime }}^2}{2}+{\mathcal{G}}_2\frac{{Y^{\prime }}^2}{2}+{\mathcal{G}}_3\frac{{Z^{\prime }}^2}{2},$$(6)

where ${\mathcal{G}}_1$, ${\mathcal{G}}_2$ and ${\mathcal{G}}_3$ are constants encompassing the shape of the galaxy, its mass density, and the inertial forces due to the rotation of the frame (Fouchard 2004). The momentum P_θ is conjugate to the angle θ; it has been introduced such that the Hamiltonian function is autonomous. Using the usual approximation ${\mathcal{G}}_2=-{\mathcal{G}}_1$, we get $${\epsilon }_{\text{G}}{\mathcal{H}}_{\text{G}}={\nu }_{\text{G}}{P}_{\theta }+{\epsilon }_{{\text{G}}_{\text{V}}}{\mathcal{H}}_{{\text{G}}_{\text{V}}}+{\epsilon }_{{\text{G}}_{\text{R}}}{\mathcal{H}}_{{\text{G}}_{\text{R}}},$$(7)

where $$\{{\epsilon }_{{\text{G}}_{\text{V}}}{\mathcal{H}}_{{\text{G}}_{\text{V}}}={\mathcal{G}}_3\frac{Z^2}{2}{\epsilon }_{{\text{G}}_{\text{R}}}{\mathcal{H}}_{{\text{G}}_{\text{R}}}={\mathcal{G}}_2\left(\frac{Y^2-X^2}{2}\cos (2\theta )-XY\sin (2\theta )\right).$$(8)

The symbols V and R are used here in reference to the vertical and radial components of the galactic tides, respectively.

The perturbations due to the planets and due to the galactic tides being both very small with respect to the Keplerian part, they act on a much longer timescale. Therefore, we use a perturbative approach to order one. The resulting Hamiltonian function is obtained by averaging $\mathcal{H}$ (Eq. (1)) over an orbital period. The momentum conjugate to the mean anomaly of the small body becomes a constant of motion, which implies the conservation of the secular semi-major axis (that we still denote a). Dropping the constant parts, the secular Hamiltonian is $$\overline{\mathcal{H}}={\nu }_{\text{G}}{P}_{\theta }+{\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}+{\epsilon }_{{\text{P}}_4}{\overline{\mathcal{H}}}_{{\text{P}}_4}+{\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}+{\epsilon }_{{\text{G}}_{\text{R}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{R}}}.$$(9)

We will now introduce explicit expressions for the small parameters: $${\epsilon }_{{\text{P}}_2}=\frac{1}{a^3}{\sum }_{i=1}^N{\mu }_i{a}_i^2,{\epsilon }_{{\text{G}}_{\text{V}}}=a^2{\mathcal{G}}_3,{\epsilon }_{{\text{P}}_4}=\frac{9}{16}\frac{1}{a^5}{\sum }_{i=1}^N{\mu }_i{a}_i^4,{\epsilon }_{{\text{G}}_{\text{R}}}=a^2{\mathcal{G}}_2.$$(10)

We note (e, I, ω, Ω) the Keplerian elements of the small body in the ecliptic reference frame, with e its eccentricity, I its inclination, ω its argument of perihelion, and Ω its longitude of ascending node. We will use the subscript G for the same quantities measured in the galactic reference frame (excepting e that does not change). Performing the required averages, the different components of Eq. (9) can be written $$\{\text{}\text{}\begin{array}{l}{\overline{\mathcal{H}}}_{{\text{P}}_2}=\frac{1-3{\cos }^{2}I}{8{(1-e^2)}^{3/2}},\hfill \\ {\overline{\mathcal{H}}}_{{\text{P}}_4}=\frac{1}{64{(1-e^2)}^{7/2}}((2+3e^2)(-3+30{\cos }^{2}I-35{\cos }^{4}I)\hfill \\ +10e^2(1-7{\cos }^{2}I){\sin }^{2}I\cos (2\omega )),\hfill \end{array}$$(11)

and $$\{{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}=\frac{{\sin }^2{I}_{\text{G}}}{4}\left(1+\frac{3}{2}{e}^2-\frac{5}{2}{e}^2\cos (2{\omega }_{\text{G}})\right),{\overline{\mathcal{H}}}_{{\text{G}}_{\text{R}}}=-\frac{1}{4}\left(1+\frac{3}{2}{e}^2\right)\cos (2{\Omega }_{\text{G}}-2\theta ){\sin }^2{I}_{\text{G}}+\frac{5}{4}{e}^2(\sin (2{\omega }_{\text{G}})\sin (2{\Omega }_{\text{G}}-2\theta )\cos I_{\text{G}}-\cos (2{\omega }_{\text{G}})\cos (2{\Omega }_{\text{G}}-2\theta )\frac{1+{\cos }^2{I}_{\text{G}}}{2}).$$(12)

We write ψ the inclination of the ecliptic plane in the galactic reference frame, and α its ascending node. Since we have neglected the precession of the ecliptic pole, the angles ψ and α are constant. The ascending node of the ecliptic can therefore be used as the origin of longitudes in the galactic frame, meaning that α ≡ 0. The corresponding conversion formulas between the two reference frames are given in Appendix A, and the values for the physical constants of the problem are gathered in Table 1. We note that the galactic and ecliptic reference frames used here are a natural choice considering the dynamics under study, but they are not the usual IAU ones (this is only a matter of origin of the longitudes).

3 The galactic Laplace plane

The explicit expressions of the small parameters (Eq. (10)) have been chosen such that the Hamiltonian functions ${\overline{\mathcal{H}}}_{{\text{P}}_2}$, ${\overline{\mathcal{H}}}_{{\text{P}}_4}$, ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$, and ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{R}}}$ have the same order of magnitude for e = 0. The secular semi-major axis rules the relative importance of the different perturbation terms. Figure 2 shows that below a ~ 600 au, the planetary perturbations dominate over the galactic tides by more than a factor 10. The situation is reversed beyond a ~ 1500 au. In between,both kinds of perturbations have the same order of magnitude (${\epsilon }_{{\text{P}}_2}$ and ${\epsilon }_{{\text{G}}_{\text{V}}}$ cross at a ~ 950 au). However, since the eccentricity appears at the denominator in ${\overline{\mathcal{H}}}_{{\text{P}}_2}$ (see Eq. (11)), we expect that the planetary perturbations always have a substantial effect in the high-eccentricity regime.

From Fig. 2, it is also clear that in the weakly perturbed intermediate regime, the planetary perturbations are dominated by the quadrupolar term, whereas the galactic tides mostly consist in their vertical component (the radial component is always smaller by one order of magnitude, see Table 1). In the remaining parts of the article, we will therefore limit the study to the simplified Hamiltonian function $$\mathcal{F}={\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}+{\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}},$$(13)

in order to draw a qualitative picture of the dynamics in the intermediate regime. This Hamiltonian has two degrees of freedom, and we will use the canonical Delaunay elements: $$\{g={\omega }_{\text{G}}h={\Omega }_{\text{G}}\text{conjugate to}\{G=L\sqrt{1-e^2}H=L\sqrt{1-e^2}\cos I_{\text{G}},$$(14)

where $L=\sqrt{\mu a}$ is a constant. We note that g only appears in ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$, whereas h only appears in ${\overline{\mathcal{H}}}_{{\text{P}}_2}$ (through cosI). In this context, the galactic coordinates are therefore the most natural coordinates to use.

First of all, we note that the Hamiltonian $\mathcal{F}$ (see Eq. (13)) is very similar to the Hamiltonian governing the secular orbital motion of a satellite perturbed by the Sun and by the J₂ flattening ofits host planet. In the quadrupolar approximation used here, the two Hamiltonian functions are even identical for small eccentricities, as shown in Appendix B. This means that the concept of “Laplace plane” introduced in the satellite case (see e.g. Tremaine et al. 2009) has its equivalent for distant trans-Neptunian objects in the galactic potential. The Laplace plane is normal to the axis around which the orbital angular momentum precesses. In other words, the orbital inclination measured with respect to the Laplace plane is almost constant, while the corresponding longitude of ascending node circulates. More specifically, a Laplace plane corresponds to a fixed points of the dynamics. The results obtained by Tremaine et al. (2009) in the satellite case remain valid here for circular orbits (e = 0 is a fixed point for the eccentricity). We have the same geometry of phase space, as shown in Fig. 3: we recognise the “circular coplanar” equilibria at Ω_G = 0 and π, among which the stable ones correspond to the classical Laplace plane (located equivalently at Ω = π and 0). We also recognise the “circular orthogonal” equilibrium for I_G = I = 90^o and Ω_G = Ω = ±π∕2. Since the phase space is a sphere, we stress that all trajectories oscillate around one of the stable Laplace equilibria. The stability of the equilibria against eccentricity growth is different from the satellite case, but this hasno consequence here considering the timescales involved (see Appendix B for details).

Low-eccentricity orbits initially lying close to the ecliptic (I_G ≈ ψ and Ω_G ≈ 0) should all precess around the classical Laplace plane. Using the expression of the Hamiltonian $\mathcal{F}$ (see Eq. (13)), the inclination of the classical Laplace plane is a root of a second order polynomial in tan I_G. Figure 4 shows the inclination of the classical Laplace plane according to the value of the secular semi-major axis a. For small values of a, this plane is very close to the ecliptic plane, whereas for large values of a, it is very close to the galactic plane. In between, the orbits of small bodies precess about an intermediary plane. This transition occurs in the region where ${\epsilon }_{{\text{P}}_2}$ and ${\epsilon }_{{\text{G}}_{\text{V}}}$ have the same order of magnitude (compare Figs. 2 and 4).

However, for nearly circular orbits, the oscillations about the classical Laplace plane in the transition regime are extremely slow compared tothe age of the solar system (see Fig. 5), meaning that in practice these orbits hardly change at all and are indeed “inert”. More precisely, the half precession period of the orbit pole exceeds the age of the solar system for semi-major axes between 350 and 13 400 au. Below 350 au, the precession is about the ecliptic pole (I ≈const.), so that a set of orbits initially lying close to the ecliptic plane never departs from it. Beyond 13 400 au, the precession is about the galactic pole (I_G ≈const.), such that orbits explore all the values of I between ψ − I_G and ψ + I_G in less than the age of the solar system. Considering a swarm of particles with initially small inclinations I, this upper limit corresponds to the transition between a disc-like and an isotropic region, even though previous authors based their criterion on a full period for very eccentric orbits (see e.g. Higuchi et al. 2007; Fouchard et al. 2018; Vokrouhlický et al. 2019). For such large values of a, the revolution period of Ω_G tends to the value obtained when neglecting the planets (noted $P_{{\Omega }^{*}}$ by Higuchi et al. 2007, see their Fig. 2).

We note that contrary to regular satellites slowly migrating from their formation region, that are expected to remain very close to their local Laplace plane (see e.g. the discussion of Polycarpe et al. 2018 about Iapetus), distant trans-Neptunian objects can be subject to “fast” changes of orbit: either a diffusion of a by planetary scattering, or an overall randomisation due to passing stars. This last mechanism is thought to have been quite efficient in the inert Oort cloud during the early stages of the solar system. This means that the orthogonal equilibrium could be populated as well, or at least, it could a priori play a role in the dynamics of distant trans-Neptunian objects. In the circular case, the orthogonal orbits are however completely frozen, as shown in Fig. 5.

As can be guessed from the expression of the Hamiltonian function, the situation is different for the eccentricity degree of freedom. Indeed, if the orbit is very eccentric (and this is the case for all known distant trans-Neptunian objects) the planetary part of the Hamiltonian is large, bending the Laplace plane towards the ecliptic (Appendix C), and resulting in much shorter timescales than shown in Fig. 5. Eccentric orbits have a much shorter revolution period of Ω_G also when neglecting the planetary perturbations (function $P_{{\Omega }^{*}}$ of Higuchi et al. 2007). The behaviour of eccentric orbits is the subject of the next section.

Fig. 2 Size of the small parameters listed in Eq. (10) with respect to the secular semi-major axis of the small body. The value of the physical parameters used are given in Table 1.

Fig. 3 Level curves of the Hamiltonian function $\mathcal{F}$ (Eq. (13)) for a circular orbit. The semi-major axis taken as parameter is a = 900 au, and the level curves are shown in black. The two graphs show the same level curves for two sets of variables (in order to avoid being misled by coordinate singularities). The dotted curves represent the inclination ψ of the ecliptic, and the blue spots show the location of the ecliptic plane. The red curve shows an example of nearly circular trajectory precessing around the normal to its local Laplace plane, obtained by numerical integration. The eccentricity slightly varies creating a deviation from the initial level curve. The oscillation period around the fixed centre is about 200 Gyrs, in accordance with Fig 5.

Fig. 4 Inclination of the classical Laplace plane with respect to the planetary plane.

Fig. 5 Period of small oscillations about the two kinds of circular Laplace equilibrium. The period of oscillations around the classical equilibrium exceeds 9 Gyrs for semi-major axes between 350 au and 13 400 au (out of the graph). The period of oscillations around the orthogonal equilibrium is linear with a, and it is everywhere dramatically long with respect to the age of the solar system (apart from very small values of a where the overall model is questionable).

4 Exploration of the dynamics

Since the Hamiltonian $\mathcal{F}$ (Eq. (13)) is composed of two parts that dominate respectively in the low- and high-semi-major axis regimes (seeFig. 2), the first step is to understand the two kinds of dynamics taken separately. Both ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$ and ${\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ are integrable, and their dynamics are well known. In this section, we recall briefly their main aspects and study the interplay between the two kinds of perturbations.

4.1 Planetary regime

If ${\epsilon }_{{\text{G}}_{\text{V}}}\ll {\epsilon }_{{\text{P}}_2}$, the dynamics is largely dominated by the planetary perturbations. Expressed in ecliptic coordinates, the Hamiltonian ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$ taken alone is trivially integrable (see Eq. (11)): the momenta are conserved, and the angles circulate with constant velocities. More specifically, e and I are constant, and the precession velocities are $$\dot{\omega }={\epsilon }_{{\text{P}}_2}\frac{3(5{\cos }^{2}I-1)}{8\sqrt{\mu a}{(1-e^2)}^2};\dot{\Omega }={\epsilon }_{{\text{P}}_2}\frac{-3\cos I}{4\sqrt{\mu a}{(1-e^2)}^2}.$$(15)

The angle ω increases for I < 63^o, decreases for 63^o < I < 117^o, and increases again for I > 117^o. In contrast, Ω decreases for I < 90^o and increases for I > 90^o. Figure 6 shows a map of these precession velocities with respect to the ecliptic inclination, as well as the places where their main integer combinations vanish. These combinations cannot be called “resonances” at this stage, because the two degrees of freedom are strictly decoupled when considering the planetary perturbations alone, but we can expect that they will have a dynamical importance in the perturbed problem (see below). Some of these combinations are mentioned by Gallardo et al. (2012) as affecting observed trans-Neptunian objects. As shown by Saillenfest et al. (2016), the hexadecapole (see Eq. (11)) and successive planetary terms only make small libration islands of ω appear at I≈ 63^o and 117^o. These islands have a maximum width of 16.4 au for the perihelion distance, which, for large semi-major axes, represents a very small variation of eccentricity. When expressed in galactic coordinates, the evolution of I_G, ω_G and Ω_G driven by ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$ are combinations of sinusoids (see Appendix A for the conversion formulas).

4.2 Planetary regime weakly perturbed by galactic tides

When the planetary perturbations dominate over the galactic tides (i.e. for small semi-major axes and/or high eccentricities), the effects of galactic tides can be studied in a perturbative approach: the planetary component ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$ (see Eq. (11)) acts as the integrable dominant part of $\mathcal{F}$, while the galactic component ${\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ (see Eq. (12)) acts as a small perturbation.

Since our dominant part ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$ is already expressed in action-angle coordinates (i.e. it does not depend on the angles), the perturbative approach is straightforward. Expressed in ecliptic coordinates, our perturbing part ${\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ is composed of several terms featuring various combinations of ω and Ω (see the complete list in Appendix A). Therefore, because of the galactic tides, such combinations become genuine resonances, whose characteristics can be obtained analytically. The procedure is detailed in Appendix D.

Figures 7 and 8 show the location and widths of all the strongest resonances (the ones that appear at first order in${\epsilon }_{{\text{G}}_{\text{V}}}\text{}$), obtained analytically. These figures are restricted to small perihelion distances, for the planets to remain by far the dominant term of the dynamics. We focus on prograde orbits, since the resonances for I > 90^o are obtained by replacing cosI by − cosI and Ω by − Ω. As shown in the top panel of Figs. 7 and 8, the libration zone of Ω has by far the largest width in inclination (it actually corresponds to the emergence of the orthogonal Laplace equilibrium, see Sect. 3). We note that the resonances ω − Ω and 2ω − Ω are the first ones to overlap when the galactic tides increase (yellow and grey areas). As shown in the bottom panel of Figs. 7 and 8, the resonances ω + Ω and 2ω + Ω are by far the largest ones in perihelion distance. The other resonances are quite small in comparison, and the libration zone of Ω even has a null width in q. The resonances ω ± 2Ω, visible in Fig. 6, do not even appear in Figs. 7 and 8: this means that they only exist at second order of ${\epsilon }_{{\text{G}}_{\text{V}}}$, and have virtually no effect in the weakly perturbed planetary regime. As explained in Appendix D, in addition to resonances, the Hamiltonian function features a term that is responsible for the emergence of the classic Laplace plane: low-inclination orbits do not precess about the ecliptic pole, as in Sect. 4.1, but about an inclined axis.

As shown by Figs. 7 and 8, when we increase the semi-major axis or the perihelion distance, the resonances become very large and overlap massively. For overly large resonances, the whole dynamical structure outlined in this section is actually destroyed: the galactic tides cannot be treated as a small perturbation anymore.

Fig. 6 Precession velocity of Ω and ω in the planetary regime. The colour represents the velocity scale from negative values in blue to positive values in red, with the same colour scale for Ω and ω. Both $\dot{\Omega }$ and $\dot{\omega }$ attain their maximum absolute value at I = 0o and 180o. The locations where the integer combinations $k\dot{\omega }+j\dot{\Omega }$ vanish (limited to ∥k∥, ∥j∥ < 3) are shown by horizontal lines. The inclination values on the left are obtained from Eq. (15), and the corresponding constant angles are written on the right.

Fig. 7 Location and widths of the strongest resonances in the planetary regime weakly perturbed by the galactic tides. The semi-major axis taken as parameter is a = 500 au. For better visibility, the perihelion distance of the resonance centre is directly used as horizontal axis. Top: location and width in inclination (filled areas). Bottom: upper and lower half widths in perihelion distance, for a centre given by the horizontal axis.

4.3 Galactic regime

If ${\epsilon }_{{\text{P}}_2}\ll {\epsilon }_{{\text{G}}_{\text{V}}}$, the dynamics is dominated by the galactic tides. The dynamics driven by ${\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ taken alone has been studied by many authors. The solutions can actually be expressed analytically in terms of elliptic integrals (see Breiter & Ratajczak 2005; Higuchi et al. 2007; Higuchi & Kokubo 2015 and references therein). The quantity $$K=\sqrt{1-e^2}\cos I_{\text{G}}$$(16)

is conserved and can be used as parameter in the Hamiltonian, which, in turns, has only one degree of freedom. Figure 9 shows the level curves of ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ for different values of K. The limit I_G = 0 or 180^o is a stable fixed point whatever the eccentricity; it coincides with the classic Laplace plane in the large-a regime (see Sect. 3), and results in a frozen orbit. Using K as parameter, this is equivalent to e² = 1 − K² (border of the forbidden regions in Fig. 9). The limit e = 0 is a fixed point with circulating Ω_G, but it is unstable for K² < 4∕5, that is, for 27^o < I_G < 153^o. For K² < 4∕5, there are two additional fixed points located at ω_G = π∕2 or 3π∕2 and $$e^2=1-\frac{\sqrt{5K^2}}{2},$$(17)

which is equivalent to the condition $$e^2=1-\frac{5}{4}{\cos }^2{I}_{\text{G}}.$$(18)

These ones are stable, still with circulating Ω_G. Finally, as already noted by Higuchi et al. (2007), the conservation of K implies that the orbit cannot become retrograde if it is prograde, and vice versa. This is the same for the planetary perturbations, even if we go beyond the quadrupolar approximation (Saillenfest et al. 2016), but this time this concerns the galactic inclination I_G, not the ecliptic one I. Moreover, Ω_G is always decreasing if I_G < 90^o and always increasing if I_G > 90^o (the period of its linear part, already mentioned in Sect. 3, is noted $P_{{\Omega }^{*}}$ by Higuchi et al. 2007).

Fig. 8 Same as Fig. 7 for a = 700 au. The hatched regions mean overlap.

Fig. 9 Level curves of the Hamiltonian function ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ (Eq. (12)). The parameters chosen are K2 < 4∕5 (left) and K2 > 4∕5 (right). Top and bottom rows: the same level curves for two sets of variables.

4.4 Intermediate non-integrable regime

In the intermediate regime (say, from a ~ 500 to 2000 au, see Fig. 2), the dynamics features two fully interacting degrees of freedom, represented by the two pairs of conjugate coordinates (g, G) and (h, H). The dynamics is chaotic in general, but can be explored through Poincaré sections, in the spirit of Li et al. (2014) and Saillenfest et al. (2017b). This method allows one to locate the regular trajectories and to determine the size of the chaotic zones. For 10 values of a, and about 15 values of the Hamiltonian spanning the different dynamical regimes for each value of a, we computed simultaneously four Poincaré sections (for increasing and decreasing g and h). We give our conclusions below and show the most representative figures of our sample.

For semi-major axes smaller than 500 au, the dynamics is dominated by the planetary perturbations, meaning that the eccentricity and ecliptic inclination do not vary much, while the angles ω and Ω circulate (see Sect. 4.1). However, the two degrees of freedom now interact, meaning that genuine resonances appear between ω and Ω (see Sect. 4.2). For a as small as 500 au, Fig. 10 shows that such resonances allow quite large variations of the perihelion distance, but at the speed of only a few au per Gyr. When varying the fixed value of the Hamiltonian, we note that the resonances ω + Ω and 2ω + Ω are by far the most prominent ones for prograde orbits (even if this is not immediately obvious with the parameters chosen to draw Fig. 10). The same holds for the resonances ω − Ω and 2ω − Ω for retrograde orbits. As expected, we also observe libration islands of ω at I≈ 63^o and 117^o, and libration islands of Ω at I≈ 90^o. Hence, the 90^o limit is not a barrier anymore for the inclination when considering the perturbed problem. Narrow chaotic regions are present (bottom graph of Fig. 10), as predicted analytically in Sect. 4.2, but this chaos acts on extremely large timescales. We also notice thin resonances that we did not mention in Sect. 4.2: using a perturbative approach, such resonances are only of order ${\epsilon }_{{\text{G}}_{\text{V}}}^2$ or more.

When we increase the semi-major axis, chaos spreads near the separatrices of the main resonances and libration islands of ω and Ω. Figure 11 shows that chaotic flips between prograde and retrograde orbits are possible for a > 600 au, but the very long timescales involved make these flips of little practical interest.

For semi-major axes between about 800 and 1100, the phase space is almost completely filled with chaos (Fig. 12). Stable trajectories only persist for very eccentric orbits, because they are governed almost entirely by the strong planetary perturbations. Moreover, the perihelion distance and the inclination evolve much faster than for a < 800 au, making their variations substantial in a duration comparable to the age of the solar system. This confirms that eccentric orbits evolve much faster than circular orbits studied in Sect. 3.

Finally, Figs. 13 and 14 show that when the semi-major axis exceeds about 1300 au, the galactic structure of the phase space emerges from the chaotic sea and progressively dominate. Resonances and libration islands of ω and Ω disappear, and the situation is now better characterised in galactic coordinates: we retrieve the structure described in Sect. 4.3 and illustrated in Fig. 9 for K² smaller or larger than 4∕5. This means that the ecliptic inclination I oscillate with a very large amplitude while the ecliptic longitude of ascending node Ω stays around 0 or π, as expected from the Laplace plane (Sect. 3). It should be noted, however, that whatever the value of the semi-major axis, there is always a chaotic region for very eccentric orbits, and a stable region for even higher e, because the denominator of ${\overline{\mathcal{H}}}_{{\text{P}}_2}$ diverges and makes the planetary perturbation dominate again (see Eq. (11)).

Fig. 10 Poincaré sections of the dynamics driven by both the planetary and galactic perturbations (Hamiltonian from Eq. (13)). The semi-major axis taken as parameter is a = 500 au. The colour scale shows the maximum orbital change rate between two successive points (see titles), and the grey zones are forbidden. Top: section for $\mathcal{F}=2\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. The islands located at e ≈ 0.8 are due to the resonance 2ω + Ω, and the islands located at e ≈ 0.2 are due to librations of ω while I ≈ 63o (Lidov–Kozai mechanism). Bottom: section for $\mathcal{F}=2\times {10}^{-8}$ au2yr−2 at ωG = 0 decreasing. The chaotic bands are due to the overlap of the resonances ω − Ω and 2ω − Ω (for I < 90o), or ω + Ω and 2ω + Ω (for I > 90o), and the islands located at I ≈ 90o are due to librations of Ω (orthogonal Laplace plane).

Fig. 11 Same as Fig. 10 for a = 600 au. Left: section for $\mathcal{F}=2\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. From top to bottom, the largest islands are due to the resonance 2ω + Ω, the libration of ω at I ≈ 63o, and the resonance 2ω − Ω. Right: section for $\mathcal{F}=5\times {10}^{-9}$ au2yr−2 at ωG = 0 decreasing. The islands located at I ≈ 90o are due to librations of Ω.

Fig. 12 Same as Fig. 10 for a = 800 au (left) and a = 1100 au (right). Left:section for $\mathcal{F}=1\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. Right: section for $\mathcal{F}=8\times {10}^{-9}$ au2yr−2 at ωG = 0 decreasing.

5 Implications for real objects

5.1 Dynamics of known inert-Oort-cloud bodies

The dynamics of the three observed members of the inert Oort cloud (Sedna, 2012 VP₁₁₃, and 2015 TG₃₈₇) are known to be stable, even though Sheppard et al. (2019) mentioned that 2015 TG₃₈₇ is at the limit of destabilisation by galactic tides (one of its numerical clones was ejected from the solar system). As expected from previous results, Poincaré sections computed in the vicinity of the orbit of 2012 VP₁₁₃ (a ≈ 270 au) show that e and I are almost constant while ω and Ω circulate. This is also the case for Sedna (a ≈ 540 au), despite its larger semi-major axis, because its orbital elements, and in particular its low inclination with respect to the ecliptic, make it unable to reach any of the features shown in Fig. 10. The case of 2015 TG₃₈₇ is more interesting (a ≈ 1190 au), because its high semi-major axis corresponds to the region where both the planetary and galactic perturbations have substantial effects. Figure 15 shows that 2015 TG₃₈₇ is not far from a chaotic zone surrounding the ω + Ω resonance. Its trajectory, however, is strictly quasi-periodic in our simplified model. Increasing the value of its semi-major axis in the uncertainty range leads to faster oscillations of q with a larger amplitude (see Fig. 15, right panel), but 2015 TG₃₈₇ is unable to reach the nearest chaotic zone. The unlikely ejection path found by Sheppard et al. (2019) involves the scattering effects of the giant planets, triggered when the perihelion distance reaches its minimum (see Fig. 1).

In addition to its few observed members, the inert Oort cloud could also contain the hypothetical Planet 9 (“P9”) proposed by Batygin & Brown (2016). Fixing its initial conditions to the nominal orbital elements adopted for instance by Fienga et al. (2016), in particular a = 700 au, q = 280 au, and I = 30^o, we obtain the left panel of Fig. 16. The situation is similar to that of 2015 TG₃₈₇, meaning that the trajectory of P9 is regular but close to a chaotic zone emerging from the ω + Ω resonance. Since P9’s orbit is still quite undetermined (if P9 ever exists), we cannot rule out the possibility that it actually lies inside the chaotic zone. We did not investigate in detail the structure of the chaos in the vicinity of P9’s orbit, but we can mention that the chaos reaches P9 if we increase its semi-major axis by only 50 au. For a = 800 au, the chaotic region extents down to q = 45 au (right graph of Fig. 16), and for a = 900 au it extents down to Neptune-crossing orbits. The long timescale involved would prevent P9 to actually encounter Neptune in 4.5 Gyrs, but its perihelion distance could anyway vary quite substantially, possibly modifying over time the characteristics of its shepherding effect on distant trans-Neptunian objects. This remains true for the updated P9 orbit obtained by Batygin et al. (2019), even though it has a somewhat smaller semi-major axis.

Fig. 13 Same as Fig. 10 for a = 1300 au. Left: section for $\mathcal{F}=1\times {10}^{-10}$ au2yr−2 at ΩG = 0 decreasing. The upper islands are due to librations of ωG with e2 ≈ 1 − K2. Right: section for $\mathcal{F}=2\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. The islands are due to librations of ωG and oscillations of e2 around $1-\sqrt{5K^2}/2$.

Fig. 14 Same as Fig. 10 for a = 2000 au. Left: section for $\mathcal{F}=7\times {10}^{-10}$ au2yr−2 at ΩG = 0 decreasing. Right: section for $\mathcal{F}=5\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. The big islands are due to librations of ωG and oscillations of e2 around $1-\sqrt{5K^2}/2$.

5.2 Evolution of a sample of objects over 4.5 Gyrs

The exploration of the inert-Oort-cloud dynamics conducted in Sect. 4 allowed us to characterise the structure of the phase space in this region, including the location of the chaotic regions. This structure should appear as imprints in the orbital distribution of a large sample of small bodies. For instance, the existence of a Laplace plane (Sect. 3) that is distinct from the ecliptic should naturally produce an accumulation of Ω near the ascending node of the galactic plane (here located at Ω = π). Moreover, we found in Sect. 4.4 that the combination ϖ = ω + Ω is among the strongest resonances in the transitional regime, which should preferentially orient ϖ around ± π∕2 minus the galactic node. However, this picture was drawn for an infinite timescale, and numerous trajectories flagged as chaotic in the Poincaré sections actually wander over the chaotic zones in a dramatically long timescale, even though no dynamical barrier prevents these trajectories from freely wandering around. In order to determine which of the dynamical structures could be discernible during a time span restricted to the age of the solar system, we monitor the evolution of a swarm of test particles over 4.5 Gyrs. The simplicity of the system under study (Eq. (13)) allows for the propagation of millions of trajectories in a reasonable amount of computation time.

Since we aim to fully explore the parameter space, rather than modelling a realistic population of trans-Neptunian objects, we do not restrict our sample to a particular distribution. Our setup is organised as follows. At first, we use a uniform distribution of semi-major axis a between 100 and 2000 au. In order to manipulate easily understandable variables, we build our initial conditions in ecliptic coordinates. The angles ω and Ω are set uniformly between 0 and 2π, and we structure our exploration in slices of perihelion distance (e.g. q ∈ [40, 60] au, [60, 80] au, etc.) and slices of ecliptic inclination cosine (e.g. cosI ∈ [0.9, 1], [0.8, 0.9], etc.). Each slice is uniformly populated by a sample of 10⁵ test particles, which are integrated numerically for 4.5 Gyrs using Hamilton’s equations of motion applied to the Hamiltonian function from Eq. (13). These propagations still do not contain the planetary scattering, active for low perihelion distances (see Fig. 1), that would add fuzziness in our distributions.

After 4.5 Gyrs, all our samples feature overdensity regions for the ecliptic angles ω and Ω at large semi-major axes. As illustrated in Fig. 17, the extension and shape of these regions are differentaccording to the sample considered, and in some cases, overdensities are noticeable even below a = 500 au. The planetary perturbations alone cannot modify the angular distributions in our samples, because they induce precession velocities that are independent of the angles (see Eq. (15)). Consequently, Fig. 17 demonstrates that the galactic tides have a noticeable effect in 4.5 Gyrs even for moderate values of the semi-major axis. It happens, however, that these overdensity regions have little to do with the dynamical mechanisms (libration zones, resonances) revealed in Sect. 4.4. In fact, as shown by Fig. 18, most of the particles follow only a small portion of the dynamical cycles involved, due to the large timescales at play. The galactic tides induce a gradient of precession velocities with respect to ω and Ω; therefore, orbits that initially precess faster catch up with orbits that precess slower, before all orbits go away on their respective dynamical paths (which can be totally different). This “phase effect” produces temporary overdensity regions, like the ones shown in Fig. 17. The gradient of precession velocities is different for each of our slices of inclination and perihelion distance, producing differing patterns. As shown in Fig. 18, the sharp patterns disappear as time goes by, replaced by actual dynamical features like resonances and libration zones. In other words, the patterns that we observe after 4.5 Gyrs are a direct relic of our initial distribution of particles. At this point, it would be tempting to conclude that the orbits of distant trans-Neptunian objects still keep a clear memory of their primordial distribution, and that all that is needed to extract the relevant dynamical patterns is to find a sufficient number of them. However, other dynamical mechanisms, like the randomisation by passing stars or the presence of distant unseen perturbers could erase the signature that we are looking for. We also stress that the patterns mentioned above (e.g. the ones appearing in Fig. 17) cannot be linked to the clustering of objects that motivates the Planet 9 hypothesis (Batygin et al. 2019), mostly because they form at too large semi-major axes values. Still, we find it a surprising coincidence that the ecliptic plane, the galactic plane, and the proposed P9 orbital plane intersect almost along the same line².

We now focus on the excursion in perihelion distance q of our samples after 4.5 Gyrs. We recall that the galactic tides produce large cycles of eccentricity and inclination (Sect. 4.3), at a rate that increases with the semi-major axis value. The planets, on the contrary, do not change the eccentricity and ecliptic inclination, but induce a precession of ω and Ω that is faster for smaller semi-major axes and smaller perihelion distances (Sect. 4.1). If this precession is fast with respect to the galactic cycles, it has the consequence of averaging to zero the galactic contribution. Hence, as a rule of thumb, the planetary perturbations block the cycles raised by the galactic tides, with an efficiency that decreases for growing semi-major axis and perihelion distance. This is indeed what we observe in Fig. 19, by comparing the behaviour of the samples with and without the planetary perturbations. This “blocking” effect is very efficient at a = 100 au (and even up to about 1500 au in the top middle panel), but almost null at a = 2000 au. It is also less efficient for high perihelion distances (right column of Fig. 19). The spreading of the distribution is damped most for small initial ecliptic inclinations (top row of Fig. 19), because this corresponds to the maximum of the planetary-induced precession velocities (see Fig. 6). Moreover, we note that two ranges of initial ecliptic inclination are much more prone to orbital variations than other ranges, as exemplified by the middle row of Fig. 19. This is particularly visible in Fig. 20, showing the value of the semi-major axis above which the perihelion distance of small bodies, starting from [40, 60] au, can get beyond 80 au in 4.5 Gyrs. The two favoured inclination ranges are approximatively I ∈ [45^o, 55^o] and [125^o, 135^o], which corresponds to the location of the two pairs of resonances (ω + Ω, 2ω + Ω) and (ω − Ω, 2ω − Ω), and the regions where they overlap (see Fig. 6). These resonances being by far the strongest ones in terms of their widths in perihelion distance (see Sects. 4.2 and 4.4), they favour variations of q.

We finally focus on the excursion in ecliptic inclination I of our samples after 4.5 Gyrs. Figure 21 shows that for small perihelion distances, the spreading of the inclination distribution is strongly damped by the planetary perturbations for initial inclinations near I = 0^o (and 180^o). This is similar to what we observed for the perihelion distance (Fig. 19). However, this time, the smallest value of the semi-major axis at which the spreading is substantial is reached for initial ecliptic inclinations I ~ 90^o (see the middle column of Fig. 21), and we observe no marked enhancement of the inclination excursions for other specific ranges of initial inclination. As before, these results are a direct consequence of the form of the galactic potential (see Sect. 4.2). Actually, the particles naturally spread in inclination as they precess about an inclined axis, corresponding to the tilt of the galactic Laplace plane (Sect. 3). For very eccentric orbits, the classic Laplace plane is severely bent towards the ecliptic (Appendix C), but the orthogonal equilibrium at I ~ 90^o produces large oscillations of the inclination (Sect. 4.2).

Fig. 15 Same as Fig. 10 for the parameters of 2015 TG387. The points of its trajectory that cross the section are represented in green. Left: nominal orbital elements given by Sheppard et al. (2019). Right: semi-major axis increased by 3σ. In both panels, the largest islands are due to the resonances ω + Ω (above) and 2ω + 3Ω (below).

Fig. 16 Same as Fig. 15 for the parameters of the Planet 9 hypothesised by Batygin & Brown (2016). Left: nominal orbital elements (see text). Right: semi-major axis increased by 100 au. In both panels, the largest islands are due to the resonances ω + Ω (above) and 2ω + 3Ω (below).

Fig. 17 Density of particles in the plane (a, ϖ = ω + Ω) after 4.5 Gyrs for two slicesof initial conditions taken as examples. Top: initial values q ∈ [40, 60] au, and cos I ∈ [0.7, 0.8]. Bottom: initial values q ∈ [40, 60] au, and cos I ∈ [−0.1, 0]. See text for the remaining initial orbital elements.

Fig. 18 Temporal evolution of orbits with a = 1800 au. The initial conditions are I = 0o, and ϖ = ω + Ω equally distributed in [0, 2π]; the initial value qi of the perihelion distance is written on the right of each graph. The vertical black line marks the 4.5 Gyrs time, at which our density maps (e.g. Fig. 17) are computed.

Fig. 19 Distribution of a few samples of particles after 4.5 Gyrs in the plane (a, q). As indicated in the titles, the first column is for an evolution with only the galactic tides (Hamiltonian ${\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$, Eq. (12)), and the second and third columns are for an evolution with both planets and galactic tides (Hamiltonian $\mathcal{F}$, Eq. (13)). Six slices of initial conditions are shown, as written in the top and right side of the graphs: the first two columns are for initial perihelion distance qi in [30, 40] au, and the thirdone in [80, 100] au, while the cosine of the initial ecliptic inclination Ii is distributed in [0.9, 1] for the top row, [0.5, 0.6] for the middle row, and [0, 0.1] for the bottom row. On each graph, the horizontal black line shows the semi-major axis of Neptune (~ 30 au) for reference.

Fig. 20 Minimum value of the semi-major axis above which particles initially sampled with q ∈ [40, 60] au spread beyond q = 80 au in 4.5 Gyrs, with respect to their initial ecliptic inclination. The horizontal bars show our twenty 0.1-width slices of cosI, connected in their centre by a full curve. The location of the four major resonances for the weakly perturbed planetary Hamiltonianare indicated by dotted black lines.

Fig. 21 Same as Fig. 19, but showing the distribution of the ecliptic inclination.

5.3 Limits of the inert region

Looking at how samples of particles initially distributed in localised regions of the orbital elements space spread under the secular action of the planets and galactic tides, we are now able to answer one of the main questions: what are the limits of the inert Oort cloud? Where in the (a, q, I) space can small bodies like Sedna and 2012 VP₁₁₃ remain efficiently fossilised since the early stages of the solar system evolution?

For eachvalue of (a, q, I), we now look in the plane (ω, Ω) for the initial condition producing the maximum variation of q or I in 4.5 Gyrs. To this end, we save the extrema q_min and q_max reached by q (resp. I) in the course of each numerical integration, and we apply an optimisation algorithm to maximise their difference Δq = q_max − q_min (resp. ΔI). We opted for the Particle Swarm Optimisation method (Poli et al. 2007) in order to limit the cases of convergence towards local maxima. The initial condition (ω, Ω) that maximises Δq is generally different from the one that maximises ΔI, so two separate optimisation procedures are needed.

Using this method, we obtain the full three dimensional structure in the (a, q, I) space of the largest orbital changes produced within our simplified model. Figures 22 and 23 show representative sections of this space inthe (a, q) and (a, I) directions. As expected, the highest orbital variations are reached for orbits with largest semi-major axes, in the regime where the galactic tides strongly dominate over the planetary perturbations (see Sect. 4). The black curve in Figs. 22 and 23 delimits the “inert” portion of the space, defined arbitrarily as Δq < 10 au or ΔI < 5^o. We see that the naive picture depicted in Introduction from previous works largely overestimates the inert region. Fig. 22 shows that orbits are truly inert only if:

• (a)

  a ≳ 500 au and the orbit is nearly circular.

• (b)

  500 ≲ a ≲ 1500 au and q is close to the planetary region.

• (c)

  a ≲ 500 au.

These three inert regions are labelled on the bottom left panel of Fig. 22. They are dynamically distinct:

• (a)

  For nearly circular orbits, we know from Sect. 3 that very large orbital variations are actually allowed by the dynamics, but that the timescale is dramatically long; this means that such orbits hardly even precess in 4.5 Gyrs (unless the semi-major axis is extremely large, see Fig. 5).

• (b)

  For small perihelion distances, the planetary perturbations produce a fast precession of the orbits, which averages out the galactic contribution; this means that such orbits have frozen q and I even when considering an infinite timescale. However, a very small perihelion distance implies that the planetary scattering, not taken into account here, is triggered (see Introduction). In case of scattering, of course, the orbit cannot be considered inert.

• (c)

  For a ≲ 500 au, the inert regions a and b merge. These orbits precess substantially, similarly to region b, even for large perihelion distances. This is where Sedna (a ≈ 540 au) and 2012 VP₁₁₃ (a ≈ 270 au) are located.

In regions b and c, Fig. 23 confirms that the border of the inert region has a complex structure that is directly linked to the main resonances between ω and Ω (see Sect. 4.2). As discussed in Sect. 5.2, this complex structure produces a very marked differential spreading of small bodiesin the space of orbital element. This structure disappears for large perihelion distances, as all resonances overlap. Sedna and 2012 VP₁₁₃ have large perihelion distances (76 au and 81 au), but not large enough to completely suppress the effect of resonances. However, due to their relatively small ecliptic inclinations (12^o and 24^o), both of them are out of any of themain resonances. This makes them true “inert” objects with precessing orbits. As expected from Sect. 5.1, on the contrary, 2015 TG₃₈₇ is out of the inert region depicted in Figs. 22 and 23 (a ≈ 1190 au, q ≈ 65 au, I ≈ 12^o). It is however close to its border.

For completeness, Appendix E gives sections of the (a, q, I) space in the (q, I) direction: we retrieve the resonance structure described analytically in Sect. 4.2.

6 Summary and conclusions

We studied the long-term orbital dynamics of small bodies in the intermediate regime between the Kuiper belt and the Oort cloud, that is, where the planetary perturbations and the galactic tides have the same order of magnitude. The two kinds of perturbations are weak in this region, and we call it the “inert Oort cloud” in reference to the few observed detached Kuiper belt objects, which have extremely stable orbits.

The problem is formally close to the case of a satellite perturbed by the J₂ flattening of its host planet and the averaged attraction from the star. As such, it possesses a tilted Laplace plane (the “galactic Laplace plane”), with a crossover located at about 1000 au. This meansthat for semi-major axes much smaller than this value (say 500 au), circular orbits precess about the ecliptic pole, whereas for semi-major axes much larger than this value (say 1500 au) they precess about the galactic pole. In between, they precess about an intermediately tilted pole. In this regime, however, the precession period for circular orbits counts in hundreds of Gyrs, meaning that these orbits hardly change at all in practice.

These dramatically long timescales are greatly reduced for eccentric orbits. The dynamics is integrable in the small and large semi-major axis regimes, when one kind of perturbation strongly dominates the other one. Between about 800 and 1100 au, however, the phase space is almost completely filled with chaos, from very eccentric down to nearly circular orbits. The chaotic diffusion timescales are quite large, but they decrease with the semi-major axis value. For semi-major axes as small as 800 au, the joint action of planets and galactic tides can produce a chaotic diffusion of perihelion distance q over tens of astronomical units in a few billion years. Even though frozen orbital regions do exist (this is the case for Sedna and 2012 VP₁₁₃), we conclude that this region is far from being inert, contrary to what one could expect from the weakness of the perturbations.

In 4.5 Gyr, the galactic tides have noticeable effects down to semi-major axes of about 500 au. At 2000 au, the orbital excursions induced can exceed 400 au in perihelion distance and 80^o in inclination. Interestingly, the largest changes of perihelion distance are reached for ecliptic inclinations I in the ranges [45^o, 55^o] and [125^o, 135^o]. These ranges are delimited by the two pairs of strong resonances (ω + Ω, 2ω + Ω) and (2ω − Ω, ω − Ω), which ease perihelion variations. When monitoring swarms of particles over 4.5 Gyrs, we also observe accumulations of orbital angles in localised zones, for semi-major axes larger than about 500 au. Indeed, due to the long timescales at play, particles do not have the time to go away on their respective dynamical paths, but they rather spread in a non uniform manner, creating (temporary) overdensities. Such accumulations are a relic of the initial distribution of small bodies, but they have little observational consequences at this stage, considering the very distant objects involved.

In conclusion, when mapping the truly “inert” region (Δq < 10 au and ΔI < 5^o over 4.5 Gyrs), we find that it is remarkably small. The precise limits of the inert region can be found in Figs. 22 and 23. It is either composed of: (a) nearly circular orbits with a ≳ 500 au, (b) orbits with 500 ≲ a ≲ 1500 au and perihelion close to the planetary region, or (c) orbits with a ≲ 500 au (as long as they are unaffected by mean-motion resonances). Moreover, orbits are truly inert only if theirperihelion distance is high enough to avoid planetary scattering; in region b, this only leaves a thin inert zone. Out of the inert region, the excursions mentioned above in perihelion distance and inclination, as well as the angular accumulations, are direct effects of the galactic tides. They are quite noticeable after 4.5 Gyrs, and can therefore be decisive when classifying observed bodies as “detached” or not, or when monitoring samples of them, as is done for P9 simulations (see e.g. Batygin et al. 2019 and references therein). Hence, we advocate including the galactic tides in numerical simulations of trans-Neptunian object with semi-major axis larger than 500 au.

Fig. 22 Limits of the inert region in the (a, q) plane. Each column corresponds to a different value of the initial ecliptic inclination (see titles). The colour scale represents the maximum possible orbital variations in 4.5 Gyrs: the top row shows the variation of ecliptic inclination, and the bottom row shows the variations of perihelion distance (see labels on the right). The black level corresponds to a variation of 5o in inclination (top row) or 10 au in perihelion distance (bottom row). Below the black level, the region can be considered inert. See text for the white symbols.

Fig. 23 Same as Fig. 22 but in the (a, I) plane. Each column corresponds to a different value of the initial perihelion distance (see titles).

Acknowledgements

We thank Vacheslav Emel’yanenko for inspiring discussions about this manuscript, as well as for his independentverifications of our results. We also thank Auriane Égal for her help about the Particle Swarm Optimisation method. This work was supported by the Programme National de Planétologie (PNP) of CNRS/INSU, co-funded by CNES.

Appendix A Conversion formulas between the ecliptic and the galactic reference frames

The galacticreference frame used here is defined with the third axis perpendicular to the galactic plane and the first axis directed towards the ascending node of the ecliptic³. We write (I_G, ω_G, Ω_G) the Keplerian elements of the small body measured in this frame. The ecliptic reference frame is defined with the third axis perpendicular to the ecliptic and the same first axis as the galactic reference frame. We write (I, ω, Ω) the Keplerian elements of the small body measured in this reference frame. Passing from one frame to another corresponds to a rotation of ± ψ around the first axis, ψ being the inclination of the ecliptic measured in the galactic reference frame. We obtain $$\cos I_{\text{G}}=\cos \psi \cos I-\sin \psi \cos \Omega \sin I,$$(A.1) $$\{\cos {\Omega }_{\text{G}}\sin I_{\text{G}}=\cos \psi \cos \Omega \sin I+\sin \psi \cos I,\sin {\Omega }_{\text{G}}\sin I_{\text{G}}=\sin \Omega \sin I,$$(A.2)

and $$\{\cos {\omega }_{\text{G}}\sin I_{\text{G}}=\cos \psi \cos \omega \sin I-\sin \psi (\sin \omega \sin \Omega -\cos \omega \cos \Omega \cos I),\sin {\omega }_{\text{G}}\sin I_{\text{G}}=\cos \psi \sin \omega \sin I+\sin \psi (\cos \omega \sin \Omega +\sin \omega \cos \Omega \cos I).$$(A.3)

The inverse formulas are obtained by replacing ψ by − ψ. Using these expressions, the Hamiltonian $\mathcal{F}$ in Eq. (13) can be written both in ecliptic or in galactic coordinates. More precisely, we have $${\overline{\mathcal{H}}}_{{\text{P}}_2}=\frac{1-3{\cos }^{2}I}{8{(1-e^2)}^{3/2}}=\frac{-1}{16{(1-e^2)}^{3/2}}[(3C^2-1)(3{\cos }^2{I}_{\text{G}}-1)+12CS\cos I_{\text{G}}\sin I_{\text{G}}\cos ({\Omega }_{\text{G}})+3S^2{\sin }^2{I}_{\text{G}}\cos (2{\Omega }_{\text{G}})],$$(A.4)

and $${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}=\frac{{\sin }^2{I}_{\text{G}}}{4}\left(1+\frac{3}{2}{e}^2-\frac{5}{2}{e}^2\cos (2{\omega }_{\text{G}})\right)=\frac{-1}{32}[2(3e^2+2)(2C^2{\cos }^{2}I+S^2{\sin }^{2}I-2)-8CS(3e^2+2)\cos I\sin I\cos (\Omega )+5S^2{e}^2{(\cos I+1)}^2\cos (2\omega +2\Omega )+20CS{e}^2(\cos I+1)\sin I\cos (2\omega +\Omega )+10(3C^2-1)e^2{\sin }^{2}I\cos (2\omega )+20CS{e}^2(\cos I-1)\sin I\cos (2\omega -\Omega )+5S^2{e}^2{(\cos I-1)}^2\cos (2\omega -2\Omega )+2S^2(3e^2+2){\sin }^{2}I\cos (2\Omega )],$$(A.5)

where C ≡ cos ψ and S ≡ sin ψ.

Appendix B Hamiltonian function for the satellite case

For comparison purpose, we present here the Hamiltonian function describing the secular evolution of a satellite perturbed by the Sun and the oblateness of its host planet in the quadrupolar approximation. Such a Hamiltonian can be written $$\overline{\mathcal{K}}={\epsilon }_{\text{J}}{\overline{\mathcal{K}}}_{\text{J}}+{\epsilon }_{\odot }{\overline{\mathcal{K}}}_{\odot },$$(B.1)

where $${\epsilon }_{\text{J}}=\frac{1}{a^3}(2{\mu }_{\text{P}}{J}_2{R}_{\text{P}}^2),{\epsilon }_{\odot }=a^2\left(\frac{3{\mu }_{\odot }}{2a_{\odot }^3{(1-e_{\odot }^2)}^{3/2}}\right),$$(B.2)

and $$\{{\overline{\mathcal{K}}}_{\text{J}}=\frac{1-3{\cos }^{2}I}{8{(1-e^2)}^{3/2}},{\overline{\mathcal{K}}}_{\odot }=\frac{{\sin }^2{I}_{\text{G}}}{4}\left(1+\frac{3}{2}{e}^2-\frac{5}{2}{e}^2\cos (2{\omega }_{\text{G}})\right)-\frac{1}{4}{e}^2.$$(B.3)

In these expressions, μ_P, J₂, and R_P are the gravitational parameter, the flattening coefficient, and the equatorial radius of the host planet, respectively, whereas μ_⊙, a_⊙, and e_⊙ are the gravitational parameter, the semi-major axis, and the eccentricity of the Sun, respectively. We keep the same notations as in the rest of the article (see e.g. Appendix A) in order to emphasize the similarities with the distant trans-Neptunian case. This time, however, the indexless orbital elements are measured with respect to the equatorial plane of the host planet, and the G index refers to the orbital plane of the Sun.

The overall Hamiltonian function in Eq. (B.1) should be compared with Eq. (13). It is well known that the averaged quadrupolar effect of inner bodies has the same form as a J₂ flattening of the central body (see e.g. Tremaine et al. 2009). As such, ε_J is strictly equivalent to the parameter ${\epsilon }_{{\text{P}}_{\text{4}}}$ used above (see Eq. (10)), and ${\overline{\mathcal{K}}}_{\text{J}}$ is identical to ${\overline{\mathcal{H}}}_{{\text{P}}_{\text{4}}}$ (see Eq. (11)). Furthermore, we see here that ε_⊙ has the same a² multiplier as ${\epsilon }_{{\text{G}}_{\text{V}}}$ (see Eq. (10)), and that ${\overline{\mathcal{K}}}_{\odot }$ has nearly the same form as ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ (see Eq. (12)), apart from the additional term − e²∕4. The two problems are thereforenot strictly equivalent, unless e = 0. Hence, the results obtained by Tremaine et al. (2009) for strictly circular orbits remain valid in the trans-Neptunian case, like the location of the equilibrium points (see Fig. 3), and their stability against inclination variation.

However, the stability of the circular equilibrium points (reusing the nomenclature of Tremaine et al. 2009)against eccentricity growth are different. This can be shown using linearised equations around the e = 0 equilibrium (and a set of variables that are not singular for circular orbits). In the satellite case, the circular coplanar equilibrium is stable for any a as long as the obliquity of the host planet is smaller than 68.875^o (see Tremaine et al. 2009). In the trans-Neptunian case, on the contrary, using the inclination of the ecliptic from Table 1 (which is smaller than 68.875^o), we find that the circular coplanar equilibrium is unstable against eccentricity growth for a between 875 and 1509 au. This instability, however, is absolutely unable to affect real bodies because it acts on a dramatically long timescale (see Fig. B.1). Finally, the stability of the circular orthogonal equilibrium against eccentricity growth is also different from the satellite case. Indeed, in the trans-Neptunian case, we obtain the stability condition $$a^5⩽\frac{3}{4}\frac{{\sum }_{i=1}^N{\mu }_i{a}_i^2}{{\mathcal{G}}_3},$$(B.4)

which is twice the analogous limit obtained in the satellite case. This corresponds to a semi-major axis of about 904 au. As before, though, the timescales involved here have no physical relevance.

Fig. B.1 Period of small oscillations about the two kinds of circular Laplace equilibrium. Top: coplanar equilibrium; bottom: orthogonal equilibrium. In the grey zone, the equilibrium point is unstable against eccentricity growth; accordingly, the oscillation period of eccentricity is replaced by the period T after which the eccentricity is multiplied by exp(2π) ≈ 535 (black curve).

Appendix C Laplace plane for eccentric orbits

Strictly speaking, the Laplace plane is defined for circular orbits (see Sect. 3 and Tremaine et al. 2009). For eccentric orbits, the two degrees of freedom are fully coupled and the orbit does not precess around a fixed pole. However, one can still get an idea of the geometry of the phase space with e > 0 by plotting the level curves of the Hamiltonian in the (I_G, Ω_G) plane for different values of (e, ω_G). As can be guessed from the expression of $\mathcal{F}$ (see Eq. (13)), we obtain strictly the same geometry as for e = 0 (Fig. 3), but where the equilibrium points at Ω_G = 0 and π, denoting the classic Laplace plane, are shifted in inclination. For example, Figs. C.1and C.2 show the inclination of this “instantaneous Laplace plane” with respect to (e, ω_G) for two given values of the semi-major axis. We see that the instantaneous equilibrium plane is close to the classic Laplace plane, except for very eccentric orbits, where it is severely bent towards the ecliptic (that it reaches for e → 1). This property can be seen in Fig. C.3, where we plot the inclination of this plane averaged over ω_G (“eccentric Laplace plane”) as a function of the semi-major axis. As shown in Appendix D, the eccentric Laplace plane is not only informative, but it has a real dynamical meaning in the weakly perturbed planetary regime.

Fig. C.1 Inclination of the classic Laplace plane with respect to the ecliptic in the eccentric case. It is obtained by considering fixed values of (e, ωG). The semi-major axis taken as parameter is a = 1000 au. The two panels show the same level curves for two sets of variables. Except in the e = 0 case, this inclination is only “instantaneous” because e and ωG actually vary. Dark colours are low inclinations, and light colours are high inclinations, as shown by the labelled levels. The thick black level shows the inclination value that is equal to the one obtained in the circular case. The mean inclination, obtained by averaging over ωG, corresponds to the vertical line ωG = π∕4 on the top panel, or the diagonal line on the bottom panel.

Fig. C.2 Same as Fig. C.1 for a = 2000 au.

Fig. C.3 Inclination of the “eccentric Laplace plane” with respect to the ecliptic. See text for a proper definition. The classic Laplace plane for a circular orbit is shown in black, while red curves show their eccentric counterparts, drawn for fixed values of the perihelion distance (see labels).

Since the known distant trans-Neptunian objects are very eccentric, we expect from Fig. C.3 that they precess about an axis that is quite closer to the ecliptic pole than predicted by the circular case shown in Fig. 4.

Appendix D Dynamics in the weakly perturbed planetary regime

We consider the dynamical system with Hamiltonian $\mathcal{F}={\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}+{\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$, where the expressions for each part can be found in Eqs. (A.4) and (A.5). If the planetary perturbations dominate over the galactic tides, ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$ acts as the integrable dominant part, whereas ${\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ acts as a small perturbation.

Luckily, the unperturbed part is already expressed in action-angle coordinates. This means that, neglecting terms in $\mathcal{O}({\epsilon }_{{\text{G}}_{\text{V}}}^2)$, the long-term behaviour of the system is simply given by the average of $\mathcal{F}$ over the non-resonant angles. This allows us to investigate the effects of each term one by one, and to study the structure of the flow in the vicinity of the resonances:

• (i)

  The first term of ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ does not include the angles; its acts therefore only as a small modulation of the precession velocities $\dot{\omega }$ and $\dot{\Omega }$ governed by ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$ (see Eq. (15) and Fig. 6).

• (ii)

  The second term of ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ is factored by cosΩ. Strictly speaking, this term cannot be called “resonant”, because it features no separatrix. It actually corresponds to the emergence of the classic Laplace plane: the orbit does not precess exactly about the ecliptic pole, as it would for ${\epsilon }_{{\text{G}}_{\text{V}}}=0$, but about a tilted pole. By averaging over ω, the momentum G becomes a constant of motion, and we retrieve exactly the “eccentric Laplace plane” introduced in Appendix. C, that rules the dynamics of the variables (I, Ω).

• (iii)

  All the remaining terms of ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ correspond to resonances and libration zones for ω and Ω. Assuming that there is a single resonance, the resonant angle can be taken as a new independent variable by a linear canonical change of coordinate (unimodular matrix). For instance, we consider the resonance ω + Ω. From the ecliptic Delaunay coordinates (P_ω, P_Ω, ω, Ω), the corresponding change of coordinates is $$\{\sigma =\omega +\Omega \gamma =\Omega \{\Sigma =P_{\omega }\Gamma =P_{\Omega }-P_{\omega }.$$(D.1)

After averaging over the circulating angle γ, we end up with one constant of motion Γ, and one degree of freedom (Σ, σ). Table D.1 gives the constants associated to all resonances that appear at first order in ${\epsilon }_{{\text{G}}_{\text{V}}}$ (that is, the ones directly appearing in the expression of ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ in Eq. (A.5)). The resonance centre is mostly governed by the unperturbed Hamiltonian ${\epsilon }_{{\text{P}}_2}{\overline{\mathcal{H}}}_{{\text{P}}_2}$: it is fixed in inclination (see Fig. 6), but goes from e = 0 to e = 1 when varying the value of the constant quantity given in Table D.1.

Since the resonances are thin in the weakly perturbed problem, we can use the pendulum approximation. This amounts to using a Taylor expansion of the Hamiltonian around the resonance centre Σ₀, keeping terms up to degree 2 for the unperturbed part, and upto degree 0 for the perturbation. The resulting Hamiltonian has the form $${\mathcal{F}}_{\text{res}}=\alpha {(\Sigma -{\Sigma }_0)}^2+\beta \cos 2\sigma \text{},$$(D.2)

which is a pendulum of centre Σ₀ and half width $\sqrt{2|\beta /\alpha |}$. Using the constants of motion in Table D.1, the widths can then either be expressed in terms of the eccentricity or in terms of the inclination.

Figures 7 and 8 show the location and widths of all the resonances that appear at first order in ${\epsilon }_{{\text{G}}_{\text{V}}}\text{}$, computed analytically using this perturbative approach. In these figures, instead of using the value of the constant quantities as parameters, we directly use the perihelion distance of the resonance centre (this allows us to draw all resonances on a single graph). In the pendulum approximation, the upper and lower half widths are equal when they are expressed in canonical coordinates (i.e. Σ), but, as shown by Figs. 7 and 8, this is not necessarily the case in I nor in q (and we indeed observe asymmetric resonances in Sect. 4.4).

The hexadecapolar planetary term (see Eq. (11)) can be easily incorporated using this perturbative approach; it only slightly changes the widths of the ω libration island.

Appendix E Maximum orbital variations reachable in 4.5 Gyrs from the secular action of the planets and galactic tides

In Sect. 5.3, we map the maximum possible variation for q and I in 4.5 Gyrs, according to the location in the (a, q, I) space. Figures 22 and 23 show sections of this space in the (a, q) and (a, I) directions. For completeness, Fig. E.1 shows here the maximum orbital variations in the (q, I) direction. Only a restricted portion of the horizontal axis is shown in order to ease the comparison with Figs. 7 and 8. We retrieve the structure of the resonances studied in Sect. 4.2, except that Fig. E.1 also incorporates a timescale information. Indeed, for large initial perihelion distances, the precession due to the planets is so slow that the trajectories do not have enough time in 4.5 Gyrs to explore the full width of the resonance (or the full extent of the chaotic region in case of resonance overlap). In Fig. E.1, this extension of timescale produces a drop of the orbital variations inside the resonances. However, for even higher perihelion distances, this drop is compensated by the overall increase of the galactic potential, which tends to shorten the timescale (see the background colour gradient, for instance along the horizontal line I_i = 15^o), while producing the generalised chaotic region studied in Sect. 4.4.

Fig. E.1 Maximum possible orbital variations produced in 4.5 Gyrs in the (q, I) plane. Each column corresponds to a different value of the semi-major axis (see titles). This figure is to be compared to Figs. 7 and 8, showing the location and widths of the main resonances obtained analytically.

References

All Tables

All Figures

Fig. 1 Naiveview of the inert Oort cloud. It is defined as the region where neither the planets nor the galactic tides have substantial effects on the orbit of small bodies. In this schematic view, the planet scattering makes small bodies move horizontally, whereas the galactic tides and the isolated mean-motion resonances with planets (labelled “planet resonances” on the graph) make them move vertically. The orbital inclination of small bodies is not considered in this picture, though it is known to play a role as well (Saillenfest et al. 2017a).

In the text

	Fig. 2 Size of the small parameters listed in Eq. (10) with respect to the secular semi-major axis of the small body. The value of the physical parameters used are given in Table 1.
In the text

Fig. 3 Level curves of the Hamiltonian function $\mathcal{F}$ (Eq. (13)) for a circular orbit. The semi-major axis taken as parameter is a = 900 au, and the level curves are shown in black. The two graphs show the same level curves for two sets of variables (in order to avoid being misled by coordinate singularities). The dotted curves represent the inclination ψ of the ecliptic, and the blue spots show the location of the ecliptic plane. The red curve shows an example of nearly circular trajectory precessing around the normal to its local Laplace plane, obtained by numerical integration. The eccentricity slightly varies creating a deviation from the initial level curve. The oscillation period around the fixed centre is about 200 Gyrs, in accordance with Fig 5.

In the text

	Fig. 4 Inclination of the classical Laplace plane with respect to the planetary plane.
In the text

Fig. 5 Period of small oscillations about the two kinds of circular Laplace equilibrium. The period of oscillations around the classical equilibrium exceeds 9 Gyrs for semi-major axes between 350 au and 13 400 au (out of the graph). The period of oscillations around the orthogonal equilibrium is linear with a, and it is everywhere dramatically long with respect to the age of the solar system (apart from very small values of a where the overall model is questionable).

In the text

Fig. 6 Precession velocity of Ω and ω in the planetary regime. The colour represents the velocity scale from negative values in blue to positive values in red, with the same colour scale for Ω and ω. Both $\dot{\Omega }$ and $\dot{\omega }$ attain their maximum absolute value at I = 0o and 180o. The locations where the integer combinations $k\dot{\omega }+j\dot{\Omega }$ vanish (limited to ∥k∥, ∥j∥ < 3) are shown by horizontal lines. The inclination values on the left are obtained from Eq. (15), and the corresponding constant angles are written on the right.

In the text

Fig. 7 Location and widths of the strongest resonances in the planetary regime weakly perturbed by the galactic tides. The semi-major axis taken as parameter is a = 500 au. For better visibility, the perihelion distance of the resonance centre is directly used as horizontal axis. Top: location and width in inclination (filled areas). Bottom: upper and lower half widths in perihelion distance, for a centre given by the horizontal axis.

In the text

	Fig. 8 Same as Fig. 7 for a = 700 au. The hatched regions mean overlap.
In the text

	Fig. 9 Level curves of the Hamiltonian function ${\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$ (Eq. (12)). The parameters chosen are K2 < 4∕5 (left) and K2 > 4∕5 (right). Top and bottom rows: the same level curves for two sets of variables.
In the text

Fig. 10 Poincaré sections of the dynamics driven by both the planetary and galactic perturbations (Hamiltonian from Eq. (13)). The semi-major axis taken as parameter is a = 500 au. The colour scale shows the maximum orbital change rate between two successive points (see titles), and the grey zones are forbidden. Top: section for $\mathcal{F}=2\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. The islands located at e ≈ 0.8 are due to the resonance 2ω + Ω, and the islands located at e ≈ 0.2 are due to librations of ω while I ≈ 63o (Lidov–Kozai mechanism). Bottom: section for $\mathcal{F}=2\times {10}^{-8}$ au2yr−2 at ωG = 0 decreasing. The chaotic bands are due to the overlap of the resonances ω − Ω and 2ω − Ω (for I < 90o), or ω + Ω and 2ω + Ω (for I > 90o), and the islands located at I ≈ 90o are due to librations of Ω (orthogonal Laplace plane).

In the text

	Fig. 11 Same as Fig. 10 for a = 600 au. Left: section for $\mathcal{F}=2\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. From top to bottom, the largest islands are due to the resonance 2ω + Ω, the libration of ω at I ≈ 63o, and the resonance 2ω − Ω. Right: section for $\mathcal{F}=5\times {10}^{-9}$ au2yr−2 at ωG = 0 decreasing. The islands located at I ≈ 90o are due to librations of Ω.
In the text

	Fig. 12 Same as Fig. 10 for a = 800 au (left) and a = 1100 au (right). Left:section for $\mathcal{F}=1\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. Right: section for $\mathcal{F}=8\times {10}^{-9}$ au2yr−2 at ωG = 0 decreasing.
In the text

	Fig. 13 Same as Fig. 10 for a = 1300 au. Left: section for $\mathcal{F}=1\times {10}^{-10}$ au2yr−2 at ΩG = 0 decreasing. The upper islands are due to librations of ωG with e2 ≈ 1 − K2. Right: section for $\mathcal{F}=2\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. The islands are due to librations of ωG and oscillations of e2 around $1-\sqrt{5K^2}/2$.
In the text

	Fig. 14 Same as Fig. 10 for a = 2000 au. Left: section for $\mathcal{F}=7\times {10}^{-10}$ au2yr−2 at ΩG = 0 decreasing. Right: section for $\mathcal{F}=5\times {10}^{-9}$ au2yr−2 at ΩG = 0 decreasing. The big islands are due to librations of ωG and oscillations of e2 around $1-\sqrt{5K^2}/2$.
In the text

	Fig. 15 Same as Fig. 10 for the parameters of 2015 TG387. The points of its trajectory that cross the section are represented in green. Left: nominal orbital elements given by Sheppard et al. (2019). Right: semi-major axis increased by 3σ. In both panels, the largest islands are due to the resonances ω + Ω (above) and 2ω + 3Ω (below).
In the text

	Fig. 16 Same as Fig. 15 for the parameters of the Planet 9 hypothesised by Batygin & Brown (2016). Left: nominal orbital elements (see text). Right: semi-major axis increased by 100 au. In both panels, the largest islands are due to the resonances ω + Ω (above) and 2ω + 3Ω (below).
In the text

	Fig. 17 Density of particles in the plane (a, ϖ = ω + Ω) after 4.5 Gyrs for two slicesof initial conditions taken as examples. Top: initial values q ∈ [40, 60] au, and cos I ∈ [0.7, 0.8]. Bottom: initial values q ∈ [40, 60] au, and cos I ∈ [−0.1, 0]. See text for the remaining initial orbital elements.
In the text

	Fig. 18 Temporal evolution of orbits with a = 1800 au. The initial conditions are I = 0o, and ϖ = ω + Ω equally distributed in [0, 2π]; the initial value qi of the perihelion distance is written on the right of each graph. The vertical black line marks the 4.5 Gyrs time, at which our density maps (e.g. Fig. 17) are computed.
In the text

Fig. 19 Distribution of a few samples of particles after 4.5 Gyrs in the plane (a, q). As indicated in the titles, the first column is for an evolution with only the galactic tides (Hamiltonian ${\epsilon }_{{\text{G}}_{\text{V}}}{\overline{\mathcal{H}}}_{{\text{G}}_{\text{V}}}$, Eq. (12)), and the second and third columns are for an evolution with both planets and galactic tides (Hamiltonian $\mathcal{F}$, Eq. (13)). Six slices of initial conditions are shown, as written in the top and right side of the graphs: the first two columns are for initial perihelion distance qi in [30, 40] au, and the thirdone in [80, 100] au, while the cosine of the initial ecliptic inclination Ii is distributed in [0.9, 1] for the top row, [0.5, 0.6] for the middle row, and [0, 0.1] for the bottom row. On each graph, the horizontal black line shows the semi-major axis of Neptune (~ 30 au) for reference.

In the text

Fig. 20 Minimum value of the semi-major axis above which particles initially sampled with q ∈ [40, 60] au spread beyond q = 80 au in 4.5 Gyrs, with respect to their initial ecliptic inclination. The horizontal bars show our twenty 0.1-width slices of cosI, connected in their centre by a full curve. The location of the four major resonances for the weakly perturbed planetary Hamiltonianare indicated by dotted black lines.

In the text

	Fig. 21 Same as Fig. 19, but showing the distribution of the ecliptic inclination.
In the text

Fig. 22 Limits of the inert region in the (a, q) plane. Each column corresponds to a different value of the initial ecliptic inclination (see titles). The colour scale represents the maximum possible orbital variations in 4.5 Gyrs: the top row shows the variation of ecliptic inclination, and the bottom row shows the variations of perihelion distance (see labels on the right). The black level corresponds to a variation of 5o in inclination (top row) or 10 au in perihelion distance (bottom row). Below the black level, the region can be considered inert. See text for the white symbols.

In the text

	Fig. 23 Same as Fig. 22 but in the (a, I) plane. Each column corresponds to a different value of the initial perihelion distance (see titles).
In the text

	Fig. B.1 Period of small oscillations about the two kinds of circular Laplace equilibrium. Top: coplanar equilibrium; bottom: orthogonal equilibrium. In the grey zone, the equilibrium point is unstable against eccentricity growth; accordingly, the oscillation period of eccentricity is replaced by the period T after which the eccentricity is multiplied by exp(2π) ≈ 535 (black curve).
In the text

Fig. C.1 Inclination of the classic Laplace plane with respect to the ecliptic in the eccentric case. It is obtained by considering fixed values of (e, ωG). The semi-major axis taken as parameter is a = 1000 au. The two panels show the same level curves for two sets of variables. Except in the e = 0 case, this inclination is only “instantaneous” because e and ωG actually vary. Dark colours are low inclinations, and light colours are high inclinations, as shown by the labelled levels. The thick black level shows the inclination value that is equal to the one obtained in the circular case. The mean inclination, obtained by averaging over ωG, corresponds to the vertical line ωG = π∕4 on the top panel, or the diagonal line on the bottom panel.

In the text

	Fig. C.2 Same as Fig. C.1 for a = 2000 au.
In the text

	Fig. C.3 Inclination of the “eccentric Laplace plane” with respect to the ecliptic. See text for a proper definition. The classic Laplace plane for a circular orbit is shown in black, while red curves show their eccentric counterparts, drawn for fixed values of the perihelion distance (see labels).
In the text

	Fig. E.1 Maximum possible orbital variations produced in 4.5 Gyrs in the (q, I) plane. Each column corresponds to a different value of the semi-major axis (see titles). This figure is to be compared to Figs. 7 and 8, showing the location and widths of the main resonances obtained analytically.
In the text