© M. Saillenfest et al. 2021

1 Introduction

The obliquity of a planet is the angle between its spin axis and the normal to its orbit. In the protoplanetary disc, giant planets are expected to form with near-zero obliquities (Ward & Hamilton 2004; Rogoszinski & Hamilton 2020). After the formation of Saturn, some dynamical mechanism must therefore have tilted its spin axis up to its current obliquity of 26.7°.

Ward & Hamilton (2004) showed that Saturn is currently located very close to a secular spin-orbit resonance with the nodal precession mode of Neptune. This resonance strongly affects Saturn’s spin axis today, and it offers a tempting explanation for its current large obliquity. For years, the scenarios that were most successful in reproducing Saturn’s current obliquity through this resonance invoked the late planetary migration (Hamilton & Ward 2004; Boué et al. 2009; Vokrouhlický & Nesvorný 2015; Brasser & Lee 2015). However, Saillenfest et al. (2021) have recently shown that this picture is incompatible with the fast tidal migration of Titan detected by Lainey et al. (2020) in two independent sets of observations – assuming that this migration is not specific to the present epoch but went on over a substantial interval of time. Indeed, satellites affect the spin-axis precession of their host planets (see e.g. Ward 1975; Tremaine 1991; Laskar et al. 1993; Boué & Laskar 2006). Since the effect of a satellite depends on its orbital distance, migrating satellites induce a long-term drift in the planet’s spin-axis precession velocity. In the course of this drift, large obliquity variations can occur if a secular spin-orbit resonance is encountered (i.e. if the planet’s spin-axis precession velocity becomes commensurate with a harmonic of its orbital precession). Because of this mechanism, dramatic variations in the Earth’s obliquity are expected to take place in a few billion years from now, as a result of the Moon’s migration (Néron de Surgy & Laskar 1997). Likewise, Jupiter’s obliquity is likely steadily increasing today and could exceed 30° in the next billions of years, as a result of the migration of the Galilean satellites (Lari et al. 2020; Saillenfest et al. 2020).

A significant migration of Saturn’s satellites implies that, contrary to previous assumptions, Saturn’s spin-axis precession velocity was much smaller in the past, precluding any resonance with an orbital frequency. The same conclusion could also hold for Jupiter (Lainey et al. 2009; Lari et al. 2020). In fact, Saillenfest et al. (2021) have shown that Titan’s migration itself is likely responsible for the resonant encounter between Saturn’s spin axis and the nodal precession mode of Neptune. Their results indicate that this relatively recent resonant encounter could explain the current large obliquity of Saturn starting from a small value, possibly less than 3°. This new paradigm solves the problem of the low probability of reproducing both the orbits and axis tilts of Jupiter and Saturn during the late planetary migration (Brasser & Lee 2015). However, it revokes the concomitant constraints on the parameters of the late planetary migration.

The findings of Saillenfest et al. (2021) have been obtained through backward integrations from Saturn’s current spin orientation, and by exploring migration histories for Titan in the vicinity of the nominal scenario of Lainey et al. (2020). However, observation uncertainties and our lack of knowledge about the past evolution of Titan’s migration rate still allow for a large variety of migration histories, and one can wonder whether the dramatic influence of Titan is a generic result or whether it is restricted to the range of parameters explored by Saillenfest et al. (2021). Moreover, even though backward numerical integrations do prove that Titan’s migration is able to raise Saturn’s obliquity, a statistical picture of the possible trajectories that could have been followed is still missing. In this regard, the likelihood offollowing a given dynamical pathway would be quite valuable, because it could be used as a constraint to the parameters of the model, in the spirit of Boué et al. (2009), Brasser & Lee (2015), and Vokrouhlický & Nesvorný (2015).

For these reasons, we aim to explore the outcomes given by all conceivable migration timescales for Titan, and to perform a statistical search for Saturn’s past obliquity. This will provide the whole region of the parameter space allowing Titan’s migration to be responsible for Saturn’s large obliquity, with the corresponding probability. Finally, since Titan still migrates today, Saturn’s obliquity could suffer from further large variations in the future, in the same way as Jupiter (Saillenfest et al. 2020). Therefore, we also aim to extend previous analyses to the future dynamics of Saturn’s spin axis.

Our article is organised as follows. In Sect. 2, we recall the dynamical model used by Saillenfest et al. (2021) and discuss the range of acceptable values for the physical parameters of Saturn and its satellites. Sections 3 and 4 are dedicated to the past spin-axis dynamics of Saturn: after having explored the parameter space and quantified the importance of non-adiabaticity, we perform Monte Carlo experimentsto search for the initial conditions of Saturn’s spin axis. In Sect. 5, we present our results about the obliquity values that will be reached by Saturn in the future. Finally, our conclusions are summarised in Sect. 6.

2 Secular dynamics of the spin axis

2.1 Equations of motion

In the approximation of rigid rotation, the spin-axis dynamics of an oblate planet subject to the lowest-order term of the torque from the Sun is given for instance by Laskar & Robutel (1993) or Néron de Surgy & Laskar (1997). Far from spin-orbit resonances, and due to the weakness of the torque, the long-term evolution of the spin axis is accurately described by the secular Hamiltonian function (i.e. averaged over rotational and orbital motions). This Hamiltonian can be written (1)

where the conjugate coordinates are X = cosε (cosine of obliquity) and − ψ (minus the precession angle). The Hamiltonian in Eq. (1) explicitly depends on time t through the orbital eccentricity e of the planet and through the functions (2)

In these expressions, q = ηcos Ω and p = ηsinΩ, where η ≡ sin(I∕2), and I and Ω are the orbital inclination and the longitude of ascending node of the planet, respectively. The quantity α is called the precession constant. It depends on the spin rate of the planet and on its mass distribution, through the formula: (3)

where is the gravitational constant, m_⊙ is the mass of the sun, ω is the spin rate of the planet, a is its semi-major axis, J₂ is its second zonal gravity coefficient, and λ is its normalised polar moment of inertia. The parameters J₂ and λ can be expressed as (4)

where A, B, and C are the equatorial and polar moments of inertia of the planet, M is its mass, and R_eq is its equatorial radius.

The precession rate of the planet is increased if it possesses massive satellites. Far-away satellites increase the torque exerted by the sun on the equatorial bulge of the planet, whereas close-in satellites artificially increase the oblateness and the rotational angular momentum of the planet (Boué & Laskar 2006). In the close-in regime, an expression for the effective precession constant has been derived by Ward (1975). It has been generalised by French et al. (1993) who included the effect of the non-zero orbital inclinations of the satellites, as they oscillate around their local ‘Laplace plane’ (see e.g. Tremaine et al. 2009). The effective precession constant is obtained by replacing J₂ and λ in Eq. (3) by the effective values: (5)

where m_k, a_k and n_k are the mass, the semi-major axis, and the mean motion of the kth satellite, ε is the obliquity of the planet, and L_k is the inclination of the Laplace plane of the kth satellite with respect to the planet’s equator. For regular satellites, L_k lies between 0 (close-in satellite) and ε (far-away satellite). The formulas of French et al. (1993) given by Eq. (5) are valid whatever the distance of the satellites, and they closely match the general precession solutions of Boué & Laskar (2006). We can also verify that the small eccentricities of Saturn’s major satellites do not contribute substantially to and λ′, allowing us to neglect them.

Because of its large mass, Titan is by far the satellite that contributes most to the value of α. Therefore, even thoughits Laplace plane is not much inclined, taking its inclination into account changes the whole satellites’ contribution by several percent¹. Tremaine et al. (2009) give a closed-form expression for L_k in the regime m_k ≪ M, where all other satellites j with a_j < a_k are also taken into account. The values obtained for Titan (L₆ ≈ 0.62°) and Iapetus (L₈ ≈ 16.03°) are very close to those found in the quasi-periodic decomposition of their ephemerides (see e.g. Vienne & Duriez 1995). The inclinations L_k of the other satellites of Saturn do not contribute substantially to the value of α.

Even though the value of α computed using Eq. (5) yields an accurate value of the current mean spin-axis precession velocity of Saturn as , it cannot be directly used to propagate the dynamics using the Hamiltonian function in Eq. (1), because α would itself be a function of X, which contradicts the Hamiltonian formulation. For this reason, authors generally assume that α only weakly depends on ε, such that the satellite’s contributions can be considered to be fixed while ε varies according toHamilton’s equations of motion (see e.g. Ward & Hamilton 2004; Boué et al. 2009; Vokrouhlický & Nesvorný 2015; Brasser & Lee 2015). In our case, Titan largely dominates the satellite’s contribution, and it is almost in the close-in regime (L₆ ≪ ε). We can therefore use the same trick as Saillenfest et al. (2021) and replace Eq. (5) by (6)

where only Titan is considered (k = 6), in the close-in regime (L₆ = 0), and where its mass m₆ has been slightly increased () so as to produce the exact same value of α today using Eq. (6) instead of Eq. (5). This slight increase in Titan’s mass has no physical meaning; it is only used here to provide the right connection between λ and today’s value of α. This point is further discussed in Sect. 2.3.

2.2 Orbital solution

The Hamiltonian function in Eq. (1) depends on the orbit of the planet and on its temporal variations. In order to explore the long-term dynamics of Saturn’s spin axis, we need an orbital solution that is valid over billions of years. In the same way as Saillenfest et al. (2020), we use the secular solution of Laskar (1990) expanded in quasi-periodic series: (7)

where ϖ is Saturn’s longitude of perihelion. The amplitudes E_k and S_k are real constants, and the angles θ_k and ϕ_k evolve linearly overtime t with frequencies μ_k and ν_k: (8)

See Appendix A for the complete orbital solution of Laskar (1990).

The series in Eq. (7) contain contributions from all the planets of the Solar System. In the integrable approximation, the frequency of each term corresponds to a unique combination of the fundamental frequencies of the system, usually noted g_j and s_j. In the limit of small masses, small eccentricities and small inclinations (Lagrange-Laplace secular system), the z series only contains the frequencies g_j, while the ζ series only contains the frequencies s_j (see e.g. Murray & Dermott 1999 or Laskar et al. 2012). This is not the case in more realistic situations. Table 1 shows the combinations of fundamental frequencies identified for the largest terms of Saturn’s ζ series obtained by Laskar (1990).

