© M. Saillenfest et al. 2022

1 Introduction

The spin axis of Uranus is inclined by about 98° with respect to its orbit normal. This large inclination is particularly puzzling as all other giant planets in the Solar System have obliquity values smaller than 30°. The most widely accepted explanation for this discrepancy is a series of impacts with large planetesimals that probably occurred at the end of the accretion phase (see e.g. Morbidelli et al. 2012; Izidoro et al. 2015; Salmon & Canup 2022). In this picture, the two ice giants Uranus and Neptune would have acquired random spin-axis orientations, whereas runaway gas accretion would have driven the obliquities of Jupiter and Saturn towards 0° (Ward & Hamilton 2004; Hamilton & Ward 2004). Differing impact histories may also explain why Uranus and Neptune have qualitatively different satellite systems and show dissimilar heat fluxes (Reinhardt et al. 2020). However, as argued by Rogoszinski & Hamilton (2020), this scenario suffers from several weaknesses. Uranus and Neptune have strikingly similar masses, radii, and spin rates; they also feature similar atmospheric dynamics (Aurnou et al. 2007), and both have non-dipolar magnetic fields which are likely produced through the same internal convection processes (see e.g. Bailey & Stevenson 2021, for a review). If the final formation stages of Uranus and Neptune had been dominated by giant impacts, one would expect them to present a much larger diversity of properties (see e.g. Chau et al. 2021), perhaps in the spirit of the terrestrial planets of the Solar System. Even though the rotation state of both Uranus and Neptune can be reproduced individually by a finely tuned collision (see e.g. Slattery et al. 1992; Reinhardt et al. 2020, and the insightful discussions of Rufu & Canup 2022), similar spin rates, above all, are not what one would expect as the outcome of random collisions. Instead, the strong similarities between Uranus and Neptune may indicate that their final properties result from smoother processes, and that gas accretion supplied enough angular momentum to dominate their rotation states, similarly to Jupiter and Saturn. This idea is supported by the fact that the spin rates of all four giant planets are a fraction of their critical rotation speeds in a consistent way with formation theories (Machida et al. 2008; Ward & Canup 2010; Lovelace et al. 2011; Batygin 2018; Bryan et al. 2018; Dong et al. 2021; Dittmann 2021). In this unified picture, we would expect the obliquities of Uranus and Neptune to have been small after their formations, and to have evolved to their observed values (about 98° and 30°) at later stages.

Today, the orbits and spin axes of Uranus and Neptune are stable and far from any kind of strong perturbation. Therefore, it seems natural to assume that their tilting occurred early, when the Solar System was still evolving. It is now believed that the last large-scale changes in the Solar System happened during a phase of instability that followed the dissipation of the protoplanetary disc, as the giant planets were migrating in a swarm of planetesimals (Tsiganis et al. 2005). It was initially proposed that this event could be related to a phase of resurging collisions recorded in the lunar craters (called the late heavy bombardment; see Gomes et al. 2005). Numerous refinements were then brought to the instability model in order to properly reproduce the structure of the Solar System as we observe it today. The most recent advances favour a very early migration and instability (within 10 Myr, or even immediately after dispersal of the gas disc), and no relation with the late heavy bombardment (see Clement et al. 2018, Clement et al. 2021; Nesvorný et al. 2018; Morgan et al. 2021).

The instability itself did not exactly place the giant planets on their current orbits. After the instability, the scattering of the remaining planetesimals produced a last grasp of migration, during which Uranus and Neptune still migrated over a few astronomical units. In this much depleted environment, the migration of Uranus and Neptune was slower, in accordance with observations: in order to reproduce the observed Kuiper belt, Neptune probably migrated during this last stage with an e-folding timescale larger than 10 Myr and possibly extending to about 100 Myr at the very end of the migration (see Nesvorný 2015, Nesvorný 2018 for a review). Yet, as noted by Ida et al. (2000), Neptune’s migration could not have been arbitrarily slow; otherwise, mean-motion resonances in the Kuiper belt would be overcrowded. Putting everything together, Uranus and Neptune most likely finished their migration from several tens to a few hundreds of millions of years at most after the dispersal of the gas disc. This timespan puts strong constraints on their tilting scenarios if they are required to happen during the planetary migration.

Apart from collisions, possible tilting mechanisms involve some kind of spin-orbit coupling (Boué & Laskar 2010; Quillen et al. 2018; Millholland & Batygin 2019; Rogoszinski & Hamilton 2020, 2021). This means that at some point, the spin-axis precession motion of the planet enters into resonance with some harmonics stemming from its orbital dynamics. If this happens while the orbits of the planets are still being reshaped, then the resonance properties change over time. Previous authors found that during planetary migration, there exist several ways to displace the resonant equilibrium configuration of the spin axis from a small obliquity ε ≈ 0° to over 90° (and even up to 180° in the case described by Quillen et al. 2018). However, because Uranus and Neptune are today so far from the Sun, their spin-axis precession timescales are extremely large and count in hundreds of millions of years. In order to yield reasonable tilting probabilities, the displacement of the equilibrium must be very slow (i.e. adiabatic) compared to the libration inside the resonance, which itself is slower than the spin-axis precession motion. Hence, for Uranus placed near its current location, authors invariably determined that the timescale for a tilting up to 90° or more greatly exceeds the expected duration of the planetary migration (see e.g. Quillen et al. 2018).

In order to speed up the spin-axis precession of Uranus and catch a strong resonance, we can imagine that Uranus formed much closer to the Sun than it is now (Rogoszinski & Hamilton 2021), or that it once had a big satellite (Boué & Laskar 2010) or a massive disc around it (Harris & Ward 1982; Rogoszinski & Hamilton 2020). Indeed, close-in planets precess faster, and satellites or circumplanetary discs have the ability to magnify the spin-axis precession rate of their host planet (see e.g. Tremaine 1991; Boué & Laskar 2006; Millholland & Batygin 2019). However, authors found that even with these refinements, the timescale required for a purely adiabatic tilting of Uranus hardly fits in the limited amount of time offered by planetary migration (see e.g. Rogoszinski & Hamilton 2021). Moreover, the late planetary migration has certainly not been perfectly smooth, but “grainy” (as a result of the scattering of large planetesimals; see e.g. Nesvorný & Vokrouhlický 2016). The level of smoothness of the planetary migration further complicates the adiabatic tilting process, as a sudden change in the semi-major axis of a given planet may release Uranus out of resonance, in which case the tilting would stop.

Yet, the tilting of Uranus may not have been produced through an adiabatic process: instead of moving the stable equilibrium point, the quickest way to incline a planet via a resonance is to do so over a single large-amplitude libration cycle (a “resonance kick”). In this case, the resonance itself must be extremely wide, which puts new constraints on the orbital dynamics of the planets. In Appendix A, we show that the orbital inclination I of the planet with respect to the invariable plane must verify tan $I\ge 3-2\sqrt{2}$ in order to allow for trajectories going from ε ≈ 0° to over 90° during a single libration inside a secular spin-orbit resonance. Therefore, irrespectively of the secular spin-orbit resonance involved (e.g. s₇ or s₈, as considered by previous authors), a large orbital inclination of more than about 10° is required for the resonant harmonic. This limit quite accurately matches the numerical experiments of Rogoszinski & Hamilton (2020). However, in order to yield a reasonable probability of reproducing Uranus’s current state (instead of a single finely tuned pathway), the inclination needed is actually much larger than 10°. In fact, the most successful tilting scenarios found in previous studies involve, altogether, a massive satellite (or disc), a resonance drift due to planetary migration, and a phase of very large orbital inclination (see Boué & Laskar 2010; Rogoszinski & Hamilton 2020). Unfortunately, as discussed by Vokrouhlický & Nesvorný (2015), large inclination values for Uranus and/or Neptune are an improbable transient feature of Solar System migration models, and none was reported in the extensive simulations of Nesvorný & Morbidelli (2012). Moreover, the mass required for the satellite in order to reach a strong secular spinorbit resonance was found to be implausibly large (e.g. 0.01 Uranus’s mass), and some additional mechanism was needed to remove the satellite after the tilting. Hence, the timescale issue remains a major problem for previous “collisionless” scenarios for the tilting of Uranus.

In the alternative migration scenario recently proposed by Lu & Laughlin (2022), the large-inclination harmonic is provided by the hypothetical Planet 9 (see e.g. Batygin et al. 2019). Even though this scenario does not alleviate the need for a massive circumplanetary disc, the large inclination of Planet 9 allows for dramatic resonance kicks that can easily produce the required level of obliquity excitation for Uranus. However, much uncertainty remains about the migration history (and existence) of Planet 9, and it is not clear whether its large inclination actually predates its scattering away, as assumed by the authors. This scenario also suffers from the same caveats as previous works regarding Uranus’s disc: a massive equatorial disc needs to be maintained and continuously fed during the whole tilting phase, which may not be realistic (see Rogoszinski & Hamilton 2020); besides, the accretion of matters from the disc onto Uranus would slow down the resonance crossing phenomenon, possibly beyond the timescale of planetary migration.

As the spin axis of Neptune is much less inclined, its current orientation is easier to explain. Possible scenarios involve a very early tilt during planetary formation (Martin et al. 2020) and/or a resonance crossing during the late planetary migration (Rogoszinski & Hamilton 2020). If Neptune was surrounded by a massive circumplanetary disc at that time, even a moderate resonance kick would be enough. For this reason, we may consider that the case of Neptune is (at least partially) solved.

Independently of Uranus and Neptune, Saillenfest & Lari (2021) have recently described a mechanism through which a migrating satellite can tilt an initially small-obliquity planet up to over ε = 90° with almost no modification to its spin rate. Their study was prompted by the works of Lainey et al. (2009, 2020), who discovered that the orbital expansion of satellites around Jupiter and Saturn are much faster than previously thought, indicating that efficient mechanisms of tidal dissipation can take place not only in terrestrial planets, as it was known for the Earth-Moon system, but also in giant planets. Such a fast satellite migration has been found to deeply affect the spin-axis dynamics of Jupiter and Saturn (see Lari et al. 2020; Saillenfest et al. 2020, 2021a,b). In fact, while exploring the future evolution of Saturn and Titan, Saillenfest & Lari (2021) showed that if a satellite substantially migrates inwards or outwards, then its influence on the spin-axis precession of its host planet changes over time in a way that is intimately linked to the properties of its equilibrium orbital plane, called the “Laplace plane” (Tremaine et al. 2009). The Laplace plane is well defined everywhere in the parameter space, except at a specific point, called S₁, where it degenerates into a continuum of equilibria. When it comes to the spin-axis dynamics of the planet, S₁ appears as a singularity located at an obliquity ε = 90°, towards which all secular spin-orbit resonances converge if the satellite is massive enough. Over the migration of its satellites, a planet therefore has a chance of being captured in resonance and brought to extreme obliquities. Depending on the migration rate of the satellite, this whole process can take place over the entire lifetime of the planetary system. Importantly, S₁ is surrounded by an unstable region which may lead to the destruction of the satellite once the tilting is achieved (Tremaine et al. 2009; Saillenfest & Lari 2021). New satellites could then be formed from the leftover debris disc (Crida & Charnoz 2012).

In this context, it appears tempting to revisit the case of Uranus’s obliquity and see whether this new satellite migration-destabilisation mechanism could solve at least part of the mystery. Since tidal migration of satellites is a process that never stops, the tilting would not need to be restricted to a short timespan, contrary to previous models, and the resonance involved would not need to be very strong. As a result, this new scenario would neither require a particularly massive satellite, nor a particularly high orbital inclination for the planet. This tilting mechanism was actually proven to work for Saturn on its current orbit, and it is expected to work as well for any giant planet of the Solar System, should it have the adequate long-range migrating satellite (Saillenfest & Lari 2021).

In this article, we aim to assess the feasibility of this tilting mechanism for Uranus and to determine the main parameters that would be required: the mass of the satellite, its migration range and timescale, etc. Due to the multi-scale nature of the mechanism, the current structure of the Solar System actually offers many usable constraints, such as its age, the orbital dynamics of the planets, and the properties of Uranus’s current satellite system.

The basic mechanism considered here is composed of two stages. The resonance capture and adiabatic tilting up to the unstable zone is studied in Sect. 2. We recall the basic mechanism through which a migrating satellite can tilt its host planet (Sect. 2.1), and conduct a preliminary analysis of the parameters needed in the case of Uranus: the mass of the satellite (Sect. 2.2), its migration range (Sect. 2.3), and the resonance capture probabilities (Sect. 2.4). The final destabilisation phase is studied in Sect. 3. We show that the satellite triggers a double synergistic destabilisation involving the planet’s spin axis (Sect. 3.1) and we explore the range of possible outcomes numerically (Sect. 3.2); we then quantify the probability of reproducing Uranus’s current spin state (Sect. 3.3) and determine the conditions under which the satellite collides into the planet (Sect. 3.4). In Sect. 4, we discuss whether the parameters that we obtain appear realistic in the broader context of satellite formation theories and tidal migration mechanisms. We conclude in Sect. 5.

2 Adiabatic tilting up to the unstable zone

2.1 Tilting mechanism

As shown by Saillenfest & Lari (2021), the influence of a given regular satellite on the spin-axis precession of its host planet critically depends on its non-dimensional “mass parameter” η, defined by

$$\eta =\frac{1}{2}\frac{m}{M}\frac{r_{\text{M}}^2}{J_2{R}_{\text{eq}}^2},$$(1)

where m is the mass of the satellite, M is the mass of the planet, J₂ is its second zonal gravity coefficient, and R_eq is a normalising radius (that we take equal to the planet’s equatorial radius). The characteristic length r_M is the distance at which the equilibrium orbital plane of the satellite lies exactly halfway between the equatorial plane and the orbital plane of the planet (the index M stands for “midpoint”). Its definition is closely related to the Laplace radius¹ r_L introduced by Tremaine et al. (2009), as

$$r_{\text{M}}^2=2r_{\text{L}}^5=2\frac{M}{m_{\odot }}{J}_2{R}_{\text{eq}}^2{a}_{\odot }^3{\left(1-e_{\odot }^2\right)}^{3/2}.$$(2)

In this expression, m_⊙ is the mass of the star, and a_⊙ and e_⊙ are the semi-major axis and the eccentricity of the planet on its orbit around the star. The value of r_M for Uranus is about 53 R_eq. We neglect for now the influence of inner satellites in the effective value of J₂ (see e.g. Tremaine et al. 2009). This means either that Uranus possessed a single satellite at that time, or that other satellites were too small and/or too close to Uranus to substantially contribute. We go back to this point in Sect. 4.

For a small mass ratio m/M, the free spin-axis precession frequency of the planet taking into account the influence of its satellite can be approximated by

$${\text{Ω}}_0=p\left(\cos \epsilon +\eta \frac{a^2}{r_{\text{M}}^2}\frac{\sin \left(2\epsilon 2I_{\text{L}}\right)}{2\sin \left(\epsilon \right)}\right),$$(3)

where ε is the planet’s obliquity, a is the semi-major axis of the satellite, and I_L is the inclination of the satellite’s local Laplace plane with respect to the planet’s equator. The leading factor p is the characteristic spin-axis precession rate of the planet (without satellite); it relates to the classic “precession constant” α as $p=\alpha {\left(1-e_{\odot }^2\right)}^{-3/2}$, and it can be written as

$$p=\frac{3}{2}\frac{\mathcal{G}{m}_{\odot }}{\omega a_{\odot }^3{\left(1-e_{\odot }^3\right)}^{3/2}}\frac{J_2}{\lambda },$$(4)

where 𝒢 is the gravitational constant, ω is the spin rate of the planet, and λ is its normalised polar moment of inertia. Equation (3) is valid as long as the satellite oscillates around its local circular Laplace equilibrium on a timescale that is much smaller than the spin-axis precession motion of the planet, which is usually well verified in practice (Saillenfest & Lari 2021). A good approximation for I_L is given by the formula

$$I_{\text{L}}=\frac{\pi }{2}+\frac{1}{2}\text{atan2}\left[-\sin \left(2\epsilon \right),-r_{\text{M}}^5/a^5-\cos \left(2\epsilon \right)\right],$$(5)

which becomes exact in the limit m/M → 0 (Tremaine et al. 2009). It should be noted that even a satellite with a small mass ratio m/M can have a large mass parameter η, and therefore a strong influence in Eq. (3); see Table 1 of Saillenfest & Lari (2021) for examples.

Because of orbital perturbations (e.g. mutual gravitational interactions with other planets), the orbital plane and eccentricity of the planet² are not fixed but vary on a secular timescale. Assuming that the orbit of the planet is long-term stable, its secular motion can be expressed (at least locally) in convergent quasi-periodic series:

$$\begin{array}{lll}z\hfill & =\hfill & e_{\odot }\exp \left(i{\varpi }_{\odot }\right)={\displaystyle \sum _k{E}_k\exp \left(i{\theta }_k\right),}\hfill \\ \zeta \hfill & =\hfill & \sin \frac{I_{\odot }}{2}\exp \left(i{\text{Ω}}_{\odot }\right)={\displaystyle \sum _k{S}_k\exp \left(i{\varphi }_k\right),}\hfill \end{array}$$(6)

where m_⊙, I_⊙, and Ω_⊙ are the planet’s longitude of pericentre, orbital inclination, and longitude of ascending node measured in an inertial frame. The amplitudes E_k and S_k are real constants and the angles θ_k and ϕ_k evolve linearly over time t with frequencies µ_k and v_k:

$$\begin{array}{lll}{\theta }_k\left(t\right)={\mu }_{k}t+{\theta }_k^{\left(0\right)}\hfill & \text{and}\hfill & {\varphi }_k\left(t\right)v_{k}t+{\varphi }_k^{\left(0\right)}\hfill \end{array},$$(7)

where the index k runs over all terms that have non-negligible amplitudes. In Appendix B, we give the orbital solution of Laskar (1990) for Uranus with amplitudes down to 10⁻⁹. This solution can be considered valid (at least in a qualitative point of view) since the end of the planetary migration. In the integrable approximation, the frequency of each term corresponds to a unique combination of the fundamental frequencies of the system, usually noted 𝒢_j and s_j. Table 1 shows the combinations of fundamental frequencies identified for the twenty largest terms of Uranus’s ζ series. In contrast to previous works, in which authors studied the effect of large inclination variations for Uranus or Neptune during early stages of their evolution, we stress that the orbit of Uranus is now very steady. Over billions of years, its inclination barely exceeds 1° with respect to the invariable plane of the Solar System. Yet, as we show below, dramatic obliquity variations can still occur if we leave the system evolve for a sufficient amount of time.

As a result of the motion of the planet’s orbital plane, its long-term spin–axis dynamics is shaped by secular spinorbit resonances, that is, by resonances between the unperturbed spin-axis precession frequency Ω₀ in Eq. (3) and the forcing frequencies µ_k and v_k appearing in Eq. (7). As detailed by Saillenfest et al. (2019b), the largest resonances are those of order 1 in the amplitudes {S_k}, for which the resonance angle is σ = ψ + ϕ_j, where ψ is the precession angle of the planet’s spin axis and j is a given index in the ζ series³. If the planet is trapped in a secular spin-orbit resonance while its satellite migrates, then its mean obliquity ε evolves together with a along a level curve of Ω₀, such that the relation Ω₀ + v_j ≈ 0 is maintained.

Saillenfest & Lari (2021) have shown that for a mass parameter η ⩾ 2, all level curves of Ω₀ with values between p and pη/2 connect ε = 0° to the singular point S₁, which has coordinates (a, ε) = (r_M, 90°); this property is illustrated in Fig. 1. In other words, any secular spin-orbit resonance with a forcing term having a negative frequency v_j such that

$$p⩽\left|v_j\right|⩽p\frac{\eta }{2}$$(8)

