spin state of Uranus if its ancient satellite is single and has a mass smaller
than $m/M\approx 1.7\times 10^{-3}$. Taking into account a system of
additional inner satellites resembling the current ones would allow to
decrease this minimum limit by roughly a factor $2$. Below $m/M\approx
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\approx 4.4\times
10^{-4}$ at most.

###  4.3 Forming Uranus’ current satellites in a debris disc

Many previous works have been devoted to the formation of Uranus’ 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’ 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\approx 4.4\times 10^{-4}$. This minimum mass is
about four times the total mass of Uranus’ 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\times 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\approx 4.4\times 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 $m_{\mathrm{disc}}$, the characteristic
lifetime of the disc
$T_{\mathrm{disc}}=m_{\mathrm{disc}}/\dot{m}_{\mathrm{disc}}$ is

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

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

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

The current migration rate of Uranus’ 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; Ćuk 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$
Gyrs, 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 away999A 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.. However, these distances are in line with the
measurement of the current orbital expansion of Saturn’s satellites (Lainey et
al., 2020). In any case, if we assume that Uranus had a distant ancient
satellite subject to substantial tidal migration (which is our basic
hypothesis throughout this paper), then some mechanism of efficient energy
dissipation must exist in the interior of Uranus. Regardless of the precise
nature of this mechanism, the dissipation efficiency needed to move Oberon to
its current location (measured by the equivalent constant parameter $k_{2}/Q$;
see Sect. 2.3) is smaller than that needed to move the distant satellite over
$10$ $R_{\mathrm{eq}}$, so if we assume that the latter is possible, then the
former should follow.

