Chaos in the inert Oort cloud

TL;DR

This paper investigates the long-term chaotic dynamics of small bodies in the inert Oort cloud, revealing the influence of galactic tides and resonances on orbital evolution over billions of years.

Contribution

It characterizes the transition between planetary and galactic influences, identifying chaotic zones and the role of resonances in the inert Oort cloud's dynamics.

Findings

Chaotic regions dominate between 800 and 1100 au.

Orbital angles cluster after 4.5 Gyrs due to slow precession.

Galactic tides affect orbits down to about 500 au.

Abstract

Context: Distant trans-Neptunian objects are subject to planetary perturbations and galactic tides. The former decrease with the distance, while the latter increase. In the intermediate regime where they have the same order of magnitude (the 'inert Oort cloud'), both are weak, resulting in very long evolution timescales. To date, three observed objects can be considered to belong to this category. Aims: We aim to provide a clear understanding of where this transition occurs, and to characterise the long-term dynamics of small bodies in the intermediate regime: relevant resonances, chaotic zones (if any), and timescales at play. Results: There exists a tilted equilibrium plane (Laplace plane) about which orbits precess. The dynamics is integrable in the low and high semi-major axis regimes, but mostly chaotic in between. From 800 to 1100 au, the chaos covers almost all the eccentricity range. The diffusion timescales are large, but not to the point of being indiscernible in a 4.5 Gyrs duration: the perihelion distance can actually vary from tens to hundreds of au. Orbital variations are favoured in specific ranges of inclination corresponding to well-defined resonances. Starting from uniform distributions, the orbital angles cluster after 4.5 Gyrs for semi-major axes larger than 500 au, because of a very slow differential precession. Conclusions: Even if it is characterised by very long timescales, the inert Oort cloud is much less inert than it appears. Orbits can be considered inert over 4.5 Gyrs only in small portions of the space of orbital elements, which include (90377) Sedna and 2012VP113. Effects of the galactic tides are discernible down to semi-major axes of about 500 au. We advocate including the galactic tides in simulations of distant trans-Neptunian objects, especially when studying the formation of detached bodies or the clustering of orbital elements.

Figures

Tables

Table 1. Table 1: Values used for the physical constants of the problem.

Table 2. Table 2: Resonance centres and constants of motion arising from the dynamics in the vicinity of each resonance appearing at first order in  ε G V subscript 𝜀 subscript G V \varepsilon_{\mathrm{G}_{\mathrm{V}}} .

Equations

$$\mathcal{H}={\mathcal{H}}_0+{\epsilon }_{\mathrm{P}}{\mathcal{H}}_{\mathrm{P}}+{\epsilon }_{\mathrm{G}}{\mathcal{H}}_{\mathrm{G}}.$$

$${\mathcal{H}}_0=-\frac{\mu }{2a},$$

$${\epsilon }_{\mathrm{P}}{\mathcal{H}}_{\mathrm{P}}={\epsilon }_{{\mathrm{P}}_0}{\mathcal{H}}_{{\mathrm{P}}_0}+{\epsilon }_{{\mathrm{P}}_2}{\mathcal{H}}_{{\mathrm{P}}_2}+{\epsilon }_{{\mathrm{P}}_4}{\mathcal{H}}_{{\mathrm{P}}_4}+\mathcal{O}({\epsilon }_{{\mathrm{P}}_6}).$$