As explained by Saillenfest et al. (2019), at first order in the amplitudes S_k and E_k, secular spin-orbit resonant angles can only be of the form σ_p = ψ + ϕ_p, where p is a given index in the ζ series. Resonances featuring terms of the z series only appear at third order and beyond. For the giant planets of the Solar System, the existing secular spin-orbit resonances are small and isolated from each other, and only first-order resonances play a substantial role (see e.g. Saillenfest et al. 2020).

Figure 1 shows the location and width of every first-order resonance for the spin-axis of Saturn in an interval of precession constant α ranging from 0″ yr⁻¹ to 2″ yr⁻¹. Because of the chaotic dynamics of the Solar System (Laskar 1989), the fundamental frequencies related to the terrestrial planets (e.g. s₁, s₂, s₄, and γ appearing in Table 1) could vary substantially over billions of years (Laskar 1990). However, they only marginally contribute to Saturn’s orbital solution and none of them takes part in the resonances shown in Fig. 1. Our secular orbital solution for Saturn can therefore be considered valid since the late planetary migration, which presumably ended at least 4 Gyr ago (see e.g. Nesvorný & Morbidelli 2012; Deienno et al. 2017; Clement et al. 2018). For this reason, we consider in all this article a maximum timespan of 4 Gyr in the past. As shown by Saillenfest et al. (2021), this timespan is more than enough for Saturn to relax to its primordial obliquity value. Our results are therefore independent of this choice, unless one considers a much slower migration rate for Titan than observed today. This last case is discussed in Sect. 3.2.

Fig. 1 Location and width of every first-order secular spin-orbit resonance for Saturn. Each resonant angle is of the form σp = ψ + ϕp where ϕp has frequency νp labelled on the graph according to its index in the orbital series (see Table 1 and Appendix A). For a given value of the precession constant α, the intervalof obliquity enclosed by the separatrix is shown in pink, as computed using the formulas of Saillenfest et al. (2019). The green bar shows Saturn’s current obliquity and the range for its precession constant considered in this article, as detailed in Sects. 2.3 and 2.4.

2.3 Precession constant

As shown by the Hamiltonian function in Eq. (1), the precession constant α is a key parameter of the spin-axis dynamics of a planet. The physical parameters of Saturn that enter into its expression (see Eq. (3)) are all very well constrained from observations, except the normalised polar moment of inertia λ. Indeed, the gravitational potential measured by spacecrafts only provides differences between the moments of inertia (e.g. the coefficient J₂). In order to obtain the individual value of a single moment of inertia, one would need to detect the precession of the spin axis or the Lense–Thirring effect, as explained for instance by Helled et al. (2011). Such measurements are difficult considering the limited timespan covered by space missions. To our knowledge, the most accurate estimate of Saturn’s polar motion, including decades of astrometric observations and Cassini data, is given by French et al. (2017). However, their estimate is still not accurate enough to bring anydecisive constraint on Saturn’s polar moment of inertia. Moreover, since the observed polar motion of Saturn is affected by many short-period harmonics, it cannot be directly linked to the secular spin-axis precession rate discussed in this article. Removing short-period harmonics from the observed signal would require an extensive modelling that is not yet available. Even though some attempts to compute a secular trend from Saturn’s spin-axis observations have been reported (as the unpublished results of Jacobson cited by Vokrouhlický & Nesvorný 2015), we must still rely on theoretical values of λ.

As pointed out by Saillenfest et al. (2020), one must be careful about the normalisation used for λ. Here, we adopt R_eq = 60268 km by convention and we renormalise each quantity in Eqs. (4) and (5) accordingly. Many different values of λ can be found in the literature. Under basic assumptions, Jeffreys (1924) obtained a value of 0.198. This value is smaller thanother estimates found in the literature, even though it is marginally compatible with the calculations of Hubbard & Marley (1989), who gave λ = 0.22037 with a 10% uncertainty. The latter value and its uncertainty have been reused by many authors afterwards, including French et al. (1993) and Ward & Hamilton (2004). Later on, Helled et al. (2009) obtained values of λ ranging between 0.207 and 0.210. From an exploration of the parameter space, Helled (2011) then found λ ∈ [0.200, 0.205], but the normalisation used in this article is ambiguous². The computations of Nettelmann et al. (2013) yielded yet another range for λ, estimated to lie in [0.219, 0.220]. Among the alternative models proposed by Vazan et al. (2016), values of λ are found to range between 0.222 and 0.228. Finally, Movshovitz et al. (2020) used a new fitting technique supposed to be less model-dependent, and obtained λ ∈ [0.2204, 0.2234] at the 3σ error level (assuming that their values are normalised using R_eq, which is not specified in the article). In the review of Fortney et al. (2018) focussing on the better knowledge of Saturn’s internal structure brought by the Cassini mission, the authors go back to a value of λ equal to 0.22 ± 10%. A value of 0.22 is also quoted in the review of Helled (2018).

Here, instead of relying on one particular estimate of λ, we turn to the exploration of the whole range of values given in the literature, which is slightly larger than λ ∈ [0.200, 0.240]. The spin velocity of Saturn is taken from Archinal et al. (2018) and its J₂ from Iess et al. (2019). For consistency with Saturn’s orbital solution (Sect. 2.2), we take its mass and secular semi-major axis from Laskar (1990).

In order to compute and λ′ in Eq. (5), we need the masses and orbital elements of Saturn’s satellites. We take into account the eight major satellites of Saturn and use the masses of the SAT427 numerical ephemerides³. These ephemerides are then digitally filtered in order to obtain the secular semi-major axes. The inclination L_k of the Laplace plane of each satellite is computed using the formula of Tremaine et al. (2009). Taking λ into its exploration interval, the current value of Saturn’s precession constant, computed from Eqs. (3) and (5), ranges from 0.747″ yr⁻¹ to 0.894″ yr⁻¹. The corresponding adjusted mass of Titan in Eq. (6) is . Similar results are obtained when using the more precise values of L_k given by theconstant terms of the full series of Vienne & Duriez (1995) and Duriez & Vienne (1997).

Becauseof tidal dissipation, satellites migrate over time. This produces a drift of the precession constant α on a timescale that is much larger than the precession motion (i.e. the circulation of ψ). The long-term spin-axis dynamics of a planet with migrating satellites is described by the Hamiltonian in Eq. (1) where α is a slowly-varying function of time. Since Titan is in the close-in regime, its outward migration produces an increase in α. The migration rate of Titan recently measured by Lainey et al. (2020) supports the tidal theory of Fuller et al. (2016), through which the time evolution of Titan’s semi-major axis can be expressed as (9)

where a₀ is Titan’s current mean semi-major axis, t₀ is Saturn’s current age, and b is a real parameter (see Lainey et al. 2020). Even though Eq. (9) only provides a crude model for Titan’s migration, the parameter b can be directly linked to the observed migration rate, independently of whether Eq. (9) is valid or not⁴. Equation (9) implies that Titan’s current tidal timescale t_tide = a₆∕ȧ₆ relates to b as b = t₀∕t_tide. Considering a 3σ error interval, the astrometric measurements of Lainey et al. (2020) yield values of b ranging in [0.18, 1.71], while their radio-science experiments yield values ranging in [0.34, 0.49]. For the long-term evolution of Saturn’s satellites, they adopt a nominal value of b₀ = 1∕3, which roughly matches the observed migration of all satellites studied. Using this nominal value, we obtain a drift of the precession constant α as depicted in Fig. 2. Taking b as parameter, a migration n times faster for Titan is obtained by using in Eq. (9) a parameter b = n b₀. The corresponding evolution of Titan’s semi-major axis is illustrated in Fig. 3.

As mentioned by Saillenfest et al. (2020), other parameters in Eq. (3) probably slightly vary over billions of years, such as the spin velocity of Saturn or its oblateness. We consider that the impact of their variations is small compared to the effect of Titan’s migration (see Fig. 2) and contained within our exploration range. Moreover, all satellites, and not only Titan, migrate over time. However, being Titan so much more massive, its fast migration is by far the dominant cause of the drift of α. Since its exact migration rate is still uncertain (see Fig. 3), this justifies our choice to only include Titan in Eq. (6), while the use of its slightly increased mass yet allows us to obtain the right value of today’s precession constant α, as if all satellites were included.

Fig. 2 Evolution of the effective precession constant of Saturn due to the migration of Titan (adapted from Saillenfest et al. 2021). The top and bottom green curves correspond to the two extreme values of the normalised polar moment of inertia λ considered in this article. They appear into α through Eq. (3). Both curves are obtained using the nominal value b = 1∕3 in Eq. (9). Today’s interval corresponds to the one shown in Fig. 1; it is independent of the value of b considered. The blue line shows Neptune’s nodal precession mode, which was higher before the end of the late planetary migration.

2.4 Current spin orientation

The initial orientation of Saturn’s spin axis is taken from the solution of Archinal et al. (2018) averaged over short-period terms. Withrespect to Saturn’s secular orbital solution (see Sect. 2.2), this gives an obliquity ε = 26.727° and a precession angle ψ = 6.402° at time J2000. The uncertainty on these values is extremely small compared to the range of α considered (see Sect. 2.3).

Fig. 3 Time evolution of Titan’s semi-major axis for different migration rates. The pink and blue intervals show the 3σ uncertainty ranges of astrometric and radio-science measurements, respectively (Lainey et al. 2020). The coloured curves are obtained by varying the parameter b in Eq. (9).

3 The past obliquity of Saturn: exploration of the parameter space

3.1 Overview of possible trajectories

From the results of their backward numerical integrations, Saillenfest et al. (2021) find that Saturn can have evolved through distinct kinds of evolution, which had previously been described by Ward & Hamilton (2004). These different kinds of evolution are set by the outcomes of the resonant encounter between Saturn’s spin-axis precession and the nodal precession mode of Neptune (term ϕ₃ in Table 1 and largest resonance in Fig. 1). The four possible types of past evolution are illustrated in Fig. 4 for b = b₀. They are namely:

