© The Authors 2023

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; Cuk 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; Cuk et al. 2020).

The latest low-order commensurability to have occurred was the 5/3 MMR between Ariel and Umbriel (e.g. Peale 1988; Cuk 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.

Ćuk 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₀ (Uranus) surrounded by two point-mass bodies m₁, m₂ ≪ m₀ (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=-{\displaystyle \sum _{k=1}^2\frac{𝒢m_0{m}_k}{r_k}\left[1+J_2{\left(\frac{R}{r_k}\right)}^2{P}_2\left({\widehat{r}}_k\cdot s\right)\right]-\frac{𝒢m_1{m}_2}{\left|r_2-r_1\right|},}$$(1)

where 𝒢 is the gravitational constant; J₂, R, and s are the second order gravity field, the radius, and the spin unit vector of the central body, respectively; r_k is the position vector of m_k with respect to the centre of mass of m₀ (planetocentric coordinates); r_k = |r_k| is the norm; ř_k = r_k/r_k is the unit vector; and P₂(x) = (3x² − 1)/2 is the Legendre polynomial of degree two. We neglected terms in (R/r_k)³ (quadrupolar approximation for the oblateness). The Hamiltonian of the problem, ℋ, 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₀, 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) $$\mathscr{H}={\mathscr{H}}_K+{\mathscr{H}}_O+{\mathscr{H}}_S+{\mathscr{H}}_R+{\mathscr{H}}_F+\frac{{\text{Θ}}^2}{2C},$$(2)

where $${\mathscr{H}}_k=-{\displaystyle \sum _{k=1}^2\frac{𝒢m_0{m}_k}{2a_k}}$$(3)

is the Keplerian part and a_k are the semi-major axes, $${\mathscr{H}}_O=-{\displaystyle \sum _{k=1}^2\frac{𝒢m_0{m}_k}{2a_k}{J}_2{\left(\frac{R}{a_k}\right)}^2\left(1-3I_k^2{\sin }^2\left({\lambda }_k-{\text{Ω}}_k\right)\right)}$$(4)

is the contribution from the oblateness of the central body, λ_k are the mean longitudes, Ω_k are the longitudes of the nodes, $$\begin{array}{l}{\mathscr{H}}_S=-\frac{𝒢m_1{m}_2}{8a_2}[4b_{\frac{1}{2}}^{\left(0\right)}\left(\alpha \right)-\alpha b_{\frac{3}{2}}^{\left(1\right)}\left(\alpha \right)(I_1^2+I_2^2\hfill \\ -2I_1{I}_2\cos \left({\text{Ω}}_2-{\Omega }_1\right)]\hfill \end{array}$$(5)

is the secular part, $b_s^{\left(j\right)}$ are Laplace coefficients, α = a₁ /a₂, and $$\begin{array}{l}{\mathscr{H}}_R=-\frac{𝒢m_1{m}_2}{8a_2}[I_1^2\cos \left(\left(p+2\right){\lambda }_2-p{\lambda }_1-2{\text{Ω}}_1\right)\hfill \\ -2I_1{I}_2\cos \left(\left(p+2\right){\lambda }_2-p{\lambda }_1-{\text{Ω}}_2-{\text{Ω}}_1\right)\hfill \\ +I_2^2\cos \left(\left(p+2\right){\lambda }_2-p{\lambda }_1-2{\text{Ω}}_2\right)]\alpha b_{\frac{3}{2}}^{\left(p+1\right)}\left(\alpha \right)\hfill \end{array}$$(6)

is the contribution from the second order resonant terms (q = 2). The ℋ_S and ℋ_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 ℋ_F corresponds to the remaining terms of the disturbing function that depend on other combinations of the angles λ₁, λ₂, Ω₁, and Ω₂ that do not appear in the expressions of ℋ_S or ℋ_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 $$\text{Θ=}C\omega ,$$(7)

which corresponds to the rotational angular momentum, ω = θ is the angular velocity, and θ 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 (Θ, θ) and Poincaré variables for the orbits (Λ_k, λ_k; Φ_k, −Ω_k) with $${\Lambda }_k={\beta }_k\sqrt{{\mu }_k{a}_k},$$(8) $${\text{Φ}}_k={\text{Λ}}_k\left(1-\cos I_k\right)\approx \frac{1}{2}{\Lambda }_k{I}_k^2,$$(9)

whereß_k = m₀m_k/(m₀ + m_k) and µ_k = 𝒢(m₀ + 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,$$(10)

which is present in all terms of the resonant Hamiltonian (Eq. (6)). Each term corresponds to a resonant combination: $${\phi }_1=\sigma -{\text{Ω}}_1,{\phi }_1=\sigma -{\text{Ω}}_2,{\phi }_3=\frac{1}{2}\left({\phi }_1+{\phi }_1\right).$$(11)

We further introduce the angles $$\gamma =\frac{p}{2}\left({\lambda }_1-{\lambda }_2\right)={\lambda }_2-\sigma \text{\hspace{1em}and\hspace{1em}}\vartheta \text{=}\theta -\sigma ,$$(12)

such that $$\left[\begin{array}{c}\sigma \\ \gamma \\ {\phi }_1\\ {\phi }_1\\ \vartheta \end{array}\right]\equiv 𝒮\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ -{\text{Ω}}_1\\ -{\text{Ω}}_2\\ \theta \end{array}\right],$$(13)

with $$𝒮=\left[\begin{array}{ccccc}-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{array}\right],$$(14)

which gives for the conjugated actions (e.g., Goldstein 1950) $$\left[\begin{array}{c}\sum \\ \text{Γ}\\ {\widetilde{\text{Φ}}}_1\\ {\widetilde{\text{Φ}}}_2\\ \tilde{\Theta }\end{array}\right]={\left({𝒮}^{-1}\right)}^T\left[\begin{array}{c}{\text{Λ}}_1\\ {\text{Λ}}_2\\ {\text{Φ}}_1\\ {\text{Φ}}_2\\ \Theta \end{array}\right].$$(15)

The new set of canonical variables that uses the resonant angles is then given by $$\begin{array}{ll}\text{\hspace{0.17em}\hspace{0.17em}Σ}={\text{Λ}}_1+{\text{Λ}}_2-{\text{Φ}}_1-{\text{Φ}}_2+\text{Θ,}\hfill & \sigma =\left(1+\frac{p}{2}\right){\lambda }_2-\frac{p}{2}{\lambda }_1\hfill \\ \text{\hspace{0.17em}\hspace{0.17em}Γ=}\left(1+\frac{2}{p}\right){\text{Λ}}_1+{\text{Λ}}_2,\hfill & \gamma ={\lambda }_2-\sigma \hfill \\ {\widetilde{\text{Φ}}}_1={\text{Φ}}_1,\hfill & {\phi }_1=\sigma -{\text{Ω}}_1\hfill \\ {\widetilde{\text{Φ}}}_2={\text{Φ}}_2,\hfill & {\phi }_2=\sigma -{\text{Ω}}_2\hfill \\ \widetilde{\text{Θ}}\text{=Θ,}\hfill & \vartheta =\theta -\sigma .\hfill \end{array}$$(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{{\text{Λ}}_k^2}{{\beta }_k^2{\mu }_k},$$(17) $$I_k\approx \sqrt{\frac{2{\text{Φ}}_k}{{\text{Λ}}_k}}.$$(18)

It is important to note, however, that the Λ_k are no longer actions of the resonant variables, and they must be obtained as (Eq. (16)) $${\text{Λ}}_1={\text{Γ}}_1-\frac{p}{2}\left({\text{Φ}}_1+{\text{Φ}}_2\right),$$(19) $${\text{Λ}}_2={\text{Γ}}_2+\left(1+\frac{p}{2}\right)\left({\text{Φ}}_1+{\text{Φ}}_2\right),$$(20)

with $${\text{Γ}}_1=\frac{p}{2}\text{Γ}\left(1-\text{Δ}\right),$$(21) $${\text{Γ}}_2=\frac{p}{2}\text{Γ}\left(1-\left(1+\frac{2}{p}\right)\text{Δ}\right),$$(22)

and $$\text{Δ}=\left(\text{Σ}-\text{Θ}\right)/\text{Γ}.$$(23)

In the approximation of low inclinations, Φ_k ≪ Λ_k (Eq. (18)), and so we also have Φ_k ≪ Γ_k, allowing us to write $${\text{Λ}}_1^{\alpha }\approx {\text{Γ}}_1^{\alpha }\left[1-\alpha \frac{p}{2}\frac{{\text{Φ}}_1+{\text{Φ}}_2}{{\text{Γ}}_1}\right],$$(24) $${\text{Λ}}_2^{\alpha }\approx {\text{Γ}}_2^{\alpha }\left[1-\alpha \left(1+\frac{p}{2}\right)\frac{{\text{Φ}}_1+{\text{Φ}}_2}{{\text{Γ}}_2}\right],$$(25)

and $$I_k\approx \sqrt{\frac{2{\text{Φ}}_k}{{\text{Γ}}_k}}.$$(26)

The angle θ does not appear in the expression of the Hamil-tonian (Eq. (2)), and so ϑ does not appear either (Eq. (16)). As a consequence, the conjugated action, Θ, is a constant of motion. According to the d’Alembert rule (conservation of the angular momentum), the remaining angles must be combined as $$\begin{array}{cc}\cos \left(k_1{\lambda }_1+k_2{\lambda }_2+k_3{\text{Ω}}_1+k_4{\text{Ω}}_2\right),& k_i\in \mathbb{Z},\end{array}$$(27)

where k₁ + k₂ + k₃ + k₄ = 0, or, using the resonant angles (Eq. (16)), $$\cos \left(\left(k_1+\frac{2}{p}{k}_1+k_2\right)\gamma -k_3{\phi }_1+\left(k_1+k_2+k_3\right){\phi }_2\right).$$(28)

Thus, the angle σ also does not appear in the Hamiltonian, and its conjugated action, Σ, is a constant of motion. We note that, $$\begin{array}{l}\sum =\left({\text{Λ}}_1-{\text{Φ}}_1\right)+\left({\text{Λ}}_1-{\text{Φ}}_2\right)+\text{Θ}\hfill \\ ={\beta }_1\sqrt{{\mu }_1{a}_1}\cos I_1{\beta }_2\sqrt{{\mu }_2{a}_2}\cos I_{2}C\omega \hfill \end{array}$$(29)

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 φ₁ and φ₂ present a slower variation than the angle γ, that is,${\dot{\phi }}_1,{\dot{\phi }}_2\ll \dot{\gamma }$. Therefore, we can construct the resonant secular normal form of the Hamiltonian (to the first order in m_k/m₀) by averaging over γ: $$\overline{\mathscr{H}}={\langle \mathscr{H}\rangle }_{\gamma }=\frac{1}{2p\pi }{\displaystyle {\int }_0^{2p\pi }\mathscr{H}d\gamma .}$$(30)

As a result, 〈ℋ_F〉_γ = 0, and since γ no longer appears in the expression of the averaged Hamiltonian, the conjugated variable, Γ (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, (Φ₁, φ₁) and (Φ₂, φ₂), and three parameters, Σ, Γ, and Θ. The auxiliary quantities Γ₁ (Eq. (21)) and Γ2 (Eq. (22)) are also constant.

The resonant secular Hamiltonian then finally reads as follows: $$\begin{array}{l}\overline{\mathscr{H}}=\left({\mathcal{K}}_a+\mathcal{S}\right)\left({\text{Φ}}_1+{\text{Φ}}_2\right)+{\mathcal{K}}_b{\left({\text{Φ}}_1+{\text{Φ}}_2\right)}^2\hfill \\ +\left({\mathcal{O}}_a+{\mathcal{S}}_b\right){\text{Φ}}_1+\left({\mathcal{O}}_b+{\mathcal{S}}_c\right){\text{Φ}}_2\hfill \\ +{\mathcal{S}}_d\sqrt{{\text{Φ}}_1}\sqrt{{\text{Φ}}_2}\cos \left({\phi }_1-{\phi }_2\right)\hfill \\ +{ℛ}_a{\text{Φ}}_1\cos \left(2{\phi }_1\right)+{ℛ}_b{\text{Φ}}_2\cos \left(2{\phi }_2\right)\hfill \\ +{ℛ}_c\sqrt{{\text{Φ}}_1}\sqrt{{\text{Φ}}_2}\cos \left({\phi }_1+{\phi }_2\right),\hfill \end{array}$$(31)

where K stands for the Keplerian coefficients (Eq. (3)), O for the oblateness coefficients (Eq. (4)), S for secular the coefficients (Eq. (5)), and R for the resonant coefficients (Eq. (6)). The Keplerian part needs to be expanded to the second order in Φk (fourth order in the inclinations), because 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 (Φ_k, φ_k) may experience some singularities when Φ_k = 0 (because of the terms in 𝒮_d and R_c). To remove this problem, we can perform a second canonical change of variables to rectangular coordinates (Φ_k, φ_k) → (y_k, iȲ_k), where $$y_k=\sqrt{{\text{Φ}}_k{\text{e}}^{\text{i}{\phi }_k}},$$(32)

and Ȳ_k is the complex conjugate of y_k. From Eq. (26), we have $$y_k\approx I_k\sqrt{\frac{{\text{Γ}}_k}{2}{\text{e}}^{\text{i}{\phi }_k}},$$(33)

and so these variables are proportional to the inclinations. The resonant secular Hamiltonian (Eq. (31)) now reads as follows: $$\begin{array}{c}\overline{\mathscr{H}}=\left({\mathcal{K}}_a+{\mathcal{S}}_a\right)\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)+{\mathcal{K}}_b{\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)}^2\\ +\left({\mathcal{O}}_a+{\mathcal{S}}_b\right)y_1{\overline{y}}_1+\left({\mathcal{O}}_b+{\mathcal{S}}_c\right)y_2{\overline{y}}_2+\frac{{\mathcal{S}}_d}{2}\left(y_1{\overline{y}}_2+{\overline{y}}_1{y}_2\right)\\ +\frac{{ℛ}_a}{2}\left(y_1^2+{\overline{y}}_1^2\right)+\frac{{ℛ}_b}{2}\left(y_2^2+{\overline{y}}_2^2\right)+\frac{{ℛ}_c}{2}\left(y_1{y}_2+{\overline{y}}_1{\overline{y}}_2\right).\end{array}$$(34)

The conservative equations of motion are simply obtained as $$\frac{d{y}_k}{dt}=\text{i}\frac{\partial \overline{\mathscr{H}}}{\partial {\overline{y}}_k},$$(35)

yielding $$\begin{array}{c}{\dot{y}}_1=\text{i}[\left({\mathcal{K}}_a+{\mathcal{S}}_a\right)y_1+2{\mathcal{K}}_b\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)y_1\\ +\left({\mathcal{O}}_a+{\mathcal{S}}_b\right)y_1+\frac{{\mathcal{S}}_d}{2}{y}_2+{ℛ}_a{\overline{y}}_1+\frac{{ℛ}_c}{2}{\overline{y}}_2]\end{array}$$(36)

and $$\begin{array}{c}{\dot{y}}_2=\text{i}[\left({\mathcal{K}}_a+{\mathcal{S}}_a\right)y_2+2{\mathcal{K}}_b\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)y_2\\ +\left({\mathcal{O}}_b+{\mathcal{S}}_c\right)y_2+\frac{{\mathcal{S}}_d}{2}{y}_1+{ℛ}_a{\overline{y}}_2+\frac{{ℛ}_c}{2}{\overline{y}}_1].\end{array}$$(37)