Yet, unless the migration rate of Oberon is much larger than that of Titan
(which cannot be firmly ruled out yet, but may seem unlikely), we cannot
expect it to have migrated up its current position in less than a couple
billion years. This means that the end of the tilting mechanism of Uranus and
the destruction of its ancient satellite cannot be arbitrarily recent if one
wants to explain the formation of its current satellites in the debris disc
formed. For instance, one could imagine that the migration of the ancient
satellite and tilting of Uranus took $2$ Gyrs; 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$ Gyrs. 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’ 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 (Ćuk 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 in both
scenarios, and trying to link it to the mechanism discussed here 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’ equator. Hence, independently of whether the current
satellites of Uranus formed before or after the tilting, invoking a
substantial outward migration for them appears 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, like giant impacts, as being a main formation
ingredient. In this view, the small to moderate obliquity values of Jupiter
($3^{\circ}$), Saturn ($27^{\circ}$), and Neptune ($30^{\circ}$) 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^{\circ}$) is more challenging. Tilting Uranus without affecting
much its spin rate would require a slow process involving the spin-axis
precession motion. Such a process hardly fits the timespan offered by some
given stage 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^{\circ}$ to $80^{\circ}$ (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’ spin-axis precession rate has
frequency $\nu_{15}=g_{7}-g_{8}+s_{7}$ (where $g_{j}$ and $s_{j}$ are the
apsidal and nodal precession modes of the Solar System planets numbered from
$1$ for Mercury to $8$ for Neptune). In order to produce a resonance with
$\nu_{15}$, Uranus must have had a satellite with minimum mass $m/M=3.5\times
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 $\nu_{3}=s_{8}$. Triggering a resonance with $\nu_{3}$ and tilting
Uranus up to a large obliquity would require a satellite with minimum mass
$m/M=4.4\times 10^{-4}$, that is, $m\approx 4\times 10^{22}$ 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 $\varepsilon$ between
$0^{\circ}$ and $20^{\circ}$. 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 in 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\times 10^{-4}$ is roughly four times the total mass of the
current satellites of Uranus; this means that the satellite involved does not
exist any more and must have been destructed during the second phase of the
evolution. More massive ancient satellites allow for reaching other
resonances. The strongest resonance is $\nu_{2}=s_{7}$ which can be reached
for a satellite mass $m/M=2.2\times 10^{-3}$. Other promising resonance
candidates are $\nu_{7}=g_{5}-g_{7}+s_{7}$ and $\nu_{5}=-g_{5}+g_{6}+s_{6}$,
reachable for masses $m/M=1.3\times 10^{-3}$ and $1.7\times 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’ 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’ 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\omega+\Omega$ located at equatorial inclination $I_{\mathrm{Q}}\approx
56^{\circ}$, and to its symmetric counterpart $2\omega-\Omega$ at
$I_{\mathrm{Q}}\approx 124^{\circ}$. This resonance is well known for
artificial satellites (see e.g. Daquin et al., 2022). Since it vanishes for a
planet obliquity $\varepsilon=90^{\circ}$, satellite collisions can only occur
when the planet’s obliquity is substantially larger than $90^{\circ}$ (as
Uranus today), or substantially smaller than $90^{\circ}$. Interestingly, the
current state of Uranus is near the maximum of our histograms in the
$\varepsilon>90^{\circ}$ 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\times 10^{-4}$, the satellite
generally collides into the planet before the obliquity reaches $80^{\circ}$.
Hence, it seems hard to reproduce the current state of Uranus through a
capture in resonance $\nu_{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 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\times 10^{-3}$,
the probability of success through a capture in resonance $\nu_{2}$ can be as
large as $86\%$. For slightly smaller masses ($m/M\approx 1.5$ to $2\times
10^{-3}$), resonances $\nu_{5}$ and $\nu_{7}$ also demonstrate large success
ratios, that 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 $\nu_{2}$ is special, because it is a secular spin-orbit resonance
with Uranus’ own nodal precession mode $s_{7}$. 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 $\nu_{2}$ produces very strong kicks that can quickly pump
the obliquity up to Uranus’ current value. For this reason, our simulations
show that as long as the satellite is initialised at the right distance
($a_{0}\approx 50\,R_{\mathrm{eq}}$), the tilting mechanism via resonance
$\nu_{2}$ can operate down to very small migration ranges (only a few
$R_{\mathrm{eq}}$, instead of $10\,R_{\mathrm{eq}}$ in the general case). When
studying the formation of the main satellites of Jupiter, Saturn, and Uranus
in gaseous circumplanetary discs, Mosqueira & Estrada (2003a, b) found that a
distance of about $50\,R_{\mathrm{eq}}$ is precisely where one would expect a
massive satellite to form around Uranus. The current main satellites of Uranus
(Miranda, Ariel, Umbriel, Titania, and Oberon) are much further in and they
may seem peculiarly small as compared to other satellites of the giant
planets. From these argument, one may argue at least that forming a satellite
with mass ratio $m/M\approx 2\times 10^{-3}$ at a distance of $a_{0}\approx
50\,R_{\mathrm{eq}}$ does not seem too unreasonable in the context of
satellite formation theories.

Any successful tilting scenario must not only reproduce the current spin state
of Uranus, but also allow for the existence of its observed satellite system.
Today’s main satellites of Uranus have prograde orbits, and they show signs of
having been formed in a tidal disc (Crida & Charnoz, 2012). We discussed the
possibility that they formed from the debris disc produced by the destruction
of Uranus’ 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\approx 3.4\,R_{\mathrm{eq}}$) up to their current distance (Oberon is at
$a\approx 23\,R_{\mathrm{eq}}$). Such a large migration range appears
unrealistic in the context of traditional tidal models, but it could be
achievable on a billion-year timescale if a mechanism similar to the ‘tidal
lock’ recently proposed for Titan is also at play for the satellites of Uranus
(Fuller et al., 2016; Lainey et al., 2020).Yet, if this is the case, the
satellites would have no particular reason to follow the peculiar distribution
described by Crida & Charnoz (2012) as an evidence of their origin in a tidal
disc. The possible relationship between our proposed ancient satellite and the
current ones is therefore not a simple question.

Moreover, in our subset of simulations that successfully reproduce the spin
state of Uranus, the disruption of the ancient satellite almost always occur
when the latter is retrograde ($I_{\mathrm{Q}}\approx 124^{\circ}$).
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 details 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 did not investigate yet the effects of tidal dissipation within the
satellite. 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 the internal composition of Uranus’ 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\approx 4.4\times 10^{-4}$; _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
reproducing the current state of Uranus via the migration of an ancient
satellite. This scenario leaves room for refinements. One promising variant
would be to invoke two satellites: a big distant one (as that described here),
plus a small one located much closer to Uranus. In this configuration, the
inner satellite would be much more prone to tidal migration than the outer
one; and while the inner satellite migrates outwards over a few
$R_{\mathrm{eq}}$, the outer one would feel a change in the effective $J_{2}$
parameter of Uranus, which would increase its characteristic radius
$r_{\mathrm{M}}$. As a result, the whole tilting mechanism would be triggered
all the same, even though the distant big satellite would not need to migrate
at all. This variant scenario exemplifies how it would be possible to relax
constraints both on the large energy dissipation required in the interior of
Uranus (as the migration of close-in satellites requires much less energy
dissipation that 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 bulk of our successful runs peaks at
Uranus’ 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’
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.

## References

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

##  Appendix A Tilting Uranus during the planetary migration

If we restrict the precession of the planet’s orbital plane to a single
dominant harmonic with amplitude $S_{0}=\sin(I_{0}/2)$ and angle $\phi_{0}$
evolving at frequency $\nu_{0}$ (see Eq. 6), the Hamiltonian function
governing the secular spin-axis dynamics of the planet is simplified into

$\displaystyle\mathcal{H}$
$\displaystyle=-\frac{\alpha}{2}\frac{X^{2}}{(1-e_{\odot}^{2})^{3/2}}+\nu_{0}\Phi_{0}$
(13) $\displaystyle+\nu_{0}(1-\cos I_{0})X-\nu_{0}\sin
I_{0}\sqrt{1-X^{2}}\cos(\psi+\phi_{0})\,,$

where $\Phi_{0}$ is the momentum conjugate to $\phi_{0}$, added to obtain an
autonomous Hamiltonian system. In this expression, $\alpha$ is the precession
constant, $e_{\odot}$ is the eccentricity of the planet, $X$ is the cosine of
obliquity, and $\psi$ is the precession angle (see e.g. Saillenfest et al.
2019b). Passing to the resonant canonical coordinates $X$ and
$-\sigma=-\psi-\phi_{0}$ and dropping unnecessary constant parts, the
Hamiltonian function can be rewritten

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

using the rescaled timescale $\mathrm{d}\tau=p\mathrm{d}t$, where

$\gamma=-\frac{\nu_{0}\cos I_{0}}{p}\,,\quad\quad\beta=-\frac{\nu_{0}\sin
I_{0}}{p}\,,$ (15)

and $p=\alpha(1-e_{\odot}^{2})^{-3/2}$ is the characteristic spin-axis
precession frequency of the planet. Equation (14) 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 $\nu_{0}<0$, so that $\gamma$ and $\beta$ are both
positive with these notations. Under the hypothesis that Uranus’ obliquity is
small initially, the fastest tilting is obtained when Uranus acquires its
current obliquity during a single cycle of $\sigma$ (‘resonance kick’).
Assuming that the planet’s orbit is stable enough to consider that $\gamma$
and $\beta$ 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 $\varepsilon=0^{\circ}$ to more than
$\varepsilon=90^{\circ}$ first requires that

$\beta\geqslant|\gamma-1/2|\,,$ (16)

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. 14). Moreover, in
order for these two points to belong to the same trajectory, the separatrix
surrounding Cassini state 2 must either contain the North pole of the sphere
or be inexistent. As shown by Saillenfest et al. (2019b), this condition can
be written

$8\beta^{2}\geqslant 1-20\gamma-8\gamma^{2}+(1+8\gamma)^{3/2}\,.$ (17)

The two curves defined by Eqs. (16) and (17) connect at
$\gamma=\gamma_{\mathrm{crit}}$, where
$\gamma_{\mathrm{crit}}=(1+\sqrt{2})/4$. For $\gamma<\gamma_{\mathrm{crit}}$,
Eq. (17) is the more stringent condition, whereas for
$\gamma>\gamma_{\mathrm{crit}}$, Eq. (16) is the more stringent condition.
According to Eq. (15), $\gamma$ and $\beta$ are linked through
$\beta=\gamma\tan I_{0}$. Therefore, the two previous conditions can be
formulated as a minimum value for $I_{0}$ instead of a minimum value for
$\beta$. The minimum inclination of the planet allowing for its tilting from
$0^{\circ}$ to more than $90^{\circ}$ as a function of $\gamma$ is shown in
Fig. 21. The global minimum of this function is located at
$\gamma=\gamma_{\mathrm{crit}}$, and its value is given by $\tan
I_{\mathrm{min}}=3-2\sqrt{2}$, that is, $I_{\mathrm{min}}\approx
9.74^{\circ}$.

Figure 21: Minimum orbital inclination of a planet such that its obliquity can
evolve from $0^{\circ}$ to more than $90^{\circ}$ over a single libration in a
secular spin-orbit resonance. Conditions 1 and 2 are given in Eqs. (16) and
(17), respectively. The minimum inclination is reached at
$\gamma_{\mathrm{crit}}=(1+\sqrt{2})/4$ and its value is given by $\tan
I_{\mathrm{min}}=3-2\sqrt{2}$.

The phase portrait of the Hamiltonian at the critical point is shown in Fig.
22. We see that it corresponds to the case where the separatrix enclosing
Cassini state 2 goes exactly from $\varepsilon=0^{\circ}$ to
$\varepsilon=90^{\circ}$.

Figure 22: Phase portrait of Colombo’s top Hamiltonian when the planet has the
minimum orbital inclination required for going from an obliquity
$\varepsilon=0^{\circ}$ to $\varepsilon\geqslant 90^{\circ}$ 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 $\varepsilon$ is the
tilt from the $z$-axis and the resonance angle $\sigma$ is the polar angle
measured in the $xy$-plane. The parameters are $\gamma=(\sqrt{2}+1)/4$ and
$\beta=(\sqrt{2}-1)/4$, respectively. The Cassini states are labelled and
showed with coloured dots. The thick black level is the separatrix of the
resonance.

The problem is the same if Uranus’ spin axis resonates with its own orbital
precession (i.e. $\nu_{0}=s_{7}$) or with a harmonic stemming from an other
planet (e.g. Neptune, as investigated by Rogoszinski & Hamilton 2021, that is,
$\nu_{0}=s_{8}$): 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
$\gamma$ and $\beta$). As discussed in Sect. 1, this alternative strongly
increases the tilting timescale, because the drift must be adiabatic enough in
order for Uranus to remain inside the resonance and gradually follow its
drift.

