© M. Saillenfest and G. Lari 2021

1 Introduction

A secular spin–orbit resonance occurs when the spin-axis precession rate of a planet becomes commensurate with a frequency (or a combination of frequencies) that appears in its orbital precession. Secular spin–orbit resonances were first studied individually and linked to Cassini’s laws (Colombo 1966; Peale 1969; Ward 1975; Henrard & Murigande 1987). The effect on the spin-axis dynamics of a whole multi-harmonic orbital precession spectrum has also been investigated, and the overlap of several secular spin–orbit resonances has been identified as responsible for large chaotic regions in the inner Solar System (Ward 1973, 1982; Laskar & Robutel 1993; Néron de Surgy & Laskar 1997; Laskar et al. 2004). More recently, higher-order resonances have been characterised in a systematic way, and their relation to the emergence of chaos has been assessed (Li & Batygin 2014; Saillenfest et al. 2019b). In fact, secular spin–orbit resonances are found to rule the long-term spin-axis dynamics of planets not only in the Solar System but also in extrasolar systems (see e.g. Atobe et al. 2004; Deitrick et al. 2018; Shan & Li 2018; Millholland & Laughlin 2019).

As shown by Ward & Hamilton (2004), Saturn is today very close to or inside a secular spin–orbit resonance with the nodal orbital precession mode of Neptune, noted s₈. The current large 26.7° obliquity of Saturn probably results from this resonance (Hamilton & Ward 2004). It was first thought that the resonance trapping occurred more than four billion years ago during the late planetary migration (Boué et al. 2009; Brasser & Lee 2015; Vokrouhlický & Nesvorný 2015). However, this would require Saturn’s satellites to not have migrated much since this event, which contradicts the fast migration of the satellites measured by Lainey et al. (2020). Instead, Saillenfest et al. (2021a) have shown that the migration of Saturn’s satellites, and in particular of Titan, is likely responsible for the resonance encounter. The resonant interaction therefore began more recently than previously thought, perhaps about one billion years ago. Using Monte Carlo simulations, Saillenfest et al. (2021b) show that in order to reproduce Saturn’s current state, the most likely dynamical pathway is a gradual tilting starting from a few degrees before the resonance encounter. Since a near-zero primordial obliquity is also what is expected from planetary formation theories (see e.g. Ward & Hamilton 2004, Rogoszinski & Hamilton 2020, and references therein), and even though primordial non-zero obliquities are not totally excluded (Millholland & Batygin 2019; Martin & Armitage 2021), this scenario appears quite promising.

If Saturn did follow its expected pathway (and was not affected by an accidental major impact; see e.g. Li & Lai 2020), then Saturn should still be trapped inside the resonance today. As Titan continues migrating, Saturn should therefore continue to follow the drift of the resonance centre in the future. Because of this mechanism, Saturn may reach very large obliquity values in the next few billions of years (Saillenfest et al. 2021b).

Yet, the final outcome of this dynamical mechanism remains unknown. The obliquity of a planet cannot increase forever, and there must exist some kind of dynamical barrier, either on the planet’s or on the satellite’s side, that would halt the tilting at some point. Even though this final outcome may not be directly relevant for Saturn and Titan because of the large timescales involved (see below), its generic nature makes it important even from the point of view of pure celestial mechanics as other planets and exoplanets may have been affected. Hence, we aim to characterise the full tilting mechanism in a general way, without anyassumption about the distance of the satellite, and with a special focus on the case of Saturn and Titan.

Since Titan is still far from its Laplace radius today and may only reach it after billions of years of continuous orbital expansion (if not tens of billions of years), previous studies have been restricted to the close-in satellite regime, in which Titan’s orbit lies in Saturn’s equatorial plane (see e.g. Goldreich 1966). Out of this regime, regular satellites are known to oscillate around their local Laplace plane, which is inclined halfway between the equator and the orbital plane of the host planet (Laplace 1805). Most importantly, Tremaine et al. (2009) found that the satellite is unstable in the vicinity of its Laplace radius if the planet’s obliquity is larger than about 69° (for a circular satellite) or 71° (for an eccentric satellite). According to Saillenfest et al. (2021b), Saturn’s obliquity may exceed these thresholds in a few gigayears from now, depending on Titan’s migration rate. Hence, if Titan happens to be located near its Laplace radius at this stage of the evolution, a non-trivial dynamics is expected for the system. The simulations of Tremaine et al. (2009) and Tamayo et al. (2013) revealed that, in some ranges of parameters, strong chaotic transitions in the satellite’s eccentricity and inclination are possible.Tamayo et al. (2013) drew a parallel between the eccentricity increase of the satellite and the ZKL mechanism (for ‘von Zeipel–Lidov–Kozai’; see Ito & Ohtsuka 2019), in which large eccentricity oscillations occur while the satellite’s argument of pericentre oscillates around a fixed value. Beyond some eccentricity threshold, the effect of planetary oblateness re-initiates the apsidal precession, with the result of averaging to zero the solar torque and stopping the eccentricity increase¹. More recently, Speedie & Zanazzi (2020) performed an extensive numerical exploration of the stability of particles initially located near their local Laplace plane. Their study confirms the secular instabilities reported by Tremaine et al. (2009) and Tamayo et al. (2013), and their fully unaveraged model allows for other kinds of instability to appear, driven by the evection and ‘ivection’ resonances².

These previous results show that studying Saturn’s tilting mechanism in a general way requires one to keep an eye on both the satellite’s and planet’s dynamics. Out of the close-in satellite regime, and a fortiori if the satellite becomes unstable, the model used by Saillenfest et al. (2021a,b) for Saturn’s spin-axis dynamics is invalidated. As a first step before developing a complete numerical model, our goal in this article is to establish a qualitative understanding of what happens to the planet and its satellite when the secular spin–orbit resonance leads them to their ultimate large-obliquity regime, where previous models fail.

In Sect. 2, we revisit the concept of Laplace surface introduced by Tremaine et al. (2009); we go further in the analytical characterisation of the equilibria and focus on the large-obliquity regime. In Sect. 3, we describe the influence of the satellite’s dynamics on the spin-axis motion of the planet. We provide simplified formulas that allow the planet’s obliquity evolution to be described as a function of the satellite’s properties. In Sect. 4, we apply our findings to Saturn and Titan and explore their coupled dynamics as Titan migrates outwards. We conclude in Sect. 5 and present some further applications of our results in other contexts.

2 Orbital motion of the satellite

In order to investigate the way satellites interact with the spin axis of their host planet, we must get a clear understanding of their orbital dynamics as well. In this section, we first consider a massless satellite orbiting an oblate planet, which has itself a fixed orbit around the star (or an orbit that can be regarded as fixed over the interval of time considered).

2.1 Equations of motion

The Hamiltonian function describing the orbital motion of the massless satellite around the planet can be written , where is the Keplerian part and gathers the orbital perturbations. The parameter ϵ ≪ 1 stresses that the orbital perturbations are small; neglecting , the long-term behaviour of the satellite is described by the secular Hamiltonian , which can be written (1)

where comes from the planet’s oblateness and comes from the star’s gravitational attraction. The secular semi-major axis a of the satellite is a constant of motion and a parameter of . Considering that a is much larger than the planet’s equatorial radius R_eq and much smaller than the star’s semi-major axis a_⊙ in its orbit around the planet, both terms of Eq. (1) can be expanded in Legendre polynomials. As Tremaine et al. (2009), we first limit the expansion to the quadrupolar approximation, which amounts to neglecting and . This leads us to define the two fixed parameters of Eq. (1) as (2)

and the two parts of the Hamiltonian function as (3)

In these expressions, μ_P and J₂ are the gravitational parameter and the second zonal gravity coefficient of the host planet, and μ_⊙, a_⊙, and e_⊙ are the gravitational parameter, the semi-major axis, and the eccentricity of the star, respectively. The usual orbital elements of the satellite are written (e, I, ω, Ω), and we use the index Q for quantities measured with respect to the planet’s equator, and the index C for quantities measured with respectto the planet’s orbital plane (that we improperly call the ‘ecliptic’). Since Eq. (3) is obtained through a multi-polar development only, it is valid for arbitrary eccentricity and inclination of the planet and of its satellite (Laskar & Boué 2010).

The Hamiltonian function in Eq. (3) has been averaged over the star’s orbital motion. As stressed by Tremaine et al. (2009), this approximation requires not only that a ≪ a_⊙, but even that a ≪ r_H, where r_H is the Hill radius of the planet. In particular, Eq. (3) does not contain the evection resonance, which can have important effects for far-away satellites (see e.g. Frouard et al. 2010; Speedie & Zanazzi 2020). In Sect. 4, we verify that Eq. (3) provides a fair approximation of Titan’s orbital dynamics.

The coordinates of the satellite measured with respect to the equator (Q) or to the ecliptic (C) are linked through the obliquity ε of the planet via the relations given in Appendix A. These relations can be used to express in terms of the equator or ecliptic coordinates only. Noting C = cosε and S = sinε, the Hamiltonian in Eq. (3) can be written equivalently as (4)

where δ_C ≡ Ω_C − Ω_P and Ω_P is the ascending node of the planet’s equator measured along the ecliptic. Likewise, the Hamiltonian in Eq. (3) can be written equivalently as (5)

where δ_Q ≡ Ω_Q − Ω_⊙ and Ω_⊙ is the ascending node of the star measured along the equator of the planet.

If the planet’s axis of figure has a fixed orientation in space, then Ω_P is a constant angle, and both the ecliptic and equatorial reference frames are inertial. This is equivalent to considering that the spin-axis precession of the planet is infinitely slow compared to the timescales relevant for the satellite. The validity of this hypothesis will be discussed in Sect. 3. For now, we consider that Ω_P is constant and examine the dynamical system described by Eq. (3), expanding on the work of Tremaine et al. (2009).

First of all, the parameters k_P and k_⊙ in Eq. (2) make appear a characteristic length called ‘Laplace radius’ defined by (6)

We also introduce a critical radius r_M, already used by Goldreich (1966), that we define by (7)

As noticed by Tamayo et al. (2013), it is more natural to use r_M as a reference radius than the conventional r_L of Tremaine et al. (2009). The symbol M stands here for ‘midpoint’ and the dynamical meaning of r_M will appear clear below³. Using this definition, we can rewrite as (8)

where a change of timescale could be used to remove the leading constant factor. In order to investigate the dynamics of a slowly migrating satellite, it is more convenient to introduce a timescale that does not involve its semi-major axis. As shown below, the frequency κ, defined as (9)

naturally appears in the dynamics, and it is therefore a good choice of characteristic timescale. We define the corresponding period as τ = 2π∕κ. Hence, we can describe the full variety of trajectories of the satellite by the only two parameters a∕r_M and ε, and their evolution timescale is provided by the period τ. Table 1 lists these parameters for various satellites in the Solar System. We also include the case of distant trans-Neptunian objects perturbed by the galactic tides, which have an almost identical dynamics (see Saillenfest et al. 2019a).

In order to study the orbital dynamics of the satellite, it is more suitable to use a set of coordinates that are not singular for circular and/or zero-inclination orbits, as the usual rectangular coordinates (10)

Alternatively, one can use a vectorial formulation as Tremaine et al. (2009).

2.2 The Laplace states

By writing down the equations of motion in a non-singular set of coordinates, we see that the condition e = 0 is an equilibrium point for the satellite whatever its other orbital elements. Moreover, linear stability analysis shows that the eccentricity and inclination degrees of freedom are decoupled in the vicinity of e = 0. Assuming that the satellite’s eccentricity is zero (for instance, if it has been damped at the time of its formation in a circumplanetary disc), we can therefore study the evolution of its inclination degree of freedom in a decoupled way.

Figure 1 shows examples of trajectories for the satellites obtained by plotting the level curves of for e = 0. The dynamics of the satellite is described by the direction of its orbital angular momentum; since the dynamics actually lie on a sphere, any planar representation of the trajectories has coordinate singularities. In Fig. 2a, we show the same phase portrait as Fig. 1 on the sphere. The system being secular, it is independent of whether the orbits and spins are prograde or retrograde. This is traduced by the invariance of the phase space to the transformations (δ_Q, I_Q) → (π + δ_Q, π − I_Q) and (ε, I_Q) → (π − ε, π − I_Q). Three kinds of equilibrium points can be seen, which we label P₁, P₂, and P₃. In the work of Tremaine et al. (2009), the points P₁ and P₃ are called ‘circular coplanar Laplace equilibria’, and the point P₂ is called ‘circular orthogonal Laplace equilibrium’. This denomination clearly reflects the geometry of these configurations. For the sake of succinctness, we call them ‘Laplace states’ 1, 2, and 3 (in reference to the famous ‘Cassini states’ described below).

The geometry of the phase portraits for any value of the parameters can be described by the location of the equilibrium points and the shape of the separatrix. The respective locations of the Laplace states when varying the parameters are illustrated in Figs. 3 and 4. Because of the symmetries mentioned above, each equilibrium point has a twin obtained by the transformation (δ_Q, I_Q) → (π + δ_Q, π − I_Q) that corresponds to the same Laplace state with reversed orbital motion.

In the space of parameters, there is a critical point, that we call S₁, defined by (11)

At point S₁, the Laplace states P₁ and P₃ and the separatrix degenerate into an equilibrium circle spanning allvalues of inclination (see Fig. 2b). All points of this circle are stable equilibrium configurations in which the linearised problem has zero eigenfrequency for inclination variations (we note it ). As shown by Figs. 3 and 4, going through point S₁ by smoothly changing the parameters inverts the locations of P₁ and P₃. We stress that in Fig. 4 the apparent jumps of P₁ (for a < r_M) and P₃ (for a > r_M) are only coordinate singularities in which P₁ or P₃ smoothly pass through the pole of the sphere (see Fig. 2). On the contrary, the jump observed at point S₁ (a = r_M) is a real singularity.

Another singularity occurs for a null or 180° obliquity. We call S₂ the corresponding region of the parameter space, defined by (12)

In region S₂, the Laplace states P₂ and P₃ and the separatrix degenerate into an equilibrium circle spanning all values of δ_Q (see Fig. 2c). All points of this circle are stable equilibrium configurations in which the linearised problem has zero eigenfrequency for inclination variations (we note it ). The regions S₁ and S₂ of the parameter space can be visualised in Fig. 5.

Apart from regions S₁ (in which P₁ becomes singular) and S₂ (in which P₂ becomes singular), the phase space keeps the same topology whatever the parameters a∕r_M > 0 and ε. This means that the Laplace states are smoothly transported by a continuous change of parameters, and they keep their stability nature against inclination variations. On Fig. 2a, such a continuous change of parameter would simply produce the rotation of the sphere around the x-axis and the narrowing or widening of the black separatrix. More precisely, Fig. 4 shows that for a < r_M, varying the obliquity produces an oscillation of P₁ around the pole and P₃ remains near 90°; for a > r_M, on the contrary, varying the obliquity makes P₁ and P₃ roll all over the sphere. The opposite behaviour would be obtained by representing the ecliptic inclination I_C instead of I_Q. The location and stability nature of the Laplace states play a fundamental role in the combined dynamics of the satellite’s orbit and the planet’s spin axis. For this reason, we review here their basic properties and go deeper than previous works in their analytic characterisation.

