# Effect of the inclination in the passage through the 5/3 mean motion resonance between Ariel and Umbriel

Sérgio R. A. Gomes, Alexandre C. M. Correia

► To cite this version:

Sérgio R. A. Gomes, Alexandre C. M. Correia. Effect of the inclination in the passage through the 5/3 mean motion resonance between Ariel and Umbriel. Astronomy & Astrophysics - A&A, 2023, 674, 10.1051/0004-6361/202346101. insu-04849170

HAL Id: insu-04849170

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

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
# Effect of the inclination in the passage through the 5/3 mean motion resonance between Ariel and Umbriel

Sérgio R. A. Gomes<fnref n="1" />, and Alexandre C. M. Correia<fnref n="1" />,<fnref n="2" />

Received 7 February 2023 / Accepted 7 April 2023

## ABSTRACT

The orbits of the main satellites of Uranus are expected to slowly drift away owing to tides raised in the planet. As a result, the 5/3 mean motion resonance between Ariel and Umbriel was likely encountered in the past. Previous studies have shown that, in order to prevent entrapment in this resonance, the eccentricities of the satellites must be larger than ~0.01 at the epoch, which is hard to explain. On the other hand, if the satellites experience some temporary capture and then escape, the inclinations rise to high values that are not observed today. We have revisited this problem both analytically and numerically focussing on the inclination, using a secular two-satellite model with circular orbits. We show that if the inclination of Umbriel was around 0.15° at the time of the 5/3 resonance encounter, capture can be avoided in about 60% of the cases. Moreover, after the resonance crossing, the inclination of Umbriel drops to a mean value around 0.08°, which is close to the presently observed one. The final inclination of Ariel is distributed between 0.01° and 0.25° with a nearly equal probability, which includes the present mean value of 0.02°.

**Key words.** celestial mechanics – planets and satellites: dynamical evolution and stability – planets and satellites: individual: Uranus – planets and satellites: individual: Ariel – planets and satellites: individual: Umbriel

## 1. Introduction

The origin of the Uranian satellite system is not yet completely understood and it is still under debate (e.g. Pollack et al. 1991; Szulágyi et al. 2018; Ishizawa et al. 2019; Inderbitzi et al. 2020; Ida et al. 2020; Rufu & Canup 2022). Despite the formation mechanism of the main satellites, their orbits slowly drift away owing to tides raised in the planet (e.g. Peale 1988; Tittemore & Wisdom 1988, 1989, 1990; Pollack et al. 1991; Ćuk et al. 2020). Because the orbits do not all evolve at the same pace, they may have encountered several mean motion resonances (MMRs) in their way since their formation about 4.5 Gyr ago. At present, there is no orbital commensurability between these satellites, but MMRs are often invoked to explain some anomalous observations, such as resurfacing events (e.g. Dermott et al. 1988; Peale 1988; Tittemore 1990), the current relatively large eccentricities (~0.001) of all satellites (e.g. Squyres et al. 1985; Smith et al. 1986; Peale 1988), or the high inclination of Miranda (~4.3°; e.g. Tittemore & Wisdom 1989, 1990; Malhotra & Dermott 1990; Verheylewegen et al. 2013; Ćuk et al. 2020).

The latest low-order commensurability to have occurred was the 5/3 MMR between Ariel and Umbriel (e.g. Peale 1988; Ćuk et al. 2020), which is consistent with recent geologic activity observed in Ariel (e.g. Zahnle et al. 2003; Cartwright et al. 2020). The current free eccentricities still observed also support this possibility (e.g. Jacobson 2014). A detailed study on the passage through this resonance was carried out by Tittemore & Wisdom (1988) using a secular resonant two-satellite planar model with small eccentricities. It was shown that if the 5/3 MMR is approached with eccentricities of both satellites smaller than ~0.01, a long-term capture in resonance is

certain. Therefore, in order to evade or skip this resonance, at least one of the eccentricities must have been close to 0.01 at the time. However, Tittemore & Wisdom (1988) recognised that it is very difficult to justify such high initial eccentricity values. Indeed, mutual perturbations between the largest Uranian satellites cannot explain such high values alone (e.g. Dermott & Nicholson 1986; Laskar 1986). Moreover, tidal friction within the satellites circularise the orbits in a short timescale (e.g. Squyres et al. 1985), and so any large eccentricity remnant left by a prior MMR crossing is expected to be quickly eroded.

Ćuk et al. (2020) revisited the passage through the 5/3 MMR using a *N*-body numerical integrator, which includes the five main satellites, non-planar orbits, and spin evolution. They started their simulations with the current eccentricities and adopted nearly zero inclination for all satellites (< 0.1°). They confirm that low initial eccentricities translate into capture in resonance, which can nevertheless be broken after some time due to some chaotic excitation of the eccentricities. They also observed that all five moons had their inclinations excited during the resonance entrapment. The effect on the inclination is particularly significant in Miranda because it has the smallest mass; however, after leaving the resonance, the remaining four moons were all left with inclination values higher than the current ones. In particular, Ariel and Umbriel acquired inclinations around 1°, which cannot be damped to the current observed mean values of ~0.02° and 0.08°, respectively. Ćuk et al. (2020) then propose that the inclination of Umbriel can be lowered if its node is involved in a secular spin-orbit resonance with the spin of Oberon (not observed today). However, even if this mechanism works, it would fail to damp the inclinations of the remaining satellites.

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

<footnote n="2">² IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>
Cuk et al. (2020) suggest that the 5/3 MMR between Ariel and Umbriel can be responsible for the current high inclination of Miranda instead of the crossing of the 3/1 MMR between Miranda and Umbriel (Tittemore & Wisdom 1989, 1990). However, the high inclination values also acquired by Ariel and Umbriel seem difficult to conciliate with the present observed system. A more plausible scenario is that Miranda’s inclination is excited during the Miranda and Umbriel 3/1 MMR crossing, and then Ariel and Umbriel skip the 5/3 MMR without being captured to prevent any damage on the inclinations. The problem with this scenario is that for initial near zero inclination values for both Ariel and Umbriel, at least one of these satellites must have an initial eccentricity larger than 0.01 to avoid entrapment in the 5/3 MMR (Tittemore & Wisdom 1988), which is also challenging to explain.

Since the inclination appears to play an important role in the passage of the 5/3 MMR between Ariel and Umbriel, here we intend to look at this problem with a different perspective. Instead of assuming a coplanar model with low initial eccentricities for Ariel and Umbriel as in Tittemore & Wisdom (1988), we assume a circular model with low inclinations using a similar approach to Tittemore & Wisdom (1989) for the study of the 3/1 MMR between Miranda and Umbriel. In Sect. 2, we introduce a secular resonant two-satellite circular model with low inclinations using complex Cartesian coordinates and describe all dynamical features present in the two degree-of-freedom phase space for different ratios of the semi-major axes. In Sect. 3, we add tidal effects to the equations of motion using a constant time-lag model and then compute the capture probabilities in the 5/3 MMR for different initial inclination values. In Sect. 4, we perform a large number of numerical simulations to estimate the possible outcomes of the passage through the 5/3 MMR as a function of the initial inclinations. Finally, in the last section, we summarise and discuss our results.

## 2. Resonant secular dynamics

In this section we study the conservative dynamics of a three-body system involved in a second order $(p+q)/p$ MMR, hence $q=2$. For the 5/3 MMR, we additionally have $p=3$, but our model is valid for any $p$ value. We further assume that the orbits are circular (zero eccentricities) and have low inclinations.

### 2.1. Hamiltonian

We consider an oblate central body of mass $m_0$ (Uranus) surrounded by two point-mass bodies $m_1, m_2 \ll m_0$ (satellites), where the subscript 1 refers to the inner orbit (Ariel) and the subscript 2 refers to the outer orbit (Umbriel). The potential energy of the system is given by (e.g. Smart 1953)

$$U = - \sum_{k=1}^{2} \frac{Gm_0 m_k}{r_k} \left[ 1 + J_2 \left( \frac{R}{r_k} \right)^2 P_2(\hat{\mathbf{r}}_k \cdot \mathbf{s}) \right] - \frac{Gm_1 m_2}{|\mathbf{r}_2 - \mathbf{r}_1|}, \quad (1)$$

where $G$ is the gravitational constant; $J_2, R$, and $\mathbf{s}$ are the second order gravity field, the radius, and the spin unit vector of the central body, respectively; $\mathbf{r}_k$ is the position vector of $m_k$ with respect to the centre of mass of $m_0$ (planetocentric coordinates); $r_k = |\mathbf{r}_k|$ is the norm; $\hat{\mathbf{r}}_k = \mathbf{r}_k/r_k$ is the unit vector; and $P_2(x) = (3x^2 - 1)/2$ is the Legendre polynomial of degree two. We neglected terms in $(R/r_k)^3$ (quadrupolar approximation for the oblateness). The Hamiltonian of the problem, $H$, was then obtained by adding the orbital and rotational kinetic energies to Eq. (1).

### 2.1.1. Expansion in elliptical elements

We can expand the Hamiltonian in elliptical elements. To the first order in the mass ratios, $m_k/m_0$, zeroth order in the eccentricities, and second order in the inclinations, $I_k$ (with respect to the equatorial plane of the central body), we have (e.g. Murray & Dermott 1999)

$$H = H_K + H_O + H_S + H_R + H_F + \frac{\Theta^2}{2C}, \quad (2)$$

where

$$H_K = - \sum_{k=1}^{2} \frac{Gm_0 m_k}{2a_k} \quad (3)$$

is the Keplerian part and $a_k$ are the semi-major axes,

$$H_O = - \sum_{k=1}^{2} \frac{Gm_0 m_k}{2a_k} J_2 \left( \frac{R}{a_k} \right)^2 \left( 1 - 3I_k^2 \sin^2 (\lambda_k - \Omega_k) \right) \quad (4)$$

is the contribution from the oblateness of the central body, $\lambda_k$ are the mean longitudes, $\Omega_k$ are the longitudes of the nodes,

$$H_S = - \frac{Gm_1 m_2}{8a_2} \left[ 4 b_{\frac{1}{2}}^{(0)}(\alpha) - \alpha b_{\frac{3}{2}}^{(1)}(\alpha) \left( I_1^2 + I_2^2 - 2I_1 I_2 \cos(\Omega_2 - \Omega_1) \right) \right] \quad (5)$$

is the secular part, $b_s^{(j)}$ are Laplace coefficients, $\alpha = a_1/a_2$, and

$$H_R = - \frac{Gm_1 m_2}{8a_2} \left[ I_1^2 \cos \left( (p+2)\lambda_2 - p\lambda_1 - 2\Omega_1 \right) \right. \\ \left. - 2I_1 I_2 \cos \left( (p+2)\lambda_2 - p\lambda_1 - \Omega_2 - \Omega_1 \right) \right. \\ \left. + I_2^2 \cos \left( (p+2)\lambda_2 - p\lambda_1 - 2\Omega_2 \right) \right] \alpha b_{\frac{3}{2}}^{(p+1)}(\alpha) \quad (6)$$

is the contribution from the second order resonant terms ($q=2$). The $H_S$ and $H_R$ terms arise solely from the direct part of the disturbing function (last term in Eq. (1)), because in the circular approximation the indirect part does not have secular or second order resonant terms. The term in $H_F$ corresponds to the remaining terms of the disturbing function that depend on other combinations of the angles $\lambda_1, \lambda_2, \Omega_1$, and $\Omega_2$ that do not appear in the expressions of $H_S$ or $H_R$.

Finally, for the last term in the Hamiltonian (Eq. (2)), corresponding to the rotational kinetic energy, $C$ is the principal moment of inertia of the central body and

$$\Theta = C\omega, \quad (7)$$

which corresponds to the rotational angular momentum, $\omega = \dot{\theta}$ is the angular velocity, and $\theta$ is the rotation angle. In the conservative case, the rotational kinetic energy is constant and could be dropped from the Hamiltonian. Nevertheless, when we include tidal dissipation (Sect. 3), there are angular momentum exchanges between the spin and the orbits, and this term cannot be neglected.
### 2.1.2. Action-angle resonant variables

We can rewrite the Hamiltonian (2) using a set of canonical action-angle variables. For that purpose, we adopted Andoyer variables for the rotation $(\Theta, \theta)$ and Poincaré variables for the orbits $(\Lambda_k, \lambda_k; \Phi_k, -\Omega_k)$ with

$$ \Lambda_k = \beta_k \sqrt{\mu_k a_k}, \qquad (8) $$

$$ \Phi_k = \Lambda_k (1 - \cos I_k) \approx \frac{1}{2} \Lambda_k I_k^2, \qquad (9) $$

where $\beta_k = m_0 m_k / (m_0 + m_k)$ and $\mu_k = G(m_0 + m_k)$. Since we aim to study the dynamics of the system near the $(p+2)/p$ MMR, we introduce the near resonant angle

$$ \sigma = \left(1 + \frac{p}{2}\right) \lambda_2 - \frac{p}{2} \lambda_1, \qquad (10) $$

which is present in all terms of the resonant Hamiltonian (Eq. (6)). Each term corresponds to a resonant combination:

$$ \varphi_1 = \sigma - \Omega_1, \quad \varphi_2 = \sigma - \Omega_2, \quad \varphi_3 = \frac{1}{2}(\varphi_1 + \varphi_2). \qquad (11) $$

We further introduce the angles

$$ \gamma = \frac{p}{2}(\lambda_1 - \lambda_2) = \lambda_2 - \sigma \quad \text{and} \quad \vartheta = \theta - \sigma, \qquad (12) $$

such that

