Issue 
A&A
Volume 677, September 2023



Article Number  A83  
Number of page(s)  17  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202244982  
Published online  08 September 2023 
Dynamics of transNeptunian objects near the 3/1 meanmotion resonance with Neptune
^{1}
CFisUC, Departamento de Física, Universidade de Coimbra,
3004516
Coimbra, Portugal
email: acor@uc.pt
^{2}
Instituto de Telecomunicações, Universidade de Aveiro,
3810193
Aveiro, Portugal
^{3}
IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université,
77 Av. DenfertRochereau,
75014
Paris, France
Received:
15
September
2022
Accepted:
14
June
2023
The complex classification of transNeptunian objects (TNOs) that are captured in meanmotion resonances (MMRs) and the constraint of their multiple origins are two significant open problems concerning the Solar System. The casebycase study of the different MMRs and their characteristics provide information about their origin and dynamics, which helps us to understand the early stages of the Solar System evolution. In this paper, we study the dynamics of the detected TNOs close to a 3/1 MMR with Neptune. We initially use a semianalytic threebody model to investigate the coplanar secular dynamics of these objects and find the stationary points. We then use surface sections and stability maps to analyse the nonaveraged dynamics. These methods allow us to isolate the different stability regions and determine the extent of the chaotic regions. We show that stability maps are an extremely powerful tool for studying the resonant dynamics when they are computed in terms of the resonant angle. We then use these maps to study the nonplanar threebody problem and the full dynamics in the presence of planetary perturbations. We confirm that TNOs near the 3/1 MMR regions can exist at very high inclinations. In the framework of the threebody problem, many of these objects can also be stable outside the 3/1 MMR owing to a Kozai secular resonance. However, when we take into account the perturbations of the four giant planets, the Kozai regions disappear and only the 3/1 MMR region remains, with eccentricities e ≲ 0.5.
Key words: celestial mechanics / Kuiper belt: general / planets and satellites: dynamical evolution and stability
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The physical mechanisms that lead to the presently observed characteristics of transNeptunian objects (TNOs) include phenomena that occurred in the past, such as the outward migration of Neptune initially proposed by Fernandez & Ip (1984). One of the most accepted scenarios is the Nice model (Tsiganis et al. 2005; Gomes et al. 2005; Morbidelli et al. 2007), where a great instability between the orbits of the giant planets scattered and ejected tens of thousands of planetesimals from their initial locations, which explains the modernday absence of a dense TNO population. Important factors in the formation of the TNOs are the composition of the primordial disc; the locations where the migration of Neptune started and ended; and the migration speed.
Numerical simulations of an outward migration of Neptune in a primordial planetesimal Kuiper Belt composed of a coldgas disc and particles with random semimajor axis between 20 and 80 au provide the capture of objects in a large number of meanmotion resonances (MMRs), with an overpopulation of the 2/1 and 3/2 MMRs (Hahn & Malhotra 2005). Indeed, several TNOs are currently observed close to MMR with Neptune (e.g. Elliot et al. 2005; Lawler et al. 2019; Bannister et al. 2018; Crompvoets et al. 2022). These can be resonant, but they can also be scattering, detached, or classical objects. Distinguishing between these possibilities is a delicate task, but the outcome has important implications for understanding the composition and the past evolution of these bodies (Gladman et al. 2008). For instance, by comparing the results of five different Neptune migration models with the high pericentre distance of detected TNOs close to 5/2 and 3/1 MMRs, Lawler et al. (2019) showed that a grainy and slow Neptune migration provides the best match to the observational data.
The complexity of MMRs has been exhaustively studied. Many authors have concentrated their efforts on analysing the topology of the problem and the equilibrium points. For the planar case, equilibrium solutions have been found for many MMRs (e.g. Message 1958; Taylor 1983a,b). Using a semianalytical secular model, Beauge (1994) demonstrated the existence of asymmetrical equilibria for the 2/1 and 3/1 MMRs, but their absence for the 3/2 and 4/3 MMRs. Voyatzis et al. (2018) determined the equilibrium points using periodic orbits, while Lan & Malhotra (2019) studied them with Poincaré surface sections. For the nonplanar case, equilibrium solutions were also obtained and studied (e.g. Ichtiaroglou et al. 1989; Hadjidemetriou 1993; Voyatzis & Kotoulas 2005; Gallardo 2019; Namouni & Morais 2020). Finally, for some resonances, it has been shown that the influence of the remaining giant planets (Jupiter, Saturn and Uranus) does not significantly modify the position of the resonant equilibria in the parameter space (e.g. Saillenfest & Lari 2017; Saillenfest et al. 2017; Malhotra et al. 2018; Lei et al. 2022). Nevertheless, perturbations from the giant planets can be important for the longterm evolution and capture inside the MMRs and for the formation of other TNO classes.
Hahn & Malhotra (2005) showed that a stirredup Kuiper Belt, which should extend at least until a = 55 au, is capable of capturing objects in the 3/1 MMR with Neptune. Bannister et al. (2018) report seven TNOs captured in a 3/1 MMR with Neptune, while Crompvoets et al. (2022) confirm a total of 12 objects inside this resonance. However, the objects with perihelion distances smaller than 36 au, which includes most of the TNOs close to the 3/1 MMR, are inside a chaotic region that provokes a diffusion in the orbital elements (Duncan et al. 1995; Gladman et al. 2002; Fernández et al. 2004). As a result, these bodies may be scattered and do not remain captured in resonance for a long period. The main characteristics of their dynamics should then be determined to evaluate whether or not these bodies are really trapped in the 3/1 MMR. Estimation of the amount of time spent by TNOs inside this resonance and determination of the mechanisms that lead these bodies to escape are also important, in order to compare the number of bodies expected in this region according to the current formation models with numbers derived from observations.
In this paper, we focus our analysis on the TNOs that are found close to the 3/1 MMR with Neptune, which are also known as threetinos. In Sect. 2, we present a dynamical analysis of these objects carried out in order to ascertain whether or not they are locked in resonance. In Sect. 3, we use a semianalytic model to describe the secular planar dynamics around the 3/1 MMR and determine the stationary points of the problem. We then investigate how the secular model is modified in the nonaveraged planar case using surface sections (Sect. 4) and stability maps (Sect. 5). These maps are then used to investigate the nonplanar dynamics and the perturbations from the remaining giant planets. Finally, we discuss the main results obtained and our conclusions in Sect. 6.
2 Dynamics of TNOs near the 3/1 MMR
In this section, we present a dynamical analysis of the already known TNOs close to a 3/1 MMR with Neptune. Using the IAU Minor Planet Center database^{1}, we identified 82 bodies with heliocentric semimajor axes between 61 and 64 au, which corresponds to orbital period ratios with Neptune of between 2.9 and 3.1. When converting to barycentric coordinates, the number of objects in this range drops to 76.
In Fig. 1, we plot the 76 selected objects in a diagram showing semimajor axis versus eccentricity, where each colour signifies a different inclination range. We note that most of the bodies have high eccentricities, in the range 0.3–0.5, and have moderate and high inclinations of up to 40°. Therefore, these objects are part of the scattering disc and detached populations. Seven objects have a sufficiently high eccentricity to cross the orbit of Neptune. Four of them, 2005 OE, 2010 GW_{64},2012 HD_{2}, and 2013 HS_{150}, have very eccentric orbits, very high inclinations, and a perihelion distance of approximately 3 and 10 au, and as a result are not classified as TNOs; they belong to the Damocloid class objects (Jewitt 2005), and are not considered in our analysis.
2.1 Frequency analysis
We first identified the bodies that can presently be captured inside the 3/1 MMR. To this end, we numerically integrated all TNOs over 100 Myr, taking into account the perturbations from the four giant planets^{2}. We used the symplectic integrator SABA4 of Laskar & Robutel (2001) and an integration step size of 0.1 yr. We then performed a frequency analysis of the heliocentric orbital elements (Laskar 1990, 1993) to determine the fundamental frequencies (n,v_{1},v_{2}) of the TNOs, where n is the main frequency associated with the mean longitude (i.e. it is the average mean motion, which corresponds to the mean value of the mean motion over 100 Myr), and v1 and v2 are the main secular frequencies associated with the longitude of the perihelion and node, respectively. Finally, we checked for combinations of these frequencies such that (1)
where n′ is Neptune’s average mean motion and k is an integer. The exact equality almost never occurs numerically, and so we consider that the previous condition is satisfied whenever (2)
When such a combination exists, the associated resonant angle (3)
is expected to librate around some equilibrium value, where λ, ϖ, and Ω are the mean longitude, the longitude of the perihelion, and the longitude of the node of the TNO, respectively, and λ′ is the mean longitude of Neptune.
In Table 1, we list the 15 objects where it was possible to identify such a combination, all for k = 0 (in Fig. 1 we mark these objects with a cross). We obtain the centre of libration, the libration period, and the amplitude of libration. The frequency analysis also allows us to determine the additional forcing frequencies and amplitudes. When the libration frequency can be clearly isolated from other perturbations, we classify the object as resonant stable (labelled with an ‘S’ in Table 1). However, in some cases the libration angle displays strong chaotic behaviour, and the resonance may be unstable (labelled with a ‘C’ in Table 1).
Finally, we additionally identify the main frequency of ϖ and Ω. In a decoupled problem, we expect that and . However, we observe that for one object, namely 2011US_{411}, both ϖ and Ω precess with the same fundamental frequency. As a result, the argument of the perihelion, (4)
also librates. In this case, all resonant angles σ_{j≠k} librate, because (Eq. (3)) (5)
Therefore, in Table 1, we identify this case with the notation ‘all’ for the resonant angle.
Crompvoets et al. (2022) report that objects 2013 UB_{17} and 2015 VM_{166} are also securely in a 3/1 MMR with Neptune, but these authors only integrate the orbits for 10 Myr. These two TNOs are included in our selection of 76 objects near the 3/1 MMR, but we were not able to confirm their longterm resonant nature. Indeed, these objects are initially resonant, but after a few million years they display chaotic evolution. We use a different method to analyse the orbits and also integrate for a ten times longer timespan. As a result, we are able to detect the longterm chaotic diffusion of the orbits more precisely. The initial orbital elements that we use may also slightly differ from those in previous studies, which can lead to a different orbital evolution. These differences show that the classification of a TNO as a resonant object is not an easy task, and must be done with caution.
In Fig. 2, we plot the orbital evolution of some TNOs near the 3/1 MMR that better illustrate the diversity of the observed behaviours. We first show an example of the many TNOs that are not captured in resonance, namely 2001 KG_{77}. We observe circulation for all resonant angles (Eq. (3)) and moderate chaotic evolution of the orbital elements. We then show three examples of captured TNOs: 2011US_{411} appears to be the most regular object in the family, as σ_{0} and ω both librate with small amplitude; 2015 KY_{173} has a nearly coplanar orbit and also a small libration amplitude, but only the angle σ_{0} librates; and finally, 2003 LG7 shows strong chaotic behaviour, although it remains in libration throughout the entire simulation. This latter object librates initially with large amplitude around 180°, but for some amount of time it librates around 100° or 260° with nearly half of the initial amplitude.
Fig. 1 Barycentric orbital elements distribution of the TNOs that are close to the 3/1 MMR with Neptune. We show 76 bodies in the range [61, 64] au, with each colour corresponding to a different inclination range. The resonant TNOs are marked with a cross (Table 1), and the remaining nonresonant with a dot. The orange vertical line gives the nominal resonant semimajor axis a ≈ 62.55 au. The black solid and dashed lines give the overlap with Neptune and Uranus orbits, respectively. 
TNOs in a 3/1 MMR with Neptune for at least 100 Myr.
2.2 Capture in resonance
Fifteen TNOs were identified in the 3/1 MMR with Neptune (Table 1), clearly suggesting some past orbital evolution. We therefore performed some quick numerical simulations of the outer migration of Neptune in the late stages of the evolution of the Solar System. Our simulations start after the occurrence of the instability due to the crossing of the 2/1 MMR resonance between Jupiter and Saturn, when Neptune had already evolved into an outer orbit with respect to Uranus (Tsiganis et al. 2005). We do not intend here to find the best migration model that reproduces the present distribution of bodies around the 3/1 MMR; our goal is simply to understand the dynamical states acquired by the TNOs when they cross this resonance. Therefore, we only include Neptune and the TNOs in the simulations.
Here, we adopt the migration model from Beaugé et al. (2006), which considers a Stoke’s drag force acting on Neptune’s orbit with an exponential decay of the semimajor axis with τ = 10^{7} yr. We fix the initial value of Neptune’s semimajor axis at 20 au, and vary the initial semimajor axes of the TNOs in the range 45–62 au with a step size of 1 au. The initial eccentricities and inclinations are in the ranges 0–0.25 and 0°–25°, respectively, with step sizes of 0.025 and 2.5°. The distribution of the angles λ, ϖ, and Ω is random in the interval [0°, 360°]. We therefore integrate a total of 2178 test particles over 40 Myr. Neptune reaches the present semimajor axis, a ≈ 30.15 au, around 18 Myr. Figure 3 shows some examples of TNOs captured in the 3/1 MMR.
We observe that 79 test particles were trapped in the 3/1 MMR at the end of the simulation (~4%). Most of the captures occurred for the resonant angle σ_{0}, but 12 also involved the angle ω, and 2 occurred for the angle σ_{2} alone. There were captures for all initial semimajor axis and inclination values, while for the initial eccentricity, capture was only observed when e > 0.1. For the TNOs captured in resonance, the final semimajor axis ends around the present value of 62.5 au, and the initial inclinations remain nearly unchanged. However, the initial eccentricity is always excited by the resonance and the final values range between 0.1 and 0.5. Figure 4 shows the distribution of the final eccentricity of the TNOs captured in the 3/1MMR as a function of the initial semimajor axis. We observe that higher eccentricities are only possible for the TNOs with smaller initial semimajor axes, because they are trapped earlier and therefore the eccentricity has more time to grow.
From our simple numerical simulation, we conclude that the outward migration of Neptune can explain most of the presently observed dynamical properties of the threetinos (Table 1). The resonant angle is mainly σ_{0}, the final eccentricities are mostly within 0.1 < e < 0.5, and the inclinations are random. We also observe that, in some simulations, σ_{0} librates with small amplitude around either 100° or 260°, while in a few it librates with large chaotic amplitudes around 180°. This last case is only observed for initial semimajor axes larger than 55 au. An interesting feature present in most captured simulations is that before the libration centre settles at 100^{°} or 260^{°} with small amplitude, the TNOs initially librate around 180^{°} with large amplitude. Therefore, the TNOs currently observed in this configuration were likely captured in the last stages of Neptune’s migration.
Fig. 2 Orbital evolution of some representative TNOs near the 3/1 resonance over 100 Myr. We show the semimajor axis, the eccentricity, the inclination, the resonant angle σ_{0}, and the argument of the pericentre (from top to bottom). The initial conditions are taken from Table 1 and we take into account the perturbations from the four giant planets. We show the barycentric osculating orbital elements. 
3 Secular dynamics
In this section, we describe the dynamics of the TNOs in the frame of the restricted threebody problem (Sun, Neptune, TNO), where the TNO is a test particle in an exterior orbit around Neptune. The Hamiltonian of the system can be written as (6)
where 𝒢 is the gravitational constant, M and m′ are the masses of the Sun and Neptune, respectively, r and r′ are the positions of the TNO and Neptune with respect to the Sun, and Δ = r – r′. The first two terms correspond to the Keplerian motion, ℋ_{k}, while the term in m′ corresponds to the interaction term, ℋ_{p}.
3.1 Expansion in elliptical elements
We first compute the resonant Hamiltonian by expanding ℋ_{p} in elliptical elements. We assume that the eccentricity of Neptune and its inclination to the reference plane are both zero (e′ = I′ = 0). For simplicity, we truncate the series up to the third order in the eccentricity of the TNO, e, and in the semimajor axis ratio, α = a′/a, where a and a′ are the semimajor axes of the TNO and Neptune, respectively. We further perform the variable change of the mean longitudes (λ, λ′) → (σ_{0}, λ′), and average over the fast angle λ′. We get^{3} (7)
where s = sin^{2} (I/2), and I is the inclination of the TNO. For the 3/1 MMR, we have α ≈ 0.5, and so we should have included more terms in the expansion of α. For instance, terms in α^{4} include the secular contribution cos 2ω and the resonant angle cos σ_{1}. A complete expansion of the resonant Hamiltonian in terms of Laplace coefficients can be found in Namouni & Morais (2018). However, our goal here is simply to identify the main resonant angles. We see that for I = 0, the only resonant angle is σ_{0}, while for e = 0 the only resonant angle is σ_{2}. We then search for a canonical change of variables that uses these two angles.
Fig. 3 Orbital evolution and capture in the 3/1 MMR of three TNOs during the migration of Neptune from 20 to 30 au. We show the semimajor axis, the eccentricity, the inclination, the resonant angle σ_{0} the argument of the pericentre, and the ratio between the mean motion of Neptune and the TNO (from top to bottom). 
Fig. 4 Distribution of the final eccentricity of the TNOs captured in the 3/1 MMR as a function of the initial semimajor axis. We also show the initial inclination of these objects with colour classes. 
3.2 Actionangle variables
For the TNO, we first adopt the Delaunay canonical variables: (8)
Neptune moves on an unperturbed orbit, with only one degree of freedom, given by λ′ = n′(t − t_{0}), where , and t is the time. We assume the conjugate action of λ′ to be Λ′. The Hamiltonian of the system becomes (9)
We introduce the new angles using a linear transformation (10)
which for the new actions gives (e.g. Goldstein 1950) (12)
The new set of canonical variables that uses the resonant angles is then given by (13)
The introduction of the angle v is important, because this angle does not appear in the expression of the Hamiltonian for e′ = 0 (circular orbit for Neptune). As a consequence, N is a constant of motion (a circular orbit for the planet introduces an additional symmetry in the problem). We further average over the fast angle and finally get for the resonant Hamiltonian (14)
Here, N is a parameter that depends on the initial conditions, and so the average resonant problem with e′ = 0 has only two degrees of freedom.
3.3 Semianalytical model
The interaction term, , can be obtained using series expansions, as in Eq. (7). However, as the eccentricity of the threetinos is usually high (Table 1), a large number of terms is required. An alternative and accurate way is to compute numerically over one orbital period of the TNO (e.g. Beauge 1994): (15)
The integration must be taken over λ instead of λ′ to ensure that the system attains the initial configuration. We also note that all orbital variables are fixed at the given values (S_{0}, S_{2}, σ_{0}, σ_{2};N) except λ and λ′ = 3λ − σ_{0} − 2ϖ. In order to numerically evaluate the integral Eq. (15), it is more convenient to use the eccentric anomalies E and E′, as the eccentricity of the TNO can reach high values (e.g. Pichierri et al. 2017). The position vector r = (x, y, z) is given by (16) (17) (18)
The eccentric anomaly can be related with the mean longitude through the Kepler equation: λ − ϖ = E − e sin E. The average perturbation (Eq. (15)) is then computed as (19)
Fig. 5 Level curves of the parameter N in the (a, e) plane (Eq. (20)). The vertical line marks the position of the semimajor axis a/a′ = 3^{2/3}, which corresponds to the nominal 3/1 MMR. 
3.4 Planar motion
The 3D restricted averaged resonant problem with e′ = 0 has two degrees of freedom (Eq. (14)). If we assume that the orbit of the TNO is on the same plane as Neptune, that is, I = 0, we get S_{2} = (G − H)/2 = 0 (Eq. (8)). The planar 2D problem therefore has only one degree of freedom, σ_{0}, and becomes integrable. This allows us to understand the dynamics very easily from the level curves of the resonant Hamiltonian Eq. (14), which can be obtained using the semianalytical approximation (Eq. (15)).
In the planar case, the parameter N simplifies as (Eq. (13)) (20)
that is, the semimajor axis and the eccentricity are related. In Fig. 5, we plot the level curves of constant N in the (a, e) plane for the 3/1 MMR. We observe that for e > 0.1, the eccentricity is approximately constant for all semimajor axes in the vicinity of the resonance. For simplicity, we can parametrise the different types of Hamiltonians (Eq. (14)) in terms of e instead of N. More precisely, in the following plots, we compute N from expression Eq. (20) using a given initial eccentricity and the nominal resonant semimajor axis a = 3^{2/3}a′ ≈ 62.55 au.
Another consequence of e being approximately constant for semimajor axes near the 3/1 MMR is a small amplitude for the eccentricity oscillations. In Fig. 6, we show the level curves of the Hamiltonian Eq. (14) for e ≈ 0.3 in the plane (e cos σ_{0}, e sin σ_{0}). We observe that the resonant island is indeed very narrow, which may explain why the capture probabilities in the 3/1 MMR in Sect. 2.2 were relatively small (~4%).
In order to better understand the topology of the 3/1 MMR, it is then more convenient to study the dynamics for the semimajor axis. In Fig. 7, we plot the level curves of the Hamiltonian (14) in the plane (c 0, a) for different Nvalues corresponding to different initial eccentricities, from e ≈ 0.1 to 0.99. We can observe how the equilibrium points and resonance libration areas evolve with the eccentricity. For e < 0.13, there is only a small libration zone around a stable equilibrium point centred at σ_{0} = 180° (first resonance area). There is also an unstable equilibrium point centred at σ_{0} = 0°. For 0.13 < e < 0.52, the equilibrium point at σ_{0} = 180° becomes unstable, and two resonant areas appear centred at asymmetric σ_{0} values, given by +σ_{r}. A wider libration region that encircles the three previous equilibrium points is also present, which explains why in the migration simulations, just after being captured in resonance, σ_{0} initially librates around 180° before settling around ±σ_{r} with small amplitude (Fig. 3). The width of the resonant areas increases with the eccentricity and the ±σ_{r} equilibria move away from 180°. For 0.52 < e < 0.97, the equilibrium point at σ_{0} = 0° becomes stable, giving rise to a second resonance area. As the eccentricity increases, the ±σ_{r} values approach 180° again and the width of the first resonant area decreases, while the size of the second resonant area around 0° increases. Finally, for e > 0.97, the asymmetric equilibria at ±σ_{r} merge with the symmetric equilibrium at σ_{0} = 180°, which becomes stable again.
In Fig. 7, we observe that the positions of the resonant equilibria (σ_{r}, a_{r}) change with eccentricity; they can be obtained by minimising the Hamiltonian Eq. (14), (21)
For a given eccentricity, we first compute the corresponding N value using expression (20) with a/a′ = 3^{2/3}. We then fix σ_{0} = 180° and vary the semimajor axis to find the a_{r} value that minimises the Hamiltonian. We finally fix the semimajor axis at a = a_{r} and vary σ_{0} to find the σ_{r} value that minimises the Hamiltonian. In principle, the pair (σ_{0}, a) that we find with this method could be away from the pair (σ_{r}, a_{r}) that minimises the Hamiltonian (Eq. (21)), but we find that it is already a very good approximation (see also Pichierri et al. 2017). The results are shown in Fig. 8. We observe that the a_{r} value is always close to the Keplerian value given by a/a′ = 3^{2/3}. Concerning σ_{r}, this latter moves away from 180^{°} as the eccentricity increases, with a maximum difference around e = 0.55, for σ_{r} ≈ 60°. For higher eccentricities, σ_{r} approaches 180^{°} again.
Our results are in perfect agreement with those obtained by Voyatzis et al. (2018), who used periodic orbits, and with those of Lan & Malhotra (2019), who constructed Poincaré surfaces. This shows that the semianalytic secular Hamiltonian Eq. (15) is also a good approach to represent the dynamics of the resonant threebody planar problem. Indeed, in Fig. 8 we additionally show the results of the numerical integrations of the orbital evolution of three TNOs taken from Table 1, whose evolution is also shown in Fig. 2. We observe that the angle σ_{0} of all these TNOs oscillates around one or both equilibrium values.
Fig. 6 Level curves of the averaged Hamiltonian Eq. (14) in the plane (e cos σ_{0}, e sin σ_{0}) for N = 15.43 au^{2}/yr (corresponding to e ≈ 0.3 and a/a′ = 3^{2/3}). The contours of the narrow libration area (separatrix) corresponding to the 3/1 MMR is given in red. 
4 Surface sections
In the previous section, we draw the level curves of the planar secular resonant Hamiltonian (Fig. 7). We observe that there are hyperbolic points and a separatrix surrounding the resonant equilibria, which can trigger chaotic motion in the nonaveraged problem. Therefore, following Lan & Malhotra (2019), here we compute the Poincaré surface sections to get a more realistic picture of the planar dynamics for the 3/1 MMR. The surface of section is crossed at the moment of the TNO passage of the pericentre, which corresponds to λ = ϖ and dλ/dt > 0.
For the construction of the surface sections, we integrate the Newtonian threebody equations for a massless TNO for 2 Myr. The surface sections are computed in the vicinity of the periodic solutions parametrised by N (Eq. (20)), such that they can be compared with the level curves from Sect. 3.4. For a given Nvalue, the initial conditions are determined for different values of a and σ_{0}. However, in this case, we need to start the simulations with a set of osculating coordinates. To get a better comparison with the secular problem (Eq. (14)), we additionally use a finite response filter (Carpino et al. 1987). This procedure consists in digitally filtering the output of the numerical integration by removing the fast frequencies with corresponding periods below 480 yr. This choice is sufficient to obtain similar curves to the approximation using the semianalytic model (Fig. 7).
In Fig. 9, we show the surface sections in the plane (σ_{0}, a) for the same Nvalues as in Fig. 7, corresponding to different eccentricities, from e ≈ 0.1–0.99. We observe that the surface sections provide similar results to the semianalytic secular model. However, as expected, there is one important difference: as the eccentricity increases, the trajectories near the separatrix become chaotic. As a consequence, for e ≥ 0.7, stable motion is only possible around one of the equilibrium points. For smaller values, we can still have stable trajectories around both equilibria. However, the presence of moderate chaos near the separatrix can trigger temporary excursions to different libration centres, as previously observed for some TNOs (Fig. 2) and in the migration simulations (Fig. 3).
5 Stability maps
The semianalytic secular model (Sect. 3) provides a good determination of the stationary points of the Hamiltonian, but the correct dynamics in the 3/1 MMR can only be accessed using a complete threebody model, such as the one obtained with the surface sections (Sect. 4). Above, we restrict our analysis to the plane, although most threetinos lie on inclined orbits (Table 1). Surface sections are very useful to describe problems with two degrees of freedom (such as the planar nonaveraged problem), but their use in studying problems with three or more is less straightforward. Therefore, to study the global dynamics, in this section we adopt a method based on stability maps. To this end, we use the frequency analysis method (Laskar 1990, Laskar 1993) to map the diffusion of the orbital evolution of the TNO.
We consider a grid of initial conditions and integrate the Newtonian threebody equations of motion for a time T. We then perform a frequency analysis of the mean longitude of the TNO, λ, using the software TRIP (Gastineau & Laskar 2011) over the time intervals [0, T/2] and [T/2, T], and determine the main frequency in each interval, f_{1} and f_{2}, respectively. The stability of the orbit is measured by the index (22)
which estimates the stability of the mean motion longdistance diffusion (Dumas & Laskar 1993). The larger δ, the more unstable the orbital motion of the TNO. For stable motion, we have δ ~ 0, while δ ≪ 1 if the motion is weakly perturbed, and δ ~ 1 when the motion is irregular. It is difficult to determine the precise value of δ for which the motion is stable or unstable, but a threshold of stability δ_{s} can be estimated such that most of the trajectories with δ < δ_{s} are stable (for more details see Couetdic et al. 2010).
The diffusion index depends on the considered time interval. Here, we integrate the equations of motion for T = 2 Myr, because this interval is able to capture the main characteristics of the dynamics regarding the resonant frequency, which lies within the range 10^{4} − 10^{5} yr (Table 1). With this time interval, we estimate that δ_{s} ~ 10^{−4}. The diffusion index δ is represented by a logarithmic colour scale calibrated such that blue and green correspond to stable trajectories (δ ≪ δ_{S}), while orange and red correspond to chaotic motion (δ ≫ δ_{S}). Finally, as for the surface sections (Sect. 4), we apply the finite response filter (Carpino et al. 1987), to cut off the fast frequencies with periods below 480 years, which corresponds to a numerical average over the fastfrequency terms of the Hamiltonian (Eq. (6)). However, in this case, the filter is only used to correct the values of the initial osculating orbital elements. The diffusion index (Eq. (22)) is still obtained with the frequency analysis of the full integration, because the main frequency of λ is the mean motion (a fast frequency).
Fig. 7 Level curves of the average Hamiltonian Eq. (14) for the 3/1 MMR in the plane σ_{0}, a). For each panel, N is constant and parametrised by a/a′ = 3^{2/3} and the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). The separatrix of the resonant region is given in red. 
Fig. 8 Equilibrium points for the 3/1 MMR with Neptune as a function of eccentricity. We show the equilibrium value for the semimajor axis (top) and the equilibrium value for the resonant angle (bottom). The red dashed line corresponds to the nominal value of the semimajor axis a/a′ = 3^{2/3} for comparison. 
Fig. 9 Surface sections in the plane (σ_{0}, a) for the threebody restricted planar problem close to the 3/1 MMR. For each panel, N is constant and parametrised by a/a′ = 3^{2/3} and the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). 
5.1 Planar case
We plot the stability maps in the plane (σ_{0}, a) using a grid of initial conditions with 200 × 200 points. As in the previous sections, for each grid, we fix the value of N. Then, for each grid point (σ_{0}, a), we compute the eccentricity from N using expression (20). In the restricted planar problem with e′ = 0, the dynamics does not depend on the angle v (Eq. (14)), and therefore we put v = 0. The initial position of Neptune is also irrelevant, and so we additionally set λ′ = 0, which implies that λ = (v + λ′)/3 = 0. For the longitude of the pericentre, we have ϖ = (3λ − λ′ − σ_{0})/2 = −σ_{0}/2. Using expressions (16)–(18), we finally get the initial condition for each grid point to use in a Newtonian threebody code.
In Fig. 10, we show the stability maps in the plane (σ_{0}, a) for the same Nvalues as in Figs. 7 and 9, corresponding to different eccentricities, from e ≈ 0.1 to 0.99. We observe that the stability maps provide similar results to the semianalytic secular model (Fig. 7) and to the surface sections (Fig. 9). We are no longerable to determine the exact position of the stationary points of the problem, but the different dynamical regimes of interest become quite visible, in particular the distinction between stable (blue and green) and unstable regions (red).
We observe that the circulation regions outside the 3/1 MMR are completely stable for e ≤ 0.3, but they become unstable for e ≥ 0.5. Indeed, TNOs located around 62.5 au with e ≥ 0.5 will cross the orbit of Neptune (Fig. 1), and can therefore only be stable if protected by the 3/1 MMR. Interestingly, the resonant areas that encircle a single asymmetric equilibrium point ±σ_{r} or both are also distinguishable. For e ≤ 0.5, these regions are both stable, but in the first case there is less diffusion (dark blue regions) than in the second case (light blue regions). This explains why the libration amplitude of the TNOs in this region presents some chaotic behaviour (Figs. 2 and 3).
For e = 0.1, we also note that the libration resonant region presents a higher level of diffusion than the surrounding circulation zone. This may explain why in the simulations with the migration of Neptune (Sect. 2.2) we never observed any capture in resonance for TNOs with e < 0.1. On the other hand, for e > 0.7, we observe that the second resonant region emerging around σ = 0° is also completely stable.
We conclude that the stability map method presented in this section gives similar results to the more traditionally used methods presented in Sects. 3 and 4. However, we believe the former has some practical advantages: it is easier to implement than semianalytical approximations or surface sections;itclearly distinguishes stability regions from the unstable ones and quantifies the orbital diffusion; and, more importantly, it is not limited by the number of degrees of freedom in the problem. As a result, it can be used to study the nonplanar dynamics (Sect. 5.2), or the dynamics in the presence of the remaining planets (Sect. 5.3).
Fig. 10 Stability maps in the plane (σ_{0}, a) for the threebody restricted planar problem close to the 3/1 MMR. For each panel, N is constant and parametrised by a/a′ = 3^{2/3} and the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 
5.2 Nonplanar case
We now consider that the inclination of the TNOs is no longer zero (I ≠ 0), as many of these objects are observed with moderate and high inclinations (Fig. 1). In this case, the secular problem has two degrees of freedom (Sect. 3), which prevents us from building an integrable model. It is also difficult to study the secular dynamics, because given the high eccentricity of the TNOs (Fig. 1), we need to include a large number of terms in the expansion of the interacting Hamiltonian (Eq. (7)). This is why in the planar model we used a semianalytical approach to accurately compute the secular perturbations. As explained in Sect. 5.1, stability maps are adapted to study the planar secular dynamics, and so here we adopt the same approach to study the nonplanar case.
For a better comparison with the planar case, we again plot the stability maps in the plane (σ_{0}, a) using a grid of initial conditions with 200 × 200 points. Neptune is on a circular orbit and lies in the plane of reference. In the restricted nonplanar problem with e′ = 0, the dynamics does not depend on the angle v (Eq. (14)), and we can set v = 0 and N as constants. Therefore, for each grid we still fix the value of N, which is now given by (Eq. (13)) (23)
that is, it depends also on the inclination. N is parametrised for the nominal semimajor axis a/a′ = 3^{2/3}, and using a given initial eccentricity and inclination of the TNO. Then, for each grid point (σ_{0}, a), we compute the eccentricity from N using expression (23) for the same initial value of the inclination. This choice can be contested, because for each Nvalue we have different pairs (e, I) that correspond to the same grid point (σ_{0}, a). However, as we recalculate the eccentricity in the planar case, we continue in this way in order to be able to make a comparison. Moreover, this is not a real problem, because as we use a threebody integration to produce the stability maps, the inclination is free to vary. Concerning the remaining orbital parameters, we set λ′ = 0, as in the planar case. As v = 0, we still have λ = (v + λ′)/3 = 0, and for the longitude of the pericentre, ϖ = (3λ − λ′ − σ_{0})/2 = −σ_{0}/2. However, the choice of Ω is not irrelevant, because we now have one additional degree of freedom. For simplicity, we initially chose Ω = 0°, but we also explore different initial values for this parameter. The initial condition for each grid point to use in the Newtonian threebody code is then finally obtained using expressions (16)–(18).
In Fig. 11, we show the stability maps for I = 30° and different eccentricities going from e ≈ 0.1 to 0.99. These maps are to be compared with those obtained for the planar motion (I = 0^{°}) shown in Fig. 10. Stability maps for I = 60^{°} and I = 85° are also drawn in Appendix A, but there are no significant differences with respect to the case for I = 30°. We nevertheless notice that, for all eccentricity values, the amplitude of the resonant areas decreases as we increase the inclination. Therefore, capture in the 3/1 MMR resonance becomes more difficult for higher inclination values, at least for the resonant angle σ_{0}.
For nonplanar systems, we observe that the main resonant structures of the planar case are still present, although some new interesting dynamics also emerge. For e ≈ 0.1, in the planar case there was a small single resonant island around σ_{0} = 180°. When we increase the inclination, the equilibrium point at σ_{0} = 180^{°} becomes unstable, and splits into two new resonant areas centred at asymmetric σ_{0} = ±σ_{r} values, as happens in the planar case for e > 0.13. For e = 0.7, we observe the opposite behaviour. In the planar case, there are three stability regions: two asymmetric regions centred at σ_{0} = ±σ_{r} and another symmetric region centred at σ_{0} = 0^{°}. For I = 30^{°}, the resonant regions merge in a single wider resonant area centred at σ_{0} = 0^{°}. A similar behaviour occurs for e = 0.99. In the planar case, there were two resonant areas centred at σ_{0} = 0° and 180°, and for I = 30° the two regions merge in a single resonant area around σ_{0} = 0^{°}. Interestingly, for e = 0.9 the three different resonant regions survive in the inclined case.
Another striking and very interesting feature occurs when we study the dynamics of the inclined TNOs. In the planar case, all trajectories with e ≥ 0.5 outside the 3/1 MMR are unstable, because they cross the orbit of Neptune. However, for I = 30^{°}, large stability regions are observed outside the resonance. A closer analysis of these orbits shows that the TNOs are protected by the Kozai resonance, that is, a secular resonance between the precession of the node and the precession of the pericentre, which also prevents close encounters with Neptune. For trajectories trapped inside the Kozai resonance, the angle ω can librate around 0°, 90°, 180°, or 270° (Eq.(4)). For trajectories outside this resonance but close to the separatrix, the eccentricity and the inclination can undergo large variations with maxima and minima in phase with the libration centres. The Kozai dynamics beyond Neptune has been studied in the nonresonant (Gallardo et al. 2012; Saillenfest et al. 2016) and resonant cases (Saillenfest & Lari 2017; Saillenfest et al. 2017; Lei et al. 2022), with findings confirming the existence of these different families of libration centres as a function of the eccentricity and inclination of the TNOs.
Due to our choice of initial conditions (Ω = 0°), we have ω = ϖ = −σ_{0}/2 for the stability maps, and so libration for ω around 0° and 180° occurs for stable regions that cluster around σ_{0} = 0°, while libration around 90° and 270° occurs for stable regions that cluster around σ_{0} = 180°. Figure 12 shows two different examples of trajectories with a = 63.1 au as a function of time that are stabilised by the Kozai resonance (marked with a black dot in Fig. 11), one for e = 0.5 with σ_{0} = 200° (ω = 260°), and another for e = 0.7 with σ_{0} = 300° (ω = 210°). In the first example, the TNO is outside the Kozai resonance, but its eccentricity decreases from 0.5 to nearly 0.4 whenever ω = 0° and ω = 180°. As a result, when the line of nodes between the orbital planes of the TNO and Neptune is aligned with the pericentre, the eccentricity is small enough to prevent close encounters with Neptune (see Fig. 1). In the second example, the TNO is trapped inside the Kozai resonance. Therefore, although the eccentricity oscillates around 0.7, close encounters with Neptune are expected when ω ≈ +85°; however, they never occur because we always have ω ≈ 180°.
For a better understanding of the two kinds of Kozai protection mechanism, Fig. 13 shows the evolution of the eccentricity and inclination as a function of ω for the previous two trajectories in Fig. 12. In the case with e = 0.5 (in red), we observe a clear correlation between the eccentricity minima and the alignment of the line of nodes with the pericentre (ω = 0° and 180°). In the case with e = 0.7 (in blue), we observe that ω ∈ [l45º,205º], and therefore never gets too close to ω ≈ +85°. We also observe a correlation between eccentricity and inclination. As the semimajor axis is approximately constant (Fig. 12), these variations are given by the ‘Kozai integral’, that is, (Eq. (23)).
The position of the Kozai resonance depends on the angle Ω (Eq. (4)), whose initial value was fixed at zero in Fig. 11. This choice should not impact the 3/1 MMR, as σ_{0} does not depend on Ω (Eq. (3)), but it will impact the location of the Kozai regions as a function of σ_{0}, because when Ω ≠ 0 we have ω = −σ_{0}/2 − Ω. Therefore, in Fig. 14 we plot the stability maps for I = 30° and for the eccentricities e ≈ 0.5 and 0.7, but now using^{4} Ω = 45°, 90° and 135°. As expected, the 3/1 MMR islands are still present at the same locations, but the Kozai stable regions are shifted by −2Ω.
We additionally observe that the 3/1 MMR dynamics cannot be completely dissociated from the Kozai dynamics. For e ~ 0.5, the Kozai areas seem to ‘compress’ the size of the libration width in the 3/1 MMR, at least for the outer regions, which encircle both the two asymmetric +σ_{r} equilibria. For e ≈ 0.7, the interaction dynamics is even richer. While for Ω = 0° the 3/1 MMR asymmetric islands centred at +σ_{r} were merged with the symmetric island centred at σ_{0} = 0°, for Ω = 90° we return to the configuration with three resonant centres, as in the planar case. For the intermediate Ω = 45° and Ω = 135°, the unstable regions overlap only one side of the asymmetric resonant islands, which becomes much more unstable than the side ‘protected’ by the Kozai mechanism. As a result, the asymmetric island within the Kozai region merges with the symmetric island, while the other asymmetric island remains isolated.
Fig. 11 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 30° and Ω = 0°. For each panel, N is constant and parametrised by the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). The black dots correspond to the initial conditions whose longterm evolution is shown in more detail in Fig. 12. 
Fig. 12 Orbital evolution with time for two TNOs near the Kozai resonance over 100 Myr for initial I = 30°. We show the semimajor axis, the eccentricity, the inclination, the resonant angle σ_{0}, and the argument of the pericentre (from top to bottom). The initial conditions correspond to the black dots in Fig. 11; i.e. a = 63.1 au, e = 0.5, and σ_{0} = 200° (top), and a = 63.1 au, e = 0.7, and σ_{0} = 300° (bottom). 
Fig. 13 Orbital evolution as a function of ω for two TNOs near the Kozai resonance. We show the inclination (top) and the eccentricity (bottom). The black dots correspond to the initial conditions taken from Figs. 11 and 12, i.e. a = 63.1 au and I = 30°, with e = 0.5 and σ_{0} = 200° (in red), and with e = 0.7 and σ_{0} = 300° (in blue). 
5.3 Perturbations by the giant planets
The stability maps also allow us to easily explore the dynamics in the presence of the remaining planets. As in previous studies, we take into account only the gravitational interactions with the giant planets (Jupiter, Saturn, Uranus and Neptune), because the perturbations from the rocky planets are not particularly important (e.g. Kuchner et al. 2002; Saillenfest & Lari 2017). For simplicity, we place the four planets initially in circular orbits and in the same reference plane (I′ = 0), and all the initial longitudes are taken equal to zero. We explore the dynamics of the TNO in the plane (Ω, a) and obtain the initial orbital elements in the same way as described in Sect. 5.2.
In Fig. 15 (top), we show the stability maps in the coplanar case (I = 0°) and for different eccentricities, e = 0.3, 0.5, and 0.7. These maps are to be compared with those obtained for the threebody problem shown in Fig. 10. For e = 0.3 and 0.5, the dynamics around the 3/1 MMR remains mostly unchanged; we only note an increase in the diffusion index levels, which is normal, because the mean motion of the TNO undergoes additional perturbations from the remaining planets. For e = 0.3, the nonresonant regions remain stable, but they become unstable for e = 0.5 because of close encounters with Neptune, as was found to be the case for the threebody problem (Fig. 10). However, for e = 0.7, we observe that all trajectories become unstable, including those that were trapped in the 3/1 MMR, because these TNOs now cross the orbit of Uranus. Indeed, in the frequency analysis presented in Sect. 2.1, the four TNOs presently observed in this region (see dashed line in Fig. 1) are completely unstable and are destroyed within a few million years. These objects merit further investigation in order to understand how they evolved into the present position, but this is out of the scope of the present paper.
In Fig. 15 (bottom), we show the stability maps for I = 30° and for the same previous eccentricities, e = 0.3, 0.5, and 0.7. These maps are to be compared with those obtained for the threebody problem shown in Fig. 11. For e = 0.3, the dynamics around the 3/1 MMR is also similar to the threebody problem, apart from the increase in the diffusion index levels. However, for e = 0.5, there is an important difference: the 3/1 MMR is still present, but the nonresonant stable areas now disappear as in the planar case. The stability achieved outside the 3/1 MMR in the threebody problem for e > 0.5 is possible due to the presence of the Kozai resonance (Sect. 5.2). When we include the remaining planets, they introduce more secular frequencies that modify the precession of the node and the precession of the pericentre of the TNO. As a result, the corresponding longitudes no longer evolve at the same pace and the argument of the pericentre (Eq. (4)) no longer librates or is in phase with the eccentricity oscillations, which leads to a subsequent close encounter with Neptune and ejection of the TNO. For e = 0.7, as in the planar case, we observe that all trajectories become unstable, because they cross the orbit of Uranus.
When we compare the planar dynamics (Fig. 15 top) with the dynamics at I = 30° (Fig. 15 bottom) in the presence of the four giant planets, we conclude that, as opposed to the threebody problem, there are no significant differences. The key ingredients to assure stable orbits for the threetinos family is to be captured in the 3/1 MMR with e ≲ 0.5 (higher eccentricities are scattered by Uranus). Orbits with e ≲ 0.4 can also be stable outside the 3/1 MMR, because they no longer cross the orbit of Neptune. These results are in perfect agreement with most of the currently observed TNOs in the area (Fig. 1).
Fig. 14 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 30° and Ω = 45°, 90°, and 135° from left to right. For each panel, N is constant and parametrised by the initial eccentricities e = 0.5 (top) and e = 0.7 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 
Fig. 15 Stability maps in the plane (σ_{0}, a) close to the 3/1 MMR in the presence of the four giant planets with I = 0° (top) and I = 30° (bottom). Each panel is parametrised by a different initial eccentricity, e = 0.3 (left), e = 0.5 (middle), and e = 0.7 (right). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 
6 Conclusions
A large number of TNOs are observed near a 3/1 MMR with Neptune, the majority of them with significant inclinations and eccentricities between 0.3 and 0.5 (Fig. 1). In this paper, we present a study of the dynamics of these objects. We first selected those that can be trapped in exact resonance and identify the main libration angle, σ_{0}, using frequency analysis.
We initially built a secular planar model that can be reduced to a problem with a single degree of freedom. We performed a canonical change of variables to use σ_{0} as the free coordinate. This allows us to follow the evolution of the centres of the equilibria with eccentricity and to understand the main dynamical structures in the problem. We then studied how the secular model undergoes modification in the nonaveraged planar case using surface sections, but also using stability maps. We show that stability maps are able to capture all the dynamical structures of the problem and additionally provide a clear view of the chaotic diffusion and unstable zones. Therefore, we finally use these maps to study the nonplanar dynamics and in the presence of perturbations from the remaining giant planets.
In the frame of the threebody problem, we observe that when we increase the inclination, the dynamics in the 3/1 MMR does not change significantly with respect to the planar case. There are nevertheless some resonant islands that split (for e = 0.1) or merge (for e = 0.7). The amplitude of the libration areas around c 0 also decreases as we increase the inclination. Therefore, capture in the 3/1 MMR for the resonant angle σ_{0} becomes more difficult for higher inclination values. The most striking difference concerns the secular dynamics. In the planar case, all TNOs with e ≥ 0.5 outside the 3/1 MMR are unstable, because they cross the orbit of Neptune. In the inclined case, large stability regions are observed, which correspond to trajectories trapped in a Kozai secular resonance. The presence of the Kozai resonance also renders the dynamics in 3/1 MMR more stable, leading to a merger of the libration regions of previously isolated resonant areas.
Using stability maps, we finally show that the inclusion of the remaining giant planets in the model also does not significantly change the global dynamics of the 3/1 MMR. The most important difference is that TNOs with e > 0.5 cannot survive, as these trajectories cross the orbit of Uranus. In addition, the Kozai stable areas disappear even for TNOs with e ≲ 0.5, because the perturbations from the planets modify the precession rates of the node and the pericentre of the TNO. Indeed, we confirm that the orbits of the three TNOs observed in this region (Fig. 1) are completely unstable and are destroyed in a few million years.
We focused the present study on the 3/1 MMR with Neptune, but our approach is completely general and can therefore be extended to other MMRs beyond Neptune or to the asteroid belt. Stability maps have already been used to study the dynamics of test particles in the Solar System near MMR (e.g. Robutel & Laskar 2001; Gallardo 2019; ForgácsDajka et al. 2023; Giuppone et al. 2023). However, these maps are usually plotted in the plane (a, e), which only provides the location of the resonances and their width as a function of eccentricity. We show here that a complete understanding of the full dynamics near MMR is only possible when we additionally plot these maps as a function of the resonant angle. Indeed, complex structures such as multiple resonant islands and Kozai interactions only become clear in this projection.
Acknowledgements
We thank Helena Morais for discussions. We are grateful to the referee JeanMarc Petit for his insightful comments. This work was supported by CFisUC strategic project (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FISAST/7002/2020), SmartGlow (POCI010247FEDER069733), PHOBOS (POCI010145FEDER029932), and ENGAGE SKA (POCI010145FEDER022217), funded by COMPETE 2020 and FCT, Portugal. We acknowledge the Laboratory for Advanced Computing at University of Coimbra (https://www.uc.pt/lca) for providing HPC resources to perform the stability maps with high resolution. This research has made use of data provided by the International Astronomical Union’s Minor Planet Center.
Appendix A Stability maps for inclined TNOs
Here we provide the stability maps for I = 60° and I = 85° and different eccentricities going from e ≈ 0.1 to 0.99 (we only consider the threebody problem). These maps are to be compared with those obtained for planar motion (I = 0°) shown in Fig. 10 and for I = 30° shown in Fig. 11. As we increase the inclination, we observe that the global dynamics does not change much, but the amplitude of the σ_{0} libration area decreases.
Fig. A.1 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 60° and Ω = 0°. For each panel, N is constant and parametrised by the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, 0.9, 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 
Fig. A.2 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 85° and Ω = 0°. For each panel, N is constant and parametrised by the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, 0.9, 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 
References
 Bannister, M.T., Gladman, B.J., Kavelaars, J.J., et al. 2018, ApJS, 236, 18 [CrossRef] [Google Scholar]
 Beauge, C. 1994, Celest. Mech. Dyn. Astron., 60, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Beaugé, C., Michtchenko, T.A., & FerrazMello, S. 2006, MNRAS, 365, 1160 [Google Scholar]
 Carpino, M., Milani, A., & Nobili, A.M. 1987, A&A, 181, 182 [NASA ADS] [Google Scholar]
 Couetdic, J., Laskar, J., Correia, A.C.M., Mayor, M., & Udry, S. 2010, A&A, 519, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crompvoets, B.L., Lawler, S.M., Volk, K., et al. 2022, Plan. Sci. J., 3, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Dumas, H.S., & Laskar, J. 1993, Phys. Rev. Lett., 70, 2975 [NASA ADS] [CrossRef] [Google Scholar]
 Duncan, M.J., Levison, H.F., & Budd, S.M. 1995, AJ, 110, 3073 [NASA ADS] [CrossRef] [Google Scholar]
 Elliot, J.L., Kern, S.D., Clancy, K.B., et al. 2005, AJ, 129, 1117 [NASA ADS] [CrossRef] [Google Scholar]
 Fernandez, J.A., & Ip, W.H. 1984, Icarus, 58, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Fernández, J.A., Gallardo, T., & Brunini, A. 2004, Icarus, 172, 372 [CrossRef] [Google Scholar]
 ForgácsDajka, E., Kővári, E., Kovács, T., Kiss, C., & Sándor, Z. 2023, ApJS, 266, 5 [CrossRef] [Google Scholar]
 Gallardo, T. 2019, Icarus, 317, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Gallardo, T., Hugo, G., & Pais, P. 2012, Icarus, 220, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194 [Google Scholar]
 Giuppone, C., Rodríguez, A., Alencastro, V., Roig, F., & Gallardo, T. 2023, Celest. Mech. Dyn. Astron., 135, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Gladman, B., Holman, M., Grav, T., et al. 2002, Icarus, 157, 269 [Google Scholar]
 Gladman, B., Marsden, B.G., & Vanlaerhoven, C. 2008, in The Solar System Beyond Neptune, eds. M.A. Barucci, H. Boehnhardt, D.P. Cruikshank, A. Morbidelli, & R. Dotson (Tucson: University of Arizona Press), 43 [Google Scholar]
 Goldstein, H. 1950, Classical Mechanics (Reading: AddisonWesley) [Google Scholar]
 Gomes, R., Levison, H.F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466 [CrossRef] [Google Scholar]
 Hadjidemetriou, J.D. 1993, Celest. Mech. Dyn. Astron., 56, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Hahn, J.M., & Malhotra, R. 2005, AJ, 130, 2392 [NASA ADS] [CrossRef] [Google Scholar]
 Ichtiaroglou, S., Katopodis, K., & Michalodimitrakis, M. 1989, J. Astrophys.Astron., 10, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Jewitt, D. 2005, AJ, 129, 530 [NASA ADS] [CrossRef] [Google Scholar]
 Kuchner, M.J., Brown, M.E., & Holman, M. 2002, AJ, 124, 1221 [NASA ADS] [CrossRef] [Google Scholar]
 Lan, L., & Malhotra, R. 2019, Celest. Mech. Dyn. Astron., 131, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1993, Phys. D Nonlinear Phenom., 67, 257 [Google Scholar]
 Laskar, J., & Robutel, P. 2001, Celest. Mech. Dyn. Astron., 80, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Lawler, S.M., Pike, R.E., Kaib, N., et al. 2019, AJ, 157, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Lei, H., Li, J., Huang, X., & Li, M. 2022, AJ, 164, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Malhotra, R., Lan, L., Volk, K., & Wang, X. 2018, AJ, 156, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Message, P.J. 1958, AJ, 63, 443 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Tsiganis, K., Crida, A., Levison, H.F., & Gomes, R. 2007, AJ, 134, 1790 [NASA ADS] [CrossRef] [Google Scholar]
 Namouni, F., & Morais, M.H.M. 2018, MNRAS, 474, 157 [CrossRef] [Google Scholar]
 Namouni, F., & Morais, M.H.M. 2020, MNRAS, 493, 2854 [NASA ADS] [CrossRef] [Google Scholar]
 Pichierri, G., Morbidelli, A., & Lai, D. 2017, A&A, 605, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Robutel, P., & Laskar, J. 2001, Icarus, 152, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Saillenfest, M., & Lari, G. 2017, A&A, 603, A79 [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. 2017, Celest.Mech. Dyn. Astron., 127, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, D.B. 1983a, Celest. Mech., 29, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, D.B. 1983b, Celest. Mech., 29, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H.F. 2005, Nature, 435, 459 [CrossRef] [Google Scholar]
 Voyatzis, G., & Kotoulas, T. 2005, Planet. Space Sci., 53, 1189 [NASA ADS] [CrossRef] [Google Scholar]
 Voyatzis, G., Tsiganis, K., & Antoniadou, K.I. 2018, Celest. Mech. Dyn. Astron., 130, 29 [NASA ADS] [CrossRef] [Google Scholar]
https://www.minorplanetcenter.net/ (data from March 8, 2023).
The ephemerides of the giant planets at JD2460000.5 were taken from the JPL Horizons System (https://ssd.jpl.nasa.gov/horizons).
All Tables
All Figures
Fig. 1 Barycentric orbital elements distribution of the TNOs that are close to the 3/1 MMR with Neptune. We show 76 bodies in the range [61, 64] au, with each colour corresponding to a different inclination range. The resonant TNOs are marked with a cross (Table 1), and the remaining nonresonant with a dot. The orange vertical line gives the nominal resonant semimajor axis a ≈ 62.55 au. The black solid and dashed lines give the overlap with Neptune and Uranus orbits, respectively. 

In the text 
Fig. 2 Orbital evolution of some representative TNOs near the 3/1 resonance over 100 Myr. We show the semimajor axis, the eccentricity, the inclination, the resonant angle σ_{0}, and the argument of the pericentre (from top to bottom). The initial conditions are taken from Table 1 and we take into account the perturbations from the four giant planets. We show the barycentric osculating orbital elements. 

In the text 
Fig. 3 Orbital evolution and capture in the 3/1 MMR of three TNOs during the migration of Neptune from 20 to 30 au. We show the semimajor axis, the eccentricity, the inclination, the resonant angle σ_{0} the argument of the pericentre, and the ratio between the mean motion of Neptune and the TNO (from top to bottom). 

In the text 
Fig. 4 Distribution of the final eccentricity of the TNOs captured in the 3/1 MMR as a function of the initial semimajor axis. We also show the initial inclination of these objects with colour classes. 

In the text 
Fig. 5 Level curves of the parameter N in the (a, e) plane (Eq. (20)). The vertical line marks the position of the semimajor axis a/a′ = 3^{2/3}, which corresponds to the nominal 3/1 MMR. 

In the text 
Fig. 6 Level curves of the averaged Hamiltonian Eq. (14) in the plane (e cos σ_{0}, e sin σ_{0}) for N = 15.43 au^{2}/yr (corresponding to e ≈ 0.3 and a/a′ = 3^{2/3}). The contours of the narrow libration area (separatrix) corresponding to the 3/1 MMR is given in red. 

In the text 
Fig. 7 Level curves of the average Hamiltonian Eq. (14) for the 3/1 MMR in the plane σ_{0}, a). For each panel, N is constant and parametrised by a/a′ = 3^{2/3} and the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). The separatrix of the resonant region is given in red. 

In the text 
Fig. 8 Equilibrium points for the 3/1 MMR with Neptune as a function of eccentricity. We show the equilibrium value for the semimajor axis (top) and the equilibrium value for the resonant angle (bottom). The red dashed line corresponds to the nominal value of the semimajor axis a/a′ = 3^{2/3} for comparison. 

In the text 
Fig. 9 Surface sections in the plane (σ_{0}, a) for the threebody restricted planar problem close to the 3/1 MMR. For each panel, N is constant and parametrised by a/a′ = 3^{2/3} and the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). 