• Type 1: for λ ≤ 0.220, Saturn went past the resonance through its hyperbolic point.

• Type 2: for λ ∈ (0.220, 0.224) ∪ (0.237, 0.241), Saturn was captured recently by crossing the separatrix of the resonance and followed the drift of its centre afterwards.

• Type 3: for λ ∈ [0.224, 237], the separatrix of the resonance appeared around Saturn’s trajectory resulting in a 100%-sure capture at low obliquity. Saturn followed the drift of its centre afterwards.

• Type 4: for λ ≥ 0.241, Saturn did not reach yet the resonance.

Figure 5 shows the current oscillation interval of Saturn’s spin axis in all four cases. Trajectories of Type 3 are those featuring the smallest libration amplitude of the resonant angle σ₃ and allowing for the smallest past obliquity of Saturn. Type 4 is ruled out by our uncertainty range for λ.

During its past evolution, Saturn also crossed a first-order secular spin-orbit resonance with the term ϕ₁₉ which has frequency g₇ − g₈ + s₇ (see Table 1). As shown in Fig. 4, however, this did not produce any noticeable change in obliquity for Saturn. Indeed, since this resonance is very small, the oscillation timescale of σ₁₉ = ψ + ϕ₁₉ inside the resonance is dramatically longer than the duration of the resonance crossing. This results in a short non-adiabatic crossing. The difference of oscillation timescales of σ₃ and σ₁₉ can be appreciated in Fig. 6. It explains why these two resonances have a so dissimilar influence on Saturn’s spin-axis dynamics. This phenomenon has been further discussed by Saillenfest et al. (2020) in the case of Jupiter.

Fig. 4 Examples illustrating the four different types of past obliquity evolution of Saturn. Each graph shows a 4-Gyr numerical trajectory (black dots) computed for Titan’s nominal migration rate and for a given value of Saturn’s normalised polar moment of inertia specifiedin title. Today’s location of Saturn is represented by the big green spot; the vertical error bar corresponds to our full exploration interval of λ. The red curves show the centre of first-order secular spin-orbit resonances (Cassini state 2) and the coloured areas represent their widths (same as Fig. 1). The top large area is the resonance with ϕ3 and the bottom thin area is the resonance with ϕ19 (see Table 1). The separatrices of the ϕ3 resonance are highlighted in blue. Going forwards in time, the trajectories go from bottom to top.

3.2 Adiabaticity of Titan’s migration

If the drift of α over time was perfectly adiabatic (i.e. infinitely slow compared to the oscillations of σ₃), the outcomeof the dynamics would not depend on the exact migration rate of Titan; the latter would only affect the evolution timescale. In the vicinity of Titan’s nominal migration rate, Saillenfest et al. (2021) show that the drift of α is almost an adiabatic process. Here, we extend the analysis to a larger interval of migration rates in order to determine the limits of the adiabatic regime.

Figure 7 shows Saturn’s obliquity 4 Gyr in the past obtained by backward numerical integrations for different migration rates of Titan and using values of λ finely sampled in its exploration interval. Migration rates comprised between the red and magenta curves are compatible with the astrometric measurements of Lainey et al. (2020), and migration rates comprised between theblue and green curves are compatible with their radio-science experiments (same colour code as in Fig. 3). As argued by Saillenfest et al. (2021), Titan’s migration may have been sporadic, in which case b would vary with time and the result would roughly correspond to a mix of several panels in Fig. 7. However, because of our current lack of knowledge about tidal dissipation mechanisms, refined evolution scenarios would only be speculative at this stage.

The blue curve of Fig. 7 confirms that the nominal migration rate of Lainey et al. (2020) is close to the adiabatic regime, since smaller rates give very similar results (see the curves for a migration two times and four times slower). Non-adiabatic signatures are only substantial in the grey areas, that is, for recently captured trajectories that crossed the resonance separatrix (evolution Type 2). Indeed, the teeth-shaped structures are due to ‘phase effects’, meaning that the precise outcome depends on the value of the resonant angle σ₃ during the separatrix crossing. For smaller migration rates, these structures pack together and tend to a smooth interval (that would be reached for perfect adiabaticity). If the migration of Titan is very slow, however, our 4-Gyr backward integrations stop while Saturn is still close to the resonance, or even inside it. The curves obtained for b ≲ 1∕7 b₀ have not enough time to completely relax from their initial shape shown in Fig. 5. This means that if, as argued by Saillenfest et al. (2021), Titan is responsible for Saturn’s current large obliquity, its migration cannot have been arbitrarily slow. Historical tidal models used to predict very small migration rates, as in the top left panel of Fig. 7. Such small migration rates are unable to noticeably affect Saturn’s obliquity over the age of the Solar System. This explains why previous studies considered that Saturn’s precession constant remained approximatively constant since the late planetary migration (Boué et al. 2009; Brasser & Lee 2015; Vokrouhlický & Nesvorný 2015). Figure 7 shows that for λ ∈ [0.200, 0.240], near-zero past obliquities can be achieved only if b ≳ 1∕16 b₀, that is, if Titan migrated by at least 1 R_eq after the late planetary migration. This condition is definitely achieved in the whole error ranges given by Lainey et al. (2020), provided that Titan’s migration did go on over a significant amount of time. Assuming that b = b₀, Titan should have migrated at least during several hundreds of million years before today in order for its semi-major axis to have changed by more than 1 R_eq. On the contrary, no substantial obliquity variation could be produced if Titan only began migrating very recently (less than a few hundreds of million years) and always remained unmoved before that. As mentioned by Saillenfest et al. (2021), this extreme possibility appears unlikely but cannot be ruled out yet.

When we increase Titan’s migration rate above its nominal value, Fig. 7 shows that the adiabatic nature of the drift of α is gradually destroyed. For b = 3b₀, phase effects become very strong and distort the whole picture. The magenta curve (which marks the limit of the 3σ error bar of Lainey et al. 2020) shows that the non-adiabaticity allows for a past obliquity of Saturn equal to exactly 0°. Such a null value is obtained when the oscillation phase of σ₃ brings Saturn’s obliquity to zero exactly together with Cassini state 2. This configuration can only happen for finely tuned values of the parameters, which is why putting a primordial obliquity ε ≈ 0° as a prerequisite puts so strong constraints on the parameter range allowed (Brasser & Lee 2015; Vokrouhlický & Nesvorný 2015).

If the resonance crossing is too fast, however, the resonant angle σ₃ has not enough time to oscillate before escaping the resonance. As a result, Saturn’s spin-axis can only follow the drift of the resonance centre during a very limited amount of time, and only a moderate obliquity kick is possible. As discussed in Sect. 3.1, this is what happens for the thin resonance with ϕ₁₉. In Fig. 7, the effect of overly fast crossings is clearly visible for b ≳ 11b₀. Beyond this approximate limit, all trajectories in our backward integrations cross the resonance separatrix, which means that trajectories of Type 3 are impossible and no small past obliquity can be obtained.

Figure 8 summarises all values of Saturn’s past obliquity obtained in our backward integrations as a function of Titan’s migration rate and Saturn’s polar moment of inertia. Non-adiabaticity is revealed by the coloured waves, denoting phase effects. As expected, the waves disappear for b ≲ b₀: this is the adiabatic regime (see Fig. 4 of Saillenfest et al. 2021 for a zoom-in view). For very small migration rates, however, Titan would not have time in 4 Gyr to migrate enough to produce substantial effects on Saturn’s obliquity. This is why the dark-blue region in Fig. 8 does not reach b = 0. For too fast migration rates, on the contrary, the resonance crossing is so brief that it can only produce a small obliquity kick. In particular no past obliquity smaller than 5° is obtained for b ≳ 10 b₀. This migration rate can therefore be considered as the largest one allowing Titan to be held responsible for Saturn’s current large obliquity.

Fig. 5 Current dynamics of Saturn’s spin axis according to its normalised polar moment of inertia λ. The value of λ (top horizontal axis) is linked to the current precession constant of Saturn (bottom horizontal axis) through Eqs. (3) and (5). The black interval shows the ‘instantaneous’ oscillation range of Saturn’s spin axis (i.e. without drift of α) obtained by numerical integration. The resonant angle is σ3 = ψ + ϕ3 (see Sect. 2). The green line shows Saturn’s current obliquity and resonant angle. The background colour indicates the type of past evolution as labelled in the top panel (see text for the numbering).

Fig. 6 Period of small oscillations about the resonance centre for a resonance with ϕ3 or ϕ19. The resonant angles are σ3 = ψ + ϕ3 and σ19 = ψ + ϕ19, respectively.Dashed curves are used for oscillations about Cassini state 2 before the separatrix appears. The appearance of the separatrix is marked by a blue dot.

3.3 Extreme phase effects

As can be guessed from the thinness of the spikes visible in some panels of Fig. 7, the variety of outcomes obtained for trajectories that cross the resonance separatrix (i.e. Types 1 and 2) depend on the resolution used for sampling the parameter λ. The deepestspikes denote the strongest phase effects; they correspond to trajectories that reach the resonance almost exactly at its hyperbolic equilibrium point (called Cassini state 4: see e.g. Saillenfest et al. 2019 for phase portraits⁵). Since the resonance island slowly drifts as α varies over time, extreme phase effects can be produced when the hyperbolic point drifts away just as the trajectory gets closer to it, maintaining the trajectory on the edge between capture and non-capture into resonance. This kind of borderline trajectory is more common for strongly non-adiabatic drifts (i.e. the spikes in Fig. 7 are wider for larger b), because a faster drift of the resonance means that trajectories need to follow less accurately the separatrix in order to ‘chase’ the hyperbolic point at the same pace as it gets away. If the drift of the resonance is too fast, however, trajectories are outrun by the resonance and strong phase effects are impossible. This is visible in the last panel of Fig. 7 (for b = 20 b₀), in which the spikes are noticeably smoothed.