$$ \begin{bmatrix} \sigma \\ \gamma \\ \varphi_1 \\ \varphi_2 \\ \vartheta \end{bmatrix} \equiv S \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ -\Omega_1 \\ -\Omega_2 \\ \theta \end{bmatrix}, \qquad (13) $$

with

$$ S = \begin{bmatrix} -p/2 & 1+p/2 & 0 & 0 & 0 \\ p/2 & -p/2 & 0 & 0 & 0 \\ -p/2 & 1+p/2 & 1 & 0 & 0 \\ -p/2 & 1+p/2 & 0 & 1 & 0 \\ p/2 & -1-p/2 & 0 & 0 & 1 \end{bmatrix}, \qquad (14) $$

which gives for the conjugated actions (e.g., Goldstein 1950)

$$ \begin{bmatrix} \Sigma \\ \Gamma \\ \tilde{\Phi}_1 \\ \tilde{\Phi}_2 \\ \tilde{\Theta} \end{bmatrix} = (S^{-1})^T \begin{bmatrix} \Lambda_1 \\ \Lambda_2 \\ \Phi_1 \\ \Phi_2 \\ \Theta \end{bmatrix}. \qquad (15) $$

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

$$ \begin{aligned} \Sigma &= \Lambda_1 + \Lambda_2 - \Phi_1 - \Phi_2 + \Theta, & \sigma &= \left(1 + \frac{p}{2}\right) \lambda_2 - \frac{p}{2} \lambda_1 \\ \Gamma &= \left(1 + \frac{2}{p}\right) \Lambda_1 + \Lambda_2, & \gamma &= \lambda_2 - \sigma \\ \tilde{\Phi}_1 &= \Phi_1, & \varphi_1 &= \sigma - \Omega_1 \\ \tilde{\Phi}_2 &= \Phi_2, & \varphi_2 &= \sigma - \Omega_2 \\ \tilde{\Theta} &= \Theta, & \vartheta &= \theta - \sigma. \end{aligned} \qquad (16) $$

### 2.1.3. Conserved quantities and average

We can rewrite the Hamiltonian (2) using the resonant canonical variables (Eq. (16)). For the actions, we can replace the semi-major axes and the inclinations using relations (8) and (9)

$$ a_k = \frac{\Lambda_k^2}{\beta_k^2 \mu_k}, \qquad (17) $$

$$ I_k \approx \sqrt{\frac{2\Phi_k}{\Lambda_k}}. \qquad (18) $$

It is important to note, however, that the $\Lambda_k$ are no longer actions of the resonant variables, and they must be obtained as (Eq. (16))

$$ \Lambda_1 = \Gamma_1 - \frac{p}{2}(\Phi_1 + \Phi_2), \qquad (19) $$

$$ \Lambda_2 = \Gamma_2 + \left(1 + \frac{p}{2}\right)(\Phi_1 + \Phi_2), \qquad (20) $$

with

$$ \Gamma_1 = \frac{p}{2} \Gamma(1 - \Delta), \qquad (21) $$

$$ \Gamma_2 = -\frac{p}{2} \Gamma \left(1 - \left(1 + \frac{2}{p}\right)\Delta\right), \qquad (22) $$

and

$$ \Delta = (\Sigma - \Theta) / \Gamma. \qquad (23) $$

In the approximation of low inclinations, $\Phi_k \ll \Lambda_k$ (Eq. (18)), and so we also have $\Phi_k \ll \Gamma_k$, allowing us to write

$$ \Lambda_1^\alpha \approx \Gamma_1^\alpha \left[ 1 - \alpha \frac{p}{2} \frac{\Phi_1 + \Phi_2}{\Gamma_1} \right], \qquad (24) $$

$$ \Lambda_2^\alpha \approx \Gamma_2^\alpha \left[ 1 + \alpha \left(1 + \frac{p}{2}\right) \frac{\Phi_1 + \Phi_2}{\Gamma_2} \right], \qquad (25) $$

and

$$ I_k \approx \sqrt{\frac{2\Phi_k}{\Gamma_k}}. \qquad (26) $$

The angle $\theta$ does not appear in the expression of the Hamiltonian (Eq. (2)), and so $\vartheta$ does not appear either (Eq. (16)). As a consequence, the conjugated action, $\Theta$, is a constant of motion. According to the d'Alembert rule (conservation of the angular momentum), the remaining angles must be combined as

$$ \cos(k_1\lambda_1 + k_2\lambda_2 + k_3\Omega_1 + k_4\Omega_2), \quad k_i \in \mathbb{Z}, \qquad (27) $$

where $k_1 + k_2 + k_3 + k_4 = 0$, or, using the resonant angles (Eq. (16)),

$$ \cos\left(\left(k_1 + \frac{2}{p}k_1 + k_2\right)\gamma - k_3\varphi_1 + \left(k_1 + k_2 + k_3\right)\varphi_2\right). \qquad (28) $$

Thus, the angle $\sigma$ also does not appear in the Hamiltonian, and its conjugated action, $\Sigma$, is a constant of motion. We note that,

$$ \begin{aligned} \Sigma &= (\Lambda_1 - \Phi_1) + (\Lambda_2 - \Phi_2) + \Theta \\ &= \beta_1 \sqrt{\mu_1 a_1} \cos I_1 + \beta_2 \sqrt{\mu_2 a_2} \cos I_2 + C\omega \end{aligned} \qquad (29) $$
<table><caption>**Table 1.** Present physical and orbital properties of the Uranian system (Thomas 1988; Jacobson 2014).</caption><thead><tr><th></th><th>m (M<sub>☉</sub> × 10<sup>-10</sup>)</th><th>R (km)</th><th>⟨T<sub>rot</sub>⟩ (day)</th><th>J<sub>2</sub>(×10<sup>-3</sup>)</th><th>C/(m<sub>0</sub>R<sup>2</sup>)</th></tr></thead><tbody><tr><td>Uranus</td><td>436562.8821</td><td>25559.</td><td>0.7183</td><td>3.5107</td><td>0.2296</td></tr><tr><td>Satellite</td><td>m (M<sub>☉</sub> × 10<sup>-10</sup>)</td><td>R (km)</td><td>⟨T<sub>orb</sub>⟩ (day)</td><td>⟨a⟩ (R)</td><td>⟨e⟩ (×10<sup>-3</sup>)</td></tr><tr><td>Ariel</td><td>6.291561</td><td>578.9</td><td>2.479971</td><td>7.468180</td><td>1.22</td></tr><tr><td>Umbriel</td><td>6.412118</td><td>584.7</td><td>4.133904</td><td>10.403550</td><td>0.0167</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.0796</td></tr></tbody></table>

does indeed correspond to the projection of the total angular momentum of the system along the direction of **s**, which must be conserved.

Near the MMR, the two resonant angles $\varphi_1$ and $\varphi_2$ present a slower variation than the angle $\gamma$, that is, $\dot{\varphi}_1, \dot{\varphi}_2 \ll \dot{\gamma}$. Therefore, we can construct the resonant secular normal form of the Hamiltonian (to the first order in $m_k/m_0$) by averaging over $\gamma$:

$$ \tilde{\mathcal{H}} = \langle \mathcal{H} \rangle_{\gamma} = \frac{1}{2p\pi} \int_{0}^{2p\pi} \mathcal{H} \, d\gamma. \qquad (30) $$

As a result, $\langle \mathcal{H}_F \rangle_{\gamma} = 0$, and since $\gamma$ no longer appears in the expression of the averaged Hamiltonian, the conjugated variable, $\Gamma$ (Eq. (16)), also becomes a constant of motion. We thus reduced a problem with five degrees of freedom initially to a problem with only two degrees of freedom, $(\Phi_1, \varphi_1)$ and $(\Phi_2, \varphi_2)$, and three parameters, $\Sigma$, $\Gamma$, and $\Theta$. The auxiliary quantities $\Gamma_1$ (Eq. (21)) and $\Gamma_2$ (Eq. (22)) are also constant.

The resonant secular Hamiltonian then finally reads as follows:

$$ \begin{aligned} \tilde{\mathcal{H}} ={}& (\mathcal{K}_a + \mathcal{S}_a)(\Phi_1 + \Phi_2) + \mathcal{K}_b(\Phi_1 + \Phi_2)^2 \\ & + (\mathcal{O}_a + \mathcal{S}_b)\Phi_1 + (\mathcal{O}_b + \mathcal{S}_c)\Phi_2 \\ & + \mathcal{S}_d \sqrt{\Phi_1} \sqrt{\Phi_2} \cos(\varphi_1 - \varphi_2) \\ & + \mathcal{R}_a \Phi_1 \cos(2\varphi_1) + \mathcal{R}_b \Phi_2 \cos(2\varphi_2) \\ & + \mathcal{R}_c \sqrt{\Phi_1} \sqrt{\Phi_2} \cos(\varphi_1 + \varphi_2), \end{aligned} \qquad (31) $$

where $\mathcal{K}$ stands for the Keplerian coefficients (Eq. (3)), $\mathcal{O}$ for the oblateness coefficients (Eq. (4)), $\mathcal{S}$ for secular the coefficients (Eq. (5)), and $\mathcal{R}$ for the resonant coefficients (Eq. (6)). The Keplerian part needs to be expanded to the second order in $\Phi_k$ (fourth order in the inclinations), because $\mathcal{K}_b$ is much larger than the remaining coefficients. The explicit expression of all these coefficients is given in Appendix A.

### 2.1.4. Complex rectangular coordinates

The equations of motion expressed in the variables $(\Phi_k, \varphi_k)$ may experience some singularities when $\Phi_k = 0$ (because of the terms in $\mathcal{S}_d$ and $\mathcal{R}_c$). To remove this problem, we can perform a second canonical change of variables to rectangular coordinates $(\Phi_k, \varphi_k) \rightarrow (y_k, i\bar{y}_k)$, where

$$ y_k = \sqrt{\Phi_k} e^{i\varphi_k}, \qquad (32) $$

and $\bar{y}_k$ is the complex conjugate of $y_k$. From Eq. (26), we have

$$ y_k \approx I_k \sqrt{\frac{\Gamma_k}{2}} e^{i\varphi_k}, \qquad (33) $$

and so these variables are proportional to the inclinations. The resonant secular Hamiltonian (Eq. (31)) now reads as follows:

$$ \begin{aligned} \tilde{\mathcal{H}} ={}& (\mathcal{K}_a + \mathcal{S}_a)(y_1\bar{y}_1 + y_2\bar{y}_2) + \mathcal{K}_b(y_1\bar{y}_1 + y_2\bar{y}_2)^2 \\ & + (\mathcal{O}_a + \mathcal{S}_b)y_1\bar{y}_1 + (\mathcal{O}_b + \mathcal{S}_c)y_2\bar{y}_2 + \frac{\mathcal{S}_d}{2}(y_1\bar{y}_2 + \bar{y}_1y_2) \\ & + \frac{\mathcal{R}_a}{2}(y_1^2 + \bar{y}_1^2) + \frac{\mathcal{R}_b}{2}(y_2^2 + \bar{y}_2^2) + \frac{\mathcal{R}_c}{2}(y_1y_2 + \bar{y}_1\bar{y}_2). \end{aligned} \qquad (34) $$

The conservative equations of motion are simply obtained as

$$ \frac{dy_k}{dt} = i \frac{\partial \tilde{\mathcal{H}}}{\partial \bar{y}_k}, \qquad (35) $$

yielding

$$ \dot{y}_1 = i \left[ (\mathcal{K}_a + \mathcal{S}_a) y_1 + 2\mathcal{K}_b (y_1\bar{y}_1 + y_2\bar{y}_2) y_1 + (\mathcal{O}_a + \mathcal{S}_b) y_1 + \frac{\mathcal{S}_d}{2} y_2 + \mathcal{R}_a \bar{y}_1 + \frac{\mathcal{R}_c}{2} \bar{y}_2 \right] \qquad (36) $$

and

$$ \dot{y}_2 = i \left[ (\mathcal{K}_a + \mathcal{S}_a) y_2 + 2\mathcal{K}_b (y_1\bar{y}_1 + y_2\bar{y}_2) y_2 + (\mathcal{O}_b + \mathcal{S}_c) y_2 + \frac{\mathcal{S}_d}{2} y_1 + \mathcal{R}_a \bar{y}_2 + \frac{\mathcal{R}_c}{2} \bar{y}_1 \right]. \qquad (37) $$

### 2.2. Dynamical evolution

The values of $m_k$, $J_2$, and $R$ are relatively well determined for the Uranian system (Table 1). Therefore, to compute the coefficients $\mathcal{K}$, $\mathcal{O}$, $\mathcal{S}$, and $\mathcal{R}$ appearing in the Hamiltonian (34), we only need to know the values of the parameters $\Gamma_1$ and $\Gamma_2$ (see Appendix A), which in turn depend on the parameters $\Gamma$ and $\Delta$ (Eqs. (21) and (22)).

If we neglect the $\mathcal{O}$ coefficients (oblateness coefficients), it is possible to eliminate the dependence in $\Gamma$, because we get $\tilde{\mathcal{H}} \propto \Gamma^{-2}$ (see Appendix A and Delisle et al. 2012, Sect. 2.1.4). Although here we cannot use this simplification, the conservative dynamics is still not very sensitive to the $\Gamma$ parameter (e.g. Tittemore & Wisdom 1988, 1989), and so we fixed it at the reference value<fnref n="1" />

$$ \Gamma = 2.6684 \times 10^{-12} M_{\odot} \text{au}^2 \text{yr}^{-1}. \qquad (38) $$