2.2 Dynamical evolution

The values of mk, J2, and R are relatively well determined for the Uranian system (Table 1). Therefore, to compute the coefficients K, O, S, and R appearing in the Hamiltonian (34), we only need to know the values of the parameters Γ1 and Γ2 (see Appendix A), which in turn depend on the parameters Γ and Δ (Eqs.(21) and (22)).

If we neglect the O coefficients (oblateness coefficients), it is possible to eliminate the dependence in Γ, because we get $\overline{\mathscr{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 Γ parameter (e.g. Tittemore & Wisdom 1988, 1989), and so we fixed it at the reference value¹ $$\text{Γ=2}.6684\times {10}^{-12}{M}_{\odot }{\text{au}}^2{\text{yr}}^{-1}.$$(38)

The dynamics of the 5/3 MMR essentially depends on the Δ parameter (Eq. (23)), which measures the proximity to the resonance. Following Delisle et al. (2012), we write $$\delta =\frac{\text{Δ}}{{\text{Δ}}_r}-1,$$(39)

where Δ_r is the value of Δ at the planar (I_k = 0) nominal resonance, that is (Eq. (23) and (29)), $${\text{Δ}}_r=\left({\text{Λ}}_{1,r}+{\text{Λ}}_{2,r}\right)/{\text{Γ}}_r,$$(40)

when n₁/n₂ = 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 Λ_k,r: $$\begin{array}{ccc}{\text{Λ}}_{2,r}=ϵ{\left(\frac{5}{3}\right)}^{1/3}{\text{Λ}}_{1,r},& \text{with}& ϵ=\frac{{\beta }_2{\mu }_2^{2/3}}{{\beta }_1{\mu }_1^{2/3}}\approx \frac{m_2}{m_1}.\end{array}$$(41)

From the expression of Γr (Eq. (16)), we additionally have $${\text{Λ}}_{2,r}={\text{Γ}}_r-\frac{5}{3}{\text{Λ}}_{1,r},$$(42)

which can be combined with expression (41) to give $${\text{Λ}}_{1,r}={\left(\frac{5}{3}+ϵ{\left(\frac{5}{3}\right)}^{1/3}\right)}^{-1}{\text{Γ}}_r$$(43) $${\text{Λ}}_{2,r}={\left(1+{ϵ}^{-1}{\left(\frac{5}{3}\right)}^{2/3}\right)}^{-1}{\text{Γ}}_r,$$(44)

and finally (Eq. (40)) $${\text{Δ}}_r=\left(1+ϵ{\left(\frac{5}{3}\right)}^{1/3}\right){\left(\frac{5}{3}+ϵ{\left(\frac{5}{3}\right)}^{1/3}\right)}^{-1}.$$(45)

2.2.1 Equilibrium points

The equilibrium points correspond to stationary solutions of the Hamiltonian. They can be found by solving $$\begin{array}{ccc}\frac{\partial \overline{\mathscr{H}}}{\partial y_1}=0& \text{and}& \frac{\partial \overline{\mathscr{H}}}{\partial y_2}=0,\end{array}$$(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 + iy_k,i, we get $$y_{2,r}=-2\frac{{\mathcal{K}}_a+2{\mathcal{K}}_b\left({\text{Φ}}_1+{\text{Φ}}_2\right)+{\mathcal{O}}_a+{\mathcal{S}}_a+{\mathcal{S}}_b+{ℛ}_a}{{\mathcal{S}}_d+{ℛ}_c}{y}_{1,r},$$(47) $$y_{1,r}=-2\frac{{\mathcal{K}}_a+2{\mathcal{K}}_b\left({\text{Φ}}_1+{\text{Φ}}_2\right)+{\mathcal{O}}_b+{\mathcal{S}}_a+{\mathcal{S}}_c+{ℛ}_b}{{\mathcal{S}}_d+{ℛ}_c}{y}_{2,r},$$(48) $$y_{2,i}=-2\frac{{\mathcal{K}}_a+2{\mathcal{K}}_b\left({\text{Φ}}_1+{\text{Φ}}_2\right)+{\mathcal{O}}_a+{\mathcal{S}}_a+{\mathcal{S}}_b+{ℛ}_a}{{\mathcal{S}}_d-{ℛ}_c}{y}_{1,i},$$(49) $$y_{1,i}=-2\frac{{\mathcal{K}}_a+2{\mathcal{K}}_b\left({\text{Φ}}_1+{\text{Φ}}_2\right)+{\mathcal{O}}_b+{\mathcal{S}}_a+{\mathcal{S}}_c+{ℛ}_b}{{\mathcal{S}}_d-{ℛ}_c}{y}_{2,i}.$$(50)

Replacing expressions (48) and (50) into expressions (47) and (49), respectively, yields that equilibria arise for $$\begin{array}{ccc}{y}_{1,r}=y_{2,r}=0& \text{or}& y_{1,i}=y_{2,i}=0.\end{array}$$(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)) $$\begin{array}{ccc}{y}_{1,i}=\pm \sqrt{{\text{Φ}}_1}& \text{and}& y_{2,i}=\pm \sqrt{{\text{Φ}}_2},\end{array}$$(52)

which we replaced in expressions (49) and (50) to determine the possible values of Φ1 and Φ2, $${\text{Φ}}_1={\text{Φ}}_2=0$$(53)

or $${\text{Φ}}_1=\frac{\left({ℛ}_c-{\mathcal{S}}_d\right){ℰ}_{\pm }-2\left({\mathcal{K}}_a+{\mathcal{O}}_a+{\mathcal{S}}_a+{\mathcal{S}}_b-{ℛ}_a\right)}{4{\mathcal{K}}_b\left(1+{ℰ}_{\pm }^2\right)},$$(54) $${\text{Φ}}_2={ℰ}_{\pm }^2{\text{Φ}}_1,$$(55)

with $$\begin{array}{l}{ℰ}_{\pm }=\frac{{\mathcal{O}}_a-{\mathcal{O}}_b-{ℛ}_a+{ℛ}_b+{\mathcal{S}}_b-{\mathcal{S}}_c}{{ℛ}_c-{\mathcal{S}}_d}\hfill \\ \pm \frac{\sqrt{{\left({\mathcal{O}}_a-{\mathcal{O}}_b-{ℛ}_a+{ℛ}_b+{\mathcal{S}}_b-{\mathcal{S}}_c\right)}^2+{\left({ℛ}_c-{\mathcal{S}}_d\right)}^2}}{{ℛ}_c-{\mathcal{S}}_d}.\hfill \end{array}$$(56)

The equilibrium point at Φ1 = Φ2 = 0 is always present, although it can be stable or unstable. The remaining equilibria only exist for some δ values. In Fig. 1, we show the evolution of the equilibrium points as a function of δ. 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 δ values far from zero, there is only one equilibrium point at I_k = 0, which is stable (in blue). For δ = 2.068 × 10⁻⁶, 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 δ = −2.540 × 10⁻⁷, which is very close to the resonance nominal value δ = 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 different values of δ (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 visible (Eq. (52)). In Fig. 2, we show the energy levels for three values of δ, which are representative of the three different equilibrium 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₁ sin φ₁, ƒ₂ sin φ₂) with cos φ₁ = cos φ₂ = 0.

For δ = 5 × 10⁻⁶ > 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 δ values), all trajectories are outside the 5/3 MMR.

For δ = 0 (Fig. 2b), the system is at the nominal resonance (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 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 δ = −2 × 10⁻⁶ < 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 δ values, but the central circulation region becomes larger, while the resonant islands become thinner.

Fig. 1 Evolution of the equilibrium points as a function of δ. 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).

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 δ = −2 × 10⁻⁶. We thus adopted this value to draw the surface 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₁ sin φ₁, I₁ cos φ₁) with cos φ₂ = 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, $\overline{\mathscr{H}}<{\mathscr{H}}_0$, the separatrix corresponds to a transition level, $\overline{\mathscr{H}}={\mathscr{H}}_0$, while the largest energies occur in the libration region, $\overline{\mathscr{H}}>{\mathscr{H}}_0$. The inner circulation region is delimited by the levels $0<\overline{\mathscr{H}}<{\mathscr{H}}_0$, where $\overline{\mathscr{H}}=0$ corresponds to the energy of the equilibrium point with y₁ = y₂ = 0 (blue point in Fig. 2c). For this energy range, there are four families, while for the remaining energies only two families exist.