In order to investigate the outcomes of extreme phase effects, one can look for the exact tip of the spikes in Fig. 7 by a fine tuning of λ. For b = b₀ (central panel), a tuning of λ at the 10⁻¹⁵ level shows that Type 2 trajectories all feature a minimum past obliquity of about 10°, as illustrated in Fig. 9. This minimum value is the same for each spike, and zooming in in Fig. 9 shows that we do reach the bottom of the spikes. For Type 1 trajectories (i.e. λ < 0.220 in the central panel of Fig. 7), we managed to find past obliquities of about 28° at the tip of the spikes, but using extended precision arithmetic may allow one to obtain even smaller values (possibly down to 10° as for Type 2 trajectories). The width of these spikes (Δλ < 10⁻¹⁵) would however make them absolutely invisible in Figs. 7 and 8. In fact, the level of fine tuning required here is so extreme that such trajectories are unlikely to have any physical relevance. They are yet possible solutions in a mathematical point of view. Some examples are given in Appendix B.

These findings can be compared to previous studies, even though previous studies relied on a different tilting scenario. For a non-adiabatic drift of the resonance and a past obliquity fixed to 1.5°, Boué et al. (2009) found that if Saturn is not currently inside the resonance (i.e. if λ < 0.220), an extremely narrow but non-zero range of initial conditions is able to reproduce Saturn’s current orientation, with a probability less than 3 × 10⁻⁸. Using a smaller set of simulations, Vokrouhlický & Nesvorný (2015) did not even find a single of these trajectories. In light of our results, we argue that these unlikely trajectories are produced through the ‘extreme phase effects’ described here. The vanishingly small probability of producing such trajectories is confirmed in Sect. 4.

Fig. 7 Past obliquity of Saturn for different migration rates of Titan. The top and bottom horizontal axes are the same as in Fig. 5 and the horizontal green line shows Saturn’s current obliquity. For a given value of the normalised polar moment of inertia λ (top horizontal axis), the curve width shows the oscillation range of obliquity 4 Gyr in the past obtained by backward numerical integration. The migration rates are labelled on each panel as a fraction of the nominal rate of Lainey et al. (2020). The four coloured curves correspond to the migration rates illustrated in Fig. 3. The grey stripes in the central panel highlight trajectories of Type 2 (same as in Fig. 5). The value of b in the top left panel corresponds to a current quality factor Q equal to 5000 (see Lainey et al. 2020).

Fig. 8 Past obliquity of Saturn as a function of Titan’s migration velocity and Saturn’s polar moment of inertia. Each panel of Fig. 7 corresponds here to a vertical slice. The colour scale depicts the minimum obliquity of the oscillation range, and the white curve highlights the 5° level. The 3σ uncertainty ranges of Lainey et al. (2020) yield today approximately b∕b0 ∈ [1∕2, 5] for the astrometric measurements and b∕b0 ∈ [1, 3∕2] for the radio-science experiments (see Fig. 3).

Fig. 9 Zoom-in view of the central panel of Fig. 7. We use a red curve to highlight the bottom limit of the blue interval, otherwise the narrowness of the spikes makes them invisible (the width of spike d is Δλ ≈ 10−14). This graph can be compared to Fig. 3 of Saillenfest et al. (2021), where such level of fine tuning is not shown due to its questionable physical relevance. See Appendix B for examples of trajectories.

Fig. 10 Snapshots of a Monte Carlo experiment computed for λ = 0.204 and Titan’s nominal migration rate. This experiment features 101 values of initial obliquity between 0° and 60°, for which 240 values of initial precession angle are regularly sampled in [0, 2π). The value of α reached by the trajectories at the time of the snapshot is labelled on each panel. Each trajectory is represented by a small dot which is coloured according to the variation range of the resonant angle σ3 (obtained bya 0.5-Gyr numerical integration with constant α). The horizontal green line shows the current obliquity of Saturn. At the beginning of the propagations, all trajectories are coloured red (since σ3 circulates), and distributed along a diagonal line. Then, as α increases over time, trajectories are dispersed off the diagonal according to the four types of trajectories depicted in Fig. 4 and labelled in the penultimate panel.

4 Monte Carlo search for initial conditions

In Sect. 3, the past behaviour of Saturn’s spin axis has been investigated using backward numerical integrations. If we now consider the space of all possible orientations for Saturn’s primordial spin axis, each dynamical pathway has a given probability of being followed. A large subset of trajectories (those of Types 1 and 2) go through the separatrix of the large resonance with ϕ₃. Separatrix crossings are known to be chaotic events (see e.g. Wisdom 1985), and since Saturn’s orbital evolution is not restricted to its 3rd harmonic, the separatrix itself appears as a thin chaotic belt (see e.g. Saillenfest et al. 2020). Therefore, we can wonder whether the chaotic divergence of trajectories during separatrix crossings could lead to some kind of time-irreversibility in our numerical solutions (see e.g. Morbidelli et al. 2020), especially in the non-adiabatic regime, which has not been studied by Saillenfest et al. (2021). These aspects can be investigated through a Monte Carlo search for the initial conditions of Saturn’s spin axis.

4.1 Capture probability

Our first experiment is designed as follows: for a given set of parameters (b, λ), values of initial obliquity are regularly sampled between 0° and 60°. Then, for each of those, we regularly sample values of initial precession angle ψ ∈ [0, 2π), and all trajectories are propagated forwards in time starting at 4 Gyr in the past (i.e. after the late planetary migration) up to today’s epoch. Figure 10 shows snapshots of this experiment for λ = 0.204 and Titan’s nominal migration rate (b = b₀). The first snapshot is taken about 20 million years after the start of the integrations, and the last snapshot is taken at today’s epoch. Changing the value of λ produces a shift of Saturn’s precession constant α but no strong variation in its drift rate (see Fig. 2). Moreover, since this drift is almost an adiabatic process for b = b₀ (see Sect. 3.2), a small change of drift rate does not modify the statistical outcome of the dynamics but only its timescale. For these reasons, a snapshot in Fig. 10 taken at a given time t for λ = 0.204 is undistinguishable from a snapshot taken at a slightly different time for another value . More precisely, if we introduce a function of time f_λ(t) such that t→α = f_λ(t), an indistinguishable snapshot is obtained for a polar moment of inertia at a time . Hence, the only parameter that matters here is the value of the precession constant α reached by the trajectories. This is why the panels of Fig. 10 are labelled by α instead of t: this way they are valid for any value of λ.

Before reaching the neighbourhood of the resonance with ϕ₃, Fig. 10 shows that all trajectories only slightly oscillate around their initial obliquity value (compare the first two snapshots, taken for two very different values of α). Then, as α continues to increase, the trajectories are gradually divided between the four possible types of evolution listed in Sect. 3.1. All trajectories with initial obliquity smaller than about 10° are captured in the resonance and lifted to high obliquities (Type 3: blue dots). Trajectories with a larger initial obliquity can either be captured (Type 2: green and orange dots) or go past the resonance through its hyperbolic point (Type 1: lowermost red dots).

Assuming that Saturn’s primordial precession angle ψ is a random number uniformly distributed in [0, 2π), the probability of capture in resonance is given by the fraction of points ending up in the pencil-shaped structure of Fig. 10. The result is shown in Fig. 11, in which we increased the resolution for better statistical significance. Assuming perfect adiabaticity, each outcome can be modelled as a probabilistic event ruled by analytical formulas (see Henrard & Murigande 1987; Ward & Hamilton 2004; Su & Lai 2020). As shown by Fig. 11, non-adiabaticity tends to smooth the probability profile and to reduce the interval of 100%-sure capture. Forgrowing initial obliquity, the probability of Type 2 trajectories (i.e. capture) decreases, favouring Type 1 trajectories instead (i.e. crossing without capture).

Fig. 11 Capture probability of Saturn in secular spin-orbit resonance with ϕ3 as a functionof its primordial obliquity. For each initial obliquity (401 values between 0 and 60°), 720 values of initial precession angle are uniformly sampled in [0, 2π) and propagated forwards in time starting from − 4 Gyr and until every trajectory has reached the resonance. This experiment is repeated with two different migration laws for Titan (see labels). The result is virtually independent of the value chosen for λ. The fractionof captured trajectories (coloured curves) is compared to the perfect adiabatic case (black curve) computed with the analytical formulas of Henrard & Murigande (1987).

4.2 Loose success criteria: probing all dynamical pathways

Over all possible trajectories, we now look for those matching Saturn’s actual spin-axis dynamics today. We first consider ‘loose success criteria’, for which a run is judged successful if: (i) Saturn’s current obliquity ε = 26.727° lies within the final spin-axis oscillation interval, and (ii) the libration amplitude of the resonant angle σ₃ lies within 5° of the actual amplitude shown in Fig. 5. These criteria are not chosen to be very strict in order to probe all dynamical pathways in the neighbourhood of Saturn’s spin-axis orientation, including some that could have been missed by the backward propagations of Sect. 3. Our results are depicted in Fig. 12. We closely retrieve the predictions of backward numerical integrations, in particular for trajectories of Type 3. Narrowing the target interval leads to an even better match. For trajectories of Type 2 (grey background), we obtain a larger spread of initial obliquities because our success criteria do not include any restriction on today’s phase of the resonant angle σ₃, but only on its oscillation amplitude. The results shown in Fig. 12 are therefore less shaped by ‘phase effects’ discussed in Sect. 3.2. Since a slight change in Titan’s migration rate would result in a phase shift, we can interpret Fig. 12 as encompassing different migration rates around the nominal observed rate. The results presented in Fig. 12 are therefore more general than those obtained using backward numerical integrations. In accordance with Fig. 11, the success ratio for Type 2 trajectories sharply decreases for increasing initial obliquity (colour gradient), because most initial conditions lead to a resonance crossing without capture. Moreover, we do not detect trajectories as extreme as those presented in Sect. 3.3 (i.e. with an initial obliquity of about 10° all over the width of Zone 2), because they require initial conditions that are too specific for our sampling; the probability of obtaining one is indeed negligible. Finally, for trajectories of Types 1 and 4, which are today out of the resonance, our ‘loose success criteria’ are extremely permissive, since the variation amplitude of σ₃ is 2π for all trajectories (see Fig. 5). This explains why Fig. 12 shows large intervals of black dots. These intervals can be spotted in Fig. 10, where they appear as the whole range of red dots that are pierced by the green horizontal line.