<footnote n="1">¹ The value of $\Gamma$ was obtained by reversing the orbital tidal evolution of the system (Eq. (73)) until the nominal resonance was achieved (Eq. (77)) and starting with the present semi-major axes (see Sect. 3.3.2).</footnote>
The dynamics of the 5/3 MMR essentially depends on the
$\Delta$ parameter (Eq. (23)), which measures the proximity to the
resonance. Following Delisle et al. (2012), we write

$$
\delta = \frac{\Delta}{\Delta_r} - 1, \tag{39}
$$

where $\Delta_r$ is the value of $\Delta$ at the planar ($I_k = 0$) nominal
resonance, that is (Eq. (23) and (29)),

$$
\Delta_r = (\Lambda_{1,r} + \Lambda_{2,r}) / \Gamma_r, \quad (40)
$$

when $n_1/n_2 = 5/3$, where $n_k$ is the mean motion of the satellite
with mass $m_k$. At the nominal resonance, using Kepler's third
law, we have the following relation for the $\Lambda_{k,r}$:

$$
\Lambda_{2,r} = \epsilon \left(\frac{5}{3}\right)^{1/3} \Lambda_{1,r}, \quad \text{with} \quad \epsilon = \frac{\beta_2 \mu_2^{2/3}}{\beta_1 \mu_1^{2/3}} \approx \frac{m_2}{m_1}. \qquad (41)
$$

From the expression of $\Gamma_r$ (Eq. (16)), we additionally have

$$
\Lambda_{2,r} = \Gamma_r - \frac{5}{3} \Lambda_{1,r}, \quad (42)
$$

which can be combined with expression (41) to give

$$
\Lambda_{1,r} = \left(\frac{5}{3} + \epsilon\left(\frac{5}{3}\right)^{1/3}\right)^{-1} \Gamma_r \quad (43)
$$

$$
\Lambda_{2,r} = \left(1 + \epsilon^{-1} \left(\frac{5}{3}\right)^{2/3}\right)^{-1} \Gamma_r, \quad (44)
$$

and finally (Eq. (40))

$$
\Delta_r = \left(1 + \epsilon\left(\frac{5}{3}\right)^{1/3}\right)\left(\frac{5}{3} + \epsilon\left(\frac{5}{3}\right)^{1/3}\right)^{-1}. \quad (45)
$$

2.2.1. Equilibrium points

The equilibrium points correspond to stationary solutions of the
Hamiltonian. They can be found by solving

$$
\frac{\partial \bar{H}}{\partial y_1} = 0 \quad \text{and} \quad \frac{\partial \bar{H}}{\partial y_2} = 0, \qquad (46)
$$

that is, finding the roots of Eqs. (36) and (37). Splitting these
equations into their real and imaginary parts, $y_k = y_{k,r} + i y_{k,i}$,
we get

$$
y_{2,r} = -2 \frac{\mathcal{K}_a + 2\mathcal{K}_b (\Phi_1 + \Phi_2) + O_a + S_a + S_b + \mathcal{R}_a}{S_d + \mathcal{R}_c} y_{1,r}, \quad (47)
$$

$$
y_{1,r} = -2 \frac{\mathcal{K}_a + 2\mathcal{K}_b (\Phi_1 + \Phi_2) + O_b + S_a + S_c + \mathcal{R}_b}{S_d + \mathcal{R}_c} y_{2,r}, \quad (48)
$$

$$
y_{2,i} = -2 \frac{\mathcal{K}_a + 2\mathcal{K}_b (\Phi_1 + \Phi_2) + O_a + S_a + S_b - \mathcal{R}_a}{S_d - \mathcal{R}_c} y_{1,i}, \quad (49)
$$

$$
y_{1,i} = -2 \frac{\mathcal{K}_a + 2\mathcal{K}_b (\Phi_1 + \Phi_2) + O_b + S_a + S_c - \mathcal{R}_b}{S_d - \mathcal{R}_c} y_{2,i}. \quad (50)
$$

Replacing expressions (48) and (50) into expressions (47)
and (49), respectively, yields that equilibria arise for

$$
y_{1,r} = y_{2,r} = 0 \quad \text{or} \quad y_{1,i} = y_{2,i} = 0. \tag{51}
$$

A more deep analysis shows that stable equilibria can only occur
when the real roots are null. We then focussed on these roots
to determine the exact position of the stable equilibria. Since
$y_{1,r} = y_{2,r} = 0$, we have (Eq. (32))

$$
y_{1,i} = \pm \sqrt{\Phi_1} \quad \text{and} \quad y_{2,i} = \pm \sqrt{\Phi_2}, \qquad (52)
$$

which we replaced in expressions (49) and (50) to determine the
possible values of $\Phi_1$ and $\Phi_2$,

$$
\Phi_1 = \Phi_2 = 0 \tag{53}
$$

or

$$
\Phi_1 = \frac{(\mathcal{R}_c - \mathcal{S}_d)\varepsilon_{\pm} - 2(\mathcal{K}_a + \mathcal{O}_a + \mathcal{S}_a + \mathcal{S}_b - \mathcal{R}_a)}{4\mathcal{K}_b(1 + \varepsilon_{\pm}^2)}, \quad (54)
$$

$$
\Phi_2 = \epsilon_{\pm}^2 \Phi_1, \quad (55)
$$

with

$$
\epsilon_{\pm} = \frac{O_a - O_b - \mathcal{R}_a + \mathcal{R}_b + S_b - S_c}{\mathcal{R}_c - S_d} \pm \frac{\sqrt{(O_a - O_b - \mathcal{R}_a + \mathcal{R}_b + S_b - S_c)^2 + (\mathcal{R}_c - S_d)^2}}{\mathcal{R}_c - S_d}. \quad (56)
$$

The equilibrium point at $\Phi_1 = \Phi_2 = 0$ is always present, although it can be stable or unstable. The remaining equilibria only exist for some $\delta$ values. In Fig. 1, we show the evolution of the equilibrium points as a function of $\delta$. We rescaled the $y_{k,i}$ by $\sqrt{\Gamma_k/2}$, such that we can translate the different equilibria in terms of inclinations (Eq. (33)).

For positive $\delta$ values far from zero, there is only one equi-
librium point at $I_k = 0$, which is stable (in blue). For $\delta =
2.068 \times 10^{-6}$, there is a first bifurcation in the equilibria: two
new stable equilibrium points appear at non-zero inclination
(in green), while the point at $I_k = 0$ becomes unstable. For
$\delta = -2.540 \times 10^{-7}$, which is very close to the resonance nominal
value $\delta = 0$ (Eq. (40)), a second bifurcation arises: two additional
unstable equilibrium points appear at a non-zero inclination (in
red), while the point at $I_k = 0$ becomes stable again.

2.2.2. Energy levels

For a better understanding of the dynamics, we can look at the
energy levels of the resonant secular Hamiltonian (34) for differ-
ent values of $\delta$ (Eq. (45)). Since our problem has two degrees
of freedom, and thus four dimensions, we need to plot these
levels on section planes. We chose the plane ($y_{1,i}$, $y_{2,i}$) with
$y_{1,r} = y_{2,r} = 0$, in order for stable equilibria to become visi-
ble (Eq. (52)). In Fig. 2, we show the energy levels for three
values of $\delta$, which are representative of the three different equi-
librium possibilities appearing in Fig. 1. Again, we rescaled $y_{k,i}$
by $\sqrt{\Gamma_k/2}$, and so we actually show the energy levels in the plane
($I_1 \sin \varphi_1$, $I_2 \sin \varphi_2$) with $\cos \varphi_1 = \cos \varphi_2 = 0$.

For $\delta = 5 \times 10^{-6} > 0$ (Fig. 2a), there is a single equilibrium point ($y_{1,i} = 0, y_{2,i} = 0$) at the centre (in blue). It corresponds to a fixed point of the Hamiltonian (34), which is surrounded by a circulating region. Therefore, in this case (and for higher $\delta$ values), all trajectories are outside the 5/3 MMR.

For $\delta = 0$ (Fig. 2b), the system is at the nominal reso-
nance (Eq. (39)). In this case, the equilibrium point at the centre
($y_{1,i} = 0, y_{2,i} = 0$) is still present (in red), but it is now unstable.
Indeed, there is a separatrix in a tilted ‘8’ shape emerging from
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Evolution of the equilibrium points as a function of $\delta$. The green lines represent stable points inside the resonance (in a libration region), the red lines represent hyperbolic points (unstable), and the blue lines represent stable fixed points (in a circulation region).</figcaption></figure>

this point that surrounds two additional stable equilibrium points (in green). Trajectories inside the separatrix that encircle the stable points are in libration and correspond to orbits inside the 5/3 MMR. Trajectories outside the separatrix are in circulation.

Finally, for $\delta = -2 \times 10^{-6} < 0$ (Fig. 2c), several equilibria exist. There are two hyperbolic points (in red) from which a separatrix with two ‘banana’ shapes emerges. This separatrix delimits the phase space in libration and circulation regions. There are two stable points (in green), one inside each banana island. Trajectories that move around these points are in libration and correspond to orbits inside the 5/3 MMR. The point $(y_{1,i} = 0, y_{2,i} = 0)$ at the centre (in blue) is again stable and inside a small circulation region. Trajectories outside the separatrix are also in circulation. This kind of phase space persists for smaller $\delta$ values, but the central circulation region becomes larger, while the resonant islands become thinner.

### 2.2.3. Surface sections

The energy levels from previous the section allowed us to identify the different regions of the phase space, but a priori they do not correspond to trajectories followed by the system. Indeed, since our problem has four dimensions, they show the trajectories when they cross the section plane with $y_{1,r} = y_{2,r} = 0$, which only remain constant for the equilibrium points. An alternative projection consists of fixing only $y_{1,r} = 0$ (or $y_{2,r} = 0$) together with a constant energy, that is, to draw Poincaré surface sections. This projection is less restrictive, and therefore allows us to distinguish between periodic (fixed points), quasi-periodic (closed curves), and chaotic trajectories. For that purpose, we used the modified Hénon method (Henon & Heiles 1964; Palaniyandi 2009).

In Fig. 2, we observe that the more diverse dynamics occurs for $\delta = -2 \times 10^{-6}$. We thus adopted this value to draw the surface

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Energy level curves in the plane $(I_1 \sin \varphi_1, I_2 \sin \varphi_2)$ with $\cos \varphi_1 = \cos \varphi_2 = 0$, for $\delta = 5 \times 10^{-6}$ (top), $\delta = 0$ (middle), and $\delta = -2 \times 10^{-6}$ (bottom). Stable equilibria are coloured in green (resonance) and blue, while unstable equilibria are coloured in red, as well as the level curves that correspond to the separatrix.</figcaption></figure>

sections. Since the Hamiltonian (34) is a four-degree function of $y_k$, the intersection of the constant energy manifold by a plane may have up to four roots (families). Each family corresponds
to a different dynamical behaviour, and so we must plot one of them at a time. However, the families are symmetric and actually we only need to show two of them. We chose to represent the families with the positive roots (that we dub 1 and 2).

In Fig. 3, we show a set of surface sections for the motion of Ariel in the plane ($y_{1,i}$, $y_{1,r}$) with $y_{2,r} = 0$. We rescaled $y_k$ again by $\sqrt{\Gamma_k/2}$, and so we actually show the surface sections in the plane ($I_1 \sin \varphi_1$, $I_1 \cos \varphi_1$) with $\cos \varphi_2 = 0$. Each panel corresponds to a different energy value, which coincide with the energy levels that are shown in Fig. 2c (we adopted the same colour code). The lowest energies occur in the circulation regions, $\tilde{\mathcal{H}} < \mathcal{H}_0$, the separatrix corresponds to a transition level, $\tilde{\mathcal{H}} = \mathcal{H}_0$, while the largest energies occur in the libration region, $\tilde{\mathcal{H}} > \mathcal{H}_0$. The inner circulation region is delimited by the levels $0 < \tilde{\mathcal{H}} < \mathcal{H}_0$, where $\tilde{\mathcal{H}} = 0$ corresponds to the energy of the equilibrium point with $y_1 = y_2 = 0$ (blue point in Fig. 2c). For this energy range, there are four families, while for the remaining energies only two families exist.

For $\tilde{\mathcal{H}} \ll 0$ (Fig. 3a), only family 1 exists, and we observe that the system is always quasi-periodic, corresponding to trajectories in the outer circulation region. As the energy increases, two islands appear, corresponding to trajectories that are in the libration region (in resonance). Initially, the motion in these new regions is also quasi-periodic. However, as the energy approaches the threshold $\tilde{\mathcal{H}} = 0$ (Figs. 3b,c), some chaotic regions appear in the transition between the circulation and libration regions. For $0 < \tilde{\mathcal{H}} < \mathcal{H}_0$ (family 1), the chaotic regions increase, while the resonant islands shrink (Figs. 3d,e), until they completely disappear for $\tilde{\mathcal{H}} = \mathcal{H}_0$ (Fig. 3f). For this specific energy range, we also needed to plot family 2. Close to $\tilde{\mathcal{H}} = 0$, we observe quasi-periodic motion in the inner circulation region, but, as we approach $\tilde{\mathcal{H}} = \mathcal{H}_0$, this area is also completely replaced by a chaotic region (Figs. 3j–l). Finally, for $\tilde{\mathcal{H}} > \mathcal{H}_0$, we observe that the chaotic region progressively vanishes and it is replaced by quasi-periodic motion in the libration region (Figs. 3g–i). In this energy range, we only have family 1 and trajectories in the outer circulation region also do not exist. Moreover, there is also a forbidden region at the centre of each panel that grows with the energy value while the libration areas shrink.