$$\{\begin{array}{rl}{\epsilon }_{{\mathrm{P}}_0}{\mathcal{H}}_{{\mathrm{P}}_0}& =-\frac{1}{r}\sum _{i=1}^N{\mu }_i,\\ {\epsilon }_{{\mathrm{P}}_2}{\mathcal{H}}_{{\mathrm{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 }_{{\mathrm{P}}_4}{\mathcal{H}}_{{\mathrm{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.\end{array}$$

$$\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),$$

$${\epsilon }_{\mathrm{G}}{\mathcal{H}}_{\mathrm{G}}={\nu }_{\mathrm{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},$$

$${\epsilon }_{\mathrm{G}}{\mathcal{H}}_{\mathrm{G}}={\nu }_{\mathrm{G}}{P}_{\theta }+{\epsilon }_{{\mathrm{G}}_{\mathrm{V}}}{\mathcal{H}}_{{\mathrm{G}}_{\mathrm{V}}}+{\epsilon }_{{\mathrm{G}}_{\mathrm{R}}}{\mathcal{H}}_{{\mathrm{G}}_{\mathrm{R}}},$$

$$\{\begin{array}{rl}{\epsilon }_{{\mathrm{G}}_{\mathrm{V}}}{\mathcal{H}}_{{\mathrm{G}}_{\mathrm{V}}}& ={\mathcal{G}}_3\frac{Z^2}{2}\\ {\epsilon }_{{\mathrm{G}}_{\mathrm{R}}}{\mathcal{H}}_{{\mathrm{G}}_{\mathrm{R}}}& ={\mathcal{G}}_2\left(\frac{Y^2-X^2}{2}\cos (2\theta )-XY\sin (2\theta )\right).\end{array}$$

$$\bar{\mathcal{H}}={\nu }_{\mathrm{G}}{P}_{\theta }+{\epsilon }_{{\mathrm{P}}_2}{\bar{\mathcal{H}}}_{{\mathrm{P}}_2}+{\epsilon }_{{\mathrm{P}}_4}{\bar{\mathcal{H}}}_{{\mathrm{P}}_4}+{\epsilon }_{{\mathrm{G}}_{\mathrm{V}}}{\bar{\mathcal{H}}}_{{\mathrm{G}}_{\mathrm{V}}}+{\epsilon }_{{\mathrm{G}}_{\mathrm{R}}}{\bar{\mathcal{H}}}_{{\mathrm{G}}_{\mathrm{R}}}.$$

$$\begin{array}{rlr}{\epsilon }_{{\mathrm{P}}_2}& =\frac{1}{a^3}\sum _{i=1}^N{\mu }_i{a}_i^2& ,\\ {\epsilon }_{{\mathrm{G}}_{\mathrm{V}}}& =a^2{\mathcal{G}}_3& ,\end{array}\begin{array}{rl}{\epsilon }_{{\mathrm{P}}_4}& =\frac{9}{16}\frac{1}{a^5}\sum _{i=1}^N{\mu }_i{a}_i^4,\\ {\epsilon }_{{\mathrm{G}}_{\mathrm{R}}}& =a^2{\mathcal{G}}_2.\end{array}$$

$$\left\{\begin{aligned} \bar{\mathcal{H}}_{\mathrm{P}_{2}}&=\frac{1-3\cos^{2}I}{8(1-e^{2})^{3/2}}\,,\\ \bar{\mathcal{H}}_{\mathrm{P}_{4}}&=\frac{1}{64(1-e^{2})^{7/2}}\Bigg{(}(2+3e^{2})(-3+30\cos^{2}I-35\cos^{4}I)\\ &\hskip 65.44142pt+10e^{2}(1-7\cos^{2}I)\sin^{2}I\cos(2\omega)\Bigg{)}\,,\end{aligned}\right.$$

$$\left\{\begin{aligned} \bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}&=\frac{\sin^{2}I_{\mathrm{G}}}{4}\left(1+\frac{3}{2}e^{2}-\frac{5}{2}e^{2}\cos(2\omega_{\mathrm{G}})\right)\,,\\ \bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{R}}}&=-\frac{1}{4}\left(1+\frac{3}{2}e^{2}\right)\cos(2\Omega_{\mathrm{G}}-2\theta)\sin^{2}I_{\mathrm{G}}\\ &+\frac{5}{4}e^{2}\Bigg{(}\sin(2\omega_{\mathrm{G}})\sin(2\Omega_{\mathrm{G}}-2\theta)\cos I_{\mathrm{G}}\\ &\hskip 28.45274pt-\cos(2\omega_{\mathrm{G}})\cos(2\Omega_{\mathrm{G}}-2\theta)\frac{1+\cos^{2}I_{\mathrm{G}}}{2}\Bigg{)}\,.\end{aligned}\right.$$

$$\mathcal{F}={\epsilon }_{{\mathrm{P}}_2}{\bar{\mathcal{H}}}_{{\mathrm{P}}_2}+{\epsilon }_{{\mathrm{G}}_{\mathrm{V}}}{\bar{\mathcal{H}}}_{{\mathrm{G}}_{\mathrm{V}}},$$

$$\{\begin{array}{rl}g& ={\omega }_{\mathrm{G}}\\ h& ={\Omega }_{\mathrm{G}}\end{array}\text{conjugate to}\{\begin{array}{rl}G& =L\sqrt{1-e^2}\\ H& =L\sqrt{1-e^2}\cos I_{\mathrm{G}}\end{array},$$

$$\dot{\omega }={\epsilon }_{{\mathrm{P}}_2}\frac{3(5{\cos }^{2}I-1)}{8\sqrt{\mu a}(1-e^2{)}^2};\dot{\Omega }={\epsilon }_{{\mathrm{P}}_2}\frac{-3\cos I}{4\sqrt{\mu a}(1-e^2{)}^2}.$$

$$K=\sqrt{1-e^2}\cos I_{\mathrm{G}}$$

$$e^2=1-\frac{\sqrt{5K^2}}{2},$$

$$e^2=1-\frac{5}{4}{\cos }^2{I}_{\mathrm{G}}.$$

$$\cos I_{\mathrm{G}}=\cos \psi \cos I-\sin \psi \cos \Omega \sin I,$$

$$\{\begin{array}{rl}\cos {\Omega }_{\mathrm{G}}\sin I_{\mathrm{G}}& =\cos \psi \cos \Omega \sin I+\sin \psi \cos I,\\ \sin {\Omega }_{\mathrm{G}}\sin I_{\mathrm{G}}& =\sin \Omega \sin I,\end{array}$$

$$\{\begin{array}{rl}\cos {\omega }_{\mathrm{G}}\sin I_{\mathrm{G}}& =\cos \psi \cos \omega \sin I\\ & -\sin \psi (\sin \omega \sin \Omega -\cos \omega \cos \Omega \cos I),\\ \sin {\omega }_{\mathrm{G}}\sin I_{\mathrm{G}}& =\cos \psi \sin \omega \sin I\\ & +\sin \psi (\cos \omega \sin \Omega +\sin \omega \cos \Omega \cos I).\end{array}$$

$${\bar{\mathcal{H}}}_{{\mathrm{P}}_2}$$

$$=\frac{-1}{16(1-e^2{)}^{3/2}}$$

$$+12CS\cos I_{\mathrm{G}}\sin I_{\mathrm{G}}\cos ({\Omega }_{\mathrm{G}})$$

$$\displaystyle+3\,S^{2}\sin^{2}I_{\mathrm{G}}\,\cos(2\Omega_{\mathrm{G}})\Bigg{]}\,,$$

$${\bar{\mathcal{H}}}_{{\mathrm{G}}_{\mathrm{V}}}$$

$$=\frac{-1}{32}$$

$$-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 )$$

$$\displaystyle+2\,S^{2}(3e^{2}+2)\sin^{2}I\,\cos(2\Omega)\Bigg{]}\,,$$

$$\bar{\mathcal{K}}={\epsilon }_{\mathrm{J}}{\bar{\mathcal{K}}}_{\mathrm{J}}+{\epsilon }_{\odot }{\bar{\mathcal{K}}}_{\odot },$$

$$\varepsilon_{\mathrm{J}}=\frac{1}{a^{3}}\big{(}2\mu_{\mathrm{P}}J_{2}R_{\mathrm{P}}^{2}\big{)}\hskip 14.22636pt,\hskip 14.22636pt\varepsilon_{\odot}=a^{2}\left(\frac{3\mu_{\odot}}{2a_{\odot}^{3}(1-e_{\odot}^{2})^{3/2}}\right)\,,$$

Full text

*Aims. *We aim to provide a clear understanding of where this transition occurs, and to characterise the long-term dynamics of small bodies in the intermediate regime: relevant resonances, chaotic zones (if any), and timescales at play.

*Methods. *The different regimes are explored analytically and numerically. We also monitored the behaviour of swarms of particles during 4.5 Gyrs in order to identify which of the dynamical features are discernible in a realistic amount of time.

*Results. *There exists a tilted equilibrium plane (Laplace plane) about which orbits precess. The dynamics is integrable in the low and high semi-major axis regimes, but mostly chaotic in between. From about 800 to 1100 astronomical units (au), the chaos covers almost all the eccentricity range. The diffusion timescales are large, but not to the point of being indiscernible in a 4.5 Gyrs duration: the perihelion distance can actually vary from tens to hundreds of au. Orbital variations are damped near the ecliptic (where previous studies focussed), but favoured in specific ranges of inclination corresponding to well-defined resonances. Moreover, starting from uniform distributions, the orbital angles cluster after 4.5 Gyrs for semi-major axes larger than 500 au, because of a very slow differential precession.

*Conclusions. *Even if it is characterised by very long timescales, the inert Oort cloud mostly features chaotic regions; it is therefore much less inert than it appears. Orbits can be considered inert over 4.5 Gyrs only in small portions of the space of orbital elements, which include (90377) Sedna and 2012VP113. Effects of the galactic tides are discernible down to semi-major axes of about 500 au. We advocate including the galactic tides in simulations of distant trans-Neptunian objects, especially when studying the formation of detached bodies or the clustering of orbital elements.

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., 2017b, a, 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, 2017b), showed that the Lidov-Kozai mechanism raised by the giant planets 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 VP113 (Trujillo & Sheppard, 2014), and 2015 TG387 (Sheppard et al., 2019). After the discovery of 2012 VP113, 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 implications of 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.

2 Unified model of planets and galactic tides

We consider a 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 part111Even though the orbits of distant bodies are essentially barycentric, and not heliocentric, the Hamiltonian function is simpler when expressed in heliocentric coordinates. This is not a problem because we then use average coordinates, in which barycentric and heliocentric elements are equivalent (the wobbles of the Sun are averaged out)., a perturbation due to the planets, and a perturbation due to the galactic tides:

[TABLE]

Expressed using Keplerian elements, the two-body part is

[TABLE]

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

We assume that the small body never goes inside the orbits of the planets. The Hamiltonian $\varepsilon_{\mathrm{P}}\mathcal{H}_{\mathrm{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 than the 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

[TABLE]

These terms correspond to the monopole (index [math]), 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:

[TABLE]

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 $\mu_{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^{\prime},Y^{\prime},Z^{\prime})$ the coordinates of the small body in an analogous reference frame, but for which at any time the $X^{\prime}$ 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

[TABLE]

where the time derivative of $\theta$ is a constant and corresponds to the angular velocity of the galactic centre seen from the Sun. In the following, we write it $\nu_{\mathrm{G}}$. In the quadrupolar approximation, the Hamiltonian function describing the orbital perturbations of the small body from the galactic tides can be written

[TABLE]

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_{\theta}$ is conjugate to the angle $\theta$; it has been introduced such that the Hamiltonian function is autonomous. Using the usual approximation $\mathcal{G}_{2}=-\mathcal{G}_{1}$, we get

[TABLE]

where

[TABLE]

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

[TABLE]

We will now introduce explicit expressions for the small parameters:

[TABLE]

We note $(e,I,\omega,\Omega)$ the Keplerian elements of the small body in the ecliptic reference frame, with $e$ its eccentricity, $I$ its inclination, $\omega$ its argument of perihelion, and $\Omega$ its longitude of ascending node. We will use the subscript $\mathrm{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

[TABLE]

and

[TABLE]

We write $\psi$ the inclination of the ecliptic plane in the galactic reference frame, and $\alpha$ its ascending node. Since we have neglected the precession of the ecliptic pole, the angles $\psi$ and $\alpha$ are constant. The ascending node of the ecliptic can therefore be used as the origin of longitudes in the galactic frame, meaning that $\alpha\equiv 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 $\bar{\mathcal{H}}_{\mathrm{P}_{2}}$, $\bar{\mathcal{H}}_{\mathrm{P}_{4}}$, $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$, and $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{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\sim 600$ au, the planetary perturbations dominate over the galactic tides by more than a factor 10. The situation is reversed beyond $a\sim 1500$ au. In between, both kinds of perturbations have the same order of magnitude ($\varepsilon_{\mathrm{P}_{2}}$ and $\varepsilon_{\mathrm{G}_{\mathrm{V}}}$ cross at $a\sim 950$ au). However, since the eccentricity appears at the denominator in $\bar{\mathcal{H}}_{\mathrm{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

[TABLE]

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:

[TABLE]

where $L=\sqrt{\mu a}$ is a constant. We note that $g$ only appears in $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$, whereas $h$ only appears in $\bar{\mathcal{H}}_{\mathrm{P}_{2}}$ (through $\cos I$). 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_{2}$ flattening of its 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 $\Omega_{\mathrm{G}}=0$ and $\pi$, among which the stable ones correspond to the classical Laplace plane (located equivalently at $\Omega=\pi$ and [math]). We also recognise the ‘circular orthogonal’ equilibrium for $I_{\mathrm{G}}=I=90^{\text{o}}$ and $\Omega_{\mathrm{G}}=\Omega=\pm\pi/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 has no consequence here considering the timescales involved (see Appendix B for details).

Low-eccentricity orbits initially lying close to the ecliptic ($I_{\mathrm{G}}\approx\psi$ and $\Omega_{\mathrm{G}}\approx 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_{\mathrm{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 $\varepsilon_{\mathrm{P}_{2}}$ and $\varepsilon_{\mathrm{G}_{\mathrm{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 to the 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$ au and $13400$ au. Below $350$ au, the precession is about the ecliptic pole ($I\approx const$), so that a set of orbits initially lying close to the ecliptic plane never departs from it. Beyond $13400$ au, the precession is about the galactic pole ($I_{\mathrm{G}}\approx const$), such that orbits explore all the values of $I$ between $\psi-I_{\mathrm{G}}$ and $\psi+I_{\mathrm{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 $\Omega_{\mathrm{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 $\Omega_{\mathrm{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.

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 (see Fig. 2), the first step is to understand the two kinds of dynamics taken separately. Both $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{P}_{2}}$ and $\varepsilon_{\mathrm{G}_{\mathrm{V}}}\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{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 $\varepsilon_{\mathrm{G}_{\mathrm{V}}}\ll\varepsilon_{\mathrm{P}_{2}}$, the dynamics is largely dominated by the planetary perturbations. Expressed in ecliptic coordinates, the Hamiltonian $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{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

[TABLE]

The angle $\omega$ increases for $I<63^{\text{o}}$, decreases for $63^{\text{o}}<I<117^{\text{o}}$, and increases again for $I>117^{\text{o}}$. In contrast, $\Omega$ decreases for $I<90^{\text{o}}$ and increases for $I>90^{\text{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 $\omega$ appear at $I\approx 63^{\text{o}}$ and $117^{\text{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_{\mathrm{G}}$, $\omega_{\mathrm{G}}$ and $\Omega_{\mathrm{G}}$ driven by $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{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 $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{P}_{2}}$ (see Eq. 11) acts as the integrable dominant part of $\mathcal{F}$, while the galactic component $\varepsilon_{\mathrm{G}_{\mathrm{V}}}\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$ (see Eq. 12) acts as a small perturbation.

Since our dominant part $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{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 $\varepsilon_{\mathrm{G}_{\mathrm{V}}}\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$ is composed of several terms featuring various combinations of $\omega$ and $\Omega$ (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 $\varepsilon_{\mathrm{G}_{\mathrm{V}}}$), 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^{\text{o}}$ are obtained by replacing $\cos I$ by $-\cos I$ and $\Omega$ by $-\Omega$. As shown in the top panel of Figs. 7-8, the libration zone of $\Omega$ 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 $\omega-\Omega$ and $2\omega-\Omega$ are the first ones to overlap when the galactic tides increase (yellow and grey areas). As shown in the bottom panel of Figs. 7-8, the resonances $\omega+\Omega$ and $2\omega+\Omega$ are by far the largest ones in perihelion distance. The other resonances are quite small in comparison, and the libration zone of $\Omega$ even has a null width in $q$. The resonances $\omega\pm 2\Omega$, visible in Fig. 6, do not even appear in Figs. 7-8: this means that they only exist at second order of $\varepsilon_{\mathrm{G}_{\mathrm{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-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.

4.3 Galactic regime

If $\varepsilon_{\mathrm{P}_{2}}\ll\varepsilon_{\mathrm{G}_{\mathrm{V}}}$, the dynamics is dominated by the galactic tides. The dynamics driven by $\varepsilon_{\mathrm{G}_{\mathrm{V}}}\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{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

[TABLE]

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 $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$ for different values of $K$. The limit $I_{\mathrm{G}}=0$ or $180^{\text{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^{2}=1-K^{2}$ (border of the forbidden regions in Fig. 9). The limit $e=0$ is a fixed point with circulating $\Omega_{\mathrm{G}}$, but it is unstable for $K^{2}<4/5$, that is, for $27^{\text{o}}<I_{\mathrm{G}}<153^{\text{o}}$. For $K^{2}<4/5$, there are two additional fixed points located at $\omega_{\mathrm{G}}=\pi/2$ or $3\pi/2$ and

[TABLE]

which is equivalent to the condition

[TABLE]

These ones are stable, still with circulating $\Omega_{\mathrm{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_{\mathrm{G}}$, not the ecliptic one $I$. Moreover, $\Omega_{\mathrm{G}}$ is always decreasing if $I_{\mathrm{G}}<90^{\text{o}}$ and always increasing if $I_{\mathrm{G}}>90^{\text{o}}$ (the period of its linear part, already mentioned in Sect. 3, is noted $P_{\Omega^{*}}$ by Higuchi et al., 2007).

4.4 Intermediate non-integrable regime

In the intermediate regime (say, from $a\sim 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. (2017a). 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 $\omega$ and $\Omega$ circulate (see Sect. 4.1). However, the two degrees of freedom now interact, meaning that genuine resonances appear between $\omega$ and $\Omega$ (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 $\omega+\Omega$ and $2\omega+\Omega$ 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 $\omega-\Omega$ and $2\omega-\Omega$ for retrograde orbits. As expected, we also observe libration islands of $\omega$ at $I\approx 63^{\text{o}}$ and $117^{\text{o}}$, and libration islands of $\Omega$ at $I\approx 90^{\text{o}}$. Hence, the $90^{\text{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 $\varepsilon_{\mathrm{G}_{\mathrm{V}}}^{2}$ or more.

When we increase the semi-major axis, chaos spreads near the separatrices of the main resonances and libration islands of $\omega$ and $\Omega$. 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 $\omega$ and $\Omega$ 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^{2}$ 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 $\Omega$ stays around [math] or $\pi$, 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 $\bar{\mathcal{H}}_{\mathrm{P}_{2}}$ diverges and makes the planetary perturbation dominate again (see Eq. 11).

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$ VP113, and $2015$ TG387) are known to be stable, even though Sheppard et al. (2019) mentioned that $2015$ TG387 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$ VP113 ($a\approx 270$ au) show that $e$ and $I$ are almost constant while $\omega$ and $\Omega$ circulate. This is also the case for Sedna ($a\approx 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$ TG387 is more interesting ($a\approx 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$ TG387 is not far from a chaotic zone surrounding the $\omega+\Omega$ 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$ TG387 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^{\text{o}}$, we obtain the left panel of Fig. 16. The situation is similar to that of $2015$ TG387, meaning that the trajectory of P9 is regular but close to a chaotic zone emerging from the $\omega+\Omega$ 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.

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 $\Omega$ near the ascending node of the galactic plane (here located at $\Omega=\pi$). Moreover, we found in Sect. 4.4 that the combination $\varpi=\omega+\Omega$ is among the strongest resonances in the transitional regime, which should preferentially orient $\varpi$ around $\pm\pi/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$ au and $2000$ au. In order to manipulate easily understandable variables, we build our initial conditions in ecliptic coordinates. The angles $\omega$ and $\Omega$ are set uniformly between [math] and $2\pi$, and we structure our exploration in slices of perihelion distance (e.g. $q\in[40,60]$ au, $[60,80]$ au, etc.) and slices of ecliptic inclination cosine (e.g. $\cos I\in[0.9,1]$, $[0.8,0.9]$, etc.). Each slice is uniformly populated by a sample of $10^{5}$ 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 $\omega$ and $\Omega$ at large semi-major axes. As illustrated in Fig. 17, the extension and shape of these regions are different according 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 $\omega$ and $\Omega$; 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 line333Measured on the ecliptic, the longitude of ascending node of the proposed P9 is close to the longitude of descending node of the galactic plane. The galactic Laplace plane mentioned in Sect. 3 is therefore rotated by $180^{\text{o}}$ with respect to the Laplace plane raised by P9, and it has nothing to do with the clustering of orbital planes mentioned for instance by Batygin et al. (2019)..

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 $\omega$ and $\Omega$ 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\in[45^{\text{o}},55^{\text{o}}]$ and $[125^{\text{o}},135^{\text{o}}]$, which corresponds to the location of the two pairs of resonances $(\omega+\Omega,2\omega+\Omega)$ and $(\omega-\Omega,2\omega-\Omega)$, 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^{\text{o}}$ (and $180^{\text{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\sim 90^{\text{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\sim 90^{\text{o}}$ produces large oscillations of the inclination (Sect. 4.2).

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$ VP113 remain efficiently fossilised since the early stages of the solar system evolution?

For each value of $(a,q,I)$, we now look in the plane $(\omega,\Omega)$ for the initial condition producing the maximum variation of $q$ or $I$ in $4.5$ Gyrs. To this end, we save the extrema $q_{\mathrm{min}}$ and $q_{\mathrm{max}}$ reached by $q$ (resp. $I$) in the course of each numerical integration, and we apply an optimisation algorithm to maximise their difference $\Delta q=q_{\mathrm{max}}-q_{\mathrm{min}}$ (resp. $\Delta 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 $(\omega,\Omega)$ that maximises $\Delta q$ is generally different from the one that maximises $\Delta 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 in the $(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-23 delimits the ‘inert’ portion of the space, defined arbitrarily as $\Delta q<10$ au or $\Delta I<5^{\text{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:

1. a)

$a\gtrsim 500$ au and the orbit is nearly circular. 2. b)

$500\lesssim a\lesssim 1500$ au and $q$ is close to the planetary region. 3. c)

$a\lesssim 500$ au.

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

1. 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). 2. 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. 3. c)

For $a\lesssim 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\approx 540$ au) and $2012$ VP113 ($a\approx 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 $\omega$ and $\Omega$ (see Sect. 4.2). As discussed in Sect. 5.2, this complex structure produces a very marked differential spreading of small bodies in the space of orbital element. This structure disappears for large perihelion distances, as all resonances overlap. Sedna and $2012$ VP113 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^{\text{o}}$ and $24^{\text{o}}$), both of them are out of any of the main resonances. This makes them true ‘inert’ objects with precessing orbits. As expected from Sect. 5.1, on the contrary, $2015$ TG387 is out of the inert region depicted in Fig. 22-23 ($a\approx 1190$ au, $q\approx 65$ au, $I\approx 12^{\text{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_{2}$ 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 means that 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$ VP113), 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$ Gyrs, 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^{\text{o}}$ in inclination. Interestingly, the largest changes of perihelion distance are reached for ecliptic inclinations $I$ in the ranges $[45^{\text{o}},55^{\text{o}}]$ and $[125^{\text{o}},135^{\text{o}}]$. These ranges are delimited by the two pairs of strong resonances $(\omega+\Omega,2\omega+\Omega)$ and $(2\omega-\Omega,\omega-\Omega)$, 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 ($\Delta q<10$ au and $\Delta I<5^{\text{o}}$ over $4.5$ Gyrs), we find that it is remarkably small. The precise limits of the inert region can be found in Fig. 22-23. It is either composed of: a) nearly circular orbits with $a\gtrsim 500$ au, b) orbits with $500\lesssim a\lesssim 1500$ au and perihelion close to the planetary region, or c) orbits with $a\lesssim 500$ au (as long as they are unaffected by mean-motion resonances). Moreover, orbits are truly inert only if their perihelion 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.

Acknowledgements.

We thank Vacheslav Emel’yanenko for inspiring discussions about this manuscript, as well as for his independent verifications 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 galactic reference 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 ecliptic444This is the same reference frame as used by Higuchi et al. (2007); the minus sign in their Eq. (22) is a typographical error.. We write $(I_{\mathrm{G}},\omega_{\mathrm{G}},\Omega_{\mathrm{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,\omega,\Omega)$ the Keplerian elements of the small body measured in this reference frame. Passing from one frame to another corresponds to a rotation of $\pm\psi$ around the first axis, $\psi$ being the inclination of the ecliptic measured in the galactic reference frame. We obtain

[TABLE]

[TABLE]

and

[TABLE]

The inverse formulas are obtained by replacing $\psi$ by $-\psi$. Using these expressions, the Hamiltonian $\mathcal{F}$ in Eq. (13) can be written both in ecliptic or in galactic coordinates. More precisely, we have

[TABLE]

and

[TABLE]

where $C\equiv\cos\psi$ and $S\equiv\sin\psi$.

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

[TABLE]

where

[TABLE]

and

[TABLE]

In these expressions, $\mu_{\mathrm{P}}$, $J_{2}$, and $R_{\mathrm{P}}$ are the gravitational parameter, the flattening coefficient, and the equatorial radius of the host planet, respectively, whereas $\mu_{\odot}$, $a_{\odot}$, and $e_{\odot}$ 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 $\mathrm{G}$ index refers to the orbital plane of the Sun.

The overall Hamiltonian function in Eq. (24) 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_{2}$ flattening of the central body (see e.g. Tremaine et al. 2009). As such, $\varepsilon_{\mathrm{J}}$ is strictly equivalent to the parameter $\varepsilon_{\mathrm{P_{4}}}$ used above (see Eq. 10), and $\bar{\mathcal{K}}_{\mathrm{J}}$ is identical to $\bar{\mathcal{H}}_{\mathrm{P_{4}}}$ (see Eq. 11). Furthermore, we see here that $\varepsilon_{\odot}$ has the same $a^{2}$ multiplier as $\varepsilon_{\mathrm{G}_{\mathrm{V}}}$ (see Eq. 10), and that $\bar{\mathcal{K}}_{\odot}$ has nearly the same form as $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$ (see Eq. 12), apart from the additional term $-e^{2}/4$. The two problems are therefore not 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^{\text{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^{\text{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. 24). 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

[TABLE]

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.

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_{\mathrm{G}},\Omega_{\mathrm{G}})$ plane for different values of $(e,\omega_{\mathrm{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 $\Omega_{\mathrm{G}}=0$ and $\pi$, denoting the classic Laplace plane, are shifted in inclination. For example, Figs. 25 and 26 show the inclination of this ‘instantaneous Laplace plane’ with respect to $(e,\omega_{\mathrm{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\rightarrow 1$). This property can be seen in Fig. 27, where we plot the inclination of this plane averaged over $\omega_{\mathrm{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.

Since the known distant trans-Neptunian objects are very eccentric, we expect from Fig. 27 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}=\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{P}_{2}}+\varepsilon_{\mathrm{G}_{\mathrm{V}}}\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$, where the expressions for each part can be found in Eqs. (22) and (23). If the planetary perturbations dominate over the galactic tides, $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{P}_{2}}$ acts as the integrable dominant part, whereas $\varepsilon_{\mathrm{G}_{\mathrm{V}}}\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{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}(\varepsilon_{\mathrm{G}_{\mathrm{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 $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{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 $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{P}_{2}}$ (see Eq. 15 and Fig. 6).

ii) The second term of $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$ is factored by $\cos\Omega$. 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 $\varepsilon_{\mathrm{G}_{\mathrm{V}}}=0$, but about a tilted pole. By averaging over $\omega$, 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,\Omega)$.

iii) All the remaining terms of $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$ correspond to resonances and libration zones for $\omega$ and $\Omega$. 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 $\omega+\Omega$. From the ecliptic Delaunay coordinates $(P_{\omega},P_{\Omega},\omega,\Omega)$, the corresponding change of coordinates is

[TABLE]

After averaging over the circulating angle $\gamma$, we end up with one constant of motion $\Gamma$, and one degree of freedom $(\Sigma,\sigma)$. Table 2 gives the constants associated to all resonances that appear at first order in $\varepsilon_{\mathrm{G}_{\mathrm{V}}}$ (that is, the ones directly appearing in the expression of $\bar{\mathcal{H}}_{\mathrm{G}_{\mathrm{V}}}$ in Eq. 23). The resonance centre is mostly governed by the unperturbed Hamiltonian $\varepsilon_{\mathrm{P}_{2}}\bar{\mathcal{H}}_{\mathrm{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 2.

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 $\Sigma_{0}$, keeping terms up to degree $2$ for the unperturbed part, and up to degree [math] for the perturbation. The resulting Hamiltonian has the form

[TABLE]

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

Figures 7 and 8 show the location and widths of all the resonances that appear at first order in $\varepsilon_{\mathrm{G}_{\mathrm{V}}}$, 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. $\Sigma$), but, as shown by Figs. 7-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 $\omega$ 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-23 show sections of this space in the $(a,q)$ and $(a,I)$ directions. For completeness, Fig. 28 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. 28 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. 28, 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^{\text{o}}$), while producing the generalised chaotic region studied in Sect. 4.4.