# Dynamics of trans-Neptunian objects near the 3/1 mean-motion resonance with Neptune

Alan J. Alves-Carmo, Timothée Vaillant, Alexandre C. M. Correia

► To cite this version:

Alan J. Alves-Carmo, Timothée Vaillant, Alexandre C. M. Correia. Dynamics of trans-Neptunian objects near the 3/1 mean-motion resonance with Neptune. Astronomy & Astrophysics - A&A, 2023, 677, 10.1051/0004-6361/202244982. insu-04849153

HAL Id: insu-04849153

https://insu.hal.science/insu-04849153v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamics of trans-Neptunian objects near the 3/1 mean-motion resonance with Neptune

Alan J. Alves-Carmo<fnref n="1" />,<fnref n="2" />, Timothée Vaillant<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, and Alexandre C. M. Correia<fnref n="1" />,<fnref n="3" />

Received 15 September 2022 / Accepted 14 June 2023

## ABSTRACT

The complex classification of trans-Neptunian objects (TNOs) that are captured in mean-motion resonances (MMRs) and the constraint of their multiple origins are two significant open problems concerning the Solar System. The case-by-case 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 semi-analytic three-body 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 non-averaged 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 non-planar three-body 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 three-body 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 \le 0.5$.

**Key words.** celestial mechanics – Kuiper belt: general – planets and satellites: dynamical evolution and stability

## 1. Introduction

The physical mechanisms that lead to the presently observed characteristics of trans-Neptunian 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 modern-day 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 cold-gas disc and particles with random semi-major axis between 20 and 80 au provide the capture of objects in a large number of mean-motion 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 semi-analytical 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 non-planar 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 long-term evolution and capture inside the MMRs and for the formation of other TNO classes.

Hahn & Malhotra (2005) showed that a stirred-up 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

<footnote n="1">¹ CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal
e-mail: acor@uc.pt</footnote>

<footnote n="2">² Instituto de Telecomunicações, Universidade de Aveiro, 3810-193 Aveiro, Portugal</footnote>

<footnote n="3">³ IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>
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 semi-analytic 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 non-averaged planar case using surface sections (Sect. 4) and stability maps (Sect. 5). These maps are then used to investigate the non-planar 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<fnref n="1" />, we identified 82 bodies with heliocentric semi-major 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 semi-major 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 GW64, 2012 HD2, and 2013 HS150, 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<fnref n="2" />. We used the symplectic integrator SABA4 of Laskar & Robutel (2001) and an integration step size

<figure><img src="image_2.png" /><figcaption>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 non-resonant with a dot. The orange vertical line gives the nominal resonant semi-major axis $a \approx 62.55$ au. The black solid and dashed lines give the overlap with Neptune and Uranus orbits, respectively.</figcaption></figure>

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 $v_1$ and $v_2$ 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

$$3n - n' + (k - 2)v_1 - k v_2 = 0, \quad (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

$$|3n - n' + (k - 2)v_1 - k v_2|/n' < 10^{-7}. \quad (2)$$

When such a combination exists, the associated resonant angle

$$\sigma_k \equiv 3\lambda - \lambda' + (k - 2)\omega - k\Omega \quad (3)$$

is expected to librate around some equilibrium value, where $\lambda, \omega$, and $\Omega$ are the mean longitude, the longitude of the perihelion, and the longitude of the node of the TNO, respectively, and $\lambda'$ 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

<footnote n="1">¹ https://www.minorplanetcenter.net/ (data from March 8, 2023).</footnote>