From the analysis of the surface sections, we conclude that the dynamics of the 5/3 MMR between Ariel and Umbriel is very rich and depends on the energy of the system. In fact, the energy depends on the value of the inclinations (Eq. (34)), given by the variables $y_1$ and $y_2$ (Eq. (33)). Therefore, the value of the inclinations of Ariel and Umbriel when the system encounters the resonance can trigger completely different behaviours. For $\tilde{\mathcal{H}} < 0$, the motion is quasi-periodic, either in circulation or libration. Near the separatrix, $\tilde{\mathcal{H}} \sim \mathcal{H}_0$, the motion is mainly chaotic. Finally, for $\tilde{\mathcal{H}} \gg \mathcal{H}_0$, the motion is again quasi-periodic, but only possible in libration with a small amplitude around the high inclination stable equilibrium points (Fig. 1).

## 3. Tidal evolution

The resonant dynamics presented in the previous section is conservative and thus the average semi-major axes remain constant. However, the orbits of the Uranian satellites are expected to evolve because of tidal interactions. The tidal contributions to the orbital and spin evolution can be obtained by considering an additional tidal potential (e.g. Darwin 1880; Kaula 1964).

### 3.1. Tidal potential energy

Tides arise from differential and inelastic deformations of an extended body (e.g. Uranus) owing to the gravitational effect of a perturber (e.g. Ariel or Umbriel). The resulting distortion gives rise to a tidal bulge, which modifies the gravitational potential of the extended body. The dissipation of the mechanical energy of tides inside the body introduces a time delay, $\tau$, between the initial perturbation and the maximal deformation. As the perturber interacts with the additional potential field, the amount of tidal potential energy is given by

