# The large obliquity of Saturn explained by the fast migration of Titan

Melaine Saillenfest, Giacomo Lari, Gwenaël Boué

► To cite this version:

Melaine Saillenfest, Giacomo Lari, Gwenaël Boué. The large obliquity of Saturn explained by the fast migration of Titan. Nature Astronomy, 2021, 5, pp.345-349. 10.1038/s41550-020-01284-x. hal-03256804

HAL Id: hal-03256804

https://hal.science/hal-03256804v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The large obliquity of Saturn explained by the fast migration of Titan

Melaine Saillenfest<fnref n="1" />✉, Giacomo Lari² and Gwenaël Boué<fnref n="1" />

The obliquity of a planet is the tilt between its equator and its orbital plane. Giant planets are expected to form with near-zero obliquities<fnref n="1" />,². After the formation of Saturn, some dynamical mechanism must therefore have tilted Saturn up to its current obliquity of 26.7°. This event is traditionally thought to have happened more than 4 Gyr ago during the late planetary migration³⁻⁵ because of the crossing of a resonance between the spin-axis precession of Saturn and the nodal orbital precession mode of Neptune⁶. Here, we show that the fast tidal migration of Titan for which the measurement is reported in ref. ⁷ is incompatible with this scenario, and that it offers a new explanation for Saturn's current obliquity. A substantial migration of Titan would prevent any early resonance, which would invalidate previous constraints on the late planetary migration that were set by the tilting of Saturn⁸⁻¹⁰. We propose instead that the resonance was encountered more recently, about 1 Gyr ago, and forced Saturn's obliquity to increase from a small value (possibly less than 3°) to its current state. This scenario suggests that Saturn's normalized polar moment of inertia lies between 0.224 and 0.237. Our findings bring out a new paradigm for the spin-axis evolution of Saturn, Jupiter¹¹ and possibly giant exoplanets in multiple systems, whereby obliquities are not settled once for all but evolve continuously as a result of the migration of their satellites.

We investigate whether the spin-axis dynamics of Saturn could have been influenced by the migration of its satellites. The torque applied by the Sun on Saturn's equatorial bulge produces a precession of its spin axis with a mean frequency of $\dot{\psi} = \alpha(1 - e^2)^{-3/2} \cos \epsilon$, where $e$ is the eccentricity of Saturn, $\epsilon$ is its obliquity, and $\alpha$ is called its precession constant¹². The value of $\alpha$ incorporates the extra torques produced by Saturn's satellites<fnref n="1" />,¹³. It also depends on several orbital and physical parameters of Saturn that are all well known except Saturn's normalized polar moment of inertia $\lambda$. Estimates found in the literature have a broad range of $\lambda \in [0.200, 0.240]$, but they are model-dependent and a proper uncertainty range is difficult to define (Methods). This whole interval gives a current precession constant $\alpha_0$ that ranges from 0.747 to 0.894" yr⁻<fnref n="1" />.

We take into account the migration of Saturn's satellites by changing their contributions to $\alpha$. The measurements reported in ref. ⁷ unveiled a fast migration for Titan and a similar tidal timescale $t_{\text{tide}} = a/\dot{a}$ for all six of the satellites of Saturn that were studied. Both these findings support the 'resonance locking' tidal theory of ref. ¹⁴. With this theory, an approximate formula was obtained⁷ for the long-term trend of Titan's semi-major axis $a$ over time, in which Titan's tidal timescale evolves in the long run with a current value of $t_{\text{tide}} \approx 13$ Gyr, in accordance with observations. We use this formula to estimate the variation of $\alpha$ over time, and as Titan is the satellite

that contributes by far the most, we include other satellites through a slight increase in Titan's mass (Methods). Titan's mean semi-major axis is currently $a=20.275$ radii of Saturn ($R_{\text{eq}}$), but ref. ⁷ suggests that during the late planetary migration, more than 4 Gyr ago, $a$ was smaller than $10 R_{\text{eq}}$. As illustrated in Fig. 1, this variation translates into a large-scale drift of Saturn's precession constant $\alpha$. Therefore, contrary to previous assumptions⁸⁻¹⁰, we find that $\alpha$ cannot be considered as roughly constant over time.