Figure 12 does not feature unexpected dynamical paths that could have been missed by our backward integrations, even though signatures of chaos are visible in the sparse spreads of coloured dots. From this close match, we conclude that the chaos is not strong enough here to significantly mingle the trajectories and to produce a substantial phenomenon of numerical irreversibility. As one can point out, separatrix crossings would have been irreversible if, in order to predict the different outcomes, we used the adiabatic invariant theory instead of numerical integrations (see Henrard & Murigande 1987; Ward & Hamilton 2004; Su & Lai 2020). Indeed, in the adiabatic invariant theory, the resonant angle is assumed to oscillate infinitely faster than the drift of α and phase effects are modelled as probabilistic events (Henrard 1982, 1993). This probabilistic modelling of chaos explains why separatrix crossings are not reversible when using this theory.

Fig. 12 Brute-force search for Saturn’s past obliquity using the loose success criteria: probing all dynamical pathways. We use Titan’s nominal migration law (b = b0). Among all trajectories evenly sampled in the space of the normalised polar moment of inertia λ (top horizontal axis, 101 values), of the initial obliquity (vertical axis, 101 values), and of the initial precession angle (240 values between 0 and 2π), we only keep those matching Saturn’s spin axis today according to our loose success criteria (see text). Each point is coloured according to the number of successful runs among the 240 initial precession angles; the success ratio is written below the colour bar. A point is not drawn if no successful trajectory is found. In the back, the blue interval shows the past obliquity of Saturn obtained by backward numerical integration (same as Fig. 7 for b = b0), showing theconsistency between backward and forward integrations in time. The background stripes and their labels have the same meaning as in Fig. 5.

4.3 Strict success criteria: relative likelihood of producing Saturn’s current state

In order to compare the likelihood of producing Saturn’s current state in the space of all possible initial conditions, our loose success criteria are not enough. Independently of whether Saturn is inside or outside the resonance today, its spin-axis precession is not uniform, which means that the phase of Saturn’s spin-axis motion at a given time is not uniformly distributed and must therefore be taken into account, too. Moreover, we saw in Sect. 3.2 that out of the strict adiabatic regime, phase effects (that are deliberately ignored by our loose success criteria) do matter to reproduce Saturn’s current spin-axis orientation; actually, the very notion of ‘libration’ loses its meaning when the drift of α is not adiabatic, since the resonance is distorted before σ₃ has time to perform a single cycle. For these reasons, we now define ‘strict success criteria’, for which a run is judged successful if: (i) today’s obliquity ε lies within 0.5° of the true value, and (ii) today’s precession angle ψ lies within 5° of the true value. These criteria are very narrow, but still within reach of our millions of numerical propagations. The result is shown in Fig. 13 for Titan’s nominal migration rate. As expected, the points are more sparse than in Fig. 12 and the success ratios are smaller. Assuming that Saturn’s primordial precession angle is a random number uniformly distributed between 0 and 2π, the colour gradient in Fig. 13 is a direct measure of the likelihood to reproduce Saturn’s current state. Type 3 trajectories are greatly favoured: they feature the maximum likelihood, which is about ten times the likelihood of Type 1 trajectories. The region with maximum likelihood is for past obliquities between about 2° and 7°, and current precession constant α between about 0.76 and 0.79″ yr⁻¹ (red box).

As already discussed by Ward & Hamilton (2004) and Hamilton & Ward (2004), there are two reasons why Type 3 trajectories are the most likely: first, they have a 100% chance of being captured inside the resonance (whereas Types 1 and 2 both have a non-zero probability of failure, see Fig. 11); second, Type 3 trajectories oscillate today with a small amplitude inside the resonance, which means that all of them feature a similar value of the precession angle ψ, imposed by the resonance relation σ₃ ~ 0. On the contrary, other types of trajectories either feature a large oscillation amplitude of σ₃ (Type 2) or circulation of σ₃ (Types 1 and 4); therefore, they only sweep over Saturn’s actual orientation once in a while, and matching it today would only be a low-probability ‘coincidental’ event⁶. As shown by Fig. 13, the least favoured trajectories are those of Type 2, especially for high initial obliquities, because of the strong decrease in capture probability (see Fig. 11).

In order to explore all migration rates and bring further constraints on the model parameters, we now turn toa second Monte Carlo experiment, with the following approach: assuming that Saturn was indeed tilted as a result of Titan’s migration, we look for the possible values of the parameters (b, λ) allowed, with their respective likelihood. This approach is similar to those used in previous studies (e.g. Vokrouhlický & Nesvorný 2015).

The notion of likelihood associated with this second experiment deserves some comments. Since Saturn’s spin axis performed many precession revolutions in 4 Gyr and since it was initially not locked in resonance, a tiny error in the model rapidly spreads over time into a uniform probability distribution of the precession angle ψ in [0, 2π). This is the reason why, in absence of any mechanism able to maintain ψ in a preferred direction, it is legitimate to consider a uniform initial distribution for ψ, as people usually do (and as we already did above). Establishing a prior distribution for ε, instead, is more hazardous: we know that near-zero values are expected from formation models, but small primordial excitations cannot be excluded. Such excitations could be attributed to the phase of planetesimal bombardment at the end of Saturn’s formation or by abrupt resonance crossings stemming from the dissipation of the protoplanetary and/or circumplanetary discs (see e.g. Millholland & Batygin 2019). Therefore, we arbitrarily consider here values of initial obliquity ε ≲ 5°, which leaves room for a few degrees of primordial obliquity excitation. This choice is somewhat guided by the 3°-obliquity of Jupiter, a part of which could possibly be primordial (Ward & Canup 2006; Vokrouhlický & Nesvorný 2015). Jupiter is located today near a secular spin-orbit resonance with s₇ (see Table 1), but contrary to Saturn, its satellites did not migrate enough yet to substantially increase its obliquity (Saillenfest et al. 2020); however, in order to ascertain possible values for Jupiter’s primordial obliquity, the effect of the past migration of the Galilean satellites would need to be studied. We choose to use a uniform random distribution of ε, resulting in a non-uniform distribution of spin-axis directions over the unit sphere that favours small obliquities⁷. The influence of our arbitrary choice of Saturn’s initial obliquity is discussed below.

In practice, our setup is the following: over a grid of point (b, λ) of the parameter space, we perform each time 2400 numerical integrations starting from random initial conditions (ε, ψ) with ε ≤ 5° and ψ ∈ [0, 2π). All trajectories are then propagated from − 4 Gyr up to today’s epoch, and we only keep trajectories matching Saturn’s current spin-axis orientation according to our strict success criteria. Figure 14 shows the result of this experiment. Again, we closely retrieve the predictions of backward integrations from Sect. 3, confirming the reversible nature of the dynamics, and helping us to interpret the patterns obtained. The wavy structure at 3b₀ ≲ b ≲ 5b₀ resembles tosome extent the successful matches of Vokrouhlický & Nesvorný (2015), reminding us that the basic dynamical ingredients are the same, even though the mechanism producing the resonance encounter in their study is different (their Fig. 7 is rotated clockwise). Unsurprisingly, the highest concentrations of matching trajectories in Fig. 14 are located in the regions where backward propagations result in near-zero primordial obliquities (compare with Fig. 8). The maximum likelihood thus favours slightly non-adiabatic migration rates, for b lying roughly between 3b₀ and 6b₀. According to Lainey et al. (2020), such values are consistent with the 3σ uncertainty ranges of Titan’s current migration rate obtained from astrometric measurements (b∕b₀ ∈ [1∕2, 5]), but not with the uncertainty ranges given by radio-science experiments (b∕b₀ ∈ [1, 3∕2]). However, successful trajectories with substantial likelihood are anyway found in a very large interval of migration rates, which extends much farther than the uncertainty range of Lainey et al. (2020). We can therefore not bring any decisive constraint on Titan’s migration history that would be tighter than those obtained from observations. Yet, the tilting of Saturn would impose strong constraint on Saturn’s polar moment of inertia. As already visible in the figures of Saillenfest et al. (2021), tilting Saturn from ε ≲ 5° in the adiabatic regime would require that λ lies between about 0.228 and 0.235. Figure 14 shows that allowing for non-adiabatic effects (b ≳ 3b₀) widens this range to about [0.224, 0.239].

Interestingly, Fig. 14 features three trajectories affected by an ‘extreme phase effect’ (see Sect. 3.3), visible as the three isolated grey points at λ ≈ 0.205, 0.210, and 0.215. These trajectories are of Type 1 (i.e. currently out of the resonance) and fit our strict success criteria. The existence of these points recalls that such trajectories are extremely rare (we found only three over millions of trials), but yet possible, as previously reported by Boué et al. (2009). They correspond to the narrow dark edges of the blue stripes in the top portion of Fig. 8. The complete trajectory producing the leftmost of these points can be found in Appendix B.

As mentioned above, the likelihood measure depicted in Fig. 14 is conditioned by our assumptions about the initial value of ε. The influence of these assumptions can be investigated from our large simulation set. Figure 15 shows the statistics restricted to the lower- and higher-obliquity halves of the distribution. Restricting the initial obliquity to ε ≲ 2.5° suppresses most successful matches from the adiabatic regime (b ≲ 3b₀), as one could have guessed from previous figures. On the contrary, restricting the statistics to the upper half of the distribution (ε ∈ [2.5°, 5°]) greatly shifts the point of maximum likelihood towards the adiabatic regime. Further experiments are provided in Appendix C with initial obliquity values up to 10°. These experiments show that the adiabatic and non-adiabatic regimes are roughly equally likely if one considers an isotropic distribution of initial spin orientations (with ε ≲ 5° or ε ≲ 10°) instead of a distribution favouring small initial obliquities as in Fig. 14. Unsurprisingly, the adiabatic regime and Titan’s nominal migration rate are the most likely if one considers initial obliquity values as 2° ≲ ε ≲ 7° (i.e. in the red box of Fig. 13).