P₁ is stable to inclination variations. As illustrated by Fig. 3, it corresponds to an orbit lying on the equator plane for close-in satellites (a ≪ r_M), and on the ecliptic plane for far-away satellites (a ≫ r_M). In between, P₁ corresponds to an intermediate tilt between the equator and the ecliptic. As a result of eccentricity and inclination damping, P₁ is expected to be the birth place of regular satellites formed in a circumplanetary disc. P₁ is therefore particularly important in satellite dynamics studies; for this reason, it is called ‘classical Laplace equilibrium’ by Tremaine et al. (2009). For δ_Q = 0 (dark blue colour in the figures), the inclination of P₁ is given by one of the two solutions of the equation⁴ (13)

the second solution being the inclination of P₃. We note them I_Q1 and I_Q3. Their closed form expressions can be written (14)

where . At a = r_M, the Laplace state P₁ lies exactly halfway between the equator and the ecliptic planes (i.e. I_Q1 = ε∕2 for ε < 90°). This is why we use the index M, for ‘midpoint’, introduced in Eq. (7). Interestingly, this midpoint does not depend on the value of the obliquity ε, but only onthe distance of the satellite. The curve described by Eq. (13) and illustrated in Fig. 3, however, is not exactly symmetric with respect to a = r_M. Its inflexion point F (for ‘flex’) is reached at radius r_F(ε), defined by (15)

and illustrated in Fig. 6. The distance between r_F and r_M is a way to quantify the asymmetry of I_Q1 as a function of a. We stress that all level curves in Fig. 6 converge at S₁. Through a smooth variation of parameters, the satellite can therefore reach the singular point S₁ from any orbital inclination between 0° and 180°. This property has important consequences for the spin-axis dynamics of the host planet, as we discuss in Sect. 3. We also show that r_F divides the close-in and far-away satellite regimes considered in previous works.

Tremaine et al. (2009) give a compact expression for the frequency of small-amplitude oscillations around P₁, which can be written as (16)

where I_Q1 is the equatorial inclination at P₁ given in Eq. (14). A negative value of means that the equilibrium point is stable. As expected, is negative all over the parameter space. For a zero-obliquity planet, Eq. (16) simplifies to (17)

It shows that the timescale parameter κ defined in Eq. (9) is the oscillation frequency around P₁ for ε = 0 and at a radius a = r_M. See Appendix B for the limit value of ξ₁ in the regions of parameter space where Eq. (16) looks undefined. As a summary, Fig. 7 shows the libration period around P₁ in the whole parameter space. The stability properties of P₂ and P₃ are not crucial for the dynamics of a regular satellite, but they can play a role if the satellite becomes unstable during its orbital migration (see Sect. 4). For this reason, a brief description of P₂ and P₃ is provided in Appendix B.

From P₃ emerges the separatrix that divides the regions of oscillations around P₁ and around P₂. Noting , the extent of the separatrix can be expressed as (18)

in which the two solutions are the minimum and maximum value of I_Q along the separatrix (see Fig. 1). From Eq. (18), we deduce that the width of the island surrounding P₂ is zero in parameter region S₂ (as expected from Fig. 2), and that it increases for growing a and for decreasing cos²ε. This has important consequences for the emergence of chaos discussed in Sect. 4. At the singular point S₁, the island covers the whole sphere.

Apart from the singular regions S₁ and S₂, the continuous behaviour of the Laplace states all over the parameter space is crucial for the long-term satellite dynamics, because if some physical mechanism induces a slow change of parameters (e.g. if the satellite migrates, or if the planet’s spin axis is gradually tilted), then the satellite would adiabatically follow the equilibrium point around which it oscillates, while conserving the phase space area J spanned by its trajectory. If the system never transits through point S₁ (for oscillations around P₁) or S₂ (for oscillations around P₂), then this adiabatic drift can go on as long as the phase space area delimited by the separatrix is wide enough to contain J. This last condition is always verified if J = 0, that is, if the satellite lies exactly on a stable Laplace state.

Fig. 1 Level curves of the Hamiltonian function for a circularorbit. The parameters are a∕rM = 1.1 and ε = 40°. The separatrix is shown by a thicker black curve. The coloured dots represent the three kinds of equilibrium points (‘Laplace states’), labelled as in the text. A dark colour is used for points P1 and P3 lying at δQ = 0 and for P2 lying at δQ = π∕2. A light colour is used for the symmetric equilibrium point that corresponds to the same Laplace state with reversed orbital motion. Figure 2a shows the same phase portrait plotted on the sphere.

Fig. 2 Level curves of the Hamiltonian function with e = 0 plotted on the sphere. The z-axis is along the spin axis of the planet. The x-axis is along the intersection of the equatorial and ecliptic planes (i.e. the line joining both equinoxes of the planet) and directed towards the ascending node of the star. The equator xy plane and the ecliptic plane are highlighted by the two outer grey circles. A point on the sphere represents the tip of the orbital angular momentum of the satellite, which has coordinates (x, y, z) = (sinIQ sinδQ, −sinIQ cosδQ, cos IQ). The colour code is the same as in Fig. 1. Panel a: same parameters as Fig. 1. Panel b: parameter region S1, defined by a∕rM = 1 and ε = 90°. The magenta curve is made of an infinity of stable equilibria resulting from the merging of P1 with the separatrix emerging from P3. Panel c: parameter region S2, defined by a∕rM > 0 and ε = 0° (or 180°). The brown curve is made of an infinity of stable equilibria resulting from the merging of P2 with the separatrix emerging from P3.

Fig. 3 Location of the Laplace states as a function of a∕rM. Each panel features a fixed value of ε (see labels). The level IQ = ε is shown by a horizontal dotted line. The colour code is the same as in Fig. 1. The black vertical line in the central panel shows the degenerate equilibrium circle produced by the merging of P1 and P3.

Fig. 4 Location of the Laplace states as a function of ε. Each panel features a fixed value of a∕rM (from top to bottom: 0.01, 0.95, 1, 1.05, and 20, respectively). Colours have the same meaning as in Fig. 3.

Fig. 5 Summary of the regions of the parameter space governing the dynamics of a massless satellite on a nearly circular orbit. The phase space features three equilibrium points (‘Laplace states’), among which P1 and P2 are stable to inclination variations and P3 is unstable. At point S1 of the parameter space, P1 and P3 degenerate into the equilibrium circle that is stableto inclination variations. In region S2 of the parameter space, P2 and P3 degenerate into the equilibrium circle that is stableto inclination variations. In regions E1, E2, and E3, respectively, P1, P2, and P3 are unstable to eccentricity variations. In region S2 ∪E2, is unstable to eccentricity variations. At point S1, is stable to eccentricity variations provided that the equatorial inclination of the satellite verifies cos2 IQ ≤ 1∕5.

Fig. 6 Inclination of the Laplace state P1 as a function of the parameters (see Eq. (14)). Some level curves are labelled in black. The critical radii and rM and rF are represented in white. At a = rM, P1 lies exactly halfway between the equator and the ecliptic (i.e. IQ = ε∕2 for a prograde spin). The abrupt transition of IQ from 0° to 180° at ε = 90° is a coordinate singularity. The S1 point is a real singularity (see text).

Fig. 7 Period of small inclination oscillations around the Laplace state P1 as a function of the parameters. The period is given by Eq. (16) and represented here in unit of the characteristic period τ = 2π∕κ. Some level curves are highlighted and labelled in black. The singular point S1 is indicated in red.

2.3 Stability to variations in eccentricity

Up to now, we assumed that the eccentricity e of the satellite is zero, which is an equilibrium point. Since the linearised system in the vicinity of e = 0 produces a decoupling between eccentricity and inclination variations, the previous analysis is valid up to order and it neglects . Since the condition e = 0 for a real satellite is never exactly verified, we must consider the stability of the Laplace states to eccentricity variations, that is, we must determine whether a small non-zero offset of eccentricity remains small or grows big over time. As before, we focus on the Laplace state P₁; the description of P₂, P₃, and the degenerate circles and are provided in Appendix B.

The eigenvalues of the eccentricity linear sub-system inform us about their stability against eccentricity growth. Tremaine et al. (2009) give a compact expression for the frequency of small-amplitude eccentricity oscillations around P₁. It can be written (19)

where I_Q1 is the equatorial inclination at P₁ given by Eq. (14). See Appendix B for the limit value of μ₁ in the regions of parameter space where Eq. (19) looks undefined. A negative value of means that the equilibrium point is stable to eccentricity variations. As noted by Tremaine et al. (2009), P₁ is stable in all the parameter space except in a small closed region resembling a cardioid. We call E₁ this region of the parameter space; it can be visualised in Figs. 5 and 8.

The boundary of E₁ is given by the roots of , which have a closed-form analytical expression. We first define two critical radii r₁ and r₂ as (20)

As shown in Fig. 8, the radii r₁ and r₂ mark the leftmost and rightmost limits of E₁, and the boundary of E₁ has a cusp at the singular point S₁. Noting , the boundary of E₁ can be expressed piecewise as (21)

for r₁ ≤ a ≤ r₂, and (22)

for r_M ≤ a ≤ r₂. Equation (22) corresponds to the cusp portion of the curve, and the two portions meet at a = r₂ (see Fig. 8). The obliquity ε ≈ 68.875° quoted by Tremaine et al. (2009) as the minimum value where P₁ can be unstable is reached at . It has actually the following closed-form: (23)

Interestingly, does not go to zero at the singular point S₁, but is discontinuous (see Appendix B). Moreover, the value of at ε = 90° and is the largest (positive) value that can ever reach in the whole parameter space: it is therefore the most unstable location of P₁ to eccentricity variations. This explains the numerical results of Tamayo et al. (2013), who note that for Uranus, whose obliquity is not far from 90°, the radius r_M is the approximatelocation at which the eccentricity grows most rapidly. This also explains why they find that the instability is more violent if the satellite reaches E₁ while migrating inwards rather than outwards (see Fig. 8).

The stability properties of P₂ and P₃ to eccentricity variations are given in Appendix B. We show that they are unstable in the regions E₂ and E₃, respectively, illustrated in Fig. 5. We note that E₂ entirely contains the region E₁; hence, in region E₁ all Laplace states are unstable to at least eccentricity or inclination variations.

Fig. 8 Region E1 of the parameter space where the Laplace state P1 is unstable to eccentricity variations. The background colour shows the normalised value of , with blue for negative values (stable) and red for positive values (unstable). The border of E1 is highlighted with a black contour, obtained using the closed-form expression in Eqs. (21) and (22). Region E1 is entirely contained inside the dashed brown lines, whose left and right limits r1 and r2 are given in Eq. (20), and whose bottom and top limits are given by Eq. (23). The Laplace state P1 is singular atpoint S1, in which is discontinuous.

2.4 Eccentric Laplace states

Along the boundaries of the regions E₁ and E₂, where the Laplace states P₁ and P₂ become unstable to eccentricity variations, Tremaine et al. (2009) show that they both bifurcate into equilibrium configurations with an eccentric orbit. We call these configurations P and P. Likewise, we show in Appendix C that along the two boundaries of the E₃ region (the V-shaped boundary for a < r_L and the drop-like boundary for a > r_L; see Fig. 5 and Appendix B), the Laplace state P₃ bifurcates into two eccentric equilibria that we call P and P, respectively.

In the space of parameters, there exist stable regions for all of these eccentric equilibria. Therefore, satellites reaching the unstable regions E₁, E₂, or E₃ via a smooth parameter change are not bound to destabilise; they can instead bifurcate to a stable eccentric configuration. The properties of the eccentric equilibria are recalled in Appendix C; we provide formulas that can be used to easily compute their locations as a function of the parameters. In our case, we are mostly interested in the equilibrium P, because it bifurcates from the classic Laplace state P₁ in which regular satellites are formed.

At equilibrium P, the orbital angles of the satellite are ω_Q = π∕2 mod π and δ_Q = 0 (or δ_Q = π for the twinequilibrium with reversed orbital motion). The equatorial inclination of the satellite at P can be written as (24)

where we define v as (25)

in which e is the satellite’s eccentricity at equilibrium. We note that the inclination in Eq. (24) has the same form as I_Q1 in Eq. (14), but where is replaced by v. The behaviour of as a function of the parameters is therefore very similar to I_Q1 except that the non-zero eccentricity of the satellite acts like a modified orbital distance. For e = 0, the definitions of and I_Q1 coincide. As shown in Appendix C, the eccentricity at equilibrium can be computed in the general case as a three-dimensional surface with an explicit parametric representation.

The eccentricity and inclination of the satellite at equilibrium P are shown in Fig. 9 as a function of the parameters. We recognise the cardioid-like boundary of the E₁ region. Since Fig. 9 is the projection of a complex three-dimensional surface, a portion of this surface has been cut off for the purpose of the figure. The removed portion of the surface connects to the grey line near the centre of the figure (see the colour discontinuity), and it can be visualised in Fig. 10. Along the cutting line, the right portion of the three-dimensional surface turns round to higher semi-major axes again.

Along the three-dimensional curve defined by (26)

the inclination of the satellite given in Eq. (24) is undefined. This is a real singularity, where P does not exist. Indeed, the curve is the eccentric continuation of the singular point S₁, at which P₁ and P₃ are degenerate. In Fig. 9, the curve is visible between the points labelled S₁ and S. Along this line, the orbital inclination has two different limits (different from 0° and 180°) according to whether the system tends to ε = 90° from below or from above. As shown in Appendix C, the point S is the location where pierces the three-dimensional surface of equilibrium. Noting , the point S has coordinates and ε = 90°. By comparing Figs. 9 and 6, we see that S can be seen as the eccentric counterpart of S₁, where inclination level curves converge. This property will be important in Sect. 3.

Contrary to the circular case, the eccentricity and inclination degrees of freedom are fully coupled at P. Therefore, in the vicinity of P, the eccentricity and inclination of the satellite both vary according to two distinct eigenfrequencies (plus their opposite). The periods of these two oscillation modes in the stable regions are shown in Fig. 11. By comparing with Fig. 7, we see that the oscillation timescale near P has the same order of magnitude as in the circular case. Moreover, we note that one frequency tends to zero at point S, in the same way as tends to zero at S₁. Along the boundary of the E₁ region, the two eigenfrequencies tend to the oscillation frequencies ξ₁ and μ₁ around P₁ given at Eqs. (16) and (19), confirming that the eccentric equilibrium P bifurcates from the circular equilibrium P₁.

As stressed by Tremaine et al. (2009), the eccentric equilibrium P is stable near its bifurcation from P₁ and in the central region of Fig. 11. These properties will be important for the future evolution of Titan described in Sect. 4. On the top tube-like portion of the equilibrium surface (not shown in Fig. 11), we show in Appendix C that P is mostly unstable, even though a small stable region exists at very high eccentricities.

Before concluding this section, we stress that the linear instability of a Laplace state does not necessarily mean that the satellite’s trajectory is chaotic, and it gives no information about the amount of eccentricity and inclination increase suffered by the satellite. Interestingly, the simulations of Tremaine et al. (2009), Tamayo et al. (2013), and Speedie & Zanazzi (2020) reveal more chaos than expected in the E₁ region, even where the eccentric equilibrium P should theoretically be stable. In the case of trans-Neptunian objects perturbed by the galactic tides (see Table 1), Saillenfest et al. (2019a) find that at a ≈ r_M the phase space is covered by chaos, allowing for transitions between circular and quasi-parabolic orbits. The emergence of violent chaos in the orbit of Titan is confirmed numerically in Sect. 4. But before speaking of chaos, we must first understand the mechanism through which Titan is brought into the unstable region. In the next section, we see that it results from an interplay between the dynamics of Titan’s orbit and Saturn’s spin axis.