would allow for a tilting of the planet between ε = 0° and ε = 90°. If the resonance is large and/or if there are large neighbouring resonances, the planet may even go beyond 90°. For retrograde resonances (i.e. positive frequency v_j), the condition in Eq. (8) allows for a tilting between ε = 180° and ε = 90°. Equation (8) can easily be transposed to higher-order secular spin-orbit resonances: v_j should simply be replaced by a combination of frequencies v_k and µ_k (see Saillenfest et al. 2019b for the variety of possible resonances). However, in the case of Uranus, resonances of order 2 and higher are very small and only a few of them may possibly allow for a capture.

Fig. 1 Level curves of the spin-axis precession rate of the planet as a function of its obliquity and the distance of its satellite. The simplified formula in Eq. (3) is used here with η = 20. Curves in the pink region connect ε = 0° to S1. Curves in the blue region connect ε = 0° to ε = 0° again. Curves in the green region never go to ε = 0° but they connect to S1. The mirror level curves exist for ε > 90° with reversed precession motion. The precession rate along the red level is Ω0 = pη/2. The precession rate along the dark green level is Ω0 = p; it connects (a, ε) = (0, 0°) to S1. The dashed blue curve is the ridge line separating the close and far satellite regimes; it has expression ${\left(a/r_{\text{M}}\right)}^5=\left[\sqrt{{\cos }^2\left(2\epsilon \right)+24}-\cos \left(2\epsilon \right)\right]/6$. See Fig. 17 of Saillenfest & Lari (2021) for examples with other values of η.

2.2 Minimum mass for the satellite

The simplified condition in Eq. (8) offers a quick way to assess whether a given satellite can produce large obliquity variations of its host planet, or which are the optimum parameters of a hypothetical satellite in order to produce a dramatic tilting. In the case of Uranus, p is much smaller than any frequency ∣v_k∣ appearing in its ζ series; therefore, the condition for a satellite to allow for a capture into resonance and tilting up to 90° is simply η ⩾ 2∣ν_★∣/p, where ν_★ is the frequency of the closest resonance (i.e. that with the lowest frequency in absolute value). This condition directly gives the minimum mass required for the satellite that we compute below.

The largest uncertainty in Uranus’s precession constant comes from its polar moment of inertia λ. Helled et al. (2010) also warn readers about the poor accuracy of the spin rate ω measured by Voyager 2, whose value is usually cited (see e.g. Yoder 1995; Archinal et al. 2018); they provide a new estimate for the solid-body rotation period of Uranus that differs from the Voyager value by several percent. As for λ, Hubbard & Marley (1989) give λ = 0.22680 obtained from their interior models using a normalising radius R_eq = 25 559 km. Guided by previous works, we assume λ = 0.23 for simplicity and consider a 10% uncertainty on the product ωλ, used to account for model dependency and the uncertainty on ω (the impact of this choice is commented below). This exploration interval entirely contains the more recent estimates of Nettelmann et al. (2013) and Neuenschwander & Helled (2022), whose values for λ range from 0.21840 to 0.22670, depending on the model used. We obtain values of the frequency parameter p of Uranus ranging from 0.0075 to 0.0092″ yr⁻¹, with minor contribution to the uncertainty range coming from other parameters.

Among all terms in Uranus’s ζ series computed by Laskar (1990), the negative frequency with smallest absolute value is ν₁₅ = 𝒢₇ − 𝒢₈ + s₇ (see Table 1 and Appendix B). Even though this term is very small and unlikely to produce a resonance capture (see below), we can use it to compute the minimum mass ever that a satellite should have in order to tilt Uranus via a secular spin-orbit resonance after the end of the planetary migration. According to Eq. (8) and with our 10% uncertainty range on the value of ωλ, we obtain a minimum mass parameter η_min ranging from about 130 to 150. From Eq. (1), this gives a minimum mass m_min for the satellite lying between about 3.2 × 10⁻⁴ and 3.8 × 10⁻⁴ the mass of Uranus. Two comments can be made about this estimate: Firstly, we note that it does not strongly depend on the poorly known value of ωλ; since we look here for an order-of-magnitude estimate, we adopt its central value for definitiveness and use it in the rest of the article. Secondly, the minimum mass obtained here is not absurdly high as compared to the mass ratios of regular satellites found in the Solar System (see e.g. Murray & Dermott 1999) or obtained in simulations of satellite formation around Uranus and Neptune (Szulágyi et al. 2018). For comparison, a satellite with mass m/M ≈ 3.5 × 10⁻⁴, that is, m ≈ 3 × 10²² kg, would be smaller than Jupiter’s moon Europa. Yet, it is still about ten times the mass of Titania or Oberon, which are today the most massive satellites of Uranus. Hence, if Uranus has been tilted with the mechanism described here, then it should have involved an ancient satellite which has now disappeared.

The need for a hypothetical past satellite does not seem to be a problem because the satellite involved is likely to be destructed anyway during the final stage of the tilting mechanism (see Saillenfest & Lari 2021). We come back to this point below. Now that the need for a hypothetical ancient satellite of Uranus is established, we must determine what could have been its physical and orbital properties in order for the tilting scenario to be plausible. Apart from the very small ν₁₅ term, Uranus’s ζ series features several low-frequency terms that may be good candidates for secular spin-orbit resonances. Table 2 lists the eight closest resonances with the corresponding minimum masses m_min for the hypothetical satellite computed through Eq. (8). For completeness, we include in this table some terms that would require quite a massive satellite to trigger a resonance, even though this possibility appears more unlikely in the context of satellite formation theories (see the discussions in Sect. 4).

2.3 Migration range and velocity

In order to allow for a resonance capture, the migration of the satellite must be slow compared to the oscillations of the resonance angle σ = ψ + ϕ_j, so that the parameter change is close to the adiabatic regime (see e.g. Su & Lai 2020; Saillenfest et al. 2020, 2021b). Among the terms listed in Table 2, this condition may not be verified for all realistic sets of parameters. For each resonance, the oscillation period near the resonance centre (called Cassini state 2) can be computed by modelling the influence of the satellite as an enhanced precession constant α for Uranus (see e.g. French et al. 1993; Ward & Hamilton 2004; Saillenfest & Lari 2021) and by applying the analytical formulas of Saillenfest et al. (2019b) or Su & Lai (2020). As described by Boué & Laskar (2006), the spin-axis precession rate of a planet increases over time if it hosts an outward-migrating close-in satellite, or an inward-migrating far-away satellite. When the satellite reaches some threshold distance, a separatrix appears around Cassini state 2, after which point the planet may be captured in the secular spin-orbit resonance considered. The libration period when the separatrix appears is given in Table 2 for each candidate resonance. An adiabatic resonance encounter could possibly have taken place for Uranus if $T_{\text{lib}}^{\left(\text{(sep)}\right)}$ is much smaller than the age of the Solar System. In this regard, the v₉ and v₁₅ resonances are very unlikely candidates, as adiabaticity would impose the satellite to hardly migrate at all over the Solar System lifetime, in which case the tilting of Uranus is impossible. Indeed, these two terms have the smallest amplitudes in our list, and the resonance width and oscillation frequency scale as the square root of the amplitude of the term (see e.g. Atobe et al. 2004; Li & Batygin 2014). The situation is different for the other candidate resonances in Table 2: we see that Uranus would have enough time to oscillate at least several dozens of time inside the resonance over the Solar System age, potentially allowing for an adiabatic capture and tilting to extreme obliquities. We quantify this point below.

Equation (3) provides a qualitative understanding of how the spin-axis precession rate varies over the space of parameters, and we used it to obtain a quick estimate for the minimum mass of the satellite m_min required for each candidate resonance. However, one may recall that Eq. (3) is approximate and strictly valid only for small satellites. As argued by Saillenfest & Lari (2021), the complete formula of French et al. (1993) is more precise and gives satisfactory results up to larger satellite masses; it can then be used to compute the locations of the resonances and their widths for any distance of the satellite. In this case, there is no closed-form expression for the minimum mass of the satellite but we can compute it in the following way: For a given mass m, the formula of French et al. (1993) is used to obtain an effective precession constant α′. This precession constant is then injected into the equations of Saillenfest et al. (2019b) in order to compute the locations and widths of all resonances in the plane (a, ε). Since a^’ depends on the obliquity, the equations are implicit and must be solved numerically. We repeat this procedure for different masses (e.g. using the bisection method) until we obtain the minimum mass needed for a given resonance to exist somewhere in the plane (a, ε) and go from ε = 0° to 90°. The minimum masses obtained in this way are listed in Table 2; we call them m_k, where k is the index of the resonant term in Uranus’s ζ series. When m_k is large, we see that our simplified analytical expression systematically underestimates its value (i.e. m_min ≲ m_k), even though it gives the correct order of magnitude.

Figure 2 shows the locations and widths of all resonances reachable by the system for different masses of the satellite. It also features the region E₁ where the classic Laplace plane of the satellite is unstable (black striped area; see Tremaine et al. 2009 and Saillenfest & Lari 2021). If the system enters the region E₁, the satellite first transfers to an eccentric equilibrium, and then it suffers from wild orbital changes, that are the most extreme in the vicinity of the singular point S₁. The behaviour of the system in this region is investigated in Sect. 3. As expected, resonances are more numerous for a more massive satellite, because their respective minimum mass criteria are verified simultaneously. When we increase the mass of the satellite, resonances appear near a = r_M (blue zone in Fig. 1), then touch the singular point S₁ (red curve in Fig. 1), at which point we define the mass m_k. Then, if we increase the mass further, resonances expand towards smaller and larger semi-major axes (pink region in Fig. 1). Figure 2 shows that the v₉ and v₁₅ resonances are indeed very thin, and that Uranus also possesses two extremely weak retrograde resonances (v₁₉ and v₃₂). Since m₂ is only slightly smaller than m8 (see Table 2), the v8 resonance appears in panel d but does not connect to S₁ (as in the blue zone in Fig. 1). From these examples, we deduce that in order to tilt Uranus from 0° to 90°, its hypothetical ancient satellite must migrate at least over a distance Δa ≈ 0.2 r_M, either inwards or outwards. This semi-major axis range represents a physical distance of about 10 Uranus’s radii. In order for the satellite to go through this distance in less than the age of the Solar System, its mean migration velocity must be larger than about 6 cm yr⁻¹. Compared to the measured expansion rates of the Moon from the Earth (4 cm yr⁻¹ ; see e.g. Williams & Boggs 2016), Ganymede from Jupiter (11 cm yr⁻¹; see Lainey et al. 2009) and Titan from Saturn (11 cm yr⁻¹; see Lainey et al. 2020), this migration velocity appears quite realistic. Hence, even though the specific dissipation mechanisms acting in the interior of Uranus are essentially unknown yet, at least we are assured that both terrestrial and gaseous planets are able to generate the appropriate range of satellite migration velocities.

It should be noted, however, that the satellite migration needs to be sustained at least up to a distance of a ≳ 0.8 r_M in order for the tilting mechanism to fully operate (see Fig. 2). Since the characteristic radius r_M is quite large for Uranus (r_M ≈ 53 R_eq), this may require the satellite to be fairly big to raise a substantial tidal bulge inside Uranus even at this distance. Assuming that the satellite migrates from a = 0.8 r_M to a = _rM, the equivalent constant parameter k₂/Q of energy dissipation within Uranus can be estimated from classical theories (see e.g. Goldreich & Soter 1966; Efroimsky & Lainey 2007). As shown by Fig. 3, the mean energy dissipation within Uranus is required to be much higher than assumed in historical works, even if the whole process takes as long as the age of the Solar System (for comparison, Goldreich & Soter 1966 gave Q ≳ 10⁵). This new scenario for the tilting of Uranus must therefore be viewed in the context of the fast satellite migration predicted by Fuller et al. (2016) and measured by Lainey et al. (2020) for the satellites of Saturn. Within this change of paradigm, ‘effective’ quality factors Q of order unity are realistic for distant satellites around gaseous planets. For instance, Fuller et al. (2016) predicted Q ≈ 20 for Titan, and Q ≈ 1 for Callisto if a similar resonance-locking mechanism is currently at play in the Jupiter–Callisto system as well. Effective quality factors can also be smaller than unity, because they are not directly linked to a physical tidal lag as they would be in classical theories (see Fuller et al. 2016). The definition of an equivalent constant Q here is therefore only used for comparison purpose.

Alternatively, a chain of mean-motion resonances among several satellites might offer another way to produce a fast satellite migration up to large distances (as it is the case for Ganymede; see e.g. Lari et al. 2020). But again, these aspects intimately depend on the specific dissipation mechanism at play and we can hardly give any more quantitative argument at this stage.

Fig. 2 Locations and widths of first-order secular spin-orbit resonances for different masses of the hypothetical ancient satellite. In each panel, the mass labelled mj is the minimum satellite mass in order to allow for Uranus to reach the vj resonance and be tilted from 0° to 90° (see Table 2). The extent of all resonances is shown in pink and the centre of the vj resonance is highlighted with a red curve. The approximate ridge line separating the close and far satellite regimes is shown by the dashed blue line (same as in Fig. 1). The black striped area is the region E1 where the satellite’s classic Laplace plane is unstable. The red dot is the singular point S1. Resonances are labelled by the frequency vj of the forcing term (arrows).

Fig. 3 Quality factor of tidal energy dissipation within Uranus in order for a satellite to migrate from a = 0.8 rM to a = rM. The ratio k2/Q is computed from classical formulas and plotted here assuming a Love number k2 = 0.1 for Uranus (Gavrilov & Zharkov 1977). The mass of the satellite is labelled above the curves. For comparison, the dotted blue line shows the case of a satellite migration from a = 0.96 rM to a = rM (i.e. a migration range reduced to Δa = 2 Req), as obtained in the numerical experiments of Sect. 3.3.

2.4 Tilting efficiency

The few closest resonances reachable by Uranus and its hypothetical satellite have diverse properties, and apart from general considerations about the level of adiabaticity required, it is not obvious which resonances are more prone to tilting Uranus than others. For this reason, we turn to a numerical exploration of the capture mechanism.

Our setting is similar to the one used by Saillenfest & Lari (2021) for Saturn: we implement the classic equations of motion for the secular spin-axis dynamics (see e.g. Laskar & Robutel 1993; Néron de Surgy & Laskar 1997), in which we modify the precession constant α according to the formula of French et al. (1993). The satellite is assumed to lie on its local Laplace plane at all time and to adiabatically follow its drift as the obliquity is changing. This approximation is valid as long as the orbital timescale of the satellite is much shorter than the motion of the planet’s spin axis, which is well verified in practice. However, it assumes that the satellite’s orbital equilibrium is stable, which is not true in region E₁. Modelling the satellite’s influence as a modified α also requires that the variations of α are slow compared to the spin-axis dynamics, which is not the case in the vicinity of S₁. For these reasons, this approach is very efficient for the stage of capture in secular spin-orbit resonance and adiabatic tilting, but it fails when the system reaches the unstable zone. We use it here to conduct a statistical analysis of the capture process and to perform millions of numerical integrations in a reasonable amount of time. As for the destabilisation stage, it is studied in Sect. 3 below using a self-consistent model.

The tilting mechanism is driven by the migration of the satellite, which runs on a billion-year timescale. Tidal migration of satellites is a complex phenomenon that deeply depends on the internal structure of their host planet, and there is much ongoing research on this topic (see e.g. Auclair-Desrotour et al. 2014; Fuller et al. 2016; Lin & Ogilvie 2021). Satellite migration around terrestrial or gaseous planets are known to involve very different physical mechanisms, and little is known about the diversity of satellite evolutions that can be produced around ice giants like Uranus and Neptune. Instead of relying on one migration scenario, we explore a large range of possibilities for the variations of the satellite’s semi-major axis a over time t. For defmitiveness, we adopt a power law of the type

$$a\left(t\right)=a_0{\left(\frac{t_0+t}{t_0}\right)}^{t_0/\tau },$$(9)