##  Appendix B Orbital solution for Uranus

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

The orbital solution is expressed in the variables $z$ and $\zeta$ as
described in Eqs. (6) and (7). In Tables 4 and 5, we give all terms of the
solution in the J2000 ecliptic and equinox reference frame.

Because of the chaotic dynamics of the Solar System (Laskar 1989, 1990), the
fundamental frequencies related to the terrestrial planets (e.g. $s_{1}$,
$s_{2}$, and $\gamma$ appearing in Table 1) could vary noticeably over
billions of years (Hoang et al. 2021). However, they only marginally
contribute to Uranus’ 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 4: Quasi-periodic decomposition of Uranus’ eccentricity and longitude of
perihelion (variable $z$).

$\begin{array}[]{rrrr}\hline\cr\hline\cr k&\mu_{k}\
(^{\prime\prime}\,\text{yr}^{-1})&E_{k}\times 10^{9}&\theta_{k}^{(0)}\
(^{\text{o}})\\\ \hline\cr 1&4.24882&37351497&210.67\\\
2&3.08952&29110713&121.36\\\ 3&0.66708&1640609&73.98\\\
4&28.22069&1540208&308.11\\\ 5&1.93168&423140&208.85\\\
6&5.40817&411981&301.79\\\ 7&2.97706&75968&306.81\\\
8&27.06140&53352&218.76\\\ 9&4.89647&47827&115.82\\\
10&52.19257&39627&45.83\\\ 11&-20.88236&35859&203.93\\\
12&1.82121&28918&154.92\\\ 13&0.77840&28310&65.10\\\
14&-19.72306&27851&293.24\\\ 15&28.86795&10506&32.64\\\
16&27.57346&9228&43.74\\\ 17&-0.49216&8225&164.74\\\ 18&29.37998&7279&38.31\\\
19&53.35188&7234&314.93\\\ 20&3.60029&4594&301.39\\\
21&-56.90922&4109&224.11\\\ 22&5.59644&3787&110.35\\\
23&7.45592&3521&200.24\\\ 24&17.91550&2648&335.34\\\
25&76.16447&2387&143.03\\\ 26&5.47449&1426&275.01\\\
27&51.03334&1325&316.29\\\ 28&5.71670&1082&120.52\\\ 29&5.35823&582&94.89\\\
30&6.93423&536&168.90\\\ 31&7.05595&423&358.70\\\ 32&17.36469&409&303.95\\\
33&7.34103&297&207.85\\\ 34&4.36906&221&220.83\\\ 35&5.99227&214&113.56\\\
36&7.57299&207&11.47\\\ 37&5.65485&206&39.22\\\ 38&5.23841&171&272.97\\\
39&17.08266&114&359.38\\\ 40&16.81285&98&93.77\\\ 41&6.82468&97&194.53\\\
42&7.20563&61&143.91\\\ 43&17.63081&42&10.70\\\ 44&19.01870&41&39.75\\\
45&17.15752&39&145.02\\\ 46&7.71663&38&93.52\\\ 47&17.81084&33&238.56\\\
48&18.18553&30&237.28\\\ 49&17.72293&25&228.46\\\ 50&16.52731&24&311.91\\\
51&18.01611&21&224.83\\\ 52&17.47683&21&80.26\\\ 53&17.55234&17&17.65\\\
54&18.46794&15&183.15\\\ 55&16.26122&15&238.89\\\ 56&18.08627&12&176.17\\\
\hline\cr\end{array}$

101010This solution has been directly obtained from Laskar (1990) as explained
in the text. The phases $\theta_{k}^{(0)}$ are given at time J2000.

Table 5: Quasi-periodic decomposition of Uranus’ inclination and longitude of
ascending node (variable $\zeta$).