For $\overline{\mathscr{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 ℋ = 0 (Figs. 3b,c), some chaotic regions appear in the transition between the circulation and libration regions. For $0<\overline{\mathscr{H}}<{\mathscr{H}}_0$(family 1), the chaotic regions increase, while the resonant islands shrink (Figs. 3d,e), until they completely disappear for $\overline{\mathscr{H}}>{\mathscr{H}}_0$ (Fig. 3f). For this specific energy range, we also needed to plot family 2. Close to $\overline{\mathscr{H}}=0$, we observe quasi-periodic motion in the inner circulation region, but, as we approach $\overline{\mathscr{H}}={\mathscr{H}}_0$, this area is also completely replaced by a chaotic region (Figs. 3j−1). Finally, for $\overline{\mathscr{H}}={\mathscr{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 y1 and y2 (Eq. (33)). Therefore, the value of the inclinations of Ariel and Umbriel when the system encounters the resonance can trigger completely different behaviours. For $\overline{\mathscr{H}}=0$, the motion is quasi-periodic, either in circulation or libration. Near the separatrix, $\overline{\mathscr{H}}~0$, the motion is mainly chaotic. Finally, for $\overline{\mathscr{H}}\gg {\mathscr{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).

Fig. 2 Energy level curves in the plane(I1 sin φ1,I2 sin φ2) with cos φ1 = cos φ2 = 0, for δ = 5 × 10−6 (top), δ = 0 (middle), and δ = −2 × 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.

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, τ, 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{𝒢m_k^2}{R}{\left(\frac{R}{r_k}\right)}^3{\left(\frac{R}{r{\prime }_k}\right)}^3{P}_2\left({\widehat{r}}_k\cdot \widehat{r}{\prime }_k\right),$$(57)

where k₂ is the elastic second Love number for potential of the body, r_k = r_k (t) is the position of the pertuber at a time t, and ${r^{\prime }}_k=r_k\left(t-\tau \right)$ 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, ${r^{\prime }}_k$, as parameters (Mignard 1979). The tidal Hamiltonian then reads (Eqs. (2) and (57)) as follows: $${\mathscr{H}}_t=\mathscr{H}+U_1+U_2.$$(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{array}{l}{U}_k=-k_2\frac{𝒢m_k^2{R}^5}{4a_k^{3}a{\prime }_k^3}[1+3\cos \left(2\theta -2\theta \prime -2{\lambda }_k+2\lambda {\prime }_k\right)\hfill \\ -\frac{3}{2}{I}_k^2(1-\cos \left(2{\lambda }_k-2{\text{Ω}}_k\right)+\cos \left(2\theta -2\theta \prime -2{\lambda }_k+2\lambda {\prime }_k\right)\hfill \\ -\cos \left(2\theta -2\theta \prime +2\lambda {\prime }_k-2{\text{Ω}}_k\right))\hfill \\ -\frac{3}{2}I{\prime }_k^2(1-\cos \left(2\lambda {\prime }_k-2\text{Ω}{\prime }_k\right)+\cos \left(2\theta -2\theta \prime -2{\lambda }_k+2\lambda {\prime }_k\right)\hfill \\ -\cos \left(2\theta -2\theta \prime -2{\lambda }_k+2\text{Ω}{\prime }_k\right))\hfill \\ +3I_{k}I{\prime }_k(\cos \left(\theta -\theta \prime -2{\lambda }_k+2\lambda {\prime }_k+{\text{Ω}}_k-\text{Ω}{\prime }_k\right)\hfill \\ +\cos \left(\theta -\theta \prime -{\text{Ω}}_k+\text{Ω}{\prime }_k\right)\hfill \\ -\cos \left(\theta -\theta \prime -2{\lambda }_k+{\text{Ω}}_k+\text{Ω}{\prime }_k\right)\hfill \\ -\cos \left(\theta -\theta \prime +2\lambda {\prime }_k-{\text{Ω}}_k-\text{Ω}{\prime }_k\right)].\hfill \end{array}$$(59)

We first performed the canonical change of variables that uses the resonant angles (Eq. (16)), and then changed to the complex Cartesian coordinates (Eq. (32)), to get $$\begin{array}{l}{U}_k=-k_2\frac{𝒢m_k^2{R}^5{\beta }_k^{12}{\mu }_k^6}{4{\text{Γ}}_k^7\text{Γ}{\prime }_k^7}[{\text{Γ}}_k\text{Γ}{\prime }_k-3\text{Γ}{\prime }_k{\overline{y}}_k{y}_k-3{\text{Γ}}_k\overline{y}{\prime }_{k}y{\prime }_k\hfill \\ -6p_k\left(\text{Γ}{\prime }_k{y}_1{\overline{y}}_1+{\text{Γ}}_{k}y{\prime }_1\overline{y}{\prime }_1+\text{Γ}{\prime }_k{y}_2{\overline{y}}_2+{\text{Γ}}_{k}y{\prime }_2\overline{y}{\prime }_2\right)\hfill \\ +3({\text{Γ}}_k\text{Γ}{\prime }_k-\text{Γ}{\prime }_k{\overline{y}}_k{y}_k-{\text{Γ}}_k\overline{y}{\prime }_{k}y{\prime }_k\hfill \\ -6p_k\left(\text{Γ}{\prime }_k{y}_1{\overline{y}}_1+{\text{Γ}}_{k}y{\prime }_1\overline{y}{\prime }_1+\text{Γ}{\prime }_k{y}_2{\overline{y}}_2+{\text{Γ}}_{k}y{\prime }_2\overline{y}{\prime }_2\right))\hfill \\ \times \cos \left(2\left(\vartheta -\vartheta \prime -q_k\left(\gamma -\gamma \prime \right)\right)\right)\hfill \\ +\frac{3}{2}{\text{Γ}}_k\left(\overline{y}{\prime }_k^2+y{\prime }_k^2\right)\cos \left(2q_k\gamma \prime \right)\hfill \\ +\frac{3}{2}\text{Γ}{\prime }_k\left({\overline{y}}_k^2+y_k^2\right)\cos \left(2q_k\gamma \right)\hfill \\ +3\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left({\overline{y}}_{k}y{\prime }_k+\overline{y}{\prime }_k{y}_k\right)\cos \left(\vartheta -\vartheta \prime \right)\hfill \\ -3\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left({\overline{y}}_k\overline{y}{\prime }_k+y_{k}y{\prime }_k\right)\cos \left(\vartheta -\vartheta \prime +2q_k\gamma \prime \right)\hfill \\ -3\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left({\overline{y}}_k\overline{y}{\prime }_k+y_{k}y{\prime }_k\right)\cos \left(\vartheta -\vartheta \prime +2q_k\gamma \right)\hfill \\ +\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left({\overline{y}}_{k}y{\prime }_k+\overline{y}{\prime }_k{y}_k\right)\cos \left(\vartheta -\vartheta \prime -2q_k\left(\gamma -\gamma \prime \right)\right)\hfill \\ +\frac{3}{2}\text{Γ}{\prime }_k\left(\overline{y}{\prime }_k^2+y_k^2\right)\cos \left(2\left(\vartheta -\vartheta \prime +q_k\gamma \prime \right)\right)\hfill \\ +\frac{3}{2}{\text{Γ}}_k\left(\overline{y}{\prime }_k^2+y{\prime }_k^2\right)\cos \left(2\left(\vartheta -\vartheta \prime -q_k\gamma \right)\right)\hfill \\ -\frac{3}{2}\text{iΓ}{\prime }_k\left(\overline{y}{\prime }_k^2-y{\prime }_k^2\right)\sin \left(2q_k\gamma \prime \right)\hfill \\ -\frac{3}{2}\text{iΓ}{\prime }_k\left({\overline{y}}_k^2-y{\prime }_k^2\right)\sin \left(2q_k\gamma \right)\hfill \\ +3\text{i}\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left(\overline{y}{\prime }_k{y}_k-{\overline{y}}_{k}y{\prime }_k\right)\sin \left(\vartheta -\vartheta \prime \right)\hfill \\ +3\text{i}\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left({\overline{y}}_k\overline{y}{\prime }_k-y_{k}y{\prime }_k\right)\sin \left(\vartheta -\vartheta \prime -2q_k\gamma \prime \right)\hfill \\ -3\text{i}\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left({\overline{y}}_k\overline{y}{\prime }_k-y_{k}y{\prime }_k\right)\sin \left(\vartheta -\vartheta \prime -2q_k\gamma \right)\hfill \\ -3\text{i}\sqrt{{\text{Γ}}_k\text{Γ}{\prime }_k}\left(\overline{y}{\prime }_k{y}_k-{\overline{y}}_{k}y{\prime }_k\right)\sin \left(\vartheta -\vartheta \prime +2q_k\left(\gamma -\gamma \prime \right)\right)\hfill \\ -\frac{3}{2}\text{iΓ}{\prime }_k\left({\overline{y}}_k^2-y_k^2\right)\sin \left(2\left(\vartheta -\vartheta \prime +q_k\gamma \prime \right)\right)\hfill \\ +\frac{3}{2}{\text{iΓ}}_k\left(\overline{y}{\prime }_k^2-y{\prime }_k^2\right)\sin \left(2\left(\vartheta -\vartheta \prime -q_k\gamma \right)\right)],\hfill \end{array}$$(60)

with p₁ = −p/2, p₂ = 1 + p/2, q₁ = 1 + 2/p, and q₂ = 1.

We note that σ does not appear in the expression of U_k. Therefore, in the presence of tides the parameter Σ (Eq. (29)) remains conserved. The fast angle γ 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 γ′ is considered as a parameter that can later cancel with γ (see Eq. (64)).

Fig. 3 Poincaré surfaces of section for Ariel in the plane (I1 sin φ1,I1 cos φ1) with cos φ2 = 0 and δ = −2 × 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 ℋ0 = 1.06 × 10−19 M⊙ au2 yr−2.

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=\text{i}{\displaystyle \sum _{j=1}^2\frac{\partial U_j}{\partial {\overline{y}}_k},\dot{\text{Γ}}=-{\displaystyle \sum _{j=1}^2\frac{\partial U_j}{\partial \gamma },\dot{\text{Θ}}=-{\displaystyle \sum _{j=1}^2\frac{\partial U_j}{\partial \vartheta }.}}}$$(61)

In principle, we should also write the equations for γ and $\dot{\gamma }$ and$\dot{\vartheta }$, 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, τ. This model is widely used and provides very simple expressions for the tidal interactions, because it can be made linear (e.g. Mignard 1979): $$\lambda {\prime }_k\approx {\lambda }_k-n_k\tau \text{\hspace{1em}and\hspace{1em}}\theta \prime \approx \theta -\omega \tau .$$(62)

It follows for the remaining primed quantities that $$y{\prime }_k\approx y_k-\text{i}{y}_k\left(p_2{\overline{n}}_2+p_1{\overline{n}}_1\right)\tau ,$$(63) $$\gamma \prime \approx \gamma -p_1\left({\overline{n}}_2-{\overline{n}}_1\right)\tau ,$$(64) $$\vartheta \prime \approx \vartheta +\left(p_2{\overline{n}}_2+p_1{\overline{n}}_1-\omega \right)\tau ,$$(65)

with $${\overline{n}}_k={\beta }_k^3{\mu }_k^2/{\text{Γ}}_k^3.$$(66)

We then substituted expressions (63) to (65) into the equations of motion (61) and averaged over the fast angle γ (as in Eq. (30)) to finally get the secular equations for the tidal evolution $${\dot{y}}_1=-\frac{3}{2}\frac{{𝒟}_1}{{\text{Γ}}_1^{13}}\left(2\text{i}p+{\overline{n}}_1\tau \right)y_1+3\text{i}\left(p+2\right)\frac{{𝒟}_2}{{\text{Γ}}_2^{13}}{y}_1,$$(67) $${\dot{y}}_2=-\frac{3}{2}\frac{{𝒟}_2}{{\text{Γ}}_1^{13}}(-2\text{i}\left(p+2\right)+{\overline{n}}_2\tau )y_2-3\text{i}p\frac{{𝒟}_1}{{\text{Γ}}_2^{13}}{y}_2,$$(68) $$\begin{array}{l}\dot{\text{Γ}}\text{=3}\frac{{𝒟}_1}{{\text{Γ}}_1^{13}}\left(1+2/p\right)[\left({\text{Γ}}_1+6p\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)-y_1{\overline{y}}_1\right)\omega \hfill \\ -\left({\text{Γ}}_1+6p\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right){\overline{n}}_1\right)]\tau \hfill \\ +3\frac{{𝒟}_2}{{\text{Γ}}_1^{13}}[{\text{Γ}}_1-6\left(p+2\right)\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)-y_2{\overline{y}}_2)\omega \hfill \\ -\left({\text{Γ}}_2-6\left(p+2\right)\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)\right){\overline{n}}_2]\tau ,\hfill \end{array}$$(69) $$\begin{array}{l}\dot{\text{Θ}}\text{=3}\frac{{𝒟}_1}{{\text{Γ}}_1^{13}}[\left({\text{Γ}}_1+6p\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)-y_1{\overline{y}}_1\right)\omega \hfill \\ -\left({\text{Γ}}_1+6p\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)-y_2{\overline{y}}_2\right){\overline{n}}_1]\tau \hfill \\ -3\frac{{𝒟}_2}{{\text{Γ}}_1^{13}}[{\text{Γ}}_2-6\left(p+2\right)\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)-y_2{\overline{y}}_2)\omega \hfill \\ -\left({\text{Γ}}_2-6\left(p+2\right)\left(y_1{\overline{y}}_1+y_2{\overline{y}}_2\right)-y_2{\overline{y}}_2\right){\overline{n}}_2]\tau ,\hfill \end{array}$$(70)