where a₀ is the initial location of the satellite, and t₀ and τ are parameters that have the dimension of time. In the following, we adopt t₀ = 4.5 Gyr and vary the parameter τ, that we call the “migration timescale” (as indeed $\tau ={\left[a/\dot{a}\right]}_{t=0}$. We stress that Eq. (9) is purely ad hoc, and that one could have chosen a linear law instead. The study of more complex migration laws (e.g. involving fast and slow regimes) is beyond the scope of this article; however, if the migration is slow enough to remain in the adiabatic regime, only the total migration range matters, and not its precise rate.

At this point, it should be recalled that, even if we do not model it in detail, the tidal migration of a satellite occurs because of a transfer of angular momentum between the planet’s rotation and the satellite’s orbit. The outward migration of the satellite is therefore accompanied by a decrease in the planet’s rotation rate ω (along with a change in its equilibrium shape). For the largest masses m_k listed in Table 2, and assuming that the increase in the satellite’s angular momentum is fully compensated by a change in ω, the migration of a satellite from a ≈ 0.8 r_M to 1 r_M would produce a change in ω of several percent. Such a variation is of the same order of magnitude as the error on the value of ω measured by Voyager 2 (Helled et al. 2010), and much less than our assumed uncertainty on the value of λ (see Sect. 2.2). Moreover, as the planet slowly cools down on a gigayear timescale after its formation, it is expected to contract slightly (see Bodenheimer & Pollack 1986; D’Angelo et al. 2021), with the result of increasing its spin rate. This opposite effect would further reduce the tidal spin down due to satellites. For these reasons, we neglect the variations of ω and J₂ in our dynamical analysis. We also consider that their changes are small enough so as to present no contradiction with our basic assumption that the spin rate of Uranus is primordial (see Sect. 1).

We are looking for trajectories that go deep inside the unstable region for the satellite. This way, the obliquity not only grows large, but we are also assured that the satellite is destabilised and potentially removed, as required to reproduce the current state of Uranus. Examples of such trajectories are given in Fig. 4. Panel a shows a capture in the v₃ resonance; in this example, the satellite has the minimum mass allowing Uranus to be captured in this resonance and tilted all the way up to 90°. The thin v₁₅ resonance is crossed as well but it hardly affects the trajectory at all. Panel b shows a case where the satellite has mass m₅ but is captured in the v₇ resonance. Indeed, for such a large mass, resonances are numerous and there are many possibilities to tilt the planet. As another example, panel c shows a case where the resonance v₅ is crossed without capture but the planet is then captured and tilted in the v₂ resonance.

However, capture and extreme tilting are not guaranteed even if the satellite has the required mass. There are four cases in which the mechanism can fail to incline the planet enough to reach the unstable zone. The first case is when the satellite does not migrate enough over the Solar System age so as to reach a ≈ r_M. This can happen if it formed too close to the planet (i.e. too far from the Laplace radius) or if it migrates too slowly. The three other cases are illustrated in Fig. 5. In panel a, the planet has a large initial obliquity, which means that it must inevitably cross the resonance separatrix instead of entering smoothly from below (see e.g. Saillenfest et al. 2020). The outcome of separatrix crossings is sensitive to the initial precession phase, and the planet has therefore a given probability of being captured or not (see e.g. Henrard & Murigande 1987; Ward & Hamilton 2004; Hamilton & Ward 2004; Su & Lai 2020). In our case, the capture probability decreases for growing initial obliquity, and panel a is a typical example of failed capture. The numerical experiments of Saillenfest et al. (2021b) show that a fast tidal migration for the satellite tends to smooth the probability profile and to reduce the interval of obliquity leading to a 100%-sure capture. In panel b, the planet is captured but evolves very close to the separatrix; as the resonance width decreases, the trajectory is then pushed out of the resonance before having reached a large obliquity. In panel c, the migration of the satellite is so fast that the resonance is crossed without altering much the obliquity. This is a typical example of non-adiabatic crossing, in which the resonance angle has not enough time to oscillate before the trajectory is already gone.

In order to quantify the ability of each candidate resonance to tilt Uranus and destroy the satellite, we turn to Monte Carlo experiments. In these experiments, the satellite is given the minimum mass m_k required for a given resonance to exist and go to 90° (see Table 2). The satellite is initialised at a₀ = 0.8 r_M, so as to focus on the effects of resonance v_k and minimise the migration range required for the tilting (see Fig. 2). Over a grid of initial obliquity ε₀ and migration timescale τ, we sample random values of initial spin-axis precession angle ψ ∈ [0,2π) and integrate all trajectories numerically until the satellite goes beyond $a^5/r_{\text{M}}^5=\left(10\sqrt{22}-4\right)/39$ (which is the rightmost limit of the unstable zone; see Saillenfest & Lari 2021), and for a maximum timespan of 4 Gyr. This integration timespan is chosen so that the entire tilting mechanism can take place between the formation of Uranus and today, with possibly some time left for the clearing of the debris disc after the satellite’s destruction (see Sects. 3 and 4 below). A statistic is drawn from each sample of initial precession angles. We consider that the system goes deep enough inside the unstable region if the obliquity grows over ε = 80°. This threshold does not mean that the obliquity cannot grow more than 80°, but the model used here fails in this regime, so we prefer to use a conservative limit. The behaviour of the system after destabilisation will be studied in Sect. 3.

Figure 6 shows the result of our experiments for the few closest resonances reachable by Uranus. Even though near-zero initial obliquities are expected from formation models, a variety of different mechanisms may have produced small primordial obliquity excitations (see e.g. Millholland & Batygin 2019; Martin et al. 2020; Rogoszinski & Hamilton 2020, 2021); for this reason, we explore a wide range of values for ε₀. Panel b for m = m₃, is a typical illustration of the different possible behaviours: for τ ≳ 13 Gyr the migration is too slow to tilt the planet in only 4 Gyr, so the success ratio is zero. For τ ≲ 2 Gyr the migration is too fast for an adiabatic capture in the v₃ resonance (as in Fig. 5c). For ε₀ ≳ 20° the trajectories cross the separatrix, leading to an inefficient probabilistic capture (as in Figs. 5a,b). This leaves us with an optimal parameter region with initial obliquity ε₀ ≲ 20° and migration timescale 2 ≲ τ ≲ 13 Gyr, for which the success probability is close to 1. The same structure is observed in panels d and e for resonances v₇ and v₅. In panels a and c, on the contrary, the resonances v₁₅ and v₉ are so thin that the resonance encounter is strongly non-adiabatic in the whole parameter region explored and no trajectory can possibly achieve the tilting of Uranus in the required time interval. In panel f, the resonance v₂ is so large that it produces a 100%-sure capture up to large initial obliquities and down to very small migration timescales; moreover, due to the large amplitude of obliquity oscillations (see Fig. 4c), trajectories can go beyond ε = 80° much before the centre of the resonance does so, which enlarges the yellow region towards the right side of the graph.

From this preliminary analysis, we conclude that a promising candidate resonance is v₃ (i.e. a resonance with the nodal orbital precession mode of Neptune, s₈): it requires a minimum mass for the satellite of only m/M ≈ 4 × 10⁻⁴ and has a high tilting efficiency for a large range of initial obliquities (0^o ⩽ ε₀ ≲ 20°). For larger masses, m/M ≳ 1 to 2 × 10⁻³ the nodal precession spectrum of Uranus allows for several other possibilities: the resonances with v₇= g₅ − 𝒢₇ + s₁ and v₅ = −𝒢₅ + g₆ + s₆ are quite good candidates, but in a more limited range of initial obliquities (ε₀ ≲ 10° and 15°, respectively). As already argued by previous authors (Boué & Laskar 2010; Rogoszinski & Hamilton 2020), the most powerful resonance in terms of tilting efficiency is v₂ = s₇ (i.e. a resonance with the nodal orbital precession mode of Uranus itself). Indeed, it has a 100% capture ratio up to large initial obliquities for Uranus (ε₀ ≲ 45°) and for any conceivable migration rate for the satellite. However, it requires quite massive satellite with m/M ≳ 2.3 × 10⁻³ : even though this value is not absurdly large, it may appear less realistic according to satellite formation theories (see the discussions in Sect. 4).

We are now assured that Uranus can be efficiently tilted to large obliquities through the migration of a hypothetical ancient satellite. The next step is to understand the behaviour of the system once it enters the unstable region. This is the purpose of the next section.

Fig. 4 Examples of capture in secular spin-orbit resonance and tilting to large obliquities. The satellite’s mass and the resonances are labelled as in Fig. 2. Numerical trajectories go from left to right (black dots) over a timespan of billions of years. The migration timescale τ of the satellite is 9 Gyr for panel a and 6 Gyr for panels b and c. Integrations are stopped in the hatched region, where the numerical model used here fails.

Fig. 5 Examples of resonant interactions without tilting all the way to the unstable region. Symbols are the same as in Fig. 4 and the satellite has mass m3 is all panels. The migration timescale τ of the satellite is 9, 2, and 0.2 Gyr for panels a, b, and c, respectively.

Fig. 6 Probability of tilting Uranus as a function of its initial obliquity and the migration velocity of its satellite. On each panel, parameters are sampled on a regular grid, and for each set of parameters, 240 numerical integrations are performed with random initial spin-axis phases ψ ∈ [0, 2π). The colour map shows the fraction of the 240 trajectories that go beyond ε = 80° during the integration. The satellite is initialised at a0 = 0.8 rM and propagated with the migration law in Eq. (9). Numerical integrations are stopped after 4 Gyr or when the satellite goes beyond the limit of the unstable zone. The mass of the satellite is labelled on each panel using the notation of Table 2. The top horizontal axis shows the time required for the satellite to reach a = rM. It can be converted into an effective tidal quality factor Q for Uranus (see Fig. 3).

3 Coupled destabilisation of the planet’s spin and satellite’s orbit

We saw that several secular spin-orbit resonances allow for a 100%-sure capture probability and tilting to ε ≳ 80° in large regions of the parameter space. Such trajectories go deep into the unstable region (hatched zone in Figs. 2, 4, 5); this is a necessary condition, but is not sufficient to reproduce the current state of Uranus. Indeed, assuming that Uranus was indeed tilted by the mechanism presented here, its 98° obliquity implies that the system not only reached the unstable region, but that some additional dynamical mechanism pushed it much further up, beyond the 90° critical point.

The simplified models used in Sect. 2 are not valid when the satellite becomes unstable. In order to investigate what happens to the system at the end of the tilting mechanism, we must use models that include the self-consistent interactions between the planet’s spin and satellite’s orbit. For this purpose, the secular model described by Correia et al. (2011) for the hierarchical planetary case can be transposed to the interactions of a planet and its satellite subject to the attraction of their host star. Due to the highly hierarchical nature of the problem, it is enough to expand all cross interactions to quadrupolar order. The Hamiltonian function is then averaged over the fast orbital and rotational angles. Equations are written in a vectorial form and integrated numerically. In our case, the dynamical variables are the rotational angular momentum of the planet G, the orbital angular momentum of the satellite G₁, and the eccentricity vector of the satellite e₁. The orbital angular momentum of the planet-satellite barycentre around the star and its eccentricity vector (written G₂ and e₂ using the notations of Correia et al. 2011) are taken as quasi-periodic functions of time, varying according to the full orbital solution of Laskar (1990). This way, all planets of the Solar System are included in the orbital motion of the planet-satellite barycentre (but we neglect their direct attraction on the satellite and on the planet’s equatorial bulge). A full description of this model is given in Appendix C. Tidal dissipation in the planet is mimicked by applying a slow variation to the satellite’s semi-major axis as in Eq. (9). In this first exploration of the dynamics, we do not include tidal dissipation within the satellite and do not track its rotational dynamics. This amounts to assuming that the energy dissipated within the satellite always remains negligible as compared to that dissipated in the planet, and that the eccentricity damping of the satellite is much slower than its dynamical evolution. The latter point can be verified a priori: for realistic dissipative parameters (see e.g. Murray & Dermott 1999), the eccentricity damping timescale of the satellite at a distance a ≳ 0.8r_M would count in gigayears. In contrast, Fig. 7 shows that as soon as the system enters (even moderately deep) into the unstable region, the increase in the satellite’s eccentricity unfolds in less than a few hundred thousand years⁴.

Fig. 7 Timescale of the satellite’s eccentricity increase when the system becomes unstable. The colour scale shows the time needed for the eccentricity to be multiplied by a factor 100, as computed from the equations linearised at the unstable equilibrium point (see Tremaine et al. 2009; Saillenfest & Lari 2021).

3.1 Destabilisation mechanism

Figure 8 shows a typical example of simulation obtained with the coupled model. The migration timescale τ is chosen to be in the middle of the parameter region producing a 100%-sure resonance capture (see Fig. 6). The first portion of the trajectory is identical to what was predicted in Sect. 2 using the simplified dynamical model (see e.g. Fig. 4c for a similar satellite mass). Indeed, the oscillations of the satellite around the Laplace equilibrium are very fast compared to the motion of the planet’s spin axis; the previous assumption that the satellite instantly follows the Laplace plane is therefore a very good approximation. The smooth increase in the satellite’s orbital inclination while its longitude of node closely oscillates around zero is a direct manifestation of this adiabatic dynamics. Once the system reaches the border of the unstable zone, however, the dynamics brutally change: Fig. 8 shows that the satellite’s eccentricity increases very quickly, which releases the planet from the secular spin-orbit resonance (σ₂ circulates again). From this point on, the satellite’s orbit starts to evolve chaotically and it can reach almost any eccentricity and inclination value, while switching sporadically between being prograde and retrograde with respect to the planet’s spin. Interestingly, these strong chaotic orbital variations affect the planet’s obliquity, which evolves pretty much as a random walk. At this stage of the evolution, there is no particular dynamical barrier anymore at ε = 90°; in the example of Fig. 8, we see that the planet’s spin can very well become retrograde and reach the current obliquity of Uranus (dotted line).

The cause of the satellite’s destabilisation is well understood (see Tremaine et al. 2009; Saillenfest & Lari 2021); the reason for the random walk in obliquity is less obvious. In order to investigate this mechanism further, we continued the numerical integration of Fig. 8 for one more gigayear with different physical parameters. The results are shown in Fig. 9. Because of the strongly chaotic nature of the dynamics, the slightest change in initial condition soon leads to a strong divergence of trajectories. In this regard, the evolutions shown in panels a, b, and c are indistinguishable: from a qualitative point of view, they may just as well be three different realisations of the same chaotic process. Hence, from panels a, b, and c, it seems clear that the chaotic obliquity variations do not depend on the residual migration of the satellite. Indeed, the destabilisation of the system is an irreversible process, and the wild orbital variations of the satellite persist even if its tidal migration stops or is reversed. Changes in the satellite’s migration are expected to take place in reality when the satellite becomes unstable; these changes can be due to tidal dissipation within the satellite (not modelled here), to the break of the tidal resonant link (in case of a migration driven by a mechanism similar to that described by Fuller et al. 2016), or to the reversed effect of energy dissipation when the satellite becomes retrograde. The absence of a visible influence of the satellite’s migration properties on the chaotic dynamics guarantees that our idealised model still yields a qualitatively relevant description of the dynamics when the system goes deep into the unstable region.

Conversely, panels d and g of Fig. 9 show that the orbital nodal precession of the planet plays an important role in the chaotic behaviour of the obliquity. If the planet’s orbit is held fixed (panel g), its obliquity still oscillates as a result of the orbital variations of the satellite, but the large-scale obliquity kicks are suppressed. Panel d shows that one single harmonic of nodal precession is enough to restore the sporadic jumps in obliquity (and term k = 2 is the largest). Other panels in Fig. 9 show that the mass of the satellite is directly related to the amplitude of the erratic variations in obliquity: the diffusion is slower for a less massive satellite, and totally suppressed for m = 0. Panel i confirms that if the satellite is removed (e.g. because of a collision; see below), the planet’s obliquity instantly freezes.

From these comparisons, we deduce that the large-scale chaotic variations in the planet’s obliquity are due to a double synergistic destabilisation: (i) the satellite’s orbit is strongly unstable because of the combined attractions of the planet’s equatorial bulge and of the star, as previously described by Tremaine et al. (2009) and Saillenfest & Lari (2021); (ii) by a partial conservation of angular momentum, the wild orbital variations of the satellite produce fast oscillations and some diffusion in the planet’s obliquity (direct effect), but also large-amplitude variations in the planet’s spin-axis precession frequency. Because of these frequency variations, the planet goes in and out of the nearby secular spin-orbit resonances, resulting in strong additional obliquity kicks. These obliquity kicks reinforce the satellite’s instability and make it explore an even wider range of orbital elements. Since all secular spin-orbit resonances converge to the singular point S₁ (see e.g. Fig. 2), all available resonances are actually “nearby” and contribute to the obliquity kicks. The largest of them is the v₂ resonance, that is, a resonance with the nodal orbital precession mode of Uranus itself. Fig. 9d shows that it contributes most to the obliquity kicks.

Fig. 8 Example of simulation using the fully coupled secular model. Black and blue are used for quantities related to the planet’s spin axis and satellite’s orbit, respectively. The orbital elements of the satellite are measured with respect to the equator and equinox of date (we indicate this reference frame by the index Q). In this example, the secular spin-orbit resonance angle is σ2 = ψ + ϕ2. The mass of the satellite is m/M = 2.2 × 10−3; it is initialised close to the local Laplace equilibrium at a = 37 Req and made migrating outwards with timescale τ = 6.5 Gyr.

Fig. 9 Chaotic diffusion of obliquity for different physical parameters. The black trajectory is the same as in Fig. 8. At t = 4 Gyr (vertical bar), parameters are instantly changed and the system is integrated for 1 Gyr more (red part). Panel a: no parameter change. Panel b:the satellite no longer migrates. Panel c: the satellite migrates inwards. Panel d: the planet’s orbital variations are restricted to the 2nd harmonic of ζ. Panel e: the mass of the satellite is halved. Panel f : the mass of the satellite is zero. Panel g: the planet’s orbit is held fixed. Panel h: the planet’s orbit is held fixed and the satellite’s mass is halved. Panel i: the planet’s orbit is held fixed and the satellite’s mass is zero.

3.2 Exploration of the parameter space

The previous example shows that the double destabilisation mechanism, while putting an end to the stable adiabatic drift inside the resonance, is able to pump the obliquity high enough to reproduce the current state of Uranus. We must now investigate the conditions under which the double destabilisation is triggered and estimate its efficiency. To this end, we turn to Monte Carlo experiments.

In our first experiment, the migration timescale of the satellite is set to τ = 7.5 Gyr and the initial obliquity of the planet is ε₀ = 3° for all numerical integrations. This puts the system in the parameter region producing a 100%-sure capture and tilting in most resonances. As shown in Fig. 6, the choice of ε₀ is not critical, as the same efficiency can be obtained in a large range of initial obliquities. The initial semi-major axis of the satellite is a₀ = 42 R_eq (that is, about 0.8 r_M). The satellite is started close to the circular Laplace equilibrium, with an eccentricity e = 10⁻⁴ and an inclination I_LP = 10⁻⁴ rad with respect to the local Laplace plane. For a given value of the satellite’s mass m, we draw a random sample of initial conditions for the planet’s spin-axis precession angle ψ ∈ [0, 2π) and for the satellite’s argument of pericentre e [0, 2π) and longitude of ascending node e [0, 2π) measured with respect to the local Laplace plane. All initial conditions are propagated numerically for 4 Gyr using the fully coupled secular model, and the maximum obliquity value reached over the full numerical integration is recorded.

The result is shown in the panel a of Fig. 10. We can identify the effects of each individual resonance which progressively appears as we take a larger mass for the satellite, and whose bottom extremity sweeps over the satellite’s initial location at a₀ ≈ 0.8 r_M (see Fig. 2 for the geometry of the resonances). Beyond a given mass value, the resonance still exists, but its bottom extremity is closer to the planet than 0.8 r_M, so the system does not encounter it during the outward migration of the satellite. The minimum mass values estimated previously (see Table 2) correspond quite accurately to the top of the obliquity peaks; their locations are marked by a small line in panel c. As expected, resonances v₁₅ and v₉ only produce tiny obliquity bumps because they are too thin to allow for an adiabatic capture in less than the age of the Solar System (see Sect. 2.3). We also checked that for a very fast satellite migration (τ ≲ 3 Gyr), the peak corresponding to resonance v₇ is greatly reduced and barely reaches ε = 30° because the probability of adiabatic capture falls to zero (see Fig. 6). The minimum mass estimate for resonance v₈ appears quite far from the corresponding obliquity peak in Fig. 10, because resonance v₈ is very distorted by the large neighbour resonance v₂. In fact, resonance v₂ is so large that it almost completely overlaps resonance v₈; it is also surrounded by several higher-order resonances that have a non-negligible width and contribute to the production of chaos (see Saillenfest et al. 2019b). By carefully analysing the trajectories in Fig. 10, we found that some of them feature captures in second-order resonances. An example of capture in a resonance with critical angle σ = 2ψ + ϕ₂ + ϕ₅ is given in Appendix D; this example confirms that the dynamics in the neighbourhood of resonance v₂ are very rich and not limited to interactions with harmonic ϕ₂ taken in isolation.

Panel a of Fig. 10 reveals that the double destabilisation mechanism described above can routinely reproduce Uranus’s obliquity through a capture into resonances v₇, v₅, v₂, v₈, or v₆. The minimum mass of the satellite for which successful trajectories are obtained is m/M ≈ 1.2 × 10⁻³ (i.e. m ≈ m₇). For smaller masses, no trajectory in our sample goes beyond the current obliquity of Uranus. As expected from the previous sections, captures into resonance v₃ for m/M ≈ 4.4 × 10⁻⁴ are perfectly stable and adiabatic, and trajectories do reach the unstable region, but then they barely manage to go over ε = 90°. This is because the satellite is too small to produce a strong back reaction on the spin axis when it becomes unstable, and the v3 resonance is quite isolated from other strong resonances. Both these aspects contribute to slow down the chaotic diffusion in obliquity, and one would need to greatly expand the integration timespan in order to give enough time for trajectories to diffuse high enough (see Sect. 3.1).

Moreover, the satellite can actually not go on kicking and jumping all over the place for a long time, because sooner or later it reaches an orbit that makes it collide into the planet. Indeed, our simulations show that the secular destabilisation of the satellite is strong enough to destroy it, similarly to asteroid disruption around white dwarfs (O’Connor et al. 2022). In panel a of Fig. 10, all integrations are continued for 4 Gyr even if the pericentre distance q of the satellite becomes smaller than 1 R_eq (this situation is unphysical). For comparison, panel b shows the same Monte Carlo experiment, but for which integrations are stopped in case of satellite-planet collision (i.e. if ever q ⩽ 1 R_eq). The fraction of collisional trajectories is shown in panel d (black histogram): we see that the transition from 0% to nearly 100% collisions is very steep for all resonances. This means that if the system goes deep in the unstable zone, collision becomes inescapable. Because of the existence of collisional trajectories, the destabilisation of the satellite can be thought of as the start of a countdown that will eventually put an end to the evolution of the system. In order to reproduce the current state of Uranus, we need that the erratic motion of the satellite brings the planet’s obliquity high enough before the collision. As shown in panel b, the maximum obliquity values reached before collision are substantially smaller than what we observed in panel a. A satellite with mass m/M ≈ 1.2 × 10⁻³ now only marginally allows for Uranus’s obliquity to be reached.

The effect of the satellite’s impact on the spin rate and obliquity of the planet can be estimated from the conservation of angular momentum (see Appendix F and the related discussion in Sect. 4.2 below). For the range of masses considered in Fig. 10, the relative change in spin rate due to the impact is a few percent at most, and the change in spin direction is no larger than a few degrees.

Depending on the cohesive strength of the satellite and the rate of its eccentricity increase, the satellite may not even reach collision but be torn apart by tidal forces farther away from Uranus. Similarly to the amount of tidal dissipation within the satellite (that we chose to ignore for now), the exact Roche limit for the satellite depends on its internal structure and composition and can only be speculative at this stage. For definitiveness, we adopt a limit of 2 R_eq, which is about halfway between the fluid and rigid Roche limits (assuming that the satellite is homogeneous and has a density similar to Uranus’s current regular satellites; see e.g. Canup & Esposito 1995; Crida & Charnoz 2012; Hesselbrock & Minton 2019). In fact, the satellite is not expected to be completely torn apart in only one passage, but rather be progressively eroded and form a torus of material. If the satellite is differentiated, its denser core may survive much longer than its outer layers, possibly down to collision (Canup 2010). Some fraction of material may also be lost in space (see e.g. Malamud & Perets 2020; Brouwers et al. 2022). In any case, the remaining debris would reorganise into a thin equatorial disc confined below the fluid Roche limit (see e.g. Hyodo et al. 2017 and discussions in Sect. 4.3). Hence, if the satellite is tidally disrupted, there should actually be a continuous (yet fast) transition regime from the solid satellite, whose influence on the planet’s spin axis is strong, to the confined debris disc, whose effect on the planet’s spin axis is negligible. Here, we do not model these phenomena, and consider that the satellite instantly disappears when crossing 2 R_eq. This approximation may seem quite crude but the results obtained are still informative, as they can be compared to the pure collisional case described above: a more realistic behaviour of the system probably lies somewhere between these two extremes. The maximum obliquity reached by the planet before the satellite crosses the Roche limit is shown in panel c of Fig. 10. Obliquity values are globally lower than in panel b, but the picture is not very different in a qualitative point of view. Indeed, almost all trajectories that cross the Roche limit reach collision shortly afterwards (see panel d).

Fig. 10 Maximum obliquity reached in the course of numerical integrations using a fully coupled secular model. For each value of the satellite’s mass (bottom axis), 96 numerical integrations are performed with random initial conditions as described in the text. The initial semi-major axis of the satellite is a0 = 42 Req and it is made migrating according to Eq. (9) with timescale τ = 7.5 Gyr. Panel a: the maximum obliquity reached over a full duration of 4 Gyr is shown by a small dot. The horizontal dotted line shows Uranus’s current obliquity. Panel b: same as panel a, except that integrations are stopped if the satellite collides into the planet. If a collision occurs, the dot on the graph shows the maximum obliquity reached before the collision; the time of reaching this maximum may be different from the collision time. Panel c: same as panel b, except that integrations are stopped if the satellite goes below the Roche limit (taken equal to 2 Req). The mass estimates from Table 2 are labelled (small vertical lines). Panel d: histogram showing the fraction of trajectories featuring a collision (black). The magenta bars on top show the small fraction of trajectories for which the satellite goes below the Roche limit but does not collide into the planet within the integration timespan.

3.3 Probability of reproducing the current state of Uranus

It is now clear that starting from a small axis tilt, the current obliquity of Uranus can be reached through the tidal migration of an ancient satellite, by a capture in secular spin-orbit resonance followed by a global instability. To reproduce the current state of Uranus, however, the planet must not only reach Uranus’s obliquity, but stop at its value. In order to estimate the probability of successful trajectories, we now examine the configuration of the system at collision, that is, at the time the satellite is expected to be destructed and leave the planet’s obliquity in a fossilised state. Panel a of Fig. 11 presents the same Monte Carlo simulations as in Fig. 10, but showing the obliquity of the planet at the time of collision (instead of the maximum obliquity reached). The intervals of mass devoid of data points are regions in which no trajectory in our sample features a collision. Obliquity values are sensibly lower than in Fig. 10b, which means that the satellite generally does not collide when the planet’s obliquity is at its maximum, but somewhat later on. We see that collisions occur when the spin axis is either plainly prograde or retrograde, but almost no collision occurs at all while the obliquity is ε ≈ 90°. The explanation for this bizarre phenomenon is given in Sect. 3.4. The colour code in Fig. 11 distinguishes trajectories ending with ε < 90° (prograde family) from those ending with ε > 90° (retrograde family). Interestingly, the bulk of the retrograde family is located right at the obliquity value of Uranus (see histogram on the right). For this reason, we define a run as “successful” to reproduce the current state of Uranus when it ends up in the retrograde family.

Panel b of Fig. 11 shows the probability of success as a function of mass (i.e. it shows the fraction of trajectories in our sample that feature a collision and for which the planet’s obliquity at collision is larger than 90°). In the most favourable region, probabilities reach 60 to 70% (the highest peak is 73%). In other words, reproducing the current state of Uranus is the most likely outcome in this range of parameters. The resonance responsible for the obliquity increase in this region is ν₂ (partially overlapping with resonance ν₈). This is not a coincidence: resonance ν₂ is the largest secular spin-orbit resonance that appears in the orbital series of Uranus.

Panels c and d of Fig. 11 present the same Monte Carlo experiment, but where the satellite is supposed to be instantly destructed when it goes below the Roche limit. The same kind of structure as in panel a is observed, but the division between the prograde and retrograde families is less marked. The probability to end up retrograde is quite smaller: it reaches 20% to 30% in the most favourable region (the highest peak is 29%). Again, this can be explained by the coupled destabilisation mechanism: as the satellite is removed earlier than in panel a, the system is offered less time to substantially excite the obliquity.

Because of the geometry of secular spin-orbit resonances, the intervals of mass showing the maximum success probability in Fig. 11 are specific to the choice of initial semi-major axis a₀ of the satellite. The choice a₀ ≈ 0.8 r_M has been made to roughly minimise the distance over which the satellite must migrate in order to reach the unstable region (see Sect. 2.3). However, similar success probabilities can be obtained if the satellite starts closer to the planet and migrates further out. In this case, the same resonances are reachable up to larger satellite masses. In order to clarify this point, we repeated our Monte Carlo experiment for different initial locations of the satellite. Figure 12 shows the probability of ending in the retrograde family as a function of both the satellite’s mass m and its initial distance a0. The effect of each individual resonance is clearly visible. As expected, resonance ν₂ for m/M > 2 × 10⁻³ produces the largest fraction of successful trajectories, but resonances ν₅ and ν₇ also generate substantial probabilities of success in particular regions of the parameter space. The range of mass having the highest success ratio increases when we decrease the initial distance of the satellite; this is a direct consequence of the geometry of secular spin-orbit resonances. Below a given initial distance, however, the satellite has not enough time in 4 Gyr to reach the unstable region, so the success probability is zero. Larger migration velocities are needed to counteract this problem; they allow one to initialise the satellite closer to the planet (compare panels a, b, and c).

If the satellite is initialised closer to the planet, the system can catch a resonance that brings it to higher obliquities before triggering the instability (see e.g. resonance v₃ in Fig. 2b, to be compared with the same resonance in Fig. 2a). This can be a way to increase the success ratio for small satellites. And indeed, because of this property, panel c in Fig. 12 shows that resonance v3 is able to generate a small fraction of successful trajectories if the satellite is started at a₀ ≲ 35 R_eq and migrates fast (see the violet shade for m/M ≈ 7 × 10⁻⁴). The success ratios, however, remain very small in resonance v₃.

In summary, Fig. 12 reveals that the regions of the parameter space allowing the current state of Uranus to be reproduced are quite vast. Through a capture in resonance v₂, the probability of obtaining a final retrograde spin axis (with an obliquity similar to Uranus’s) can be larger than 80%. Other promising resonances are v₅ and v₇, with success ratios up to 50%. Thanks to the large width of resonance v₂, the tilting mechanism can also operate down to very small migration ranges: this is the case in the uppermost portion of the graphs, for initial semi-major axes a₀ ≳ 50 R_eq. To further illustrate this point, Fig. 13 shows the time elapsed before collision and the migration range covered by the satellite in these Monte Carlo experiments. In some narrow range of parameters, we see that the current state of Uranus can be reproduced even if the satellite migrates over distances as small as 2 R_eq. This reduced migration range increases the effective quality factor Q that would be required for Uranus (see Fig. 3).

Fig. 11 Final value of the planet’s obliquity in the simulations of Fig. 10. Panels a and b: the satellite is removed when it collides into the planet. Trajectories in which the satellite does not collide within the integration timespan are not shown. The main resonances are labelled. The obliquity histogram is given on the right. Obliquity values larger than 90° are highlighted in red; their ratio over all trajectories (including those without collision) is shown in panel b as a function of mass. Panels c and d: the satellite is removed when it goes below the Roche limit. Trajectories in which it never does are not shown. Obliquity values larger than 90° are highlighted in blue.

Fig. 12 Probability of destructing the satellite and producing a retrograde obliquity as a function of the mass of the satellite (bottom axis), its initial semi-major axis (vertical axis), and its migration timescale (labels). Each pixel features 96 numerical integrations over 4 Gyr with random initial angles (same setting as in Fig. 10). In the upper panels, the satellite is removed when it collides into the planet (q < Req). In the lower panels, the satellite is removed when it goes below the Roche limit (q < 2 Req). White pixels mean that the satellite survived in all 96 simulations. Coloured pixels show the fraction of simulations in which the satellite is destructed and the planet’s spin axis ends in the retrograde family. Resonances are labelled in panel a.

Fig. 13 Collision time and migration distance covered by the satellite in the simulations of Fig. 12. Top panels: integration time at collision. Bottom panels: distance covered by the satellite before collision. Similar graphs are obtained when considering the time of crossing the Roche limit (with only slightly smaller integration times and migration distances).

3.4 Collision keyholes

The peculiar distribution of obliquity in Fig. 11, showing accumulations in a prograde and a retrograde spin “families” and much depletion around ε ≈ 90°, deserves a more in-depth analysis. We write I_Q the orbital inclination of the satellite measured with respect to the planet’s equator and I_C its inclination measured with respect to the planet’s orbital plane (Q and C stand for “equator” and “ecliptic”). When examining the orbital elements of the satellite at the time of its collision, we find that the values of I_Q are localised in two very narrow bands that can be visualised in Fig. 14. These two bands are localised at I_Q ≈ 55° and 125°; we see that they are less narrow when the satellite crosses the Roche limit than when the satellite impacts the planet. In contrast, we find no preferential orientation of the satellite with respect to the planet’s orbital plane, with inclinations I_C covering almost the entire range from 0° to 180°. This seems to indicate that the pathways to collision (‘collision keyholes’) involve predominantly planet-satellite interactions.

Figure 15 shows an example of the track left by the satellite in the plane (e, I_Q) as the system wanders in the chaotic zone. Even though the obliquity of the planet covers a large range of values during the evolution (see Fig. 8), the two collision keyholes are extremely localised in inclination I_Q: we clearly see that no collision can possibly occur outside of I_Q ≈ 55° and 125°, which explains the existence of the narrow bands in Fig. 14. We performed various numerical tests, and found that the same chaotic structure as in Fig. 15 is produced even in the most simple case of a massless satellite around a planet with fixed orbit and fixed spin axis. This may not be surprising, because collisions occur on a very short timescale compared to the orbital and spin-axis motions of the planet, and we only consider relatively small satellites here. The nature of the collision keyholes should therefore have a simple dynamical origin that we investigate here.

The Hamiltonian function describing the secular motion of a massless satellite around a planet with fixed orbit and spin axis can be written

$$\mathscr{H}=k_{\text{p}}{\mathscr{H}}_{\text{p}}+k_{\odot }{\mathscr{H}}_{\odot },$$(10)

where the first term comes from the interaction of the satellite with the equatorial bulge of the planet, and the second term contains the orbital perturbation produced by the star. The explicit expression of the constant factors {k_p, k_⊙} and of the Hamiltonian functions {ℋ_P,ℋ_⊙}expanded at quadrupolar order are given in Appendix E. Because ℋ_P contains a factor (1 − e²)^−3/2, the first term in Eq. (10) becomes strongly dominant whenever the eccentricity grows large. Indeed, the contribution of k_P ℋ_P produces a fast precession of the satellite’s argument of pericentre ω_Q and longitude of ascending node Ω_Q, at rates

$$\begin{array}{lll}{\dot{\omega }}_{\text{Q}}=k_{\text{P}}=\frac{5{\cos }^2{I}_Q-1}{\sqrt{\mathcal{G}Ma}{\left(1-e^2\right)}^2}\hfill & ;\hfill & {\dot{\text{Ω}}}_{\text{Q}}=k_{\text{P}}=\frac{-2\cos I_Q}{\sqrt{\mathcal{G}Ma}{\left(1-e^2\right)}^2},\hfill \end{array}$$(11)

with the result of averaging out the angular dependency in k_⊙ℋ_⊙. As a consequence, the eccentricity and equatorial inclination of the satellite become constants of motion when the eccentricity grows large. Because of this property, the eccentricity of the satellite can never grow arbitrarily close to 1: the contribution from k_Pℋ_P makes it freeze at some point and leaves it with no alternative but decreasing again.

In order to determine the orbital configurations in which the eccentricity can grow the largest, we must focus on the intermediate regime in which k_Pℋ_P is dominant but the contribution from k_⊙ℋ_⊙ is not yet negligible. If we adopt a perturbative approach and consider that the two degrees of freedom are weakly coupled, any combination ξ = iω_Q + jΩ_Q can become a resonant angle, where (i, j) ∈ ℤ² \ {(0, 0)}. The possible resonances at first order in the perturbation are given in the left column of Table 3. Their nominal location is obtained by imposing ξ = 0 from Eq. (11), which results in one specific value for the inclination I_Q, that we note I₀.

The properties of these resonances can be studied with the method described in Appendix E. In the vicinity of each resonance, the averaged system admits a constant quantity K that links the eccentricity and the inclination of the satellite, in the same way as the Kozai constant (see e.g. Lidov 1962; Kozai 1962; Saillenfest et al. 2016). The expression of K associated with each resonance is given in the right column of Table 3. Because of this constant quantity, the resonance width ΔI_Q in inclination has an equivalent width Δe in eccentricity (except for the “resonance” ξ = Ω_Q for which the eccentricity itself is constant; see Table 3). The widths of each individual resonance is shown in Fig. 16 as a function of the distance of the satellite a and obliquity of the planet ε.

For a fixed value of obliquity ε, the extension of the resonances can be visualised in the space of the satellite’s orbital elements, as illustrated in Fig. 17. As the resonances all feature a distinct constant quantity K, they lie in a different space and we cannot directly judge whether they are overlapping or not in the sense of Chirikov (1979). However, Fig. 17 gives a precise idea of the variations of orbital elements expected for the satellite under the influence of each resonance. For a ≲ 0.2 r_M, resonances are small and well separated from each other for any obliquity value. When the satellite is close to a = r_M, on the contrary, and especially if the obliquity is large, resonances allow for extreme variations in eccentricity and inclination, in accordance with what we observed in previous simulations⁵. Moreover, it is clear from Fig. 17 that no resonance at all is able to substantially increase the eccentricity over large intervals of equatorial inclination located around I_Q = 0°, 90°, and 180°. This property explains the peculiar shape of the black region in Fig. 15; the two spikes are produced by the two triplets of resonances that have the largest Δe (see the horizontal dashed lines).

In Fig. 16, we see that the resonance that produces by far the largest variations in eccentricity is 2ω_Q + Ω_Q (which becomes 2ω_Q − Ω_Q if the orbit is retrograde). This resonance, however, is very narrow in inclination: a few degrees at most. Its nominal location is I_Q ≈ 56° (or I_Q ≈ 124° for a retrograde orbit), as given in Table 3; this location closely matches the collision keyholes in Figs. 14 and 15. We deduce that the secular resonance 2ω_Q + Ω_Q (when the orbit is prograde) or 2ω_Q − Ω_Q (when the orbit is retrograde) is the main responsible for the collision of the satellite into the planet. This resonance is well known for artificial satellites; it has been studied in depth by Daquin et al. (2022) in the more complicated case of lunisolar perturbations. In our case, Figs. 16 and 17 reveal that the widths of this resonance abruptly drop to zero for obliquity values ε = 0°, 90°, and 180°; this explains why no collision at all occurs around these values of obliquity, as noticed previously in Fig. 11. One can verify this property by running a simulation with the planet’s obliquity fixed to ε = 90°: we obtain a graph similar to Fig. 15, but in which the two spikes are smoother, like eroded, and the collision limit is out of reach for the satellite. The strongest resonance that persists for ε = 90° is ξ = ω_Q − Ω_Q located at I_Q ≈ 46° (or equivalently ω_Q − Ω_Q at I_Q ≈ 134°), but it is not quite large enough for the satellite to reach collision (see Fig. 16).

It should be noted, however, that the simplified model used here to describe the secular resonances predicts an equal probability for collisions at I_Q ≈ 56° and 124°. Equal probabilities are indeed obtained for small satellites, but more massive satellites happen to collide much more frequently in the retrograde keyhole (see Fig. 14). A careful examination of Fig. 15 reveals that the retrograde peak is indeed slightly bigger. This subtlety cannot be explained using a simplified restricted model.

Another property of Fig. 14 needs to be clarified: the extreme majority of satellite destructions (collisions of Roche crossings) that occur when ε > 90° (coloured dots) go through the retrograde keyhole at I_Q ≈ 124° and almost never in the prograde keyhole, whatever the mass of the satellite. In other words, in almost every trajectory that successfully reproduces the current state of Uranus, the ancient satellite collides into Uranus while being retrograde with respect to its equator. This property can be explained by the “countdown to collision” mentioned earlier: As the satellite wanders about in the chaotic zone, it will reach sooner or later a collision keyhole. Because of the partial conservation of the planet-satellite angular momentum, the planet’s spin axis diffuses much more efficiently from ε < 90° to ε > 90° (or back again) when the satellite is highly retrograde with respect to its equator (this can be verified in Fig. 8). During the short interval while ε ≈ 90°, collision keyholes are closed and the satellite is temporarily safe from destruction. However, when the obliquity becomes substantially larger than 90° and the collision keyholes open again, the satellite has a much larger probability of colliding right away, while still being retrograde with respect to the planet’s equator, than going back prograde and then collide. Figures 18 and 19 present an example of trajectory for both cases. Figure 18 shows the most common situation, in which the satellite is destructed right away after having become retrograde and driven the obliquity to values larger than 90°. Figure 19 shows the very improbable situation (much less than 1% of cases) in which the satellite goes back prograde before being destructed. One can note that these two trajectories have been obtained using two very different migration timescales; this highlights the large variety of parameters from which the current state of Uranus can be reproduced via the mechanism described here.

Fig. 14 Same as Fig. 11, but showing the orbital inclination of the satellite at the time of its destruction. Coloured dots show obliquity values larger than 90°; grey dots show obliquity values smaller than 90° (same as in Fig. 11).

Fig. 15 Chaotic region explored by the satellite. The black dots show the track left during 4 Gyr by a single numerical integration. The physical parameters are the same as in Fig. 8. The blue and red vertical lines show the Roche and collision limits, respectively. The white horizontal dashed lines show the locations of the secular resonances that have the largest width in eccentricity (see text).

Fig. 16 Widths of first-order secular resonances for the satellite as a function of the parameters. Only prograde resonances are shown. The widths of retrograde resonances are obtained by inverting the sign of ΩQ in the titles. For each resonance, the inclination and eccentricity of the satellite are linked through a constant quantity K given in Table 3. The fixed value of K used to plot this figure is given in Appendix E. The colour gradient shows the width of the resonance as measured in inclination (top row) and in eccentricity (bottom row). Because of its overly large ΔIQ, the width displayed in the top right panel is shown divided by 10 (as illustrated in Fig. 17, it reaches 180° at a = rM and ε = 90°).

Fig. 17 Location and width of first-order secular resonances for the satellite. Each panel is plotted for one fixed value of the planet’s obliquity ε from 0° (left) to 90° (right). Resonances are labelled in the leftmost panel by their respective critical angle ξ (omitting the Q indexes). The coloured intervals show the extremal values of equatorial inclination IQ spanned by the resonance separatrix. Inside each resonance, a conserved quantity links IQ to the eccentricity e, as shown by the colour gradient.

4 Discussions

4.1 Forming Uranus’s ancient satellite

In Sect. 3, we have seen that when trying to reproduce the current spin state of Uranus via the migration of a single ancient satellite, large probabilities of success can be obtained for a satellite mass m/M ≈ 1.7 × 10⁻³ (i.e. about the mass of Jupiter’s moon Ganymede, m ≈ 15 × 10²² kg) or larger. In order to achieve the tilting of Uranus in a billion-year timescale without invoking extravagant migration rates, this satellite should be formed during the early stages of the Solar System evolution and start its outward migration at an initial distance a₀ ranging from about 40 to 50 R_eq. We must now discuss whether the existence a such a satellite around Uranus could appear realistic in the context of satellite formation theories.

For postulating that Uranus and Neptune had small initial obliquities, our main hypothesis is that they formed in a similar manner as Jupiter and Saturn, with the infall of gas dominating their final spin state (see Sect. 1). With this hypothesis in mind, it seems natural to consider that a primordial generation of satellites formed around Uranus and Neptune through a mechanism similar to the Galilean satellites around Jupiter or Titan around Saturn. Mosqueira & Estrada (2003a,b) modelled the formation of the main satellites of Jupiter, Saturn, and Uranus in a gaseous circumplanetary disc. According to their work, Uranus today “lacks” a big satellite, that should have been formed at a distance of 57 Req and would have been Uranus’s analogue of Ganymede or Titan. The authors interpreted this discrepancy as the signature of a substantial inward migration of Uranus’s protosatellites in the gaseous disc; they did not notice that a satellite at this distance around Uranus is unlikely to survive because it is close to the unstable zone (see Tremaine et al. 2009; Saillenfest & Lari 2021). Even though major improvements have been made since then in satellite formation theories (see e.g. Batygin & Morbidelli 2020), a large distant satellite does seem to lack in Uranus’s system when compared to the systems of Jupiter and Saturn. This missing satellite may resemble what we are looking for, but some caution is needed. On the one hand, it is now admitted that planets and satellites have substantially migrated and been reorganised after their formation (see e.g. Tsiganis et al. 2005; Nesvorný & Morbidelli 2012; Lainey et al. 2009, 2020; Lari et al. 2020). Hence, trying to adjust formation scenarios to today’s locations of the planets and their satellites is probably inadequate. On the other hand, we assumed so far that at the time of Uranus’s tilting, the hypothetical ancient satellite was not affected by the attraction of additional close-in satellites. To lowest order, the gravitational attraction of inner moons is equivalent to increasing the planet’s J₂ (Tremaine et al. 2009). If we take the current satellites of Uranus into account at their current locations, the effective value of J₂ felt by the hypothetical distant satellite would be multiplied by a factor k ≈ 5.4. As a result, the Laplace radius is multiplied by k^1/5 ≈ 1.4 (see Eq. (2)) and r_M becomes 75R_eq instead of 53R_eq. This new value still does not seem too unrealistic, but it does not match any more that proposed by Mosqueira & Estrada (2003a,b). In order to avoid this increase in distance for the hypothetical ancient satellite, the current moons of Uranus should either have been formed after the tilting (see Sect. 4.3 below), or have been located closer to Uranus than they are today.

The total satellites-to-planet mass ratios for Jupiter, Saturn, and Uranus are today a few times 10⁻⁴. Canup & Ward (2006) argued that such similar mass ratios can be explained naturally if the satellites of all giant planets formed in a similar way within a circumplanetary disc of gas and dust, with multiple generations of satellites being formed and lost by migrating through the disc. Under the widely accepted assumption that Triton is a captured object, Rufu & Canup (2017) then showed that a primordial satellite system with mass ratio roughly equal to 10⁻⁴ or less is also expected to have existed around Neptune. This general picture seems to contradict the existence of a primordial satellite with mass m/M ≳ 10⁻³ around Uranus, as proposed here. However, even though the typical ratio of 10⁻⁴ was found by Canup & Ward (2006) to only weakly depend on the poorly known parameters involved, we have no guarantee that satellites around Jupiter, Saturn, and Uranus did form in roughly similar external conditions. Mass ratios larger than 10⁻³ can actually be obtained from the formulas of Canup & Ward (2006) by slightly tweaking the poorly known parameters. Moreover, the assumption that the regular satellites of all four giant planets have been shaped through the same physical processes is itself debated. Circumplanetary discs around the Solar System giant planets probably had much more diverse histories than initially envisioned by Canup & Ward (2006), and additional phenomena, as the truncation of the circumplanetary disc due to the planet’s magnetosphere, are expected to have played an important role (see e.g. Sasaki et al. 2010; Batygin & Morbidelli 2020).

Recently, Szulágyi et al. (2018) breathed new life into the possibility that the satellites of Uranus and Neptune formed early in a gaseous circumplanetary disc. The range of masses that they obtain in the framework of their population synthesis appears in line with our proposed ancient satellite. Yet, in this paper, the formation of the disc relies on a heat sink forcing the gas to not exceed a temperature of 100 K in the neighbourhood of the planet, which may not be realistic. Besides, the satellite population obtained strongly depends on the flux of mass injected within the disc during planetary formation, which is an unknown quantity. For these reasons, it seems to us that the mass distribution obtained by Szulágyi et al. (2018) cannot be used as an argument in favour of given mass values in our scenario: using a higher mass flux would increase the mass of the satellites obtained, and vice versa. What we conclude from Szulágyi et al. (2018) and previous works is that satellites as massive as what we need here can be formed routinely from circumplanetary discs that appear ‘realistic’ for Uranus and Neptune. As a general rule, very big satellites may also exist around giant planets (with masses as large as m/M ≈ 10⁻²), but they are probably not a generic outcome of formation processes within gaseous circumplanetary discs. Recent 3D hydrodynamic simulations seem to suggest that, unlike Jupiter and Saturn, planets similar to Uranus and Neptune do not develop significantly massive discs unless the equation of state of the gas is taken to be isothermal (see Fung et al. 2019). For this reason, the formation of a massive satellite around Uranus may have been unlikely.

Unusually big satellites may also be formed during the protoplanetary disc phase as captured coorbital proto-cores (Hansen 2019). However, those are expected to have large orbital inclinations, with most likely values near 180°, in contrast with the scenario proposed here in which the satellite starts close to its local (prograde) Laplace plane. Hence the formation scenario of Hansen (2019) could be invoked in this context only if it is associated with an efficient mechanism of inclination damping. If instead the satellite remains highly retrograde, then tidal dissipation within the planet would result in an inward (rather than outward) migration. As discussed in Sect. 4.2, the mechanism for tilting Uranus can work just as well, or even better, if the satellite migrates inwards, however there is no evidence showing that a powerful inward migration of retrograde satellites is possible. The capture of a large body during the planetesimal-driven planetary migration is also a possibility for Uranus to have acquired a big ancient satellite (see e.g. Agnor & Hamilton 2006; Li & Christou 2020); however, captured objects typically have large eccentricities and inclinations (Nesvorný et al. 2007), in contrast to the satellite needed here for the tilting of Uranus.

Fig. 18 Example of simulation that reproduces the current state of Uranus. We assume that the satellite is instantly removed when it goes below the Roche limit (red point). The mass of the satellite is m/M = 2:2 × 103 and it migrates with a timescale τ = 7:5 Gyr.

Fig. 19 Same as Fig. 18, but where the satellite crosses the Roche limit while having a prograde orbit with respect to the planet’s equator. The mass of the satellite is m/M = 2 × 10−3 and it migrates with a timescale τ = 3 Gyr.

4.2 Smallest ancient satellite allowing for Uranus to be tilted

Even though a mass ratio of m/M ≈ 1.7 × 10⁻³ does not appear absurdly high in the context of satellite formation theories, smaller satellites are generally thought to be more likely for a primordial formation in a gaseous circumplanetary disc. For comparison, the mass ratios for the Galilean satellites range from about 3 to 8 × 10⁻⁵, and the mass ratio for Titan is about 2 × 10⁻⁴. The possibility for Uranus and Neptune to have had a massive circumplanetary disc during the late stages of their formation is itself much debated (see e.g. Szulágyi et al. 2018; Reinhardt et al. 2020 and Sect. 4.1 above). As the formation of small satellites requires less material and imposes weaker constraints on their hypothetical gaseous birth disc, it appears natural to try to minimise as much as possible the mass required for the ancient satellite in our scenario.

Because of the characteristics of Uranus’s orbital motion, we can be assured from strong dynamical grounds that no tilting can be achieved for a single ancient satellite with mass smaller than m/M ≈ 4.4 × 10⁻⁴ (see Sect. 2). With this minimum mass, a satellite allows for a capture of Uranus in secular spin-orbit resonance with v₃ and adiabatic tilting up to large obliquity values. When the system reaches the unstable zone and the satellite goes wild, however, a mass of m/M ≈ 4.4 × 10⁻⁴ is insufficient to produce a strong back reaction on the planet, and therefore the satellite is generally destructed before the planet’s spin axis has time to diffuse much higher up (see Sect. 3). As a result, the planet’s obliquity is left fossilised to values that are still close to the boundary of the unstable zone (i.e. from about 75° to 80° using the parameters of Fig. 11), and far from Uranus’s current value. We can think of different ways to solve this problem and increase the success rate for small satellites.

The first and most obvious possibility would be that the ancient satellite was not alone around Uranus. As explained above, the addition of close-in satellites increases the effective J₂ of Uranus (see e.g. Tremaine et al. 2009). This increase pushes away the Laplace radius of our hypothetical ancient satellite, but it also decreases the mass that it needs to have in order for Uranus to reach a given secular spin-orbit resonance. If Uranus’s J₂ is increased by a factor k by a set of inner satellites, then r_M is enlarged by a factor k^1/5 and the mass needed to reach a given resonance is divided by k^2/5 (see equations in Sect. 2.1). It is difficult to judge a priori which satellite configuration would be most appropriate here, and how would inner satellites behave when their distant neighbour becomes unstable. To get a rough idea of the effect of inner satellites, we note that the current satellite system of Uranus would push the critical radius r_M to 75 R_eq and divide the mass needed for the ancient satellite by roughly a factor two. In particular, the strong resonance v₂ showing the highest success probabilities in our scenario would be reachable by an ancient satellite with mass m/M ≈ 10⁻³ instead of 2 × 10⁻³. This might increase the plausibility of our scenario. Yet, a set of inner satellites would probably not help to increase the success probabilities obtained in other resonances, because the issue does not come from the tilting phase (resonance v₃ is perfect for that) but from the destabilisation phase (during which small satellites do not produce a strong enough back reaction on the planet’s obliquity). A variant of this possibility would be to have a succession of several migrating satellites that reach the unstable region one after another: the first one would produce the resonance capture and tilting to about 80°, and the other ones would just reach the unstable region (with the obliquity being already high) and produce some chaotic obliquity diffusion before being disrupted. If they sum up, each of these limited obliquity changes may allow the current obliquity of Uranus to be reached. However, in any case, the exact behaviour of a set of satellites is hard to predict and one would need to explore these scenarios self-consistently using numerical simulations.

The second possibility would be to delay the instability of the satellite such that the obliquity can climb higher up inside the resonance before the destabilisation phase is triggered. As the instability originates from an eccentricity increase (see Tremaine et al. 2009; Saillenfest & Lari 2021), such a delay could be produced through tidal dissipation inside the satellite, whose main effect is to damp eccentricity. Yet, dissipation within the satellite should not be too strong either, otherwise the instability could be completely suppressed, or the migration of the satellite could be reversed before the instability has been properly triggered. As mentioned earlier, the instability is so strong that it is not expected to be suppressed for realistic dissipation parameters of the satellite (see Fig. 7 and the related discussion). However, the rate of energy dissipation within the satellite depends on its interior properties which can only be speculative at this stage. Because of the wide range of possible behaviours that can be obtained when varying the dissipation rate, we leave the exploration of this possibility for future works.

The third possibility would be to follow a secular spin-orbit resonance that is connected higher up to the unstable region, so that the instability is triggered at a larger value of obliquity. The different panels in Fig. 2 show that a given resonance can reach the unstable zone at different obliquity values, depending on the exact mass of the satellite. This change in the resonance shape is the reason why resonances in Fig. 12 show a higher success ratio for smaller a₀ and larger mass (i.e. in the lower end of the coloured bands). For resonance v₃, however, a non-negligible fraction of successful trajectories can be obtained only for a₀ ≲ 33 R_eq (see panels c and f of Fig. 12). Such small values of a₀ require a very fast tidal migration in order for the system to reach the unstable zone in less than the age of the Solar System. For even smaller initial semi-major axes, the migration rate of the satellite would need to be so large that adiabatic capture in resonance v₃ would be endangered (see Fig. 6). This third solution in order to increase the success rate of small satellites is therefore not conclusive. Alternatively, as shown in Fig. 2, very late instabilities would be produced for a satellite migrating inwards from an initial semi-major axis a₀ > r_M. In order for the satellite to migrate inwards, it should be retrograde with respect to the planet’s spin motion. Retrograde massive satellites with I_Q ≈ 180° can be formed through the mechanism described by Hansen (2019). However, it is not clear whether retrograde satellites could be able to trigger tidal dissipation mechanisms that are as efficient as those produced by prograde satellites (e.g. by a process analogous to that described by Fuller et al. 2016). This possibility is therefore quite speculative at this stage.

The fourth possibility would be to invoke an additional mechanism after the adiabatic drift in resonance v₃ in order to produce the final missing tilt Δε of about 20° and reach ε = 98°. The first idea that comes to mind is to take into account the effect of the satellite collision in the final obliquity value of the planet. The satellite destabilisation may also trigger a chain reaction involving other satellites and lead to several impacts. Assuming a perfect merger, the effect of a satellite collision into the planet can be estimated from the conservation of their total angular momentum. Indeed, the spin-axis motion of the planet and the destabilisation of its satellite are due to the secular effect of external perturbers (namely, the Sun and other planets); as such, external perturbations act on long timescales (at least tens of thousands of years; see Table 1 of Saillenfest & Lari 2021). During the very limited time interval of a collision, the planet and its satellite can therefore be considered as an isolated system. The details of computations are given in Appendix F. Considering that the ancient satellite is destabilised by the mechanism described above, the maximum obliquity change of the planet due to the impact of its satellite is, in radians,

$$\text{Δ}\epsilon \approx 35\frac{m}{M}.$$(12)

This means that in order to produce a 10° extra tilt, the ancient satellite must have at least a mass m/M ≈ 5 × 10⁻³, that is, bigger than Mercury. In order to produce 15°, it must have the size of Mars. These very large masses do not seem realistic in the context of satellite formation, and they clearly do not help to increase the success rate of small satellites. The effects of a chain reaction involving several satellites are difficult to predict a priori, but because of the conservation of total angular momentum, the total mass of the satellites would need to be quite large anyway, and in the same order of magnitude as in Eq. (12).

Alternatively, given the small size of the obliquity gap to be filled, one can think of a subsequent impact with a rogue minor planet, once the tilting over ε ≈ 80° has been achieved through the satellite migration. We consider here the scenario of a small late impact that only marginally alters the spin rate and obliquity (quite differently from giant impact scenarios mentioned in Sect. 1, for which the impact itself fully settles the spin rate and obliquity of the planet). Since collisions capable of a tilting over Δε ≈ 20° can still vary a planet’s spin rate by as much as a factor of two (Rogoszinski & Hamilton 2020), the most likely size of impactors must be determined with the condition that it should not significantly change the planet’s spin rate. We accomplish this by using the code developed by Rogoszinski & Hamilton (2020, 2021) to calculate a planet’s final spin state after collisions. The planet grows by summing the angular momenta of impactors and the planet, and the planet’s final tilt and spin are extracted from its resulting angular momentum vector. For simplicity, we assume that the impactors originate from within the Solar System, they travel on trajectories that are parallel to the plane of the planet’s orbit, and all their mass is absorbed upon impact. These approximations are sufficient because we assume small orbital inclinations, and less than 0.1 M_⊙ of debris is expected to be ejected from the system after an Earth-mass strike (Kegerreis et al. 2018, Kegerreis et al. 2019). The impactor hits at a random location on the planet’s surface and over all possible spin-axis precession phases, so we run this code for a half million iterations to generate probability distributions of the planet’s spin state. Finally, since the impactors are likely to travel on low-eccentricity orbits prograde to Uranus’s, we follow the approach of Hamilton & Burns (1994) and sample the relative speed of impactors between 0 and 0.4 times the Keplerian speed of a circular orbit at the distance of Uranus (6.8 km s⁻¹ ). The values obtained are much less than Uranus’s escape speed (21.4 km s⁻¹), so we include gravitational focusing (i.e. trajectories are aimed closer to the planet’s centre). Spin and obliquity must be handled simultaneously: we calculate the likelihood of generating Uranus’s current spin state by measuring the fraction of instances that are within both ±5° of Uranus current obliquity and ±10% of its current spin rate. Figure 20 shows the results obtained for initial obliquities between 0° and 90° and impactor masses from 0 to 1 M_⊕. The most likely impactor that can generate Uranus’s current spin state from an initial tilt of 85° has a mass of about 0.1 M_⊕. Again, this is close to the size of Mars. Even though pebble accretion models posit a surplus of Mars-sized cores in the formation region of the giant planets (see e.g. Levison et al. 2015), and especially beyond Saturn (Izidoro et al. 2015), these planetesimals are thought to have been accreted and scattered quickly during planetary formation. After billions of years of evolution, as required by our tilting mechanism, the structure of the Solar System was most likely very similar to what it is now. The most massive objects that may possibly impact Uranus today are trans-Neptunian dwarf planets, whose small masses lead to a negligible probability of producing the small missing tilt that we are looking for (Eris is 0.003 M_⊕; compare with Fig. 20). Hence, even though planetesimal impacts probably happened during the late stages of planetary formation and possibly produced a small initial obliquity excitation (which would not be in contradiction with our tilting scenario; see Fig. 6), impacts of external bodies do not allow us to decrease the mass of the satellite needed to tilt Uranus in our scenario.

In conclusion of this subsection, it seems difficult to reproduce the current spin state of Uranus if its ancient satellite is single and has a mass smaller than m/M ≈ 1.7 × 10⁻³. Taking into account a system of additional inner satellites resembling the current ones would allow this lower limit to be decreased by roughly a factor two. Below m/M ≈ 10⁻³, however, our simulations show that satellites hardly manage to excite the planet’s obliquity high enough. In order to increase further the success rate of small satellites, another possibility would be to include the effects of tidal energy dissipation within the satellites. This could help to decrease the minimum mass for a single satellite down to m/M ≈ 4.4 × 10⁻⁴ at most.

Fig. 20 Probability of reproducing Uranus’s current spin state as a function of the mass of a single impact strike. The planet is initialised with the current spin rate of Uranus, considered to be primordial. The mass of the impactor is given in Earth masses (M⊕).

4.3 Forming Uranus’s current satellites in a debris disc

Many previous works have been devoted to the formation of Uranus’s current satellites in the context of the giant impact hypothesis (see e.g. Morbidelli et al. 2012; Ida et al. 2020; Rufu & Canup 2022; Salmon & Canup 2022; Woo et al. 2022). If our new tilting scenario is correct, then it should also not be in contradiction with the existence of Uranus’s current satellite system. As discussed above, strong dynamical arguments show that Uranus cannot be tilted through the scenario considered here if its ancient satellite is single and has a mass smaller than m/M ≈ 4.4 × 10⁻⁴. This minimum mass is about four times the total mass of Uranus’s current satellite system. Assuming that the ancient satellite goes below the Roche limit during the final instability (see Sect. 3), then it would be torn apart into a large amount of debris. The remaining pieces of material would then collide with each other and rapidly reorganise into an equatorial disc confined inside the Roche limit (see Morbidelli et al. 2012; Hyodo et al. 2017).

Since the current regular satellites of Uranus seem to have been formed in such a tidal disc (Crida & Charnoz 2012; Hesselbrock & Minton 2019), it is tempting to link them to the disruption of our hypothetical ancient massive satellite. If the ancient satellite had mass m/M = 4.4 × 10⁻⁴ and was initially the only regular satellite of Uranus, then about 75% of its mass should have fallen onto the planet or been ejected, while the remaining 25% served as building blocks for the current satellite system of Uranus. These proportions are very similar to those obtained by Hesselbrock & Minton (2017) in their simulations of the formation of Phobos from a debris disc. Hesselbrock & Minton (2017) show that the precise formation efficiency of the new generation of satellites depends on the exact radius at which the ancient satellite is torn apart (which itself depends on its cohesive strength): if the ancient satellite is destructed closer to the planet, then the viscous spreading of the debris disc results in more material to be lost by falling onto the planet. Hence, in the context of the tilting mechanism described above, the minimum mass m/M ≈ 4.4 × 10⁻⁴ for the ancient satellite is compatible with the formation of all current regular satellites of Uranus from the debris disc, but more massive ancient satellites are not ruled out.

When the disc spreads beyond the fluid Roche limit, the pieces of debris aggregate into new satellites through the mechanism described by Crida & Charnoz (2012). For a disc mass m_disc, the characteristic lifetime of the disc T_disc = m_disc/ṁ_disc is

$$T_{\text{disc}}\approx 0.0425{\left(\frac{M}{m_{\text{disc}}}\right)}^2{T}_{\text{Roche}},$$

where T_Roche is the orbital period at the fluid Roche limit. Assuming that the disc contains just enough material to form all the current regular satellites of Uranus, one obtains a characteristic lifetime of about 4000 yr. This timescale is reduced to 200 yr if the disc has an initial mass m_disc/M ≈ 4.4 × 10⁻⁴. The formation of the new generation of satellites is therefore extremely fast at the beginning, and it slows down asymptotically as m_disc ∝ t^−1/2. When the disc mass has decreased below a given threshold, other phenomena can take over: Estrada & Durisen (2021) show that micro-meteoroid bombardment can slow down ring particles and force small-mass rings to drift inwards and eventually vanish relatively quickly. These short evolution timescales would explain why no massive disc remains around Uranus today, even in the hypothesis that the destruction of its ancient satellite (and the end of the tilting mechanism) occurred recently in the history of the Solar System, perhaps a few hundreds of million years ago. The current inner sparse ring system of Uranus may be a remnant of the old massive debris disc, and the innermost satellites of Uranus may be evolved residuals of the last generation of satellites formed from the disc. These innermost satellites are known to be unstable on short timescales (see e.g. Duncan & Lissauer 1997; French & Showalter 2012; French et al. 2015), and it is probable that they continuously form new ring particles and new small satellites through mutual collisions and reaccretion (Tiscareno et al. 2013).

After their formation in the debris disc, the regular satellites of Uranus must have migrated up to their current location. Oberon is the farthest regular satellite of Uranus and it is located today at a semi-major axis of a ≈ 23 R_eq. This distance is to be compared to the synchronous orbit around Uranus, which lies somewhat below 3.5 R_eq (the migration from the fluid Roche limit at 2.5 R_eq to the synchronous orbit is produced through Lindblad torques from the disc while it is massive enough; see Hesselbrock & Minton 2019).

The current migration rate of Uranus’s satellites is essentially unconstrained today from direct observations (Lainey 2016), and even though theoretical arguments can be used to give bounds to the unknown parameters (Tittemore & Wisdom 1990; Cuk et al. 2020), our lack of knowledge about the dissipation mechanisms at play in the interior of Uranus leaves room for a large range of possible scenarios. For instance, previous arguments based on mean-motion resonance encounters may be invalidated if satellites are undergoing a mechanism similar to the “tidal resonance locking” of Fuller et al. (2016), for which the effective quality factor Q of the planet evolves over time in such a way that semi-major axis ratios are constant (Crida 2020). In fact, as shown by Lainey et al. (2020), the migration history of satellites around giant planets can result to be spectacularly different from what one would expect from constant-Q models. Based on the observed migration rate of Saturn’s satellites, and assuming that the tidal mechanism of Fuller et al. (2016) was not triggered long after their formation, then Titan, for instance, would have migrated across a radial distance of 10 radii of Saturn in 4 Gyr, that is, more than 600 000 km. In contrast, if we suppose that Oberon was formed from a debris disc at the Roche limit of Uranus, then it would need to have migrated over about 20 radii of Uranus before today, which are about 500 000 km. These large numbers may seem extraordinary high when viewed in the context of previous works. For instance, Morbidelli et al. (2012) argued that a tidal disc around Uranus cannot have been the birth place of its current satellites because they are too far away⁶. However, these distances are in line with the measurement of the current orbital expansion of Saturn’s satellites (Lainey et al. 2020). In any case, if we assume that Uranus had a distant ancient satellite subject to substantial tidal migration (which is our basic hypothesis throughout this paper), then some mechanism of efficient energy dissipation must exist in the interior of Uranus. Regardless of the precise nature of this mechanism, the dissipation efficiency needed to move Oberon to its current location (measured by the equivalent constant parameter k₂/Q; see Sect. 2.3) is smaller than that needed to move the distant satellite over 10 R_eq, so if we assume that the latter is possible, then the former should follow.

Yet, unless the migration rate of Oberon is much larger than that of Titan (which cannot be firmly ruled out yet, but may seem unlikely), we cannot expect it to have migrated up its current position in less than a couple billion years. This means that the end of the tilting mechanism of Uranus and the destruction of its ancient satellite cannot be arbitrarily recent if one wants to explain the formation of its current satellites in the debris disc formed. For instance, one could imagine that the migration of the ancient satellite and tilting of Uranus took 2 Gyr; then it was followed by the formation of its current satellites in the debris disc; and then the migration of its current satellites from the synchronous orbit to their current location took an extra 2 Gyr. Of course, the exact timing of these two phases is highly speculative.

The creation of a tidal disc that would be suitable for the formation of the current satellites of Uranus is itself not straightforward. Section 3.4 reveals that, for dynamical reasons, the massive ancient satellite is much more likely to collide into Uranus through the retrograde keyhole, whereas the current satellites of Uranus are prograde with respect to its spin axis. Yet, several factors could lead to the formation of a prograde tidal disc even if the satellite collides on a retrograde orbit. For instance, the satellite is not expected to be destructed at once but to form a torus of debris, and all pieces of debris still orbit in the highly unstable region before they are damped to the equator plane, so whether the final tidal disc should be prograde or retrograde is perhaps not as obvious as it may appear. Impacts of debris into the planet could also lead to the production of prograde ejecta.

Instead of being directly torn apart by tides, the destabilisation of the ancient satellite may also have triggered a collision cascade in a pre-existing system of inner moons. In this case, a large portion of the total satellite mass would need to be ejected or collide into the planet anyway, in order to reproduce the current mass of Uranus’s system. Furthermore, dramatic satellite collisions beyond the Roche limit are not expected to produce a tidal disc, but to lead to a rapid reaccretion of debris into new moons without substantial spreading (Hyodo & Charnoz 2017). This means that if satellite-satellite collisions occurred far from the Roche limit, then the current satellites of Uranus would need to be the remnants of the collision cascade, which seems to be in contradiction with their spacing and mass distribution that closely match those expected from formation in a tidal disc (Crida & Charnoz 2012). Given the large eccentricity and inclination of the massive ancient satellite when it is destabilised, one would also expect a substantial orbital excitation from such large satellite-satellite collisions. Even though eccentricities can be quickly damped through tidal dissipation, this is not the case of inclinations, which have a much longer damping timescale. Today, only Miranda has a substantial orbital inclination, and its origin is well explained as a consequence of a past 5:3 resonance between Ariel and Umbriel (Cuk et al. 2020). Previous studies obtained similar inclination excitations via the evolution of Miranda through the 3:1 mean-motion resonance with Umbriel (Tittemore & Wisdom 1989; Malhotra & Dermott 1990; Verheylewegen et al. 2013). Depending on the actual migration rate of the satellites, this inclination increase may be quite recent, and trying to link it to the destabilisation of our big ancient moon would be doubtful. A way to protect the current satellites of Uranus from the destabilising action of their distant neighbour would be to place them closer to the planet at the time of the tilting, because this would lock them more tightly to Uranus’s equator. Hence, independently of whether the current satellites of Uranus formed before or after the tilting, invoking a substantial outward migration for them seems inescapable.

5 Summary and conclusion

Giant planets are thought to form in a protoplanetary disc made of gas and dust. During their last formation stages, gas accretion dominates their final spin state and gives them a primordial near-zero obliquity. The similarities between the four giant planets of the Solar System – their spin rates in particular − suggest such a common formation mechanism, and may disfavour stochastic processes, such as giant impacts, as being a main formation ingredient. In this view, the small to moderate obliquity values of Jupiter (3°), Saturn (27°), and Neptune (30°) can be explained by post-formation events (see e.g. Ward & Hamilton 2004; Hamilton & Ward 2004; Ward & Canup 2006; Vokrouhlický & Nesvorný 2015; Rogoszinski & Hamilton 2020; Saillenfest et al. 2021a). Explaining the extreme obliquity of Uranus (98^°) is more challenging. Tilting Uranus without affecting its spin rate much would require a slow process involving its spin-axis precession motion. Such a process hardly fits the timespan offered by the different stages of the early Solar System evolution (e.g. dissipation of the gas disc or late planetary migration). This is why previous studies investigating the tilting of Uranus as a post-formation event faced a timescale issue (Boué & Laskar 2010; Quillen et al. 2018; Rogoszinski & Hamilton 2021).

In contrast, the tidal migration of satellites is an everlasting process that can slowly change the orbit of natural satellites over billions of years. Massive satellites are known to affect the spin-axis motion of planets (see e.g. Tremaine 1991; French et al. 1993; Boué & Laskar 2006), and over the lifetime of the Solar System the migration of the main satellites of Jupiter and Saturn can potentially produce dramatic obliquity changes to their host planet (Saillenfest et al. 2020, 2021b). In this article, we examined whether this mechanism could also apply to Uranus and possibly explain its extreme obliquity.

The mechanism considered here involves secular spin-orbit resonances, that is, resonances between the precession of the spin axis of the planet and some harmonics appearing in its orbital precession. It is made of two different phases as satellites migrate: (i) a capture in resonance and steady obliquity increase, and (ii) a violent destabilisation of the satellite when the obliquity exceeds 70° to 80° (Saillenfest & Lari 2021). Both phases can be used to determine the properties of the satellite that would be needed to reproduce the current spin state of Uranus.

Phase 1 As we are studying a very slow process and the orbit of Uranus has been stable for billions of years, we need to consider the current orbital dynamics of Uranus as a forcing to its spin-axis motion. This puts strong constraints on the possible secular spin-orbit resonances involved. The orbital forcing term that is closest to Uranus’s spin-axis precession rate has frequency ν₁₅ = 𝒢₇ − 𝒢₈ + s₇ (where 𝒢_j and s_j are the apsidal and nodal precession modes of the Solar System planets numbered from 1 for Mercury to 8 for Neptune). In order to produce a resonance with ν₁₅, Uranus must have had a satellite with minimum mass m/M = 3.5 × 10⁻⁴ (assuming that Uranus had a single satellite at that time); however, this resonance is too weak to allow for an adiabatic capture and tilting in less than the age of the Solar System. The second closest forcing term has frequency v₃ = s₈. Triggering a resonance with v₃ and tilting Uranus up to a large obliquity would require a satellite with minimum mass m/M = 4.4 × 10⁻⁴, that is, m ≈ 4 × 10²² kg (smaller than Jupiter’s moon Europa). Simulations show that there is a 100% probability for Uranus to be captured and tilted in this resonance in a wide interval of satellite migration rates and for primordial obliquities ε between 0° and 20°. In order to complete the tilting phase, the satellite must migrate over a range of about 10 radii of Uranus during its history; covering this distance in less than the age of the Solar System would require the satellite to migrate at least 6 cm yr⁻¹ on average. This migration rate is orders of magnitude higher than what can be inferred from historical models of tidal dissipation within Uranus, but quite comparable to the 11 cm yr⁻¹ measured for Titan (see Lainey et al. 2020). A mechanism similar to that described by Fuller et al. (2016) could be a viable explanation for a fast satellite expansion around Uranus.

A mass of m/M = 4.4 × 10⁻⁴ is roughly four times the total mass of the current satellites of Uranus; this means that the satellite involved does not exist anymore and it must have been destructed during the second phase of the evolution. More massive ancient satellites allow for other resonances to be reached. The strongest resonance is ν₂ = s₇ which can be reached for a satellite mass m/M = 2.2 × 10⁻³. Other promising resonance candidates are ν₇ = 𝒢₅ − 𝒢₇ + s₇ and ν₅ = −𝒢₅ + 𝒢₆ + s₆, which can be reached for masses m/M = 1.3 × 10⁻³ and 1.7 × 10⁻³, respectively (i.e. between the masses of Jupiter’s moons Callisto and Ganymede). These mass estimates hold for a single ancient satellite. If Uranus had a system of additional inner moons at that time, then the mass needed would be reduced. For instance, if we place the hypothetical ancient satellite together with Uranus’s current moon system, all mass estimates given above should be divided by roughly a factor two.

Phase 2 We explored the destabilisation phase using a coupled secular model including, altogether, the satellite’s orbit and the planet’s spin-axis dynamics. When the system reaches the unstable zone, the satellite’s orbit can be excited very quickly and reach almost any value of eccentricity and inclination. If the satellite is massive enough, its wild orbital variations produce a back reaction on the planet’s spin axis: the planet sweeps over various secular spin-orbit resonances and undergoes erratic obliquity kicks. These kicks allow the planet to go much deeper in the unstable region and reach obliquity values similar to Uranus’s today. Sooner or later, the satellite reaches the Roche limit (or collides into the planet), which puts an end to the chaotic dynamics and fossilises the planet’s spin in its last state. Satellite collisions occur through the action of the secular resonance 2ω + Ω located at equatorial inclination I_Q ≈ 56°, and to its symmetric counterpart 2ω − Ω at I_Q ≈ 124°. This resonance is well known for artificial satellites (see e.g. Daquin et al. 2022). Since it vanishes for a planet obliquity ε = 90°, satellite collisions can only occur when the planet’s obliquity is substantially larger than 90° (as Uranus today), or substantially smaller than 90°. Interestingly, the current state of Uranus is near the maximum of our histograms in the ε > 90° side. This may be in favour of the mechanism proposed here.

Larger resonances and more massive satellites produce stronger obliquity kicks. Our simulations show that for m/M = 4.4 × 10⁻⁴, the satellite generally collides into the planet before the obliquity reaches 80°. Hence, it seems hard to reproduce the current state of Uranus through a capture in resonance ν₃, even if this resonance demonstrates a great capture and tilting efficiency. We discussed other phenomena that may produce the missing extra tilt for such small ancient satellites (e.g. a subsequent collision with a minor planet) but found none to be satisfactory. More massive satellites are much more successful in reproducing the current spin state of Uranus. If the satellite has a mass m/M between 2 and 3 × 10⁻³, the probability of success through a capture in resonance ν₂ can be as large as 86%. For slightly smaller masses (m/M ≈ 1.5 to 2 × 10⁻³), resonances ν₅ and ν₇ also demonstrate large success ratios, which reach 50% in some ranges of parameters. As before, these mass estimates can be reduced if we assume that Uranus already had a set of additional inner moons at that time; their survival, however, is far from being guaranteed.

Resonance ν₂ is special because it is a secular spin-orbit resonance with Uranus’s own nodal precession mode s₇. Due to its large width, it produces large-amplitude oscillations of the obliquity that can destabilise the satellite very early on, after only a short passage through Phase 1. Moreover, resonance ν₂ produces very strong kicks that can quickly pump the obliquity up to Uranus’s current value. For this reason, our simulations show that as long as the satellite is initialised at the right distance (a₀ ≈ 50 R_eq), the tilting mechanism via resonance ν₂ can operate down to very small migration ranges (only a few R_eq, instead of 10 R_eq in the general case). When studying the formation of the main satellites of Jupiter, Saturn, and Uranus in gaseous circumplanetary discs, Mosqueira & Estrada (2003a,b) found that a distance of about 50 R_eq is precisely where one would expect a massive satellite to form around Uranus. The current main satellites of Uranus (Miranda, Ariel, Umbriel, Titania, and Oberon) are much further in and they may seem peculiarly small as compared to other satellites of the giant planets. From these argument, one may argue at least that forming a satellite with mass ratio m/M ≈ 2 × 10⁻³ at a distance of a₀ ≈ 50 R_eq does not seem too unreasonable in the context of satellite formation theories.

Any successful tilting scenario must not only reproduce the current spin state of Uranus, but also allow for the existence of its observed satellite system. Today’s main satellites of Uranus have prograde orbits, and they show signs of having been formed in a tidal disc (Crida & Charnoz 2012). We have discussed the possibility that they formed from the debris disc produced by the destruction of Uranus’s hypothetical ancient satellite, after the latter went below the Roche limit. The formation of the new generation of satellites is expected to be very fast compared to other dynamical processes at play. Then, the newborn satellites would need to migrate outwards from the synchronous orbit (a ≈ 3.4 R_eq) up to their current distance (Oberon is at a ≈ 23 R_eq). Such a large migration range appears unrealistic in the context of traditional tidal models, but it could be achievable on a billion-year timescale if a mechanism similar to the ‘tidal lock’ recently proposed for Titan is also at play for the satellites of Uranus (Fuller et al. 2016; Lainey et al. 2020).Yet, if this is the case, the satellites would have no particular reason to follow the peculiar distribution described by Crida & Charnoz (2012) as evidence of their origin in a tidal disc. The possible relationship between our proposed ancient satellite and the current ones is therefore not a simple question.

Moreover, in our subset of simulations that successfully reproduce the spin state of Uranus, the disruption of the ancient satellite almost always occur when the latter is retrograde (I_Q ≈ 124°). Therefore, our simplified secular model does not provide a straightforward way to produce a prograde tidal disc from its debris, and one would need to model the disruption process in more detail to assess the properties of the disc – including tidal dissipation within the satellite, orbital and collisional evolution of the debris, and the possibility of additional ejecta from the planet.

These issues are solved if one considers that the current satellites of Uranus are primordial and formed in the same gaseous disc as our hypothetical ancient satellite (i.e. following the scenario of Mosqueira & Estrada 2003a,b). It is not clear, however, whether the current satellites would have survived during the violent unstable phase of their massive neighbour, and whether they would have kept their orbits as unexcited as they are today.

Eccentricities can be subsequently damped by tidal dissipation, but inclinations have much longer damping timescales. A possible solution to this problem would be to consider that the current satellites were closer to Uranus during the tilting and more tightly locked to its equatorial plane, and that they substantially migrated afterwards. As these scenarios involve several moons, one would need to investigate them using N-body simulations to get definite answers.

Our results are based on numerical explorations performed with secular models. Being fast and valid for arbitrary values of eccentricity and inclination, these models allow us to explore a vast range of parameters. During the first phase of the evolution, we are guaranteed that they give a qualitatively accurate picture. During the second phase, on the contrary, it may seem questionable to average the dynamics over the orbit of the satellite even when its eccentricity grows close to 1 and its precession becomes very fast. More work is now needed to determine the exact ultimate fate of the satellite and the properties of the tidal disc that it would create. In this extreme dynamical regime, unaveraged numerical integrations are probably required.

We have not investigated the effects of tidal dissipation within the satellite yet. Once again, tidal dissipation would affect the second phase of the evolution, when the satellite’s eccentricity strongly increases. The exact outcome of the destabilisation process is hard to predict a priori because it involves three competing effects: tidal dissipation within the planet (which produces the outward satellite migration), dynamical destabilisation (which increases the eccentricity), and tidal dissipation within the satellite (which circularises the orbit and produces an inward migration). The net effect on the satellite strongly depends on the relative timescales of these three competing effects, which are related to the specific mechanisms of energy dissipation at play. As dissipation processes inside ice giant planets are essentially unknown today, and since the internal composition of Uranus’s ancient satellite is highly speculative, we did not try to obtain a definite answer here. However, our results can serve as a robust starting point for future experiments. According to the level of tidal dissipation within the satellite, we can expect three different behaviours: (i) If dissipation is exceptionally strong, the destabilisation of the satellite would be inhibited and our scenario would fail to reproduce the spin state of Uranus; (ii) if dissipation is moderate, the destabilisation may be delayed but not suppressed, and this could increase the fraction of successful runs for a single small satellite, down to a mass m/M ≈ 4.4 × 10⁻⁴; and (iii) if dissipation is small, the picture outlined here would remain unchanged.

In this article, we have settled the basic ingredients that would allow for the current state of Uranus to be reproduced via the migration of an ancient satellite. This scenario can be refined. One promising variant would be to invoke two satellites: a big distant one (as that described here), plus a small one located much closer to Uranus. In this configuration, the inner satellite would be much more prone to tidal migration than the outer one; and while the inner satellite migrates outwards over a few R_eq, the outer one would feel a change in the effective J₂ parameter of Uranus, which would increase its characteristic radius r_M. As a result, the whole tilting mechanism would be triggered all the same, even though the distant big satellite would not need to migrate at all. This variant scenario exemplifies how it would be possible to relax constraints both on the large energy dissipation required in the interior of Uranus (as the migration of close-in satellites requires much less energy dissipation than the migration of distant ones), and on the large mass needed for the ancient satellite (as the combined action of several moderately massive satellites can produce the same effect). The exploration of this and other possible variant scenarios are left for future works.

Despite the unanswered questions mentioned above, this new picture for the tilting of Uranus appears quite promising to us. To our knowledge, this is the first time that a single mechanism is able to both tilt Uranus (phase 1) and fossilise its spin axis in its final state (phase 2) without invoking a giant impact or other external phenomena. The distribution of our successful runs peaks at Uranus’s location, which appears as a natural outcome of the dynamics. This picture also seems appealing as a generic phenomenon: Jupiter today is about to begin the tilting phase (Saillenfest et al. 2020), Saturn may be halfway in (Saillenfest et al. 2021b), and Uranus would have completed the final stage, with the destruction of its satellite.

Confronting this picture to observations is not an easy task. Part of the answer would be given by a measure of the tidal migration rate of Uranus’s current satellites. A high rate would indicate that they formed substantially closer to Uranus: this would give the possibility that they formed from the debris of the ancient satellite, or that they were protected from its wild destabilisation phase.

Acknowledgements

We thank the anonymous referee for his/her thorough analysis of our manuscript and inspiring comments about this new tilting mechanism.

Appendix A Tilting Uranus during the planetary migration

If we restrict the precession of the planet’s orbital plane to a single dominant harmonic with amplitude S₀ = sin(I₀/2) and angle ϕ₀ evolving at frequency v₀ (see Eq. 6), the Hamiltonian function governing the secular spin-axis dynamics of the planet is simplified into

$$\begin{array}{lll}\mathscr{H}\hfill & =\hfill & -\frac{\alpha }{2}\frac{X^2}{{\left(1-e_{\odot }^2\right)}^{3/2}}+v_0{\text{Φ}}_0\hfill \\ \hfill & +\hfill & v_0\left(1-\cos I_0\right)X-v_0\sin I_0\sqrt{1-X^2}\cos \left(\psi +{\varphi }_0\right),\hfill \end{array}$$(A.1)

where Φ₀ is the momentum conjugate to ϕ₀, added to obtain an autonomous Hamiltonian system. In this expression, α is the precession constant, e_⊙ is the eccentricity of the planet, X is the cosine of obliquity, and ψ is the precession angle (see e.g. Saillenfest et al. 2019b). Passing to the resonant canonical coordinates X and −σ = −ψ − ϕ₀ and dropping unnecessary constant parts, the Hamiltonian function can be rewritten

$$\mathscr{H}=-\frac{1}{2}{X}^2+\gamma X+\beta \sqrt{1-X^2}\cos \sigma ,$$(A.2)

using the rescaled timescale dτ = pdt, where

$$\begin{array}{ll}\gamma =-\frac{v_0\cos I_0}{p},\hfill & \beta =-\frac{v_0\sin I_0}{p},\hfill \end{array}$$(A.3)

and $p=\alpha {\left(1-e_{\odot }^2\right)}^{-3/2}$ is the characteristic spin-axis precession frequency of the planet. Equation (A.2) is the Hamiltonian of Colombo’s top problem, which has been widely studied (Colombo 1966; Henrard & Murigande 1987; Saillenfest et al. 2019b; Haponiak et al. 2020; Su & Lai 2020, 2022). Its equilibrium points are generally called the ‘Cassini states’, following Peale (1969).

Since mutual planetary perturbations dominantly produce a retrograde nodal precession, we have v₀ < 0, so that γ and β are both positive with these notations. Under the hypothesis that Uranus’s obliquity is small initially, the fastest tilting is obtained when Uranus acquires its current obliquity during a single cycle of σ (‘resonance kick’). Assuming that the planet’s orbit is stable enough to consider that γ and β do not vary much over this cycle, the existence of such an extreme trajectory puts strong constraints on the parameters. The existence of a trajectory connecting ε = 0° to more than ε = 90° first requires that

$$\beta ⩾\left|\gamma -1/2\right|,$$(A.4)

which is the condition for which a trajectory passing at X = 0 has the same Hamiltonian value has a trajectory passing at X = 1 (see Eq. (A.2)). Moreover, in order for these two points to belong to the same trajectory, the separatrix surrounding Cassini state 2 must either contain the north pole of the sphere or be inexistent. As shown by Saillenfest et al. (2019b), this condition can be written

$$8{\beta }^2⩾1-20\gamma -8{\gamma }^2+{\left(1+8\gamma \right)}^{3/2}.$$(A.5)

The two curves defined by Eqs. (A.4) and (A.5) connect at γ = γ_crit, where ${\gamma }_{\text{crit}}=\left(1+\sqrt{2}\right)/4$. For γ < γ_crit, Eq. (A.5) is the more stringent condition, whereas for γ > γ_crit, Eq. (A.4) is the more stringent condition. According to Eq. (A.3), γ and β are linked through β = γ tan I₀. Therefore, the two previous conditions can be formulated as a minimum value for I₀ instead of a minimum value for β. The minimum inclination of the planet allowing for its tilting from 0° to more than 90° as a function of γ is shown in Fig. A.1. The global minimum of this function is located at γ = γ_crit, and its value is given by tan $I_{\min }=3-2\sqrt{2},$ that is, I_min ≈ 9.74°.

Fig. A.1 Minimum orbital inclination of a planet such that its obliquity can evolve from 0° to more than 90° over a single libration in a secular spin-orbit resonance. Conditions 1 and 2 are given in Eqs. (A.4) and (A.5), respectively. The minimum inclination is reached at ${\gamma }_{\text{crit}}=\left(1+\sqrt{2}\right)/4$ and its value is given by tan $I_{\min }=3-2\sqrt{2}$.

The phase portrait of the Hamiltonian at the critical point is shown in Fig. A.2. We see that it corresponds to the case where the separatrix enclosing Cassini state 2 goes exactly from ε = 0° to ε = 90°.

The problem is the same if Uranus’s spin axis resonates with its own orbital precession (i.e. v₀ = s₇) or with a harmonic stemming from an other planet (e.g. Neptune, as investigated by Rogoszinski & Hamilton 2021, that is, v₀ = s₈): this harmonic must anyway fulfil the conditions described above, which means that the planet must have a large inclination in order to allow for a fast tilting. Imposing small inclination values, instead, means that Uranus cannot acquire a large obliquity during a single libration, and that we need an additional drift of the resonance centre (i.e. a variation of γ and β). As discussed in Sect. 1, this alternative strongly increases the tilting timescale, because the drift must be adiabatic enough in order for Uranus to remain inside the resonance and gradually follow its drift.

Appendix B Orbital solution for Uranus

The secular orbital solution of Laskar (1990) is obtained by multiplying the normalised proper modes $z_i^{•}$ and ${\zeta }_i^{•}$ (Tables VI and VII of Laskar 1990) by the matrix $\tilde{S}$ corresponding to the linear part of the solution (Table V of Laskar 1990). In the series obtained, the terms with the same combination of frequencies are then merged together, resulting in 56 terms in eccentricity and 60 terms in inclination. This forms the secular part of the orbital solution of Uranus, which is what is required by our averaged models.

The orbital solution is expressed in the variables z and ζ as described in Eqs. (6) and (7). In Tables B.1 and B.2, we give all terms of the solution in the J2000 ecliptic and equinox reference frame.

Fig. A.2 Phase portrait of Colombo’s top Hamiltonian when the planet has the minimum orbital inclination required for going from an obliquity ε = 0° to ε > 90° over a single libration. The top and bottom panels show the same phase portrait using two sets of coordinates. In the bottom panel, the obliquity ε is the tilt from the z-axis and the resonance angle σ is the polar angle measured in the xy-plane. The parameters are $\gamma =\left(\sqrt{2}+1\right)/4$ and $\beta =\left(\sqrt{2}-1\right)/4$, respectively. The Cassini states are labelled and showed with coloured dots. The thick black level is the separatrix of the resonance.

Because of the chaotic dynamics of the Solar System (Laskar 1989, 1990), the fundamental frequencies related to the terrestrial planets (e.g. s₁, s₂, and γ appearing in Table 1) could vary noticeably over billions of years (Hoang et al. 2021). However, they only marginally contribute to Uranus’s orbital solution, and none of them takes part in the resonances studied in this article. We can therefore safely consider that this orbital solution is valid on a billion-year timescale, at least in a qualitative point of view.

Appendix C Coupled model for the satellite and the spin axis of its host planet

In Sect. 3, we investigate the effect of the satellite destabilisation on the coupled dynamics of the satellite’s orbit and the planet’s spin axis. This investigation requires a self-consistent coupled model that also incorporates the orbital variations of the planet due to perturbations from all other Solar System planets, in order for the secular spin-orbit resonances to be present. In this section, we describe the model that we use for this purpose.

Appendix C.1 Unaveraged Hamiltonian function

Our setting is similar to that of Correia et al. (2011). We consider three bodies with masses m_i and positions x_i measured in an inertial reference frame. The momenta conjugate to the positions x_i are X_i = m_iẋ_i. We use the index 0 for the planet, index 1 for the satellite, and index 2 for the star. The planet is assumed to be an extended rigid body. The position x₀ of the planet is that of its centre of mass. We write (I, J, K) the basis vectors associated with the principal axes of inertia of the planet; the unitary vectors (I, J, K) are attached to the rigid planet and rotate with it. In this system of coordinates, the inertia tensor I of the planet is diagonal and writes:

$$I={\left(\begin{array}{lll}A\hfill & 0\hfill & 0\hfill \\ 0\hfill & B\hfill & 0\hfill \\ 0\hfill & 0\hfill & C\hfill \end{array}\right)}_{\left(I,J,K\right)}\text{\hspace{0.17em}\hspace{0.17em}where\hspace{0.17em}}A⩽B⩽C.$$(C.1)

Following Boué & Laskar (2006), the Hamiltonian function describing the dynamics of the system, as built from its total energy, can be written

$$\mathscr{H}={\mathscr{H}}_{\text{E}}+{\mathscr{H}}_{\text{N}}+{\mathscr{H}}_1,$$(C.2)

where ℋ_E corresponds to the Eulerian free rigid rotation of the planet, ℋ_N describes the Newtonian attraction of three mass points, and ℋ₁ contains the interactions of the non-spherical component of the planet with the satellite and the star mass points. The Eulerian part ℋ_E is

$${\mathscr{H}}_{\text{E}}=\frac{1}{2}{G}^{\text{T}}{ℐ}^{-1}G=\frac{{\left(G\cdot I\right)}^2}{2A}+\frac{{\left(G\cdot J\right)}^2}{2B}+\frac{{\left(G\cdot K\right)}^2}{2C},$$(C.3)

where G is the rotational angular momentum of the planet. The Newtonian part ℋ_N can be written

$${\mathscr{H}}_{\text{N}}={\displaystyle \sum _{i=0}^2\frac{{\Vert X_i\Vert }^2}{2m_i}-{\displaystyle \sum _{0⩽i⩽j⩽2}\frac{\mathcal{G}{m}_i{m}_j}{\Vert x_i-x_j\Vert },}}$$(C.4)

where 𝒢 is the gravitational constant. Because of the hierarchical nature of the system, we use the following set of Jacobi coordinates:

$$\{\begin{array}{l}{r}_0=\frac{m_0{x}_0+m_1{x}_1+m_2{x}_2}{m_0+m_1+m_2},\hfill \\ r_1=x_1-x_0,\hfill \\ r_2=x_2-\frac{m_0{x}_0+m_1{x}_1}{m_0+m_1},\hfill \end{array}$$(C.5)

with conjugate momenta

$$\{\begin{array}{l}{p}_0={\beta }_0{\dot{r}}_0=X_0+X_1+X_2,\hfill \\ p_1={\beta }_1{\dot{r}}_1=X_1-m_1\frac{X_0+X_1}{m_0+m_1},\hfill \\ p_2={\beta }_2{\dot{r}}_2=X_2-m_2\frac{X_0+X_1+X_2}{m_0+m_1+m_2},\hfill \end{array}$$(C.6)

where the reduced masses are defined as

$$\{\begin{array}{l}{\beta }_0=m_0+m_1+m_2,\hfill \\ {\beta }_1=\frac{m_0{m}_1}{m_0+m_1},\hfill \\ {\beta }_2=\frac{\left(m_0+m_1\right)m_2}{m_0+m_1+m_2}.\hfill \end{array}$$(C.7)

The vector r₀ describes the position of the barycentre of the whole system in the inertial reference frame, the vector r₁ is the position of the satellite with respect to the planet, and the vector r₂ is the position of the star with respect to the barycentre of the planet and its satellite. In this system of coordinates, the Newtonian part of the Hamiltonian function can be rewritten

$${\mathscr{H}}_{\text{N}}=\frac{{\Vert {\text{p}}_0\Vert }^2}{2{\beta }_0}+{\mathscr{H}}_{\text{K}}+{\mathscr{H}}_{\text{M}},$$(C.8)

where ℋ_K is the Hamiltonian describing two decoupled Keplerian motions:

$${\mathscr{H}}_{\text{N}}={\displaystyle \sum _{i=1}^2\left(\frac{{\Vert p_i\Vert }^2}{2{\beta }_i}-\frac{{\mu }_i{\beta }_i}{\Vert r_i\Vert }\right)}=-\frac{{\mu }_1{\beta }_1}{2a_1}-\frac{{\mu }_2{\beta }_2}{2a_2},$$(C.9)

in which we define µ₁ = 𝒢(m₀ + m₁) and µ₂ = 𝒢(m₀ + m₁ + m₂), and a_i are the associated semi-major axes. The Hamiltonian function ℋ_M contains the mutual Newtonian perturbations:

$${\mathscr{H}}_{\text{M}}=\left(\frac{\mathcal{G}\left(m_0+m_1\right)m_2}{\Vert r_2\Vert }-\frac{\mathcal{G}{m}_0{m}_2}{\Vert x_2-x_0\Vert }\right)-\frac{\mathcal{G}{m}_1{m}_2}{\Vert x_2-x_1\Vert }.$$(C.10)

We take the hierarchical nature of system into account by noting that

$$\begin{array}{lll}{x}_2-x_0\hfill & =\hfill & \frac{m_1}{m_0+m_1}{r}_1+r_2,\hfill \\ x_2-x_1\hfill & =\hfill & \frac{m_0}{m_0+m_1}{r}_1+r_2,\hfill \end{array}$$(C.11)

where ∥r₁∥ ≪ ∥|r₂∥, and by performing a multipolar expansion of ℋ_M using the Legendre polynomials. This gives

$${\mathscr{H}}_{\text{M}}=\frac{\mathcal{G}{\beta }_1{m}_2}{2r_2^3}\left(r_1^2-3\frac{\left(r_1\cdot r_2\right)}{r_2^2}\right)+O\left({\left[\frac{r_1}{r_2}\right]}^3\right),$$(C.12)

in which we write r_i ≡ ∥r_i∥.

The interaction part ℋ_I is obtained from the gravitational potential of a mass element within the planet interacting with the outer two point-mass bodies. Writing y the position of the mass element measured with respect to the planet’s barycentre and y_i ≡ x_i − x₀, we again take advantage of the hierarchical nature of the system by noting that ∥y∥ ≪ ∥y_i∥ ∀i = 1,2 and that the planet is almost spherical. By using a multipolar development of 1/∥y_i - y∥ and integrating over the planet’s volume, we get

$$\begin{array}{lll}{\mathscr{H}}_1\hfill & =\hfill & {\displaystyle \sum _{i=1}^2[\frac{\mathcal{G}{m}_i}{2{\Vert y_i\Vert }^5}[\left(2A-B-C\right){\left(y_i\cdot I\right)}^2+\left(2B-A-C\right){\left(y_i\cdot J\right)}^2}\hfill \\ \hfill & +\hfill & \left(2C-A-B\right){\left(y_i\cdot K\right)}^2]+O\left({\left[\frac{R_{\text{eq}}}{\Vert y_i\Vert }\right]}^3\right)],\hfill \end{array}$$(C.13)

where R_eq is the equatorial (i.e. largest) radius of the planet. Replacing y₁ and y₂ by their expressions as a function of r₁ and r₂, we finally obtain

$$\begin{array}{lll}{\mathscr{H}}_1\hfill & =\hfill & {\displaystyle \sum _{i=1}^2[\frac{\mathcal{G}{m}_i}{2r_i^5}[\left(2A-B-C\right){\left(r_i\cdot I\right)}^2+\left(2B-A-C\right){\left(r_i\cdot J\right)}^2}\hfill \\ \hfill & +\hfill & \left(2C-A-B\right){\left(r_i\cdot K\right)}^2],\hfill \\ \hfill & +\hfill & O\left({\left[\frac{R_{\text{eq}}}{r_1}\right]}^3\right)+O\left({\left[\frac{R_{\text{eq}}}{r_2}\right]}^3\right)+O\left({\left[\frac{r_1}{r_2}\right]}^3\right)+O\left(\frac{r_1}{r_2}{\left[\frac{R_{\text{eq}}}{r_2}\right]}^3\right).\hfill \end{array}$$(C.14)

We notice that r₀ does not appear in the Hamiltonian function, thanks to the reduction of the total barycentre. As a result, the isolated constant term ∥p₀∥²/(2β₀) in Eq. (C.8) can be dropped from the Hamiltonian function.

Appendix C.2 Secular system

The Hamiltonian function obtained above can be decomposed into

$$\mathscr{H}={\mathscr{H}}_0+ϵ{\mathscr{H}}_1,$$(C.15)

where ℋ₀ = ℋ_E + ℋ_K is the dominant integrable part and ϵℋ₁ = ℋ_M + ℋ_I is a perturbation (we introduce a factor ϵ ≪ 1 to emphasise the smallness of the perturbation). Assuming that there is no resonance involving the mean motions and/or the planet’s spin rate, the long-term dynamics of the system at first order in ≪ is given by the classic averaging method.

In order to average ≪ℋ₁ over the fast angles of the rotational motion of the planet, one would need to express ℋ_E in action-angle coordinates. Yet, because of the planet’s relaxation to hydrostatic equilibrium, we know that the planet’s free rotational motion is close to that of an axisymmetric ellipsoid, for which the Andoyer angles 𝒢 and ℓ circulate in an independent manner (i.e. the motion is composed of a main rotation about the constant rotational angular momentum G and a slower rotation about the principal axis K). Assuming that none of the two corresponding fast frequencies⁷ is involved in a resonance, we can therefore bypass the standard procedure and average the Hamiltonian function over 𝒢 and ℓ separately (see Boué & Laskar 2006; Vaillant et al. 2019). Expressing the vectors I, J, and K in terms of the Andoyer angles, we obtain

$$\begin{array}{lll}{\langle {\left(r_i\cdot I\right)}^2\rangle }_{g,\ell }\hfill & =\hfill & {\langle {\left(r_i\cdot J\right)}^2\rangle }_{g,\ell }\hfill \\ \hfill & =\hfill & \frac{1+{\cos }^{2}J}{4}{r}_i^2-\frac{1}{2}\left(1-\frac{3}{2}{\sin }^{2}J\right){\left(r_i\cdot w\right)}^2,\hfill \\ {\langle {\left(r_i\cdot K\right)}^2\rangle }_{g,\ell }\hfill & =\hfill & \frac{{\sin }^{2}J}{2}{r}_i^2+\left(1-\frac{3}{2}{\sin }^{2}J\right){\left(r_i\cdot w\right)}^2,\hfill \end{array}$$(C.16)

where w = G/𝒢 and J is the angle between G and K, such that

$${\langle {\mathscr{H}}_1\rangle }_{g,\ell }=-\frac{2C-A-B}{2}\left(1-\frac{3}{2}{\sin }^{2}J\right){\displaystyle \sum _{i=1}^2\frac{\mathcal{G}{m}_i}{2r_i^3}\left[1-3\frac{{\left(r_i\cdot w\right)}^2}{r_i^2}\right].}$$(C.17)

Upon averaging, J becomes a constant angle, and the planet’s rotational dynamics is reduced to the orientation of the unitary angular momentum vector w.

It remains to average the perturbation Hamiltonian over the fast orbital angles, which are the mean anomalies M₁ and M₂. The required integrals are performed using the classical formulas (see e.g. Appendix A of Boué & Laskar 2006). Dropping the constant parts, the secular Hamiltonian function is finally

$$\overline{\mathscr{H}}={\overline{\mathscr{H}}}_{\text{M}}+{\overline{\mathscr{H}}}_1,$$(C.18)

where

$${\overline{\mathscr{H}}}_{\text{M}}=\frac{\gamma }{6}\left[1-6e_1^2-3\left(1-e_1^2\right){\left(k_1\cdot k_2\right)}^2+15{\left({\text{e}}_1\cdot k_2\right)}^2\right],$$(C.19)

and

$${\overline{\mathscr{H}}}_1={\displaystyle \sum _{i=1}^2\frac{{\alpha }_i}{6}}\left[1-3{\left(k_i\cdot w\right)}^2\right],$$(C.20)

in which we have defined

$$\gamma =\frac{3\mathcal{G}{\beta }_1{m}_2{a}_1^2}{4a_2^3{\left(1-e_2^2\right)}^{3/2}},$$(C.21)

and

$${\alpha }_i=\frac{3\mathcal{G}{m}_0{m}_i{R}_{\text{eq}}^2{J}_1}{2a_i^3{\left(1-e_i^2\right)}^{3/2}}\left(1-\frac{3}{2}{\sin }^{2}J\right),$$(C.22)

where

$$J_2=\frac{2C-A-B}{2m_0{R}_{\text{eq}}^2}.$$(C.23)

In these expressions, e₁ and e₂ are the eccentricity vectors of the satellite around the planet and of the planet-satellite barycentre around the star⁸ (i.e. the vector pointing to the pericentre and whose norm is the eccentricity), and k₁ and k₂ are the unitary angular momentum vectors of the two orbits.

Following Correia et al. (2011), we express the equations of motion in terms of the (non-canonical) vectorial elements G, e₁, and G₁, where G₁ is the orbital angular momentum of the satellite:

$$G_1={\beta }_1\sqrt{{\mu }_1{a}_1\left(1-e_1^2\right)}{k}_1.$$(C.24)

We obtain

$$\begin{array}{lll}\dot{G}\hfill & =\hfill & {\displaystyle \sum _{i=1}^2{\alpha }_i\left(k_i\cdot w\right)\left(k_i\times w\right),}\hfill \\ {\dot{e}}_1\hfill & =\hfill & \frac{{\alpha }_1}{2\Vert G_1\Vert }\left[2\left(k_i\cdot w\right)\left(e_1\times w\right)-\left(1-5{\left(k_1\cdot w\right)}^2\right)\left(k_1\times e_1\right)\right]\hfill \\ \hfill & +\hfill & \frac{\gamma \left(1-e_1^2\right)}{\Vert G_1\Vert }[2\left(k_1\times e_1\right)+\left(k_1\cdot k_2\right)\left(e_1\times k_2\right)\hfill \\ \hfill & \hfill & -5\left(e_1\cdot k_2\right)\left(k_1\times k_2\right)],\hfill \\ {\dot{G}}_1\hfill & =\hfill & {\alpha }_1\left(k_1\cdot w\right)\left(k_1\times w\right)+\gamma \left(1-e_1^2\right)\left(k_1\cdot k_2\right)\left(k_1\times k_2\right)\hfill \\ \hfill & -\hfill & 5\gamma \left(e_1\cdot k_2\right)\left(e_1\times k_2\right).\hfill \end{array}$$(C.25)

The elements e₂ and k₂ of the planet-satellite barycentre around the star require a specific treatment. Indeed, the three-body system that we consider here is actually not isolated, but subject to perturbations from additional planets. Their direct interactions with the satellite’s orbit and the planet’s equatorial bulge are negligible; however, they are the main contributors to the secular motion of the planet-satellite barycentre around the star. This contribution is essential because it gives rise to the secular spin-orbit resonances. For this reason, we do not integrate the temporal evolution of e₂ and k₂ but take them as known functions of time obtained from earlier works. Given the small size of the satellites considered in this article, we can safely assume that its influence did not alter noticeably the motion of the barycentre of Uranus and its satellite system around the Sun. Hence, we compute e₂ and k₂ from the secular orbital solution of Laskar (1990) described in Appendix B. The conversion from the variables z and ζ to the vectors e₂ and k₂ is straightforward. This solution acts as an indirect forcing term to the dynamics of the satellite and the spin axis of the planet.

The constant value of Uranus’s rotational angular momentum 𝒢 is quite uncertain, as both its spin rate ω and moment of inertia are poorly known (see the discussion in Sect 2.2). For this reason, it is enough to use the approximate expression for a principal axis rotator $G\approx C\omega =\lambda \omega m_0{R}_{\text{eq}}^2$, where λ is the normalised polar moment of inertia. We adopt the nominal value of λω described in Sect. 2.2. We also neglect the sin² J term in Eq. (C.22) which is not known for Uranus but expected to be extremely small⁹. The migration of the satellite is modelled by applying a slow drift in a1 in the equations of motion according to the migration law in Eq. (9). As discussed in Sect. 2.4, we neglect the simultaneous decrease in Uranus’s spin rate ω, whose contribution is smaller than our uncertainty on the value of 𝒢.

Appendix D High-order secular spin-orbit resonances

For a given planet, the strongest secular spin-orbit resonances are of first order in the planet’s orbital inclination (see e.g. Saillenfest et al. 2019b). First-order resonances have resonant angles of the form σ_k = ψ + ϕ_k, where ψ is the spin-axis precession angle and ϕ_k is the circulating angle of the kth harmonic in the planet’s nodal precession spectrum (see Eq. 6). The secular spin-axis dynamics, however, is not restricted to first-order resonances. At second order in inclination, resonances involve two different harmonics with resonant angles of the form σ_j,k = 2ψ +ϕ_j+ ϕ_k. At order three, there exist five different kinds of resonances whose properties are described by Saillenfest et al. (2019b).

In the case of Uranus on its current orbit, resonances are rare and isolated from each other. Yet, harmonic k = 2 has quite a large amplitude (see Table 1), and its frequency is not far from those of harmonics k = 5 and 8. This suggests that high-order resonances could play a substantial role in the dynamics in the neighbourhood of resonance σ₂. And indeed, we found that some fraction of the trajectories shown in Fig. 10 for satellite masses in the neighbourhood of m = m₂ is produced by stable captures in second-order resonances. Figure D.1 shows an example of such trajectories: σ₂ = ψ + ϕ₂ first starts to oscillate, but it is soon replaced by σ₂₅ = 2ψ + ϕ₂ + ϕ₅, while both σ₂ and σ₅ circulate. The nearby resonance σ₈ does not play any role in this example.

Fig. D.1 Example of stable capture in a secular spin-orbit resonance of second order. Black and blue are used for quantities related to the planet’s spin axis and satellite’s orbit, respectively. The numerical integration is performed with the fully coupled secular model presented in Sect. 3. The mass of the satellite is m/M = 2 × 10−3. The satellite is initialised close to the local circular Laplace equilibrium at a = 42 Req and it is made migrating outwards over Δa = 9 Req.

Appendix E Secular resonances of a massless satellite

In Sect. 3.4, we investigate the origin of the peculiar dynamical pathways through which the satellite can collide into the planet. To this end, we simplify the problem and consider the dynamics of a massless satellite around a planet having a frozen orbit and spin-axis orientation. The secular Hamiltonian function of the satellite can be written ℋ = k_Pℋ_P + k_⊙ℋ_⊙. Expanding both contributions to quadrupole order, Saillenfest & Lari (2021) define the two constant coefficients as

$$\begin{array}{lll}{k}_{\text{P}}=\frac{3}{4}\frac{\mathcal{G}M}{a}{J}_2\frac{R_{\text{eq}}^2}{a^2}\hfill & \text{and}\hfill & k_{\odot }=\frac{3}{8}\frac{\mathcal{G}{m}_{\odot }}{a_{\odot }{\left(1-e_{\odot }^2\right)}^{3/2}}\frac{a^2}{a_{\odot }^2}\hfill \end{array},$$(E.1)

and the two Hamiltonian functions as

$${\mathscr{H}}_{\text{P}}=\frac{1-3{\cos }^2{I}_Q}{3{\left(1-e^2\right)}^{3/2}},$$(E.2)

and

$$\begin{array}{lll}{\mathscr{H}}_{\odot }\hfill & =\hfill & -\frac{1}{8}[8e^2+2\left(3e^2+2\right)\left(2C^2{\cos }^2{I}_Q+S^2{\sin }^2{I}_Q-2\right)\hfill \\ \hfill & \hfill & +8CS\left(3e^2+2\right)\cos I_Q\sin I_Q\cos {\text{Ω}}_Q\hfill \\ \hfill & \hfill & +5S^2{e}^2{\left(\cos I_Q+1\right)}^2\cos \left(2{\omega }_Q+2{\text{Ω}}_Q\right)\hfill \\ \hfill & \hfill & -20CS{e}^2\left(\cos I_Q+1\right)\sin I_Q\cos \left(2{\omega }_Q+{\text{Ω}}_Q\right)\hfill \\ \hfill & \hfill & +10\left(3C^2-1\right)e^2{\sin }^2{I}_Q\cos \left(2{\omega }_Q\right)\hfill \\ \hfill & \hfill & -20CS{e}^2\left(\cos I_Q-1\right)\sin I_Q\cos \left(2{\omega }_Q+{\text{Ω}}_Q\right)\hfill \\ \hfill & \hfill & +5S^2{e}^2{\left(\cos I_Q-1\right)}^2\cos \left(2{\omega }_Q-2{\text{Ω}}_Q\right)\hfill \\ \hfill & \hfill & +2S^2\left(3e^2+2\right){\sin }^2{I}_Q\cos \left(2{\text{Ω}}_Q\right)].\hfill \end{array}$$(E.3)

In this expression, C = cos ε and S = sin ε are constant, and the orbital angles of the satellite (I_Q, ω_Q, Ω_Q) are measured with respect to the equator and equinox of the planet. All other quantities that appear in the Hamiltonian function are described in the main text.

Because of the factor (1 − e²)^−3/2 in Eq. (E.2), k_Pℋ_P is the dominant contribution to the dynamics whenever the eccentricity grows large. We consider the intermediate regime in which k_Pℋ_P is the dominant integrable part of the Hamiltonian and k_⊙ℋ_⊙ is a small perturbation. Written in terms of the Delaunay elements, k_⊙ℋ_⊙ is expressed in action-angle coordinates, so we can directly apply a perturbative method. At first order in the perturbation, the possible resonances are those directly appearing in Eq. (E.3); for each of them, the resonance angle ξ is a linear combination of ω_Q and Ω_Q. These resonances are listed in Table 3 of the main text.

Using a perturbative approach, we can study the orbital dynamics of the satellite in the vicinity of each resonance. The method is the following (see e.g. Saillenfest et al. 2019a, Talu et al. 2021): The angle ξ is first taken as coordinate by a canonical change of variables (ω_Q, Ω_Q) → (ξ,γ). The conjugate momenta of ξ and γ are written Ξ and Γ. In the vicinity of the resonance, ξ is a slow angle and γ is a fast angle. At first order in the perturbation, the long-term dynamics is described by averaging the Hamiltonian function over γ. The Hamiltonian system obtained is independent of γ, which means that Γ is a constant of motion and can be used as parameter. We end up with a one-degree-of-freedom system and any possible trajectory for the satellite can be represented as a level curve of the Hamiltonian function. In practice, Γ is rewritten as a constant quantity K, whose expression for each resonance is given in Table 3. By virtue of the constant nature of K, it is equivalent to represent the level curves of the Hamiltonian function in the (e, ξ) plane or in the (I_Q, ξ) plane.

When it emerges from chaos and reaches the weakly perturbed regime, the satellite can have any value of K. For definitiveness, we consider here for each resonance the value of K such that e = 0 for I_Q = I₀, where I₀ is the nominal inclination of the resonance (see Table 3). Even though this represents a loss of generality, we later see that studying just one value of K is enough to understand the correct hierarchy between the various resonances. The case ξ = Ω_Q is special, because the constant quantity is the eccentricity e itself. By putting e = 0, we obtain the Hamiltonian function that describes the circular Laplace equilibria (see Tremaine et al. 2009). The ‘resonance’ in this case is the libration island surrounding the orthogonal equilibrium (called Laplace state P₂ by Saillenfest & Lari 2021), for which the minimum and maximum value of I_Q along the separatrix are given by

Fig. E.1 Geometry of the phase portrait describing the orbital dynamics of a satellite in the vicinity of one of the secular resonances listed in Table 3 (except the last one). The eccentricity e and equatorial inclination IQ of the satellite are linked through a constant quantity K. The origin e = 0 of the graph corresponds to the nominal inclination I0 of the resonance.

$${\cos }^2{I}_{\text{Q}}=\frac{1+u-\sqrt{1+u^2+2u\cos \left(2\epsilon \right)}}{2u},$$(E.4)

where $u=r_{\text{M}}^4/a^5$. All other resonances in Table 3 produce very similar phase portraits. Since some resonances are very large, we do not use the pendulum approximation but solve for the exact resonant Hamiltonian function. The structure of the phase space is illustrated in Fig. E.1. The width Δe of the separatrix gives the maximum eccentricity reachable by the satellite inside the resonance. It is linked to an analogous width ΔI_Q through the constant quantity K. The widths of each resonance as a function of the parameters are shown in Figs. 16 and 17 of the main text.

Appendix F Effect of a satellite collision on the obliquity

In Sect. 2.4, we observe that large regions of the parameter space lead to a resonance capture of Uranus. As the satellite continues to migrate outwards, the system converges towards a configuration for which the satellite is highly unstable. The effects of this destabilisation are investigated in Sect. 3; they can lead to a direct collision of the satellite into the planet. In this section, we investigate what would be the effect of a satellite-planet collision on the final spin rate and obliquity of the planet.

As explained in Sect. 4.2, the total angular momentum L of the planet and its satellite can be considered constant during the collision. The constant vector L can be written

$$L=L_{\text{P}}+L_{\text{S}},$$(F.1)

where L_P is the contribution of the planet’s spin and L_S is the contribution of the satellite’s orbit. After the removal of the satellite, the obliquity of Uranus is not expected to vary at all due to external perturbations because Uranus is released far from any kind of spin-orbit resonance. The magnitude of L can therefore be computed using the current parameters of Uranus. A set of inner satellites can also be included; however, the current satellites of Uranus only contribute to 1% of L = ∥L∥, which is much smaller than our uncertainty on Uranus’s moment of inertia. To leading order, ∥L_P∥ ∝ ω and ∥L_P∥ ≫ ∥L_S∥, so the relative change in the planet’s spin rate is

$$\frac{\text{Δ}\omega }{\omega }\approx \frac{L-L_{\text{P}}}{L_{\text{P}}}\approx \frac{L_{\text{S}}}{L_{\text{P}}}\cos I_{\text{Q}},$$(F.2)

where L_P = ∥L_P∥, L_S = ∥L_S∥, and I_Q is the orbital inclination of the satellite with respect to the planet’s equator at the collision time (i.e. the angle between L_P and L_S). Using the gyroscopic approximation for the spin of the planet, the ratio L_S/L_P can be written as

$$\frac{L_{\text{S}}}{L_{\text{P}}}=\frac{m}{M}\frac{n}{\omega }\frac{a^2}{R_{\text{eq}}^2}\frac{\sqrt{1-e^2}}{\lambda },$$(F.3)

where n and e are the satellite’s mean motion and eccentricity, and other quantities are described in Sect. 2.1. The collision is expected to occur when the semi-major axis of the satellite is a ≈ r_M as defined in Eq. (2). Moreover, the collision condition requires that a(1 − e) ≲ R_eq which gives a minimum bound for the eccentricity value at the collision time. When using the physical parameters of Uranus in Eq. (F.3), we obtain a value of about 35 m/M, where the leading factor varies from 32 to 39 when assuming a 10% uncertainty on the value of ωλ as in Sect. 2.2. The order of magnitude for the relative change in the planet’s spin rate is therefore

$$\frac{\text{Δ}\omega }{\omega }\approx 35\frac{m}{M}\cos I_{\text{Q}}.$$(F.4)

The variation in the planet’s spin-axis orientation due to the impact is quantified by the angle φ between L and L_P. From Eq. (F.1), it can be expressed as

$$\sin \phi =\frac{L_{\text{S}}}{L}\sin I_Q,$$(F.5)

that is,

$$\sin \phi =35\frac{m}{M}\sin I_Q.$$(F.6)

Under favourable longitude of node of the satellite at the time of collision, the spin-axis deviation φ directly equals the obliquity variation Δε. As shown in Sect. 3.4, satellite collisions only occur when I_Q ≈ 55° or I_Q ≈ 125°. Therefore, for the mass range discussed throughout this article, the change in spin rate due to the satellite’s impact is a few percent at most, and the change in obliquity is no larger than a few degrees.

The maximum spin-axis deviation φ would be obtained for a grazing impact with I_Q = 90°. Hence, an upper bound for the obliquity change in radians due to the satellite collision is Δε ≈ 35 m/M.

References