This discussion shows how important is the prior chosen for the initial conditions. Assumption biases are unavoidable and were also present in previous studies: Boué et al. (2009) assumed ε = 1.5°; Vokrouhlický & Nesvorný (2015) assumed ε = 0.1° (with respect to the orbit averaged over all angles but ϕ₃); and Brasser & Lee (2015) assumed ε ≈ 0.05° (with respect to the invariable plane, i.e. the orbit averaged over all ϕ_k). As shown by our results, leaving room for a few degrees of extra primordial excitation, or even only 0.5°, in any of thosestudies could have greatly enhanced the chances of success. As recalled above, a few degrees of primordial obliquity excitation are plausible and could be explained in different ways. In this regard, the most general overview of our findings is given by Fig. 8, since it does not presuppose any initial obliquity for Saturn, and Fig. 13 shows the respective likelihood of each dynamical pathway, still with no assumption about the initial obliquity.

Fig. 13 Same as Fig. 12, but using the strict success criteria: comparing the relative likelihood of producing Saturn’s current state. As in Fig. 12, each point of the graph is made of 240 simulations with initial ψ ∈ [0, 2π). The red rectangle highlights the region featuring the highest success ratios.

Fig. 14 Distribution of the solutions starting from a low primordial obliquity and matching our strict success criteria. For each set (b, λ) of the parameters, 2400 values of initial obliquity ε and precession angle ψ are drawn from a uniform random distribution in (ε, ψ) ∈ [0°, 5°] × [0, 2π). Coloured dots show the parameter sets (b, λ) for which at least one successful trajectory was found; the success ratio is written below the colour bar. Light-grey crosses mean that no successful trajectory was found over our 2400 initial conditions. The black contours show the 5°-level obtained through backward numerical integrations (same as Fig. 8), showing the consistency between backward and forward integrations in time. The black lines in the top portion show the approximate location of the border of the blue stripes in Fig. 8, where extreme phase effects can happen; the corresponding ranges of parameters are so narrow that they are missed by the resolution of Fig. 8 (see Sect. 3.3).

Fig. 15 Same as Fig. 14, but for statistics based on a sub-sample of simulations. (a) Initial conditions in (ε, ψ) ∈ [0°, 2.5°] × [0, 2π). (b) Initial conditions in (ε, ψ) ∈ [2.5°, 5°] × [0, 2π). In both panels, each point is made of about 1200 initial conditions extracted from the simulations from Fig. 14.

5 The future obliquity of Saturn

Since Titan goes on migrating today, Saturn’s obliquity is likely to continuously vary over time. Hence, we can wonder whether it could reach large values, in the same way as Jupiter (Saillenfest et al. 2020). In order to explore the future obliquity dynamics of Saturn, we propagate Saturn’s spin-axis from today up to 5 Gyr in the future.

Figure 16 shows the summary of our results for finely sampled values of λ and b. Contrary to Fig. 8, we restrict here our sampling to b < 3b₀ because for larger migration rates, Titan goes beyond the Laplace radius during the integration timespan (a₆ ≈ 40 R_eq) and the close-satellite approximation used in Eq. (6) is invalidated. Faster migration rates are anyway disfavoured by the 3σ uncertainty range of the radio-science experiments of Lainey et al. (2020).

The top portion of Fig. 16 features trajectories of Type 1. Such trajectories are currently above the resonance with ϕ₃ (see Fig. 4) and they go farther away from it as α continues to increase. The increase in α makes them cross the resonances with ϕ₅₁, with ϕ₁₄, and with ϕ₁₅ (see Fig. 1 and Table 1). Being very small, these resonances are crossed quickly and they donot produce noticeable obliquity variations in Fig. 16. This explains why the top portion of the figure is coloured almost uniformly with an obliquity value approximatively equal to today’s. For b ≈ 3b₀ (the fastest migration presented in Fig. 16), trajectories reach the lower fringe of the strong resonance with ϕ₄ at the end of the integration, but they do not actually reach it.

The middle portion of Fig. 16 features trajectories of Types 2 and 3. Such trajectories are currently inside the resonance with ϕ₃ (see Fig. 4) and they follow its centre as α increases. After 5 Gyr from now, Saturn can therefore reach very large obliquity values, provided that Titan goes on migrating as predicted by Lainey et al. (2020). For migration rates lying in the error range of the radio-science experiments of Lainey et al. (2020), Saturn’s obliquity can grow as large as 65°. As noticed by Saillenfest et al. (2021), the resonance width increases up to α ≈0.971″ yr⁻¹, but decreases beyond. The trajectories featuring a large libration amplitude and a large increase in α have therefore a risk of being expelled out of the resonance, as described by Su & Lai (2020). After a careful examination, we found that expulsion out of resonance only occurs for the largest migration velocities and in a tiny interval of λ located at the very edge of the brightly-coloured region of Fig. 16. An example of such a trajectory is presented inFig. 17. The expelled trajectories reach slightly smaller values of obliquity than if they continued to follow the resonance centre; however, this behaviour concerns such a small range of parameters (which is almost undistinguishable in Fig. 16) that it is very unlikely to have any consequence for Saturn.

The bottom portion of Fig. 16 features trajectories of Type 4, which are ruled out by our uncertainty range for λ ∈ [0.200, 0.240]. Such trajectories did not reach yet the resonance today (see Fig. 4), but they will in the future as α continues toincrease. The resonance encounter can either lead to a capture (like Type 2 trajectories) or to a permanent decrease in obliquity (like Type 1 trajectories). The outcome is determined by the phase of σ₃ at crossing time, which depends on the migration velocity. This explains why the two possible outcomes are organised in Fig. 16 as narrow bands that are close to each other for a slow migration and spaced for a fast migration. In a perfect adiabatic regime, the bands would be so close to each other that the outcome could be modelled as a probabilistic event.

Fig. 16 Future obliquity of Saturn as a function of Titan’s migration velocity and Saturn’s polar moment of inertia. The axes are the same as in Fig. 8. The 3σ uncertainty ranges of Lainey et al. (2020) yield approximately b∕b0 ∈ [1∕2, 5] for the astrometric measurements and b∕b0 ∈ [1, 3∕2] for the radio-science experiments. Some level curves are shown in red.

Fig. 17 Example of Type 2 trajectory that is expelled out of the resonance. It has been obtained for λ = 0.221076 and a migration rate b∕b0 = 3. The integration runs from today up to 5 Gyr in the future.

6 Discussion and conclusion

Since giant planets are expected to form with near-zero obliquities, some mechanism must have tilted Saturn after its formation. Saillenfest et al. (2021) have shown that the fast migration of Titan measured by Lainey et al. (2020) may be responsible for the current large obliquity of Saturn. Through an extensive set of numerical simulations, we further investigated the long-term spin-axis dynamics of Saturn and determined the variety of laws for Titan’s migration compatible with this scenario.

Saturn is located today near a strong secular spin-orbit resonance with the nodal precession mode of Neptune (Ward & Hamilton 2004). As Titan migrates over time, it produces a drift in Saturn’s spin-axis precession velocity, which led Saturn to encounter this resonance. The continuous migration of Titan shifts the resonance centre over time, which can force Saturn’s obliquity to vary. Through this mechanism, Saturn’s obliquity can have grown from a small to a large value provided that: (i) Titan migrated over a large enough distance to substantially shift the resonance centre, and (ii) Titan migrated slowly enough for Saturn to adiabatically follow the resonance shift. The first condition is met if Titan migrated over a distance of at least one radius of Saturn after the late planetary migration, more than 4 Gyr ago. Assuming that Titan’s migration is continuous, this requires migration velocities larger than about n ≈ 0.06 times the nominal rate given by Lainey et al. (2020). For comparison, astrometric measurements predict n ≳ 0.5. The second condition is met if Titan’s migration velocity does not exceed n ≈ 10 times the nominal rate, while astrometric measurements predict n ≲ 5. Therefore, the scenario proposed by Saillenfest et al. (2021) is realistic over the whole range of migration rates obtained from observations. It even allows for more complex scenarios in which Titan would alternate between fast and slow migration regimes.

For the largest migration rates of Titan allowed by observational uncertainty ranges, non-adiabatic effects are quite pronounced, but not to the point of preventing Saturn from following the resonance centre. Interestingly, non-adiabaticity even allows for an exactly zero value for Saturn’s primordial obliquity. Zero values are however disfavoured by the error range of radio-science experiments, which yield most likely primordial obliquities between 2° and 7°.

Our Monte Carlo experiments do not reveal a strong chaotic mixing of trajectories, even though borderline separatrix-crossing trajectories do exhibit a noticeable chaotic spreading. All possible dynamical paths fall into the four types of trajectories obtained by Saillenfest et al. (2021) through backward numerical integrations, and we detected no substantial numerical irreversibility. For Titan’s nominal migration rate, our experiments show that all trajectories with initial obliquity smaller than about 10° are captured inside the resonance with a 100% probability. Such trajectories can match Saturn’s current orientation if its normalised polar moment of inertia λ lies in about [0.224, 0.237], as previously reported. Interestingly, small past obliquities ε ≲ 10° in our Monte Carlo experiments also feature the highest likelihood of reproducing Saturn’s current spin-axis orientation, surpassing high-obliquity alternatives by a factor of about ten. Yet, other values of λ cannot be completely ruled out; they would mean that Saturn’s past obliquity was larger or similar as today and one would need to find another explanation for its large value.

In the future, the still ongoing migration of Titan is expected to produce dramatic effects on Saturn’s obliquity provided that Saturn is currently located inside the resonance, that is, if λ lies in about [0.220, 0.241]. Depending on the precise migration rate of Titan, Saturn’s obliquity would then range between 55° and 65° after 5 Gyr from now, and we even obtain values exceeding 75° when considering the full 3σ uncertainty of the astrometric measurements of Lainey et al. (2020). For smaller values of λ, Saturn’s obliquity is not expected to change much in the future because the migration of Titan pushes it away from the resonance. No strong obliquity variations would be expected either if Titan’s migration rate strongly drops in the future (i.e. if Titan is released out of the tidal resonance-locking mechanism of Fuller et al. 2016), but to our knowledge, there is no evidence showing that it could be the case.

