Issue 
A&A
Volume 654, October 2021



Article Number  A83  
Number of page(s)  32  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202141467  
Published online  15 October 2021 
Future destabilisation of Titan as a result of Saturn’s tilting
^{1}
IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Université de Lille,
75014
Paris,
France
email: melaine.saillenfest@obspm.fr
^{2}
Department of Mathematics, University of Pisa,
Largo Bruno Pontecorvo 5,
56127
Pisa,
Italy
Received:
4
June
2021
Accepted:
6
July
2021
Context. As a result of Titan’s migration and Saturn’s probable capture in secular spin–orbit resonance, recent works show that Saturn’s obliquity could be steadily increasing today and may reach large values in the next billions of years. Satellites around highobliquity planets are known to be unstable in the vicinity of their Laplace radius, but the approximations used so far for Saturn’s spin axis are invalidated in this regime.
Aims. We aim to investigate the behaviour of a planet and its satellite when the satellite crosses its Laplace radius while the planet is locked in secular spin–orbit resonance.
Methods. We expand on previous works and revisit the concept of Laplace surface. We use it to build an averaged analytical model that couples the planetary spinaxis and satellite dynamics.
Results. We show that the dynamics is organised around a critical point, S_{1}, at which the phasespace structure is singular, located at 90° obliquity and near the Laplace radius. If the spinaxis precession rate of the planet is maintained fixed by a resonance while the satellite migrates outwards or inwards, then S_{1} acts as an attractor towards which the system is forced to evolve. When it reaches the vicinity of S_{1}, the entire system breaks down, either because the planet is expelled from the secular spin–orbit resonance or because the satellite is ejected or collides into the planet.
Conclusions. Provided that Titan’s migration is not halted in the future, Titan and Saturn may reach instability between a few gigayears and several tens of gigayears from now, depending on Titan’s migration rate. The evolution would destabilise Titan and drive Saturn towards an obliquity of 90°. Our findings may have important consequences for Uranus. They also provide a straightforward mechanism for producing transiting exoplanets with a faceon massive ring, a configuration that is often put forward to explain some superpuff exoplanets.
Key words: planets and satellites: dynamical evolution and stability / planets and satellites: formation / celestial mechanics
© M. Saillenfest and G. Lari 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
A secular spin–orbit resonance occurs when the spinaxis 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 spinaxis dynamics of a whole multiharmonic 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, higherorder 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 longterm spinaxis 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_{8}. 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 nearzero 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 nonzero 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 closein 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 nontrivial 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 reinitiates the apsidal precession, with the result of averaging to zero the solar torque and stopping the eccentricity increase^{1}. 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^{2}.
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 closein satellite regime, and a fortiori if the satellite becomes unstable, the model used by Saillenfest et al. (2021a,b) for Saturn’s spinaxis 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 largeobliquity 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 largeobliquity regime. In Sect. 3, we describe the influence of the satellite’s dynamics on the spinaxis 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 longterm 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 semimajor 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 semimajor 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_{2} are the gravitational parameter and the second zonal gravity coefficient of the host planet, and μ_{⊙}, a_{⊙}, and e_{⊙} are the gravitational parameter, the semimajor 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 multipolar 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 faraway 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 spinaxis 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^{3}. 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 semimajor 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 transNeptunian 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 zeroinclination 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 nonsingular 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_{1}, P_{2}, and P_{3}. In the work of Tremaine et al. (2009), the points P_{1} and P_{3} are called ‘circular coplanar Laplace equilibria’, and the point P_{2} 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_{1}, defined by (11)
At point S_{1}, the Laplace states P_{1} and P_{3} 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_{1} by smoothly changing the parameters inverts the locations of P_{1} and P_{3}. We stress that in Fig. 4 the apparent jumps of P_{1} (for a < r_{M}) and P_{3} (for a > r_{M}) are only coordinate singularities in which P_{1} or P_{3} smoothly pass through the pole of the sphere (see Fig. 2). On the contrary, the jump observed at point S_{1} (a = r_{M}) is a real singularity.
Another singularity occurs for a null or 180° obliquity. We call S_{2} the corresponding region of the parameter space, defined by (12)
In region S_{2}, the Laplace states P_{2} and P_{3} 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_{1} and S_{2} of the parameter space can be visualised in Fig. 5.
Apart from regions S_{1} (in which P_{1} becomes singular) and S_{2} (in which P_{2} 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 xaxis 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_{1} around the pole and P_{3} remains near 90°; for a > r_{M}, on the contrary, varying the obliquity makes P_{1} and P_{3} 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_{1} is stable to inclination variations. As illustrated by Fig. 3, it corresponds to an orbit lying on the equator plane for closein satellites (a ≪ r_{M}), and on the ecliptic plane for faraway satellites (a ≫ r_{M}). In between, P_{1} corresponds to an intermediate tilt between the equator and the ecliptic. As a result of eccentricity and inclination damping, P_{1} is expected to be the birth place of regular satellites formed in a circumplanetary disc. P_{1} 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_{1} is given by one of the two solutions of the equation^{4} (13)
the second solution being the inclination of P_{3}. 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_{1} 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_{1}. Through a smooth variation of parameters, the satellite can therefore reach the singular point S_{1} from any orbital inclination between 0° and 180°. This property has important consequences for the spinaxis dynamics of the host planet, as we discuss in Sect. 3. We also show that r_{F} divides the closein and faraway satellite regimes considered in previous works.
Tremaine et al. (2009) give a compact expression for the frequency of smallamplitude oscillations around P_{1}, which can be written as (16)
where I_{Q1} is the equatorial inclination at P_{1} 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 zeroobliquity planet, Eq. (16) simplifies to (17)
It shows that the timescale parameter κ defined in Eq. (9) is the oscillation frequency around P_{1} for ε = 0 and at a radius a = r_{M}. See Appendix B for the limit value of ξ_{1} in the regions of parameter space where Eq. (16) looks undefined. As a summary, Fig. 7 shows the libration period around P_{1} in the whole parameter space. The stability properties of P_{2} and P_{3} 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_{2} and P_{3} is provided in Appendix B.
From P_{3} emerges the separatrix that divides the regions of oscillations around P_{1} and around P_{2}. 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_{2} is zero in parameter region S_{2} (as expected from Fig. 2), and that it increases for growing a and for decreasing cos^{2}ε. This has important consequences for the emergence of chaos discussed in Sect. 4. At the singular point S_{1}, the island covers the whole sphere.
Apart from the singular regions S_{1} and S_{2}, the continuous behaviour of the Laplace states all over the parameter space is crucial for the longterm 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_{1} (for oscillations around P_{1}) or S_{2} (for oscillations around P_{2}), 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.
Parameters and dynamical timescales of some satellites and their host planets in the Solar System.
Fig. 1 Level curves of the Hamiltonian function for a circularorbit. The parameters are a∕r_{M} = 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 P_{1} and P_{3} lying at δ_{Q} = 0 and for P_{2} 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 zaxis is along the spin axis of the planet. The xaxis 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) = (sinI_{Q} sinδ_{Q}, −sinI_{Q} cosδ_{Q}, cos I_{Q}). The colour code is the same as in Fig. 1. Panel a: same parameters as Fig. 1. Panel b: parameter region S_{1}, defined by a∕r_{M} = 1 and ε = 90°. The magenta curve is made of an infinity of stable equilibria resulting from the merging of P_{1} with the separatrix emerging from P_{3}. Panel c: parameter region S_{2}, defined by a∕r_{M} > 0 and ε = 0° (or 180°). The brown curve is made of an infinity of stable equilibria resulting from the merging of P_{2} with the separatrix emerging from P_{3}. 
Fig. 3 Location of the Laplace states as a function of a∕r_{M}. Each panel features a fixed value of ε (see labels). The level I_{Q} = ε 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 P_{1} and P_{3}. 
Fig. 4 Location of the Laplace states as a function of ε. Each panel features a fixed value of a∕r_{M} (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 P_{1} and P_{2} are stable to inclination variations and P_{3} is unstable. At point S_{1} of the parameter space, P_{1} and P_{3} degenerate into the equilibrium circle that is stableto inclination variations. In region S_{2} of the parameter space, P_{2} and P_{3} degenerate into the equilibrium circle that is stableto inclination variations. In regions E_{1}, E_{2}, and E_{3}, respectively, P_{1}, P_{2}, and P_{3} are unstable to eccentricity variations. In region S_{2} ∪E_{2}, is unstable to eccentricity variations. At point S_{1}, is stable to eccentricity variations provided that the equatorial inclination of the satellite verifies cos^{2} I_{Q} ≤ 1∕5. 
Fig. 6 Inclination of the Laplace state P_{1} as a function of the parameters (see Eq. (14)). Some level curves are labelled in black. The critical radii and r_{M} and r_{F} are represented in white. At a = r_{M}, P_{1} lies exactly halfway between the equator and the ecliptic (i.e. I_{Q} = ε∕2 for a prograde spin). The abrupt transition of I_{Q} from 0° to 180° at ε = 90° is a coordinate singularity. The S_{1} point is a real singularity (see text). 
Fig. 7 Period of small inclination oscillations around the Laplace state P_{1} 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 S_{1} 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 nonzero offset of eccentricity remains small or grows big over time. As before, we focus on the Laplace state P_{1}; the description of P_{2}, P_{3}, and the degenerate circles and are provided in Appendix B.
The eigenvalues of the eccentricity linear subsystem inform us about their stability against eccentricity growth. Tremaine et al. (2009) give a compact expression for the frequency of smallamplitude eccentricity oscillations around P_{1}. It can be written (19)
where I_{Q1} is the equatorial inclination at P_{1} given by Eq. (14). See Appendix B for the limit value of μ_{1} 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_{1} is stable in all the parameter space except in a small closed region resembling a cardioid. We call E_{1} this region of the parameter space; it can be visualised in Figs. 5 and 8.
The boundary of E_{1} is given by the roots of , which have a closedform analytical expression. We first define two critical radii r_{1} and r_{2} as (20)
As shown in Fig. 8, the radii r_{1} and r_{2} mark the leftmost and rightmost limits of E_{1}, and the boundary of E_{1} has a cusp at the singular point S_{1}. Noting , the boundary of E_{1} can be expressed piecewise as (21)
for r_{1} ≤ a ≤ r_{2}, and (22)
for r_{M} ≤ a ≤ r_{2}. Equation (22) corresponds to the cusp portion of the curve, and the two portions meet at a = r_{2} (see Fig. 8). The obliquity ε ≈ 68.875° quoted by Tremaine et al. (2009) as the minimum value where P_{1} can be unstable is reached at . It has actually the following closedform: (23)
Interestingly, does not go to zero at the singular point S_{1}, 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_{1} 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_{1} while migrating inwards rather than outwards (see Fig. 8).
The stability properties of P_{2} and P_{3} to eccentricity variations are given in Appendix B. We show that they are unstable in the regions E_{2} and E_{3}, respectively, illustrated in Fig. 5. We note that E_{2} entirely contains the region E_{1}; hence, in region E_{1} all Laplace states are unstable to at least eccentricity or inclination variations.
Fig. 8 Region E_{1} of the parameter space where the Laplace state P_{1} 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 E_{1} is highlighted with a black contour, obtained using the closedform expression in Eqs. (21) and (22). Region E_{1} is entirely contained inside the dashed brown lines, whose left and right limits r_{1} and r_{2} are given in Eq. (20), and whose bottom and top limits are given by Eq. (23). The Laplace state P_{1} is singular atpoint S_{1}, in which is discontinuous. 
2.4 Eccentric Laplace states
Along the boundaries of the regions E_{1} and E_{2}, where the Laplace states P_{1} and P_{2} 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_{3} region (the Vshaped boundary for a < r_{L} and the droplike boundary for a > r_{L}; see Fig. 5 and Appendix B), the Laplace state P_{3} 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_{1}, E_{2}, or E_{3} 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_{1} 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)
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 nonzero 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 threedimensional 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 cardioidlike boundary of the E_{1} region. Since Fig. 9 is the projection of a complex threedimensional 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 threedimensional surface turns round to higher semimajor axes again.
Along the threedimensional 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_{1}, at which P_{1} and P_{3} are degenerate. In Fig. 9, the curve is visible between the points labelled S_{1} 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 threedimensional 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_{1}, 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_{1}. Along the boundary of the E_{1} region, the two eigenfrequencies tend to the oscillation frequencies ξ_{1} and μ_{1} around P_{1} given at Eqs. (16) and (19), confirming that the eccentric equilibrium P bifurcates from the circular equilibrium P_{1}.
As stressed by Tremaine et al. (2009), the eccentric equilibrium P is stable near its bifurcation from P_{1} 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 tubelike 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_{1} region, even where the eccentric equilibrium P should theoretically be stable. In the case of transNeptunian 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 quasiparabolic 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 threedimensional 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 threedimensional surface. The colour is the same as in Fig. 9, top panel. In Fig. 9, the top tubelike 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 S_{1} and S are labelled, and the cutting line of the threedimensional surface is shown by a grey line (same as Fig. 9). Threedimensional figures are provided in Appendix C, where the geometry of the stable regions is clearer. 
3 Spinaxis 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 spinaxis 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 spinaxis precession frequency to a fixed value.
A selfconsistent 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_{1}, 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 spinaxis 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 spinaxis 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_{2} 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 longterm stable, so that its secular orbital motion can be expressed (at least locally) in convergent quasiperiodic 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 quasiperiodic 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 longterm stable planetary system, then the Hamiltonian in Eq. (27) can be considered as a perturbation to the unperturbed Hamiltonian . In this setting, the longterm spinaxis dynamics of the planet is shaped by resonances between the unperturbed spinaxis precession frequency and the forcing frequencies appearing in Eq. (33). In the Solar System, the orbital precession motions of the terrestrial planets contain numerous largeamplitude 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 spinaxis 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). Secondorder resonances involve two terms in the series, and thirdorder resonances involve three. Eccentricitydriven 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 spinaxis precession frequency is forced to remain approximatively constant, equal to a combination of frequencies ν_{k}.
In the vicinity of a firstorder secular spin–orbit resonance, the Hamiltonian function reduces to the wellknown ‘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 spinaxis precession frequency of the planet. Contrary to Saillenfest et al. (2019b), we do not expand the eccentricity variations of the planet in quasiperiodic 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_{2} and C_{3} are present (see Fig. 13a). If γ^{2∕3} + β^{2∕3} < 1, then all four Cassini states are present, and C_{2} 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 spinaxis 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 ε_{2} = 0° and ε_{3} = 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_{2} 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_{1} and C_{4} 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_{2} becomes vanishingly narrow and it merges with the Cassini states C_{2} and C_{4}, 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 nonzero libration amplitudes).
Hence, as a summary, if we increase the parameter α from 0 to ∞, the Cassini state C_{2} gradually passes from ε_{2} ≈ 0° (without separatrix) to ε_{2} = 90° (inside the resonance separatrix). These properties will be important below.
Fig. 12 Spinaxis dynamics in the vicinity of a firstorder 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 spinaxis precession motion. In Sect. 3.1, this property allowed us to treat the orbital variations as a forcing term in the spinaxis 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_{1} 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 spinaxis precession timescale of the planet as T = 2π∕p, where has been introduced in Sect. 3.1. Without satellite, the free spinaxis 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_{1}, this timescale condition may be violated at the border of region E_{1} or in the extreme vicinity of the singular point S_{1}, 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 spinaxis precession frequency of their host planet. If the satellites oscillate about the Laplace state P_{1} with τ ≪ T, French et al. (1993) give an elegant expression for their contribution: one must simply replace J_{2} and λ in Eq. (29) by the effective values (39)
where m_{k}, a_{k} and n_{k} are the mass, the semimajor 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^{5}. Using these expressions, the free spinaxis precession frequency Ω_{0} 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 closein regime considered previously by Saillenfest et al. (2020, 2021b). It only requires that the satellites oscillate around the Laplace state P_{1}, which is what we expect for any regular satellite (see e.g. Tremaine et al. 2009), unless it reaches the highobliquity unstable region E_{1} (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 selfconsistent, because the satellite is considered to be massless when dealing with the orbital dynamics (Sect. 2) but massive when computing its longterm 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^{−4} of Saturn’s, the inclination L_{k} and the precession rate Ω_{0} obtained through Eqs. (14) and (39) are very close to those obtained using the selfconsistent (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 satellitetoplanet 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 zeroobliquity planet simplifies to (40)
Likewise, for ε = 90° the parameter becomes (41)
However, since Ω_{0} is factored by cos ε, the spinaxis precession frequency of the planet is in any case zero for ε = 90°, except at a = r_{M}, where the Laplace state P_{1} 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 fastspinning planet like Saturn, the oblateness coefficient J_{2} 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_{2} 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 lowmass satellite can have a large mass parameter η. For instance, Titan has a mass of m∕M ≈ 10^{−4} but a mass parameter η ≈ 10, and Titan greatly affects Saturn’s spinaxis 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 spinaxis 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 spinaxis 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 closein and faraway 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 Ω_{0} divides the closein and faraway satellite regimes, characterised by the wellknown power laws in a^{2} and a^{−3}, respectively. The magnitude of Ω_{0} 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 Ω_{0} 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 spinaxis precession frequency of the planet is intimately linked to the properties of the Laplace state P_{1} of its satellite. We further analyse this relation below.
Figure 16 shows the value of Ω_{0} at its maximum (i.e. at a = r_{F}) as a function of the planet’s obliquity. The global maximum of Ω_{0} is reached at ε = 0°, for which . It has value (44)
For an obliquity ε → 90°, the maximum of Ω_{0} is reached at a = r_{M}, that is, at the singular point S_{1} of the satellite’s dynamics (see Sect. 2). Figure 16 shows that the maximum value of Ω_{0} has two different limits at point S_{1} according to whether the obliquity tends to 90° from above or from below. These limits are Ω_{0} = ± Ω_{S}, where (45)
We note that these limits are nonzero 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 nonzero finite value when cosε tends to zero (α′ actually goes to infinity). We deduce that when the system approaches the singular point S_{1}, 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 Ω_{0} 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 spinaxis 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π∕Ω_{0}, instead of the frequency Ω_{0} 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 selfconsistent models.
For a very small mass parameter η ≪ 1 (e.g. for Deimos), the effect of the satellite is negligible and the level curves of Ω_{0} 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 Ω_{0} 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 Ω_{0} 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 Ω_{0}: instead, the level curves must go around the ridge line, and they converge at the singular point S_{1} (see Fig. 17). More precisely, all level curves with frequency values Ω_{0} ≤ Ω_{S} converge to S_{1}. For a large mass parameter η, this condition is verified for almost all level curves of Ω_{0}. This property is related to the orbital plane of the satellite, which can reach S_{1} from any inclination (see Fig. 6).
Figure 17 shows that for a large enough value of η, numerous level curves of Ω_{0} connect the singular region S_{2} (i.e. ε = 0) to the singular point S_{1}. 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_{1} 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 Ω_{0} going from ε = 0° to ε = 90° verify p ≤ Ω_{0} ≤ Ω_{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 Ω_{0} = 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_{1} is surrounded by the unstable region E_{1} (see Sect. 2.3), we expect the satellite’s orbit to become unstable before actually reaching S_{1}. 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_{1}, 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 xyplane. The obliquity ε is the tilt from the zaxis, and the resonance angle σ is the polar angle measured in the xyplane. 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 Spinaxis 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 spinaxis precession frequency of the planet as a function of its obliquity. The curve is obtained by injecting a = r_{F} 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 Spinaxis 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 spinaxis 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 S_{1} and the ridge line a = r_{F} are indicated in red. 
Largest terms of Saturn’s inclination series in the J2000 ecliptic and equinox reference frame.
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 longterm 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 quasiperiodic 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 spinaxis 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^{−1}. 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 firstorder secular spin–orbit resonances, even when considering the full quasiperiodic solution of Laskar (1990). The largest resonance is with ν_{3} = s_{8} (see Table 2). Other resonances can be identified in Table A.2 of Saillenfest et al. (2021b): there are two prograde resonances (ν_{19} and ν_{51}) and two retrograde resonances (ν_{28} and ν_{44}). As revealed by their high index in the orbital series, these four resonances are very small. This explains why no relevant second or higherorder resonance can possibly affect Saturn. Then, we know that all precession frequencies Ω_{0} verifying p ≤ Ω_{0} ≤ Ω_{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^{−1} and Ω_{S} ≈ 1.20″ yr^{−1}, so this condition is met by all five resonances mentioned above (even though ν_{51} is right at the limit, since ν_{51}≈ Ω_{S}).
Figure 18 shows the location and width of the firstorder secular spin–orbit resonances reachable by Saturn as a function of the distance of Titan. The full effect of Titan on Saturn’s spinaxis 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_{1}. 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_{1} before reaching the singularity S_{1}; 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_{1} 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^{−1}, Saturn would be unable to reach ν_{14} and ν_{15} for any distance of Titan, even though the trajectory in Fig. 17 of Saillenfest et al. (2021b) appears at the same height as ν_{14} and ν_{15} on the graph. Moreover, as shown in Sect. 3.2, α can tend toinfinity while the spinaxis precession frequency remains finite, which is quite counterintuitive. 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 spinaxis precession.
Assuming that Saturn follows the centre of the resonance with s_{8} (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 closesatellite 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_{1}, where the libration period of σ_{3} 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 spinaxis precession period is maintained to a constant value as it enters the resonance.
Titan is located today at a mean semimajor axis of a_{0} ≈ 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_{0} 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 radioscience 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 radioscience experiments. Figures 18 and 19 show that the system will first reach the unstable region E_{1} 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 radioscience 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_{1}. For this reason, Fig. 19 can only provide a qualitative picture of the evolution of the system. Because of the intricate and multitimescale nature of the dynamics, building a selfconsistent numerical model for the evolution of Saturn and Titan is challenging: the system involves the orbital dynamics of Titan and Saturn, the spinaxis 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 multiharmonic 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 spinaxis 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 firstorder 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 semimajor 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 = r_{F} dividing the closein and faraway satellite regimes is represented by a dashed blue line. The region E_{1} where the satellite’s Laplace state P_{1} is unstable isshown in black. The singular point S_{1} 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 s_{8}. 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 selfexplanatory. 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 semimajor axis of Titan and the obliquity of Saturn evolve over time as slowvarying parameters. As a first approximation, their evolution law is provided by the topleft panel of Fig. 19 (blue curve). Therefore, in our simulations we make the obliquity of Saturn and the semimajor 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_{1}. 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_{1} (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_{1} (see the bottom left panel of Fig. 19).
For comparison, we also performed direct unaveraged numerical integrations of the restricted threebody problem including Saturn’s oblateness. In order to reproduce the drift in Titan’s semimajor 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 spinaxis 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 spinaxis 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 spinaxis motion should instead be integrated as well in a selfconsistent 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_{1} and the vicinity of the singular point S_{1}. 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 Longterm dynamics of Titan as it migrates away from Saturn. Saturn is assumed to remain locked in secular spin–orbit resonance with s_{8} 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 semimajor 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 P_{1}, 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 threebody problem including Saturn’s oblateness. The output is then digitally filtered to remove its shortperiod 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 spinaxis 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 firstorder secular spin–orbit resonances, and the blue hatched area is the region E_{1} where the Titan’s circular equilibrium is unstable (same as Fig. 18). 
4.3 Saturn’s spinaxis
Using the Hamiltonian function in Eq. (27), our setting is similar to the migration simulations of Saillenfest et al. (2021b), except that both the semimajor 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 closein or faraway 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_{1}, 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 longterm 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_{1}, where the inclination variations of Titan as a function of obliquity are extremely sharp (and discontinuous exactly at S_{1}). But as shown in Sect. 4.2, Titan is expected anyway to be destabilised before actually reaching S_{1}. Therefore, we stress that this model is not selfconsistent; 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 illdefined near the singular point S_{1} (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 ν_{51} 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_{1} (hatched region). Eventually, Saturn crosses the s_{8} resonance again when Titan passes in the farsatellite 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_{8} (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_{1} 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_{1}, 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 selfconsistent numerical model; these questions are left for future works.
In panels e and f, Saturn did not reach yet the resonance with s_{8} 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_{0} 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 darkblue 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 largeamplitude 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 = r_{M}. 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_{8}. 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_{1}, P_{2}, and P_{3}) are organised around two critical regions S_{1} and S_{2} in the parameter space, in which thepairs (P_{1}, P_{3}) or (P_{2}, P_{3}), respectively, are degenerate. In particular, S_{1} is defined as the point where the planet’s obliquity is ε = 90° and the satellite’s semimajor axis is a = 2^{1∕5}r_{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_{1} and S_{2} 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_{1}). 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 spinaxis precession rate Ω_{0} 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 Ω_{0} 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 Ω_{0} in the plane (a, ε). When the satellite’s inclination is fully taken into account (i.e. without the usual closein or faraway approximations), the singularity S_{1} in its orbital dynamics is transferred to the spinaxis precession rate of the planet. As a result, if the satellite is massive enough, most level curves of Ω_{0} converge towards the singularity S_{1}. 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_{8} resonance, and confirmed using numerical simulations.
In the vicinity of S_{1}, however, several kinds of instability are expected to happen. Firstly, S_{1} lies at the border of a region in the parameter space where the circular Laplace state P_{1} 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_{1}, where the instability is strongest, before being destabilised.
Secondly, all neighbouring secular spin–orbit resonances converge towards S_{1}. 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_{1}.
Thirdly, the widths of all secular spin–orbit resonances decrease in the vicinity of S_{1}. 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 selfconsistent model coupling the orbital motions of Saturn and Titan, the spinaxis dynamics of Saturn torqued by the Sun and by Titan, and the multiharmonic 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 resonancelocking 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 ‘superpuff’ exoplanets thought to possess a massive faceon 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)
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_{1}, P_{2}, and P_{3}. Even though the equilibrium P_{1} 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_{1}
P_{1} is stable to inclination variations (see Fig. 1). For δ_{Q} = 0, the equatorial inclination I_{Q1} of the satellite at P_{1} is given by Eq. (14). Inclination oscillations around the Laplace state P_{1} have frequency ξ_{1} 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_{1} (i.e. a = r_{M}), since we retrieve the eigenfrequency along the degenerate stable circle . For a very close satellite, ξ_{1} 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_{1} is stable to eccentricity variations, except in the region E_{1} described in Sect. 2.3. Eccentricity oscillations around P_{1} have frequency μ_{1} given in Eq. (19). For a zeroobliquity 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 nearzero obliquity, the oscillations of the inclination and eccentricity of a satellite around the classical Laplace state P_{1} 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_{1} produces a discontinuity of . This could have been expected since P_{1} and P_{3} 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, μ_{1} 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^{2}I_{Q} ≤ 1∕5, that is, if 63° ≲ I_{Q} ≲ 117°. This corresponds to the region where the J_{2}induced precession of ω_{Q} is negative.
Appendix B.2 Circular equilibrium P_{2}
P_{2} 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_{2} corresponds to an orbit that is perpendicular to both the equator and the ecliptic planes. The frequency of smallamplitude inclination oscillations around P_{2} is given by (B.7)
As expected, the value of is always negative, and it tends to zero in the region S_{2} of the parameter space (i.e. ε = 0° or 180°), since we retrieve the eigenfrequency along the degenerate stable circle .
The frequency of smallamplitude eccentricity oscillations around P_{2} is given by (B.8)
We deduce that P_{2} is stable to eccentricity variations only for close enough satellites that verify a ≤ r_{3}, where we define (B.9)
The region where P_{2} 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_{2}, it is only stable for closeenough satellites.
Appendix B.3 Circular equilibrium P_{3}
P_{3} is unstable to inclination variations (see Fig. 1). The squared eigenvalue of the linearised problem at P_{3} 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_{3} exists, and tends to 0 at the singular regions S_{1} and S_{2}, where we retrieve the eigenfrequencies and along the degenerate stable circles and .
Even though P_{3} 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_{3} is stable to eccentricity variations for closeenough satellites and in a region resembling a water drop. The squared eigenvalue of the linearised problem at P_{3} is given by Eq. (19) where I_{Q1} must be replaced by I_{Q3} = I_{Q1} ± 90°. We call E_{3} the unstable region where . Its border can be expressed piecewise as (B.11)
for r_{3} ≤ 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_{3} is given in Eq. (B.9).
Fig. B.1 Region E_{3} of the parameter space where the Laplace state P_{3} 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 E_{3} is obtained using the closedform expression in Eqs. (B.11) and (22). The droplike 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_{1}, P_{2}, and P_{3} are described. Wherever they are defined, all three of them are stable to eccentricity variations, except in the regions E_{1}, E_{2}, and E_{3}, 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_{1} and P_{2} where these equilibria become unstable (i.e. at the borders of the E_{1} and E_{2} 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_{3} at the border of region E_{3}; 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_{1} region defined at Eqs. (21) and (22); this shows that P indeed bifurcates away from P_{1} 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 rootfinding algorithms to reverse the expression for η. We define the obliquity solutions ε_{+} and ε_{−} as (C.2)
The values of cos^{2}ε_{+} and cos^{2}ε_{−} 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^{2}ε_{+} and cos^{2}ε_{−} 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 closedform 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 threedimensional 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_{1} where the circular equilibrium P_{1} 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 tubelike portion of the threedimensional surface has been cut off.
Along the threedimensional curve defined in Eq. (26) and depicted by a magenta curve in Figs. C.1 and C.2, the threedimensional surface has a cusp. As discussed in the main text, the curve is the eccentric continuation of the singular point S_{1} described in Sect. 2.2, where the equilibrium point P_{1} is singular. We note that pierces the threedimensional 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 closedform 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 threedimensional 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^{2}), where the full threedimensional 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 semimajor 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^{6}.
Fig. C.1 Location of the eccentric equilibrium P as a function of the parameters. The colour shades show the two closedform 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 threedimensional 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 tubelike 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_{3} defined in Eq. (B.9); this shows that P indeed bifurcates away from P_{2} where the latter becomes unstable. Hence, P exists for a ≥ r_{3} and the eccentricity of the satellite at P is (C.7)
The bifurcation of P from P_{2} can be visualised as a threedimensional 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 threedimensional surface describing the eccentricity of the satellite. As in Fig. C.2, the surface is seen from two viewing angles and the top tubelike portion has been cut for better readability. 
Fig. C.4 Same as Fig. C.3 but replacing the coordinates a∕r_{M} and e by r_{M} ∕a and e^{2}, respectively.Contrary to previous figures, the top tubelike 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)
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_{3}, we note that and are equal to the oscillation frequencies ξ_{2} and μ_{2} around P_{2} given in Eqs. (B.7) and (B.8), confirming that the eccentric equilibrium P bifurcates away from the circular equilibrium P_{2}.
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 semimajor axis a∕r_{M} 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_{3} is unstable to eccentricity variations in the region E_{3} of the parameter space. In its illustrations in Figs. 5 and B.1, we see that the border of E_{3} can be divided into two components: a small droplike boundary for a > r_{L} (whose expression is given by Eq. (22)), and a Vshaped 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_{3} bifurcates into two distinct eccentric equilibria. We call them P and P, respectively, for the eccentric equilibria bifurcating away from the droplike boundary and from the Vshaped 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 droplike boundary of the E_{3} region; this shows that P indeed bifurcates away from P_{3} 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 closedform 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^{2}ε_{+} and cos^{2}ε_{−} 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 threedimensional surface illustrated in Fig. C.7. As expected, the bottom section of this surface (i.e. at e = 0) coincides with the droplike boundary of the region E_{3} where the circular equilibrium P_{3} 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 threedimensional 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 threedimensional 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 threedimensional surface. Interestingly, a small stable region exists on the top portion of the threedimensional 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_{3} 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 loweccentricity portion of the threedimensional 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^{2}), where the full threedimensional 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)
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 Vshaped boundary of the region E_{3} (see Eq. (B.11)), but where u is replaced by w. This shows that P indeed bifurcates away from P_{3} 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 threedimensional domain (C.14)
is the eccentric continuation of the singular region S_{2} described in Sect. 2, where P_{2} and P_{3} 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_{3} and the two curves that define can be visualised in Fig. C.10. We see that is the contact region between the two threedimensional 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 cos^{2} ε_{+} 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 tubelike 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 threedimensional surface describing the eccentricity of the satellite. As in Fig. C.7, the surface is seen from two viewing angles and the top tubelike portion has been cut for better readability. 
Fig. C.9 Same as Fig. C.8 but replacing the coordinates a∕r_{M} and e by r_{M} ∕a and e^{2}, respectively.Contrary to previous figures, the top tubelike 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 Selfconsistent formula for the spinaxis precession rate
In Eq. (42), we give a simplified expression for the spinaxis 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 closedform 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 selfconsistent. The satellite is considered massless to compute its equilibrium inclination (Laplace state P_{1}), 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 longterm 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 selfconsistent 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 spinaxis 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, highorder 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 timeevolution 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)
are constant, and the physical parameters of the system are contained in the coefficients (D.4)
In this selfconsistent 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 masslesscase 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 selfconsistent model of Boué & Laskar (2006), the precession rate in Eq. (D.1) cannot be expressed in terms of a few macroparameters, 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 nonnegligible (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 satellitetoplanet mass ratios. For the EarthMoon system, whose mass ratio is as large as about 0.01, a selfconsistent precession model would be required to get precise quantitative results.
Fig. D.1 Same as Fig. 17, but using the selfconsistent 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
 Akinsanmi, B., Santos, N. C., Faria, J. P., et al. 2020, A&A, 635, A8 [Google Scholar]
 Atobe, K., Ida, S., & Ito, T. 2004, Icarus, 168, 223 [Google Scholar]
 Boué, G., & Laskar, J. 2006, Icarus, 185, 312 [NASA ADS] [CrossRef] [Google Scholar]
 Boué, G., & Laskar, J. 2010, ApJ, 712, L44 [Google Scholar]
 Boué, G., Laskar, J., & Kuchynka, P. 2009, ApJ, 702, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Brasser, R., & Lee, M. H. 2015, AJ, 150, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Colombo, G. 1966, AJ, 71, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Ćuk, M., Hamilton, D. P., Lock, S. J., & Stewart, S. T. 2016, Nature, 539, 402 [Google Scholar]
 Deitrick, R., Barnes, R., Quinn, T. R., et al. 2018, AJ, 155, 60 [Google Scholar]
 Fouchard, M. 2004, MNRAS, 349, 347 [Google Scholar]
 French, R. G., Nicholson, P. D., Cooke, M. L., et al. 1993, Icarus, 103, 163 [CrossRef] [Google Scholar]
 Frouard, J., Fouchard, M., & Vienne, A. 2010, A&A, 515, A54 [Google Scholar]
 Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P. 1966, Rev. Geophys. Space Phys., 4, 411 [Google Scholar]
 Hamilton, D. P., & Ward, W. R. 2004, AJ, 128, 2510 [NASA ADS] [CrossRef] [Google Scholar]
 Haponiak, J., Breiter, S., & Vokrouhlický, D. 2020, Celest. Mech. Dyn. Astron., 132, 24 [CrossRef] [Google Scholar]
 Henrard, J., & Murigande, C. 1987, Celest. Mech., 40, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Ito, T., & Ohtsuka, K. 2019, Monographs on Environment, Earth Planets, 7, 1 [Google Scholar]
 Kreyche, S. M., Barnes, J. W., Quarles, B. L., et al. 2020, Planet. Sci. J., 1, 8 [Google Scholar]
 Lainey, V., Gomez Casajus, L., Fuller, J., et al. 2020, Nat. Astron., 4, 1053 [CrossRef] [Google Scholar]
 Laplace, P. S. 1805, Traité de Mécanique Céleste, 4 (Courcier, Paris) [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [Google Scholar]
 Laskar, J., & Robutel, P. 1993, Nature, 361, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Correia, A. C. M., Gastineau, M., et al. 2004, Icarus, 170, 343 [Google Scholar]
 Li, G., & Batygin, K. 2014, ApJ, 790, 69 [Google Scholar]
 Li, J., & Lai, D. 2020, ApJ, 898, L20 [Google Scholar]
 Martin, R. G., & Armitage, P. J. 2021, ApJ, 912, L16 [Google Scholar]
 Millholland, S., & Batygin, K. 2019, ApJ, 876, 119 [CrossRef] [Google Scholar]
 Millholland, S., & Laughlin, G. 2019, Nat. Astron., 3, 424 [CrossRef] [Google Scholar]
 Murray, C. A. 1989, A&A, 218, 325 [NASA ADS] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge University Press) [Google Scholar]
 Néron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [Google Scholar]
 Peale, S. J. 1969, AJ, 74, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Piro, A. L., & Vissapragada, S. 2020, AJ, 159, 131 [Google Scholar]
 Rogoszinski, Z., & Hamilton, D. P. 2020, ApJ, 888, 60 [CrossRef] [Google Scholar]
 Rogoszinski, Z., & Hamilton, D. P. 2021, Planet. Sci. J., 2, 78 [Google Scholar]
 Saillenfest, M., Fouchard, M., Ito, T., & Higuchi, A. 2019a, A&A, 629, A95 [Google Scholar]
 Saillenfest, M., Laskar, J., & Boué, G. 2019b, A&A, 623, A4 [Google Scholar]
 Saillenfest, M., Lari, G., & Courtot, A. 2020, A&A, 640, A11 [EDP Sciences] [Google Scholar]
 Saillenfest, M., Lari, G., & Boué, G. 2021a, Nat. Astron., 5, 345 [Google Scholar]
 Saillenfest, M., Lari, G., Boué, G., & Courtot, A. 2021b, A&A, 647, A92 [Google Scholar]
 Shan, Y., & Li, G. 2018, AJ, 155, 237 [Google Scholar]
 Souami, D., & Souchay, J. 2012, A&A, 543, A133 [Google Scholar]
 Speedie, J., & Zanazzi, J. J. 2020, MNRAS, 497, 1870 [Google Scholar]
 Tamayo, D., Burns, J. A., Hamilton, D. P., & Nicholson, P. D. 2013, AJ, 145, 54 [Google Scholar]
 Touma, J., & Wisdom, J. 1998, AJ, 115, 1653 [Google Scholar]
 Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706 [NASA ADS] [CrossRef] [Google Scholar]
 Vokrouhlický, D., & Nesvorný, D. 2015, ApJ, 806, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Ward, W. R. 1973, Science, 181, 260 [Google Scholar]
 Ward, W. R. 1975, AJ, 80, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Ward, W. R. 1982, Icarus, 50, 444 [Google Scholar]
 Ward, W. R., & Canup, R. M. 2006, ApJ, 640, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Ward, W. R., & Hamilton, D. P. 2004, AJ, 128, 2501 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, W., & Fabrycky, D. 2019, AAS J. submitted [arXiv:1904.02290] [Google Scholar]
The same mechanism was described by Saillenfest et al. (2019a) as a protection mechanism for inner Oort cloud objects against the action of galactic tides.
The ‘ivection’ resonance mentioned by Speedie & Zanazzi (2020) is order zero in eccentricity; it should not be confused with the mixedtype ‘eviction’ resonance described by Touma & Wisdom (1998), which is order two in eccentricity and can only be triggered once the eccentricity is high enough. See the preprint of Xu & Fabrycky (2019) for more details.
Ćuk et al. (2016) and Speedie & Zanazzi (2020) go one step further and redefine r_{L} by adding the factor 2 in Eq. (6). We rather prefer to introduce a different symbol, because using differing definitions for the classic ‘Laplace radius’ is misleading when it comes to comparison with previous works.
There seems to be a typographical error in Eqs. (22) and (23) of Tremaine et al. (2009): for both equations, the first equality is correct but not the second one. We give the correct expression in Eq. (13).
Ward & Hamilton (2004) give another expression for and λ′ in their Eqs. (2) and (3). When the results are compared to the selfconsistent theory of Boué & Laskar (2006), the expression of Ward & Hamilton (2004) appears to be erroneous. We suspect that Eq. (2) of Ward & Hamilton (2004) actually contains a typographical error. This likely error seems to have been propagated in Eq. (44) of Millholland & Laughlin (2019).
Due to their similar properties, the equilibria P and P have been studied together by Tremaine et al. (2009). Their Fig. 5 features therefore both P (for inclinations smaller than 90°) and P (for inclinations larger than 90°).
All Tables
Parameters and dynamical timescales of some satellites and their host planets in the Solar System.
Largest terms of Saturn’s inclination series in the J2000 ecliptic and equinox reference frame.
All Figures
Fig. 1 Level curves of the Hamiltonian function for a circularorbit. The parameters are a∕r_{M} = 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 P_{1} and P_{3} lying at δ_{Q} = 0 and for P_{2} 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 zaxis is along the spin axis of the planet. The xaxis 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) = (sinI_{Q} sinδ_{Q}, −sinI_{Q} cosδ_{Q}, cos I_{Q}). The colour code is the same as in Fig. 1. Panel a: same parameters as Fig. 1. Panel b: parameter region S_{1}, defined by a∕r_{M} = 1 and ε = 90°. The magenta curve is made of an infinity of stable equilibria resulting from the merging of P_{1} with the separatrix emerging from P_{3}. Panel c: parameter region S_{2}, defined by a∕r_{M} > 0 and ε = 0° (or 180°). The brown curve is made of an infinity of stable equilibria resulting from the merging of P_{2} with the separatrix emerging from P_{3}. 

In the text 
Fig. 3 Location of the Laplace states as a function of a∕r_{M}. Each panel features a fixed value of ε (see labels). The level I_{Q} = ε 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 P_{1} and P_{3}. 

In the text 
Fig. 4 Location of the Laplace states as a function of ε. Each panel features a fixed value of a∕r_{M} (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 P_{1} and P_{2} are stable to inclination variations and P_{3} is unstable. At point S_{1} of the parameter space, P_{1} and P_{3} degenerate into the equilibrium circle that is stableto inclination variations. In region S_{2} of the parameter space, P_{2} and P_{3} degenerate into the equilibrium circle that is stableto inclination variations. In regions E_{1}, E_{2}, and E_{3}, respectively, P_{1}, P_{2}, and P_{3} are unstable to eccentricity variations. In region S_{2} ∪E_{2}, is unstable to eccentricity variations. At point S_{1}, is stable to eccentricity variations provided that the equatorial inclination of the satellite verifies cos^{2} I_{Q} ≤ 1∕5. 

In the text 
Fig. 6 Inclination of the Laplace state P_{1} as a function of the parameters (see Eq. (14)). Some level curves are labelled in black. The critical radii and r_{M} and r_{F} are represented in white. At a = r_{M}, P_{1} lies exactly halfway between the equator and the ecliptic (i.e. I_{Q} = ε∕2 for a prograde spin). The abrupt transition of I_{Q} from 0° to 180° at ε = 90° is a coordinate singularity. The S_{1} point is a real singularity (see text). 

In the text 
Fig. 7 Period of small inclination oscillations around the Laplace state P_{1} 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 S_{1} is indicated in red. 

In the text 
Fig. 8 Region E_{1} of the parameter space where the Laplace state P_{1} 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 E_{1} is highlighted with a black contour, obtained using the closedform expression in Eqs. (21) and (22). Region E_{1} is entirely contained inside the dashed brown lines, whose left and right limits r_{1} and r_{2} are given in Eq. (20), and whose bottom and top limits are given by Eq. (23). The Laplace state P_{1} is singular atpoint S_{1}, in which is discontinuous. 

In the text 
Fig. 9 Eccentricity and inclination of the satellite at the eccentric equilibrium P. The threedimensional 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 threedimensional surface. The colour is the same as in Fig. 9, top panel. In Fig. 9, the top tubelike 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 S_{1} and S are labelled, and the cutting line of the threedimensional surface is shown by a grey line (same as Fig. 9). Threedimensional figures are provided in Appendix C, where the geometry of the stable regions is clearer. 

In the text 
Fig. 12 Spinaxis dynamics in the vicinity of a firstorder 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 xyplane. The obliquity ε is the tilt from the zaxis, and the resonance angle σ is the polar angle measured in the xyplane. 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 Spinaxis 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 spinaxis precession frequency of the planet as a function of its obliquity. The curve is obtained by injecting a = r_{F} 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 Spinaxis 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 spinaxis 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 S_{1} and the ridge line a = r_{F} are indicated in red. 

In the text 
Fig. 18 Location and width of every firstorder 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 semimajor 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 = r_{F} dividing the closein and faraway satellite regimes is represented by a dashed blue line. The region E_{1} where the satellite’s Laplace state P_{1} is unstable isshown in black. The singular point S_{1} 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 s_{8}. 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 selfexplanatory. 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 Longterm dynamics of Titan as it migrates away from Saturn. Saturn is assumed to remain locked in secular spin–orbit resonance with s_{8} 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 semimajor 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 P_{1}, 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 threebody problem including Saturn’s oblateness. The output is then digitally filtered to remove its shortperiod 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 spinaxis 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 firstorder secular spin–orbit resonances, and the blue hatched area is the region E_{1} 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 = r_{M}. 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 E_{3} of the parameter space where the Laplace state P_{3} 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 E_{3} is obtained using the closedform expression in Eqs. (B.11) and (22). The droplike 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 closedform 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 threedimensional 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 tubelike 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 threedimensional surface describing the eccentricity of the satellite. As in Fig. C.2, the surface is seen from two viewing angles and the top tubelike portion has been cut for better readability. 

In the text 
Fig. C.4 Same as Fig. C.3 but replacing the coordinates a∕r_{M} and e by r_{M} ∕a and e^{2}, respectively.Contrary to previous figures, the top tubelike 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 semimajor axis a∕r_{M} 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 cos^{2} ε_{+} 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 tubelike 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 threedimensional surface describing the eccentricity of the satellite. As in Fig. C.7, the surface is seen from two viewing angles and the top tubelike portion has been cut for better readability. 

In the text 
Fig. C.9 Same as Fig. C.8 but replacing the coordinates a∕r_{M} and e by r_{M} ∕a and e^{2}, respectively.Contrary to previous figures, the top tubelike 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 selfconsistent 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.