A planet undergoes a secular spin-orbit resonance when there is a commensurability between the precession frequency of its spin axis and the frequency of one or several harmonics that appear in its orbital precession¹²,¹⁵,¹⁶. The orbital motion of the giant planets is now very steady¹⁷, such that the nodal orbital precession mode of Neptune, $s_8 \approx -0.6919"$ yr⁻<fnref n="1" />, is almost constant. By contrast, using the nominal migration law of Titan, values of $\alpha$ range between 0.314 and 0.376" yr⁻<fnref n="1" /> after the end of the late planetary migration (Fig. 1). These values are much smaller than $s_8$ or any other orbital frequency of the planets¹⁷, which prevents Saturn from having an early encounter with a low-order secular spin-orbit resonance.

This result is not specific to the migration law used for Titan. It arises for any realistic law that is compatible with the measurements reported in ref. ⁷, unless we are currently witnessing a transient episode of exceptionally fast migration for Titan and coincidentally similar $t_{\text{tide}}$ for most major satellites of Saturn, which appears unlikely. Moreover, as the planets were more tightly packed during their late migration, their orbital precession frequencies were all higher than those of today⁹, whereas $\alpha$ was not much affected by the slight change in Saturn's semi-major axis. This stretches even further the mismatch between the two kinds of precession frequencies. The ring of Saturn, even if it was more massive than it is today¹⁸, did not contribute substantially, and a massive circumplanetary disc could have produced only a small obliquity kick (Methods). Consequently, Saturn's obliquity remained essentially unaltered at that time. Provided that Titan did migrate substantially before today (that is, over a few $R_{\text{eq}}$ or more), this result would invalidate previous constraints on the late planetary migration that were set by the tilting of Saturn; the two events are probably unrelated.

Instead, the fast migration of Titan implies that Saturn reached the resonance with $s_8$ much more recently, as a result of the increase in $\alpha$. To explore possible trajectories of Saturn's spin axis, we propagate it numerically backwards in time for $\lambda$ finely sampled in its exploration interval. We use the same setting as in ref. ¹¹, namely, the orbit of Saturn evolves according to the full series provided in ref. ¹⁷ and the precession constant $\alpha$ is taken as a slowly varying parameter. Numerical experiments show that the system is not chaotic enough to produce noticeable numerical irreversibility, such that backward propagation does reveal the only dynamical pathway that corresponds to a given state of Saturn today. We limit our

<footnote n="1">¹IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Université de Lille, Paris, France. ²Department of Mathematics, University of Pisa, Pisa, Italy. ✉e-mail: melaine.saillenfest@obspm.fr</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 1** | Evolution of the effective precession constant $\alpha$ of Saturn due to the tidal migration of Titan. For a given value of $\lambda$, $\alpha$ is propagated through time by using Titan's nominal migration law given in ref. 7. The green error bar shows the current value $\alpha_0$ obtained by considering Saturn's normalized polar moment of inertia $\lambda$ lying in [0.200, 0.240]. The blue line shows the value of the nodal precession mode of Neptune, $s_8$; it has been nearly constant since the end of the late planetary migration$^{17}$, and was higher before. The spin-axis precession rate of Saturn scales as $\alpha \cos \varepsilon$.

propagations to 4 Gyr in the past, that is, confidently after the late planetary migration. We find that the resonance with $s_8$ was reached about 1 Gyr ago (Fig. 2). Before this event, the obliquity of Saturn oscillated with a small amplitude around a fixed value. The crossing of the thin resonance with $g_7 - g_8 + s_7$ hardly affected it at all. During the late planetary migration that occurred before our 4 Gyr timespan, the planetary precession frequencies were all higher than those of today, which resulted in an even steadier obliquity value.

Distinct dynamical pathways could have been followed by Saturn according to its current precession constant $\alpha_0$, whose value depends on Saturn's polar moment of inertia $\lambda$. Values of $\lambda$ ranging between about 0.224 and 0.237 imply that Saturn currently oscillates with a small amplitude inside the resonance. In this case, when Saturn's precession constant became greater than $\alpha \approx 0.696'' \text{yr}^{-1}$, the separatrices of the resonance appeared around its trajectory, which led to capture with 100% certainty (Fig. 2a). We expect this kind of evolution if Saturn reached the resonance with a small obliquity$^1$. Other values of $\lambda$ imply either that Saturn was captured recently with an already large obliquity (Fig. 2b), or that Saturn went past the resonance through its hyperbolic point (Fig. 2c). In our summary of Saturn's past obliquity values (Fig. 3), the U-shaped portion of the curve corresponds to values of $\alpha_0$ that put Saturn today inside the resonance with $s_8$. If Saturn is currently out of resonance, previous studies (even though they relied on a different scenario) also reported the impossibility (or negligible probability) of producing Saturn's current obliquity from an initially low value$^{6,8-10}$. The bulk of successful runs in previous studies correspond to values of $\alpha_0$ that put Saturn today in the vicinity of the resonance centre, which is also what we observe here in the context of Titan's tidal migration. Using the nominal migration law of Titan, a minimum past obliquity of Saturn equal to about 2.5° is reached for $\alpha_0 \approx 0.775'' \text{yr}^{-1}$ (that is, $\lambda \approx 0.231$) (Fig. 3).