$$U_k = -k_2 \frac{G m_k^2 R^5}{R} \left(\frac{R}{r_k}\right)^3 \left(\frac{R}{r'_k}\right)^3 P_2(\hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}'_k), \quad (57)$$

where $k_2$ is the elastic second Love number for potential of the body, $\mathbf{r}_k = \mathbf{r}_k(t)$ is the position of the pertuber at a time $t$, and $\mathbf{r}'_k = \mathbf{r}_k(t - \tau)$ is its position when it exerts the perturbation. In this work we solely consider the deformations raised on Uranus by its satellites, since we assumed the satellites as point masses (Sect. 2). This choice is fully justified because for circular orbits and synchronous satellites the tides raised by the planet on its satellites can be neglected (e.g. Correia 2009).

Although tidal effects do not preserve the mechanical energy, it is possible to extend the Hamiltonian formalism from Sect. 2 by considering the primed quantities, $\mathbf{r}'_k$, as parameters (Mignard 1979). The tidal Hamiltonian then reads (Eqs. (2) and (57)) as follows:

$$\mathcal{H}_t = \mathcal{H} + U_1 + U_2. \quad (58)$$

As in Sect. 2.1.1, we first expanded $U_k$ in elliptical elements. To the first order in the mass ratios, zeroth order in the eccentricities, and second order in the inclinations, we have

$$\begin{aligned} U_k = & -k_2 \frac{G m_k^2 R^5}{4 a_k^3 a_k'^3} \left[ \right. \\ & \left. \begin{aligned} & 1 + 3 \cos(2\theta - 2\theta' - 2\lambda_k + 2\lambda_k') \\ & - \frac{3}{2} I_k^2 \left( 1 - \cos(2\lambda_k - 2\Omega_k) + \cos(2\theta - 2\theta' - 2\lambda_k + 2\lambda_k') \right. \\ & \left. \quad - \cos(2\theta - 2\theta' + 2\lambda_k' - 2\Omega_k) \right) \\ & - \frac{3}{2} I_k'^2 \left( 1 - \cos(2\lambda_k' - 2\Omega_k') + \cos(2\theta - 2\theta' - 2\lambda_k + 2\lambda_k') \right. \\ & \left. \quad - \cos(2\theta - 2\theta' - 2\lambda_k + 2\Omega_k') \right) \\ & + 3 I_k I_k' \left( \cos(\theta - \theta' - 2\lambda_k + 2\lambda_k' + \Omega_k - \Omega_k') \right. \\ & \left. \quad + \cos(\theta - \theta' - \Omega_k + \Omega_k') \right. \\ & \left. \quad - \cos(\theta - \theta' - 2\lambda_k + \Omega_k + \Omega_k') \right. \\ & \left. \quad - \cos(\theta - \theta' + 2\lambda_k' - \Omega_k - \Omega_k') \end{aligned} \right]. \end{aligned} \quad (59)$$

We first performed the canonical change of variables that uses the resonant angles (Eq. (16)), and then changed to the complex
<figure><img src="image_4.png" /><figcaption>Family 2</figcaption></figure>

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

**Fig. 3.** Poincaré surfaces of section for Ariel in the plane ($I_1 \sin \varphi_1, I_1 \cos \varphi_1$) with $\cos \varphi_2 = 0$ and $\delta = -2 \times 10^{-6}$. Each panel was obtained with a different energy value, corresponding to the energy levels shown in Fig. 2c (we adopted the same colour code), and $H_0 = 1.06 \times 10^{-19} M_{\odot} \text{au}^2 \text{yr}^{-2}$.
Cartesian coordinates (Eq. (32)), to get

$$
\begin{align}
U_k = & -k_2 \frac{\mathcal{G} m_k^2 R^5 \beta_k^{12} \mu_k^6}{4 \Gamma_k^7 \Gamma_k^{7'}} \left[ \begin{aligned}[t]
& \Gamma_k \Gamma_k' - 3 \Gamma_k' \bar{y}_k y_k - 3 \Gamma_k \bar{y}_k' y_k' \\
& - 6 p_k \left( \Gamma_k' y_1 \bar{y}_1 + \Gamma_k y_1' \bar{y}_1' + \Gamma_k' y_2 \bar{y}_2 + \Gamma_k y_2' \bar{y}_2' \right) \\
& + 3 \left( \Gamma_k \Gamma_k' - \Gamma_k' \bar{y}_k y_k - \Gamma_k \bar{y}_k' y_k' \right) \\
& - 6 p_k \left( \Gamma_k' y_1 \bar{y}_1 + \Gamma_k y_1' \bar{y}_1' + \Gamma_k' y_2 \bar{y}_2 + \Gamma_k y_2' \bar{y}_2' \right) \\
& \times \cos(2(\theta - \vartheta' - q_k(\gamma - \gamma'))) \\
& + \frac{3}{2} \Gamma_k \left( \bar{y}_k'^2 + y_k'^2 \right) \cos(2q_k \gamma') \\
& + \frac{3}{2} \Gamma_k' \left( \bar{y}_k^2 + y_k^2 \right) \cos(2q_k \gamma) \\
& + 3 \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k y_k' + \bar{y}_k' y_k \right) \cos(\theta - \vartheta') \\
& - 3 \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k \bar{y}_k' + y_k y_k' \right) \cos(\theta - \vartheta' + 2q_k \gamma') \\
& - 3 \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k \bar{y}_k' + y_k y_k' \right) \cos(\theta - \vartheta' - 2q_k \gamma) \\
& + 3 \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k y_k' + \bar{y}_k' y_k \right) \cos(\theta - \vartheta' - 2q_k(\gamma - \gamma')) \end{aligned} \right] \tag{60} \\
& + \frac{3}{2} \Gamma_k' \left( \bar{y}_k^2 + y_k^2 \right) \cos(2(\theta - \vartheta' + q_k \gamma')) \\
& + \frac{3}{2} \Gamma_k \left( \bar{y}_k'^2 + y_k'^2 \right) \cos(2(\theta - \vartheta' - q_k \gamma)) \\
& - \frac{3}{2} i \Gamma_k \left( \bar{y}_k^2 - y_k'^2 \right) \sin(2q_k \gamma') \\
& - \frac{3}{2} i \Gamma_k' \left( \bar{y}_k^2 - y_k^2 \right) \sin(2q_k \gamma) \\
& + 3i \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k' y_k - \bar{y}_k y_k' \right) \sin(\theta - \vartheta') \\
& + 3i \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k \bar{y}_k' - y_k y_k' \right) \sin(\theta - \vartheta' + 2q_k \gamma') \\
& - 3i \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k \bar{y}_k' - y_k y_k' \right) \sin(\theta - \vartheta' - 2q_k \gamma) \\
& - 3i \sqrt{\Gamma_k \Gamma_k'} \left( \bar{y}_k' y_k - \bar{y}_k y_k' \right) \sin(\theta - \vartheta' - 2q_k(\gamma - \gamma')) \\
& - \frac{3}{2} i \Gamma_k' \left( \bar{y}_k^2 - y_k^2 \right) \sin(2(\theta - \vartheta' + q_k \gamma')) \\
& + \frac{3}{2} i \Gamma_k \left( \bar{y}_k^2 - y_k'^2 \right) \sin(2(\theta - \vartheta' - q_k \gamma)) \Biggr], \nonumber
\end{align}
$$

with $p_1 = -p/2$, $p_2 = 1 + p/2$, $q_1 = 1 + 2/p$, and $q_2 = 1$.

We note that $\sigma$ does not appear in the expression of $U_k$. Therefore, in the presence of tides the parameter $\Sigma$ (Eq. (29)) remains conserved. The fast angle $\gamma$ is still present in the expression of $U_k$, but at this stage we cannot perform an average as in Sect. 2.1.3 because $\gamma'$ is considered as a parameter that can later cancel with $\gamma$ (see Eq. (64)).

## 3.2. Secular equations of motion

The equations of motion are obtained from Eq. (58) using the Hamilton equations. The additional contributions from tides only

derive from the tidal potential energy $U_k$, and are given by

$$ \dot{y}_k = i \sum_{j=1}^{2} \frac{\partial U_j}{\partial \bar{y}_j}, \quad \dot{\Gamma} = - \sum_{j=1}^{2} \frac{\partial U_j}{\partial \gamma}, \quad \dot{\Theta} = - \sum_{j=1}^{2} \frac{\partial U_j}{\partial \theta}. \quad (61) $$

In principle, we should also write the equations for $\dot{\gamma}$ and $\dot{\theta}$, but these angles disappear from the equations of motion with some of the following simplifications, and so we do not need them to get a closed set for the secular evolution of the system.

To handle the expression of the primed quantities, we need to use a tidal model. For simplicity, we adopt here the weak friction model (e.g. Singer 1968; Alexander 1973), which assumes a constant and small time delay, $\tau$. This model is widely used and provides very simple expressions for the tidal interactions, because it can be made linear (e.g. Mignard 1979):

$$ \chi_k' \approx \lambda_k - n_k \tau \quad \text{and} \quad \theta' \approx \theta - \omega \tau. \quad (62) $$

It follows for the remaining primed quantities that

$$ y_k' \approx y_k - i y_k (p_2 \bar{n}_2 + p_1 \bar{n}_1) \tau, \quad (63) $$

$$ \gamma' \approx \gamma - p_1 (\bar{n}_2 - \bar{n}_1) \tau, \quad (64) $$

$$ \vartheta' \approx \vartheta + (p_2 \bar{n}_2 + p_1 \bar{n}_1 - \omega) \tau, \quad (65) $$

with

$$ \bar{n}_k = \beta_k^3 \mu_k^2 / \Gamma_k^3. \quad (66) $$

We then substituted expressions (63) to (65) into the equations of motion (61) and averaged over the fast angle $\gamma$ (as in Eq. (30)) to finally get the secular equations for the tidal evolution

$$ \dot{y}_1 = -\frac{3}{2} \frac{\mathcal{D}_1}{\Gamma_1^{13}} (2ip + \bar{n}_1\tau) y_1 + 3i(p + 2) \frac{\mathcal{D}_2}{\Gamma_2^{13}} y_1, \quad (67) $$

$$ \dot{y}_2 = -\frac{3}{2} \frac{\mathcal{D}_2}{\Gamma_2^{13}} (-2i(p+2) + \bar{n}_2\tau) y_2 - 3ip \frac{\mathcal{D}_1}{\Gamma_1^{13}} y_2, \quad (68) $$

$$
\begin{align}
\dot{\Gamma} ={}& 3 \frac{\mathcal{D}_1}{\Gamma_1^{13}} (1 + 2/p) \Biggl[ & (\Gamma_1 + 6p(y_1 \bar{y}_1 + y_2 \bar{y}_2) - y_1 \bar{y}_1) \omega \nonumber \\
               & & - (\Gamma_1 + 6p(y_1 \bar{y}_1 + y_2 \bar{y}_2)) \bar{n}_1 \Biggr] \tau \tag{69} \\
               & & + 3 \frac{\mathcal{D}_2}{\Gamma_2^{13}} \Biggl[ & (\Gamma_2 - 6(p + 2)(y_1 \bar{y}_1 + y_2 \bar{y}_2) - y_2 \bar{y}_2) \omega \nonumber \\
               & & - (\Gamma_2 - 6(p + 2)(y_1 \bar{y}_1 + y_2 \bar{y}_2)) \bar{n}_2 \Biggr] \tau, \nonumber
\end{align}
$$

$$
\begin{align}
\dot{\Theta} ={}& -3 \frac{\mathcal{D}_1}{\Gamma_1^{13}} \Biggl[ (\Gamma_1 + 6p(y_1 \bar{y}_1 + y_2 \bar{y}_2) - y_1 \bar{y}_1) \omega \nonumber \\
& - (\Gamma_1 + 6p(y_1 \bar{y}_1 + y_2 \bar{y}_2) - y_1 \bar{y}_1) n_1 \Biggr] \tau \nonumber \\
& -3 \frac{\mathcal{D}_2}{\Gamma_2^{13}} \Biggl[ (\Gamma_2 - 6(p + 2)(y_1 \bar{y}_1 + y_2 \bar{y}_2) - y_2 \bar{y}_2) \omega \nonumber \\
& - (\Gamma_2 - 6(p + 2)(y_1 \bar{y}_1 + y_2 \bar{y}_2) - y_2 \bar{y}_2) n_2 \Biggr] \tau, \tag{70}
\end{align}
$$

where

$$ \mathcal{D}_k = k_2 \mathcal{G} m_k^2 \beta_k^{12} \mu_k^6 R^5. \quad (71) $$

We note that in the expressions of $y_1$ (Eq. (67)) and $y_2$ (Eq. (68)), we have a conservative contribution (imaginary terms) and a dissipative contribution (real terms in $\tau$). The conservative contributions result from a permanent tidal deformation and only slightly modify the fundamental frequencies of the system, while the dissipative contributions modify the secular evolution.
### 3.3. Tidal constraints

Tidal dissipation induces variations in the parameters $\Gamma$ (Eq. (69)) and $\Theta$ (Eq. (70)). Then, the coefficients $K$, $O$, $S$, and $R$ appearing in the Hamiltonian (34) slowly change in time (Eqs. (21)–(23)), modifying the phase space (Fig. 2), which translates into a secular evolution of the system. We note that for the oblateness coefficients (Eqs. (A.3) and (A.4)), changes are observed not only due to $\Gamma$, but also in $J_2$, according to (e.g. Correia & Rodríguez 2013)

$$J_2 = k_f \frac{\omega^2 R^3}{3Gm_0}, \qquad (72)$$

where $\omega = \Theta/C$ (Eq. (7)) and $k_f$ is the fluid second Love number for potential. Using the present rotation rate and $J_2$ of Uranus (Table 1), we obtained $k_f = 0.356$.

#### 3.3.1. Orbital evolution

For a better understanding of the tidal evolution of the satellites' orbits, we can compute the secular evolution of the semi-major axes and inclinations from Eqs. (17) and (18) as

$$\dot{a}_k \approx 2A_k \left( \frac{\omega}{n_k} \left(1 - \frac{1}{2}I_k^2\right) - 1 \right) a_k, \qquad (73)$$

$$\dot{I}_k \approx -\frac{A_k \omega}{2 n_k} I_k, \qquad (74)$$

with

$$A_k = \frac{3Gm_k^2 R^5}{\beta_k a_k^8} k_2 \tau. \qquad (75)$$

Both Ariel and Umbriel have very low inclinations (Table 1). Combined with the fast rotation rate of Uranus ($\omega/n_k \gg 1$), we conclude that tides induce an outward migration of the satellites (Eq. (73)), and they damp their inclinations (Eq. (74)).

#### 3.3.2. Resonance encounter

We can extrapolate the past evolution of the semi-major axes and determine when the 5/3 MMR encounter may have occurred (Eq. (73)). Neglecting the effect of the inclination and assuming a constant rotation rate for Uranus ($\omega/n_k \approx \text{cte}$), we get

$$\frac{a_2^7 da_2}{m_2 (\omega/n_1 - 1)} \approx \frac{a_1^7 da_1}{m_1 (\omega/n_2 - 1)}, \qquad (76)$$

which can be integrated to obtain $a_2$ as a function of $a_1$. Starting with the present system (Table 1) and assuming that the satellites were not temporarily captured into resonance, we can move backwards until the ratio of the nominal resonance is achieved (Eq. (41))

$$a_2/a_1 \approx (5/3)^{2/3} \approx 1.4057. \qquad (77)$$

This allowed us to get the best estimate for the semi-major axes of Ariel and Umbriel at the exact 5/3 MMR (Fig. 4),

$$a_1/R = 7.3906, \quad a_2/R = 10.3891, \qquad (78)$$

and to compute other related parameters, such as $\Gamma$ (Eq. (38)).

<figure><img src="image_6.png" /><figcaption>Fig. 4. Tidal evolution of the semi-major axes of Ariel and Umbriel (Eq. (76)). The orange point marks the current observed values (Table 1). The dashed line gives the position of the nominal 5/3 MMR relation between $a_2$ and $a_1$ (Eq. (77)). The blue point gives the best estimation for the semi-major axes at the nominal resonance (Eq. (78)).</figcaption></figure>

#### 3.3.3. Evolution timescale

The orbital evolution timescale depends on $\tau$, which is related to the tidal dissipation (Eq. (62)). A more commonly used dimensionless quantity to measure the tidal dissipation is given by the quality factor (e.g. Correia & Valente 2022),

$$Q = 1/(2\omega\tau). \qquad (79)$$

As a result of studying the likelihood of resonance crossing in the Uranian system, Tittemore & Wisdom (1990) stated that the 2/1 MMR Ariel-Umbriel cannot be crossed, while the 3/1 MMR Miranda-Umbriel can be, and thus they constrained the interval of $Q$ to be between 11 000 and 39 000. Following the same approach, but using more recent measures for the satellites' masses (Jacobson 2014), we find that $Q = 8000$ is a more suitable value. From Eq. (79), we then computed $\tau \approx 0.62$ s. In fact, the exact dissipation rate depends on the product $k_2\tau$ (Eq. (75)). As in previous studies, in the present work we adopted $k_2 = 0.104$ (Gavrilov & Zharkov 1977), which translates into

$$\frac{k_2}{Q} = 1.3 \times 10^{-5} \Leftrightarrow k_2\tau = 0.064 \text{ s.} \qquad (80)$$

Using Eq. (73), we estimate that the 5/3 MMR was crossed about 640 Myr ago. We can also estimate the damping timescale for the inclinations of both satellites (Eq. (74)). We get

$$\tau_{\text{inc}} \approx \frac{4\beta_k a_k^8 n_k Q}{3Gm_k^2 R^5 k_2}, \qquad (81)$$

which yields about 180 Gyr and 1500 Gyr for Ariel and Umbriel, respectively. We hence conclude that the presently observed inclination values are likely unchanged since the system crossed the 5/3 MMR.

### 3.4. Capture probability

The behaviour of the system when it crosses the 5/3 MMR is not straightforward, since the problem has two degrees of freedom. As explained in Sect. 2.2.3, the different behaviours while crossing the resonance depend on the energy of the system, which corresponds to different inclinations. In order to get an idea of the critical inclinations that tend to trap the system in resonance or skip it, we can build a simplified one degree-of-freedom model.
The system can either be in resonance with the angle $\varphi_1$ or $\varphi_2$. Following Tittemore & Wisdom (1988), we retained only the terms associated with each angle and obtained their associated one degree-of-freedom simplified Hamiltonian (Eq. (34))

$$ \mathcal{H}_k = \mathcal{A}_k y_k \bar{y}_k + \mathcal{B}_k (y_k \bar{y}_k)^2 + \frac{C_k}{2} (y_k^2 + \bar{y}_k^2), \qquad (82) $$

with

$$ \mathcal{A}_1 = \mathcal{K}_a + O_a + S_a + S_b, \quad \mathcal{B}_1 = \mathcal{K}_b, \quad C_1 = \mathcal{R}_a, \qquad (83) $$

$$ \mathcal{A}_2 = \mathcal{K}_a + O_b + S_a + S_c, \quad \mathcal{B}_2 = \mathcal{K}_b, \quad C_2 = \mathcal{R}_b. \qquad (84) $$

In this simplified case, we could then directly apply the theory of the adiabatic invariant for resonance capture (Henrard 1982; Henrard & Lemaître 1983). In general, the phase space of a given resonance presents three different regions delimited by the separatrix (see for instance Fig. 2c): a libration region, with area $\mathcal{J}_k^R$, and two circulation regions, one outside the libration region and another encircled by the separatrix, with an area $\mathcal{J}_k^C$. Consider an initial trajectory with some energy $H_k < H_0$. As the system evolves due to tides, the phase space and the areas of each region also change. When the trajectory encounters the resonance, depending on the exact place where the separatrix is crossed, the system can either evolve into the libration region (resonance) or into the inner circulation region. We can compute the capture probability analytically, provided that the evolution is adiabatic, that is, the tidal induced variations are much slower than the conservative inclination variations. The capture probability can then be obtained by the modification of the phase space with time, that is, by the change in the areas encircled by the separatrix (e.g. Yoder 1979; Henrard 1982):

$$ P_{\text{cap}} = \frac{\dot{\mathcal{J}}_k^R}{\dot{\mathcal{J}}_k^R + \dot{\mathcal{J}}_k^C}. \qquad (85) $$

The area of the resonance region is obtained by integrating over the separatrix,

$$ \mathcal{J}_k^R = i \oint y_k d\bar{y}_k. \qquad (86) $$

The energy of the separatrix ($H_0$) can be found by finding the hyperbolic points (see Sect. 2.2.1) and computing the Hamiltonian (82) at these points,

$$ H_0 = - \frac{(\mathcal{A}_k + C_k)^2}{4\mathcal{B}_k}. \qquad (87) $$

Replacing $H_0$ into expression (82), we find, for the separatrix points, that

$$ y_k = - \frac{2\mathcal{A}_k \bar{y}_k \pm \sqrt{-2C_k (\mathcal{A}_k + 2\bar{y}_k^2 \mathcal{B}_k + C_k)^2 / \mathcal{B}_k}}{4\bar{y}_k^2 \mathcal{B}_k + 2C_k}. \qquad (88) $$

There are two solutions, the + corresponding to the branch between the libration and the outer circulation region ($\mathcal{J}_k^+$), and the - corresponding to the branch between the libration and the inner circulation region ($\mathcal{J}_k^-$). Therefore, we have $\mathcal{J}_k^+ = \mathcal{J}_k^+ + \mathcal{J}_k^-$ and $\mathcal{J}_k^+ = -\mathcal{J}_k^-$, which gives (Eq. (85))

$$ P_{\text{cap}} = \frac{\dot{\mathcal{J}}_k^+ + \dot{\mathcal{J}}_k^-}{\dot{\mathcal{J}}_k^+} = 1 + \frac{\partial \mathcal{J}_k^-}{\partial \Gamma} \bigg|_{\partial \Gamma}, \qquad (89) $$

<figure><img src="image_7.png" /><figcaption>Fig. 5. Capture probabilities in the $\varphi_1$ (top) and $\varphi_2$ (bottom) resonances. The solid line gives the theoretical approximation given by expression (89), while the dots give the results of numerical simulations. We ran 100 initial conditions with $\varphi_k$ differing by 1.8°.</figcaption></figure>

where,

$$ \frac{\partial \mathcal{J}_k^-}{\partial \Gamma} = \frac{1}{\mathcal{B}_k^2} \left( \arcsin \left( \sqrt{-\frac{C_k}{\mathcal{A}_k}} \right) + \frac{\pi}{2} \right) \left( \mathcal{B}_k \mathcal{A}_k' - \mathcal{A}_k \mathcal{B}_k' \right) + \frac{1}{\mathcal{B}_k^2} \sqrt{-\frac{C_k + \mathcal{A}_k}{C_k}} \left( \mathcal{B}_k \mathcal{C}_k' - C_k \mathcal{B}_k' \right), \qquad (90) $$

$$ \frac{\partial \mathcal{J}_k^+}{\partial \Gamma} = \frac{1}{\mathcal{B}_k^2} \left( \arcsin \left( \sqrt{-\frac{C_k}{\mathcal{A}_k}} \right) - \frac{\pi}{2} \right) \left( \mathcal{B}_k \mathcal{A}_k' - \mathcal{A}_k \mathcal{B}_k' \right) + \frac{1}{\mathcal{B}_k^2} \sqrt{-\frac{C_k + \mathcal{A}_k}{C_k}} \left( \mathcal{B}_k \mathcal{C}_k' - C_k \mathcal{B}_k' \right), \qquad (91) $$

and $\mathcal{A}_k' = \partial \mathcal{A}_k / \partial \Gamma$, $\mathcal{B}_k' = \partial \mathcal{B}_k / \partial \Gamma$, and $C_k' = \partial C_k / \partial \Gamma$. We note that since $\Gamma$ is the only time-dependent quantity appearing in $\mathcal{A}_k$, $\mathcal{B}_k$, and $C_k$ (Appendix A), $\dot{\mathcal{J}}_k = \dot{\Gamma} \partial \mathcal{J}_k / \partial \Gamma$ (we neglected the small changes in $\Theta$ from the oblateness coefficients, $O_k$).

In Fig. 5, we show the probability of capture in the $\varphi_1$ and $\varphi_2$ resonances obtained with expression (89). For some inclination values, we also show the results obtained with numerical integrations of the equations of motion derived from the simplified Hamiltonian (Eq. (82)) together with the secular tidal equations (Eqs. (67)–(70)). For each initial inclination, we ran 100 simulations where the initial angle $\varphi_k$ was uniformly sampled. The amount of simulations captured in resonance at the end of the simulation are marked with a dot. The statistical fluctuation, represented as error bars, were estimated using binomial statistics,
with the expression

$$ \Delta P = \sqrt{\frac{P_{\text{cap}}}{N} (1 - P_{\text{cap}})}, \qquad (92) $$

where N = 100 is the number of simulations (e.g. Tittemore & Wisdom 1988). We observe that there is a good agreement between the theoretical curve (Eq. (89)) and the output of the numerical simulations, that is, the adiabatic approximation holds.

In Fig. 5, we observe that for initial inclinations lower than 0.05°, the system is always captured in resonance. However, as we increase the initial inclination, the capture probability quickly decreases, it becomes ~50% for $I_k = 0.1°$, and drops to ~30% for $I_k = 0.2°$. These results suggest that a system with nearly coplanar orbits cannot escape the 5/3 MMR, but for inclinations higher than about 0.1°, it may be able to evade it.

We cannot completely rely on the conclusions obtained with the simplified Hamiltonian, mainly for two reasons. One is because the complete Hamiltonian (Eq. (34)) depends on the inclination of the other body. When we simplified the Hamiltonian (Eq. (82)) for $y_1$, we dropped all terms in $y_2$ (and vice versa), which is equivalent to setting $y_2 = 0$. However, if we set $y_2 \neq 0$, more terms appear in the Hamiltonian, leading to a different distribution in the capture probabilities. The other reason is that the complete Hamiltonian has two degrees of freedom, and so for some combinations of the inclination values, the system can be chaotic (see Sect. 2.2.3). For the trajectories crossing the chaotic regions, the final outcome is unpredictable.

# 4. Numerical simulations

We have seen that Ariel and Umbriel almost certainly encountered the 5/3 MMR at some time in the past. Depending on the inclination values of these satellites, the system may experience rather different behaviours. Near coplanar orbits are expected to become trapped in resonance, while some inclination can lead to alternative scenarios. In order to get a global view of the resonant dynamics with tides, in this section we integrate the complete set of differential Eqs. (36), (37), and (67)–(70).

## 4.1. Setup

When a resonance is crossed, we cannot perform a backward integration because we have a stochastic evolution. Therefore, we need to place the system slightly before the resonance encounter and then integrate it forwards. It is not possible to determine the exact semi-major axes prior to resonance crossing, but if the system does not spend much time in resonance, the semi-major axes should not differ much from the estimation given in Sect. 3.3.2. We still need to slightly decrease $a_1$ (or increase $a_2$) to move the system out of the nominal resonance (Eq. (78)). Since tides are stronger in Ariel, we opted to shift $a_1$ and kept $a_2$ constant:

$$ a_1/R = 7.3892, \quad a_2/R = 10.3891. \qquad (93) $$

These values of the semi-major axes allowed us to compute the initial $\Gamma$ parameter (Eq. (16)). For some given initial inclination values, we could also compute the initial value of $\Delta$ (Eq. (23)), which translates into an initial $\delta > 0$ (Eq. (39)).

The physical properties of Uranus and its satellites can be found in Table 1. The total angular momentum is conserved and can also be obtained from the present system,

$$ \Sigma = 9.367 \times 247 \times 10^{-10} M_{\odot} \text{au}^2 \text{yr}^{-1}. \qquad (94) $$

<figure><img src="image_8.png" /><figcaption>Fig. 6. Tidal evolution of the system as function of $\delta$ for initial inclinations $I_1 = I_2 = 0.01°$. The results of the numerical simulation (in black) are superimposed in the equilibria map (Fig. 1). We show the evolution for the resonant angle $\varphi_1$ (top) and $\varphi_2$ (bottom).</figcaption></figure>

The initial rotational rate, $\omega$, was obtained from $\Sigma$ using Eq. (29) with the pre-resonance semi-major axes (Eq. (93)). Finally, for the tidal dissipation, we adopted $k_2\tau = 0.064$ s (Eq. (80)).

## 4.2. Comparison with analytical estimations

In general, tidal effects are weak and only correspond to small perturbations of the conservative dynamics (adiabatic approximation). To verify that the numerical integrations do follow the theoretical predictions from Sects. 2 and 3, we initially ran a few simulations with the full equations of motion and then superimposed the output in the equilibria map as a function of $\delta$ (Fig. 1). Since $\omega > n_k$, tidal effects are expected to increase the value of $\Gamma$ (Eq. (69)) and thus decrease the value of $\delta$ (Eq. (39)). Therefore, the results of the simulations as a function of time must be read from the right to the left in the figure.

In Fig. 6, we show the results of a first experiment with very small initial inclinations for both satellites, $I_1 = I_2 = 0.01°$. Initially, when $\delta > 0$, the system is in circulation with a small amplitude around the equilibrium point at zero ($y_1 = 0, y_2 = 0$). The system encounters the resonance when $\delta \approx 0$. However, just before $\delta = 0$, two stable equilibrium points emerge (corresponding to the 5/3 MMR), and the equilibrium point at zero becomes unstable. The system is thus forced to follow one of the two resonance branches. As the system evolves and $\delta < 0$, the equilibrium point at zero becomes stable again. However, because the amplitude of oscillation is small, the system closely follows the resonant branch, which corresponds to an increase in the inclinations. We thus confirm that for initially near coplanar orbits, the system cannot avoid capture in the 5/3 MMR (Sect. 3.4).

In Fig. 7, we show the results of a second experiment with higher initial inclinations for both satellites, $I_1 = I_2 = 0.1°$. The initial evolution for $\delta > 0$ is similar to the case with lower initial inclinations (Fig. 6), except that the amplitude of oscillation
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Two examples of tidal evolution of the system as a function of $\delta$ for initial inclinations $I_1 = I_2 = 0.1^{\circ}$. The results of the numerical simulation (in black) are superimposed in the equilibria map (Fig. 1). We show the evolution for the resonant angle $\varphi_1$ (top) and $\varphi_2$ (bottom).</figcaption></figure>

is ten times larger in this case (corresponding to a system with a higher energy). As the system encounters the resonance at $\delta \approx 0$, it is not able to follow one of the resonant equilibria and it remains in a chaotic region around the separatrix (Fig. 3f). In the example on the left (Fig. 7a), after some time in the chaotic region with $\delta < 0$, the system finds a way into the libration region and follows one of the resonant branches in a quasi-periodic orbit (Fig. 3i). On the other hand, in the example on the right (Fig. 7b), the system finds an alternative path back into the circulation region around the equilibrium point at zero in a quasi-periodic orbit (Fig. 3k). We thus confirm that for initial orbits with some inclination, the system experiences a chaotic regime for some time (Sect. 2.2.3), after which it can either be captured in the 5/3 MMR, or escape it with some probability (Sect. 3.4).

### 4.3. Impact of the initial inclination

The present mean inclinations of Ariel and Umbriel (Table 1) are likely unchanged after the system quits the 5/3 MMR (Eq. (81)). However, they may have been rather different before this encounter, as the inclinations undergo some excitations while crossing the resonance (Sect. 2). The fact that the system is not trapped in resonance at present, implies that the separatrix was crossed at some point and the inclinations had to experience some chaotic oscillations (Fig. 7b). Therefore, it is impossible to simply integrate backwards and determine the exact inclination values prior to the resonance encounter. To have a more clear idea of what the system may have been, we need to perform a statistical study of its past evolution, starting with arbitrary initial inclinations for both Ariel and Umbriel and then reject those that are not coherent with the present observations.

For that purpose, we explored a mesh of initial inclinations ranging between $0.001^{\circ}$ and $0.2^{\circ}$ with a stepsize of $0.05^{\circ}$. For each pair ($I_1, I_2$), we ran 1000 simulations evenly sampled over the angle $\sigma$ (Eq. (10)) for 100 Myr, in a total of 25 000 experiments. In Table 2, we list the complete set of initial

conditions together with a summary of the outcome of the resonance crossing. The results of the simulations are presented in percentage because they can be seen as a statistical distribution of a given final evolution possibility. The total number of events of each kind can be simply obtained multiplying the percent by ten.

#### 4.3.1. Capture probability

For each run, we first evaluated whether capture in resonance occurred or not. Capture takes place whenever at least one resonant angle, $\varphi_k$ (Eq. (11)), switches from circulation to libration, and the mean motion ratio becomes approximately constant ($n_1/n_2 \approx 5/3$). The mean motion ratio criterion is usually more useful to automatically detect captures because inside chaotic regions the behaviour of the resonant angles can be quite erratic. Conversely, capture does not occur when all resonant angles remain in circulation, which introduces only a small jump in the mean motion ratios. For this last case, we consider that the resonance is skipped. In Table 2, the relative number of these events is dubbed as $P_s$ (skip probability).

The capture probability can be simply evaluated as $1 - P_s$. However, when a capture occurs, we still have to distinguish between trajectories that remain captured for long periods of time from those that are able to escape shortly after. Indeed, in many resonant lockings, the inclinations remain chaotic and the system is able to evade the resonance after some time. In theory, all trajectories evading the resonance could lead to the present system. In practice, that is not possible, because the longer the system stays in resonance, the higher the final inclinations become (Fig. 6), and they cannot be conciliated with the presently observed values (Table 1). Moreover, Ćuk et al. (2020) have shown that as long as Ariel and Umbriel stay in resonance, the inclinations of the three other large satellites of Uranus also grow, in particular that of Miranda. After only 10 Myr in resonance, is it impossible to explain the observed values since tides
<table><caption>**Table 2.** Initial conditions and summary of the numerical simulations's results of the 5/3 MMR crossing.</caption><thead><tr><th>#</th><th>I<sub>1</sub> (°)</th><th>I<sub>2</sub> (°)</th><th>P<sub>c</sub> (%)</th><th>P<sub>e</sub> (%)</th><th>P<sub>s</sub> (%)</th><th>⟨I<sub>1</sub>⟩ (°)</th><th>ΔI<sub>1</sub> (°)</th><th>⟨I<sub>2</sub>⟩ (°)</th><th>ΔI<sub>2</sub> (°)</th></tr></thead><tbody><tr><td>1</td><td>0.001</td><td>0.001</td><td>100.0</td><td>—</td><td>—</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>0.05</td><td>0.001</td><td>100.0</td><td>—</td><td>—</td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>0.10</td><td>0.001</td><td>99.2</td><td>0.8</td><td>—</td><td></td><td></td><td></td><td></td></tr><tr><td>4</td><td>0.15</td><td>0.001</td><td>100.0</td><td>—</td><td>—</td><td></td><td></td><td></td><td></td></tr><tr><td>5</td><td>0.20</td><td>0.001</td><td>100.0</td><td>—</td><td>—</td><td></td><td></td><td></td><td></td></tr><tr><td>6</td><td>0.001</td><td>0.05</td><td>100.0</td><td>—</td><td>—</td><td></td><td></td><td></td><td></td></tr><tr><td>7</td><td>0.05</td><td>0.05</td><td>100.0</td><td>—</td><td>—</td><td></td><td></td><td></td><td></td></tr><tr><td>8</td><td>0.10</td><td>0.05</td><td>97.3</td><td>2.5</td><td>0.2</td><td></td><td></td><td></td><td></td></tr><tr><td>9</td><td>0.15</td><td>0.05</td><td>82.4</td><td>10.5</td><td>7.1</td><td>0.149</td><td>0.062</td><td>0.056</td><td>0.023</td></tr><tr><td>10</td><td>0.20</td><td>0.05</td><td>80.9</td><td>9.3</td><td>9.8</td><td>0.187</td><td>0.054</td><td>0.047</td><td>0.019</td></tr><tr><td>11</td><td>0.001</td><td>0.10</td><td>90.1</td><td>9.2</td><td>0.7</td><td>0.170</td><td>0.041</td><td>0.047</td><td>0.016</td></tr><tr><td>12</td><td>0.05</td><td>0.10</td><td>83.1</td><td>13.6</td><td>3.3</td><td>0.167</td><td>0.050</td><td>0.041</td><td>0.022</td></tr><tr><td>13</td><td>0.10</td><td>0.10</td><td>50.6</td><td>28.5</td><td>20.9</td><td>0.145</td><td>0.059</td><td>0.049</td><td>0.020</td></tr><tr><td>14</td><td>0.15</td><td>0.10</td><td>38.8</td><td>19.2</td><td>42.0</td><td>0.119</td><td>0.068</td><td>0.068</td><td>0.020</td></tr><tr><td>15</td><td>0.20</td><td>0.10</td><td>36.2</td><td>15.8</td><td>48.0</td><td>0.141</td><td>0.066</td><td>0.081</td><td>0.023</td></tr><tr><td>16</td><td>0.001</td><td>0.15</td><td>39.6</td><td>26.5</td><td>33.9</td><td>0.112</td><td>0.068</td><td>0.083</td><td>0.019</td></tr><tr><td>17</td><td>0.05</td><td>0.15</td><td>38.6</td><td>25.8</td><td>35.6</td><td>0.114</td><td>0.069</td><td>0.095</td><td>0.024</td></tr><tr><td>18</td><td>0.10</td><td>0.15</td><td>35.1</td><td>25.1</td><td>39.8</td><td>0.088</td><td>0.081</td><td>0.109</td><td>0.017</td></tr><tr><td>19</td><td>0.15</td><td>0.15</td><td>31.0</td><td>23.2</td><td>45.8</td><td>0.099</td><td>0.082</td><td>0.127</td><td>0.023</td></tr><tr><td>20</td><td>0.20</td><td>0.15</td><td>30.0</td><td>15.9</td><td>54.1</td><td>0.104</td><td>0.085</td><td>0.136</td><td>0.021</td></tr><tr><td>21</td><td>0.001</td><td>0.20</td><td>32.7</td><td>26.2</td><td>41.1</td><td>0.107</td><td>0.073</td><td>0.144</td><td>0.018</td></tr><tr><td>22</td><td>0.05</td><td>0.20</td><td>31.9</td><td>28.6</td><td>39.5</td><td>0.116</td><td>0.070</td><td>0.152</td><td>0.018</td></tr><tr><td>23</td><td>0.10</td><td>0.20</td><td>30.1</td><td>23.9</td><td>46.0</td><td>0.109</td><td>0.072</td><td>0.162</td><td>0.015</td></tr><tr><td>24</td><td>0.15</td><td>0.20</td><td>27.1</td><td>22.1</td><td>50.8</td><td>0.117</td><td>0.075</td><td>0.173</td><td>0.013</td></tr><tr><td>25</td><td>0.20</td><td>0.20</td><td>27.0</td><td>15.2</td><td>57.8</td><td>0.119</td><td>0.076</td><td>0.182</td><td>0.014</td></tr></tbody></table>

**Notes.** $I_k$ is the initial inclination of each satellite; $P_c$ is the number of simulations trapped in resonance more than 10 Myr; $P_e$ is the number of simulations that were captured in resonance, but escaped in less than 10 Myr; $P_s$ is the number of simulations that skipped the resonance; and $\langle I_k \rangle$ and $\Delta I_k$ are the mean and the standard deviation, respectively, of a Rice distribution adjusted to the final results.

are unable to efficiently damp the inclinations (Eq. (81)). Therefore, in our analysis, we split the capture events into those that stay in resonance for more than 10 Myr and those that are able to evade it before that time. In Table 2, the relative number of escapes is dubbed as $P_e$ (escape probability) and the relative number of long-term entrapments is referred to as $P_c$ (capture for more than 10 Myr probability).

In most sets of the simulations shown in Table 2, all three scenarios described above (capture, escape, and skip) are simultaneously present. To better illustrate the different possibilities, in Fig. 8, we show one example of each case, corresponding to set #17 (with an initial $I_1 = 0.05^\circ$ and $I_2 = 0.15^\circ$).

In Fig. 8a, the system is permanently captured in resonance. Prior to the resonance encounter, the semi-major axes ratio follows the asymptotic evolution predicted by Eq. (73). When the system comes across the resonance (at $t \approx 11$ Myr), the semi-major axes ratio becomes constant (since $n_1/n_2 \approx 5/3$), deviating considerably from the asymptotic line. Indeed, looking at Fig. 4, when capture occurs, instead of following the black curve, the system follows the dashed line, corresponding to the nominal resonance (Eq. (77)). Shortly after being captured, the system enters into the chaotic region. The inclinations of both satellites are excited and grow steadily on average, in particular that of Ariel, which oscillates between the two resonant branches (green lines in Fig. 1). For a time $t \approx 40$ Myr, the system finds a way into the libration region and becomes quasi-periodic, as in the example shown in Fig. 7a.

In Fig. 8b, the system is temporarily captured in resonance for less than 10 Myr. As in the previous example, the semi-major

axes ratio initially follows the asymptotic evolution predicted by Eq. (73), but it switches to a constant ratio as soon as capture in resonance occurs (at $t \approx 11$ Myr). As before, during the resonant entrapment, the system enters into the chaotic region and the inclinations of both satellites start to grow. However, for a time around $t \approx 18$ Myr (that is, just 7 Myr after being captured), the system finds a way outside the chaotic region that breaks the resonant locking and returns into the circulation region around the equilibrium point at zero (blue line in Fig. 1). From that point on, the semi-major axes ratio again follows the asymptotic evolution predicted by Eq. (73), though restarting with a slightly higher ratio, and the inclinations of both satellites stabilise around a given constant mean value. The final inclination of Ariel is always higher than its initial value because of the forcing during the resonant phase. The final inclination of Umbriel is less impacted by this mechanism and it does not have a systematic trend because it oscillates around the initial value with a large amplitude during the resonance crossing.

Finally, in Fig. 8c, the system shortly skips the resonance. As usual, the semi-major axes ratio initially follows the asymptotic evolution predicted by Eq. (73), and it undergoes some perturbations during the resonance crossing. However, it quickly returns to the asymptotic evolution, though restarting with a slightly lower ratio. In this case, there is almost no chaotic evolution for the inclinations because, just after crossing the separatrix, the system directly goes into the circulation region around the equilibrium point at zero (blue line in Fig. 1). The inclination of Ariel initially briefly grows since, owing to the topology of the 5/3 MMR (Fig. 2), it is not possible to reach the inner circulation
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Three examples of a system crossing the 5/3 MMR with initial $I_1 = 0.05^\circ$ and $I_2 = 0.15^\circ$ (set #17 in Table 2). We show the ratio between the semi-major axes (top), the inclination of Ariel (middle), and the inclination of Umbriel (bottom) as a function of time. Each column corresponds to a different simulation. We show an example of a system that is permanently caught in resonance (a), one that is captured but evades the resonance in less than 10 Myr (b), and another that skips the resonance without capture (c). The blue line gives the asymptotic evolution predicted by Eq. (73), while the dashed line gives the position of the nominal resonance Eq. (77). The orange dots mark the present mean inclinations (Table 1).</figcaption></figure>

region without moving around the resonant equilibrium for a short time (less than 1 Myr). Since in the example #17 the initial inclination of Ariel is indeed relatively small, for a short moment the inclination of Ariel follows the resonant branch and has to increase (green line in Fig. 1). Nevertheless, as soon as the equilibrium point at zero becomes stable, the system moves back into the circulation region and the inclination of Ariel drops. In the example shown in Fig. 8c, the inclination of Ariel decreases to a mean value smaller than the initial one, but actually in other examples it can be anything between zero and the maximum previously attained. The brief resonant excitation only involves the angle $\phi_1$, and so it does not impact the inclination of Umbriel much. Its amplitude grows due to the mutual interactions with Ariel, but the mean value remains constant. However, as the resonance is skipped, the mean value of Umbriel inclination suddenly drops to a lower level. This reduction is always observed because the inner circulation region is confined within a region of low inclination for Umbriel, in particular for $\delta$ values very close to zero (Fig. 2c).

Figure 8c also provides an example of a simulation where the system crosses the 5/3 MMR and subsequently evolves into the presently observed configuration (orange dots, taken from Table 1). Not all simulations #17 that skip the resonance end in the present system, though a significant number is consistent with it. On the other hand, some simulations starting with different initial conditions (Table 2) can also evolve into the present state. Therefore, only a more refined analysis of the final distribution of the inclinations of both Ariel and Umbriel

can provide more insight into the occurrence likelihood of the currently observed system.

### 4.3.2. Statistics

For the set of simulations #1–#8 (Table 2), that is to say those with initial $I_2 \le 0.05^\circ$, the system is captured in resonance nearly in 100% of the cases. For the initial inclination of Umbriel, these results are thus in agreement with the predictions of the simplified model presented in Sect. 3.4. Since capture in the 5/3 MMR is not observed today, we can immediately exclude this range of initial conditions. We hence conclude that, regardless of the initial inclination of Ariel, the system requires some moderate initial inclination for Umbriel ($I_2 \ge 0.1^\circ$) to evolve into the present state.

The resonance crossing is a stochastic process and therefore the same initial inclinations with a slightly different initial resonant angle may end up in a completely different final state. In all of the other sets of simulations that we performed (#9–#25), the system can either remain captured for a long time, escape the resonance in less than 10 Myr, or simply skip it. In general, as we increase the initial inclinations of both satellites, the probability of capture in resonance decreases, in conformity with the analysis from Sect. 3.4. Interestingly, the number of trajectories temporarily captured for less than 10 Myr does not change much with the initial conditions, they occur around 20–25% of the time. This suggests that the chaotic diffusion, which characterises this transient regime, is not very sensitive to the initial
<figure><img src="image_11.png" /><figcaption>Fig. 9. Histograms for the final distribution of the inclinations of Ariel (top) and Umbriel (bottom) for different sets of initial conditions. We also plotted the best fit Rice distribution (Eq. (95)) to each histogram (red curve) and the corresponding mean inclination, $\langle I \rangle$, and standard deviation, $\Delta I$. We used the sets #9 (left), #17 (centre), and #19 (right) from Table 2.</figcaption></figure>

inclinations. Finally, as a result of the previous two outcomes, as we increase the initial inclinations of both satellites, the number of simulations that simply skip the resonance also increases.

For those systems that skip or escape the 5/3 MMR, one can ask if the final inclinations are in agreement with the present observations. However, there is not an easy answer because of the chaotic diffusion. Indeed, the final inclinations of Ariel and Umbriel never end exactly with the same values, but they rather follow some kind of statistical distribution. Considering only the trajectories that quickly evade or skip the resonance ($P_e$ and $P_s$), we can build a histogram to better understand how they are distributed for each set of simulations #9 – #25 (i.e. we only consider simulations with less than 95% of capture probability).

In Fig. 9, we show three examples of histograms for the sets #9 ($I_1 = 0.15^\circ$, $I_2 = 0.05^\circ$), #17 ($I_1 = 0.05^\circ$, $I_2 = 0.15^\circ$), and #19 ($I_1 = 0.15^\circ$, $I_2 = 0.15^\circ$). The final inclinations of Ariel and Umbriel are distributed in classes with a size of $0.015^\circ$ and the number of events in each class is normalised by the total number of trajectories that quickly evaded or skipped the resonance. We observe that for Ariel, the final inclinations are more or less evenly distributed between $0.01^\circ$ and $0.25^\circ$, while for Umbriel the final inclinations closely pile up around some mean value.

To better analyse the results in a systematic way, and since $I_k \propto |y_k|$ (Eq. (33)), we fitted a Rice distribution to each data set (Rice 1945), which describes the modulus of a random walk variable in two dimensions. This function has the form

$$f(I) = \frac{I}{\Delta I^2} B_0 \left( \frac{I \langle I \rangle}{\Delta I^2} \right) \exp \left( - \frac{I^2 + \langle I \rangle^2}{2\Delta I^2} \right), \quad (95)$$

where $\langle I \rangle$ is the mean inclination, $\Delta I$ is the standard deviation, and $B_0(x)$ is the modified Bessel function, given by

$$B_0(x) = \sum_{n=0}^{\infty} \frac{x^{2n}}{n!^2}. \quad (96)$$

In Fig. 9, for each histogram we also show the curve of the Rice distribution corresponding to the best fit parameters $\langle I \rangle$ and $\Delta I$. We verified there is a reasonably good agreement between the histogram and the derived distribution, in particular for the final inclination of Umbriel. We hence adopted the best fit parameters obtained in this way to characterise each data set.

The statistical results obtained are listed in Table 2. From a detailed analysis, for the trajectories that are not trapped in resonance, we observe the following: for Ariel, regardless of the its initial inclination value, the final inclinations are always more or less uniformly distributed<fnref n="2" /> in the interval $I_1 \in [0.01^\circ, 0.25^\circ]$; and for Umbriel, the final inclinations are gathered around the mean value with a standard deviation $\Delta I_2 \approx 0.02^\circ$, and the mean value for the final inclination of Umbriel increases with the initial inclinations of both satellites. A more subtle analysis additionally shows that the average final inclination of Umbriel is always lower than its initial value and it also depends on the initial inclination of Ariel. It is approximately given by (see Fig. 10)

$$\langle I_2 \rangle - I_2 \approx -0.06^\circ + 0.22 I_1, \quad (97)$$

<footnote n="2">² We note that by construction, the Rice distribution must be zero when the inclination is zero and it has to peak around the mean value (Eq. (95)). As the final inclinations of Ariel are more or less uniformly distributed, a step function would provide a better adjustment. Nevertheless, we kept the Rice distribution for Ariel for simplicity since it is able to correctly capture the mean value and the dispersion around it.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Variation in the average inclination of Umbriel after crossing the 5/3 MMR as a function of the initial inclination of Ariel. The points are taken from the simulations #11 to #25 in Table 2, and the line gives the adjustment of a linear regression.</figcaption></figure>

that is, to say the lower the initial inclination of Ariel is, the larger the decrease observed in the initial inclination of Umbriel.

The present average inclinations of Ariel and Umbriel are $I_1 \approx 0.02^\circ$ and $I_2 \approx 0.08^\circ$, respectively (Table 1). After crossing the 5/3 MMR, the inclination of Ariel is always approximately uniformly distributed and not very sensitive to the initial conditions. The presently observed value is thus compatible with near zero or moderate initial inclinations of both satellites ($I_k < 0.2^\circ$), but no additional constraints can be derived. However, the final inclination of Umbriel strongly depends on the initial inclination of both satellites (Eq. (97)). We hence conclude that the best configuration that reproduces the present data was obtained for initial $I_1 \approx 0.20^\circ$ and $I_2 \approx 0.10^\circ$ (#15, Table 2) or initial $I_1 \le 0.05^\circ$ and $I_2 \approx 0.15^\circ$ (#16 and #17, Table 2).

The satellites of Uranus were likely formed in a circumplanetary disk (e.g. Pollack et al. 1991; Szulágyi et al. 2018; Ishizawa et al. 2019; Inderbitzi et al. 2020), and so the initial inclinations should have been extremely small. Therefore, initial inclinations of $0.15^\circ$ or higher are very difficult to explain. However, Umbriel was most likely previously involved in a 3/1 MMR with Miranda (Tittemore & Wisdom 1989, 1990), which can excite the inclination of Umbriel to $0.15^\circ$ prior to the encounter with the 5/3 MMR. We hence conclude that the most likely scenario for the initial inclinations of Ariel and Umbriel is $I_1 \le 0.05^\circ$ and $I_2 \approx 0.15^\circ$ (#16 and #17, Table 2).

## 5. Conclusion

Ariel and Umbriel have almost certainly passed through the 5/3 MMR in the past owing to the tidal evolution of their orbits. However, the exact mechanism that allows the system to evade capture in this resonance is a mystery. For coplanar orbits (zero inclinations), the eccentricity of at least one satellite must be close to $0.01$ at the time, which is unlikely because tides are expected to quickly damp the eccentricities to nearly zero (Tittemore & Wisdom 1988). For non-coplanar orbits, the inclinations appear to grow to high values, which is unlikely because tides are very inefficient to damp the inclinations to the presently observed near zero values (Cuk et al. 2020).

To address this question, in this paper we revisited the 5/3 MMR crossing problem only focussing on the inclination. We adopted a secular resonant two-satellite circular model with low inclinations, using a Hamiltonian approach similar to Tittemore & Wisdom (1989) for the study of the 3/1 resonance between Miranda and Umbriel. However, in our model

we included the spin evolution of Uranus, we used complex Cartesian coordinates, and we adopted the total angular momentum of the system as a canonical variable, which is conserved and naturally removes one degree of freedom from the problem. We thus only needed to perform one average over a fast angle, instead of the widely used average over two fast angles (e.g. Tittemore & Wisdom 1988, 1989; Michtchenko & Ferraz-Mello 2001; Alves et al. 2016). We also implemented a Hamiltonian extension for tides based on the constant time-lag model, which provides the exact tidal evolution for all variables in the problem.

Our model is valid for any second order MMR, and thus not restricted to the 5/3 MMR between Ariel and Umbriel. For instance, it can be directly applied to the passage through the 3/1 MMR between Miranda and Umbriel. More generally, it can also be used to study a large number of planetary systems near second order resonances (in the circular approximation), such as the 3/1 MMR for HD 60532 (Laskar & Correia 2009), the 5/3 MMR for HD 33844 (Wittenmyer et al. 2016), the 7/5 MMR for HD 41248 (Jenkins et al. 2013), or the 9/7 MMR for Kepler-29 (Vissapragada et al. 2020).

Applying our model to Ariel and Umbriel, we have shown that, prior to the 5/3 MMR encounter, the system is in circulation around an equilibrium point at the origin of the coordinates ($y_1 = 0, y_2 = 0$). As the system approaches the resonance, this equilibrium becomes unstable, while two other stable symmetrical equilibria appear, corresponding to libration in resonance. Shortly after, the equilibrium at zero becomes stable again, corresponding to a new circulation region (see Fig. 2). For initial very low inclination values (low energy), the system is thus forced to follow one of the resonant equilibria (capture). However, for moderate inclinations (higher energy), the system may encounter the resonance when all the stable equilibrium possibilities are already available and directly go to the circulation region (skip). Alternatively, the system enters into a chaotic regime and can subsequently evade into the circulation (escape) or the libration regions (capture). The chaotic nature of the system as a function of the energy is clearly portrayed in a sequence of Poincaré surface sections obtained with the modified Hénon method (Fig. 3).

The crossing of the 5/3 MMR is a stochastic process, and so we performed a large number of numerical simulations covering many different combinations for the initial inclinations of Ariel and Umbriel. The results show that the initial inclination of Umbriel must have been higher than about $0.1^\circ$ to avoid a permanent capture. Moreover, in order to conciliate the output of the simulations with the presently observed system (Table 1), we find that the optimal inclinations for the satellites prior to the resonance encounter are $I_1 \le 0.05^\circ$ for Ariel and $I_2 \approx 0.15^\circ$ for Umbriel. In this configuration, about 60% of the simulations avoid capture in resonance and the final inclination distribution of Umbriel clusters around the present mean value of $0.08^\circ$. The final inclination distribution of Ariel spreads between $0.01^\circ$ and $0.25^\circ$ with a nearly equal probability, which includes the present mean value of $0.02^\circ$. The inclination damping timescales of ~180 Gyr for Ariel and ~1500 Gyr for Umbriel additionally suggest that the inclination values obtained just after crossing the resonance have nearly remained unchanged up to the present.

The satellites of Uranus were presumably formed in a circumplanetary disk, and so the primordial eccentricities and inclinations should have been extremely small. Before the encounter with the 5/3 MMR, the Uranian satellites may have crossed other MMRs (e.g. Peale 1988; Cuk et al. 2020). In general, these resonances excite the eccentricities and the inclinations of the
bodies involved in the commensurability. Tides are very efficient at damping the eccentricities, but not the inclinations. Therefore, while the remnant eccentricities are quickly eroded, the inclinations are fossilised until the next resonant encounter. For instance, the present high inclination of Miranda (~4.3°) is usually explained after a passage through the 3/1 MMR between Miranda and Umbriel that occurred several billion years ago. Interestingly, this resonance also involves Umbriel, whose inclination can therefore also be excited, though to a much lower value than that of Miranda (Tittemore & Wisdom 1989, 1990). As a consequence, prior to the encounter with the 5/3 MMR, it is reasonable to assume a near zero inclination for Ariel and an inclination around 0.15° for Umbriel.

The results obtained in this paper are very compelling, but they need to be taken with caution because our model is limited to circular orbits and two satellites. In a non-circular model, the eccentricity terms introduce other resonant angles that provide additional libration and chaotic regions to the problem (Tittemore & Wisdom 1988). Indeed, for the 3/1 MMR between Miranda and Umbriel, it was shown that the coupling between the eccentricity and inclination resonances may cause significant variations in the eccentricity evolution of Miranda (Tittemore & Wisdom 1990). Moreover, the presence of the remaining three large satellites or Uranus can also introduce three-body resonances that may further excite the eccentricities and the inclinations (Ćuk et al. 2020). Therefore, in order to fully understand how exactly the system evaded the 5/3 MMR between Ariel and Umbriel and subsequently settled into the present state, future work should take into account the effect of the eccentricities and from all five major satellites in the Uranian system.

*Acknowledgements.* This work was supported by grant SFRH/BD/143371/2019, and by projects CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), 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 the resources to perform the numerical simulations.

## References

* Alexander, M. E. 1973, Ap&SS, 23, 459

* Alves, A. J., Michtchenko, T. A., & Tadeu dos Santos, M. 2016, Celest. Mech. Dyn. Astron., 124, 311

* Cartwright, R. J., Beddingfield, C. B., Nordheim, T. A., et al. 2020, ApJ, 898, L22

* Correia, A. C. M. 2009, ApJ, 704, L1

* Correia, A. C. M., & Rodríguez, A. 2013, ApJ, 767, 128

* Correia, A. C. M., & Valente, E. F. S. 2022, Celest. Mech. Dyn. Astron., 134, 24

* Ćuk, M., El Moutamid, M., & Tiscareno, M. S. 2020, Plan. Sci. J., 1, 22

* Darwin, G. H. 1880, Philos. Trans. R. Soc. London, 171, 713

* Delisle, J.-B., Laskar, J., Correia, A. C. M., & Boué, G. 2012, A&A, 546, A71

* Dermott, S. F., & Nicholson, P. D. 1986, Nature, 319, 115

* Dermott, S. F., Malhotra, R., & Murray, C. D. 1988, Icarus, 76, 295

* Gavrilov, S. V., & Zharkov, V. N. 1977, Icarus, 32, 443

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

* Henon, M., & Heiles, C. 1964, AJ, 69, 73

* Henrard, J. 1982, Celest. Mech., 27, 3

* Henrard, J., & Lemaître, A. 1983, Celest. Mech., 30, 197

* Ida, S., Ueta, S., Sasaki, T., & Ishizawa, Y. 2020, Nat. Astron., 4, 880

* Inderbitzi, C., Szulágyi, J., Cilibrasi, M., & Mayer, L. 2020, MNRAS, 499, 1023

* Ishizawa, Y., Sasaki, T., & Hosono, N. 2019, ApJ, 885, 132

* Jacobson, R. A. 2014, AJ, 148, 76

* Jenkins, J. S., Tuomi, M., Brasser, R., Ivanyuk, O., & Murgas, F. 2013, ApJ, 771, 41

* Kaula, W. M. 1964, Rev. Geophys., 2, 661

* Laskar, J. 1986, A&A, 166, 349

* Laskar, J., & Correia, A. C. M. 2009, A&A, 496, L5

* Malhotra, R., & Dermott, S. F. 1990, Icarus, 85, 444

* Michtchenko, T. A., & Ferraz-Mello, S. 2001, Icarus, 149, 357

* Mignard, F. 1979, Moon Planets, 20, 301

* Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge, UK: Cambridge University Press)

* Palaniyandi, P. 2009, Chaos Solitons Fractals, 39, 1877

* Peale, S. J. 1988, Icarus, 74, 153

* Pollack, J. B., Lunine, J. I., & Tittemore, W. C. 1991, in Uranus, eds. J. T. Bergstrahl, E. D. Miner, & M. S. Matthews (Tucson: University of Arizona Press), 469

* Rice, S. O. 1945, Bell Syst. Tech. J., 24, 46

* Rufu, R., & Canup, R. M. 2022, ApJ, 928, 123

* Singer, S. F. 1968, Geophys. J. R. Astron. Soc., 15, 205

* Smart, W. M. 1953, Celestial Mechanics (London, New York: Longmans, Green)

* Smith, B. A., Soderblom, L. A., Beebe, R., et al. 1986, Science, 233, 43

* Squyres, S. W., Reynolds, R. T., & Lissauer, J. J. 1985, Icarus, 61, 218

* Szulágyi, J., Cilibrasi, M., & Mayer, L. 2018, ApJ, 868, L13

* Thomas, P. C. 1988, Icarus, 73, 427

* Tittemore, W. C. 1990, Icarus, 87, 110

* Tittemore, W. C., & Wisdom, J. 1988, Icarus, 74, 172

* Tittemore, W. C., & Wisdom, J. 1989, Icarus, 78, 63

* Tittemore, W. C., & Wisdom, J. 1990, Icarus, 85, 394

* Verheylewegen, E., Noyelles, B., & Lemaître, A. 2013, MNRAS, 435, 1776

* Vissapragada, S., Jontof-Hutter, D., Shporer, A., et al. 2020, AJ, 159, 108

* Wittenmyer, R. A., Johnson, J. A., Butler, R. P., et al. 2016, ApJ, 818, 35

* Yoder, C. F. 1979, Celest. Mech., 19, 3

* Zahnle, K., Schenk, P., Levison, H., & Dones, L. 2003, Icarus, 163, 263
# Appendix A: Conservative Hamiltonian coefficients

We note that

$$ \mathcal{K}_a = -\frac{p}{2} \frac{\mu_1^2 \beta_1^3}{\Gamma_1^3} + \left(1 + \frac{p}{2}\right) \frac{\mu_2^2 \beta_2^3}{\Gamma_2^3} \qquad (A.1) $$

$$ \mathcal{K}_b = -\frac{3}{2} \left( \frac{p^2}{4} \frac{\mu_1^2 \beta_1^3}{\Gamma_1^4} + \left(1 + \frac{p}{2}\right)^2 \frac{\mu_2^2 \beta_2^3}{\Gamma_2^4} \right) \qquad (A.2) $$

$$ O_a = \frac{3}{2} J_2 R^2 \left( -(p-1) \frac{\mu_1^4 \beta_1^7}{\Gamma_1^7} + (p+2) \frac{\mu_2^4 \beta_2^7}{\Gamma_2^7} \right) \qquad (A.3) $$

$$ O_b = \frac{3}{2} J_2 R^2 \left( -p \frac{\mu_1^4 \beta_1^7}{\Gamma_1^7} + (p+3) \frac{\mu_2^4 \beta_2^7}{\Gamma_2^7} \right) \qquad (A.4) $$

$$ S_a = \frac{G \mu_2 m_1 m_2 \beta_2^2}{\Gamma_2^3} \left(1 + \frac{p}{2} + \frac{\mu_2 \beta_2^2}{\mu_1 \beta_1^2} \left[ \left(1 + \frac{p}{2}\right) \frac{\Gamma_1^2}{\Gamma_2^2} + \frac{p}{2} \frac{\Gamma_1}{\Gamma_2} \right] \frac{\partial}{\partial \alpha} \right) b_{\frac{p}{2}}^{(0)}(\alpha_0) \qquad (A.5) $$

$$ S_b = \frac{G \mu_2 m_1 m_2 \beta_2^2}{4 \Gamma_1 \Gamma_2^2} \alpha_0 b_{\frac{p}{2}}^{(1)}(\alpha_0) \qquad (A.6) $$

$$ S_c = \frac{G \mu_2 m_1 m_2 \beta_2^2}{4 \Gamma_2^3} \alpha_0 b_{\frac{p}{2}}^{(1)}(\alpha_0) \qquad (A.7) $$

$$ S_d = - \frac{G \mu_2 m_1 m_2 \beta_2^2}{2 \Gamma_2^2 \sqrt{\Gamma_1 \Gamma_2}} \alpha_0 b_{\frac{p}{2}}^{(1)}(\alpha_0) \qquad (A.8) $$

$$ \mathcal{R}_a = - \frac{G \mu_2 m_1 m_2 \beta_2^2}{4 \Gamma_1 \Gamma_2^2} \alpha_0 b_{\frac{p}{2}}^{(p+1)}(\alpha_0) \qquad (A.9) $$

$$ \mathcal{R}_b = - \frac{G \mu_2 m_1 m_2 \beta_2^2}{4 \Gamma_2^3} \alpha_0 b_{\frac{p}{2}}^{(p+1)}(\alpha_0) \qquad (A.10) $$

$$ \mathcal{R}_c = \frac{G \mu_2 m_1 m_2 \beta_2^2}{2 \Gamma_2^2 \sqrt{\Gamma_1 \Gamma_2}} \alpha_0 b_{\frac{p}{2}}^{(p+1)}(\alpha_0) \qquad (A.11) $$

where $\alpha_0 \approx 0.7114$ corresponds to $\alpha = a_1/a_2$ at the nominal resonance (Eq. (77)).