Open Access
Issue
A&A
Volume 668, December 2022
Article Number A108
Number of page(s) 35
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202243953
Published online 12 December 2022

© M. Saillenfest et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 Introduction

The 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 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. s7 or s8, 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 S1, where it degenerates into a continuum of equilibria. When it comes to the spin-axis dynamics of the planet, S1 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, S1 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

(1)

where m is the mass of the satellite, M is the mass of the planet, J2 is its second zonal gravity coefficient, and Req is a normalising radius (that we take equal to the planet’s equatorial radius). The characteristic length rM 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 radius1 rL introduced by Tremaine et al. (2009), as

(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 rM for Uranus is about 53 Req. We neglect for now the influence of inner satellites in the effective value of J2 (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

(3)

where ε is the planet’s obliquity, a is the semi-major axis of the satellite, and IL 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 , and it can be written as

(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 IL is given by the formula

(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 planet2 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:

(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 Ek and Sk are real constants and the angles θk and ϕk evolve linearly over time t with frequencies µk and vk:

(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−9. 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 sj. 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 Ω0 in Eq. (3) and the forcing frequencies µk and vk appearing in Eq. (7). As detailed by Saillenfest et al. (2019b), the largest resonances are those of order 1 in the amplitudes {Sk}, 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 ζ series3. 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 Ω0, such that the relation Ω0 + vj 0 is maintained.

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

(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 vj), 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: vj should simply be replaced by a combination of frequencies vk 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.

Table 1

First twenty terms of Uranus’s inclination and longitude of ascending node in the J2000 ecliptic and equinox reference frame.

thumbnail 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 = /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 . 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 ∣vk∣ 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 Req = 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−1, 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 ν15 = 𝒢7 − 𝒢8 + s7 (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 mmin for the satellite lying between about 3.2 × 10−4 and 3.8 × 10−4 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−4, that is, m ≈ 3 × 1022 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 ν15 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 mmin 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 is much smaller than the age of the Solar System. In this regard, the v9 and v15 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 mmin 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 mk, where k is the index of the resonant term in Uranus’s ζ series. When mk is large, we see that our simplified analytical expression systematically underestimates its value (i.e. mminmk), 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 E1 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 E1, 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 S1. 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 = rM (blue zone in Fig. 1), then touch the singular point S1 (red curve in Fig. 1), at which point we define the mass mk. 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 v9 and v15 resonances are indeed very thin, and that Uranus also possesses two extremely weak retrograde resonances (v19 and v32). Since m2 is only slightly smaller than m8 (see Table 2), the v8 resonance appears in panel d but does not connect to S1 (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 rM, 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−1. Compared to the measured expansion rates of the Moon from the Earth (4 cm yr−1 ; see e.g. Williams & Boggs 2016), Ganymede from Jupiter (11 cm yr−1; see Lainey et al. 2009) and Titan from Saturn (11 cm yr−1; 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 rM in order for the tilting mechanism to fully operate (see Fig. 2). Since the characteristic radius rM is quite large for Uranus (rM ≈ 53 Req), 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 rM to a = rM, the equivalent constant parameter k2/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 ≳ 105). 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.

Table 2

Minimum mass of the satellite in order for Uranus to reach the few closest resonances and be tilted from 0° to 90°.

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

thumbnail 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 E1. 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 S1. 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

(9)

where a0 is the initial location of the satellite, and t0 and τ are parameters that have the dimension of time. In the following, we adopt t0 = 4.5 Gyr and vary the parameter τ, that we call the “migration timescale” (as indeed . 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 mk 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 rM to 1 rM 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 J2 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 v3 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 v15 resonance is crossed as well but it hardly affects the trajectory at all. Panel b shows a case where the satellite has mass m5 but is captured in the v7 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 v5 is crossed without capture but the planet is then captured and tilted in the v2 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 ≈ rM. 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 mk required for a given resonance to exist and go to 90° (see Table 2). The satellite is initialised at a0 = 0.8 rM, so as to focus on the effects of resonance vk and minimise the migration range required for the tilting (see Fig. 2). Over a grid of initial obliquity ε0 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 (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 ε0. Panel b for m = m3, 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 v3 resonance (as in Fig. 5c). For ε0 ≳ 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 ε0 ≲ 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 v7 and v5. In panels a and c, on the contrary, the resonances v15 and v9 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 v2 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 v3 (i.e. a resonance with the nodal orbital precession mode of Neptune, s8): it requires a minimum mass for the satellite of only m/M ≈ 4 × 10−4 and has a high tilting efficiency for a large range of initial obliquities (0oε0 ≲ 20°). For larger masses, m/M ≳ 1 to 2 × 10−3 the nodal precession spectrum of Uranus allows for several other possibilities: the resonances with v7= g5 − 𝒢7 + s1 and v5 = −𝒢5 + g6 + s6 are quite good candidates, but in a more limited range of initial obliquities (ε0 ≲ 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 v2 = s7 (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 (ε0 ≲ 45°) and for any conceivable migration rate for the satellite. However, it requires quite massive satellite with m/M ≳ 2.3 × 10−3 : 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.

thumbnail 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.

thumbnail 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.

thumbnail 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 G1, and the eccentricity vector of the satellite e1. The orbital angular momentum of the planet-satellite barycentre around the star and its eccentricity vector (written G2 and e2 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.8rM 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 years4.

thumbnail 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 (σ2 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 S1 (see e.g. Fig. 2), all available resonances are actually “nearby” and contribute to the obliquity kicks. The largest of them is the v2 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.

thumbnail 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.

thumbnail 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 ε0 = 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 ε0 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 a0 = 42 Req (that is, about 0.8 rM). The satellite is started close to the circular Laplace equilibrium, with an eccentricity e = 10−4 and an inclination ILP = 10−4 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 a0 ≈ 0.8 rM (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 rM, 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 v15 and v9 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 v7 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 v8 appears quite far from the corresponding obliquity peak in Fig. 10, because resonance v8 is very distorted by the large neighbour resonance v2. In fact, resonance v2 is so large that it almost completely overlaps resonance v8; 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ψ + ϕ2 + ϕ5 is given in Appendix D; this example confirms that the dynamics in the neighbourhood of resonance v2 are very rich and not limited to interactions with harmonic ϕ2 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 v7, v5, v2, v8, or v6. The minimum mass of the satellite for which successful trajectories are obtained is m/M ≈ 1.2 × 10−3 (i.e. mm7). For smaller masses, no trajectory in our sample goes beyond the current obliquity of Uranus. As expected from the previous sections, captures into resonance v3 for m/M ≈ 4.4 × 10−4 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 Req (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 Req). 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−3 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 Req, 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 Req. 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).

thumbnail 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 ν2 (partially overlapping with resonance ν8). This is not a coincidence: resonance ν2 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 a0 of the satellite. The choice a0 ≈ 0.8 rM 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 ν2 for m/M > 2 × 10−3 produces the largest fraction of successful trajectories, but resonances ν5 and ν7 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 v3 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 a0 ≲ 35 Req and migrates fast (see the violet shade for m/M ≈ 7 × 10−4). The success ratios, however, remain very small in resonance v3.

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 v2, 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 v5 and v7, with success ratios up to 50%. Thanks to the large width of resonance v2, 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 a0 ≳ 50 Req. 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 Req. This reduced migration range increases the effective quality factor Q that would be required for Uranus (see Fig. 3).

thumbnail 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.

thumbnail 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.

thumbnail 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 IQ the orbital inclination of the satellite measured with respect to the planet’s equator and IC 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 IQ are localised in two very narrow bands that can be visualised in Fig. 14. These two bands are localised at IQ ≈ 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 IC 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, IQ) 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 IQ: we clearly see that no collision can possibly occur outside of IQ ≈ 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

(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 {kp, k} and of the Hamiltonian functions {ℋP,ℋ}expanded at quadrupolar order are given in Appendix E. Because ℋP contains a factor (1 − e2)−3/2, the first term in Eq. (10) becomes strongly dominant whenever the eccentricity grows large. Indeed, the contribution of kPP produces a fast precession of the satellite’s argument of pericentre ωQ and longitude of ascending node ΩQ, at rates

(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 kPP 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 kPP 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) ∈ ℤ2 \ {(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 IQ, that we note I0.

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 ΔIQ 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 rM, resonances are small and well separated from each other for any obliquity value. When the satellite is close to a = rM, 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 simulations5. 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 IQ = 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 IQ 56° (or IQ 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 IQ 46° (or equivalently ωQ − ΩQ at IQ 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 IQ 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 IQ 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.

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

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

Table 3

Properties of the secular resonances appearing in the orbital dynamics of a massless satellite at first order in the solar perturbation.

thumbnail 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°).

thumbnail 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−3 (i.e. about the mass of Jupiter’s moon Ganymede, m ≈ 15 × 1022 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 a0 ranging from about 40 to 50 Req. 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 J2 (Tremaine et al. 2009). If we take the current satellites of Uranus into account at their current locations, the effective value of J2 felt by the hypothetical distant satellite would be multiplied by a factor k ≈ 5.4. As a result, the Laplace radius is multiplied by k1/5 ≈ 1.4 (see Eq. (2)) and rM becomes 75Req instead of 53Req. 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−4. 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−4 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−3 around Uranus, as proposed here. However, even though the typical ratio of 10−4 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−3 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−2), 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.

thumbnail 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.

thumbnail 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−3 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−5, and the mass ratio for Titan is about 2 × 10−4. 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−4 (see Sect. 2). With this minimum mass, a satellite allows for a capture of Uranus in secular spin-orbit resonance with v3 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−4 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 J2 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 J2 is increased by a factor k by a set of inner satellites, then rM is enlarged by a factor k1/5 and the mass needed to reach a given resonance is divided by k2/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 rM to 75 Req and divide the mass needed for the ancient satellite by roughly a factor two. In particular, the strong resonance v2 showing the highest success probabilities in our scenario would be reachable by an ancient satellite with mass m/M ≈ 10−3 instead of 2 × 10−3. 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 v3 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 a0 and larger mass (i.e. in the lower end of the coloured bands). For resonance v3, however, a non-negligible fraction of successful trajectories can be obtained only for a0 ≲ 33 Req (see panels c and f of Fig. 12). Such small values of a0 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 v3 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 a0 > rM. In order for the satellite to migrate inwards, it should be retrograde with respect to the planet’s spin motion. Retrograde massive satellites with IQ ≈ 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 v3 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,

(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−3, 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−1 ). The values obtained are much less than Uranus’s escape speed (21.4 km s−1), 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−3. 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−3, 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−4 at most.

thumbnail 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−4. 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−4 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−4 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 mdisc, the characteristic lifetime of the disc Tdisc = mdisc/disc is

where TRoche 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 mdisc/M ≈ 4.4 × 10−4. The formation of the new generation of satellites is therefore extremely fast at the beginning, and it slows down asymptotically as mdisct−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 Req. This distance is to be compared to the synchronous orbit around Uranus, which lies somewhat below 3.5 Req (the migration from the fluid Roche limit at 2.5 Req 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 away6. 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 k2/Q; see Sect. 2.3) is smaller than that needed to move the distant satellite over 10 Req, 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 ν15 = 𝒢7 − 𝒢8 + s7 (where 𝒢j and sj 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 ν15, Uranus must have had a satellite with minimum mass m/M = 3.5 × 10−4 (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 v3 = s8. Triggering a resonance with v3 and tilting Uranus up to a large obliquity would require a satellite with minimum mass m/M = 4.4 × 10−4, that is, m ≈ 4 × 1022 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−1 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−1 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−4 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 ν2 = s7 which can be reached for a satellite mass m/M = 2.2 × 10−3. Other promising resonance candidates are ν7 = 𝒢5 − 𝒢7 + s7 and ν5 = −𝒢5 + 𝒢6 + s6, which can be reached for masses m/M = 1.3 × 10−3 and 1.7 × 10−3, 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 IQ ≈ 56°, and to its symmetric counterpart 2ω − Ω at IQ ≈ 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−4, 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 ν3, 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−3, the probability of success through a capture in resonance ν2 can be as large as 86%. For slightly smaller masses (m/M ≈ 1.5 to 2 × 10−3), resonances ν5 and ν7 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 ν2 is special because it is a secular spin-orbit resonance with Uranus’s own nodal precession mode s7. 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 ν2 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 (a0 ≈ 50 Req), the tilting mechanism via resonance ν2 can operate down to very small migration ranges (only a few Req, instead of 10 Req 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 Req 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−3 at a distance of a0 ≈ 50 Req 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 Req) up to their current distance (Oberon is at a ≈ 23 Req). 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 (IQ ≈ 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−4; 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 Req, the outer one would feel a change in the effective J2 parameter of Uranus, which would increase its characteristic radius rM. 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 S0 = sin(I0/2) and angle ϕ0 evolving at frequency v0 (see Eq. 6), the Hamiltonian function governing the secular spin-axis dynamics of the planet is simplified into

(A.1)

where Φ0 is the momentum conjugate to ϕ0, 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 −σ = −ψϕ0 and dropping unnecessary constant parts, the Hamiltonian function can be rewritten

(A.2)

using the rescaled timescale dτ = pdt, where

(A.3)

and 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 v0 < 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

(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

(A.5)

The two curves defined by Eqs. (A.4) and (A.5) connect at γ = γcrit, where . 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 I0. Therefore, the two previous conditions can be formulated as a minimum value for I0 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 that is, Imin ≈ 9.74°.

thumbnail 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 and its value is given by tan .

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. v0 = s7) or with a harmonic stemming from an other planet (e.g. Neptune, as investigated by Rogoszinski & Hamilton 2021, that is, v0 = s8): 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 and (Tables VI and VII of Laskar 1990) by the matrix 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.

thumbnail 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 and , 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. s1, s2, 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.

Table B.1

Quasi-periodic decomposition of Uranus’s eccentricity and longitude of perihelion (variable z).

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.

Table B.2

Quasi-periodic decomposition of Uranus’s inclination and longitude of ascending node (variable ζ).

Appendix C.1 Unaveraged Hamiltonian function

Our setting is similar to that of Correia et al. (2011). We consider three bodies with masses mi and positions xi measured in an inertial reference frame. The momenta conjugate to the positions xi are Xi = mii. 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 x0 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:

(C.1)

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

(C.2)

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

(C.3)

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

(C.4)

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

(C.5)

with conjugate momenta

(C.6)

where the reduced masses are defined as

(C.7)

The vector r0 describes the position of the barycentre of the whole system in the inertial reference frame, the vector r1 is the position of the satellite with respect to the planet, and the vector r2 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

(C.8)

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

(C.9)

in which we define µ1 = 𝒢(m0 + m1) and µ2 = 𝒢(m0 + m1 + m2), and ai are the associated semi-major axes. The Hamiltonian function ℋM contains the mutual Newtonian perturbations:

(C.10)

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

(C.11)

where ∥r1∥ ≪ ∥|r2∥, and by performing a multipolar expansion of ℋM using the Legendre polynomials. This gives

(C.12)

in which we write ri ≡ ∥ri.

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 yixix0, we again take advantage of the hierarchical nature of the system by noting that ∥y∥ ≪ ∥yi∥ ∀i = 1,2 and that the planet is almost spherical. By using a multipolar development of 1/yi - y∥ and integrating over the planet’s volume, we get

(C.13)

where Req is the equatorial (i.e. largest) radius of the planet. Replacing y1 and y2 by their expressions as a function of r1 and r2, we finally obtain

(C.14)

We notice that r0 does not appear in the Hamiltonian function, thanks to the reduction of the total barycentre. As a result, the isolated constant term ∥p02/(2β0) 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

(C.15)

where ℋ0 = ℋE + ℋK is the dominant integrable part and ϵ1 = ℋ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 ≪ℋ1 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 frequencies7 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

(C.16)

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

(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 M1 and M2. 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

(C.18)

where

(C.19)

and

(C.20)

in which we have defined

(C.21)

and

(C.22)

where

(C.23)

In these expressions, e1 and e2 are the eccentricity vectors of the satellite around the planet and of the planet-satellite barycentre around the star8 (i.e. the vector pointing to the pericentre and whose norm is the eccentricity), and k1 and k2 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, e1, and G1, where G1 is the orbital angular momentum of the satellite:

(C.24)

We obtain

(C.25)

The elements e2 and k2 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 e2 and k2 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 e2 and k2 from the secular orbital solution of Laskar (1990) described in Appendix B. The conversion from the variables z and ζ to the vectors e2 and k2 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 , where λ is the normalised polar moment of inertia. We adopt the nominal value of λω described in Sect. 2.2. We also neglect the sin2 J term in Eq. (C.22) which is not known for Uranus but expected to be extremely small9. 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 σ2. And indeed, we found that some fraction of the trajectories shown in Fig. 10 for satellite masses in the neighbourhood of m = m2 is produced by stable captures in second-order resonances. Figure D.1 shows an example of such trajectories: σ2 = ψ + ϕ2 first starts to oscillate, but it is soon replaced by σ25 = 2ψ + ϕ2 + ϕ5, while both σ2 and σ5 circulate. The nearby resonance σ8 does not play any role in this example.

thumbnail 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 ℋ = kPP + k. Expanding both contributions to quadrupole order, Saillenfest & Lari (2021) define the two constant coefficients as

(E.1)

and the two Hamiltonian functions as

(E.2)

and

(E.3)

In this expression, C = cos ε and S = sin ε are constant, and the orbital angles of the satellite (IQ, ω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 − e2)−3/2 in Eq. (E.2), kPP is the dominant contribution to the dynamics whenever the eccentricity grows large. We consider the intermediate regime in which kPP 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 (IQ, ξ) 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 IQ = I0, where I0 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 P2 by Saillenfest & Lari 2021), for which the minimum and maximum value of IQ along the separatrix are given by

thumbnail 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.

(E.4)

where . 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 ΔIQ 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

(F.1)

where LP is the contribution of the planet’s spin and LS 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, ∥LP∥ ∝ ω and ∥LP∥ ≫ ∥LS∥, so the relative change in the planet’s spin rate is

(F.2)

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

(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 arM as defined in Eq. (2). Moreover, the collision condition requires that a(1 − e) ≲ Req 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

(F.4)

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

(F.5)

that is,

(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 IQ ≈ 55° or IQ ≈ 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 IQ = 90°. Hence, an upper bound for the obliquity change in radians due to the satellite collision is Δε ≈ 35 m/M.

References

  1. Agnor, C. B., & Hamilton, D. P. 2006, Nature, 441, 192 [Google Scholar]
  2. Archinal, B. A., Acton, C. H., A’Hearn, M.F., et al. 2018, Celest. Mech. Dyn. Astron., 130, 22 [NASA ADS] [CrossRef] [Google Scholar]
  3. Atobe, K., Ida, S., & Ito, T. 2004, Icarus, 168, 223 [NASA ADS] [CrossRef] [Google Scholar]
  4. Auclair-Desrotour, P., Le Poncin-Lafitte, C., & Mathis, S. 2014, A&A, 561, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Aurnou, J., Heimpel, M., & Wicht, J. 2007, Icarus, 190, 110 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bailey, E., & Stevenson, D. J. 2021, Planet. Sci. J., 2, 64 [Google Scholar]
  7. Batygin, K. 2018, AJ, 155, 178 [Google Scholar]
  8. Batygin, K., & Morbidelli, A. 2020, ApJ, 894, 143 [Google Scholar]
  9. Batygin, K., Adams, F. C., Brown, M. E., & Becker, J. C. 2019, Phys. Rep., 805, 1 [Google Scholar]
  10. Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 [Google Scholar]
  11. Boué, G., & Laskar, J. 2006, Icarus, 185, 312 [Google Scholar]
  12. Boué, G., & Laskar, J. 2010, ApJ, 712, L44 [CrossRef] [Google Scholar]
  13. Brouwers, M. G., Bonsor, A., & Malamud, U. 2022, MNRAS, 509, 2404 [NASA ADS] [Google Scholar]
  14. Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, Nat. Astron., 2, 138 [Google Scholar]
  15. Canup, R. M. 2010, Nature, 468, 943 [Google Scholar]
  16. Canup, R. M., & Esposito, L. W. 1995, Icarus, 113, 331 [NASA ADS] [CrossRef] [Google Scholar]
  17. Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834 [Google Scholar]
  18. Chandler, S. C. 1891, AJ, 11, 59 [NASA ADS] [CrossRef] [Google Scholar]
  19. Chau, A., Reinhardt, C., Izidoro, A., Stadel, J., & Helled, R. 2021, MNRAS, 502, 1647 [Google Scholar]
  20. Chirikov, B. V. 1979, Phys. Rep., 52, 263 [Google Scholar]
  21. Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J. 2018, Icarus, 311, 340 [Google Scholar]
  22. Clement, M. S., Kaib, N. A., Raymond, S. N., & Chambers, J. E. 2021, Icarus, 367, 114585 [NASA ADS] [CrossRef] [Google Scholar]
  23. Colombo, G. 1966, AJ, 71, 891 [Google Scholar]
  24. Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 111, 105 [NASA ADS] [CrossRef] [Google Scholar]
  25. Crida, A. 2020, Nat. Astron., 4, 1024 [Google Scholar]
  26. Crida, A., & Charnoz, S. 2012, Science, 338, 1196 [Google Scholar]
  27. Cuk, M., El Moutamid, M., & Tiscareno, M. S. 2020, Planet. Sci. J., 1, 22 [Google Scholar]
  28. D’Angelo, G., Weidenschilling, S.J., Lissauer, J.J., & Bodenheimer, P. 2021, Icarus, 355, 114087 [Google Scholar]
  29. Daquin, J., Legnaro, E., Gkolias, I., & Efthymiopoulos, C. 2022, Celest. Mech. Dyn. Astron., 134, 6 [NASA ADS] [CrossRef] [Google Scholar]
  30. Dittmann, A. J. 2021, MNRAS, 508, 1842 [Google Scholar]
  31. Dong, J., Jiang, Y.-F., & Armitage, P. J. 2021, ApJ, 921, 54 [NASA ADS] [CrossRef] [Google Scholar]
  32. Duncan, M. J., & Lissauer, J. J. 1997, Icarus, 125, 1 [NASA ADS] [CrossRef] [Google Scholar]
  33. Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572 [Google Scholar]
  34. Efroimsky, M., & Lainey, V. 2007, J. Geophys. Res. (Planets), 112 [Google Scholar]
  35. Estrada, P., & Durisen, R.H. 2021, in European Planetary Science Congress, EPSC2021-869 [Google Scholar]
  36. French, R. S., & Showalter, M. R. 2012, Icarus, 220, 911 [NASA ADS] [CrossRef] [Google Scholar]
  37. French, R. G., Nicholson, P. D., Cooke, M. L., et al. 1993, Icarus, 103, 163 [Google Scholar]
  38. French, R. G., Dawson, R. I., & Showalter, M. R. 2015, AJ, 149, 142 [NASA ADS] [CrossRef] [Google Scholar]
  39. Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867 [Google Scholar]
  40. Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152 [Google Scholar]
  41. Gavrilov, S. V., & Zharkov, V. N. 1977, Icarus, 32, 443 [NASA ADS] [CrossRef] [Google Scholar]
  42. Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [Google Scholar]
  43. Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466 [Google Scholar]
  44. Hamilton, D. P., & Burns, J. A. 1994, Science, 264, 550 [Google Scholar]
  45. Hamilton, D. P., & Ward, W. R. 2004, AJ, 128, 2510 [Google Scholar]
  46. Hansen, B. M. S. 2019, Sci. Adv., 5, eaaw8665 [Google Scholar]
  47. Haponiak, J., Breiter, S., & Vokrouhlický, D. 2020, Celest. Mech. Dyn. Astron., 132, 24 [Google Scholar]
  48. Harris, A. W., & Ward, W. R. 1982, Annu. Rev. Earth Planet. Sci., 10, 61 [CrossRef] [Google Scholar]
  49. Henrard, J., & Murigande, C. 1987, Celest. Mech., 40, 345 [Google Scholar]
  50. Helled, R., Anderson, J. D., & Schubert, G. 2010, Icarus, 210, 446 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hesselbrock, A. J., & Minton, D. A. 2017, Nat. Geosci., 10, 266 [Google Scholar]
  52. Hesselbrock, A. J., & Minton, D. A. 2019, AJ, 157, 30 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hoang, N. H., Mogavero, F., & Laskar, J. 2021, A&A, 654, A156 [Google Scholar]
  54. Hubbard, W. B., & Marley, M. S. 1989, Icarus, 78, 102 [Google Scholar]
  55. Hyodo, R., & Charnoz, S. 2017, AJ, 154, 34 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hyodo, R., Charnoz, S., Ohtsuki, K., & Genda, H. 2017, Icarus, 282, 195 [Google Scholar]
  57. Ida, S., Bryden, G., Lin, D. N. C., & Tanaka, H. 2000, ApJ, 534, 428 [NASA ADS] [CrossRef] [Google Scholar]
  58. Ida, S., Ueta, S., Sasaki, T., & Ishizawa, Y. 2020, Nat. Astron., 4, 880 [Google Scholar]
  59. Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220 [Google Scholar]
  60. Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364, aat2965 [Google Scholar]
  61. Izidoro, A., Morbidelli, A., Raymond, S. N., Hersant, F., & Pierens, A. 2015, A&A, 582, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Kegerreis, J. A., Teodoro, L. F. A., Eke, V. R., et al. 2018, ApJ, 861, 52 [NASA ADS] [CrossRef] [Google Scholar]
  63. Kegerreis, J. A., Eke, V. R., Gonnet, P., et al. 2019, MNRAS, 487, 5029 [Google Scholar]
  64. Konopliv, A. S., Park, R. S., Rivoldini, A., et al. 2020, Geophys. Res. Lett., 47, e90568 [NASA ADS] [CrossRef] [Google Scholar]
  65. Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
  66. Lainey, V. 2016, Celest. Mech. Dyn. Astron., 126, 145 [Google Scholar]
  67. Lainey, V., Arlot, J.-E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957 [Google Scholar]
  68. Lainey, V., Gomez Casajus, L., Fuller, J., et al. 2020, Nat. Astron., 4, 1053 [Google Scholar]
  69. Lari, G., Saillenfest, M., & Fenucci, M. 2020, A&A, 639, A40 [Google Scholar]
  70. Laskar, J. 1989, Nature, 338, 237 [Google Scholar]
  71. Laskar, J. 1990, Icarus, 88, 266 [Google Scholar]
  72. Laskar, J., & Robutel, P. 1993, Nature, 361, 608 [Google Scholar]
  73. Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322 [Google Scholar]
  74. Li, D., & Christou, A. A. 2020, AJ, 159, 184 [NASA ADS] [CrossRef] [Google Scholar]
  75. Li, G., & Batygin, K. 2014, ApJ, 790, 69 [Google Scholar]
  76. Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [Google Scholar]
  77. Lin, Y., & Ogilvie, G. I. 2021, ApJ, 918, L21 [NASA ADS] [CrossRef] [Google Scholar]
  78. Lovelace, R. V. E., Covey, K. R., & Lloyd, J. P. 2011, AJ, 141, 51 [CrossRef] [Google Scholar]
  79. Lu, T., & Laughlin, G. 2022, Planet. Sci. J., 3, 221 [Google Scholar]
  80. Machida, M. N., Kokubo, E., Inutsuka, S.-I., & Matsumoto, T. 2008, ApJ, 685, 1220 [Google Scholar]
  81. Malamud, U., & Perets, H. B. 2020, MNRAS, 493, 698 [Google Scholar]
  82. Malhotra, R., & Dermott, S. F. 1990, Icarus, 85, 444 [NASA ADS] [CrossRef] [Google Scholar]
  83. Martin, R. G., Zhu, Z., & Armitage, P. J. 2020, ApJ, 898, L26 [Google Scholar]
  84. Millholland, S., & Batygin, K. 2019, ApJ, 876, 119 [Google Scholar]
  85. Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, Icarus, 219, 737 [NASA ADS] [CrossRef] [Google Scholar]
  86. Morgan, M., Seligman, D., & Batygin, K. 2021, ApJ, 917, L8 [NASA ADS] [CrossRef] [Google Scholar]
  87. Mosqueira, I., & Estrada, P. R. 2003a, Icarus, 163, 198 [NASA ADS] [CrossRef] [Google Scholar]
  88. Mosqueira, I., & Estrada, P. R. 2003b, Icarus, 163, 232 [NASA ADS] [CrossRef] [Google Scholar]
  89. Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge University Press) [Google Scholar]
  90. Néron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [Google Scholar]
  91. Nesvorný, D. 2015, AJ, 150, 73 [CrossRef] [Google Scholar]
  92. Nesvorný, D. 2018, ARA&A, 56, 137 [CrossRef] [Google Scholar]
  93. Nesvorný, D., & Morbidelli, A. 2012, AJ, 144, 117 [Google Scholar]
  94. Nesvorný, D., & Vokrouhlický, D. 2016, ApJ, 825, 94 [CrossRef] [Google Scholar]
  95. Nesvorný, D., Vokrouhlický, D., & Morbidelli, A. 2007, AJ, 133, 1962 [CrossRef] [Google Scholar]
  96. Nesvorný, D., Vokrouhlický, D., Bottke, W. F., & Levison, H. F. 2018, Nat. Astron., 2, 878 [Google Scholar]
  97. Nettelmann, N., Helled, R., Fortney, J. J., & Redmer, R. 2013, Planet. Space Sci., 77, 143 [Google Scholar]
  98. Neuenschwander, B. A., & Helled, R. 2022, MNRAS, 512, 3124 [Google Scholar]
  99. O’Connor, C. E., Teyssandier, J., & Lai, D. 2022, MNRAS, 513, 4178 [Google Scholar]
  100. Peale, S. J. 1969, AJ, 74, 483 [Google Scholar]
  101. Quillen, A. C., Chen, Y.-Y., Noyelles, B., & Loane, S. 2018, Celest. Mech. Dyn. Astron., 130, 11 [Google Scholar]
  102. Reinhardt, C., Chau, A., Stadel, J., & Helled, R. 2020, MNRAS, 492, 5336 [Google Scholar]
  103. Rogoszinski, Z., & Hamilton, D. P. 2020, ApJ, 888, 60 [Google Scholar]
  104. Rogoszinski, Z., & Hamilton, D. P. 2021, Planet. Sci. J., 2, 78 [Google Scholar]
  105. Rufu, R., & Canup, R. M. 2017, AJ, 154, 208 [NASA ADS] [CrossRef] [Google Scholar]
  106. Rufu, R., & Canup, R. M. 2022, ApJ, 928, 123 [NASA ADS] [CrossRef] [Google Scholar]
  107. Saillenfest, M., & Lari, G. 2021, A&A, 654, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2016, Celest. Mech. Dyn. Astron., 126, 369 [NASA ADS] [CrossRef] [Google Scholar]
  109. Saillenfest, M., Fouchard, M., Ito, T., & Higuchi, A. 2019a, A&A, 629, A95 [Google Scholar]
  110. Saillenfest, M., Laskar, J., & Boué, G. 2019b, A&A, 623, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Saillenfest, M., Lari, G., & Courtot, A. 2020, A&A, 640, A11 [EDP Sciences] [Google Scholar]
  112. Saillenfest, M., Lari, G., & Boué, G. 2021a, Nat. Astron., 5, 345 [Google Scholar]
  113. Saillenfest, M., Lari, G., Boué, G., & Courtot, A. 2021b, A&A, 647, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Salmon, J., & Canup, R. M. 2022, ApJ, 924, 6 [NASA ADS] [CrossRef] [Google Scholar]
  115. Sasaki, T., Stewart, G. R., & Ida, S. 2010, ApJ, 714, 1052 [Google Scholar]
  116. Serra, D., Lari, G., Tommei, G., et al. 2019, MNRAS, 490, 766 [Google Scholar]
  117. Slattery, W. L., Benz, W., & Cameron, A. G. W. 1992, Icarus, 99, 167 [NASA ADS] [CrossRef] [Google Scholar]
  118. Su, Y., & Lai, D. 2020, ApJ, 903, 7 [Google Scholar]
  119. Su, Y., & Lai, D. 2022, MNRAS, 509, 3301 [Google Scholar]
  120. Szulágyi, J., Cilibrasi, M., & Mayer, L. 2018, ApJ, 868, L13 [Google Scholar]
  121. Talu, T., Alessi, E. M., & Tommei, G. 2021, Universe, 7, 482 [Google Scholar]
  122. Tiscareno, M. S., Hedman, M. M., Burns, J. A., & Castillo-Rogez, J. 2013, ApJ, 765, L28 [NASA ADS] [CrossRef] [Google Scholar]
  123. Tittemore, W. C., & Wisdom, J. 1989, Icarus, 78, 63 [NASA ADS] [CrossRef] [Google Scholar]
  124. Tittemore, W. C., & Wisdom, J. 1990, Icarus, 85, 394 [NASA ADS] [CrossRef] [Google Scholar]
  125. Tremaine, S. 1991, Icarus, 89, 85 [Google Scholar]
  126. Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706 [Google Scholar]
  127. Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [Google Scholar]
  128. Vaillant, T., Laskar, J., Rambaux, N., & Gastineau, M. 2019, A&A, 622, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Verheylewegen, E., Noyelles, B., & Lemaitre, A. 2013, MNRAS, 435, 1776 [Google Scholar]
  130. Vokrouhlický, D., & Nesvorný, D. 2015, ApJ, 806, 143 [Google Scholar]
  131. Ward, W. R., & Canup, R. M. 2006, ApJ, 640, L91 [Google Scholar]
  132. Ward, W. R., & Canup, R. M. 2010, AJ, 140, 1168 [Google Scholar]
  133. Ward, W. R., & Hamilton, D. P. 2004, AJ, 128, 2501 [Google Scholar]
  134. Williams, J. G., & Boggs, D. H. 2016, Celest. Mech. Dyn. Astron., 126, 89 [Google Scholar]
  135. Woo, J. M. Y., Reinhardt, C., Cilibrasi, M., et al. 2022, Icarus, 375, 114842 [NASA ADS] [CrossRef] [Google Scholar]
  136. Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants, ed. T.J. Ahrens, 1 [Google Scholar]

1

There is currently no consensus in the literature on the exact definition of the “Laplace radius”. To avoid confusion with previous works, we stick to the definition of Tremaine et al. (2009) for rL and introduce rM as a distinct characteristic length.

2

Throughout this article, the orbit of the planet-satellite barycentre around the star is generally referred to as “the orbit of the planet” for simplicity. This distinction will be important in Sect. 3, when using a self-consistent coupled model for the planet and its satellite.

3

For short, the first-order secular spin–orbit resonance with resonant angle σ = ψ + ϕj will be called the “vj resonance”, where j is a given index in Table 1.

4

The destabilisation process is complicated by the existence of the stable eccentric equilibrium (see Tremaine et al. 2009; Saillenfest & Lari 2021). Yet, we have checked that when the destabilisation is triggered, the timescale of the eccentricity increase in our simulations is indeed consistent with that shown in Fig. 7.

5

The same dynamical mechanism has been reported by Saillenfest et al. (2019a) for trans-Neptunian objects perturbed by the galactic tides; in this case, rM ≈ 1000 au.

6

A fast satellite migration is also discarded by Rufu & Canup (2022), even though it could be a promising solution to their alternative version of the coaccretion + giant impact model; see their final discussion.

7

𝒢 is associated with the proper rotation (or spin motion), and with the Eulerian free nutation, also called Chandler wobble (from Chandler 1891 who observed it for the Earth). The period of the free wobble is expected to be a week or so for rigid Jupiter and Saturn, and about 50 days for Uranus and Neptune (as computed from their spin rates and moments of inertia).

8

As only the norm of e2 appears at quadrupolar order, it does not matter whether we consider the orbit of the star around the planet-satellite barycentre or the contrary.

9

The relaxation of planets to hydrostatic equilibrium tends to reduce the angle J, and specific internal processes are needed to maintain a non-zero value. The mean value of J is measured to be about 7 × 10−7 rad for the Earth (Chandler 1891) and 3 × 10−8 rad for Mars (Konopliv et al. 2020). The Chandler wobble has not been detected yet for gaseous planets, but from the latest Juno data, the constant J for Jupiter is constrained to be smaller than 10−6 rad and compatible with zero (Iess et al. 2018; Serra et al. 2019; Durante et al. 2020). Similar constraints have been obtained for Saturn (Iess et al. 2019).

All Tables

Table 1

First twenty terms of Uranus’s inclination and longitude of ascending node in the J2000 ecliptic and equinox reference frame.

Table 2

Minimum mass of the satellite in order for Uranus to reach the few closest resonances and be tilted from 0° to 90°.

Table 3

Properties of the secular resonances appearing in the orbital dynamics of a massless satellite at first order in the solar perturbation.

Table B.1

Quasi-periodic decomposition of Uranus’s eccentricity and longitude of perihelion (variable z).

Table B.2

Quasi-periodic decomposition of Uranus’s inclination and longitude of ascending node (variable ζ).

All Figures

thumbnail 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 = /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 . See Fig. 17 of Saillenfest & Lari (2021) for examples with other values of η.

In the text
thumbnail 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).

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail 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).

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail 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).

In the text
thumbnail 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).

In the text
thumbnail 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°).

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail 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 and its value is given by tan .

In the text
thumbnail 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 and , respectively. The Cassini states are labelled and showed with coloured dots. The thick black level is the separatrix of the resonance.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.