The migration law for Titan proposed by Lainey et al. (2020) and used in this article is very simplified. Since our conclusions remain valid in a much larger interval of migration rates than allowed by the observational uncertainties, we can be confident that no major change would be produced by using different (and possibly more realistic) migrations laws, unless Titan underwent extreme variations in migration rate in the past. For instance, if Titan’s migration is not continuous and if it was only triggered very recently (less than a few hundreds of million years ago), then Saturn’s past obliquity dynamics would not have been affected. As mentioned by Saillenfest et al. (2021), this alternative is unlikely but cannot be ruled out considering our current knowledge of the tidal dissipation within Saturn.

The past and future behaviour of Saturn’s spin axis is very sensitive to its normalised polar moment of inertia λ. An accuracy of at least three digits would be required to securely assert which dynamical path was followed by Saturn and what will be the future evolution of its spin axis. Model-dependent theoretical values are not enough for this purpose, and it is still unclear what is the true uncertainty of values inferred from the Cassini data (Helled 2011; Fortney et al. 2018; Movshovitz et al. 2020). A precise value of λ would inform us about whether Saturn is currently inside the resonance (which is the most likely alternative), or outside the resonance. If Saturn is confirmed to be currently in resonance, it would imply that Titan’s past migration rate never became so fast as to eject Saturn from the resonance or to prevent its capture in the first place. However, this constraint would not be very stringent: simulations show that Saturn can be captured into resonance even if Titan’s migration rate is increased by a factor ten from the nominal measured value. If, on the contrary, Saturn turns out to be currently out of resonance, then it would imply that its primordial obliquity was high, and most probably even higher than 30°, regardless of Titan’s precise migration history. This last possibility is not what one would expect from planetary formation models, and our results show that it is also unlikely in a dynamical point of view.

Previous works reveal that numerous dynamical mechanisms can alter the obliquity of a planet (see e.g. Laskar & Robutel 1993; Correia & Laskar 2001; Quillen et al. 2018; Millholland & Batygin 2019). The fast migration of satellites and capture in a secular spin-orbit resonance offers one more alternative, and we have shown that it can result in a steady increase in obliquity, possibly lasting over the whole lifetime of the planetary system. In the broad context of exoplanets, we can therefore expect that only a few would have conserved their primordial axis tilt, whether they are close-in and likely tidally locked (Millholland & Laughlin 2019), or whether they are largely spaced and have very stable orbits like Jupiter and Saturn.

Acknowledgements

Our work greatly benefited from discussions with David Nesvorný; we thank him very much. We are also very grateful to Dan Tamayo for his in-depth review and inspiring comments. G.L. acknowledges financial support from the Italian Space Agency (ASI) through agreement 2017-40-H.0.

Appendix A Orbital solution for Saturn

The secular orbital solution of Laskar (1990) is obtained by multiplying the normalised proper modes and (Tables VI and VII of Laskar 1990) by the matrix corresponding to the linear part of the solution (Table V of Laskar 1990). In the series obtained, the terms with the same combination of frequencies are then merged together, resulting in 56 terms in eccentricity and 60 terms in inclination. This forms the secular part of the orbital solution of Saturn, which is what is required by our averaged model.

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

Appendix B Examples of trajectories featuring extreme phase effects

In Sect. 3.3, we show that trajectories crossing the separatrix can feature extreme phase effects when they reach the resonance in the vicinity of its hyperbolic point and follow its drift over time. This maintains them on the edge between capture (Type 2 trajectory) and non-capture (Type 1 trajectory).

Figure B.1 shows examples of such trajectories obtained for Titan’s nominal migration rate. These trajectories are of Type 2 (i.e. currently inside the resonance). Instead of the precession angle ψ, we plot the resonant angle σ₃ = ψ + ϕ₃, where ϕ₃ evolves as in Eq. (8). The elliptic point of the resonance (Cassini state 2) is located at σ₃ = 0, and the hyperbolic equilibrium point (Cassini state 4) is located at σ₃ = π (see e.g. Saillenfest et al. 2019). We see that passing from one spike of Fig. 9 to the next one corresponds to performing one more oscillation inside the resonance. For a purely adiabatic dynamics, all spikes would be infinitelyclose to each other, such that it would be impossible to get one specific trajectory by finely tuning λ.

Figure B.2 shows another example of extreme phase effect but for a trajectory of Type 1 (i.e. currently outside the resonance). It is obtained using a strongly non-adiabatic migration, which widens the parameter ranges allowing for extreme phase effects (see Sect. 3.3). This trajectory does not exactly match Saturn’s spin-axis orientation today, but it lies within our strict success criteria defined in Sect. 4.3: its current coordinates ε and ψ are within 0.4° and 4.9° of the actual ones, respectively. This trajectory appears in the top portion of Fig. 14 as the leftmost isolated grey point. It can be linked to the bottom of a spike in Fig. 7, that is, to one of the top blue stripes of Fig. 8. After having bifurcated away from the hyperbolic point, Fig. B.2 shows that this trajectory has performed one complete revolution of σ₃. The two other isolated grey points in Fig. 14 have performed zero and two, respectively.

Fig. B.1 Example of trajectories featuring an extreme phase effect. Left column: evolution of the obliquity, and right column: evolution of the resonant angle σ3 = ψ + ϕ3. The migration parameter is b = b0. For each row, the parameter λ used corresponds to the tip of a spike in Fig. 9 (see labels), tuned at the 10−15 level. The pink area represents the interval occupied by the resonance once the separatrix appears. The blue curve shows the locationof the hyperbolic equilibrium point (Cassini state 4). The green point shows Saturn current location (at t = 0).

Fig. B.2 Same as Fig. B.1, but for a trajectory of Type 1. This trajectory has a migration parameter b = 7.37 b0 and a normalised polar moment of inertia λ = 0.2114.

Appendix C Experiments on the initial obliquity prior

In Sect. 4.3, a Monte Carlo experiment is performed in order to look for the most likely values of Saturn’s precession constant and Titan’s migration rate. Formation models predict that Saturn’s primordial obliquity was near-zero, but the statistics obtained greatly depend on the precise distribution used as initial conditions. In this section, we investigate further this dependence with additional Monte Carlo experiments.

Figure C.1 shows the statistics obtained when assuming a uniform distribution of initial conditions over the spherical cap defined by ε ≤ 5° (i.e. with a uniform sampling of cosε instead of ε). Contrary to Fig. 14, this distribution is isotropic: it assumes that all directions over the spherical cap are equiprobable; small obliquity values are not particularly favoured. In practice, we can avoid running millions of simulations again by simply weighting the count of each run in Fig. 14 by sin ε. As illustrated in Fig. C.2, this trick allows us to mimic a uniform distribution of cos ε from a uniform distribution of ε. This method has the drawback of reducing by roughly a factor two the resolution at the high-obliquity end of the distribution (since trajectories are weighted by a factor w ≈ 2), but this is not a problem here thanks to the high number of simulations.

Interestingly, Fig. C.1 shows that a uniform distribution of initial spin directions over the sphere yields approximatively equal likelihoods for the adiabatic (b ≲ 3b₀) and non-adiabatic (b ≳ 3b₀) regimes. If we enlarge the distribution of initial conditions to ε ≤ 10°, Fig. C.3 shows that the limits between the adiabatic and non-adiabatic regimes completely vanish, leaving only one large region with roughly constant likelihood.

Figure C.4 shows the distribution of successful runs starting from initial obliquities in the range 2.5° ≤ ε ≤ 7.5°. This interval turns out to be roughly the one that favours most the adiabatic regime and Titan’s nominal migration rate, to the detriment of the non-adiabatic regime. This is not surprising, since past obliquities between about 2° and 7° are the most likely for Titan’s nominal migration rate (see Fig. 13).

Fig. C.1 Same as Fig. 14, but considering an isotropic distribution of initial spin orientation with ε ≤ 5°. It is obtained from Fig. 14 by weighting the count of each run by sin ε (see Fig. C.2).

Fig. C.2 Sampled distribution of initial obliquity for one arbitrary point of Fig. 14, made of 2400 simulations. The raw count of sampled trajectories is shown in red; the weighted count is shown in blue. Top: histogram with respectto the obliquity. Bottom: histogram with respect to the cosine of obliquity. The probability density functions (‘pdf’) are shown by the red and blue curves.

Fig. C.3 Same as Fig. 14, but for a range of initial spin orientations enlarged to 0° ≤ ε ≤ 10°. Each point is made of 2400 numerical simulations. (a) Uniform random distribution of (ε, ψ) ∈ [0°, 10°] × [0, 2π). (b) Uniform random distribution of (ε, ψ) over the spherical cap defined by ε ≤ 10°. As in Fig. C.1, panel b is obtained by weighting the count numbers of panel a. The black contours show the 5° and 10° levels obtained through backward numerical integrations (see Fig. 8).

Fig. C.4 Same as Fig. 14, but for initial obliquities uniformly distributed between 2.5° and 7.5°. It is obtained from a sub-sample of Fig. C.3a, such that each point is made of about 1200 runs. The black contours show the 5° and 10° levels obtained through backward numerical integrations (same as Fig. C.3).

References

All Tables

All Figures

Fig. 1 Location and width of every first-order secular spin-orbit resonance for Saturn. Each resonant angle is of the form σp = ψ + ϕp where ϕp has frequency νp labelled on the graph according to its index in the orbital series (see Table 1 and Appendix A). For a given value of the precession constant α, the intervalof obliquity enclosed by the separatrix is shown in pink, as computed using the formulas of Saillenfest et al. (2019). The green bar shows Saturn’s current obliquity and the range for its precession constant considered in this article, as detailed in Sects. 2.3 and 2.4.

In the text

Fig. 2 Evolution of the effective precession constant of Saturn due to the migration of Titan (adapted from Saillenfest et al. 2021). The top and bottom green curves correspond to the two extreme values of the normalised polar moment of inertia λ considered in this article. They appear into α through Eq. (3). Both curves are obtained using the nominal value b = 1∕3 in Eq. (9). Today’s interval corresponds to the one shown in Fig. 1; it is independent of the value of b considered. The blue line shows Neptune’s nodal precession mode, which was higher before the end of the late planetary migration.