Fig. 9 Eccentricity and inclination of the satellite at the eccentric equilibrium P. The three-dimensional surface of equilibrium has been cut along the grey line, as shown in Fig. 10. In the bottom panel, some level curves are plotted in white (see labels). Noting , the extremity of the grey cutting line have coordinates and at the top and bottom points, and and ε = 90° at the middle.

Fig. 10 Eccentricity at the eccentric equilibrium P seen as a three-dimensional surface. The colour is the same as in Fig. 9, top panel. In Fig. 9, the top tube-like portion of the surface has been cut off along the grey line. As detailed in Appendix C, the top portion extends up to a → ∞, where it tends to e = 1. Sections of this surface can be seen in Fig. 6 of Tremaine et al. (2009).

Fig. 11 Period of small oscillations around P as a function of the parameters. There are two oscillation modes, shown in the two panels. In the hatched area, the equilibrium is unstable. The singular points S1 and S are labelled, and the cutting line of the three-dimensional surface is shown by a grey line (same as Fig. 9). Three-dimensional figures are provided in Appendix C, where the geometry of the stable regions is clearer.

3 Spin-axis dynamics of the planet

In the previous section, the spin axis of the host planet was assumed to be fixed in an inertial frame. Actually, because of the torque applied by the star and the satellite on its equatorial bulge, the spin axis of the planet is made to slowly precess over time. In this section, we aim to get a qualitative understanding of the effect of the satellite on the spin-axis motion of its host planet, with an eye on the case where the planet is locked in a secular spin–orbit resonance, that is, where additional perturbations maintain the planet’s spin-axis precession frequency to a fixed value.

A self-consistent model for the dynamics of a satellite and the spin axis of its host planet has been derived by Boué & Laskar (2006): under the assumption that the satellite’s argument of pericentre stably circulates, they obtained a full analytical characterisation of the averaged dynamics, which was proven to be integrable. However, this model does not hold if the system is affected by additional perturbations. In particular, mutual interactions between planets result in their nodal and apsidal precession motions (see e.g. Murray & Dermott 1999), whose multiple modes and harmonics are responsible for the secular spin–orbit resonances. Besides, the assumptions of Boué & Laskar (2006) cannot apply if the system reaches the region E₁, as the satellite’s pericentre can become stationary near the equilibrium (see Sect. 2). Consequently, the model of Boué & Laskar (2006) will serve us as a reference for the ‘instantaneous’ value of the secular spin-axis precession rate of the planet, but it cannot be used (as such) to describe the dynamics inside a secular spin–orbit resonance.

In this section, we first recall the properties of secular spin–orbit resonances (Sect. 3.1), and then we study the effect of a satellite on a resonantly locked planet (Sect. 3.2).

3.1 Secular spin–orbit resonance

In the approximation of rigid rotation, the secular spin-axis dynamics of an oblate planet is ruled by the Hamiltonian function (27)

where the conjugate canonical coordinates used here are X = cosε (cosine of obliquity) and − ψ (minus the precession angle). The first part comes from the torque exerted by the star on the equatorial bulge of the planet at quadrupolarorder. It can be written (28)

where the parameter α is called the ‘precession constant’. In the absence of satellite, the expression of the precession constant is given for instance by Néron de Surgy & Laskar (1997) as (29)

where ω is the spin rate of the planet and λ is its normalised polar moment of inertia. The parameters J₂ and λ are related to the moments of inertia A, B, and C of the planet through (30)

where M is the mass of the planet. The second part of the Hamiltonian function stems from the motion of the planet’s orbital plane, produced for instance through mutual perturbations with other planets. It can be written (31)

where , , and are explicit functions of time whose expressions in terms of the planet’s classical orbital elements are given for instance by Laskar & Robutel (1993) and Néron de Surgy & Laskar (1997). If the planet’s orbit were fixed, then would be identically zero.

We assume that the orbit of the planet is long-term stable, so that its secular orbital motion can be expressed (at least locally) in convergent quasi-periodic series. Truncating the series describing its orbital inclination motion to N terms, it can be expressed as (32)

where S_k is a positive real constant and ϕ_k evolves linearly over time with frequency ν_k, that is, (33)

for any k = 1, 2...N. In Eq. (32), I and Ω are the orbital inclination and the longitude of ascending node of the planet measured in an inertial reference frame, not to be confused with the satellite’s orbital elements used in Sect. 2. As shown by Saillenfest et al. (2019b), the Hamiltonian function is proportional to the amplitudes S_k of the quasi-periodic decomposition. Therefore, if the planet is not much inclined with respect to the invariable plane of the system (i.e. S_k ≪ 1 for ν_k ≠ 0), which is what we expect in a long-term stable planetary system, then the Hamiltonian in Eq. (27) can be considered as a perturbation to the unperturbed Hamiltonian . In this setting, the long-term spin-axis dynamics of the planet is shaped by resonances between the unperturbed spin-axis precession frequency and the forcing frequencies appearing in Eq. (33). In the Solar System, the orbital precession motions of the terrestrial planets contain numerous large-amplitude harmonics, creating a collection of wide secular spin–orbit resonances which overlap with each other and create wide chaotic zones (Laskar & Robutel 1993). The orbital precession motions of the giant planets, on the contrary, are composed of many fewer strong harmonics, so that large secular spin–orbit resonances are rare and isolated from each other. Depending on their spin-axis precession frequency, the giant planets of the Solar System can therefore be captured into an isolated resonance and oscillate stably within its separatrix (see e.g. Ward & Hamilton 2004; Ward & Canup 2006; Saillenfest et al. 2020, 2021b).

Using a perturbative approach, Saillenfest et al. (2019b) have described the properties of all resonances up to order three in the amplitudes {S_k}. The largest resonances are those of order 1, for which the resonance angle is σ = ψ + ϕ_j, where j is a given index in the orbital series in Eq. (32). Second-order resonances involve two terms in the series, and third-order resonances involve three. Eccentricity-driven resonances only appear at order three and beyond. In any case, the resonance angle is a linear combination involving ψ and one or several ϕ_k. If the planet is trapped inside one of those resonances, then the resonance angle oscillates around a fixed value, which means that the spin-axis precession frequency is forced to remain approximatively constant, equal to a combination of frequencies ν_k.

In the vicinity of a first-order secular spin–orbit resonance, the Hamiltonian function reduces to the well-known ‘Colombos’s top Hamiltonian’ (Colombo 1966; Henrard & Murigande 1987). This Hamiltonian can be written (34)

where the conjugate coordinates are X = cosε and − σ = − ψ − ϕ_j (i.e. minus the resonance angle). Neglecting terms of order four and higher in the amplitudes {S_k}, the parameters γ and β in Eq. (34) are (35)

where is the characteristic spin-axis precession frequency of the planet. Contrary to Saillenfest et al. (2019b), we do not expand the eccentricity variations of the planet in quasi-periodic series: since eccentricity variations only appear at third order in the amplitudes, they are not important for our present qualitative description of the dynamics. Written as in Eq. (35), eccentricity variations simply produce slight fluctuations in the resonance parameters γ and β.

As defined in Eq. (35), for small amplitudes S_k, the parameters γ and β are both positive if ν_j < 0. This corresponds to a prograde resonance, for which the resonance centre is located at an obliquity ε ≤ 90°. An example is presented in Fig. 12. On the contrary, retrograde resonances are obtained for ν_j > 0, for which γ and β are negative. Due to symmetries, changing the sign of γ is equivalentto replacing X by − X, and changing the sign of β is equivalent to replacing σ by σ + π. Following Peale (1969), the equilibrium points are usually called ‘Cassini states’, numbered from 1 to 4, as labelled in Fig. 12.

Henrard & Murigande (1987) showed that the phase space has two different topologies according to whether γ^2∕3 + β^2∕3 is smaller or larger than 1. If γ^2∕3 + β^2∕3 > 1, then there is no resonance (i.e. no separatrix) and only the Cassini states C₂ and C₃ are present (see Fig. 13a). If γ^2∕3 + β^2∕3 < 1, then all four Cassini states are present, and C₂ becomes the resonance centre (see Fig. 13b). As noted by Saillenfest et al. (2019b), in some parameter region, the resonance contains the north pole and/or the south pole of the sphere. We stress that the resonance can be quite large even for a moderate amplitude S_j.

As detailed below, the effect of satellites on the spin-axis motion can be modelled by replacing α in Eq. (28) by an effective precession constant α′, whose value depends on the distance of the satellites. Therefore, we are interested in how the geometry of the phase space is changed when modifying the parameter α appearing in Eq. (35) via p. Analytical formulas for the locations of the equilibrium points and the separatrix are provided by Saillenfest et al. (2019b) or Haponiak et al. (2020); we recall here their behaviour when varying the parameter α. Due to symmetries, we only describe the case of a prograde resonance, for which γ and β are positive.

Figure 14 shows the generic behaviour of the system with respect to the parameter 1∕γ, which is proportionalto α. For α → 0, both γ and β tend to infinity, so only Cassini states 2 and 3 are present (see Fig. 13a). Their asymptotic locations are (36)

which, for a small amplitude S_j, represents just a small offset from ε₂ = 0° and ε₃ = 180° (see Fig. 14). In other words, for α → 0 the resonance is infinitely far away and it only contributes through a residual shift of the equilibrium points at the north and south poles of the sphere. Therefore, the obliquity is almost constant and σ circulates.

If we increase α above zero, the Cassini state C₂ is tilted away from the north pole of the sphere (see Figs. 13a and 14); then, when γ^2∕3+ β^2∕3 becomes smaller than1, the Cassini states C₁ and C₄ appear together with the separatrix. As shown in previous articles (see in particular Figs. B.1 and B.2 of Saillenfest et al. 2020), for increasing α the resonance first contains the north pole of the sphere, and then the separatrix crosses the north pole and moves down to larger obliquities. This transition is visible in Fig. 14 as the very narrow interval of 1∕γ in which the pink region touches ε = 0°.

For α →∞, the parameters γ and β both tend to zero and the location of the Cassini states tend to (37)

Moreover, for α →∞ the separatrix enclosing C₂ becomes vanishingly narrow and it merges with the Cassini states C₂ and C₄, producing the degenerate equilibrium circle shown in Fig. 13c. For large but finite values of α, we note that the resonance width goes beyond ε = 90° (there is no topological boundary at ε = 90° for non-zero libration amplitudes).

Hence, as a summary, if we increase the parameter α from 0 to ∞, the Cassini state C₂ gradually passes from ε₂ ≈ 0° (without separatrix) to ε₂ = 90° (inside the resonance separatrix). These properties will be important below.

Fig. 12 Spin-axis dynamics in the vicinity of a first-order secular spin–orbit resonance. Some level curves of the Hamiltonian function in Eq. (34) are represented in black or in red, according to whether they are outside or inside the resonance, respectively. The separatrix is shown by a thick black curve. The constant parameters are γ = 0.75 and β = 0.03. The equilibrium points (‘Cassini states’) are represented by coloured dots. Figure 13b shows the same phase portrait plotted on the sphere.

3.2 Effect of a satellite

The orbital angular momentum of the planet usually greatly exceeds its rotational angular momentum. Therefore, the planet’s orbit remains almost unaffected by the spin-axis precession motion. In Sect. 3.1, this property allowed us to treat the orbital variations as a forcing term in the spin-axis dynamics. Now, if the planet has a satellite that lies in its local Laplace plane (i.e. if it is in the Laplace state P₁ described in Sect. 2.2), then the planet’s spin axis and the satellite’s orbit rigidly precess as a whole about the planet’s orbital angular momentum (Boué & Laskar 2006). In Eq. (4), this precession would be traduced by a drift of Ω_P over time.

We define the characteristic spin-axis precession timescale of the planet as T = 2π∕p, where has been introduced in Sect. 3.1. Without satellite, the free spin-axis precession frequency of the planet is simply (38)

as obtained from Eq. (28). The characteristic timescale T is given in Table 1 for some planets of the Solar System. The value of T can be compared to the characteristic timescale τ of the satellite’s orbital dynamics. We see that τ ≪ T for all satellites listed in Table 1. This large separation of timescales justifies the approximation made by Tremaine et al. (2009) and used in Sect. 2 to consider that the planet’s equatorial reference frame is inertial despite its precession motion. For a satellite oscillating about the Laplace state P₁, this timescale condition may be violated at the border of region E₁ or in the extreme vicinity of the singular point S₁, that is, where one of its orbital oscillation frequencies tends to zero (see Fig. 7). The behaviour of satellites in this critical regime will be investigated numerically in Sect. 4.

Substantially massive satellites contribute to the spin-axis precession frequency of their host planet. If the satellites oscillate about the Laplace state P₁ with τ ≪ T, French et al. (1993) give an elegant expression for their contribution: one must simply replace J₂ and λ in Eq. (29) by the effective values (39)

where m_k, a_k and n_k are the mass, the semi-major axis, and the mean motion of the kth satellite, and L_k is the inclination of the Laplace plane of the kth satellite with respect to the planet’s equator⁵. Using these expressions, the free spin-axis precession frequency Ω₀ of the planet is still obtained from Eq. (28), but where α is replaced by an effective precession constant α′. We stress that Eq. (39) is valid whatever the distance of the satellite, and not only in the close-in regime considered previously by Saillenfest et al. (2020, 2021b). It only requires that the satellites oscillate around the Laplace state P₁, which is what we expect for any regular satellite (see e.g. Tremaine et al. 2009), unless it reaches the high-obliquity unstable region E₁ (see Sect. 2.3).

For small satellites, the value of L_k can be directly taken from Eq. (14). In this case, the model is not self-consistent, because the satellite is considered to be massless when dealing with the orbital dynamics (Sect. 2) but massive when computing its long-term influence on the planet’s spin axis. Yet, because τ ≪T (see Table 1), this approximation results to be very accurate for satellites having a small mass ratio m_k ∕M. For Titan, whose mass is about 10⁻⁴ of Saturn’s, the inclination L_k and the precession rate Ω₀ obtained through Eqs. (14) and (39) are very close to those obtained using the self-consistent (but quite complicated) model of Boué & Laskar (2006), as detailed in Appendix D. The approximation is less good for very massive satellites like the Moon (m_k∕M ≈ 0.01), but we can check that the qualitative picture described below remains valid, which means that our analysis captures the essence of the dynamics even for large satellite-to-planet mass ratios.

We note that Eq. (39) looks undefined for ε = 0 or 90°. However, computing L_k using Eq. (14), the contribution of satellites around a zero-obliquity planet simplifies to (40)

Likewise, for ε = 90° the parameter becomes (41)

However, since Ω₀ is factored by cos ε, the spin-axis precession frequency of the planet is in any case zero for ε = 90°, except at a = r_M, where the Laplace state P₁ of the satellites is undefined (see Sect. 2).

In order to keep the discussion as general as possible before focussing on particular bodies, one more approximation can be made. Indeed, even for a fast-spinning planet like Saturn, the oblateness coefficient J₂ appearing in Eq. (29) is a small parameter compared to λ (whose order of magnitude is slightly less than unity). As a result, the relative increase in J₂ produced by satellites in Eq. (39) is much larger than the relative increase in λ. This discrepancy is amplified bythe factor n_k∕ω appearing in Eq. (39), which further reduces the satellites’ contribution to λ′. Therefore, as a first approximation, the contribution of satellites to λ′ is negligible compared to their contribution to . Assuming that λ′≈ λ, and considering that the planet has a single main satellite, the free precession frequency of the planet’s spin axis simplifies to (42)

where we have introduced the ‘mass parameter’ η of the satellite, defined as (43)