The current orientation of Saturn's spin axis implies that if $\alpha$ evolved in a purely adiabatic way (that is, very slowly with respect to the obliquity oscillations), then Saturn's primordial obliquity could not have been exactly zero$^1$. However, the migration rate of Titan reported in ref. 7 is so fast that it departs somewhat from adiabaticity,

and one could expect that different migration laws would lead to different outcomes. The influence of Titan's migration law can be investigated by varying its tidal timescale $t_{\text{tide}} = a/\dot{a}$. The nominal

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

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

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

**Fig. 2** | Examples of possible obliquity evolution for Saturn. **a–c.** Each graph shows a 4 Gyr numerical trajectory of Saturn (black dots) computed for $\alpha$ evolving as in Fig. 1. Saturn lies today along the green error bar (the same bar as in Fig. 1 at $t=0$), and the value of $\lambda$ used is provided in each graph. The pink areas show all first-order secular spin-orbit resonances$^{16}$; the large area on the right is the resonance with $s_8$ and the thin area on the left is the resonance with $g_7 - g_8 + s_7$ (ref. 17). The separatrices of the $s_8$ resonance are highlighted in blue. If we assume that Saturn's primordial precession angle $\psi$ is a random number distributed uniformly between 0 and $2\pi$, trajectories that encounter the separatrix would have a given probability of being captured (**b**) or of crossing over the resonance (**c**). A given value of $\lambda$, however, corresponds to one and only one dynamical pathway such as those shown in these examples.
<figure><img src="image_7.png" /><figcaption>**Fig. 3 | Past obliquity of Saturn as a function of its current precession constant $α_0$.** For this graph, Titan migrates according to its nominal migration law as illustrated in Fig. 1. Each point of the black curve is obtained from a numerical simulation (see Fig. 2 for examples). The green line shows Saturn's current obliquity and the black dotted line highlights the 5° level.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 4 | Past obliquity of Saturn for different migration histories of Titan.** On the horizontal axis, $t_{\text{tide}} = a/\dot{a}$ is Titan's current migration timescale. The horizontal axis is not linear but scaled as $1/t_{\text{tide}}$. Between the white curves, Saturn's past obliquity is smaller than 5°. The curve of Fig. 3 corresponds here to a vertical slice at $t_{\text{tide}} \approx 13$ Gyr. With the assumption that $t_{\text{tide}} = 6.5$ Gyr instead, the resonance with $s_8$ would have been encountered about 0.4 Gyr ago, whereas with $t_{\text{tide}} = 65$ Gyr, it would have been encountered about 2.8 Gyr ago.</figcaption></figure>

law of ref. 7 corresponds today to $t_{\text{tide}} \approx 13$ Gyr, and radio-science experiments yield values ranging between 9.3 Gyr and 13.2 Gyr at the $3\sigma$ uncertainty level. As Titan's fast migration may have been episodic¹⁴, we extend our analysis beyond this interval. Our result summary (Fig. 4) shows a very horizontal-like structure in spite

of shallow wavy patterns, which indicates that Titan's migration is an almost adiabatic process. Therefore, Saturn's spin-axis dynamics depend mostly on the range over which Titan migrates and not much on its migration rate. For instance, the results obtained for $t_{\text{tide}} = 65$ Gyr are similar to those obtained in a scenario in which Titan remained unmoved for billions of years and then started to migrate at its observed rate during the last 1.5 Gyr only. This example shows that our results are valid for a large variety of migration histories that are compatible with the measurements reported in ref. 7, and even if Titan's migration was sporadic. The minimum past obliquity of Saturn is firmly below 5° in an almost unvaried range of $α_0$ values. The remaining slightly non-zero excess of obliquity can be explained by many different phenomena, such as the light bombardment phase that occurred at the end of Saturn's formation, or an abrupt resonance crossing that stemmed from the dissipation of the protoplanetary disc¹⁹. The tidal migration of Titan is therefore able to explain the current large obliquity of Saturn starting from a small, presumably primordial, value.