where $${𝒟}_k=k_2𝒢m_k^2{\beta }_k^{12}{\mu }_k^6{R}^5.$$(71)

We note that in the expressions of y1 (Eq. (67)) and y2 (Eq. 68)), we have a conservative contribution (imaginary terms) and a dissipative contribution (real terms in τ). 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 Γ (Eq. (69)) and Θ (Eq. (70)). Then, the coefficients 𝒦, O, 𝒮, and ℛ 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 Γ, but also in J₂, according to (e.g. Correia & Rodríguez 2013) $$J_2=k_{\text{f}}\frac{{\omega }^2{R}^3}{3𝒢m_0},$$(72)

where ω = Θ/C (Eq. (7)) and k_f is the fluid second Love number for potential. Using the present rotation rate and J₂ 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,$$(73) $${\dot{I}}_k\approx -\frac{A_k}{2}\frac{\omega }{n_k}{I}_k,$$(74)

with $$A_k=\frac{3𝒢m_k^2{R}^5}{{\beta }_k{a}_k^8}{k}_2\tau .$$(75)

Both Ariel and Umbriel have very low inclinations (Table 1). Combined with the fast rotation rate of Uranus (ω/n_k ≫ 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 (ω/n_k ≈ cte), we get $$\frac{a_2^{7}d{a}_2}{m_2\left(\omega /n_1-1\right)}\approx \frac{a_1^{7}d{a}_1}{m_2\left(\omega /n_2-1\right)},$$(76)

which can be integrated to obtain a₂ as a function of a₁. 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 {\left(5/3\right)}^{2/3}\approx 1.4057$$(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,a_2/R=10.3891,$$(78)

and to compute other related parameters, such as Γ (Eq. (38)).

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 a2 and a1 (Eq. (77)). The blue point gives the best estimation for the semi-major axes at the nominal resonance (Eq. (78)).

3.3.3 Evolution timescale

The orbital evolution timescale depends on τ, which is related to the tidal dissipation (Eq. (62)). A more commonly used dimen-sionless quantity to measure the tidal dissipation is given by the quality factor (e.g. Correia & Valente 2022), $$Q=1/\left(2\omega \tau \right).$$(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 11000 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 τ ≈ 0.62 s. In fact, the exact dissipation rate depends on the product k₂τ (Eq. (75)). As in previous studies, in the present work we adopted k₂ = 0.104 (Gavrilov & Zharkov 1977), which translates into $$\frac{k_2}{Q}=1.3\times {10}^{-5}\iff k_2\tau =0.064\text{\hspace{0.17em}s.}$$(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}{3𝒢m_k^2{R}^5{k}^2},$$(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 φ₁ or φ₂. 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)) $${\mathscr{H}}_k={𝒜}_k{y}_k{\overline{y}}_k+{ℬ}_k{\left(y_k{\overline{y}}_k\right)}^2+\frac{{\mathcal{𝒞}}_k}{2}\left(y_k^2{\overline{y}}_k^2\right),$$(82)

with $${𝒜}_1={𝒦}_a+{𝒪}_a+{𝒮}_a+{𝒮}_{b,}{ℬ}_1={𝒦}_b,{𝒞}_1={ℛ}_a,$$(83) $${𝒜}_1={𝒦}_a+{𝒪}_a+{𝒮}_a+{𝒮}_{c,}{ℬ}_2={𝒦}_b,{𝒞}_2={ℛ}_b,$$(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^C$, and two circulation regions, one outside the libra-tion region and another encircled by the separatrix, with an area ${\mathcal{J}}_k^C$. Consider an initial trajectory with some energy ℋ_k < ℋ₀. 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{𝒥}}_k^R}{{\dot{𝒥}}_k^R+{𝒥}_k^C}.$$(85)