The mass parameters of some satellites in the Solar System are given in Table 1. We stress that even a low-mass satellite can have a large mass parameter η. For instance, Titan has a mass of m∕M ≈ 10⁻⁴ but a mass parameter η ≈ 10, and Titan greatly affects Saturn’s spin-axis dynamics (see below). Therefore, contrary to what one could think a priori (see e.g. Li & Batygin 2014), a large mass ratio m∕M is not necessarily required to substantially alter a planet’s obliquity. Using Eq. (42), the spin-axis precession rate of the planet normalised by p is only a function of ε, η, and a∕r_M. We can therefore study its behaviour in a very generic way, as we did for the satellite’s dynamics in Sect. 2.

Figure 15 shows the spin-axis precession frequency of the planet as a function of the distance of its satellite. We retrieve the classical curve shown in Fig. 5 of Boué & Laskar (2006), with the close-in and far-away satellite regimes. For a →0 and a →∞, the precession frequency tends to the value pcosε obtained without satellite. In between, the satellite enhances the precession frequency of the planet by an increment that is proportional to its mass parameter η (see Eq. (42)). For a given obliquity, the maximum value of Ω₀ divides the close-in and far-away satellite regimes, characterised by the well-known power laws in a² and a⁻³, respectively. The magnitude of Ω₀ differs when varying the satellite’s mass parameter η, but the location of its maximum asa function of a is independent of η. As shown in Fig. 15, the maximum of Ω₀ is located somewhat below the midpoint radius a = r_M. By analysing Eq. (42), we find that the maximum of the curve is located at a = r_F, defined in Eq. (15). We recall that r_F is the inflexion point of the satellite’s inclination, illustrated in Fig. 6. We see here that the spin-axis precession frequency of the planet is intimately linked to the properties of the Laplace state P₁ of its satellite. We further analyse this relation below.

Figure 16 shows the value of Ω₀ at its maximum (i.e. at a = r_F) as a function of the planet’s obliquity. The global maximum of Ω₀ is reached at ε = 0°, for which . It has value (44)

For an obliquity ε → 90°, the maximum of Ω₀ is reached at a = r_M, that is, at the singular point S₁ of the satellite’s dynamics (see Sect. 2). Figure 16 shows that the maximum value of Ω₀ has two different limits at point S₁ according to whether the obliquity tends to 90° from above or from below. These limits are Ω₀ = ± Ω_S, where (45)

We note that these limits are non-zero and well defined. This is far from obvious when modelling the satellites by an effective precession constant α′, as one would never guess that α′cosε tends to a non-zero finite value when cosε tends to zero (α′ actually goes to infinity). We deduce that when the system approaches the singular point S₁, the notion of ‘precession constant’ loses its meaning. In this regime, classical figures drawn in the plane (ε, α), like those used by Saillenfest et al. (2021b) are inappropriate.

We are interested in the level curves of Ω₀ as a function of the parameters. Indeed, if the planet is trapped in a secular spin–orbit resonance during the migration of its satellite, its precession frequency is forced to oscillate around a constant value while the parameters a and ε vary (see Sect. 3.1). This is likely the case for Saturn and Titan (see Saillenfest et al. 2021b), and it will probably be the case for Jupiter and its satellites in the future (see Saillenfest et al. 2020). Figure 17 shows the spin-axis precession period of the planet in the parameter space for different values of the satellite’s mass parameter η. We choose here to plot the precession period |2π∕Ω₀|, instead of the frequency Ω₀ for easier comparison with the satellite’s oscillation period shown in Fig. 7. Values of τ, η and T for real bodies can be found in Table 1. As shown in Appendix D, Fig. 17 presents a very good agreement with the precession period obtained in self-consistent models.

For a very small mass parameter η ≪ 1 (e.g. for Deimos), the effect of the satellite is negligible and the level curves of Ω₀ would appear perfectly horizontal in Fig. 17. Therefore, even if the planet is trapped in a secular spin–orbit resonance, its mean obliquity would remain unaltered over the migration of its satellite. For a substantial value of η, on the contrary, Fig. 17 shows that the level curves of Ω₀ lose their horizontal shape and rearrange around the ridge line a = r_F, where the satellite’s contribution is maximum. Indeed, since the factor in front of η in Eq. (44) is close to 1∕2, the relative difference between Ω_max and Ω_S is very small for a large value of η, namely (46)

Therefore, a large value of η implies that Ω₀ is approximatively constant along the ridge line a = r_F (i.e. the two curves in Fig. 16 are nearly horizontal). Consequently, the ridge line a = r_F creates a barrier that cannot be crossed by the level curves of Ω₀: instead, the level curves must go around the ridge line, and they converge at the singular point S₁ (see Fig. 17). More precisely, all level curves with frequency values |Ω₀ |≤ Ω_S converge to S₁. For a large mass parameter η, this condition is verified for almost all level curves of Ω₀. This property is related to the orbital plane of the satellite, which can reach S₁ from any inclination (see Fig. 6).

Figure 17 shows that for a large enough value of η, numerous level curves of Ω₀ connect the singular region S₂ (i.e. ε = 0) to the singular point S₁. Therefore, if the planet is trapped in a secular spin–orbit resonance, the migration of its satellite through a = r_M forces its obliquity to increase all the way from 0° to 90°. We see that such an extreme obliquity increase can take place over a very short migration range for the satellite (e.g. in panel d if its distance decreases from a ≈ 1.05r_M to a = r_M). The theoretical limit for the obliquity reached through the mechanism described by Saillenfest et al. (2021b) is therefore 90°, or even more if the resonance is large and its width extends beyond ε = 90° (see Sect. 3.1). If ever the system manages to go through the singular point S₁ in some way, one could even imagine a scenario where the planet then picks a retrograde resonance (see e.g. Kreyche et al. 2020) and goes on tilting up to 180°.

The level curves of Ω₀ going from ε = 0° to ε = 90° verify p ≤ Ω₀ ≤ Ω_S. Therefore, the minimum mass parameter allowing the planet to tilt all the way from 0° to 90° through the migration of its satellite is obtained by putting Ω_S = p in Eq. (45), which gives η = 2. In other words, η = 2 is the minimum mass parameter for which the level curve Ω₀ = p, which starts at (a, ε) = (0, 0°), connects to (a, ε) = (r_M, 90°): see the level curve labelled ‘1’ in Fig. 17. In Table 1, the condition η ≥ 2 is verified for the Moon, Titan, and Oberon, but not for Deimos, Callisto, and Iapetus.

In practice, since the singular point S₁ is surrounded by the unstable region E₁ (see Sect. 2.3), we expect the satellite’s orbit to become unstable before actually reaching S₁. Moreover, the width of the secular spin–orbit resonance of the planet decreases as the obliquity increases (see Fig. 14), and since many level curves converge to S₁, other secular spin–orbit resonances (if any) would necessarily overlap at some point and create a chaotic region. For these two reasons, we expect the planet to be released out of resonance before actually reaching ε = 90°. This double destabilisation of the planet and of its satellite is investigated in the next section.

Fig. 13 Level curves of the Hamiltonian plotted on the sphere. The planet’s orbit lies in the xy-plane. The obliquity ε is the tilt from the z-axis, and the resonance angle σ is the polar angle measured in the xy-plane. The colour code is the same as in Fig. 12. Panel a: geometry for γ2∕3 + β2∕3 > 1. Panel b: geometry for γ2∕3 + β2∕3 < 1 (same parameters as in Fig. 12). Panel c: geometry for γ and β tending to zero.

Fig. 14 Bifurcation diagram of the system as a function of 1∕γ, which is proportional to α. The second parameter used in this example is β = 0.03γ. The Cassini states are labelled with the same colour code as previous figures.

Fig. 15 Spin-axis precession frequency of a planet orbited by a single satellite. The frequency is given by the simplified expression in Eq. (42) and represented here in unit of the characteristic frequency p, for different obliquity values (see labels). In this example, the mass parameter of the satellite is set to η = 10, which is close to Titan’s value (see Table 1).

Fig. 16 Maximum spin-axis precession frequency of the planet as a function of its obliquity. The curve is obtained by injecting a = rF from Eq. (15) in Eq. (42). The extreme values Ωmax and ΩS are given in Eqs. (44) and (45). For ε = 90°, the value at the maximum is undefined.

Fig. 17 Spin-axis precession period of a planet orbited by a single satellite. The period is given by Eq. (42) and represented here in unit of the characteristic period T = 2π∕p. For ε >90°, the spin-axis precession frequency is negative. Each panel corresponds to a given value of the satellite’s mass parameter η (see Eq. (43)), as labelled in the top right corner. Some level curves are highlighted in black. The singular point S1 and the ridge line a = rF are indicated in red.

4 Saturn and Titan

4.1 Overview of the dynamics

The Hamiltonian function in Eq. (27) explicitly 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. As in previous studies, we use the secular solution of Laskar (1990) expanded in quasi-periodic series, that is, under the form given in Eq. (32). The ten largest terms of the ζ series of Saturn are shown in Table 2, ordered by decreasing amplitude.

As explained in Sect. 3, a satellite is able to tilt a planet from ε = 0° to 90° if its mass parameter η is larger than 2. This condition is met for Titan, which has η ≈ 12.4. Furthermore, the maximum spin-axis precession frequency of the planet is reached along a = r_F, and the global maximum Ω_max is given by Eq. (44). For Saturn and Titan (see Table 1), we obtain Ω_max ≈ 1.41″ yr⁻¹. Therefore, all frequencies in Saturn’s orbital decomposition whose magnitude exceeds this value are unreachable by Saturn, whatever the distance of Titan. This only leaves a handful of possible first-order secular spin–orbit resonances, even when considering the full quasi-periodic solution of Laskar (1990). The largest resonance is with ν₃ = s₈ (see Table 2). Other resonances can be identified in Table A.2 of Saillenfest et al. (2021b): there are two prograde resonances (ν₁₉ and ν₅₁) and two retrograde resonances (ν₂₈ and ν₄₄). As revealed by their high index in the orbital series, these four resonances are very small. This explains why no relevant second- or higher-order resonance can possibly affect Saturn. Then, we know that all precession frequencies Ω₀ verifying p ≤ Ω₀ ≤ Ω_S have a level curve that connects ε = 0° to 90° (or ε = 180° to 90° for a retrograde spin). For Saturn and Titan, we have p ≈ 0.19″ yr⁻¹ and Ω_S ≈ 1.20″ yr⁻¹, so this condition is met by all five resonances mentioned above (even though ν₅₁ is right at the limit, since |ν₅₁|≈ Ω_S).

Figure 18 shows the location and width of the first-order secular spin–orbit resonances reachable by Saturn as a function of the distance of Titan. The full effect of Titan on Saturn’s spin-axis is taken into account, including its contribution in λ′ (see Eq. (39)). For each resonance, the location of the Cassini states and the separatrix are obtained using the exact analytical formulas of Saillenfest et al. (2019b); however, since Titan’s contribution to the precession constant α itself depends on the obliquity ε (see Eq. (39)), the equations become implicit and must be solved numerically (e.g. with the bisection method). As expected, all five resonances converge at the singular point S₁. If the planet is trapped in a secular spin–orbit resonance during the migration of its satellite, we see that it can behave very differently according to whether the satellite migrates outwards or inwards. For an outward migration, the system goes straight across the unstable region E₁ before reaching the singularity S₁; therefore, the satellite is expected to become gradually eccentric and eventually destabilise (see Sect. 2.4 and Tremaine et al. 2009). On the contrary, for an inward migration, the system can go very close to S₁ before being brutally destabilised at the singularity.

Figure 18 can be compared to Figs. 1 and 17 of Saillenfest et al. (2021b), drawn in terms of Saturn’s effective precession constant (the vertical and horizontal axes are inverted). Figure 17 of Saillenfest et al. (2021b) is a good example of how using a formalism with the precession constant α can be misleading when the satellite gets close to its Laplace radius. Indeed, because of the frequency cut at Ω_max ≈ 1.41″ yr⁻¹, Saturn would be unable to reach ν₁₄ and ν₁₅ for any distance of Titan, even though the trajectory in Fig. 17 of Saillenfest et al. (2021b) appears at the same height as ν₁₄ and ν₁₅ on the graph. Moreover, as shown in Sect. 3.2, α can tend toinfinity while the spin-axis precession frequency remains finite, which is quite counter-intuitive. In order to prevent misinterpretations, we advocate avoiding using α as parameter when the satellite approaches a = r_M, and using the general formula in Eq. (39) or the model of Boué & Laskar (2006) for the spin-axis precession.

Assuming that Saturn follows the centre of the resonance with s₈ (Cassini state 2) and Titan remains in its Laplace plane as it migrates (Laplace state 1), all the properties of the system can be monitored through the analytical formulas given in Sects. 2 and 3. The general behaviour of the system is presented in Fig. 19. On the top left panel, we see that the close-satellite approximation used in previous articles underestimates the obliquity increase of Saturn, even though it remains valid up to an obliquity of about 60°. As detailed in Sect. 2, we note that Titan’s inclination is very different according to whether it reaches the singular point at a = r_M from above or from below. For an outward migration, Titan reaches the singularity with a quite moderate inclination of about 15°. Yet, since the width of the separatrix enclosing Laplace state 2 tends to 180° at a = r_M, this gives an idea of the large inclination variations expected if Titan deviates from the exact equilibrium point, for instance because of a coupling with the eccentricity becoming unstable (see Sect. 2.3). The right column of Fig. 19 gives an idea of the large separation of timescales between the dynamics of Titan and of Saturn’s spin axis. The two timescales become commensurate only in the vicinity of S₁, where the libration period of σ₃ tends to zero while Titan’s libration periods tend to infinity because of the nearby instability. On the top right panel, we can appreciate how Saturn’s spin-axis precession period is maintained to a constant value as it enters the resonance.

Titan is located today at a mean semi-major axis of a₀ ≈ 0.49 r_M. If Titan goes on migrating outwards, it will eventually reach a = r_M at some time in the future. The simplified migration law for Titan provided by Lainey et al. (2020) is (47)

where t₀ is Saturn’s current age and b is a real parameter. According to this formula, Titan will reach a = r_M after an interval of time from today equal to (48)

The astrometric measurements of Lainey et al. (2020) yield values of b ranging in [0.18, 1.71], and their radio-science experiments yield values ranging in [0.34, 0.49]. We deduce that if Titan goes on migrating as expected, it will reach a = r_M between about 2.4 and 230 Gyr from now according to astrometric measurements, and in about 15 to 33 Gyr according to radio-science experiments. Figures 18 and 19 show that the system will first reach the unstable region E₁ when a ≈ 0.83 r_M, that is, between about 1.6 and 80 Gyr from now according to astrometric measurements and in about 8.7 to 17 Gyr according to radio-science experiments. These large values show that the system is unlikely to destabilise before the Sun leaves the main sequence. Yet, Saturn and Titan are not located exactly at their respective equilibrium points, but they oscillate around them with substantial amplitudes. As shown by Tamayo et al. (2013), this can speed up the destabilisation process.

