Issue 
A&A
Volume 629, September 2019



Article Number  A95  
Number of page(s)  20  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201936298  
Published online  12 September 2019 
Chaos in the inert Oort cloud
^{1}
IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Université de Lille,
75014
Paris, France
email: melaine.saillenfest@obspm.fr
^{2}
National Astronomical Observatory of Japan,
2211 Osawa, Mitaka,
Tokyo 1818588, Japan
Received:
12
July
2019
Accepted:
9
August
2019
Context. Distant transNeptunian 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 longterm 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 semimajor 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 welldefined resonances. Moreover, starting from uniform distributions, the orbital angles cluster after 4.5 Gyrs for semimajor 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 semimajor axes of about 500 au. We advocate including the galactic tides in simulations of distant transNeptunian objects, especially when studying the formation of detached bodies or the clustering of orbital elements.
Key words: celestial mechanics / Kuiper belt: general / Oort cloud
© M. Saillenfest et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 transNeptunian or Kuiper belt objects (with all their subclasses) and the Oort cloud. This distinction was made because the transNeptunian population is indeed observed on orbits lying beyond or close to Neptune, whereas the longperiod 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 transNeptunian 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 transNeptunian 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 transNeptunian 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 semimajor 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 semimajor axis (Gallardo et al. 2012), because energy kicks result in larger variations of semimajor axis if the semimajor 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 semimajor 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 meanmotion 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 semimajor axis remains at the resonance location. This mechanism, however, is only efficient for semimajor 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 semimajor axis should be smaller than 500 au for the planets to possibly have a substantial effect through meanmotion resonances. As regards the effects of the galactic tides, Fouchard et al. (2017) showed that an object with perihelion in the JupiterSaturn region, that is, below 15 au from the Sun, should have a semimajor 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 semimajor axis is smaller than 1600 au and the perihelion distance is larger than 45 au, but the semimajor 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 semimajor 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_{113} (Trujillo & Sheppard 2014), and 2015 TG_{387} (Sheppard et al. 2019). After the discovery of 2012 VP_{113}, 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 longterm 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 longterm 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 mixedtype 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 meanmotion 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 Sunbody Keplerian part^{1}, a perturbation due to the planets, and a perturbation due to the galactic tides: (1)
Expressed using Keplerian elements, the twobody part is (2)
where a is the semimajor 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 can therefore be expanded in Legendre polynomials. As explained in the Introduction, meanmotion 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 (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: (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 . We call this frame the “ecliptic” reference frame. The quantities μ_{i} and a_{i} are the gravitational parameter and the semimajor axis of the planet i, for a total of N planets. Because of their small semimajor 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 lowestorder 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 (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 (6)
where , and 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 , we get (7)
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 (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 semimajor axis (that we still denote a). Dropping the constant parts, the secular Hamiltonian is (9)
We will now introduce explicit expressions for the small parameters: (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 (11)
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).
Values used for the physical constants of the problem.
3 The galactic Laplace plane
The explicit expressions of the small parameters (Eq. (10)) have been chosen such that the Hamiltonian functions , , , and have the same order of magnitude for e = 0. The secular semimajor 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 ( and cross at a ~ 950 au). However, since the eccentricity appears at the denominator in (see Eq. (11)), we expect that the planetary perturbations always have a substantial effect in the higheccentricity 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 (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: (14)
where is a constant. We note that g only appears in , whereas h only appears in (through cosI). In this context, the galactic coordinates are therefore the most natural coordinates to use.
First of all, we note that the Hamiltonian (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 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 transNeptunian 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).
Loweccentricity 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 (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 semimajor 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 and 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 semimajor 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 disclike 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 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 transNeptunian 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 transNeptunian 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 transNeptunian 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 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 semimajor 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 (Eq. (13)) for a circular orbit. The semimajor 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 semimajor 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 (Eq. (13)) is composed of two parts that dominate respectively in the low and highsemimajor axis regimes (seeFig. 2), the first step is to understand the two kinds of dynamics taken separately. Both and 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 , the dynamics is largely dominated by the planetary perturbations. Expressed in ecliptic coordinates, the Hamiltonian 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 (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 transNeptunian 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 semimajor axes, represents a very small variation of eccentricity. When expressed in galactic coordinates, the evolution of I_{G}, ω_{G} and Ω_{G} driven by 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 semimajor axes and/or high eccentricities), the effects of galactic tides can be studied in a perturbative approach: the planetary component (see Eq. (11)) acts as the integrable dominant part of , while the galactic component (see Eq. (12)) acts as a small perturbation.
Since our dominant part is already expressed in actionangle coordinates (i.e. it does not depend on the angles), the perturbative approach is straightforward. Expressed in ecliptic coordinates, our perturbing part 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), 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 , 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: lowinclination 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 semimajor 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 and attain their maximum absolute value at I = 0^{o} and 180^{o}. The locations where the integer combinations 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 semimajor 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 , the dynamics is dominated by the galactic tides. The dynamics driven by 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 (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 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 largea 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 Ω_{G}, but it is unstable for K^{2} < 4∕5, that is, for 27^{o} < I_{G} < 153^{o}. For K^{2} < 4∕5, there are two additional fixed points located at ω_{G} = π∕2 or 3π∕2 and (17)
which is equivalent to the condition (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 by Higuchi et al. 2007).
Fig. 9
Level curves of the Hamiltonian function (Eq. (12)). The parameters chosen are K^{2} < 4∕5 (left) and K^{2} > 4∕5 (right). Top and bottom rows: the same level curves for two sets of variables. 
4.4 Intermediate nonintegrable 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 semimajor 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 or more.
When we increase the semimajor 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 semimajor 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 semimajor 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^{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 Ω stays around 0 or π, as expected from the Laplace plane (Sect. 3). It should be noted, however, that whatever the value of the semimajor axis, there is always a chaotic region for very eccentric orbits, and a stable region for even higher e, because the denominator of 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 semimajor 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 au^{2}yr^{−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 ≈ 63^{o} (Lidov–Kozai mechanism). Bottom: section for au^{2}yr^{−2} at ω_{G} = 0 decreasing. The chaotic bands are due to the overlap of the resonances ω − Ω and 2ω − Ω (for I < 90^{o}), or ω + Ω and 2ω + Ω (for I > 90^{o}), and the islands located at I ≈ 90^{o} are due to librations of Ω (orthogonal Laplace plane). 
Fig. 11
Same as Fig. 10 for a = 600 au. Left: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. From top to bottom, the largest islands are due to the resonance 2ω + Ω, the libration of ω at I ≈ 63^{o}, and the resonance 2ω − Ω. Right: section for au^{2}yr^{−2} at ω_{G} = 0 decreasing. The islands located at I ≈ 90^{o} are due to librations of Ω. 
Fig. 12
Same as Fig. 10 for a = 800 au (left) and a = 1100 au (right). Left:section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. Right: section for au^{2}yr^{−2} at ω_{G} = 0 decreasing. 
5 Implications for real objects
5.1 Dynamics of known inertOortcloud bodies
The dynamics of the three observed members of the inert Oort cloud (Sedna, 2012 VP_{113}, and 2015 TG_{387}) are known to be stable, even though Sheppard et al. (2019) mentioned that 2015 TG_{387} 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_{113} (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 semimajor 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_{387} is more interesting (a ≈ 1190 au), because its high semimajor axis corresponds to the region where both the planetary and galactic perturbations have substantial effects. Figure 15 shows that 2015 TG_{387} is not far from a chaotic zone surrounding the ω + Ω resonance. Its trajectory, however, is strictly quasiperiodic in our simplified model. Increasing the value of its semimajor axis in the uncertainty range leads to faster oscillations of q with a larger amplitude (see Fig. 15, right panel), but 2015 TG_{387} 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_{387}, 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 semimajor 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 Neptunecrossing 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 transNeptunian objects. This remains true for the updated P9 orbit obtained by Batygin et al. (2019), even though it has a somewhat smaller semimajor axis.
Fig. 13
Same as Fig. 10 for a = 1300 au. Left: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. The upper islands are due to librations of ω_{G} with e^{2} ≈ 1 − K^{2}. Right: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. The islands are due to librations of ω_{G} and oscillations of e^{2} around . 
Fig. 14
Same as Fig. 10 for a = 2000 au. Left: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. Right: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. The big islands are due to librations of ω_{G} and oscillations of e^{2} around . 
5.2 Evolution of a sample of objects over 4.5 Gyrs
The exploration of the inertOortcloud 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 transNeptunian 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 semimajor 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^{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 ω and Ω at large semimajor 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 semimajor 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 transNeptunian 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 semimajor 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^{2}.
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 semimajor 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 semimajor 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 semimajor 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 planetaryinduced 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 semimajor 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 semimajor 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 TG_{387}. The points of its trajectory that cross the section are represented in green. Left: nominal orbital elements given by Sheppard et al. (2019). Right: semimajor 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: semimajor 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 = 0^{o}, and ϖ = ω + Ω equally distributed in [0, 2π]; the initial value q_{i} 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 , Eq. (12)), and the second and third columns are for an evolution with both planets and galactic tides (Hamiltonian , 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 q_{i} in [30, 40] au, and the thirdone in [80, 100] au, while the cosine of the initial ecliptic inclination I_{i} 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 semimajor axis of Neptune (~ 30 au) for reference. 
Fig. 20
Minimum value of the semimajor 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.1width 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. 
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_{113} 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 semimajor 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 semimajor 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_{113} (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_{113} 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_{387} 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 longterm 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 meansthat for semimajor axes much smaller than this value (say 500 au), circular orbits precess about the ecliptic pole, whereas for semimajor 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 semimajor 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 semimajor axis value. For semimajor 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_{113}), 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 semimajor 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 semimajor 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 meanmotion 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 transNeptunian object with semimajor 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 5^{o} 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, cofunded 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^{3}. 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 (A.1) (A.2)
The inverse formulas are obtained by replacing ψ by − ψ. Using these expressions, the Hamiltonian in Eq. (13) can be written both in ecliptic or in galactic coordinates. More precisely, we have (A.4)
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 (B.1)
In these expressions, μ_{P}, J_{2}, 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 semimajor 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 transNeptunian 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_{2} flattening of the central body (see e.g. Tremaine et al. 2009). As such, ε_{J} is strictly equivalent to the parameter used above (see Eq. (10)), and is identical to (see Eq. (11)). Furthermore, we see here that ε_{⊙} has the same a^{2} multiplier as (see Eq. (10)), and that has nearly the same form as (see Eq. (12)), apart from the additional term − e^{2}∕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 transNeptunian 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 transNeptunian 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 transNeptunian case, we obtain the stability condition (B.4)
which is twice the analogous limit obtained in the satellite case. This corresponds to a semimajor 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 (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 semimajor 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 semimajor 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 semimajor 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.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 transNeptunian 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 , 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, acts as the integrable dominant part, whereas acts as a small perturbation.
Luckily, the unperturbed part is already expressed in actionangle coordinates. This means that, neglecting terms in , the longterm behaviour of the system is simply given by the average of over the nonresonant 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 does not include the angles; its acts therefore only as a small modulation of the precession velocities and governed by (see Eq. (15) and Fig. 6).
 (ii)
The second term of 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 , 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 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 (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 (that is, the ones directly appearing in the expression of in Eq. (A.5)). The resonance centre is mostly governed by the unperturbed Hamiltonian : 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.
Resonance centres and constants of motion arising from the dynamics in the vicinity of each resonance appearing at first order in .
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 Σ_{0}, keeping terms up to degree 2 for the unperturbed part, and upto degree 0 for the perturbation. The resulting Hamiltonian has the form (D.2)
which is a pendulum of centre Σ_{0} and half width . 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 , 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 semimajor 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
 Bannister, M. T., Shankman, C., Volk, K., et al. 2017, AJ, 153, 262 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Brown, M. E. 2016, AJ, 151, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Adams, F. C., Brown, M. E., & Becker, J. C. 2019, Phys. Rep., 805, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, J. C., Adams, F. C., Khain, T., Hamilton, S. J., & Gerdes, D. 2017, AJ, 154, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Brasser, R., Duncan, M. J., Levison, H. F., Schwamb, M. E., & Brown, M. E. 2012, Icarus, 217, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Breiter, S., & Ratajczak, R. 2005, MNRAS, 364, 1222 [NASA ADS] [CrossRef] [Google Scholar]
 Bretagnon, P. 1982, A&A, 114, 278 [NASA ADS] [Google Scholar]
 Brown, M. E., Trujillo, C., & Rabinowitz, D. 2004, ApJ, 617, 645 [NASA ADS] [CrossRef] [Google Scholar]
 Dones, L., Weissman, P. R., Levison, H. F., & Duncan, M. J. 2004, in Star Formation in the Interstellar Medium: In Honor of David Hollenbach, eds. D. Johnstone, F. C. Adams, D. N. C. Lin, D. A. Neufeeld, & E. C. Ostriker, ASP Conf. Ser., 323, 371 [NASA ADS] [Google Scholar]
 Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2016, A&A, 587, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fouchard, M. 2004, MNRAS, 349, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Fouchard, M., Rickman, H., Froeschlé, C., & Valsecchi, G. B. 2017, Icarus, 292, 218 [NASA ADS] [CrossRef] [Google Scholar]
 Fouchard, M., Higuchi, A., Ito, T., & Maquet, L. 2018, A&A, 620, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallardo, T., Hugo, G., & Pais, P. 2012, Icarus, 220, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Gladman, B., Holman, M., Grav, T., et al. 2002, Icarus, 157, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R., Gallardo, T., Fernández, J., & Brunini, A. 2005, Celest. Mech. Dyn. Astron., 91, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R. S., Soares, J. S., & Brasser, R. 2015, Icarus, 258, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Higuchi, A., & Kokubo, E. 2015, AJ, 150, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Higuchi, A., Kokubo, E., Kinoshita, H., & Mukai, T. 2007, AJ, 134, 1693 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Jílková, L., Portegies Zwart, S., Pijloo, T., & Hammer, M. 2015, MNRAS, 453, 3157 [Google Scholar]
 Kaib, N. A., Pike, R., Lawler, S., et al. 2019, AJ, 158, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, G., Naoz, S., Holman, M., & Loeb, A. 2014, ApJ, 791, 86 [Google Scholar]
 Morbidelli, A.,& Nesvorny, D. 2019, ArXiv eprints [arXiv:1904.02980] [Google Scholar]
 Murray, C. A. 1989, A&A, 218, 325 [NASA ADS] [Google Scholar]
 Poli, R., Kennedy, J., & Blackwell, T. 2007, Swarm Intell., 1, 33 [CrossRef] [Google Scholar]
 Polycarpe, W., Saillenfest, M., Lainey, V., et al. 2018, A&A, 619, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2016, Celest. Mech. Dyn. Astron., 126, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2017a, Celest. Mech. Dyn. Astron., 127, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2017b, Celest. Mech. Dyn. Astron., 129, 329 [Google Scholar]
 Sheppard, S. S., Trujillo, C. A., Tholen, D. J., & Kaib, N. 2019, AJ, 157, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Vokrouhlický, D., Nesvorný, D., & Dones, L. 2019, AJ, 157, 181 [NASA ADS] [CrossRef] [Google Scholar]
Even 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).
Measured 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^{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).
This is the same reference frame as used by Higuchi et al. (2007); the minus sign in their Eq. (22) is a typographical error.
All Tables
Resonance centres and constants of motion arising from the dynamics in the vicinity of each resonance appearing at first order in .
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 meanmotion 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 semimajor 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 (Eq. (13)) for a circular orbit. The semimajor 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 semimajor 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 and attain their maximum absolute value at I = 0^{o} and 180^{o}. The locations where the integer combinations 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 semimajor 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 (Eq. (12)). The parameters chosen are K^{2} < 4∕5 (left) and K^{2} > 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 semimajor 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 au^{2}yr^{−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 ≈ 63^{o} (Lidov–Kozai mechanism). Bottom: section for au^{2}yr^{−2} at ω_{G} = 0 decreasing. The chaotic bands are due to the overlap of the resonances ω − Ω and 2ω − Ω (for I < 90^{o}), or ω + Ω and 2ω + Ω (for I > 90^{o}), and the islands located at I ≈ 90^{o} are due to librations of Ω (orthogonal Laplace plane). 

In the text 
Fig. 11
Same as Fig. 10 for a = 600 au. Left: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. From top to bottom, the largest islands are due to the resonance 2ω + Ω, the libration of ω at I ≈ 63^{o}, and the resonance 2ω − Ω. Right: section for au^{2}yr^{−2} at ω_{G} = 0 decreasing. The islands located at I ≈ 90^{o} 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 au^{2}yr^{−2} at Ω_{G} = 0 decreasing. Right: section for au^{2}yr^{−2} at ω_{G} = 0 decreasing. 

In the text 
Fig. 13
Same as Fig. 10 for a = 1300 au. Left: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. The upper islands are due to librations of ω_{G} with e^{2} ≈ 1 − K^{2}. Right: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. The islands are due to librations of ω_{G} and oscillations of e^{2} around . 

In the text 
Fig. 14
Same as Fig. 10 for a = 2000 au. Left: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. Right: section for au^{2}yr^{−2} at Ω_{G} = 0 decreasing. The big islands are due to librations of ω_{G} and oscillations of e^{2} around . 

In the text 
Fig. 15
Same as Fig. 10 for the parameters of 2015 TG_{387}. The points of its trajectory that cross the section are represented in green. Left: nominal orbital elements given by Sheppard et al. (2019). Right: semimajor 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: semimajor 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 = 0^{o}, and ϖ = ω + Ω equally distributed in [0, 2π]; the initial value q_{i} 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 , Eq. (12)), and the second and third columns are for an evolution with both planets and galactic tides (Hamiltonian , 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 q_{i} in [30, 40] au, and the thirdone in [80, 100] au, while the cosine of the initial ecliptic inclination I_{i} 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 semimajor axis of Neptune (~ 30 au) for reference. 

In the text 
Fig. 20
Minimum value of the semimajor 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.1width 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 5^{o} 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 semimajor 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 semimajor 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.