This scenario is, however, not unique. On the one hand, our current knowledge of the tidal resonance-locking mechanism of ref. 14 leaves the possibility that we are now observing a very fast migration of Titan that otherwise remained nearly zero over all of Titan's history. This extreme possibility appears unlikely, as the resonance-locking mechanism has seemingly triggered for all of the six moons of Saturn studied in ref. 7, which favours a continuous rather than sporadic migration. On the other hand, the possibility remains that Saturn acquired its large obliquity via an early impact²⁰ and was then not much affected by the recent resonance crossing with $s_8$ (Fig. 2b,c). Ruling out this possibility would require a better knowledge of Titan's polar moment of inertia, preferably obtained in a model-independent way.

The mechanism presented in this letter is generic and also applies to Jupiter, even though Jupiter did not have enough time to be tilted substantially yet¹¹. Giant exoplanets could harbour massive satellites as well²¹. If a fast satellite migration is indeed the rule⁷, the obliquities of such planets are vulnerable to dramatic increase via the same mechanism. The existence of secular spin-orbit resonances requires only a precessing orbit, for instance due to other planetary companions¹⁵. The number of strong harmonics in the precession spectrum of a planet is larger for planetary systems that feature substantial eccentricities and inclinations. The various orbital architectures observed in exoplanetary systems should therefore allow for much more diverse secular spin-orbit resonances than for the giant planets of the Solar System, which have today rather cold orbits²². The mechanism presented in this letter could therefore greatly affect the obliquity distribution of exoplanets.

## Methods

**Long-term dynamics of Saturn's spin axis.** The long-term spin-axis dynamics of an oblate planet is ruled by the Hamiltonian function

$$ \mathcal{H}(X, -\psi, t) = -\frac{\alpha}{2} \frac{X^2}{(1-e(t)^2)^{3/2}} - \sqrt{1-X^2} (\mathcal{A}(t) \sin \psi + \mathcal{B}(t) \cos \psi) + 2XC(t), \quad (1) $$

where the conjugate coordinates are $X = \cos \epsilon$ (cosine of obliquity) and $-\psi$ (minus precession angle)¹². The quantity $\alpha$ is called the precession constant. Equation (1) explicitly depends on time $t$ through the orbital eccentricity $e$ of the planet and through the functions