Because of the instabilities described above, Saturn and Titan are not expected to exactly follow Cassini state 2 and Laplace state 1, especially when they approach the singularity point S₁. For this reason, Fig. 19 can only provide a qualitative picture of the evolution of the system. Because of the intricate and multi-timescale nature of the dynamics, building a self-consistent numerical model for the evolution of Saturn and Titan is challenging: the system involves the orbital dynamics of Titan and Saturn, the spin-axis dynamics of Saturn torqued by the Sun and by Titan, and other planets of the Solar System should be included as well in some way to produce the multi-harmonic orbital precession of Saturn. Such a numerical model should also be fast enough to be usable for gigayear propagations (while Titan’s orbital period today is only a few weeks). The design of such a model and the statistical exploration of the chaotic behaviour of Titan and Saturn are left for future works. Yet, a precise idea of the outcomes of the instability can already be obtained by mixing the two simplified models presented in Sects. 2 and 3: on the one hand we explore numerically the dynamics of Titan when Saturn’s spin-axis drifts inside the resonance, and on the other hand we explore the behaviour of Saturn while Titan migrates in the Laplace surface.

Fig. 18 Location and width of every first-order secular spin–orbit resonance reachable by Saturn over the migration of Titan, with amplitudes down to 10−8. Each resonant angle is of the form σj = ψ + ϕj where ϕj has frequency νj labelled on the graph according to its index in the orbital series (see Table 2). For a given value of Titan’s semi-major axis, the interval of obliquity enclosed by the separatrix is shown in pink. The centre of the resonance (i.e. Cassini state 2) with ϕ3 is shown by a red line and the upper and lower separatrices are highlighted in blue. The ridge line a = rF dividing the close-in and far-away satellite regimes is represented by a dashed blue line. The region E1 where the satellite’s Laplace state P1 is unstable isshown in black. The singular point S1 is shown by a red dot.

Fig. 19 Evolution of Saturn’s spin axis and Titan’s orbit while following the centre of the secular spin–orbit resonance with s8. We use blue for Saturn’s spin axis and red for Titan’s orbit. Each panel is described directly on the graph when the legend is not self-explanatory. In the right column, time is shown both in normalised units (left vertical axis) and in physical units (right vertical axis). The conversion factors are given in Table 1. All curves are obtained through the analytical formulas described in the text, assuming that Saturn closely oscillates near Cassini state 2 and Titan closely oscillates near Laplace state 1. In the top left panel, the green curve shows the location of the resonance centre according to the formula used by Saillenfest et al. (2021b) valid for a close satellite. In the bottom right panel, the black curve in the unstable region corresponds to the period needed for the eccentricity to be multiplied by exp (2π) ≈ 535.

4.2 Titan’s orbit

Using the Hamiltonian function in Eq. (1), our setting is similar to the migration simulations of Tremaine et al. (2009), except that both the semi-major axis of Titan and the obliquity of Saturn evolve over time as slow-varying parameters. As a first approximation, their evolution law is provided by the top-left panel of Fig. 19 (blue curve). Therefore, in our simulations we make the obliquity of Saturn and the semi-major axis of Titan vary simultaneously, the latter evolving according to the migration law given by Eq. (47). We tried various values of b in the full uncertainty range [0.18, 1.71], but the statistics of the simulations result to be absolutely independent of Titan’s migration velocity. Indeed, the gigayear timescale of Titan’s orbital expansion always remains extremely large as compared to the timescale of its secular dynamics (τ ≈ 4501 yr; see Table 1 and Fig. 7). Yet, when Titan reaches the unstable region, we do obtain several possible outcomes due to the intrinsic chaotic divergence of trajectories.

Figure 20 shows two examples of simulations. As expected, Titan closely follows its local Laplace plane (blue curve), until it reaches the neighbourhood of the unstable region E₁. At this point, Titan’s eccentricity increases, as the trajectory wanders in the vicinity of the stable eccentric equilibrium P (red curve). Then, P becomes unstable, as shown by the hatched region in Fig. 11, and Titan’s evolution becomes chaotic. This is where the two solutions in Fig. 20 diverge. In the case labelled ‘Solution 1’, Titan jumps right away to a trajectory reaching very high eccentricity and inclination values. In the case labelled ‘Solution 2’, Titan catches the eccentric equilibrium P when this equilibrium becomes stable again (see Fig. 11), but it eventually goes back to the unstable zone where it reaches high eccentricity and inclination values. In this example, a new transition occurs just before the end of the simulation, and Titan’s equatorial inclination starts oscillating between 0° and 180°. In both cases presented in Fig. 20, the ecliptic inclination (not shown) at the end of the integration also oscillates roughly between 0° and 180°. Such extreme inclination variations are not surprising: in the circular case, Fig. 2 shows how the level curves of the Hamiltonian pass from a horizontal structure for ε = 0° (panel c) to a vertical structure at S₁ (panel b), where the orbit can roll all the way around its nodal line. This is traduced by the extent of the separatrix that reaches 180° at S₁ (see the bottom left panel of Fig. 19).

For comparison, we also performed direct unaveraged numerical integrations of the restricted three-body problem including Saturn’s oblateness. In order to reproduce the drift in Titan’s semi-major axis, we added in its equations of motion a small additional acceleration that depends on its velocity, and Saturn’s obliquity is varied accordingly. No orbital or spin-axis precession motions for Saturn are included. In these simulations, Titan’s migration is sped up by a factor of about 500 as compared to Eq. (47), which yields reasonable computation times (a few days or so). As explained above, this acceleration factor is justified by the extremely large separation between Titan’s orbital timescale (thousands of years) and its migration timescale (billions of years). Figure 21 shows two examples of such simulations, chosen for their similarity with Fig. 20. They show that the secular model truncated at quadrupole order used throughout this article does capture the essence of the dynamics. In particular, the evection and ‘ivection’ resonances reported by Speedie & Zanazzi (2020) to produce additional unstable regions are not found to play any role in Titan’s future evolution.

Of course, when Titan’s eccentricity increases and begins to oscillate widely, our model breaks down because the influence of Titan on Saturn’s spin-axis precession is no longer given by Eq. (39). As a result, Saturn’s obliquity should no longer follow the law given in Fig. 18. In order to explore the dynamics in the unstable region, Saturn’s spin-axis motion should instead be integrated as well in a self-consistent way. Yet, Figs. 20 and 21 already give an idea of what can happen to Titan’s orbit when the system reaches the unstable region E₁ and the vicinity of the singular point S₁. We see that its eccentricity and inclination can reach almost any value. In particular, Titan goes well below its Roche radius in Figs. 20 and 21.

Fig. 20 Long-term dynamics of Titan as it migrates away from Saturn. Saturn is assumed to remain locked in secular spin–orbit resonance with s8 at all time; we take its obliquity evolution at the centre of the resonance (red curve in Fig. 18). Titan evolves according to the secular Hamiltonian function given in Eq. (1) in which its semi-major axis is a slowly varying parameter. The top and bottom panels give two outcomes of the chaotic transitions obtained by using a slightly different migration rate. In all panels, the blue curve shows the location of the circular Laplace state P1, and the red curve shows the location of the eccentric Laplace state P.

Fig. 21 Same as Fig. 20, except that the integration is performed using the unaveraged equations of the restricted three-body problem including Saturn’s oblateness. The output is then digitally filtered to remove its short-period component. Due to the chaotic divergence of trajectories, we ran several simulations with different migration velocity and chose two of them that resemble the trajectories in Fig. 20.

Fig. 22 Numerical integration of Saturn’s spin-axis dynamics for Titan remaining at all time in its circular Laplace equilibrium as it migrates away. Titan’s migration parameter is set to b = 1 (see Eq. (47)) and each panel corresponds to a given value of Saturn’s polar moment of inertia (see the red labels in Fig. 23). Trajectories are shown in black. The pink bands are Saturn’s first-order secular spin–orbit resonances, and the blue hatched area is the region E1 where the Titan’s circular equilibrium is unstable (same as Fig. 18).

4.3 Saturn’s spin-axis

Using the Hamiltonian function in Eq. (27), our setting is similar to the migration simulations of Saillenfest et al. (2021b), except that both the semi-major axis and the inclination of Titan evolve over time as slowly varying parameters. The problem with this approach, and the reason why it has not been used in previous works, is that Titan’s inclination (when it is not in the close-in or far-away regime) depends on Saturn’s obliquity. Therefore, Saturn’s effective precession constant α′ depends on the obliquity in a complicated way (see Sect. 3.2), and we lose the Hamiltonian structure described by Eq. (27).

Yet, except in the vicinity of the singular point S₁, the dependence of α′ on the obliquity is weak. Therefore, at first level of approximation, the equations of motion obtained from the Hamiltonian in Eq. (27) are still valid, and the dependence of α′ on the obliquity can be added afterwards in the equations of motion, in order to account for its long-term drift. Including Titan’s Laplace plane inclination in α′ means that for any obliquity variation of Saturn, Titan instantly moves at the new equilibrium configuration. As explained above, this approximation is justified by the large separation between the two timescales, but it fails near S₁, where the inclination variations of Titan as a function of obliquity are extremely sharp (and discontinuous exactly at S₁). But as shown in Sect. 4.2, Titan is expected anyway to be destabilised before actually reaching S₁. Therefore, we stress that this model is not self-consistent; we use it here as a quick way to assess the relevance of our analytical predictions in Sect. 4.1, and to give a first qualitative picture of the different possible trajectories for Saturn. Apart from the obliquity dependence in α′, our model is the same as that of Saillenfest et al. (2021b): the orbit of Saturn evolves according to the full series of Laskar (1990), and Titan is made to migrate outwards according to the migration law in Eq. (47). Since Saturn’s polar moment of inertia and Titan’s migration rate are not well known, we perform a large number of simulations with parameters (b, λ) sampled in their uncertainty ranges. We refer to Saillenfest et al. (2021b) for a discussion about all parameters and their uncertainties. We note that because of our rigid rotation model, λ can be considered as an effective parameter that may slightly differ from what would be obtained using a refined model with differential rotation.

Figure 22 shows examples of trajectories for six different values of Saturn’s normalised polar moment of inertia λ. Contrary to Saillenfest et al. (2021b), we do not represent the trajectories as a function of Saturn’s precession constant because this constant is ill-defined near the singular point S₁ (see Sect. 3.2). In order to compare Fig. 22 with previous works, we stress that the initial point of the trajectory isthe same on each panel; what differs here is the locations of the resonances, which are slightly shifted from one value of λ to another. Werecognise the different types of trajectories described by Saillenfest et al. (2021b):

In panel a, Saturn is currently out of the resonance and it goes farther away as Titan migrates. The two crossings of the ν₅₁ resonance do not produce substantial obliquity variations.

In panel b, Saturn currently oscillates inside the resonance with a large amplitude. When the resonance width decreases, the adiabatic invariant cannot be conserved and Saturn escapes the resonance by crossing the separatrix. In this case, Saturn reaches a large obliquity, but Titan may remain stable anyway because it does not enter deep inside the unstable zone E₁ (hatched region). Eventually, Saturn crosses the s₈ resonance again when Titan passes in the far-satellite regime, producing an obliquity kick.

In panels c and d, Saturn currently oscillates closely around the resonance centre (with a minimum libration amplitude of 30° or so; see Ward & Hamilton 2004). As shown in previous works, this configuration is the most likely in a dynamical point of view, regardless of the actual mechanism that is responsible for Saturn’s resonance encounter with s₈ (Hamilton & Ward 2004; Boué et al. 2009; Vokrouhlický & Nesvorný 2015; Saillenfest et al. 2021b). In this case, we see that Saturn is able to get very close to the singular point S₁ before beingdestabilised, because the neighbouring resonances are thin and their overlap does not produce much chaos. Afterthe chaotic transition, Saturn can either be ejected from the resonance with an obliquity of about 90° or slightly more (panel c), or it can be recaptured at once when the resonance width increases again (panel d). However, as shown in Sect. 4.2, Titan is expected to be completely destabilised before reaching S₁, so the evolution of Saturn’s obliquity should remain frozen at some point as Titan is removed (collision or ejection). The questions of where this transition happens and what is the statistical outcome of the destabilisation would require a self-consistent numerical model; these questions are left for future works.

In panels e and f, Saturn did not reach yet the resonance with s₈ today. As discussed by Saillenfest et al. (2021b), this configuration would require a value of λ that is slightly out of its expected range (namely λ ≳ 0.241 while we expect λ ∈ [0.200, 0.240]). We include it here for completeness. As Saturn encounters the resonance, it can either cross it without being captured (panel e), in which case its obliquity suffers from a small kick, or it can be captured (panel f), in which case we end up with the same kind of evolution as in panel b.

As a summary of these numerical experiments, Fig. 23 shows the maximum obliquity reached by Saturn for Titan migrating from its current location a = a₀ up to its midpoint radius a = r_M. The figure shows the results obtained in a grid made of 250 values of b and 501 values of λ. We stress that, contrary to previous works, the integration duration is not the same for each run, but depends on b (see the top horizontal axis). The bottom dark-blue region in Fig. 23 corresponds to the cases where Saturn is out of the resonance today and goes farther away as Titan migrates. The large coloured region corresponds to values of the parameters that put Saturn inside the resonance today. It is narrower for larger b because large-amplitude librations are unstable if Titan’s migration is too fast. The top region (which is out of the expected range for λ) corresponds to cases where Saturn did not reach yet the resonance today but will in the future, resulting in a capture (coloured stripe) or not (dark background). See Fig. 16 of Saillenfest et al. (2021b) for a discussion about this striped pattern.

From Fig. 23, we conclude that Saturn gets extreme obliquities in a large region of the parameter space. If Saturn is inside the resonance today, then its obliquity reaches at least about 75°. This limit isrobust in spite of our simplified model because it is reached before encountering the unstable region (see Fig. 22b). The preferred parameter range for Saturn in previous studies (Boué et al. 2009; Vokrouhlický & Nesvorný 2015; Saillenfest et al. 2021b) is precisely the one producing the maximum obliquity increase (about 91° in our simplified model); it is also the one producing the maximum destabilisation for Titan (see Sect. 4.2). However, according tothe exact migration rate of Titan, the system may not reach the instability by the end of the Sun’s main sequence.

Fig. 23 Maximum obliquity reached by Saturn for Titan migrating from its current location up to a = rM. Titan is assumed to remain at all time in its circular Laplace equilibrium. The integration time (top horizontal axis) depends on Titan’s migration parameter b (bottom horizontal axis). The noisy grey curve is the 90° obliquity level. The red labels show the location of the six examples of trajectory presented in Fig. 22.

5 Summary and discussions

Titan is observed to migrate away from Saturn much faster than previously thought (Lainey et al. 2020). This migration is likely responsible for Saturn’s current axis tilt of 26.7° (Saillenfest et al. 2021a), in which case Saturn should still be trapped today in secular spin–orbit resonance with Neptune’s nodal precession mode s₈. Since Titan goes on migrating today, Saturn’s obliquity is expected to increase in the future (Saillenfest et al. 2021b). In this article, we investigated the final outcome of this mechanism, and the behaviour of the system when Titan will cross its Laplace radius.

The intricate nature of the dynamics requires the evolution of both the satellite’s orbit and the planet’s spin axis to be studied. Building on the work of Tremaine et al. (2009), we have shown that the three circular equilibria for the satellite (dubbed here Laplace states P₁, P₂, and P₃) are organised around two critical regions S₁ and S₂ in the parameter space, in which thepairs (P₁, P₃) or (P₂, P₃), respectively, are degenerate. In particular, S₁ is defined as the point where the planet’s obliquity is ε = 90° and the satellite’s semi-major axis is a = 2^1∕5r_L (where r_L is the Laplace radius defined by Tremaine et al. 2009).

We found that all three circular equilibria bifurcate to eccentric equilibrium configurations (noted , , , and ) in some regions of the parameter space. The location of all eccentric equilibria can be expressed with explicit parametric representations. The critical regions S₁ and S₂ both have an eccentric continuation, in which the pairs or , respectively, are degenerate.