In the text

	Fig. 3 Time evolution of Titan’s semi-major axis for different migration rates. The pink and blue intervals show the 3σ uncertainty ranges of astrometric and radio-science measurements, respectively (Lainey et al. 2020). The coloured curves are obtained by varying the parameter b in Eq. (9).
In the text

Fig. 4 Examples illustrating the four different types of past obliquity evolution of Saturn. Each graph shows a 4-Gyr numerical trajectory (black dots) computed for Titan’s nominal migration rate and for a given value of Saturn’s normalised polar moment of inertia specifiedin title. Today’s location of Saturn is represented by the big green spot; the vertical error bar corresponds to our full exploration interval of λ. The red curves show the centre of first-order secular spin-orbit resonances (Cassini state 2) and the coloured areas represent their widths (same as Fig. 1). The top large area is the resonance with ϕ3 and the bottom thin area is the resonance with ϕ19 (see Table 1). The separatrices of the ϕ3 resonance are highlighted in blue. Going forwards in time, the trajectories go from bottom to top.

In the text

Fig. 5 Current dynamics of Saturn’s spin axis according to its normalised polar moment of inertia λ. The value of λ (top horizontal axis) is linked to the current precession constant of Saturn (bottom horizontal axis) through Eqs. (3) and (5). The black interval shows the ‘instantaneous’ oscillation range of Saturn’s spin axis (i.e. without drift of α) obtained by numerical integration. The resonant angle is σ3 = ψ + ϕ3 (see Sect. 2). The green line shows Saturn’s current obliquity and resonant angle. The background colour indicates the type of past evolution as labelled in the top panel (see text for the numbering).

In the text

	Fig. 6 Period of small oscillations about the resonance centre for a resonance with ϕ3 or ϕ19. The resonant angles are σ3 = ψ + ϕ3 and σ19 = ψ + ϕ19, respectively.Dashed curves are used for oscillations about Cassini state 2 before the separatrix appears. The appearance of the separatrix is marked by a blue dot.
In the text

Fig. 7 Past obliquity of Saturn for different migration rates of Titan. The top and bottom horizontal axes are the same as in Fig. 5 and the horizontal green line shows Saturn’s current obliquity. For a given value of the normalised polar moment of inertia λ (top horizontal axis), the curve width shows the oscillation range of obliquity 4 Gyr in the past obtained by backward numerical integration. The migration rates are labelled on each panel as a fraction of the nominal rate of Lainey et al. (2020). The four coloured curves correspond to the migration rates illustrated in Fig. 3. The grey stripes in the central panel highlight trajectories of Type 2 (same as in Fig. 5). The value of b in the top left panel corresponds to a current quality factor Q equal to 5000 (see Lainey et al. 2020).

In the text

Fig. 8 Past obliquity of Saturn as a function of Titan’s migration velocity and Saturn’s polar moment of inertia. Each panel of Fig. 7 corresponds here to a vertical slice. The colour scale depicts the minimum obliquity of the oscillation range, and the white curve highlights the 5° level. The 3σ uncertainty ranges of Lainey et al. (2020) yield today approximately b∕b0 ∈ [1∕2, 5] for the astrometric measurements and b∕b0 ∈ [1, 3∕2] for the radio-science experiments (see Fig. 3).

In the text

Fig. 9 Zoom-in view of the central panel of Fig. 7. We use a red curve to highlight the bottom limit of the blue interval, otherwise the narrowness of the spikes makes them invisible (the width of spike d is Δλ ≈ 10−14). This graph can be compared to Fig. 3 of Saillenfest et al. (2021), where such level of fine tuning is not shown due to its questionable physical relevance. See Appendix B for examples of trajectories.

In the text

Fig. 10 Snapshots of a Monte Carlo experiment computed for λ = 0.204 and Titan’s nominal migration rate. This experiment features 101 values of initial obliquity between 0° and 60°, for which 240 values of initial precession angle are regularly sampled in [0, 2π). The value of α reached by the trajectories at the time of the snapshot is labelled on each panel. Each trajectory is represented by a small dot which is coloured according to the variation range of the resonant angle σ3 (obtained bya 0.5-Gyr numerical integration with constant α). The horizontal green line shows the current obliquity of Saturn. At the beginning of the propagations, all trajectories are coloured red (since σ3 circulates), and distributed along a diagonal line. Then, as α increases over time, trajectories are dispersed off the diagonal according to the four types of trajectories depicted in Fig. 4 and labelled in the penultimate panel.

In the text

Fig. 11 Capture probability of Saturn in secular spin-orbit resonance with ϕ3 as a functionof its primordial obliquity. For each initial obliquity (401 values between 0 and 60°), 720 values of initial precession angle are uniformly sampled in [0, 2π) and propagated forwards in time starting from − 4 Gyr and until every trajectory has reached the resonance. This experiment is repeated with two different migration laws for Titan (see labels). The result is virtually independent of the value chosen for λ. The fractionof captured trajectories (coloured curves) is compared to the perfect adiabatic case (black curve) computed with the analytical formulas of Henrard & Murigande (1987).

In the text

Fig. 12 Brute-force search for Saturn’s past obliquity using the loose success criteria: probing all dynamical pathways. We use Titan’s nominal migration law (b = b0). Among all trajectories evenly sampled in the space of the normalised polar moment of inertia λ (top horizontal axis, 101 values), of the initial obliquity (vertical axis, 101 values), and of the initial precession angle (240 values between 0 and 2π), we only keep those matching Saturn’s spin axis today according to our loose success criteria (see text). Each point is coloured according to the number of successful runs among the 240 initial precession angles; the success ratio is written below the colour bar. A point is not drawn if no successful trajectory is found. In the back, the blue interval shows the past obliquity of Saturn obtained by backward numerical integration (same as Fig. 7 for b = b0), showing theconsistency between backward and forward integrations in time. The background stripes and their labels have the same meaning as in Fig. 5.

In the text

	Fig. 13 Same as Fig. 12, but using the strict success criteria: comparing the relative likelihood of producing Saturn’s current state. As in Fig. 12, each point of the graph is made of 240 simulations with initial ψ ∈ [0, 2π). The red rectangle highlights the region featuring the highest success ratios.
In the text

Fig. 14 Distribution of the solutions starting from a low primordial obliquity and matching our strict success criteria. For each set (b, λ) of the parameters, 2400 values of initial obliquity ε and precession angle ψ are drawn from a uniform random distribution in (ε, ψ) ∈ [0°, 5°] × [0, 2π). Coloured dots show the parameter sets (b, λ) for which at least one successful trajectory was found; the success ratio is written below the colour bar. Light-grey crosses mean that no successful trajectory was found over our 2400 initial conditions. The black contours show the 5°-level obtained through backward numerical integrations (same as Fig. 8), showing the consistency between backward and forward integrations in time. The black lines in the top portion show the approximate location of the border of the blue stripes in Fig. 8, where extreme phase effects can happen; the corresponding ranges of parameters are so narrow that they are missed by the resolution of Fig. 8 (see Sect. 3.3).

In the text

	Fig. 15 Same as Fig. 14, but for statistics based on a sub-sample of simulations. (a) Initial conditions in (ε, ψ) ∈ [0°, 2.5°] × [0, 2π). (b) Initial conditions in (ε, ψ) ∈ [2.5°, 5°] × [0, 2π). In both panels, each point is made of about 1200 initial conditions extracted from the simulations from Fig. 14.
In the text

	Fig. 16 Future obliquity of Saturn as a function of Titan’s migration velocity and Saturn’s polar moment of inertia. The axes are the same as in Fig. 8. The 3σ uncertainty ranges of Lainey et al. (2020) yield approximately b∕b0 ∈ [1∕2, 5] for the astrometric measurements and b∕b0 ∈ [1, 3∕2] for the radio-science experiments. Some level curves are shown in red.
In the text

	Fig. 17 Example of Type 2 trajectory that is expelled out of the resonance. It has been obtained for λ = 0.221076 and a migration rate b∕b0 = 3. The integration runs from today up to 5 Gyr in the future.
In the text

Fig. B.1 Example of trajectories featuring an extreme phase effect. Left column: evolution of the obliquity, and right column: evolution of the resonant angle σ3 = ψ + ϕ3. The migration parameter is b = b0. For each row, the parameter λ used corresponds to the tip of a spike in Fig. 9 (see labels), tuned at the 10−15 level. The pink area represents the interval occupied by the resonance once the separatrix appears. The blue curve shows the locationof the hyperbolic equilibrium point (Cassini state 4). The green point shows Saturn current location (at t = 0).

In the text

	Fig. B.2 Same as Fig. B.1, but for a trajectory of Type 1. This trajectory has a migration parameter b = 7.37 b0 and a normalised polar moment of inertia λ = 0.2114.
In the text

	Fig. C.1 Same as Fig. 14, but considering an isotropic distribution of initial spin orientation with ε ≤ 5°. It is obtained from Fig. 14 by weighting the count of each run by sin ε (see Fig. C.2).
In the text

	Fig. C.2 Sampled distribution of initial obliquity for one arbitrary point of Fig. 14, made of 2400 simulations. The raw count of sampled trajectories is shown in red; the weighted count is shown in blue. Top: histogram with respectto the obliquity. Bottom: histogram with respect to the cosine of obliquity. The probability density functions (‘pdf’) are shown by the red and blue curves.
In the text

Fig. C.3 Same as Fig. 14, but for a range of initial spin orientations enlarged to 0° ≤ ε ≤ 10°. Each point is made of 2400 numerical simulations. (a) Uniform random distribution of (ε, ψ) ∈ [0°, 10°] × [0, 2π). (b) Uniform random distribution of (ε, ψ) over the spherical cap defined by ε ≤ 10°. As in Fig. C.1, panel b is obtained by weighting the count numbers of panel a. The black contours show the 5° and 10° levels obtained through backward numerical integrations (see Fig. 8).

In the text

	Fig. C.4 Same as Fig. 14, but for initial obliquities uniformly distributed between 2.5° and 7.5°. It is obtained from a sub-sample of Fig. C.3a, such that each point is made of about 1200 runs. The black contours show the 5° and 10° levels obtained through backward numerical integrations (same as Fig. C.3).
In the text