<footnote n="2">² The ephemerides of the giant planets at JD2460000.5 were taken from the JPL Horizons System (https://ssd.jpl.nasa.gov/horizons).</footnote>
<table><caption>**Table 1.** TNOs in a 3/1 MMR with Neptune for at least 100 Myr.</caption><thead><tr><th>TNO</th><th>Type</th><th>Resonant<br>angle</th><th>Libration<br>centre (deg)</th><th>Libration<br>ampl. (deg)</th><th>Libration<br>period (kyr)</th><th>a<br>(au)</th><th>e<br>(deg)</th><th>I<br>(deg)</th><th>λ<br>(deg)</th><th>ω<br>(deg)</th><th>Ω<br>(deg)</th></tr></thead><tbody><tr><td>2000 SS<sub>331</sub></td><td>S</td><td>σ<sub>0</sub></td><td>78.4</td><td>40.4</td><td>19.1</td><td>62.771</td><td>0.514</td><td>2.984</td><td>272.899</td><td>272.990</td><td>146.809</td></tr><tr><td>2003 LG<sub>7</sub></td><td>C</td><td>σ<sub>0</sub></td><td>180.0</td><td>165.2</td><td>—</td><td>62.135</td><td>0.478</td><td>20.104</td><td>1.740</td><td>342.235</td><td>238.346</td></tr><tr><td>2004 VU<sub>130</sub></td><td>S</td><td>σ<sub>0</sub></td><td>180.0</td><td>158.8</td><td>44.6</td><td>62.262</td><td>0.429</td><td>8.0236</td><td>325.186</td><td>271.547</td><td>42.958</td></tr><tr><td>2011 US<sub>411</sub></td><td>S</td><td>all</td><td>70.5</td><td>31.1</td><td>18.2</td><td>62.431</td><td>0.500</td><td>22.041</td><td>17.679</td><td>9.271</td><td>24.813</td></tr><tr><td>2014 FX<sub>71</sub></td><td>C</td><td>σ<sub>0</sub></td><td>221.4</td><td>158.8</td><td>—</td><td>62.517</td><td>0.459</td><td>17.726</td><td>225.553</td><td>193.391</td><td>275.858</td></tr><tr><td>2014 SK<sub>404</sub></td><td>S</td><td>σ<sub>0</sub></td><td>76.4</td><td>29.1</td><td>23.9</td><td>62.524</td><td>0.415</td><td>26.162</td><td>254.756</td><td>243.124</td><td>128.143</td></tr><tr><td>2014 SR<sub>373</sub></td><td>S</td><td>σ<sub>0</sub></td><td>71.7</td><td>22.4</td><td>27.9</td><td>62.486</td><td>0.383</td><td>35.583</td><td>232.502</td><td>214.479</td><td>165.746</td></tr><tr><td>2015 GA<sub>55</sub></td><td>S</td><td>σ<sub>0</sub></td><td>261.2</td><td>62.5</td><td>35.6</td><td>62.625</td><td>0.387</td><td>17.958</td><td>204.417</td><td>214.371</td><td>29.045</td></tr><tr><td>2015 KY<sub>173</sub></td><td>S</td><td>σ<sub>0</sub></td><td>273.5</td><td>44.2</td><td>24.4</td><td>62.656</td><td>0.436</td><td>8.996</td><td>119.677</td><td>135.324</td><td>152.811</td></tr><tr><td>2015 RA<sub>278</sub></td><td>S</td><td>σ<sub>0</sub></td><td>86.6</td><td>43.7</td><td>23.5</td><td>62.623</td><td>0.428</td><td>11.344</td><td>140.300</td><td>112.388</td><td>202.01</td></tr><tr><td>2015 RZ<sub>277</sub></td><td>S</td><td>σ<sub>0</sub></td><td>180.0</td><td>168.5</td><td>38.1</td><td>62.247</td><td>0.445</td><td>5.827</td><td>127.120</td><td>45.786</td><td>194.182</td></tr><tr><td>2015 UH<sub>87</sub></td><td>C</td><td>σ<sub>0</sub></td><td>157.9</td><td>155.1</td><td>—</td><td>62.613</td><td>0.449</td><td>12.098</td><td>260.861</td><td>11.783</td><td>98.258</td></tr><tr><td>2015 VN<sub>166</sub></td><td>S</td><td>σ<sub>0</sub></td><td>76.4</td><td>31.4</td><td>20.2</td><td>62.716</td><td>0.481</td><td>6.030</td><td>146.592</td><td>119.125</td><td>227.357</td></tr><tr><td>2016 SO<sub>56</sub></td><td>S</td><td>σ<sub>0</sub></td><td>74.7</td><td>14.0</td><td>24.5</td><td>62.602</td><td>0.396</td><td>16.306</td><td>204.431</td><td>174.026</td><td>164.311</td></tr><tr><td>2019 GZ<sub>129</sub></td><td>S</td><td>σ<sub>0</sub></td><td>180.0</td><td>164.4</td><td>46.7</td><td>62.224</td><td>0.367</td><td>8.523</td><td>49.577</td><td>348.889</td><td>117.442</td></tr></tbody></table>

Notes. TNOs that are shown in Fig. 2 are marked with an asterisk. The ‘S’ stands for stable resonance (the libration frequency can be isolated from other perturbations), while ‘C’ indicates that the libration angle displays strong chaotic behaviour. We show the barycentric osculating orbital elements at the epoch JD2460000.5.

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 $\omega$ and $\Omega$. In a decoupled problem, we expect that $\dot{\omega} \approx v_1$ and $\dot{\Omega} \approx v_2$. However, we observe that for one object, namely 2011 US$_{411}$, both $\omega$ and $\Omega$ precess with the same fundamental frequency. As a result, the argument of the perihelion,

$$ \omega = \omega - \Omega, \qquad (4) $$

also librates. In this case, all resonant angles $\sigma_{j \neq k}$ librate, because (Eq. (3))

$$ \sigma_j = \sigma_k + (j - k) \omega. \qquad (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 long-term 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 time-span. As a result, we are able to detect the long-term 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: 2011 US$_{411}$ appears to be the most regular

object in the family, as $\sigma_0$ and $\omega$ both librate with small amplitude; 2015 KY$_{173}$ has a nearly coplanar orbit and also a small libration amplitude, but only the angle $\sigma_0$ librates; and finally, 2003 LG$_7$ 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.

## 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 semi-major axis with $\tau = 10^7$ yr. We fix the initial value of Neptune’s semi-major axis at 20 au, and vary the initial semi-major 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 $\lambda$, $\omega$, and $\Omega$ is random in the interval [0°, 360°]. We therefore integrate a total of 2178 test particles over 40 Myr. Neptune reaches the present semi-major axis, $a \approx 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
<figure><img src="image_3.png" /><figcaption>Fig. 2. Orbital evolution of some representative TNOs near the 3/1 resonance over 100 Myr. We show the semi-major axis, the eccentricity, the inclination, the resonant angle $\sigma_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.</figcaption></figure>

captures occurred for the resonant angle $\sigma_0$, but 12 also involved the angle $\omega$, and 2 occurred for the angle $\sigma_2$ alone. There were captures for all initial semi-major axis and inclination values, while for the initial eccentricity, capture was only observed when $e \ge 0.1$. For the TNOs captured in resonance, the final semi-major 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/1 MMR as a function of the initial semi-major axis. We observe that higher eccentricities are only possible for the TNOs with smaller initial semi-major 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 $\sigma_0$, the final eccentricities are mostly within $0.1 < e < 0.5$, and the inclinations are random. We also observe that, in some simulations, $\sigma_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 semi-major 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.

## 3. Secular dynamics

In this section, we describe the dynamics of the TNOs in the frame of the restricted three-body 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

$$ \mathcal{H} = \mathcal{H}_k + \mathcal{H}_p = \frac{|\dot{\mathbf{r}}|^2}{2} - \frac{G\mathcal{M}}{|\mathbf{r}|} - Gm' \left( \frac{1}{\Delta} - \frac{\mathbf{r} \cdot \mathbf{r}'}{|\mathbf{r}'|^3} \right), \quad (6) $$

where $G$ is the gravitational constant, $M$ and $m'$ are the masses of the Sun and Neptune, respectively, $\mathbf{r}$ and $\mathbf{r}'$ are the positions of the TNO and Neptune with respect to the Sun, and $\Delta = |\mathbf{r} - \mathbf{r}'|$. The first two terms correspond to the Keplerian motion, $\mathcal{H}_k$, while the term in $m'$ corresponds to the interaction term, $\mathcal{H}_p$.

### 3.1. Expansion in elliptical elements

We first compute the resonant Hamiltonian by expanding $\mathcal{H}_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 semi-major axis ratio, $\alpha = a'/a$, where $a$ and $a'$ are the semi-major axes of the TNO and Neptune, respectively. We further perform the variable change of the mean longitudes ($\lambda, \lambda'$) $\rightarrow (\sigma_0, \lambda')$, and average over the fast
<figure><img src="image_4.png" /><figcaption>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 semi-major axis, the eccentricity, the inclination, the resonant angle $\sigma_0$, the argument of the pericentre, and the ratio between the mean motion of Neptune and the TNO (from top to bottom).</figcaption></figure>

angle $\lambda'$. We get<fnref n="3" />

$$ \langle \mathcal{H}_p \rangle = -\frac{Gm'}{a} \left[ \begin{aligned} & \frac{\alpha^2}{4} \left( 1 + \frac{3}{2}e^2 \right) (1 - 6s + 6s^2) \\ & -\frac{3}{8}\alpha^{-2}e^2 (1 - 9\alpha^3) (1 - s) \cos \sigma_0 \\ & +\frac{15}{8}\alpha^3 (1 - 6e^2) s (1 - s)^2 \cos \sigma_2 \end{aligned} \right], \quad (7) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Distribution of the final eccentricity of the TNOs captured in the 3/1 MMR as a function of the initial semi-major axis. We also show the initial inclination of these objects with colour classes.</figcaption></figure>

where $s = \sin^2(I/2)$, and $I$ is the inclination of the TNO. For the 3/1 MMR, we have $\alpha \approx 0.5$, and so we should have included more terms in the expansion of $\alpha$. For instance, terms in $\alpha^4$ include the secular contribution $\cos 2\omega$ and the resonant angle $\cos \sigma_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 $\sigma_0$, while for $e = 0$ the only resonant angle is $\sigma_2$. We then search for a canonical change of variables that uses these two angles.

### 3.2. Action-angle variables

For the TNO, we first adopt the Delaunay canonical variables:

$$ \begin{aligned} L &= \sqrt{GMa}, & G &= L\sqrt{1-e^2}, & H &= G \cos I \\ l &= \lambda - \varpi, & g &= \omega, & h &= \Omega \end{aligned} \quad (8) $$

Neptune moves on an unperturbed orbit, with only one degree of freedom, given by $\lambda' = n'(t - t_0)$, where $n' = \sqrt{G(M + m')/a'^3}$, and $t$ is the time. We assume the conjugate action of $\lambda'$ to be $\Lambda'$. The Hamiltonian of the system becomes

$$ \mathcal{H} = -\frac{G^2 M^2}{2L^2} + n'\Lambda' + \mathcal{H}_p(L, G, H, l, g, h, \lambda'). \quad (9) $$

We introduce the new angles using a linear transformation

$$ \begin{bmatrix} \sigma_0 \\ \sigma_2 \\ \nu \\ \tilde{\lambda} \end{bmatrix} \equiv S \begin{bmatrix} l \\ g \\ h \\ \lambda' \end{bmatrix}, \quad (10) $$

with

$$ S = \begin{bmatrix} 3 & 1 & 1 & -1 \\ 3 & 3 & 1 & -1 \\ 3 & 3 & 3 & -1 \\ 0 & 0 & 0 & 1 \end{bmatrix}, \quad (11) $$

<footnote n="3">³ In the expression of $\cos \sigma_0$ we neglected the term in $\alpha^3$ because we already have terms in $\alpha^{-2}$.</footnote>
which for the new actions gives (e.g. Goldstein 1950)

$$ \begin{bmatrix} S_0 \\ S_2 \\ N \\ \tilde{\Lambda} \end{bmatrix} = (S^{-1})^T \begin{bmatrix} L \\ G \\ H \\ \Lambda' \end{bmatrix}. \qquad (12) $$

The new set of canonical variables that uses the resonant angles is then given by

$$ \begin{align} S_0 &= \frac{1}{2}(L - G), & \sigma_0 &= 3\lambda - \lambda' - 2\omega \nonumber \\ S_2 &= \frac{1}{2}(G - H), & \sigma_2 &= 3\lambda - \lambda' - 2\Omega \nonumber \\ N &= \frac{1}{6}(3H - L), & \nu &= 3\lambda - \lambda' \nonumber \\ \tilde{\Lambda} &= \frac{1}{3}(3\Lambda' + L), & \tilde{\lambda} &= \lambda' \tag{13} \end{align} $$

The introduction of the angle $\nu$ 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 $\tilde{\lambda}$ and finally get for the resonant Hamiltonian

$$ \langle \mathcal{H} \rangle = - \frac{\mathcal{G}^2 M^2}{18(S_0 + S_2 + N)^2} - n'(S_0 + S_2 + N) + \tilde{\mathcal{H}}_p(S_0, S_2, \sigma_0, \sigma_2; N). \qquad (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. Semi-analytical model

The interaction term, $\tilde{\mathcal{H}}_p$, 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 $\tilde{\mathcal{H}}_p$ numerically over one orbital period of the TNO (e.g. Beaugé 1994):

$$ \tilde{\mathcal{H}}_p = -Gm' \frac{1}{2\pi} \int_0^{2\pi} \left( \frac{1}{\Delta} - \frac{\mathbf{r} \cdot \mathbf{r}'}{|\mathbf{r}'|^3} \right) d\lambda. \qquad (15) $$

The integration must be taken over $\lambda$ instead of $\lambda'$ 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, \sigma_0, \sigma_2; N$) except $\lambda$ and $\lambda' = 3\lambda - \sigma_0 - 2\omega$. 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 $\mathbf{r} = (x, y, z)$ is given by

$$ x = a \left[ (\cos E - e)(\cos \omega \cos \Omega - \sin \Omega \sin \omega \cos I) - \sqrt{1 - e^2} \sin E (\sin \omega \cos \Omega + \sin \Omega \cos \omega \cos I) \right], \qquad (16) $$

$$ y = a \left[ (\cos E - e)(\cos \omega \sin \Omega + \cos \Omega \sin \omega \cos I) - \sqrt{1 - e^2} \sin E (\sin \omega \sin \Omega - \cos \omega \cos \Omega \cos I) \right], \qquad (17) $$

$$ z = a \left[ (\cos E - e) \sin \omega + \sqrt{1 - e^2} \sin E \cos \omega \right] \sin I. \qquad (18) $$

The eccentric anomaly can be related with the mean longitude through the Kepler equation: $\lambda - \omega = E - e \sin E$. The average perturbation $\tilde{\mathcal{H}}_p$ (Eq. (15)) is then computed as

$$ \tilde{\mathcal{H}}_p = -Gm' \frac{1}{2\pi} \int_0^{2\pi} \left( \frac{1}{\Delta} - \frac{\mathbf{r} \cdot \mathbf{r}'}{|\mathbf{r}'|^3} \right) (1 - e \cos E) dE. \qquad (19) $$

<figure><img src="image_6.png" /><figcaption>Fig. 5. Level curves of the parameter $N$ in the $(a, e)$ plane (Eq. (20)). The vertical line marks the position of the semi-major axis $a/a' = 3^{2/3}$, which corresponds to the nominal 3/1 MMR.</figcaption></figure>

### 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, $\sigma_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 semi-analytical approximation (Eq. (15)).

In the planar case, the parameter $N$ simplifies as (Eq. (13))

$$ N = \frac{1}{6} \sqrt{GMa} (3\sqrt{1 - e^2} - 1); \qquad (20) $$

that is, the semi-major 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 semi-major 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 semi-major axis $a = 3^{2/3}a' \approx 62.55$ au.

Another consequence of $e$ being approximately constant for semi-major 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 \approx 0.3$ in the plane $(e \cos \sigma_0, e \sin \sigma_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 semi-major axis. In Fig. 7, we plot the level curves of the Hamiltonian (14) in the plane $(\sigma_0, a)$ for different $N$-values corresponding to different initial eccentricities, from $e \approx 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 $\sigma_0 = 180^\circ$ (first resonance area). There is also an unstable equilibrium point
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Level curves of the averaged Hamiltonian Eq. (14) in the plane $(e \cos \sigma_0, e \sin \sigma_0)$ for $N = 15.43$ au$^2$/yr (corresponding to $e \approx 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.</figcaption></figure>

centred at $\sigma_0 = 0^\circ$. For $0.13 < e < 0.52$, the equilibrium point at $\sigma_0 = 180^\circ$ becomes unstable, and two resonant areas appear centred at asymmetric $\sigma_0$ values, given by $\pm\sigma_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, $\sigma_0$ initially librates around $180^\circ$ before settling around $\pm\sigma_r$ with small amplitude (Fig. 3). The width of the resonant areas increases with the eccentricity and the $\pm\sigma_r$ equilibria move away from $180^\circ$. For $0.52 < e < 0.97$, the equilibrium point at $\sigma_0 = 0^\circ$ becomes stable, giving rise to a second resonance area. As the eccentricity increases, the $\pm\sigma_r$ values approach $180^\circ$ again and the width of the first resonant area decreases, while the size of the second resonant area around $0^\circ$ increases. Finally, for $e > 0.97$, the asymmetric equilibria at $\pm\sigma_r$ merge with the symmetric equilibrium at $\sigma_0 = 180^\circ$, which becomes stable again.

In Fig. 7, we observe that the positions of the resonant equilibria ($\sigma_r, a_r$) change with eccentricity; they can be obtained by minimising the Hamiltonian Eq. (14),

$$ \frac{\partial \langle \mathcal{H} \rangle}{\partial a} = 0 \quad \text{and} \quad \frac{\partial \langle \mathcal{H} \rangle}{\partial \sigma_0} = 0. \qquad (21) $$

For a given eccentricity, we first compute the corresponding $N$ value using expression (20) with $a/a' = 3^{2/3}$. We then fix $\sigma_0 = 180^\circ$ and vary the semi-major axis to find the $a_r$ value that minimises the Hamiltonian. We finally fix the semi-major axis at $a = a_r$ and vary $\sigma_0$ to find the $\sigma_r$ value that minimises the Hamiltonian. In principle, the pair $(\sigma_0, a)$ that we find with this method could be away from the pair $(\sigma_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 $\sigma_r$, this latter moves away from $180^\circ$ as the eccentricity increases, with a

maximum difference around $e = 0.55$, for $\sigma_r \approx 60^\circ$. For higher eccentricities, $\sigma_r$ approaches $180^\circ$ 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 semi-analytic secular Hamiltonian Eq. (15) is also a good approach to represent the dynamics of the resonant three-body 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 $\sigma_0$ of all these TNOs oscillates around one or both equilibrium values.

## 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 non-averaged 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 $\lambda = \bar{\omega}$ and $d\lambda/dt > 0$.

For the construction of the surface sections, we integrate the Newtonian three-body 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 $N$-value, the initial conditions are determined for different values of $a$ and $\sigma_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 semi-analytic model (Fig. 7).

In Fig. 9, we show the surface sections in the plane $(\sigma_0, a)$ for the same $N$-values as in Fig. 7, corresponding to different eccentricities, from $e \approx 0.1-0.99$. We observe that the surface sections provide similar results to the semi-analytic 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 \ge 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 semi-analytic 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 three-body 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 non-averaged problem), but their use in studying problems with three or more is less
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Level curves of the average Hamiltonian Eq. (14) for the 3/1 MMR in the plane $(\sigma_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.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Equilibrium points for the 3/1 MMR with Neptune as a function of eccentricity. We show the equilibrium value for the semi-major axis (top) and the equilibrium value for the resonant angle (bottom). The red dashed line corresponds to the nominal value of the semi-major axis $a/a' = 3^{2/3}$ for comparison.</figcaption></figure>

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, 1993) to map the diffusion of the orbital evolution of the TNO.

We consider a grid of initial conditions and integrate the Newtonian three-body equations of motion for a time $T$. We then perform a frequency analysis of the mean longitude of the TNO, $\lambda$, 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

$$ \delta = \left| 1 - \frac{f_2}{f_1} \right|, \qquad (22) $$

which estimates the stability of the mean motion long-distance diffusion (Dumas & Laskar 1993). The larger $\delta$, the more unstable the orbital motion of the TNO. For stable motion, we have $\delta \sim 0$, while $\delta \ll 1$ if the motion is weakly perturbed, and $\delta \sim 1$ when the motion is irregular. It is difficult to determine the precise value of $\delta$ for which the motion is stable or unstable, but a threshold of stability $\delta_s$ can be estimated such that most of the trajectories with $\delta < \delta_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 $\delta_s \sim 10^{-4}$. The diffusion index $\delta$ is represented by a logarithmic colour scale calibrated such that blue and green correspond to stable trajectories ($\delta \ll \delta_s$), while orange and red correspond to chaotic motion ($\delta \gg \delta_s$). Finally, as for the surface
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Surface sections in the plane $(\sigma_0, a)$ for the three-body 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).</figcaption></figure>

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 fast-frequency 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 $\lambda$ is the mean motion (a fast frequency).

## 5.1. Planar case

We plot the stability maps in the plane $(\sigma_0, a)$ using a grid of initial conditions with $200 \times 200$ points. As in the previous sections, for each grid, we fix the value of $N$. Then, for each grid point $(\sigma_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 $\lambda' = 0$, which implies that $\lambda = (v + \lambda')/3 = 0$. For the longitude of the pericentre, we have $\omega = (3\lambda - \lambda' - \sigma_0)/2 = -\sigma_0/2$. Using expressions (16)–(18), we finally get the initial condition for each grid point to use in a Newtonian three-body code.

In Fig. 10, we show the stability maps in the plane $(\sigma_0, a)$ for the same $N$-values as in Figs. 7 and 9, corresponding to different eccentricities, from $e \approx 0.1$ to 0.99. We observe that the stability maps provide similar results to the semi-analytic secular model (Fig. 7) and to the surface sections (Fig. 9). We are no longer able

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 \le 0.3$, but they become unstable for $e \ge 0.5$. Indeed, TNOs located around 62.5 au with $e \ge 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 $\pm\sigma_r$ or both are also distinguishable. For $e \le 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 $\sigma = 0^\circ$ 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 semi-analytical approximations or surface sections; it clearly distinguishes stability regions from the unstable ones and quantifies the orbital diffusion; and, more importantly, it is not limited by
<figure><img src="image_11.png" /><figcaption>Fig. 10. Stability maps in the plane $(\sigma_0, a)$ for the three-body 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)).</figcaption></figure>

the number of degrees of freedom in the problem. As a result, it can be used to study the non-planar dynamics (Sect. 5.2), or the dynamics in the presence of the remaining planets (Sect. 5.3).

## 5.2. Non-planar case

We now consider that the inclination of the TNOs is no longer zero ($I \neq 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 semi-analytical 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 non-planar case.

For a better comparison with the planar case, we again plot the stability maps in the plane $(\sigma_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 non-planar problem with $e' = 0$, the dynamics does not depend on the angle $\nu$ (Eq. (14)), and we can set $\nu = 0$ and $N$ as constants. Therefore,

for each grid we still fix the value of $N$, which is now given by (Eq. (13))

$$N = \frac{1}{6} \sqrt{\overline{Ma}} (3 \sqrt{1 - e^2} \cos I - 1); \quad (23)$$

that is, it depends also on the inclination. $N$ is parametrised for the nominal semi-major axis $a/a' = 3^{2/3}$, and using a given initial eccentricity and inclination of the TNO. Then, for each grid point $(\sigma_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 $N$-value we have different pairs $(e, I)$ that correspond to the same grid point $(\sigma_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 three-body integration to produce the stability maps, the inclination is free to vary. Concerning the remaining orbital parameters, we set $\lambda' = 0$, as in the planar case. As $\nu = 0$, we still have $\lambda = (\nu + \lambda')/3 = 0$, and for the longitude of the pericentre, $\omega = (3\lambda - \lambda' - \sigma_0)/2 = -\sigma_0/2$. However, the choice of $\Omega$ is not irrelevant, because we now have one additional degree of freedom. For simplicity, we initially chose $\Omega = 0^{\circ}$, but we also explore different initial values for this parameter. The initial condition for each grid point to use in the Newtonian three-body code is then finally obtained using expressions (16)–(18).
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Stability maps in the plane $(\sigma_0, a)$ for the three-body restricted problem close to the 3/1 MMR with $I = 30^\circ$ and $\Omega = 0^\circ$. 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 long-term evolution is shown in more detail in Fig. 12.</figcaption></figure>

In Fig. 11, we show the stability maps for $I = 30^\circ$ and different eccentricities going from $e \approx 0.1$ to $0.99$. These maps are to be compared with those obtained for the planar motion ($I = 0^\circ$) shown in Fig. 10. Stability maps for $I = 60^\circ$ and $I = 85^\circ$ are also drawn in Appendix A, but there are no significant differences with respect to the case for $I = 30^\circ$. 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 $\sigma_0$.

For non-planar systems, we observe that the main resonant structures of the planar case are still present, although some new interesting dynamics also emerge. For $e \approx 0.1$, in the planar case there was a small single resonant island around $\sigma_0 = 180^\circ$. When we increase the inclination, the equilibrium point at $\sigma_0 = 180^\circ$ becomes unstable, and splits into two new resonant areas centred at asymmetric $\sigma_0 = \pm\sigma_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 $\sigma_0 = \pm\sigma_r$ and another symmetric region centred at $\sigma_0 = 0^\circ$. For $I = 30^\circ$, the resonant regions merge in a single wider resonant area centred at $\sigma_0 = 0^\circ$. A similar behaviour occurs for $e = 0.99$. In the planar case, there were

two resonant areas centred at $\sigma_0 = 0^\circ$ and $180^\circ$, and for $I = 30^\circ$ the two regions merge in a single resonant area around $\sigma_0 = 0^\circ$. 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 \ge 0.5$ outside the 3/1 MMR are unstable, because they cross the orbit of Neptune. However, for $I = 30^\circ$, 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 $\omega$ can librate around $0^\circ, 90^\circ, 180^\circ$, or $270^\circ$ (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 non-resonant (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.
<figure><img src="image_13.png" /></figure>

**Fig. 12.** Orbital evolution with time for two TNOs near the Kozai resonance over 100 Myr for initial $I = 30°$. We show the semi-major 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).

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 $ω ∈ [145°, 205°]$, and therefore never gets too close to $ω ≈ ±85°$. We also observe a correlation between eccentricity and inclination. As the semi-major axis is approximately constant

<figure><img src="image_14.png" /></figure>

**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).

(Fig. 12), these variations are given by the ‘Kozai integral’, that is, $\sqrt{1 - e^2} \cos I \approx \text{const}$ (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 $Ω \neq 0$ we have $ω = -σ_0/2 - Ω$. Therefore, in Fig. 14 we plot the stability maps for $I = 30°$ and for the eccentricities $e \approx 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 \approx 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 \approx 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.

### 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

<footnote>$^4$ As $σ_0 = -2(ω + Ω)$, the maps are identical for $Ω$ and $Ω + 180°$.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Stability maps in the plane $(\sigma_0, a)$ for the three-body restricted problem close to the 3/1 MMR with $I = 30^\circ$ and $\Omega = 45^\circ, 90^\circ$, and $135^\circ$ 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)).</figcaption></figure>

TNO in the plane $(\sigma_0, 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^\circ$) and for different eccentricities, $e = 0.3, 0.5$, and $0.7$. These maps are to be compared with those obtained for the three-body 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 non-resonant 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 three-body 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^\circ$ 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

three-body problem shown in Fig. 11. For $e = 0.3$, the dynamics around the 3/1 MMR is also similar to the three-body 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 non-resonant stable areas now disappear as in the planar case. The stability achieved outside the 3/1 MMR in the three-body problem for $e \ge 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^\circ$ (Fig. 15 bottom) in the presence of the four giant planets, we conclude that, as opposed to the three-body 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 \lesssim 0.5$ (higher eccentricities are scattered by Uranus). Orbits with $e \lesssim 0.4$ can also be stable outside the 3/1 MMR, because they no longer cross the orbit of Neptune.
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Stability maps in the plane $(\sigma_0, a)$ close to the 3/1 MMR in the presence of the four giant planets with $I = 0^\circ$ (top) and $I = 30^\circ$ (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)).</figcaption></figure>

These results are in perfect agreement with most of the currently observed TNOs in the area (Fig. 1).

## 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, $\sigma_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 $\sigma_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 non-averaged 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 non-planar dynamics and in the presence of perturbations from the remaining giant planets.

In the frame of the three-body 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 $\sigma_0$ also decreases as we increase the inclination. Therefore, capture in the 3/1 MMR for the resonant angle $\sigma_0$ becomes more difficult for higher inclination values. The most striking difference concerns the secular dynamics. In the planar case, all TNOs with $e \ge 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 \le 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ács-Dajka 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 Jean-Marc Petit for his insightful comments. This work was supported by CFisUC strategic project (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), SmartGlow (POCI-01-0247-FEDER-069733), PHOBOS (POCI-01-0145-FEDER-029932), and ENGAGE SKA (POCI-01-0145-FEDER-022217), 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.

## References

Bannister, M. T., Gladman, B. J., Kavelaars, J. J., et al. 2018, ApJS, 236, 18

Beauge, C. 1994, Celest. Mech. Dyn. Astron., 60, 225

Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006, MNRAS, 365, 1160

Carpino, M., Milani, A., & Nobili, A. M. 1987, A&A, 181, 182

Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, A&A, 519, A10

Crompvoets, B. L., Lawler, S. M., Volk, K., et al. 2022, Plan. Sci. J., 3, 113

Dumas, H. S., & Laskar, J. 1993, Phys. Rev. Lett., 70, 2975

Duncan, M. J., Levison, H. F., & Budd, S. M. 1995, AJ, 110, 3073

Elliot, J. L., Kern, S. D., Clancy, K. B., et al. 2005, AJ, 129, 1117

Fernandez, J. A., & Ip, W. H. 1984, Icarus, 58, 109

Fernández, J. A., Gallardo, T., & Brunini, A. 2004, Icarus, 172, 372

Forgács-Dajka, E., Kővári, E., Kovács, T., Kiss, C., & Sándor, Z. 2023, ApJS, 266, 5

Gallardo, T. 2019, Icarus, 317, 121

Gallardo, T., Hugo, G., & Pais, P. 2012, Icarus, 220, 392

Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194

Giuppone, C., Rodríguez, A., Alencastro, V., Roig, F., & Gallardo, T. 2023, Celest. Mech. Dyn. Astron., 135, 3

Gladman, B., Holman, M., Grav, T., et al. 2002, Icarus, 157, 269

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

Goldstein, H. 1950, *Classical Mechanics* (Reading: Addison-Wesley)

Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466

Hadjidemetriou, J. D. 1993, Celest. Mech. Dyn. Astron., 56, 201

Hahn, J. M., & Malhotra, R. 2005, AJ, 130, 2392

Ichtiaroglou, S., Katopodis, K., & Michalodimitrakis, M. 1989, *J. Astrophys. Astron.*, 10, 367

Jewitt, D. 2005, AJ, 129, 530

Kuchner, M. J., Brown, M. E., & Holman, M. 2002, AJ, 124, 1221

Lan, L., & Malhotra, R. 2019, Celest. Mech. Dyn. Astron., 131, 39

Laskar, J. 1990, Icarus, 88, 266

Laskar, J. 1993, *Phys. D Nonlinear Phenom.*, 67, 257

Laskar, J., & Robutel, P. 2001, Celest. Mech. Dyn. Astron., 80, 39

Lawler, S. M., Pike, R. E., Kaib, N., et al. 2019, AJ, 157, 253

Lei, H., Li, J., Huang, X., & Li, M. 2022, AJ, 164, 74

Malhotra, R., Lan, L., Volk, K., & Wang, X. 2018, AJ, 156, 55

Message, P. J. 1958, AJ, 63, 443

Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F., & Gomes, R. 2007, AJ, 134, 1790

Namouni, F., & Morais, M. H. M. 2018, MNRAS, 474, 157

Namouni, F., & Morais, M. H. M. 2020, MNRAS, 493, 2854

Pichierri, G., Morbidelli, A., & Lai, D. 2017, A&A, 605, A23

Robutel, P., & Laskar, J. 2001, Icarus, 152, 4

Saillenfest, M., & Lari, G. 2017, A&A, 603, A79

Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2016, Celest. Mech. Dyn. Astron., 126, 369

Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2017, Celest. Mech. Dyn. Astron., 127, 477

Taylor, D. B. 1983a, Celest. Mech., 29, 75

Taylor, D. B. 1983b, Celest. Mech., 29, 51

Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459

Voyatzis, G., & Kotoulas, T. 2005, Planet. Space Sci., 53, 1189

Voyatzis, G., Tsiganis, K., & Antoniadou, K. I. 2018, Celest. Mech. Dyn. Astron., 130, 29
# 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 \approx 0.1$ to 0.99 (we only consider the three-body 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 $\sigma_0$ libration area decreases.

<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Stability maps in the plane $(\sigma_0, a)$ for the three-body restricted problem close to the 3/1 MMR with $I = 60°$ and $\Omega = 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)).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. A.2.** Stability maps in the plane $(\sigma_0, a)$ for the three-body restricted problem close to the 3/1 MMR with $I = 85^\circ$ and $\Omega = 0^\circ$. 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)).</figcaption></figure>