Regular satellites like Titan form in their classical Laplace plane (i.e. at equilibrium P₁). As long as the equilibrium remains stable, they stay in its vicinity during their orbital migration, and their orbital inclination varies accordingly (see e.g. Tremaine et al. 2009). Using numerical integrations, we verified that this is indeed the case for Titan, even when taking into account its fast orbital expansion measured by Lainey et al. (2020).

As shown in previous works, satellites increase the mean spin-axis precession rate Ω₀ of their host planet, with a magnitude that depends on their orbital distance (see e.g. Boué & Laskar 2006). Consequently, if the planet is trapped in a secular spin–orbit resonance, any migration of its satellites is compensated by an obliquity change, so that Ω₀ is maintained fixed at the resonant value (usually called Cassini state 2). In other words, the planet and its satellite follow a level curve of Ω₀ in the plane (a, ε). When the satellite’s inclination is fully taken into account (i.e. without the usual close-in or far-away approximations), the singularity S₁ in its orbital dynamics is transferred to the spin-axis precession rate of the planet. As a result, if the satellite is massive enough, most level curves of Ω₀ converge towards the singularity S₁. This means that if the satellite reaches the vicinity of its Laplace radius during its migration, the obliquity of its host planet is driven towards ε = 90°. If the resonance is large, the planet can even go beyond this limit. Simplified analytical formulas give the conditions required for a full 90°-tilt. These conditions are met for Titan and Saturn and the s₈ resonance, and confirmed using numerical simulations.

In the vicinity of S₁, however, several kinds of instability are expected to happen. Firstly, S₁ lies at the border of a region in the parameter space where the circular Laplace state P₁ is unstable. As discussed in previous works (Tremaine et al. 2009; Tamayo et al. 2013), when a satellite crosses this border, it can transfer to the stable eccentric equilibrium P, but this equilibrium soon becomes unstable, too. At this point, numerical integrations show that Titan’s eccentricity and inclination completely destabilise, potentially allowing for the ejection of Titan or its collision on Saturn. In a hypothetical system in which the satellite would migrate inwards, the destabilisation is expected to be more violent, because in this case the system would smoothly reach the extreme vicinity of S₁, where the instability is strongest, before being destabilised.

Secondly, all neighbouring secular spin–orbit resonances converge towards S₁. The planet’s spin axis is therefore expected to reach a chaotic region at some point, produced by resonance overlap. In the case of Saturn, the neighbouring resonances are very thin, and numerical experiments show that the chaotic region is restricted to a very small region near S₁.

Thirdly, the widths of all secular spin–orbit resonances decrease in the vicinity of S₁. As a result, if the libration amplitude of the planet inside the resonance is too large, it can be ejected before actually reaching ε = 90°. Using a simplified numerical model, we performed a preliminary survey of Saturn’s behaviour and explored a large range for its poorly known moment of inertia. We found that Saturn can indeed be ejected from resonance with an obliquity ε ≳ 75°, depending on its libration amplitude inside the resonance. Regardless of the mechanism responsible for Saturn’s capture in secular spin–orbit resonance, previous studies show that Saturn is likely located deep inside the resonance today (Boué et al. 2009; Vokrouhlický & Nesvorný 2015; Saillenfest et al. 2021b). In this case, Saturn’s obliquity increase is maximum, and it may reach ε ≈ 91° in the future (provided that Titan is not ejected before; see above). Determining the statistical outcome of this double dynamical instability for Saturn and Titan is left for future works. It will require a self-consistent model coupling the orbital motions of Saturn and Titan, the spin-axis dynamics of Saturn torqued by the Sun and by Titan, and the multi-harmonic orbital precession of Saturn produced by the other planets of the Solar System.

Our results about Saturn and Titan rely on the assumption that Titan will go on migrating for gigayears in the future. If instead, Titan’s migration rate strongly drops before reaching the unstable zone (i.e. if Titan is released out of the tidal resonance-locking mechanism of Fuller et al. 2016) the system would remain frozen, with roughly constant obliquity for Saturn and fixed orbit for Titan. However, to our knowledge, there is no evidence showing that Titan’s migration would stop in the future (at least not before it becomes strongly unstable). The timescale required to tilt Saturn also plays an important role. According to the precise migration rate of Titan, Saturn and Titan should reach the instability region between a few gigayears and several tens of gigayears from now, provided that Titan goes on migrating as expected. Therefore, the evolution timescale may be too slow for the system to reach instability by the end of the Sun’s main sequence. Yet, even if it does not reach instability in time, Saturn would anyway get very large obliquity values, as already described in Fig. 16 of Saillenfest et al. (2021b).

Moreover, the mechanism described here is generic: it only requires a secular spin–orbit resonance and a substantially massive migrating moon. The structure of the resonance and its convergence towards ε = 90° near the moon’s Laplace radius are the same for any (exo)planet considered, as well as the mechanism of destabilisation. In particular, tilting Uranus on a gigayear timescale by a former satellite now removed by the instability is a promising mechanism that has not been explored yet (Boué & Laskar 2010; Rogoszinski & Hamilton 2020, 2021). This mechanism may also provide a dynamical explanation for the ‘super-puff’ exoplanets thought to possess a massive face-on ring (Akinsanmi et al. 2020; Piro & Vissapragada 2020): the destroyed satellite would both provide the tilting mechanism (for the obliquity to be near 90°) and the ring material. More work is now needed to assess the feasibility of this scenario to these specific systems.

Acknowledgements

We thank Ariane Courtot for her work that led us to a deeper understanding of the coupling between satellite and spin axis. M.S. also thanks Aurore Mazur for her contribution to the study of Titan’s destabilisation mechanism. We are grateful to the anonymous referee for having carefully read our manuscript and pointed out some ambiguity in terminology.

Appendix A Conversion formulas between the equator and ecliptic reference frames

The ecliptic reference frame is defined with the third axis perpendicular to the planet’s orbit; we write (I_C, ω_C, Ω_C) the Keplerian elements of the satellite measured in this frame, which are, respectively, the inclination, the argument of pericentre, and the longitude of the ascending node. The equator reference frame is defined with the third axis perpendicular to the planet’s equator; we write (I_Q, ω_Q, Ω_Q) the Keplerian elements of the satellite measured in this frame. Passing from one frame to the other corresponds to a rotation of ± ε around the mutual line of nodes of the two reference planes (i.e. the line joining both equinoxes), ε being the obliquity of the planet. We obtain (A.1)

and (A.2)

where δ_Q ≡ Ω_Q − Ω_⊙ and Ω_⊙ is the ascending node of the star measured along the equator of the planet. The first equation reverses as (A.3)

where δ_C ≡ Ω_C − Ω_P and Ω_P is the ascending node of the planet’s equator measured along the ecliptic. These three equations are all what we need to express the Hamiltonian function in Eq. (1) in terms of the equator or ecliptic coordinates only. The expressions obtained are given in Eqs. (3), (4) and (5).

Appendix B Circular Laplace equilibria

In Sect. 2, we describe the secular orbital dynamics of a massless satellite perturbed by the Sun and by the oblateness of its host planet at quadrupole order. We call the three kinds of circular equilibria the ‘Laplace states’ P₁, P₂, and P₃. Even though the equilibrium P₁ is the one that matters most for regular satellites, all equilibria can play a role in the dynamics when the satellite becomes unstable and explores wide regions in the phase space (see Sect. 4). For this reason, we recall here the stability properties of all three circular equilibria.

Appendix B.1 Circular equilibrium P₁

P₁ is stable to inclination variations (see Fig. 1). For δ_Q = 0, the equatorial inclination I_Q1 of the satellite at P₁ is given by Eq. (14). Inclination oscillations around the Laplace state P₁ have frequency ξ₁ given in Eq. (16). For a ε = 0° or 180°, it simplifies into Eq. (17). Likewise, the value obtained for ε = 90° is (B.1)

and we note that it tends to 0 at the singular point S₁ (i.e. a = r_M), since we retrieve the eigenfrequency along the degenerate stable circle . For a very close satellite, ξ₁ becomes independent of the obliquity: (B.2)

where is the satellite’s mean motion. This is the usual nodal precession frequency produced by the oblateness of the central body in absence of other perturber (see e.g. Murray & Dermott 1999).

P₁ is stable to eccentricity variations, except in the region E₁ described in Sect. 2.3. Eccentricity oscillations around P₁ have frequency μ₁ given in Eq. (19). For a zero-obliquity planet, Eq. (19) simplifies to (B.3)

which is exactly equal to the inclination eigenfrequency given by Eq. (17). Therefore, for a planet with a near-zero obliquity, the oscillations of the inclination and eccentricity of a satellite around the classical Laplace state P₁ have the same frequencies. The value obtained for ε = 90° is (B.4)

We see that does not go to zero when a → r_M and that it has two different limits; this means that going through the singularity S₁ produces a discontinuity of . This could have been expected since P₁ and P₃ are inverted(i.e. there is a discontinuity in their location as well, jumping by 90° in inclination; see Sect. 2.2).

For a very close satellite, μ₁ becomes independent of the obliquity: (B.5)

which is the same value of in Eq. (B.2). This is the usual pericentre precession frequency produced by the oblateness of the central body in absence of other perturber (see e.g. Murray & Dermott 1999).

The eccentricity eigenfrequency along the degenerate stable circle (see Fig. 2b) is given by (B.6)

The value of depends on the inclination of the satellite along the degenerate circle. We retrieve the discontinuous limits of Eq. (B.4) when a → r_M by putting I_Q = 0° or 90° in Eq. (B.6). Along the degenerate circle, the satellite is stable only if cos²I_Q ≤ 1∕5, that is, if 63° ≲ I_Q ≲ 117°. This corresponds to the region where the J₂-induced precession of ω_Q is negative.

Appendix B.2 Circular equilibrium P₂

P₂ is stable to inclination variations, and located at δ_Q = ± π∕2 and I_Q = 90° (see Fig. 1). As described by Tremaine et al. (2009), P₂ corresponds to an orbit that is perpendicular to both the equator and the ecliptic planes. The frequency of small-amplitude inclination oscillations around P₂ is given by (B.7)

As expected, the value of is always negative, and it tends to zero in the region S₂ of the parameter space (i.e. ε = 0° or 180°), since we retrieve the eigenfrequency along the degenerate stable circle .

The frequency of small-amplitude eccentricity oscillations around P₂ is given by (B.8)

We deduce that P₂ is stable to eccentricity variations only for close enough satellites that verify a ≤ r₃, where we define (B.9)

The region where P₂ is unstable to eccentricity variations is therefore (B.10)

The eccentricity eigenfrequency along the degenerate stable circle (see Fig. 2c) is also given by Eq. (B.8). Hence, as for the Laplace state P₂, it is only stable for close-enough satellites.

Appendix B.3 Circular equilibrium P₃

P₃ is unstable to inclination variations (see Fig. 1). The squared eigenvalue of the linearised problem at P₃ is given by Eq. (16) where I_Q1 must be replaced by I_Q3, which is equal to I_Q1 ± 90°. As expected, whenever P₃ exists, and tends to 0 at the singular regions S₁ and S₂, where we retrieve the eigenfrequencies and along the degenerate stable circles and .

Even though P₃ is unstable to inclination variations anywhere it exists in the parameter space, this is not the case for eccentricity variations (Tremaine et al. 2009). As illustrated in Figs. 5 and B.1, P₃ is stable to eccentricity variations for close-enough satellites and in a region resembling a water drop. The squared eigenvalue of the linearised problem at P₃ is given by Eq. (19) where I_Q1 must be replaced by I_Q3 = I_Q1 ± 90°. We call E₃ the unstable region where . Its border can be expressed piecewise as (B.11)

for r₃ ≤ a ≤ r_L, and Eq. (22) for r_L ≤ a ≤ r_M, recalling that r_L is the traditional ‘Laplace radius’ defined in Eq. (7), and r₃ is given in Eq. (B.9).

Fig. B.1 Region E3 of the parameter space where the Laplace state P3 is unstable to eccentricity variations. The colours and labels have the same meaning as in Fig. 8, and the axes have the same scale. The border of E3 is obtained using the closed-form expression in Eqs. (B.11) and (22). The drop-like portion of the contour reaches its maximum width at where the obliquity on the boundary is .

Appendix C Eccentric Laplace equilibria

The Hamiltonian function in Eq. (1) describes the secular orbital dynamics of a massless satellite perturbed by the star and by the oblateness of its host planet at quadrupole order. In Sect. 2.2, the circular equilibrium configurations, called Laplace states P₁, P₂, and P₃ are described. Wherever they are defined, all three of them are stable to eccentricity variations, except in the regions E₁, E₂, and E₃, respectively,described in Sect. 2.3 and illustrated in Fig. 5.

Tremaine et al. (2009) have shown that the system also admits eccentric equilibria, that is, configurations in which the satellite has a frozen eccentric orbit. Two of these configurations bifurcate away from P₁ and P₂ where these equilibria become unstable (i.e. at the borders of the E₁ and E₂ regions).

In this section, we recall these results and go further in the analytical characterisation of the eccentric equilibria, allowing one to compute them more easily and in a form that is more directly usable in the context of our work. Furthermore, we show that the remaining eccentric equilibria described by Tremaine et al. (2009) bifurcate away from P₃ at the border of region E₃; therefore, all eccentric equilibria (that we call below P, P, P, and P) emerge from a bifurcation of a circular Laplace equilibrium.

Appendix C.1 Eccentric equilibrium P

The eccentric equilibrium P corresponds to one of the configurations called ‘eccentric coplanar–coplanar Laplace equilibria’ by Tremaine et al. (2009). As such, the orbital angles of the satellite at P are ω_Q = π∕2 mod π and δ_Q = 0. As before, because of the symmetries of the secular problem, each eccentric equilibrium point has a twin obtained by the transformation (δ_Q, I_Q) → (π + δ_Q, π − I_Q) that corresponds to the same Laplace state with reversed orbital motion (see e.g. Fig. 1). The equatorial inclination of the satellite at P is given in Eq. (24), and its eccentricity e is obtained by solving the equation (C.1)

By injecting the expression of into Eq. (C.1), we obtain an equation that only depends on e and on the parameters of the problem, a∕r_M and ε. If we solve it for e → 0, we obtain again the boundary of the E₁ region defined at Eqs. (21) and (22); this shows that P indeed bifurcates away from P₁ where the latter becomes unstable.

When trying to solve Eq. (C.1) for as a function of and ε, we obtain an intractable polynomial of order 20. However, there exists a straightforward analytical expression for ε as a function of u and η. Even though this is not exactly what we are looking for, this expression can be used to visualise the solutions in a graphical form and to set up efficient root-finding algorithms to reverse the expression for η. We define the obliquity solutions ε₊ and ε₋ as (C.2)

where (C.3)

and (C.4)

The values of cos²ε₊ and cos²ε₋ are shown in Fig. C.1. They give the location of the eccentric equilibrium P in the whole space of parameters. We see that solutions exist for any eccentricity value. However, they lie in restricted regions of the (η, u) space, whose boundaries can be defined piecewise by analytical formulas of the form u = f(η). These formulas are directly labelled on the figure, except u₊ and u₋ whose expressions are: (C.5)

Along the curves u₊(η) and u₋(η), the variable Y cancels in Eq. (C.3), which means that cos²ε₊ and cos²ε₋ have the same value (double root). Along the remaining boundary curves in Fig. C.1 (green, orange, and magenta curves), the solution cosε represented is equal to zero.

