The past and future obliquity of Saturn as Titan migrates

Aims: It has recently been shown that the fast migration of Titan could be responsible for the current 26.7{\deg}-tilt of Saturn's spin axis. We aim to quantify the level of generality of this result and to measure the most likely sets of parameters. We also aim to determine the obliquity that Saturn will reach in the future. Methods: We explore a broad range of parameters and propagate numerically the orientation of Saturn's spin axis both backward and forward in time. Results: In the adiabatic regime, the likelihood of reproducing Saturn's current spin-axis orientation is maximised for primordial obliquities between about 2{\deg} and 7{\deg}. For a slightly faster migration than expected from radio-science experiments, non-adiabatic effects even allow for exactly null primordial obliquities. Starting from such small tilts, Saturn's spin axis can evolve up to its current state provided that: i) the semi-major axis of Titan changed by more than 5% of its current value, and ii) its migration rate does not exceed 10 times the nominal measured rate. In comparison, observational data suggest that the increase of Titan's semi-major axis exceeded 50% over 4 Gyrs, and error bars imply that the current migration rate is unlikely to be larger than 1.5 times its nominal value. Conclusions: If Titan did migrate substantially before today, tilting Saturn from a small obliquity is not only possible, but it is the most likely scenario. Saturn's obliquity is expected to be still increasing today and could exceed 65{\deg} in the future. Maximising the likelihood would also put strict constraints on Saturn's polar moment of inertia. However, the possibility remains that Saturn's primordial obliquity was already large, for instance as a result of a massive collision. The unambiguous distinction between these two scenarios would be given by a precise measure of Saturn's polar moment of inertia.


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 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 • .  have shown 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 Boué et al. 2009;Vokrouhlický & Nesvorný 2015;Brasser & Lee 2015). However, Saillenfest et al. (2021) recently showed 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;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 pre-cession 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 of 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 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 migra-Article number, page 1 of 24 arXiv:2101.10366v1 [astro-ph.EP] 25 Jan 2021 A&A proofs: manuscript no. saturnspin tion (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 of following 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 . 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 experiments to 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.

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  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 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 eccen-tricity e of the planet and through the functions and C(t) = qṗ − pq . (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: where G 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 2 is its second zonal gravity coefficient, and λ is its normalised polar moment of inertia. The parameters J 2 and λ can be expressed as 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 closein 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 2 and λ in Eq. (3) by the effective values: 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 J 2 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 though its Laplace plane is not much inclined, taking its inclination into account changes the whole satellites' contribution by several percent 1 . 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 6 ≈ 0.62 • ) and Iapetus (L 8 ≈ 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ψ = αX/(1 − e 2 ) 3/2 , 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 to Hamilton's equations of motion (see e.g. 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 6 ε). We can therefore use the same trick as Saillenfest et al. (2021) and replace Eq. (5) bỹ where only Titan is considered (k = 6), in the close-in regime (L 6 = 0), and where its mass m 6 has been slightly increased (m 6 ≈ 1.04 m 6 ) so as to produce the exact same value of α today using Eq. (6) instead of Eq. (5). This slight increase of 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

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: 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 over time t with frequencies µ k and ν k : 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 1999or 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). Notes. Due to the secular resonance (g 1 − g 5 ) − (s 1 − s 2 ), an additional fundamental frequency γ appears in term 20 (see Laskar 1990). (*) There are typographical errors in Laskar (1990) in the identification of the 14th and 16th terms.
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 −1 to 2 ·yr −1 . Because of the chaotic dynamics of the Solar System (Laskar 1989), the fundamental frequencies related to the terrestrial planets (e.g. s 1 , s 2 , s 4 , 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 Gyrs 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 Gyrs 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.

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 obser-Article number, page 3 of 24 A&A proofs: manuscript no. saturnspin 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 interval of 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.
vations, 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 2 ). 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 any decisive 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) obtains a value of 0.198. This value is smaller than other estimates found in the literature, even though it is marginally compatible with the calculations of Hubbard & Marley (1989), who give λ = 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 . 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 2 . The computations of Nettelmann et al. (2013) yield 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 by 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 2 from Iess et al. (2019). For consistency with Saturn's orbital solution (Sect. 2.2), we take its mass and secular semimajor axis from Laskar (1990).
In order to compute J 2 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 3 . 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) ism 6 ≈ 1.04 m 6 . Similar results are obtained when using the more precise values of L k given by the constant terms of the full series of Vienne & Duriez (1995) and Duriez & Vienne (1997).
Because of 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 of α. 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 where a 0 is Titan's current mean semi-major axis, t 0 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 4 . Equation (9) implies that Titan's current tidal timescale t tide = a 6 /ȧ 6 relates to b as b = t 0 /t tide . Considering a Saillenfest et al.: The past and future obliquity of Saturn as Titan migrates   For the long-term evolution of Saturn's satellites, they adopt a nominal value of b 0 = 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 0 . 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 massm 6 yet allows us to obtain the right value of today's precession constant α, as if all satellites were included.

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. With respect 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).

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 . 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 φ 3 in Table 1 and largest resonance in Fig. 1). The four possible kinds of past evolution are illustrated in Fig. 4 for b = b 0 . They are namely: 1. For λ 0.220, Saturn went past the resonance through its hyperbolic point. 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. 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. 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 σ 3 and allowing for the smallest past obliquity of Saturn. Type 4 is ruled out by our uncertainty range of λ.
During its past evolution, Saturn also crossed a first-order secular spin-orbit resonance with the term φ 19 which has frequency g 7 − g 8 + s 7 (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 σ 19 = ψ + φ 19 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 Article number, page 5 of 24 A&A proofs: manuscript no. saturnspin for Titan's nominal migration rate and for a given value of Saturn's normalised polar moment of inertia λ = C/(MR 2 eq ) specified in 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 further highlighted in blue. Going forward in time, the trajectories go from bottom to top. timescales of σ 3 and σ 19 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.

Adiabaticity of Titan's migration
If the drift of α over time was perfectly adiabatic (i.e. infinitely slow compared to the oscillations of σ 3 ), the outcome of 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) showed 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 Gyrs 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 the blue 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 σ 3 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 0 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 very first 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 0 , 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 0 , 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 0 , 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  (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). equal to exactly 0 • . Such a null value is obtained when the oscillation phase of σ 3 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 σ 3 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 φ 19 . In Fig. 7, the effect of overly fast crossings is clearly visible for b 11b 0 . Beyond this approximate limit, all trajectories in our backward integrations cross the resonance separatrix,  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.
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. Nonadiabaticity is revealed by the coloured waves, denoting phase effects. As expected, the waves disappear for b b 0 : this is the adiabatic regime (see Fig. 4 of Saillenfest et al. 2021 for a zoomin view). For very small migration rates, however, Titan did not have time in 4 Gyrs to migrate enough to produce substantial effects on Saturn's obliquity. This is why the dark-blue region of 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 0 . This migration rate can therefore be considered as the largest one allowing Titan to be held responsible for Saturn's current large obliquity.

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 deepest spikes denote the strongest phase effects; they correspond to trajectories that reached the resonance almost exactly at its hyperbolic equilibrium point (called Cassini state 4: see e.g. Saillenfest et al. 2019 for phase portraits 5 ). 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 5 There is a typographical error in Saillenfest et al. (2019): the list of the Cassini states given before Eq. (22) should read (4,2,3,1) instead of (1,2,3,4) in order to match the denomination introduced by Peale (1969).
Article number, page 7 of 24 A&A proofs: manuscript no. saturnspin 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 0 ), 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 0 (central panel), a tuning of λ at the 10 −15 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 −15 ) 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 they relied on a different tilting scenario. For a nonadiabatic 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 −8 . 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   Fig. 7. We use a red curve to highlight the bottom limit of the blue interval, otherwise the narrowness of the spikes make 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.
vanishingly small probability of producing such trajectories is confirmed in Sect. 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 φ 3 . 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.

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 Gyrs 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 0 ). The first snapshot is taken about 20 million years after the start of the integrations, whereas the last snapshot is taken at today's epoch. See Fig. 4 for one example of trajectory. Changing the value of λ produces a shift of Saturn's precession constant α but no strong variation of its drift rate (see Fig. 2). Moreover, since this drift is almost an adiabatic process for b = b 0 (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 timet for another valuẽ λ. 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 timet = f −1 λ (α). 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 φ 3 , 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;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. For growing initial obliquity, the probability of Type 2 trajectories (i.e. capture) decreases, favouring Type 1 trajectories instead (i.e. crossing without capture).

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 spinaxis oscillation interval, and ii) the libration amplitude of the resonant angle σ 3 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 σ 3 , 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 of 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 with the initial obliquity value (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 σ 3 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 can 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;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(Henrard , 1993. This probabilistic modelling of chaos explains why separatrix crossings are not reversible when using this theory.   Fig. 11. Capture probability of Saturn in secular spin-orbit resonance with φ 3 as a function of 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 forward in time starting from −4 Gyrs 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 fraction of captured trajectories (coloured curves) is compared to the perfect adiabatic case (black curve) computed with the analytical formulas of Henrard & Murigande (1987).

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 of 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 therefore must be considered, 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 σ 3 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 lies for past obliquities between about 2 • and 7 • , and current precession constant α between about 0.76 and 0.79 ·yr −1 (red box).
As already discussed by  and , 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 Type 1 and 2 both have a non-zero probability of failure, see Fig. 11); second, Type 3 trajectories oscillate today with a small amplitude Article number, page 11 of 24 A&A proofs: manuscript no. saturnspin inside the resonance, which means that all of them feature a similar precession angle ψ imposed by the resonance relation σ 3 ∼ 0. On the contrary, other types of trajectories either feature a large oscillation amplitude of σ 3 (Type 2) or circulation of σ 3 (Types 1 and 4); therefore, they only sweep over Saturn's actual orientation once in a while, and matching it today is only a low-probability 'coincidental' event 6 . As shown by Fig. 13, the least favoured trajectories are those of Type 2, especially for high initial obliquities, because of the strong decrease of capture probability (see Fig. 11).
In order to explore all migration rates and bring further constraints on the model parameters, we now turn to a 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 Gyrs 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 6 The same argument has been pointed out for Jupiter by Ward & Canup (2006) and Saillenfest et al. (2020). 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 7 (see Table 1), but contrary to Saturn, its satellites did not migrate enough yet to substantially increase its obliquity ; 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 nonuniform distribution of spin-axis directions over the unit sphere that favours small obliquities 7 . 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 Gyrs 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 0 b 5b 0 resemble to some extent the suc- 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  cessful matches of Vokrouhlický & Nesvorný (2015), reminding us that the basic dynamical ingredients are the same even though the mechanism of resonance encounter 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 0 and 6b 0 . 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 0 ∈ [1/2, 5]), but not with the uncertainty ranges given by radio-science experiments (b/b 0 ∈ [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 0 ) 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 0 ), 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 φ 3 ); and Brasser & Lee (2015) assumed ε ≈ 0.05 • (with respect to the in-variable 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 those studies 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.

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 . In order to explore the future obliquity dynamics of Saturn, we propagate Saturn's spin-axis from today up to 5 Gyrs 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 0 because for larger migration rates, Titan goes beyond the Laplace radius during the integration timespan (a 6 ≈ 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 φ 3 (see Fig. 4) and they go farther away from it as α continues to increase. The increase of α makes them cross the resonances with φ 51 , with φ 14 , and with φ 15 (see Fig. 1 and Table 1). Being very small, these resonances are crossed quickly and they do not 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 0 (the fastest migration presented in Fig. 16), trajectories reach the lower fringe of the strong resonance with φ 4 at the end of the integration, but they do not actually reach it.
The middle portion of Fig. 16 features trajectories of Types 1 and 2. Such trajectories are currently inside the resonance with φ 3 (see Fig. 4) and they follow its centre as α increases. After 5 Gyrs 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 of α 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 in Fig. 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 to increase. The resonance encounter can either lead to a capture (like Type 2 trajectories) or to a permanent decrease of obliquity (like Type 1 trajectories). The outcome is determined by the phase of σ 3 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.

Discussion and conclusion
Since 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 have further investigated the longterm 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 ). As Titan migrates over time, it produces a drift of 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 small to large values 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 at least 1 Saturn's radius after the late planetary migration, about 4 Gyrs 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 of Saturn's primordial obliquity. Zero values are however disfavoured by the error range of radioscience experiments, which yield most likely primordial obliquities values between 2 • and 7 • .
Our Monte Carlo experiments do not reveal a strong chaotic mixing of trajectories, even though borderline separatrixcrossing 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 the highest likelihood of reproducing Saturn's current spin-axis orientation, surpassing high-obliquity alternatives by a factor of about 10. 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 have dramatic effects on Saturn's obliquity provided that it 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 Gyrs 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 of migration rates 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 it would not have affected Saturn's past obliquity dynamics. 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. As shown above, 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. Modeldependent theoretical values are not enough for this purpose, and it is still unclear what are 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. 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 of 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. 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 σ 3 = ψ+φ 3 , where φ 3 evolves as in Eq. (8). The elliptic point of the resonance (Cassini state 2) is located at σ 3 = 0, and the hyperbolic equilibrium point (Cassini state 4) is located at σ 3 = π (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 infinitely close 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 σ 3 . The two other isolated grey points in Fig. 14 have performed 0 and 2, respectively.

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 depends 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 Notes. This solution has been directly obtained from Laskar (1990) as explained in the text. The phases θ (0) k are given at time J2000.
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 simulation again by simply weight-ing 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 2 the resolution at the highobliquity 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 0 ) and non-adiabatic (b 3b 0 ) 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).