$$ \left\{ \begin{array}{l} \mathcal{A}(t) = \frac{2(q+p C(t))}{\sqrt{1-p^2-q^2}}, \\ \mathcal{B}(t) = \frac{2(p-q C(t))}{\sqrt{1-p^2-q^2}}, \end{array} \right. \quad \text{and} \quad C(t) = \dot{q}\dot{p} - p\dot{q}. \quad (2) $$

In these expressions, $q = \eta \cos \Omega$ and $p = \eta \sin \Omega$, where $\eta = \sin(I/2)$ and $I$ and $\Omega$ are the orbital inclination and the longitude of ascending node of the planet, respectively. If the planet has migrating satellites, its precession constant $\alpha$ is a slowly varying function of time.
**Saturn's precession constant.** The precession constant of a planet can be written

$$ \alpha = \frac{3 G m_{\odot} J_2 + q}{2 \omega a_s^2} \frac{\lambda + \ell}{\lambda + \ell}, \quad (3) $$

where $q$ and $\ell$ are the satellites' contributions$^{1,3}$. Here, $Gm_{\odot}$ is the Sun's gravitational parameter, $\omega$ is Saturn's spin rate, $a_s$ its semi-major axis, $J_2$ its second zonal gravity coefficient and $\lambda$ its polar moment of inertia C normalized by using its mass $M$ and its equatorial radius $R_{eq}$ as $\lambda \equiv C/(MR_{eq}^2)$.

When computing $\alpha$ through equation (3), the values of $J_2$, $\lambda$, $q$ and $\ell$ must be normalized by using the same reference radius so that the final value of $\alpha$ is independent of the normalization used. Different normalization conventions are the source of much confusion in the literature. Here, we use quantities normalized by using $R_{eq}=60,268$ km; values taken from other articles are renormalized accordingly. The choice of $a_s$ in equation (3) is also critical, as it can lead to substantially different values of $\alpha$ and place Saturn today inside or outside the resonance with $s_6$. As equation (3) relates to secular theories (that is, averaged over rotational and orbital motions), one must use secular elements such as those computed in the context of ref. $^{17}$. The spin rate and current spin-axis orientation of Saturn are given by the latest IAU (International Astronomical Union) recommendations$^{2,3}$ and its $J_2$ is given in ref. $^{24}$.

Saturn's normalized polar moment of inertia $\lambda$ is still poorly constrained today. Even though many different values can be found in the literature, they are all model-dependent and generally disagree with each other. Citing ref. $^{25}$, authors often consider $\lambda=0.22$ as a central value, with a 10% error meant to account for model dependence$^{1,8,26}$. This corresponds to an interval that is slightly larger than $\lambda \in [0.200, 0.240]$. In more recent articles, the actual uncertainties of estimates of $\lambda$ still seem unclear, even when taking into account the data from the Cassini mission. When reviewing values obtained through different models, specialists still adhere to $\lambda=0.22$ with a 10% variation$^{27}$.

Values of $q$ and $\ell$ including the meaningful contributions from all satellites of Saturn have been computed by other authors$^{26}$. Here, we use updated values of $q$ and $\ell$ obtained by taking into account the non-zero equilibrium orbital inclination ('Laplace plane') for all satellites, and by using the satellites' masses and mean orbital elements that come from modern ephemerides (Jet Propulsion Laboratory SAT427 (https://ssd.jpl.nasa.gov/)). Compared to the values in ref. $^{26}$, these refined values of $q$ and $\ell$ lead to a slightly different relation between the value of Saturn's current precession constant $a_0$ and the value of Saturn's normalized polar moment of inertia $\lambda$.

**Migration of Saturn's satellites.** The contributions $q$ and $\ell$ of migrating satellites vary over time, producing a drift of $\alpha$. Through the tidal theory of ref. $^{14}$, the semi-major-axis variations of most major satellites of Saturn are approximated in ref. $^{7}$ as

$$ a(t) \approx a_0 \left(\frac{t}{t_0}\right)^B, \quad (4) $$

where $a_0$ is the satellite's current mean semi-major axis, $t_0$ is Saturn's current age, and $B \approx 1/3$ is a constant. The tidal timescale $t_{\text{tide}} = a/\dot{a}$ of the satellites can be derived from equation (4); it evolves in the long run, with a current value of $t_{\text{tide}}=t_0/B \approx 13$ Gyr.

Being so much more massive, Titan is the satellite that contributes by far the most to the value of $\alpha$. The second largest contribution comes from Iapetus, but it is only 5% of Titan's$^{28}$. Therefore, most of the drift of $\alpha$ over time is due to Titan's migration. For this reason and because of the quasi-adiabatic nature of the drift, other satellites can be modelled as a slight increase in Titan's mass, which allows us to include only Titan in equation (3). The adjusted mass of Titan is about 4% larger. A more refined model that takes into account the migration of all satellites (or variations of Saturn's physical parameters in equation (3)) would only change the evolution timescale of $\alpha$ slightly (Fig. 1), and not our conclusions. Here, the slight increase in Titan's mass is meant only to provide the right connection between $\lambda$ and today's value of $\alpha$, such that $\alpha(t=0)=a_0$.

The effect of a past massive ring on the value of $\alpha$ can be estimated by putting an annulus of material at Saturn's fluid Roche radius in equation (3). Giving it five times the mass of Titan, as suggested in ref. $^{28}$, does not substantially modify the value of $\alpha$. Our conclusions about the past small value of $\alpha$ hold even if the annulus is several tens of times more massive than Titan. Only a wide and massive long-lived circumplanetary disc could possibly have enhanced $\alpha$ enough to reach a resonance$^{19}$, but its decay would then have produced only a small obliquity kick$^{6}$.

**Orbital solution of Saturn.** The Hamiltonian function in equation (1) depends on the orbit of the planet and on its temporal variations. We need a secular orbital solution that is valid over billions of years, which is well beyond the timespan covered by ephemerides. Luckily, the orbital dynamics of the giant planets of the Solar System are almost integrable and excellent solutions have been developed. We use the secular solution of ref. $^{17}$ expanded into quasi-periodic series:

$$ \begin{aligned} z &= e \exp(i\dot{\omega}) = \sum_k E_k \exp(i\theta_k), \\ \zeta &= \eta \exp(i\Omega) = \sum_k S_k \exp(i\phi_k), \end{aligned} \quad (5) $$

where $\dot{\omega}$ is Saturn's longitude of perihelion, and $k$ runs over all terms of the series. The amplitudes $E_k$ and $S_k$ are real constants, and the angles $\theta_k$ and $\phi_k$ evolve linearly over time $t$, with frequencies $\mu_k$ and $\nu_k$:

$$ \theta_k(t) = \mu_k t + \theta_k^{(0)}, \quad \text{and} \quad \phi_k(t) = \nu_k t + \phi_k^{(0)}. \quad (6) $$

As described in ref. $^{11}$, the complete secular orbital solution of ref. $^{17}$ is obtained by multiplying the normalized proper modes $z_i^*$ and $\zeta_i^*$ (tables VI and VII of ref. $^{17}$) by the matrix $\mathcal{S}$ that corresponds to the linear part of the solution (table V of ref. $^{17}$). In the series obtained, the terms with the same combination of frequencies are then merged, which results in 56 terms in eccentricity and 60 terms in inclination, with amplitudes down to $10^{-8}$. These terms contain contributions from all the planets of the Solar System. Even though frequencies associated to the terrestrial planets could vary somewhat over billions of years as a result of their chaotic dynamics$^{17}$, they contribute only marginally to Saturn's orbital solution and none of them takes part in the resonances shown in Fig. 2. This secular orbital solution can therefore be considered valid over a billion-year timescale, since the end of the late planetary migration.

Because the orbital series of Saturn are not restricted here to the $s_8$ harmonic, our results differ to some degree from previous analytical predictions. In particular, the minimum past obliquity of Saturn is about $2.5^\circ$ for an adiabatic resonance capture (see main text). This value is smaller than the $4.5^\circ$ predicted previously$^{1}$, because of extra obliquity oscillations produced by other terms in the orbital series.

**Numerical integration.** The dynamical system described by equation (1) is propagated numerically. The precession constant $\alpha$ is made to vary according to the migration of Titan in equation (4), and the orbit of Saturn is made to vary according to the quasi-periodic series in equation (5). Instead of the variables $X$ and $\psi$ shown in equation (1), numerical calculations are made with the equivalent spherical coordinates of the spin axis

$$ \begin{cases} x = \sin(\varepsilon) \cos(\psi), \\ y = \sin(\varepsilon) \sin(\psi), \\ z = \cos(\varepsilon), \end{cases} \quad (7) $$

which are not singular for a zero obliquity. Numerical integrations are performed with a step size small enough such that $\sqrt{x^2 + y^2 + z^2}$ is conserved to within the computer's rounding-off errors.

## Data availability

Source Data for the figures are provided with this paper. The data supporting other findings of this study are available from the corresponding author upon reasonable request.

## Code availability

The parameters and equations of motion are fully described within the paper. All data can be reproduced using any standard implementation. The numerical integration scheme used is fully available from ref. $^{29}$.

Received: 5 June 2020; Accepted: 30 November 2020;

Published online: 18 January 2021

## References

1. Ward, W. R. & Hamilton, D. P. Tilting Saturn. I. Analytic model. *Astron. J.* **128**, 2501–2509 (2004).

2. Rogoszinski, Z. & Hamilton, D. P. Tilting ice giants with a spin-orbit resonance. *Astrophys. J.* **888**, 60 (2020).

3. Malhotra, R. The origin of Pluto's peculiar orbit. *Nature* **365**, 819–821 (1993).

4. Tsiganis, K., Gomes, R., Morbidelli, A. & Levison, H. F. Origin of the orbital architecture of the giant planets of the Solar System. *Nature* **435**, 459–461 (2005).

5. Nesvorný, D. & Morbidelli, A. Statistical study of the early Solar System's instability with four, five, and six giant planets. *Astron. J.* **144**, 117 (2012).

6. Hamilton, D. P. & Ward, W. R. Tilting Saturn. II. Numerical model. *Astron. J.* **128**, 2510–2517 (2004).

7. Lainey, V. et al. New tidal paradigm in giant planets supported by rapid orbital expansion of Titan. *Nat. Astron.* **4**, 1053–1058 (2020).