Instead of projecting the solutions in the (η, u) plane as in Fig. C.1, the closed-form expression in Eq. (C.2) can be used to plot the eccentricity of the satellite at equilibrium P as a function of the parameters in a parametric form. We obtain a three-dimensional surface, as illustrated in Fig. C.2. As expected, the bottom section of this surface (i.e. at e = 0) coincides with the border of the region E₁ where the circular equilibrium P₁ is unstable (see Figs. 5 and 8). In the main text, Fig. 9 shows the value of the eccentricity at P as a colour scale, as well as the value of the inclination obtained using Eq. (24). In order to draw this figure, the tube-like portion of the three-dimensional surface has been cut off.

Along the three-dimensional curve defined in Eq. (26) and depicted by a magenta curve in Figs. C.1 and C.2, the three-dimensional surface has a cusp. As discussed in the main text, the curve is the eccentric continuation of the singular point S₁ described in Sect. 2.2, where the equilibrium point P₁ is singular. We note that pierces the three-dimensional surface in Fig. C.2 and reappears on the other side (this corresponds to the dotted portion in Fig. C.1).

We did not find a closed-form expression for the boundary dividing the regions of the parameter space where P is stable from the regions where it is unstable. However, the stability nature of P as a function of the parameters can be determined numerically: the equilibrium is stable wherever the linearised system has no eigenvalue with positive real part. Figure C.3 shows the stable and unstable regions projected on the three-dimensional surface of P. Sections of this surface can be seen in Fig. 6 of Tremaine et al. (2009). In the stable regions, the oscillation frequencies of the satellite are illustrated in Fig. 11.

For completeness, Fig. C.4 illustrates the stability nature of P in the space (r_M∕a, ε, e²), where the full three-dimensional shape can be visualised, including the region where a →∞. We see that there is a small additional region where P is stable, for a large semi-major axis and an eccentricity very close to 1 (see the thin grey border on the right side of the figure). This thin stable region also appears in Fig. 5 of Tremaine et al. (2009) as the blue points lying at the tip of the middle triangle⁶.

Fig. C.1 Location of the eccentric equilibrium P as a function of the parameters. The colour shades show the two closed-form obliquity solutions given by Eq. (C.2). Some levels are highlighted in white. The definition regions of the solutions are bounded by curves of the form u = f(η), as labelled on the figures (coloured curves). Among those, u+ and u− are given in Eq. (C.5). The junction points (black labels) are given in Table C.1.

Fig. C.2 Eccentricity e of the satellite at equilibrium P as a function of the parameters. The solutions lie on a three-dimensional surface shown here from two viewing angles. This figure is obtained by stitching together the patches represented in Fig. C.1. For better comparison, the dividing curves are drawn on the surface using the same colour code. The top tube-like portion of the surface has been cut; as shown in Fig. C.1, it extends up to a → ∞ (i.e. u → 0). The equilibrium P is singular along the magenta curve, called in the text.

Appendix C.2 Eccentric equilibrium P

The eccentric equilibrium P corresponds to the configuration called ‘eccentric orthogonal–coplanar Laplace equilibrium’ by Tremaine et al. (2009). As such, the orbital inclination of the satellite at P is I_Q = I_C = 90°, while its orbital angles are ω_C = 0 mod π and δ_Q = π∕2 (or 3π∕2 for its twin equilibrium with reversed orbital motion). The eccentricity of the satellite at P is given by the relation (C.6)

where and . If we solve it for e → 0, we obtain again the critical distance a = r₃ defined in Eq. (B.9); this shows that P indeed bifurcates away from P₂ where the latter becomes unstable. Hence, P exists for a ≥ r₃ and the eccentricity of the satellite at P is (C.7)

The bifurcation of P from P₂ can be visualised as a three-dimensional surface, illustrated in Fig. C.5.

Fig. C.3 Stability of the equilibrium point P as a function of the parameters. The stable and unstable regions are painted in grey and red, respectively, and they are projected on the three-dimensional surface describing the eccentricity of the satellite. As in Fig. C.2, the surface is seen from two viewing angles and the top tube-like portion has been cut for better readability.

Fig. C.4 Same as Fig. C.3 but replacing the coordinates a∕rM and e by rM ∕a and e2, respectively.Contrary to previous figures, the top tube-like portion of the surface is seen in its totality.

As shown by Tremaine et al. (2009), P is stable wherever it exists. The eigenfrequencies of the linearised system are given by (C.8)

and (C.9)

The eccentricity and inclination of the satellite are coupled in the vicinity of P, so they their oscillation spectra both contain the two frequencies and . At the transition radius a = r₃, we note that and are equal to the oscillation frequencies ξ₂ and μ₂ around P₂ given in Eqs. (B.7) and (B.8), confirming that the eccentric equilibrium P bifurcates away from the circular equilibrium P₂.

Fig. C.5 Eccentricity e of the satellite at equilibrium P as a function of the parameters. The eccentricity is given by Eq. (C.7); it only depends on the semi-major axis a∕rM of the satellite. The equilibrium is stable wherever it exists.

Appendix C.3 Eccentric equilibria P and P

As explained in Sect. 2, the circular equilibrium P₃ is unstable to eccentricity variations in the region E₃ of the parameter space. In its illustrations in Figs. 5 and B.1, we see that the border of E₃ can be divided into two components: a small drop-like boundary for a > r_L (whose expression is given by Eq. (22)), and a V-shaped boundary for a < r_L spanning all values of obliquity ε from 0 to 180° (and whose expression is given by Eq. (B.11)). Below, we show that along these two boundaries, the circular equilibrium P₃ bifurcates into two distinct eccentric equilibria. We call them P and P, respectively, for the eccentric equilibria bifurcating away from the drop-like boundary and from the V-shaped boundary.

The eccentric equilibrium P correspondsto one of the configurations called ‘eccentric coplanar–coplanar Laplace equilibria’ by Tremaine et al. (2009). As such, the orbital angles of the satellite at P are ω_Q = π∕2 mod π and δ_Q = 0 (or δ_Q = π for the twinequilibrium with reversed orbital motion). The equatorial inclination of the satellite at P can be written as (C.10)

where v depends on the eccentricity at equilibrium, as given in Eq. (25). The inclination in Eq. (C.10) has the same form as I_Q3 in Eq. (14), but where is replaced by v. The two definitions coincide for e = 0.

The eccentricity e of the satellite at P is obtained by solving the equation (C.1), but where must be replaced by . If we solve it for e → 0, we obtain again the drop-like boundary of the E₃ region; this shows that P indeed bifurcates away from P₃ along this boundary. The resolution of the equation in the general case is very similar to what is presented in Appendix C.1 for P: we obtain a closed-form analytical expression for ε as a function of u and η, which is again given by Eqs. (C.2), (C.3), and (C.4), but with a differing domain of definition.

The values of cos²ε₊ and cos²ε₋ for the eccentric equilibrium P are shown in Fig. C.6. They give its location in the whole space of parameters. Contrary to P, we see that solutions do not exist for all eccentricity values: there is a gap between η_g and η_f where P does not exist (see the corresponding values of eccentricity in Table C.1).

The expression in Eq. (C.2) can be used to plot the eccentricity of the satellite at equilibrium P in a parametric form. We obtain the three-dimensional surface illustrated in Fig. C.7. As expected, the bottom section of this surface (i.e. at e = 0) coincides with the drop-like boundary of the region E₃ where the circular equilibrium P₃ is unstable to eccentricity variations (see Figs. 5 and B.1). As mentioned in Appendix C.1, the eccentric equilibrium P is singular along the three-dimensional curve given in Eq. (26), where is undefined. It corresponds to a degenerate equilibrium which results from the merging of P and P. Indeed, we note that is a contact region between the two three-dimensional surfaces shown in Figs. C.2 and C.7 (see the magenta curve).

Figure C.8 highlights the regions of the parameter space where P is stable, projected on the three-dimensional surface. Interestingly, a small stable region exists on the top portion of the three-dimensional surface. This stable region is visible in Fig. 5 by Tremaine et al. (2009) as the blue points in the uppermost triangle. Therefore, even though the circular equilibrium P₃ described in Sect. 2 is unstable to inclination variations in the whole space of parameters, this is not everywhere the case for its eccentric counterpart P. Yet, Fig. C.8 shows that this stable region is disconnected from the low-eccentricity portion of the three-dimensional surface, which means that a regular satellite starting with a roughly circular orbit cannot reach it by a smooth change of parameters occurring after its formation. For completeness, Fig. C.9 illustrates the stability nature of P in the space (r_M∕a, ε, e²), where the full three-dimensional shape can be visualised, including the region where a →∞. No additional stable region can be seen.

The eccentric equilibrium P corresponds to the configuration called ‘eccentric coplanar–orthogonal Laplace equilibrium’ by Tremaine et al. (2009). As such, the orbital angles of the satellite at P are ω_Q = 0 mod π and δ_Q = 0 (or δ_Q = π for the twinequilibrium with reversed orbital motion). The equatorial inclination of the satellite at P can be written as (C.11)

where we define (C.12)

in which e is the satellite’s eccentricity at equilibrium. The inclination in Eq. (C.11) has the same form as I_Q3 in Eq. (14), but where is replaced by w. The two definitions coincide for e = 0.

The eccentricity e of the satellite at P is given by the implicit equation (C.13)

that must be inverted to obtain w (and thus e) as a function of ε. We recognise the V-shaped boundary of the region E₃ (see Eq. (B.11)), but where u is replaced by w. This shows that P indeed bifurcates away from P₃ along this boundary. Through the inversion of Eq. (C.13), the inclination is a function of ε only and it does not depend on the orbital distance a∕r_M. The eccentric equilibrium P is unstable wherever it exists, except for ε = 0° or 180°, where its unstable mode tends to zero. Indeed, the three-dimensional domain (C.14)

is the eccentric continuation of the singular region S₂ described in Sect. 2, where P₂ and P₃ merge with the separatrix. Similarly, results from the merging of P and P, which degenerate into an equilibrium region . The bifurcation of P from P₃ and the two curves that define can be visualised in Fig. C.10. We see that is the contact region between the two three-dimensional surfaces shown in Figs. C.5 and C.10.

Inside the region , the degenerate equilibrium is defined by an eccentricity e given by Eq. (C.7), an inclination I_Q = 90°, an argument of pericentre ω_Q = 0 mod π, and any value for the longitude of node δ_Q. In the vicinity of , the eigenfrequencies of the system are given by , and by Eq. (C.9) for .

Fig. C.6 Same as Fig. C.1, but showing the location of the eccentric equilibrium P as a function of the parameters. On the top panel, the definition region of cos2 ε+ is very narrow, comprised between the curves u = η3(5 − 4η2) and u = u−(η).

Fig. C.7 Eccentricity e of the satellite at equilibrium P as a functionof the parameters. This figure is obtained by stitching together the patches represented in Fig. C.6. The top tube-like portion of the surface has been cut; as shown in Fig. C.6, it extends up to a → ∞ (i.e. u → 0). The equilibrium P is singular along the magenta curve, called in the text.

Fig. C.8 Stability of the equilibrium point P as a function of the parameters. The stable and unstable regions are painted in grey and red, respectively, and they are projected on the three-dimensional surface describing the eccentricity of the satellite. As in Fig. C.7, the surface is seen from two viewing angles and the top tube-like portion has been cut for better readability.

Fig. C.9 Same as Fig. C.8 but replacing the coordinates a∕rM and e by rM ∕a and e2, respectively.Contrary to previous figures, the top tube-like portion of the surface is seen in its totality.

Fig. C.10 Eccentricity e of the satellite at equilibrium P as a function of the parameters. The colour gives its stability: grey for stable and red for unstable. P is unstable in the whole parameter space except along , which results from the merging of P and P into a degenerate equilibrium.

Appendix D Self-consistent formula for the spin-axis precession rate

In Eq. (42), we give a simplified expression for the spin-axis precession rate of an oblate planet affected by a satellite. Taking advantage of the simplicity of this expression, many properties of the precession rate are given by closed-form analytical formulas, such as the location and magnitude of the maximum rate, or the shape of its level curves as a function of the parameters (see Sect. 3.2). However, as stressed in the main text, Eq. (42) is not self-consistent. The satellite is considered massless to compute its equilibrium inclination (Laplace state P₁), but massive when it comes to its influence on the planet’s spin axis. Then, the spin axis is taken as fixed when studying the satellite’s dynamics, but moving when taking into account secular spin–orbit resonances. Both these approximations are expected to be accurate for satellites with a small enough mass, because in this case the satellite’s orbit evolves on a much shorter timescale than the planet’s spin axis, and only the very long-term component of the satellite’s dynamics (mean secular motion) affects the planet’s obliquity. Yet, the validity range of our approximations are not known a priori, and since we explore a large range of mass parameters (see Fig. 17), a comparison with self-consistent models would be useful.

For this purpose, we use the model of Boué & Laskar (2006), which describes the orbital dynamics of an oblate planet and of its (massive) satellite, together with the spin-axis dynamics of the planet, considering the fully coupled equations of motion at quadrupolar order. Since the satellite is nonetheless considered to be less massive than its host planet, high-order terms in δ are neglected, where δ = m∕(m + M) with our notations. Even though the model of Boué & Laskar (2006) stands for any eccentricity of the satellite, it is averaged over the satellite’s argument of pericentre, which means that it cannot be used when the satellite oscillates around one of the stable eccentric equilibria described in Sect. 2.4. Under the approximation that the total angular momentum of the system is contained in the orbital motion of the planet (which is almost exactly verified in practice), Boué & Laskar (2006) give an analytical expression for the time-evolution of the satellite’s orbital pole and the planet’s spin axis (see their Eq. 133). Due to its high level of generality, this expression is quite complicated and cumbersome to use. However, it shows that if the satellite is placed on its circular Laplace plane, then the planet and its satellite rigidly precess as a whole at the frequency (D.1)

where (D.2)

in which (D.3)

are constant, and the physical parameters of the system are contained in the coefficients (D.4)

and (D.5)

In this self-consistent model, the equatorial inclination I_Q of the satellite at its circular Laplace equilibrium is obtained by cancelling the nutation amplitude of the solution (coefficient 𝔰 in Eq. 133 of Boué & Laskar 2006). This amounts to solving for I_Q the equation (D.6)

where Ω depends itself on I_Q through Eq. (D.1). For a given set of parameters, this equation can easily be solved numerically. For small masses, one can check that we retrieve very accurately the inclination given by the classical massless-case formula in Eq. (13).

Figure D.1 shows the precession period obtained for the parameters of Saturn and Titan. The mass of Titan is varied between each panel so as to produce the same mass parameter η as in Fig. 17; all other physical parameters are left unchanged. Because physical parameters are deeply entangled in the self-consistent model of Boué & Laskar (2006), the precession rate in Eq. (D.1) cannot be expressed in terms of a few macro-parameters, as done in Sect. 3.2. As a result, Fig. D.1 is specific to the parameters of Titan and Saturn, whereas Fig. 17 is generic.

For small satellite masses, Figs. 17 and D.1 give undistinguishable results (panels a and b). For a mass of the order of Titan’s (panel c), small differences in magnitude are noticeable near a = r_M, as shown by the slight displacements of the level curves labelled 0.2 and 0.17. For a satellite about ten times as massive as Titan (panel d), the ridge line is slightly shifted right as compared to our simplified model, but its magnitude is still very similar to that of Fig. 17. For even larger masses (not shown), the ridge line is further shifted right, and the contribution of the satellite to λ becomes non-negligible (see Eq. 39). However, the overall structure of the dynamics remains the same, with a singular point at ε = 90° towards which the level curves converge.