In the text 
Fig. 10 Stability maps in the plane (σ_{0}, a) for the threebody restricted planar problem close to the 3/1 MMR. For each panel, N is constant and parametrised by a/a′ = 3^{2/3} and the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 

In the text 
Fig. 11 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 30° and Ω = 0°. For each panel, N is constant and parametrised by the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, e = 0.9, e = 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). The black dots correspond to the initial conditions whose longterm evolution is shown in more detail in Fig. 12. 

In the text 
Fig. 12 Orbital evolution with time for two TNOs near the Kozai resonance over 100 Myr for initial I = 30°. We show the semimajor axis, the eccentricity, the inclination, the resonant angle σ_{0}, and the argument of the pericentre (from top to bottom). The initial conditions correspond to the black dots in Fig. 11; i.e. a = 63.1 au, e = 0.5, and σ_{0} = 200° (top), and a = 63.1 au, e = 0.7, and σ_{0} = 300° (bottom). 

In the text 
Fig. 13 Orbital evolution as a function of ω for two TNOs near the Kozai resonance. We show the inclination (top) and the eccentricity (bottom). The black dots correspond to the initial conditions taken from Figs. 11 and 12, i.e. a = 63.1 au and I = 30°, with e = 0.5 and σ_{0} = 200° (in red), and with e = 0.7 and σ_{0} = 300° (in blue). 

In the text 
Fig. 14 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 30° and Ω = 45°, 90°, and 135° from left to right. For each panel, N is constant and parametrised by the initial eccentricities e = 0.5 (top) and e = 0.7 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 

In the text 
Fig. 15 Stability maps in the plane (σ_{0}, a) close to the 3/1 MMR in the presence of the four giant planets with I = 0° (top) and I = 30° (bottom). Each panel is parametrised by a different initial eccentricity, e = 0.3 (left), e = 0.5 (middle), and e = 0.7 (right). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 

In the text 
Fig. A.1 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 60° and Ω = 0°. For each panel, N is constant and parametrised by the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, 0.9, 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 

In the text 
Fig. A.2 Stability maps in the plane (σ_{0}, a) for the threebody restricted problem close to the 3/1 MMR with I = 85° and Ω = 0°. For each panel, N is constant and parametrised by the initial eccentricities (from left to right): e = 0.1, e = 0.3, e = 0.5 (top), and e = 0.7, 0.9, 0.99 (bottom). The colour scale corresponds to the relative frequency diffusion index in logarithmic scale (Eq. (22)). 

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.