8. Boué, G., Laskar, J. & Kuchynka, P. Speed limit on Neptune migration imposed by Saturn tilting. *Astrophys. J.* **702**, L19–L22 (2009).

9. Brasser, R. & Lee, M. H. Tilting Saturn without tilting Jupiter: constraints on giant planet migration. *Astron. J.* **150**, 157 (2015).

10. Vokrouhlický, D. & Nesvorný, D. Tilting Jupiter (a bit) and Saturn (a lot) during planetary migration. *Astrophys. J.* **806**, 143 (2015).

11. Saillenfest, M., Lari, G. & Courtot, A. The future large obliquity of Jupiter. *Astron. Astrophys.* **640**, A11 (2020).

12. Laskar, J. & Robutel, P. The chaotic obliquity of the planets. *Nature* **361**, 608–612 (1993).
13. Ward, W. R. Tidal friction and generalized Cassini's laws in the solar system. *Astron. J.* **80**, 64–70 (1975).

14. Fuller, J., Luan, J. & Quataert, E. Resonance locking as the source of rapid tidal migration in the Jupiter and Saturn moon systems. *Mon. Not. R. Astron. Soc.* **458**, 3867–3879 (2016).

15. Millholland, S. & Laughlin, G. Obliquity-driven sculpting of exoplanetary systems. *Nat. Astron.* **3**, 424–433 (2019).

16. Saillenfest, M., Laskar, J. & Boué, G. Secular spin-axis dynamics of exoplanets. *Astron. Astrophys.* **623**, A4 (2019).

17. Laskar, J. The chaotic motion of the solar system: a numerical estimate of the size of the chaotic zones. *Icarus* **88**, 266–291 (1990).

18. Charnoz, S., Canup, R. M., Crida, A. & Dones, L. in *Planetary Ring Systems: Properties, Structure, and Evolution* (eds Tiscareno, M. S. & Murray, C. D.) 517–538 (Cambridge Univ. Press, 2018).

19. Millholland, S. & Batygin, K. Excitation of planetary obliquities through planet–disk interactions. *Astrophys. J.* **876**, 119 (2019).

20. Safronov, V. S. Sizes of the largest bodies falling onto the planets during their formation. *Sov. Astron.* **9**, 987–991 (1966).

21. Kipping, D. M. In search of exomoons. In *Proc. Frank N. Bash Symposium 2013: New Horizons in Astronomy*, 012 (PoS, 2014).

22. Kreyche, S. M. et al. Retrograde-rotating exoplanets experience obliquity excitations in an eccentricity-enabled resonance. *Planet. Sci. J.* **1**, 8 (2020).

23. Archinal, B. A. et al. Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015. *Celest. Mech. Dyn. Astron.* **130**, 22 (2018).

24. Iess, L. et al. Measurement and implications of Saturn's gravity field and ring mass. *Science* **364**, aat2965 (2019).

25. Hubbard, W. B. & Marley, M. S. Optimized Jupiter, Saturn, and Uranus interior models. *Icarus* **78**, 102–118 (1989).

26. French, R. G. et al. Geometry of the Saturn system from the 3 July 1989 occultation of 28 Sgr and Voyager observations. *Icarus* **103**, 163–214 (1993).

27. Fortney, J. J. et al. in *Saturn in the 21st Century* (eds Baines, K. H. et al.) 44–68 (Cambridge Univ. Press, 2018).

28. Hesselbrock, A. J. & Minton, D. A. An ongoing satellite–ring cycle of Mars and the origins of Phobos and Deimos. *Nat. Geosci.* **10**, 266–269 (2017).

29. Rein, H. & Spiegel, D. S. IAS15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits. *Mon. Not. R. Astron. Soc.* **446**, 1424–1437 (2015).

## Acknowledgements

We thank L. Gomez Casajus for fruitful discussions and S. Renner for his suggestions during the writing of our manuscript. G.L. acknowledges financial support from the Italian Space Agency (ASI) through agreement no. 2017-40-H.O.

## Author contributions

G.B. conceived the original idea. G.L. compiled the data. G.L. and M.S. made the computations. M.S. wrote the article. All authors participated in supervising the whole study.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-020-01284-x.

**Correspondence and requests for materials** should be addressed to M.S.

**Peer review information** *Nature Astronomy* thanks David Nesvorny and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2021