From this comparison, we conclude that the simplified model presented in Sect. 3.2 provides a very good qualitative description of the dynamics, even though the exact magnitude of the precession rate can differ somewhat for large satellite-to-planet mass ratios. For the Earth-Moon system, whose mass ratio is as large as about 0.01, a self-consistent precession model would be required to get precise quantitative results.

Fig. D.1 Same as Fig. 17, but using the self-consistent precession rate of Boué & Laskar (2006) given in Eq. (D.1). The physical parameters are those of Titan and Saturn, except that the mass of Titan is adjusted to produce the same mass parameters η as in Fig. 17. In order to ease the comparison, the ridge line and singular point (in red) are those of the simplified model, and the labelled level curves are the same as in Fig. 17 (except the level 0.021 in panel d which is added here).

References

All Tables

All Figures

Fig. 1 Level curves of the Hamiltonian function for a circularorbit. The parameters are a∕rM = 1.1 and ε = 40°. The separatrix is shown by a thicker black curve. The coloured dots represent the three kinds of equilibrium points (‘Laplace states’), labelled as in the text. A dark colour is used for points P1 and P3 lying at δQ = 0 and for P2 lying at δQ = π∕2. A light colour is used for the symmetric equilibrium point that corresponds to the same Laplace state with reversed orbital motion. Figure 2a shows the same phase portrait plotted on the sphere.

In the text

Fig. 2 Level curves of the Hamiltonian function with e = 0 plotted on the sphere. The z-axis is along the spin axis of the planet. The x-axis is along the intersection of the equatorial and ecliptic planes (i.e. the line joining both equinoxes of the planet) and directed towards the ascending node of the star. The equator xy plane and the ecliptic plane are highlighted by the two outer grey circles. A point on the sphere represents the tip of the orbital angular momentum of the satellite, which has coordinates (x, y, z) = (sinIQ sinδQ, −sinIQ cosδQ, cos IQ). The colour code is the same as in Fig. 1. Panel a: same parameters as Fig. 1. Panel b: parameter region S1, defined by a∕rM = 1 and ε = 90°. The magenta curve is made of an infinity of stable equilibria resulting from the merging of P1 with the separatrix emerging from P3. Panel c: parameter region S2, defined by a∕rM > 0 and ε = 0° (or 180°). The brown curve is made of an infinity of stable equilibria resulting from the merging of P2 with the separatrix emerging from P3.

In the text

	Fig. 3 Location of the Laplace states as a function of a∕rM. Each panel features a fixed value of ε (see labels). The level IQ = ε is shown by a horizontal dotted line. The colour code is the same as in Fig. 1. The black vertical line in the central panel shows the degenerate equilibrium circle produced by the merging of P1 and P3.
In the text

	Fig. 4 Location of the Laplace states as a function of ε. Each panel features a fixed value of a∕rM (from top to bottom: 0.01, 0.95, 1, 1.05, and 20, respectively). Colours have the same meaning as in Fig. 3.
In the text

Fig. 5 Summary of the regions of the parameter space governing the dynamics of a massless satellite on a nearly circular orbit. The phase space features three equilibrium points (‘Laplace states’), among which P1 and P2 are stable to inclination variations and P3 is unstable. At point S1 of the parameter space, P1 and P3 degenerate into the equilibrium circle that is stableto inclination variations. In region S2 of the parameter space, P2 and P3 degenerate into the equilibrium circle that is stableto inclination variations. In regions E1, E2, and E3, respectively, P1, P2, and P3 are unstable to eccentricity variations. In region S2 ∪E2, is unstable to eccentricity variations. At point S1, is stable to eccentricity variations provided that the equatorial inclination of the satellite verifies cos2 IQ ≤ 1∕5.

In the text

Fig. 6 Inclination of the Laplace state P1 as a function of the parameters (see Eq. (14)). Some level curves are labelled in black. The critical radii and rM and rF are represented in white. At a = rM, P1 lies exactly halfway between the equator and the ecliptic (i.e. IQ = ε∕2 for a prograde spin). The abrupt transition of IQ from 0° to 180° at ε = 90° is a coordinate singularity. The S1 point is a real singularity (see text).

In the text

	Fig. 7 Period of small inclination oscillations around the Laplace state P1 as a function of the parameters. The period is given by Eq. (16) and represented here in unit of the characteristic period τ = 2π∕κ. Some level curves are highlighted and labelled in black. The singular point S1 is indicated in red.
In the text

Fig. 8 Region E1 of the parameter space where the Laplace state P1 is unstable to eccentricity variations. The background colour shows the normalised value of , with blue for negative values (stable) and red for positive values (unstable). The border of E1 is highlighted with a black contour, obtained using the closed-form expression in Eqs. (21) and (22). Region E1 is entirely contained inside the dashed brown lines, whose left and right limits r1 and r2 are given in Eq. (20), and whose bottom and top limits are given by Eq. (23). The Laplace state P1 is singular atpoint S1, in which is discontinuous.

In the text

	Fig. 9 Eccentricity and inclination of the satellite at the eccentric equilibrium P. The three-dimensional surface of equilibrium has been cut along the grey line, as shown in Fig. 10. In the bottom panel, some level curves are plotted in white (see labels). Noting , the extremity of the grey cutting line have coordinates and at the top and bottom points, and and ε = 90° at the middle.
In the text

	Fig. 10 Eccentricity at the eccentric equilibrium P seen as a three-dimensional surface. The colour is the same as in Fig. 9, top panel. In Fig. 9, the top tube-like portion of the surface has been cut off along the grey line. As detailed in Appendix C, the top portion extends up to a → ∞, where it tends to e = 1. Sections of this surface can be seen in Fig. 6 of Tremaine et al. (2009).
In the text

Fig. 11 Period of small oscillations around P as a function of the parameters. There are two oscillation modes, shown in the two panels. In the hatched area, the equilibrium is unstable. The singular points S1 and S are labelled, and the cutting line of the three-dimensional surface is shown by a grey line (same as Fig. 9). Three-dimensional figures are provided in Appendix C, where the geometry of the stable regions is clearer.

In the text

Fig. 12 Spin-axis dynamics in the vicinity of a first-order secular spin–orbit resonance. Some level curves of the Hamiltonian function in Eq. (34) are represented in black or in red, according to whether they are outside or inside the resonance, respectively. The separatrix is shown by a thick black curve. The constant parameters are γ = 0.75 and β = 0.03. The equilibrium points (‘Cassini states’) are represented by coloured dots. Figure 13b shows the same phase portrait plotted on the sphere.

In the text

Fig. 13 Level curves of the Hamiltonian plotted on the sphere. The planet’s orbit lies in the xy-plane. The obliquity ε is the tilt from the z-axis, and the resonance angle σ is the polar angle measured in the xy-plane. The colour code is the same as in Fig. 12. Panel a: geometry for γ2∕3 + β2∕3 > 1. Panel b: geometry for γ2∕3 + β2∕3 < 1 (same parameters as in Fig. 12). Panel c: geometry for γ and β tending to zero.

In the text

	Fig. 14 Bifurcation diagram of the system as a function of 1∕γ, which is proportional to α. The second parameter used in this example is β = 0.03γ. The Cassini states are labelled with the same colour code as previous figures.
In the text

	Fig. 15 Spin-axis precession frequency of a planet orbited by a single satellite. The frequency is given by the simplified expression in Eq. (42) and represented here in unit of the characteristic frequency p, for different obliquity values (see labels). In this example, the mass parameter of the satellite is set to η = 10, which is close to Titan’s value (see Table 1).
In the text

	Fig. 16 Maximum spin-axis precession frequency of the planet as a function of its obliquity. The curve is obtained by injecting a = rF from Eq. (15) in Eq. (42). The extreme values Ωmax and ΩS are given in Eqs. (44) and (45). For ε = 90°, the value at the maximum is undefined.
In the text

Fig. 17 Spin-axis precession period of a planet orbited by a single satellite. The period is given by Eq. (42) and represented here in unit of the characteristic period T = 2π∕p. For ε >90°, the spin-axis precession frequency is negative. Each panel corresponds to a given value of the satellite’s mass parameter η (see Eq. (43)), as labelled in the top right corner. Some level curves are highlighted in black. The singular point S1 and the ridge line a = rF are indicated in red.

In the text

Fig. 18 Location and width of every first-order secular spin–orbit resonance reachable by Saturn over the migration of Titan, with amplitudes down to 10−8. Each resonant angle is of the form σj = ψ + ϕj where ϕj has frequency νj labelled on the graph according to its index in the orbital series (see Table 2). For a given value of Titan’s semi-major axis, the interval of obliquity enclosed by the separatrix is shown in pink. The centre of the resonance (i.e. Cassini state 2) with ϕ3 is shown by a red line and the upper and lower separatrices are highlighted in blue. The ridge line a = rF dividing the close-in and far-away satellite regimes is represented by a dashed blue line. The region E1 where the satellite’s Laplace state P1 is unstable isshown in black. The singular point S1 is shown by a red dot.

In the text

Fig. 19 Evolution of Saturn’s spin axis and Titan’s orbit while following the centre of the secular spin–orbit resonance with s8. We use blue for Saturn’s spin axis and red for Titan’s orbit. Each panel is described directly on the graph when the legend is not self-explanatory. In the right column, time is shown both in normalised units (left vertical axis) and in physical units (right vertical axis). The conversion factors are given in Table 1. All curves are obtained through the analytical formulas described in the text, assuming that Saturn closely oscillates near Cassini state 2 and Titan closely oscillates near Laplace state 1. In the top left panel, the green curve shows the location of the resonance centre according to the formula used by Saillenfest et al. (2021b) valid for a close satellite. In the bottom right panel, the black curve in the unstable region corresponds to the period needed for the eccentricity to be multiplied by exp (2π) ≈ 535.

In the text

Fig. 20 Long-term dynamics of Titan as it migrates away from Saturn. Saturn is assumed to remain locked in secular spin–orbit resonance with s8 at all time; we take its obliquity evolution at the centre of the resonance (red curve in Fig. 18). Titan evolves according to the secular Hamiltonian function given in Eq. (1) in which its semi-major axis is a slowly varying parameter. The top and bottom panels give two outcomes of the chaotic transitions obtained by using a slightly different migration rate. In all panels, the blue curve shows the location of the circular Laplace state P1, and the red curve shows the location of the eccentric Laplace state P.

In the text

Fig. 21 Same as Fig. 20, except that the integration is performed using the unaveraged equations of the restricted three-body problem including Saturn’s oblateness. The output is then digitally filtered to remove its short-period component. Due to the chaotic divergence of trajectories, we ran several simulations with different migration velocity and chose two of them that resemble the trajectories in Fig. 20.

In the text

Fig. 22 Numerical integration of Saturn’s spin-axis dynamics for Titan remaining at all time in its circular Laplace equilibrium as it migrates away. Titan’s migration parameter is set to b = 1 (see Eq. (47)) and each panel corresponds to a given value of Saturn’s polar moment of inertia (see the red labels in Fig. 23). Trajectories are shown in black. The pink bands are Saturn’s first-order secular spin–orbit resonances, and the blue hatched area is the region E1 where the Titan’s circular equilibrium is unstable (same as Fig. 18).

In the text

Fig. 23 Maximum obliquity reached by Saturn for Titan migrating from its current location up to a = rM. Titan is assumed to remain at all time in its circular Laplace equilibrium. The integration time (top horizontal axis) depends on Titan’s migration parameter b (bottom horizontal axis). The noisy grey curve is the 90° obliquity level. The red labels show the location of the six examples of trajectory presented in Fig. 22.

In the text

	Fig. B.1 Region E3 of the parameter space where the Laplace state P3 is unstable to eccentricity variations. The colours and labels have the same meaning as in Fig. 8, and the axes have the same scale. The border of E3 is obtained using the closed-form expression in Eqs. (B.11) and (22). The drop-like portion of the contour reaches its maximum width at where the obliquity on the boundary is .
In the text

Fig. C.1 Location of the eccentric equilibrium P as a function of the parameters. The colour shades show the two closed-form obliquity solutions given by Eq. (C.2). Some levels are highlighted in white. The definition regions of the solutions are bounded by curves of the form u = f(η), as labelled on the figures (coloured curves). Among those, u+ and u− are given in Eq. (C.5). The junction points (black labels) are given in Table C.1.

In the text

Fig. C.2 Eccentricity e of the satellite at equilibrium P as a function of the parameters. The solutions lie on a three-dimensional surface shown here from two viewing angles. This figure is obtained by stitching together the patches represented in Fig. C.1. For better comparison, the dividing curves are drawn on the surface using the same colour code. The top tube-like portion of the surface has been cut; as shown in Fig. C.1, it extends up to a → ∞ (i.e. u → 0). The equilibrium P is singular along the magenta curve, called in the text.

In the text

	Fig. C.3 Stability of the equilibrium point P as a function of the parameters. The stable and unstable regions are painted in grey and red, respectively, and they are projected on the three-dimensional surface describing the eccentricity of the satellite. As in Fig. C.2, the surface is seen from two viewing angles and the top tube-like portion has been cut for better readability.
In the text

	Fig. C.4 Same as Fig. C.3 but replacing the coordinates a∕rM and e by rM ∕a and e2, respectively.Contrary to previous figures, the top tube-like portion of the surface is seen in its totality.
In the text

	Fig. C.5 Eccentricity e of the satellite at equilibrium P as a function of the parameters. The eccentricity is given by Eq. (C.7); it only depends on the semi-major axis a∕rM of the satellite. The equilibrium is stable wherever it exists.
In the text

	Fig. C.6 Same as Fig. C.1, but showing the location of the eccentric equilibrium P as a function of the parameters. On the top panel, the definition region of cos2 ε+ is very narrow, comprised between the curves u = η3(5 − 4η2) and u = u−(η).
In the text

	Fig. C.7 Eccentricity e of the satellite at equilibrium P as a functionof the parameters. This figure is obtained by stitching together the patches represented in Fig. C.6. The top tube-like portion of the surface has been cut; as shown in Fig. C.6, it extends up to a → ∞ (i.e. u → 0). The equilibrium P is singular along the magenta curve, called in the text.
In the text

	Fig. C.8 Stability of the equilibrium point P as a function of the parameters. The stable and unstable regions are painted in grey and red, respectively, and they are projected on the three-dimensional surface describing the eccentricity of the satellite. As in Fig. C.7, the surface is seen from two viewing angles and the top tube-like portion has been cut for better readability.
In the text

	Fig. C.9 Same as Fig. C.8 but replacing the coordinates a∕rM and e by rM ∕a and e2, respectively.Contrary to previous figures, the top tube-like portion of the surface is seen in its totality.
In the text

	Fig. C.10 Eccentricity e of the satellite at equilibrium P as a function of the parameters. The colour gives its stability: grey for stable and red for unstable. P is unstable in the whole parameter space except along , which results from the merging of P and P into a degenerate equilibrium.
In the text

Fig. D.1 Same as Fig. 17, but using the self-consistent precession rate of Boué & Laskar (2006) given in Eq. (D.1). The physical parameters are those of Titan and Saturn, except that the mass of Titan is adjusted to produce the same mass parameters η as in Fig. 17. In order to ease the comparison, the ridge line and singular point (in red) are those of the simplified model, and the labelled level curves are the same as in Fig. 17 (except the level 0.021 in panel d which is added here).

In the text