$\begin{array}[]{rrrr}\hline\cr\hline\cr k&\nu_{k}\
(^{\prime\prime}\,\text{yr}^{-1})&S_{k}\times 10^{9}&\phi_{k}^{(0)}\
(^{\text{o}})\\\ \hline\cr 1&0.00000&13773646&107.59\\\
2&-3.00557&8871413&320.33\\\ 3&-0.69189&563042&203.96\\\
4&-26.33023&347710&307.29\\\ 5&-2.35835&299979&224.75\\\
6&-4.16482&187859&231.66\\\ 7&-1.84625&182575&224.56\\\
8&-3.11725&59252&146.97\\\ 9&-1.19906&25881&313.99\\\
10&11.50319&18941&101.01\\\ 11&10.34389&11930&11.68\\\
12&-26.97744&10362&225.10\\\ 13&-5.61755&10270&348.70\\\
14&20.96631&7346&237.78\\\ 15&-0.58033&5474&197.32\\\
16&-5.50098&3662&342.89\\\ 17&-50.30212&2748&29.83\\\
18&-7.07963&2372&273.81\\\ 19&0.46547&1575&106.88\\\
20&82.77163&1514&308.95\\\ 21&-5.85017&1238&165.47\\\
22&-7.19493&1177&105.12\\\ 23&-6.96094&1120&97.96\\\
24&-17.74818&1109&303.28\\\ 25&-5.21610&927&18.91\\\ 26&-5.37178&922&35.48\\\
27&-5.10025&903&195.38\\\ 28&58.80017&887&32.90\\\ 29&-6.84091&826&109.96\\\
30&-7.33264&810&196.75\\\ 31&34.82788&701&114.12\\\ 32&0.57829&545&283.72\\\
33&-5.96899&529&350.64\\\ 34&-18.85115&463&240.24\\\ 35&-27.48935&449&38.53\\\
36&-25.17116&425&35.94\\\ 37&-6.15490&383&89.77\\\ 38&9.18847&310&181.15\\\
39&-28.13656&282&134.07\\\ 40&-6.73842&265&44.50\\\ 41&-7.40536&245&233.35\\\
42&-7.48780&228&47.95\\\ 43&-6.56016&214&303.47\\\ 44&-8.42342&178&211.21\\\
45&18.14984&139&111.19\\\ 46&-19.40256&122&28.16\\\ 47&-17.19656&119&153.99\\\
48&-18.01114&84&62.09\\\ 49&-17.66094&80&318.93\\\ 50&-17.83857&65&109.13\\\
51&-17.54636&58&66.71\\\ 52&-17.94404&41&32.26\\\ 53&-18.59563&37&278.11\\\
54&-19.13075&15&125.90\\\ 55&-18.30007&11&78.29\\\ 56&-18.97001&6&253.36\\\
57&-18.69743&3&41.70\\\ 58&-18.77933&3&42.83\\\ 59&-18.22681&3&226.30\\\
60&-19.06544&3&230.21\\\ \hline\cr\end{array}$

111111This solution has been directly obtained from Laskar (1990) as explained
in the text. The phases $\phi_{k}^{(0)}$ are given at time J2000.

##  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 other planets
in the system, in order for the secular spin-orbit resonances to also be
present. In this section, we describe the model that we use for this purpose.

###  C.1 Unaveraged Hamiltonian function

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

$\mathcal{I}=\begin{pmatrix}A&0&0\\\ 0&B&0\\\
0&0&C\end{pmatrix}_{(\mathbf{I},\mathbf{J},\mathbf{K})}\quad\text{where
}A\leqslant B\leqslant C\,.$ (18)

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

$\mathcal{H}=\mathcal{H}_{\mathrm{E}}+\mathcal{H}_{\mathrm{N}}+\mathcal{H}_{\mathrm{I}}\,,$
(19)

where $\mathcal{H}_{\mathrm{E}}$ corresponds to the Eulerian free rigid
rotation of the planet, $\mathcal{H}_{\mathrm{N}}$ describes the Newtonian
attraction of three mass points, and $\mathcal{H}_{\mathrm{I}}$ contains the
interactions of the non-spherical component of the planet with the satellite
and the star mass points. The Eulerian part $\mathcal{H}_{\mathrm{E}}$ is

$\mathcal{H}_{\mathrm{E}}=\frac{1}{2}\mathbf{G}^{\mathrm{T}}\mathcal{I}^{-1}\mathbf{G}=\frac{(\mathbf{G}\cdot\mathbf{I})^{2}}{2A}+\frac{(\mathbf{G}\cdot\mathbf{J})^{2}}{2B}+\frac{(\mathbf{G}\cdot\mathbf{K})^{2}}{2C}\,,$
(20)

where $\mathbf{G}$ is the rotational angular momentum of the planet. The
Newtonian part $\mathcal{H}_{\mathrm{N}}$ can be written

$\mathcal{H}_{\mathrm{N}}=\sum_{i=0}^{2}\frac{\|\mathbf{X}_{i}\|^{2}}{2m_{i}}-\sum_{0\leqslant
i<j\leqslant
2}\frac{\mathcal{G}m_{i}m_{j}}{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|}\,,$ (21)

where $\mathcal{G}$ is the gravitational constant. Because of the hierarchical
nature of the system, we use the following set of Jacobi coordinates:

$\left\\{\begin{aligned}
\mathbf{r}_{0}&=\frac{m_{0}\mathbf{x}_{0}+m_{1}\mathbf{x}_{1}+m_{2}\mathbf{x}_{2}}{m_{0}+m_{1}+m_{2}}\,,\\\
\mathbf{r}_{1}&=\mathbf{x}_{1}-\mathbf{x}_{0}\,,\\\
\mathbf{r}_{2}&=\mathbf{x}_{2}-\frac{m_{0}\mathbf{x}_{0}+m_{1}\mathbf{x}_{1}}{m_{0}+m_{1}}\,,\end{aligned}\right.$
(22)

with conjugate momenta

$\left\\{\begin{aligned}
\mathbf{p}_{0}&=\beta_{0}\dot{\mathbf{r}}_{0}=\mathbf{X}_{0}+\mathbf{X}_{1}+\mathbf{X}_{2}\,,\\\
\mathbf{p}_{1}&=\beta_{1}\dot{\mathbf{r}}_{1}=\mathbf{X}_{1}-m_{1}\frac{\mathbf{X}_{0}+\mathbf{X}_{1}}{m_{0}+m_{1}}\,,\\\
\mathbf{p}_{2}&=\beta_{2}\dot{\mathbf{r}}_{2}=\mathbf{X}_{2}-m_{2}\frac{\mathbf{X}_{0}+\mathbf{X}_{1}+\mathbf{X}_{2}}{m_{0}+m_{1}+m_{2}}\,,\end{aligned}\right.$
(23)

where the reduced masses are defined as

$\left\\{\begin{aligned} \beta_{0}&=m_{0}+m_{1}+m_{2}\,,\\\
\beta_{1}&=\frac{m_{0}m_{1}}{m_{0}+m_{1}}\,,\\\
\beta_{2}&=\frac{(m_{0}+m_{1})m_{2}}{m_{0}+m_{1}+m_{2}}\,.\end{aligned}\right.$
(24)

The vector $\mathbf{r}_{0}$ describes the position of the barycentre of the
whole system in the inertial reference frame, the vector $\mathbf{r}_{1}$ is
the position of the satellite with respect to the planet, and the vector
$\mathbf{r}_{2}$ 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

$\mathcal{H}_{\mathrm{N}}=\frac{\|\mathbf{p}_{0}\|^{2}}{2\beta_{0}}+\mathcal{H}_{\mathrm{K}}+\mathcal{H}_{\mathrm{M}}\,,$
(25)

where $\mathcal{H}_{\mathrm{K}}$ is the Hamiltonian describing two decoupled
Keplerian motions:

$\mathcal{H}_{\mathrm{K}}=\sum_{i=1}^{2}\left(\frac{\|\mathbf{p}_{i}\|^{2}}{2\beta_{i}}-\frac{\mu_{i}\beta_{i}}{\|\mathbf{r}_{i}\|}\right)=-\frac{\mu_{1}\beta_{1}}{2a_{1}}-\frac{\mu_{2}\beta_{2}}{2a_{2}}\,,$
(26)

in which we define $\mu_{1}=\mathcal{G}(m_{0}+m_{1})$ and
$\mu_{2}=\mathcal{G}(m_{0}+m_{1}+m_{2})$, and $a_{i}$ are the associated semi-
major axes. The Hamiltonian function $\mathcal{H}_{\mathrm{M}}$ contains the
mutual Newtonian perturbations:

$\mathcal{H}_{\mathrm{M}}=\left(\frac{\mathcal{G}(m_{0}+m_{1})m_{2}}{\|\mathbf{r}_{2}\|}-\frac{\mathcal{G}m_{0}m_{2}}{\|\mathbf{x}_{2}-\mathbf{x}_{0}\|}\right)-\frac{\mathcal{G}m_{1}m_{2}}{\|\mathbf{x}_{2}-\mathbf{x}_{1}\|}\,.$
(27)

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

$\displaystyle\mathbf{x}_{2}-\mathbf{x}_{0}$
$\displaystyle=\frac{m_{1}}{m_{0}+m_{1}}\mathbf{r}_{1}+\mathbf{r}_{2}\,,$ (28)
$\displaystyle\mathbf{x}_{2}-\mathbf{x}_{1}$
$\displaystyle=-\frac{m_{0}}{m_{0}+m_{1}}\mathbf{r}_{1}+\mathbf{r}_{2}\,,$

where $\|\mathbf{r}_{1}\|\ll\|\mathbf{r}_{2}\|$, and by performing a
multipolar expansion of $\mathcal{H}_{\mathrm{M}}$ using the Legendre
polynomials. This gives

$\mathcal{H}_{\mathrm{M}}=\frac{\mathcal{G}\beta_{1}m_{2}}{2r_{2}^{3}}\left(r_{1}^{2}-3\frac{(\mathbf{r}_{1}\cdot\mathbf{r}_{2})^{2}}{r_{2}^{2}}\right)+\mathcal{O}\left(\left[\frac{r_{1}}{r_{2}}\right]^{3}\right)\,,$
(29)

in which we write $r_{i}\equiv\|\mathbf{r}_{i}\|$.

The interaction part $\mathcal{H}_{\mathrm{I}}$ is obtained from the
gravitational potential of a mass element within the planet interacting with
the outer two point-mass bodies. Writing $\mathbf{y}$ the position of the mass
element measured with respect to the planet’s barycentre and
$\mathbf{y}_{i}\equiv\mathbf{x}_{i}-\mathbf{x}_{0}$, we again take advantage
of the hierarchical nature of the system by noting that
$\|\mathbf{y}\|\ll\|\mathbf{y}_{i}\|\,\forall i=1,2$ and that the planet is
almost spherical. By using a multipolar development of
$1/\|\mathbf{y}_{i}-\mathbf{y}\|$ and integrating over the planet’s volume, we
get

$\displaystyle\mathcal{H}_{\mathrm{I}}$
$\displaystyle=\sum_{i=1}^{2}\Bigg{[}\frac{\mathcal{G}m_{i}}{2\|\mathbf{y}_{i}\|^{5}}\bigg{[}(2A-B-C)\left(\mathbf{y}_{i}\cdot\mathbf{I}\right)^{2}+(2B-A-C)\left(\mathbf{y}_{i}\cdot\mathbf{J}\right)^{2}$
(30)
$\displaystyle+(2C-A-B)\left(\mathbf{y}_{i}\cdot\mathbf{K}\right)^{2}\bigg{]}+\mathcal{O}\left(\left[\frac{R_{\mathrm{eq}}}{\|\mathbf{y}_{i}\|}\right]^{3}\right)\Bigg{]}\,,$

where $R_{\mathrm{eq}}$ is the equatorial (i.e. largest) radius of the planet.
Replacing $\mathbf{y}_{1}$ and $\mathbf{y}_{2}$ by their expressions as a
function of $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, we finally obtain

$\displaystyle\mathcal{H}_{\mathrm{I}}$
$\displaystyle=\sum_{i=1}^{2}\frac{\mathcal{G}m_{i}}{2r_{i}^{5}}\bigg{[}(2A-B-C)\left(\mathbf{r}_{i}\cdot\mathbf{I}\right)^{2}+(2B-A-C)\left(\mathbf{r}_{i}\cdot\mathbf{J}\right)^{2}$
(31)
$\displaystyle+(2C-A-B)\left(\mathbf{r}_{i}\cdot\mathbf{K}\right)^{2}\bigg{]}$
$\displaystyle+\mathcal{O}\left(\left[\frac{R_{\mathrm{eq}}}{r_{1}}\right]^{3}\right)+\mathcal{O}\left(\left[\frac{R_{\mathrm{eq}}}{r_{2}}\right]^{3}\right)+\mathcal{O}\left(\left[\frac{r_{1}}{r_{2}}\right]^{3}\right)+\mathcal{O}\left(\frac{r_{1}}{r_{2}}\left[\frac{R_{\mathrm{eq}}}{r_{2}}\right]^{2}\right)\,.$

We notice that $\mathbf{r}_{0}$ does not appear in the Hamiltonian function,
thanks to the reduction of the total barycentre. As a result, the isolated
constant term $\|\mathbf{p}_{0}\|^{2}/(2\beta_{0})$ in Eq. (25) can be dropped
from the Hamiltonian function.

###  C.2 Secular system

The Hamiltonian function obtained above can be decomposed into

$\mathcal{H}=\mathcal{H}_{0}+\epsilon\mathcal{H}_{1}\,,$ (32)

where $\mathcal{H}_{0}=\mathcal{H}_{\mathrm{E}}+\mathcal{H}_{\mathrm{K}}$ is
the dominant integrable part and
$\epsilon\mathcal{H}_{1}=\mathcal{H}_{\mathrm{M}}+\mathcal{H}_{\mathrm{I}}$ is
a perturbation (we introduce a factor $\epsilon\ll 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 $\epsilon$ is given by the classic averaging method.

In order to average $\epsilon\mathcal{H}_{1}$ over the fast angles of the
rotational motion of the planet, one would need to express
$\mathcal{H}_{\mathrm{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 $g$ and $\ell$ circulate in an independent manner (i.e. the
motion is composed of a main rotation about the constant rotational angular
momentum $\mathbf{G}$ and a slower rotation about the principal axis
$\mathbf{K}$). Assuming that none of the two corresponding fast
frequencies121212$g$ is associated with the proper rotation (or spin motion),
and $\ell$ 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 of moments of
inertia). is involved in a resonance, we can therefore bypass the standard
procedure and average the Hamiltonian function over $g$ and $\ell$ separately
(see Boué & Laskar 2006; Vaillant et al. 2019). Expressing the vectors
$\mathbf{I}$, $\mathbf{J}$, and $\mathbf{K}$ in terms of the Andoyer angles,
we obtain

$\displaystyle\left<(\mathbf{r}_{i}\cdot\mathbf{I})^{2}\right>_{g,\ell}$
$\displaystyle=\left<(\mathbf{r}_{i}\cdot\mathbf{J})^{2}\right>_{g,\ell}$ (33)
$\displaystyle=\frac{1+\cos^{2}J}{4}r_{i}^{2}-\frac{1}{2}\left(1-\frac{3}{2}\sin^{2}J\right)(\mathbf{r}_{i}\cdot\mathbf{w})^{2}\,,$
$\displaystyle\left<(\mathbf{r}_{i}\cdot\mathbf{K})^{2}\right>_{g,\ell}$
$\displaystyle=\frac{\sin^{2}J}{2}r_{i}^{2}+\left(1-\frac{3}{2}\sin^{2}J\right)(\mathbf{r}_{i}\cdot\mathbf{w})^{2}\,,$

where $\mathbf{w}=\mathbf{G}/G$ and $J$ is the angle between $\mathbf{G}$ and
$\mathbf{K}$, such that

$\left<\mathcal{H}_{\mathrm{I}}\right>_{g,\ell}=-\frac{2C-A-B}{2}\left(1-\frac{3}{2}\sin^{2}J\right)\sum_{i=1}^{2}\frac{\mathcal{G}m_{i}}{2r_{i}^{3}}\bigg{[}1-3\frac{(\mathbf{r}_{i}\cdot\mathbf{w})^{2}}{r_{i}^{2}}\bigg{]}\,.$
(34)

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

It remains to average the perturbation Hamiltonian over the fast orbital
angles, which are the mean anomalies $M_{1}$ and $M_{2}$. 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

$\bar{\mathcal{H}}=\bar{\mathcal{H}}_{\mathrm{M}}+\bar{\mathcal{H}}_{\mathrm{I}}\,,$
(35)

where

$\bar{\mathcal{H}}_{\mathrm{M}}=\frac{\gamma}{6}\,\bigg{[}1-6e_{1}^{2}-3(1-e_{1}^{2})(\mathbf{k}_{1}\cdot\mathbf{k}_{2})^{2}+15(\mathbf{e}_{1}\cdot\mathbf{k}_{2})^{2}\bigg{]}\,,$
(36)

and

$\displaystyle\bar{\mathcal{H}}_{\mathrm{I}}$
$\displaystyle=\sum_{i=1}^{2}\frac{\alpha_{i}}{6}\,\bigg{[}1-3(\mathbf{k}_{i}\cdot\mathbf{w})^{2}\bigg{]}\,,$
(37)

in which we have defined

$\gamma=\frac{3\mathcal{G}\beta_{1}m_{2}a_{1}^{2}}{4a_{2}^{3}(1-e_{2}^{2})^{3/2}}\,,$
(38)

and

$\alpha_{i}=\frac{3\mathcal{G}m_{0}m_{i}R_{\mathrm{eq}}^{2}J_{2}}{2a_{i}^{3}(1-e_{i}^{2})^{3/2}}\left(1-\frac{3}{2}\sin^{2}J\right)\,,$
(39)

where

$J_{2}=\frac{2C-A-B}{2m_{0}R_{\mathrm{eq}}^{2}}\,.$ (40)

In these expressions, $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$ are the
eccentricity vectors of the satellite around the planet and of the planet-
satellite barycentre around the star131313As only the norm of $\mathbf{e}_{2}$
appears at quadrupolar order, it does not matter whether we consider the orbit
of the star around the planet-satellite barycentre or the contrary. (i.e. the
vector pointing to the pericentre and whose norm is the eccentricity), and
$\mathbf{k}_{1}$ and $\mathbf{k}_{2}$ 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 $\mathbf{G}$, $\mathbf{e}_{1}$, and
$\mathbf{G}_{1}$, where $\mathbf{G}_{1}$ is the orbital angular momentum of
the satellite:

$\mathbf{G}_{1}=\beta_{1}\sqrt{\mu_{1}a_{1}(1-e_{1}^{2})}\,\mathbf{k}_{1}\,.$
(41)

We obtain

$\displaystyle\dot{\mathbf{G}}$
$\displaystyle=-\sum_{i=1}^{2}\alpha_{i}\,(\mathbf{k}_{i}\cdot\mathbf{w})(\mathbf{k}_{i}\times\mathbf{w})\,,$
(42) $\displaystyle\dot{\mathbf{e}}_{1}$
$\displaystyle=\frac{\alpha_{1}}{2\|\mathbf{G}_{1}\|}\bigg{[}2(\mathbf{k}_{1}\cdot\mathbf{w})(\mathbf{e}_{1}\times\mathbf{w})-\Big{(}1-5(\mathbf{k}_{1}\cdot\mathbf{w})^{2}\Big{)}(\mathbf{k}_{1}\times\mathbf{e}_{1})\bigg{]}$
$\displaystyle+\frac{\gamma\,(1-e_{1}^{2})}{\|\mathbf{G}_{1}\|}\bigg{[}2(\mathbf{k}_{1}\times\mathbf{e}_{1})+(\mathbf{k}_{1}\cdot\mathbf{k}_{2})(\mathbf{e}_{1}\times\mathbf{k}_{2})$
$\displaystyle\hskip
142.26378pt-5(\mathbf{e}_{1}\cdot\mathbf{k}_{2})(\mathbf{k}_{1}\times\mathbf{k}_{2})\bigg{]}\,,$
$\displaystyle\dot{\mathbf{G}}_{1}$
$\displaystyle=\alpha_{1}\,(\mathbf{k}_{1}\cdot\mathbf{w})(\mathbf{k}_{1}\times\mathbf{w})+\gamma\,(1-e_{1}^{2})(\mathbf{k}_{1}\cdot\mathbf{k}_{2})(\mathbf{k}_{1}\times\mathbf{k}_{2})$
$\displaystyle-5\gamma\,(\mathbf{e}_{1}\cdot\mathbf{k}_{2})(\mathbf{e}_{1}\times\mathbf{k}_{2})\,.$

The elements $\mathbf{e}_{2}$ and $\mathbf{k}_{2}$ 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 $\mathbf{e}_{2}$ and $\mathbf{k}_{2}$ 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
$\mathbf{e}_{2}$ and $\mathbf{k}_{2}$ from the secular orbital solution of
Laskar (1990) described in Appendix B. The conversion from the variables $z$
and $\zeta$ to the vectors $\mathbf{e}_{2}$ and $\mathbf{k}_{2}$ 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’ rotational angular momentum $G$ is quite
uncertain, as both its spin rate $\omega$ and moment of inertia are poorly
known (see the discussion in Sect 2.2). For this reason, it is enough to use
the approximate expression for a principal axis rotator $G\approx
C\omega=\lambda\omega\,m_{0}R_{\mathrm{eq}}^{2}$, where $\lambda$ is the
normalised polar moment of inertia. We adopt the nominal value of
$\lambda\omega$ described in Sect. 2.2. We also neglect the $\sin^{2}J$ term
in Eq. (39) which is not known for Uranus but expected to be extremely
small141414The 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\times
10^{-7}$ rad for the Earth (Chandler 1891) and $3\times 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).. The migration of the satellite
is modelled by applying a slow drift in $a_{1}$ 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’ spin rate $\omega$, whose
contribution is smaller than our uncertainty on the value of $G$.

##  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
$\sigma_{k}=\psi+\phi_{k}$, where $\psi$ is the spin-axis precession angle and
$\phi_{k}$ is the circulating angle of the $k$th 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 $\sigma_{j,k}=2\psi+\phi_{j}+\phi_{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 $\sigma_{2}$. And indeed, we found
that some fraction of the trajectories shown in Fig. 10 for satellite masses
in the neighbourhood of $m=m_{2}$ is produced by stable captures in second-
order resonances. Figure 23 shows an example of such trajectories:
$\sigma_{2}=\psi+\phi_{2}$ first starts to oscillate, but it is soon replaced
by $\sigma_{2,5}=2\psi+\phi_{2}+\phi_{5}$, while both $\sigma_{2}$ and
$\sigma_{5}$ circulate. The nearby resonance $\sigma_{8}$ does not play any
role in this example.

Figure 23: 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\times 10^{-3}$. The satellite is initialised close
to the local circular Laplace equilibrium at $a=42$ $R_{\mathrm{eq}}$ and it
is made migrating outwards over $\Delta a=9$ $R_{\mathrm{eq}}$.

##  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
$\mathcal{H}=k_{\mathrm{P}}\mathcal{H}_{\mathrm{P}}+k_{\odot}\mathcal{H}_{\odot}$.
Expanding both contributions to quadrupole order, Saillenfest & Lari (2021)
define the two constant coefficients as

$k_{\mathrm{P}}=\frac{3}{4}\frac{\mathcal{G}M}{a}J_{2}\frac{R_{\mathrm{eq}}^{2}}{a^{2}}\quad\text{and}\quad
k_{\odot}=\frac{3}{8}\frac{\mathcal{G}m_{\odot}}{a_{\odot}(1-e_{\odot}^{2})^{3/2}}\frac{a^{2}}{a_{\odot}^{2}}\,,$
(43)

and the two Hamiltonian functions as

$\mathcal{H}_{\mathrm{P}}=\frac{1-3\cos^{2}I_{\mathrm{Q}}}{3(1-e^{2})^{3/2}}\,,$
(44)

and

$\displaystyle\mathcal{H}_{\odot}$
$\displaystyle=-\frac{1}{8}\bigg{[}8e^{2}+2(3e^{2}+2)(2C^{2}\cos^{2}I_{\mathrm{Q}}+S^{2}\sin^{2}I_{\mathrm{Q}}-2)$
(45) $\displaystyle\quad\quad+8CS(3e^{2}+2)\cos I_{\mathrm{Q}}\sin
I_{\mathrm{Q}}\cos\Omega_{\mathrm{Q}}$
$\displaystyle\quad\quad+5S^{2}e^{2}(\cos
I_{\mathrm{Q}}+1)^{2}\cos(2\omega_{\mathrm{Q}}+2\Omega_{\mathrm{Q}})$
$\displaystyle\quad\quad-20CSe^{2}(\cos I_{\mathrm{Q}}+1)\sin
I_{\mathrm{Q}}\cos(2\omega_{\mathrm{Q}}+\Omega_{\mathrm{Q}})$
$\displaystyle\quad\quad+10(3C^{2}-1)e^{2}\sin^{2}I_{\mathrm{Q}}\cos(2\omega_{\mathrm{Q}})$
$\displaystyle\quad\quad-20CSe^{2}(\cos I_{\mathrm{Q}}-1)\sin
I_{\mathrm{Q}}\cos(2\omega_{\mathrm{Q}}-\Omega_{\mathrm{Q}})$
$\displaystyle\quad\quad+5S^{2}e^{2}(\cos
I_{\mathrm{Q}}-1)^{2}\cos(2\omega_{\mathrm{Q}}-2\Omega_{\mathrm{Q}})$
$\displaystyle\quad\quad+2S^{2}(3e^{2}+2)\sin^{2}I_{\mathrm{Q}}\cos(2\Omega_{\mathrm{Q}})\bigg{]}\,.$

In this expression, $C=\cos\varepsilon$ and $S=\sin\varepsilon$ are constant,
and the orbital angles of the satellite
$(I_{\mathrm{Q}},\omega_{\mathrm{Q}},\Omega_{\mathrm{Q}})$ are measured with
respect to the equator and equinox of the planet. All other quantities that
appear in the Hamiltonian function are described in the main text.

Because of the factor $(1-e^{2})^{-3/2}$ in Eq. (44),
$k_{\mathrm{P}}\mathcal{H}_{\mathrm{P}}$ is the dominant contribution to the
dynamics whenever the eccentricity grows large. We consider the intermediate
regime in which $k_{\mathrm{P}}\mathcal{H}_{\mathrm{P}}$ is the dominant
integrable part of the Hamiltonian and $k_{\odot}\mathcal{H}_{\odot}$ is a
small perturbation. Written in terms of the Delaunay elements,
$k_{\mathrm{P}}\mathcal{H}_{\mathrm{P}}$ 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.
(45); for each of them, the resonance angle $\xi$ is a linear combination of
$\omega_{\mathrm{Q}}$ and $\Omega_{\mathrm{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 $\xi$ is first
taken as coordinate by a canonical change of variables
$(\omega_{\mathrm{Q}},\Omega_{\mathrm{Q}})\rightarrow(\xi,\gamma)$. The
conjugate momenta of $\xi$ and $\gamma$ are written $\Xi$ and $\Gamma$. In the
vicinity of the resonance, $\xi$ is a slow angle and $\gamma$ is a fast angle.
At first order in the perturbation, the long-term dynamics is described by
averaging the Hamiltonian function over $\gamma$. The Hamiltonian system
obtained is independent of $\gamma$, which means that $\Gamma$ 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, $\Gamma$ 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,\xi)$ plane or in the
$(I_{\mathrm{Q}},\xi)$ plane.

When it emerges from chaos and reaches the weakly perturbed regime, the
satellite can have any value of $K$. For definitiveness, we consider here for
each resonance the value of $K$ such that $e=0$ for $I_{\mathrm{Q}}=I_{0}$,
where $I_{0}$ is the nominal inclination of the resonance (see Table 3). Even
though this represents a loss of generality, we will see that studying just
one value of $K$ is enough to understand the correct hierarchy between the
various resonances. The case $\xi=\Omega_{\mathrm{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 $I_{\mathrm{Q}}$
along the separatrix are given by

$\cos^{2}I_{\mathrm{Q}}=\frac{1+u-\sqrt{1+u^{2}+2u\cos(2\varepsilon)}}{2u}\,,$
(46)

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

Figure 24: 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
$I_{\mathrm{Q}}$ of the satellite are linked through a constant quantity $K$.
The origin $e=0$ of the graph corresponds to the nominal inclination $I_{0}$
of the resonance.

##  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 $\mathbf{L}$ of the
planet and its satellite can be considered constant during the collision. The
constant vector $\mathbf{L}$ can be written

$\mathbf{L}=\mathbf{L}_{\mathrm{P}}+\mathbf{L}_{\mathrm{S}}\,,$ (47)

where $\mathbf{L}_{\mathrm{P}}$ is the contribution of the planet’s spin and
$\mathbf{L}_{\mathrm{S}}$ is the contribution of the satellite’s orbit. After
the removal of the satellite, the obliquity of Uranus is not expected to vary
at all due to external perturbations because Uranus is released far from any
kind of spin-orbit resonance. The magnitude of $\mathbf{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=\|\mathbf{L}\|$, which is much smaller than our uncertainty on
Uranus’ moment of inertia. To leading order, the relative change in the
planet’s spin rate is

$\frac{\Delta\omega}{\omega}\approx\frac{L-L_{\mathrm{P}}}{L_{\mathrm{P}}}\approx\frac{L_{\mathrm{S}}}{L_{\mathrm{P}}}\cos
I_{\mathrm{Q}}\,,$ (48)

where $L_{\mathrm{P}}=\|\mathbf{L}_{\mathrm{P}}\|$,
$L_{\mathrm{S}}=\|\mathbf{L}_{\mathrm{S}}\|$, and $I_{\mathrm{Q}}$ is the
orbital inclination of the satellite with respect to the planet’s equator at
the collision time (i.e. the angle between $\mathbf{L}_{\mathrm{P}}$ and
$\mathbf{L}_{\mathrm{S}}$). Using the gyroscopic approximation for the spin of
the planet, the ratio $L_{\mathrm{S}}/L_{\mathrm{P}}$ can be written as

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

where $n$ and $e$ are the satellite’s mean motion and eccentricity, and other
quantities are described in Sect. 2.1. The collision is expected to occur when
the semi-major axis of the satellite is $a\approx r_{\mathrm{M}}$ as defined
in Eq. (2). Moreover, the collision condition requires that $a(1-e)\lesssim
R_{\mathrm{eq}}$ which gives a minimum bound for the eccentricity value at the
collision time. When using in Eq. (49) the physical parameters of Uranus, 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 $\omega\lambda$ as
in Sect. 2.2. The order of magnitude for the relative change in the planet’s
spin rate is therefore

$\frac{\Delta\omega}{\omega}\approx 35\frac{m}{M}\cos I_{\mathrm{Q}}\,.$ (50)

The variation in the planet’s spin-axis orientation due to the impact is
quantified by the angle $\varphi$ between $\mathbf{L}$ and
$\mathbf{L}_{\mathrm{P}}$. From Eq. (47), it can be expressed as

$\sin\varphi=\frac{L_{\mathrm{S}}}{L}\sin I_{\mathrm{Q}}\,,$ (51)

that is,

$\sin\varphi\approx 35\frac{m}{M}\sin I_{\mathrm{Q}}\,.$ (52)

Under favourable longitude of node of the satellite at the time of collision,
the spin-axis deviation $\varphi$ directly equals the obliquity variation
$\Delta\varepsilon$. As shown in Sect. 3.4, satellite collisions only occur
when $I_{\mathrm{Q}}\approx 55^{\circ}$ or $I_{\mathrm{Q}}\approx
125^{\circ}$. 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 $\varphi$ would be obtained for a grazing
impact with $I_{\mathrm{Q}}=90^{\circ}$. Hence, an upper bound for the
obliquity change in radians due to the satellite collision is
$\Delta\varepsilon\approx 35\,m/M$.