Issue 
A&A
Volume 640, August 2020



Article Number  A11  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202038432  
Published online  31 July 2020 
The future large obliquity of Jupiter
^{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:
18
May
2020
Accepted:
3
June
2020
Aims. We aim to determine whether Jupiter’s obliquity is bound to remain exceptionally small in the Solar System, or if it could grow in the future and reach values comparable to those of the other giant planets.
Methods. The spinaxis of Jupiter is subject to the gravitational torques from its regular satellites and from the Sun. These torques evolve over time due to the longterm variations of its orbit and to the migration of its satellites. With numerical simulations, we explore the future evolution of Jupiter’s spin axis for different values of its moment of inertia and for different migration rates of its satellites. Analytical formulas show the location and properties of all relevant resonances.
Results. Because of the migration of the Galilean satellites, Jupiter’s obliquity is currently increasing, as it adiabatically follows the drift of a secular spinorbit resonance with the nodal precession mode of Uranus. Using the current estimates of the migration rate of the satellites, the obliquity of Jupiter can reach values ranging from 6° to 37° after 5 Gyr from now, according to the precise value of its polar moment of inertia. A faster migration for the satellites would produce a larger increase in obliquity, as long as the drift remains adiabatic.
Conclusions. Despite its peculiarly small current value, the obliquity of Jupiter is no different from other obliquities in the Solar System: It is equally sensitive to secular spinorbit resonances and it will probably reach comparable values in the future.
Key words: planets and satellites: dynamical evolution and stability / planets and satellites: fundamental parameters
© M. Saillenfest et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://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
The obliquity of a planet is the angle between its spin axis and the normal to its orbit. A nonzero obliquity results in seasonal climate changes along the planet’s orbit, as occurs on Earth. In the protoplanetary disc, giant planets are expected to form with nearzero obliquities, while terrestrial planets should exhibit more random values (see e.g. Ward & Hamilton 2004; Rogoszinski & Hamilton 2020a). Yet, the planets of the Solar System all feature a large variety of obliquities. The case of Mercury is special because the strong tidal dissipation due to the proximity of the Sun now tightly maintains Mercury’s obliquity to a nearzero value (see e.g. Correia & Laskar 2010). Excluding Mercury, Jupiter is by far the planet of the Solar System that has the smallest obliquity (see Table 1). This small value seems to put Jupiter in a different category, and it appears unclear why Jupiter should be the only giant planet to indefinitely preserve its primordial obliquity.
Large obliquity changes can be produced by strong impacts. An impact with a planetarysized body is thought to have created the Moon and affected the spin axis of the Earth, which has remained unchanged ever since (Canup & Asphaug 2001; Li & Batygin 2014b). Largescale collisions have also probably participated in increasing the obliquity of Uranus (Boué & Laskar 2010; Morbidelli et al. 2012; Rogoszinski & Hamilton 2020a).
Apart from collisions, a wellknown mechanism that can modify the obliquity of a planet is a socalled “secular spinorbit resonance”, that is, a near commensurability between the frequency of precession of the spin axis and the frequency of one (or several) harmonics appearing in the precession of the orbit. This mechanism happens to be extremely common in planetary systems. The overlap of such resonances produces a large chaotic region for the spin axis of the terrestrial planets of the Solar System (see Laskar & Robutel 1993). This chaos probably had a strong influence on the early obliquity of Venus, which was then driven to its current value by the solar tides combined with its thick atmosphere (Correia & Laskar 2001, 2003; Correia et al. 2003). The Moon currently protects the Earth from large chaotic variations in its obliquity (Laskar et al. 1993; Li & Batygin 2014a), but due to tidal dissipation within the EarthMoon system, the Earth will eventually reach the chaotic region in a few gigayears from now (see Néron de Surgy & Laskar 1997). This chaotic zone also strongly affects the obliquity of Mars, which still currently wanders between 0° and more than 60° (Laskar et al. 2004a; Brasser & Walsh 2011). As shown by Millholland & Batygin (2019), secular spinorbit resonances can also take place very early in the history of a planet, that is, within the protoplanetary disc itself. More generally, secular spinorbit resonances are thought to strongly affect the obliquity of exoplanets (see e.g. Atobe et al. 2004; Brasser et al. 2014; Deitrick et al. 2018b,a; Shan & Li 2018; Millholland & Laughlin 2018, 2019; Quarles et al. 2019; Saillenfest et al. 2019; Kreyche et al. 2020).
For the giant planets of the Solar System, the secular spinorbit resonances are relatively thin today and well separated from each other. This is why it is so difficult to explain the large obliquity of Uranus by a spinorbit coupling, now that the precession of Uranus’ spin axis is far from any firstorder resonances (see e.g. Boué & Laskar 2010, Rogoszinski & Hamilton 2020a,b). Jupiter and Saturn, on the contrary, are located very close to strong resonances: Jupiter is close to resonance with the nodal precession mode of Uranus (Ward & Canup 2006), and Saturn is close to resonance with the nodal precession mode of Neptune (Ward & Hamilton 2004; Hamilton & Ward 2004; Boué et al. 2009). Therefore, the dynamics of Jupiter’s spin axis seems to be equally affected by secular spinorbit resonances as other planets in the Solar System. This was confirmed by Brasser & Lee (2015) and Vokrouhlický & Nesvorný (2015), who show that models of the late planetary migration have to be finely tuned to avoid overexciting Jupiter’s obliquity by spinorbit coupling, while tilting Saturn to its current orientation. In this regard, the spinaxis dynamics of Jupiter does not appear to be special at all, in contrast to its small obliquity value.
In this article, we aim to investigate the future longterm spinaxis dynamics of Jupiter. In particular, we want to determine whether Jupiter’s obliquity is bound to remain exceptionally small in the Solar System, or if it could grow in the future and reach values comparable to those of the other planets.
The precession motion of a planet’s spin axis depends on the physical properties of the planet (mass repartition and spin velocity), but also on external torques applied to its equatorial bulge. These torques come from the combined gravitational attraction of the Sun and of satellites (if it has any). Since the orbit of Jupiter is stable over billions of years (Laskar 1990), the direct torque from the Sun will not noticeably change in the future. However, Jupiter’s satellites are known to migrate over time because of tidal dissipation. The future longterm orbital evolution of the Galilean satellites has been recently explored by Lari et al. (2020). The solutions that they describe can therefore be used as a guide to study the future spinaxis dynamics of Jupiter. Due to their much smaller masses, the other satellites of Jupiter do not contribute noticeably to its spinaxis dynamics.
In Sect. 2, we describe our dynamical model and discuss the range of acceptable values for the physical parameters of Jupiter, in particular its polar moment of inertia. In Sect. 3, we present our results about the future spinaxis dynamics of Jupiter: We explore the outcomes given by different values of the poorly known physical parameters of Jupiter and by different migration rates for its satellites. Our conclusions are summarised in Sect. 4.
Current obliquities of the planets of the Solar System.
2 Secular dynamics of the spin axis
2.1 Equations of motion
The spinaxis dynamics of an oblate planet subject to the lowestorder term of the torque from the Sun is given for instance by Laskar & Robutel (1993) or Néron de Surgy & Laskar (1997). Far from spinorbit resonances, and due to the weakness of the torque, the longterm evolution of the spinaxis is accurately described by the secular Hamiltonian function (i.e. averaged over rotational and orbital motions). This Hamiltonian can be written (1)
where the conjugate coordinates are X (cosine of obliquity) and − ψ (minus the precession angle). The Hamiltonian in Eq. (1) depends explicitly on time t through the orbital eccentricity e and through the functions (2)
In these expressions, q = η cos Ω and p = η sin Ω, where η ≡ sin(I∕2), and I and Ω are the orbital inclination and the longitude of ascending node of the planet, respectively. If the orbit of the planet is fixed in time, its obliquity is constant and its precession angle ψ circulates with constant angular velocity . The quantity α is called the precession constant. It depends on the spin rate of the planet and of its mass distribution, through the formula: (3)
where is the gravitational constant, m_{⊙} is the mass of the Sun, ω is the spin rate of the planet, a is its semimajor axis, J_{2} is its second zonal gravity coefficient, and λ is its normalised polar moment of inertia. We retrieve the expression given for instance by Néron de Surgy & Laskar (1997) by noting that (4)
where A, B, and C are the equatorial and polar moments of inertia of the planet, M is its mass, and R_{eq} is its equatorial radius.
The precession rate of the planet is increased if it possesses massive satellites. If the satellites are far away from the planet, their equilibrium orbital plane (called Laplace plane, see Tremaine et al. 2009) is close to the orbital plane of the planet; therefore, faraway satellites increase the torque exerted by the Sun on the equatorial bulge of the planet. If the satellites are close to the planet, on the contrary, their equilibrium orbital plane coincides with the equator of the planet and precesses with it as a whole (Goldreich 1965); therefore, closein satellites artificially increase the oblateness and the rotational angular momentum of the planet. In the closein satellite regime, an expression for the effective precession constant has been derived by Ward (1975). As detailed by French et al. (1993), it consists in replacing J_{2} and λ in Eq. (3) by the effective values: (5)
where m_{k}, a_{k}, and n_{k} are the mass, the semimajor axis, and the mean motion of the kth satellite. In these expressions, the eccentricities and inclinations of the satellites are neglected. This approximation has been widely used in the literature. In the case of a single satellite, Boué & Laskar (2006) have obtained a general expression for the precession rate of a planet with an eccentric and inclined satellite, encompassing both the closein and faraway regimes. Using their article, we can verify that the Galilean satellites are in the closein regime. The Laplace plane of Callisto is inclined today by less than 1° with respect to Jupiter’s equator. The small eccentricities and inclinations of the Galilean satellites would contribute to and λ′ with terms of order and , so even if e_{k} increases up to 0.1 (a value found by Lari et al. 2020 in some cases) or if I_{k} increases up to 10°, the additional contribution to and λ′ would only be of order 10^{−4} and 10^{−6}, respectively. As we see below, this contribution is much smaller than our uncertainty on the value of λ, allowing us to stick to the approximation given by Eq. (5).
2.2 Orbital solution
The Hamiltonian given in Eq. (1) depends on the orbit of the planet and on its temporal variations. In order to explore the longterm dynamics of Jupiter’s spin axis, we need an orbital solution that is valid over billions of years. This is well beyond the timespan covered by ephemerides. Luckily, the orbital dynamics of the giant planets of the Solar System are almost integrable and excellent solutions have been developed. We use the secular solution of Laskar (1990) expanded in quasiperiodic series: (6)
where ϖ is Jupiter’s longitude of perihelion. The amplitudes E_{k} and S_{k} are real constants, and the angles θ_{k} and ϕ_{k} evolve linearly overtime t, with frequencies μ_{k} and ν_{k}: (7)
The complete orbital solution of Laskar (1990) can be found in Appendix A for amplitudes down to 10^{−8}.
The series in Eq. (6) contain contributions from all the planets of the Solar System. In the integrable approximation, the frequency of each term corresponds to a unique combination of the fundamental frequencies of the system, usually noted g_{j} and s_{j}. In the limit of small masses, small eccentricities and small inclinations (LagrangeLaplace secular system), the z series only contains the frequencies g_{j}, while the ζ series only contains the frequencies s_{j} (see e.g. Murray & Dermott 1999 or Laskar et al. 2012). This is not the case in more realistic situations, as recalled for instance by Kreyche et al. (2020) in the context of obliquity dynamics. In planetary systems featuring meanmotion resonances, the spinaxis of a planet can be affected by shifted orbital precession frequencies (Millholland & Laughlin 2019) or by secondary resonances (Quillen et al. 2017, 2018). However, this does not apply in the Solar System as it is today, even when the existing near commensurabilities (like the “great Jupiter–Saturn inequality”) are taken into account. Table 2 shows the combinations of fundamental frequencies identified for the largest terms of Jupiter’s ζ series obtained by Laskar (1990).
As explained by Saillenfest et al. (2019), at first order in the amplitudes S_{k} and E_{k}, secular spinorbit resonant angles can only be of the form σ_{p} = ψ + ϕ_{p}, where p is a given index in the ζ series. Resonances featuring terms of the z series only appear at third order and beyond. For the terrestrial planets of the Solar System, the z and ζ series converge very slowly, which implies that large resonances are very numerous. These resonances overlap massively and produce wide chaotic zones in the obliquity dynamics (see Laskar & Robutel 1993; Néron de Surgy & Laskar 1997; Correia et al. 2003; Laskar et al. 2004a). The situation is very different for the giant planets of the Solar System, for which the z and ζ series converge quickly owing to the quasiintegrable nature of their dynamics. Therefore, the secular spinorbit resonances are small and isolated from each other, and only firstorder resonances play a substantial role.
Figure 1 shows the location and width of every firstorder resonance for the spinaxis of Jupiter in an interval of precession constant α ranging from 0 to 5″ yr^{−1}. Because of the chaotic dynamics of the Solar System (Laskar 1989), the fundamental frequencies related to the terrestrial planets (e.g. s_{1}, s_{2}, and γ appearing in Table 2) could vary substantially over billions of years (Laskar 1990). However, they only marginally contribute to Jupiter’s orbital solution and none of them takes part in the resonances shown in Fig. 1. Our secular orbital solution of Jupiter can therefore be considered valid over a billionyear timescale.
First twenty terms of Jupiter’s inclination and longitude of ascending node in the J2000 equator and equinox reference frame.
2.3 Precession constant
As shown by the Hamiltonian function in Eq. (1), the precession constant α is a key parameter of the spinaxis dynamics of a planet. Among the physical parameters of Jupiter that enter into its expression (see Eq. (3)), all are very well constrained from observations except the normalised polar moment of inertia λ.
While comparing the values of λ given in the literature, one must be careful about the normalisation used. Equation (4) explicitly requires a normalisation using the equatorial radius R_{eq}, since it is linked to the value of J_{2}. However, published values of the polar moment of inertia are often normalised using the mean radius of Jupiter, which differs from R_{eq} by a factor of about 0.978. This distinction seems to have been missed by Ward & Canup (2006), who quote the nominal value given by D. R. Williams in the NASA Jupiter fact sheet^{1} as 0.254, whereas it actually translates into λ = 0.243 when it is normalised using R_{eq}. Ward & Canup (2006) also mention that “theoretical values [of λ] range from 0.255 for the extreme of a constantdensity core and massless envelope to 0.221 for a constantdensity envelope and pointmass core”. Unfortunately, these numbers are taken from a conference talk given by W. B. Hubbard in 2005 so we cannot check how they have been obtained. Since Eq. (3) is used, however, we can assume that they have been properly normalised using R_{eq}.
As is shown in Fig. 1, the spinaxis of Jupiter is located very close to a strong secular spinorbit resonance. The corresponding term of the orbital series is related to the precession mode of Uranus (term k = 4 in Table 2), and the resonant angle is σ_{4} = ψ + ϕ_{4}. As noted by Ward & Canup (2006), dissipative processes during the early planetary evolution are expected to have forced Jupiter’s spin axis to spiral down towards the centre of the resonance, called Cassini state 2. And indeed, the current value of σ_{4} is very close to zero, which has a low probability to happen if Jupiter is far from Cassini state 2 because σ_{4} would then circulate between 0° and 360°. In order to match Cassini state 2, however, Jupiter’s normalised moment of inertia should be λ ≈ 0.2365 (see Fig. 2). Since this value is not far from what is proposed in the literature, this prompted Ward & Canup (2006) to consider this value as likely for Jupiter.
As noted by Le Maistre et al. (2016), the value of λ ≈ 0.2365 corresponds to a massive core for Jupiter, and estimates obtained from models of Jupiter’s interior structure are generally higher. Helled et al. (2011) obtain values of λ ranging from 0.251 to 0.253, that were confirmed by Nettelmann et al. (2012). These values are consistent with the range of λ ∈ [0.221, 0.255] quoted above. Other studies seem to agree on even higher values: Wahl et al. (2017) and Ni (2018) present values of λ ranging between 0.2629 and 0.2644, compatible with the findings of Hubbard & Marley (1989), Nettelmann et al. (2012), and Hubbard & Militzer (2016). Finally, both low and high values are obtained by Vazan et al. (2016), who give either λ = 0.247 or λ = 0.262 for three different models. As explained by Le Maistre et al. (2016), however, all these values are modeldependent and still a matter of debate. Hopefully, the Juno mission will provide direct observational constraints soon that will help us to determine which models of Jupiter’s interior structure are the most relevant.
Here, instead of relying on one particular value of λ, we turn to the exploration of the whole range of values given in the literature, namely λ ∈ [0.220, 0.265]. The rotation velocity of Jupiter is taken from Archinal et al. (2018) and the other physical parameters are fixed to those used by Lari et al. (2020) for consistency with the satellites’ orbital evolution (see below). The corresponding value for the current precession constant of Jupiter, computed from Eqs. (3) and (5), ranges from 2.64 to 3.17″ yr^{−1}. Given this large range, using updated physical parameters (see e.g. Folkner et al. 2017; Iess et al. 2018; Serra et al. 2019) would only slightly shift the value of α within our exploration interval.
Because of tidal dissipation, satellites slowly migrate over time. This produces a drift of the precession constant α on a timescale that is much larger than the precession motion (i.e. the circulation of ψ). The longterm spinaxis dynamics of a planet with migrating satellites is described by the Hamiltonian in Eq. (1), but where α is a slowlyvarying function of time. In the EarthMoon system, the outward migration of the Moon produces a decrease of α that pushes the Earth towards a wide chaotic region (see Néron de Surgy & Laskar 1997). This decrease of α is due to the fact that the Moon is in the farsatellite regime (see Boué & Laskar 2006). The Galilean satellites, on the contrary, are in the closesatellite regime, and their outward migration produces an increase of α, as shown by Eq. (5). This increase can be quantified using the longterm orbital solution of Lari et al. (2020) depicted in Fig. 3 and interpolating between data points. The result is presented in Fig. 4 for the two extreme values of λ considered in this article, as well as for the value of λ ≈ 0.2365 proposed by Ward & Canup (2006). Despite the various outcomes of the dynamics described by Lari et al. (2020), the result on the evolution of α is almost undistinguishable from one of their simulations to another, even if the eccentricities of the satellites are taken into account in Eq. (5). Indeed, the variation of α mostly depends on the drift of the satellites’ semimajor axes, which is almost identical in every simulation of Lari et al. (2020).
Since the rate of energy dissipation between Jupiter and its satellites is not well known today, the timescale of the drift shown in Figs. 3 and 4 could somewhat contract or expand. This point is further discussed in Sect. 3. Moreover, other parameters in Eq. (3) probably slightly vary over billions of years, such as the spin velocity of Jupiter or its oblateness. We consider that the impact of their variations on the value of α is small and contained within our exploration range.
Fig. 1 Location and width of every firstorder secular spinorbit resonance for Jupiter. Each resonant angle is of the form σ_{p} = ψ + ϕ_{p} where ϕ_{p} has frequency ν_{p} labelled on the graph according to its index in the orbital series (see Table 2 and Appendix A). For a given value of the precession constant α, the interval of obliquity enclosed by the separatrix is shown in pink, as computed using the exact formulas given by Saillenfest et al. (2019). The green bar on the left shows Jupiter’s current obliquity and the range for its precession constant considered in this article, as detailed in Sects. 2.3 and 2.4. 
Fig. 2 Trajectory of Jupiter’s spin axis in the vicinity of resonance with the fourth harmonics of ζ (see Table 2). Being farther away from Jupiter’s precession frequency, the contribution of other harmonics can be averaged; their mean contribution is included here up to third order in the amplitudes (as in Eq. (17) of Saillenfest et al. 2019). Each trajectory corresponds to a level curve of the Hamiltonian, which has only one degree of freedom. The red dot shows the current location of Jupiter, and the black dot shows Cassini state 2. The red curve is the current trajectory of Jupiter’s spin axis for λ = 0.250 (top) or λ = 0.237 (bottom). 
Fig. 3 Typical evolution of the semimajor axes of the Galilean satellites obtained by Lari et al. (2020). The values are expressed in unit of Jupiter’s equatorial radius. The bump at about 1.8 Gyr is due to the capture of Callisto into resonance. 
Fig. 4 Evolution of the effective precession constant of Jupiter due to the migration of its satellites. The top and bottom curves correspond to the two extreme values of the normalised polar moment of inertia λ considered in this article. They appear into α through Eq. (3). The central curve corresponds to the value of λ that places Jupiter just near Cassini state 2 with the precession mode of Uranus (Ward & Canup 2006). 
2.4 Initial conditions
The initial orientation of the spin axis is taken from the solution of Archinal et al. (2018) averaged over shortperiod terms. At the level of precision required by our exploratory study, the refined orientation obtained by Durante et al. (2020) is undistinguishable from this nominalorientation. With respect to Jupiter’s secular orbital solution (see Sect. 2.2), this gives an obliquity ε = 3.120° and a precession angle ψ = −137.223° at time J2000. The uncertainty on these values is extremely small compared to the range of α considered (see Sect. 2.3). Since the uncertainty is smaller than the curve width of our figures, we do not consider any error bar on the initial value of ε and ψ.
3 Obliquity evolution with migrating satellites
For values of λ finely sampled in our exploration interval, the spin axis of Jupiter is numerically propagated forwards in time for 5 Gyr. By virtue of trigonometric identities, moving Jupiter’s orbit one step forwards in time using the quasiperiodic decomposition in Eq. (6) only amounts to computing a few sums and products. The trajectories obtained are shown in Fig. 5 for a few values of λ. They are projected in the plane of the obliquity and the precession constant of Jupiter, where we localise also the centres and widths of allfirstorder secular spinorbit resonances. See Appendix B for further details about the geometry of the resonances.
For values of λ smaller than about 0.228, Jupiter starts outside of the large resonance with ϕ_{4}, and the increase of its precession constant α pushes it even farther away over time. As shown by the trajectory computed for λ = 0.227, the crossing of the very thin resonance with ϕ_{19} twists the trajectory a little, but this cannot produce any large change of obliquity. Indeed the resonance with ϕ_{19} is not strong enough to capture Jupiter’s spin axis: It is crossed quickly as α increases, and Fig. 6shows that the libration period of σ_{19} = ψ + ϕ_{19} is very long. This results in a nonadiabatic crossing (see Appendix C for details). Consequently, no major obliquity variation for Jupiter can be expected in the future if λ < 0.228. However, such small values of λ seem to be ruled out by most models of Jupiter’s interior (see Sect. 2.3).
For values of λ larger than 0.228, on the contrary, Jupiter is currently located inside or below the large resonance with ϕ_{4}. As predicted, the value λ = 0.2365 results in very small oscillations around Cassini state 2. As its precession constant α slowly increases with time, Jupiter is captured into the resonance and follows the drift of its centre towards large obliquities. Indeed, the resonance with ϕ_{4} is large, and the libration period of σ_{4} = ψ + ϕ_{4} is short compared to the variation timescale of α (see Figs. 5 and 6). This results in an adiabatic capture. The various possible outcomes of adiabatic and nonadiabatic crossings of secular spinorbit resonances have recently been studied by Su & Lai (2020). However, the orbital motion is here not limited to a single harmonic, and Appendix C shows that the separatrix of the resonance is replaced by a chaotic “moat”. Properly speaking, the resonance with ϕ_{4} becomes a “true resonance” only as soon as the separatrix appears, that is, for α larger than α ≈3.04″ yr^{−1} (see Appendix B). In the whole range of values of λ > 0.228 considered in this article, the spinaxis of Jupiter is initially located close enough to Cassini state 2 to invariably end up inside the separatrix of the resonance when it appears. The capture probability is therefore 100%. None of our simulation shows a release out of resonance or a turnoff towards Cassini state 1, which could have been a possible outcome if Jupiter’s spin axis was initially located farther away from Cassini state 2 or if the drift of α was not adiabatic^{2}. Since in canonical coordinates the resonance width increases for α growing up to 4.244″ yr^{−1}, no separatrix crossing can happen, even for a large libration amplitude inside the resonance (e.g. for λ = 0.228 in Fig. 5). The maximum obliquity reached by Jupiter is therefore only limited by the finite amount of time considered.
If the Galilean satellites migrate faster than shown in Fig. 3, the obliquity reached in 5 Gyr would be larger than that presented in Fig. 5. The migration rate of the satellites is not well known. According to Lari et al. (2020), the longterm migration rate of the satellites varies by ± 15% over the uncertainty range of the parameter measured by Lainey et al. (2009). This parameter quantifies the dissipation within Jupiter at Io’s frequency. Figure 7 shows the maximum obliquity reached in 5 Gyr for λ sampled in our exploration interval and sampled in its uncertainty range. We retrieve the discontinuity at λ ≈ 0.228 discussed before, below which only small obliquity variations are possible. For λ > 0.228, as expected, we see that a fast migration and a small moment of inertia produce a fast increase of obliquity, which reaches 37° in 5 Gyr in the most favourable case of Fig. 7. On the contrary, a slow migration and a large moment of inertia produce a slow increase of obliquity, which barely reaches 6° in 5 Gyr in the most unfavourable case of Fig. 7.
Fig. 5 Future evolution of Jupiter’s spin axis projected in the plane of the obliquity and the precession constant α. Each panel corresponds to a value of the normalised polar moment of inertia of Jupiter given in title. The green bar shows the initial location of Jupiter’s spin axis according to our exploration interval of λ; the central mark is the value proposed by Ward & Canup (2006). The red curves show the centre of all firstorder secular spinorbit resonances (Cassini state 2) and the coloured areas represent their widths (same as Fig. 1). From bottom to top, the resonances are with ϕ_{6}, with ϕ_{4}, and with ϕ_{19} (see Table 2). The black dots show the numerical solutions obtained over a timespan of 5 Gyr from now; they evolve from bottom to top. According to the exact migration rate of the Galilean satellites, the timeline could somewhat contract or expand (see text). 
Fig. 6 Period of small oscillations about the resonance centre for a resonance with ϕ_{4} or ϕ_{19}. Even though complete closedform solutions exist (see Haponiak et al. 2020), the smalloscillation limit leads to handier formulas, suitable for order of magnitude estimates. The resonant angles are σ_{4} = ψ + ϕ_{4} and σ_{19} = ψ + ϕ_{19}, respectively.Dashed curves are used for oscillations about Cassini state 2 before the separatrix appears. The appearance of the separatrix is marked by a blue dot. 
Fig. 7 Maximum obliquity reached by Jupiter after 5 Gyr from now as a function of its normalised polar moment of inertia (right vertical axis) and the dissipation parameter of Jupiter at Io’s frequency (horizontal axis). The left vertical axis shows the current precession constant α of Jupiter. Some level curves are shown in red. 
4 Discussion and conclusion
Prompted by the peculiarly small value of the current obliquity of Jupiter, we studied the future longterm evolution of its spinaxis under the influence of its slowly migrating satellites.
Jupiter is located today near a strong secular spinorbit resonance with the nodal precession mode of Uranus (Ward & Canup 2006). Because of this resonance, the obliquity of Jupiter is found to be currently increasing, provided that its normalised polar moment of inertia is larger than about 0.228. Such a small value seems to be ruled out by models of Jupiter’s interior (see e.g. Helled et al. 2011; Hubbard & Militzer 2016; Wahl et al. 2017). For larger valuesof λ, the migration of the Galilean satellites induces an adiabatic drift of the precession constant α of Jupiter that pushes its spin axis inside the resonance and forces it to follow the resonance centre towards high obliquities. For the value λ ≈ 0.2365 proposed by Ward & Canup (2006), the obliquity can reach values as large as 30° in the next 5 Gyr. For the value λ ≈ 0.252 obtained byHelled et al. (2011), the obliquity reaches values ranging from about 17° to 23°. The increase is more modest for values close to λ ≈ 0.264 found by other authors, for which the maximum value of the obliquity ranges from about 6° to 17°. Hence, our main conclusion is that, contrary to Saturn, Jupiter did not have time to tilt much yet from its primordial orientation, but it will in the future and possibly a lot.
The model of tidal dissipation applied by Lari et al. (2020) to the Galilean satellites and used here to compute the drift of α is simplified. The current migration rates of satellites in the Solar System have been proved to be higher than previously thought (see Lainey et al. 2009, 2017). As discussed by Lari et al. (2020), the migration of the Galilean satellites could be even faster than considered here if ever one of the outer satellites was pushed by a resonance with the frequency of an internal oscillation of Jupiter (Fuller et al. 2016). This would result in a faster increase of Jupiter’s obliquity. This increase would be halted, however, if the satellites ever migrate so fast as to break the adiabaticity of the capture into secular spinorbit resonance. In this case, Jupiter would cross the resonance and exit without following the drift of its centre (see e.g. Ward & Hamilton 2004; Su & Lai 2020). Numerical experiments show that adiabaticity would be broken for a migration more than 110 times faster than currently estimated. Such an extremely fast migration seems unlikely. Moreover, with such a fast migration, Callisto and then Ganymede would soon go beyond the closesatellite regime (Boué & Laskar 2006): This would slow down the increase of α and possibly restore the adiabaticity of its drift. Therefore, the future increase of Jupiter’s obliquity appears to be a robust result.
The maximum obliquity that Jupiter will reach could be very large, but it depends on the precise value of Jupiter’s polar moment of inertia and on the precise migration rate of the Galilean satellites. We hope to obtain soon new estimates for these two crucial parameters, in particular from the results of the Juno and JUICE missions.
Acknowledgements
We thank Marco Fenucci for his help and his suggestions during the redaction of our manuscript. We also thank the anonymous referee for her/his valuable comments. G.L. acknowledges financial support from the Italian Space Agency (ASI) through agreement 201740H.0 in the context of the NASA Juno mission.
Appendix A Orbital solution for Jupiter
The secular orbital solution of Laskar (1990) is obtained by multiplying the normalised proper modes and (Tables VI and VII of Laskar 1990) by the matrix corresponding to the linear part of the solution (Table V of Laskar 1990). In the series obtained, the terms with the same combination of frequencies are then merged together, resulting in 56 terms in eccentricity and 60 terms in inclination. This forms the secular part of the orbital solution of Jupiter, which is what is required by our averaged model.
The orbital solution is expressed in the variables z and ζ as describedin Eqs. (6) and (7). In Tables A.1 and A.2, we give the terms of the solution in the J2000 ecliptic and equinox reference frame for amplitudes down to 10^{−8}.
Quasiperiodic decomposition of Jupiter’s eccentricity and longitude of perihelion (variable z).
Quasiperiodic decomposition of Jupiter’s inclination and longitude of ascending node (variable ζ).
Appendix B Crossing the resonance with ϕ4
Figures 1 and 5 show the location and width of all firstorder secular spinorbit resonances produced by Jupiter’s orbital solution (Appendix A). In particular, Jupiter is located very close to the large resonance with ϕ_{4}, whose frequency is the nodal precession mode of Uranus. Figures B.1 and B.2 show the geometry of the resonance with ϕ_{4} for different values of the precession constant α of Jupiter. These graphs can be understood as horizontal sections of Fig. 5, where we can locate the centre of the resonance (i.e. Cassini state 2) and the separatrix width. For easier comparison, Figs. 5 and B.1 share the same horizontal axis.
Fig. B.1 Level curves of the Hamiltonian function in the vicinity of resonance with ϕ_{4}. The resonant angle is σ_{4} = ψ + ϕ_{4}. Other terms are averaged and included up to the third order of their amplitudes (see Saillenfest et al. 2019). Each panel corresponds to a different value of the precession constant α. Equilibriumpoints (called Cassini states) are shown by black spots. The interior of the resonance is coloured red and the separatrix is shown with a thick red curve. The location and width of the resonance for continuous values of α can be seen in Fig. 5. In order to avoid being misled by coordinate singularities, Fig. B.2 shows the same level curves in a different set of coordinates. 
Fig. B.2 Same as Fig. B.1, but using polar coordinates that are not singular for an obliquity ε = 0°. The outer black circle corresponds to an obliquity ε = 40°. 
Appendix C Crossing the resonance with ϕ_{19}
As a matter of fact, Jupiter’s orbital motion is not restricted to the ϕ_{4} term. However, secular spinorbit resonances with all other terms (apart from ϕ_{19}) are located very far from the location of Jupiter, so that their effects average over time. The case of ϕ_{19} is special: even thoughit is very thin, this resonance is not far from Jupiter’s location (see the upper red curve in Fig. 5), which means that ψ + ϕ_{19} is a slow angle that cannot be averaged out.
Instead of considering only ϕ_{4}, as in Appendix B, a more rigorous model of the longterm spinaxis dynamics of Jupiter consists in averaging the Hamiltonian function over all angles except resonances with both ϕ_{4} and ϕ_{19}. For a constant value of α, this results in a twodegreeoffreedom Hamiltonian system, in which the two angle coordinates are σ_{4} = ψ + ϕ_{4} and σ_{19} = ψ + ϕ_{19}. The dynamics can then be studied using Poincaré surfaces of section. Figure C.1 shows two examples of sections. The lower island centred at σ_{19} = 0 corresponds to the thin resonance with ϕ_{19}: as expected,it is completely distorted as compared to the unperturbed separatrix (blue curve) due to the proximity of the large resonance with ϕ_{4}. It still persists, however, as a set of periodic orbits. In contrast, the large resonance with ϕ_{4} is hardly affected at all by the ϕ_{19} term, which only transforms its separatrix into a thin chaotic belt. In the left panel of Fig. C.1, we can also recognise Cassini state 1 with ϕ_{4} (for σ_{4} = π and a small obliquity), that is also visible in Fig. B.1.
Fig. C.1 Poincaré surfaces of section showing the dynamics in the vicinity of resonances with ϕ_{4} and ϕ_{19}. Each graph corresponds to a different value of α (see titles). The separatrices of the two resonances taken separately are shown with coloured curves: red for ϕ_{4} and blue for ϕ_{19}. 
We investigated whether Jupiter could be trapped into the thin resonance with ϕ_{19} and follow its resonance centre, but we found out that this can never happen. On the one hand, the current phase of σ_{19} is close to π, so that even if λ is finely tuned to place Jupiter right inside the resonance, it ends up near the separatrix, leading to an unstable resonant motion. On the other hand, as shown in Fig. 6, the libration period of σ_{19} is extremelylarge (the width and oscillation frequency both scale as the square root of the amplitude S_{19}). This means that, as α increases, the crossing of this resonance is not adiabatic. The libration periods shown in Fig. 6 should be compared to the time needed for α to go through the resonant region. According to Fig. 4, the mean increase rate of α is 0.086″ yr^{−1} Gyr^{−1} and according to Fig. 5, the resonances have a vertical width Δα ≈ 0.21″ yr^{−1} for the ϕ_{4} resonance and Δα ≈ 0.01″ yr^{−1} for the ϕ_{19} resonance (computed at the right separatrix when it appears). Therefore, the time that would be needed for α to cross the ϕ_{4} resonance is Δt ≈ 2.5 Gyr, which corresponds to many oscillation periods of σ_{4} (about 100): this is the adiabatic regime. On the contrary, the time needed for α to cross the ϕ_{19} resonance is Δt ≈ 0.1 Gyr, which corresponds to less than one oscillation period of σ_{19} (about 0.2): this is the nonadiabatic regime. As a result, a resonance capture with ϕ_{19} is extremely unlikely, even if the orbital motion of Jupiter was restricted to its 19th harmonic: Jupiter’s spin axis enters the resonance and exits before σ_{19} has time to oscillate. With suitable initial conditions, Jupiter roughly follows the resonance centre during the crossing, producing a bump in the obliquity evolution (see Fig. 5 for λ = 0.227), but nothing more can possibly happen. This kind of nonadiabatic resonance crossing is described by Ward et al. (1976), Laskar et al. (2004b), and Ward & Hamilton (2004) using Fresnel integrals.
References
 Archinal, B. A., Acton, C. H., A’Hearn, M. F., et al. 2018, Celest. Mech. Dyn. Astron., 130, 22 [Google Scholar]
 Atobe, K., Ida, S., & Ito, T. 2004, Icarus, 168, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Boué, G., & Laskar, J. 2006, Icarus, 185, 312 [NASA ADS] [CrossRef] [Google Scholar]
 Boué, G., & Laskar, J. 2010, ApJ, 712, L44 [NASA ADS] [CrossRef] [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]
 Brasser, R., & Walsh, K. J. 2011, Icarus, 213, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Brasser, R., Ida, S., & Kokubo, E. 2014, MNRAS, 440, 3685 [NASA ADS] [CrossRef] [Google Scholar]
 Canup, R. M., & Asphaug, E. 2001, Nature, 412, 708 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2003, Icarus, 163, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010, Icarus, 205, 338 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Laskar, J., & de Surgy, O. N. 2003, Icarus, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Deitrick, R., Barnes, R., Bitz, C., et al. 2018a, AJ, 155, 266 [CrossRef] [Google Scholar]
 Deitrick, R., Barnes, R., Quinn, T. R., et al. 2018b, AJ, 155, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572 [CrossRef] [Google Scholar]
 Folkner, W. M., Iess, L., Anderson, J. D., et al. 2017, Geophys. Res. Lett., 44, 4694 [NASA ADS] [CrossRef] [Google Scholar]
 French, R. G., Nicholson, P. D., Cooke, M. L., et al. 1993, Icarus, 103, 163 [CrossRef] [Google Scholar]
 Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P. 1965, AJ, 70, 5 [NASA ADS] [CrossRef] [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]
 Helled, R., Anderson, J. D., Schubert, G., & Stevenson, D. J. 2011, Icarus, 216, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Hubbard, W. B., & Marley, M. S. 1989, Icarus, 78, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Hubbard, W. B., & Militzer, B. 2016, ApJ, 820, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Konopliv, A. S., Park, R. S., & Ermakov, A. I. 2020, Icarus, 335, 113386 [CrossRef] [Google Scholar]
 Kreyche, S. M., Barnes, J. W., Quarles, B. L., et al. 2020, Planet. Sci. J., 1, 8 [CrossRef] [Google Scholar]
 Lainey, V., Arlot, J.E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Lari, G., Saillenfest, M., & Fenucci, M. 2020, A&A, 639, A40 [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J. 1989, Nature, 338, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Robutel, P. 1993, Nature, 361, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Joutel, F., & Robutel, P. 1993, Nature, 361, 615 [Google Scholar]
 Laskar, J., Correia, A. C. M., Gastineau, M., et al. 2004a, Icarus, 170, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004b, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., Boué, G., & Correia, A. C. M. 2012, A&A, 538, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Le Maistre, S., Folkner, W. M., Jacobson, R. A., & Serra, D. 2016, Planet. Space Sci., 126, 78 [CrossRef] [Google Scholar]
 Li, G., & Batygin, K. 2014a, ApJ, 790, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Li, G., & Batygin, K. 2014b, ApJ, 795, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Millholland, S., & Batygin, K. 2019, ApJ, 876, 119 [CrossRef] [Google Scholar]
 Millholland, S., & Laughlin, G. 2018, ApJ, 869, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Millholland, S., & Laughlin, G. 2019, Nat. Astron., 3, 424 [CrossRef] [Google Scholar]
 Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, Icarus, 219, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Néron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [Google Scholar]
 Nettelmann, N., Becker, A., Holst, B., & Redmer, R. 2012, ApJ, 750, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Ni, D. 2018, A&A, 613, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peale, S. J. 1969, AJ, 74, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Quarles, B., Li, G., & Lissauer, J. J. 2019, ApJ, 886, 56 [CrossRef] [Google Scholar]
 Quillen, A. C., Chen, Y.Y., Noyelles, B., & Loane, S. 2018, Celest. Mech. Dyn. Astron., 130, 11 [CrossRef] [Google Scholar]
 Quillen, A. C., NicholsFleming, F., Chen, Y.Y., & Noyelles, B. 2017, Icarus, 293, 94 [CrossRef] [Google Scholar]
 Rogoszinski, Z., & Hamilton, D. P. 2020a, ApJ, 888, 60 [CrossRef] [Google Scholar]
 Rogoszinski, Z., & Hamilton, D. P. 2020b, AAS J., submitted [arXiv:2004.14913] [Google Scholar]
 Saillenfest, M., Laskar, J., & Boué, G. 2019, A&A, 623, A4 [CrossRef] [EDP Sciences] [Google Scholar]
 Serra, D., Lari, G., Tommei, G., et al. 2019, MNRAS, 490, 1 [CrossRef] [Google Scholar]
 Shan, Y., & Li, G. 2018, AJ, 155, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Su, Y., & Lai, D. 2020, ApJ, submitted [arXiv:2004.14380] [Google Scholar]
 Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706 [NASA ADS] [CrossRef] [Google Scholar]
 Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016, ApJ, 829, 118 [Google Scholar]
 Vokrouhlický, D., & Nesvorný, D. 2015, ApJ, 806, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 [NASA ADS] [CrossRef] [Google Scholar]
 Ward, W. R. 1975, AJ, 80, 64 [NASA ADS] [CrossRef] [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]
 Ward, W. R., Colombo, G., & Franklin, F. A. 1976, Icarus, 28, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants, ed. T. J. Ahrens (Washington, DC: American Geophysical Union), 1 [Google Scholar]
There is a typographical error in Saillenfest et al. (2019): the list of the Cassini states given before Eq. (22) should read (4,2,3,1) instead of (1,2,3,4) in order to match the denomination introduced by Peale (1969).
All Tables
First twenty terms of Jupiter’s inclination and longitude of ascending node in the J2000 equator and equinox reference frame.
Quasiperiodic decomposition of Jupiter’s eccentricity and longitude of perihelion (variable z).
Quasiperiodic decomposition of Jupiter’s inclination and longitude of ascending node (variable ζ).
All Figures
Fig. 1 Location and width of every firstorder secular spinorbit resonance for Jupiter. Each resonant angle is of the form σ_{p} = ψ + ϕ_{p} where ϕ_{p} has frequency ν_{p} labelled on the graph according to its index in the orbital series (see Table 2 and Appendix A). For a given value of the precession constant α, the interval of obliquity enclosed by the separatrix is shown in pink, as computed using the exact formulas given by Saillenfest et al. (2019). The green bar on the left shows Jupiter’s current obliquity and the range for its precession constant considered in this article, as detailed in Sects. 2.3 and 2.4. 

In the text 
Fig. 2 Trajectory of Jupiter’s spin axis in the vicinity of resonance with the fourth harmonics of ζ (see Table 2). Being farther away from Jupiter’s precession frequency, the contribution of other harmonics can be averaged; their mean contribution is included here up to third order in the amplitudes (as in Eq. (17) of Saillenfest et al. 2019). Each trajectory corresponds to a level curve of the Hamiltonian, which has only one degree of freedom. The red dot shows the current location of Jupiter, and the black dot shows Cassini state 2. The red curve is the current trajectory of Jupiter’s spin axis for λ = 0.250 (top) or λ = 0.237 (bottom). 

In the text 
Fig. 3 Typical evolution of the semimajor axes of the Galilean satellites obtained by Lari et al. (2020). The values are expressed in unit of Jupiter’s equatorial radius. The bump at about 1.8 Gyr is due to the capture of Callisto into resonance. 

In the text 
Fig. 4 Evolution of the effective precession constant of Jupiter due to the migration of its satellites. The top and bottom curves correspond to the two extreme values of the normalised polar moment of inertia λ considered in this article. They appear into α through Eq. (3). The central curve corresponds to the value of λ that places Jupiter just near Cassini state 2 with the precession mode of Uranus (Ward & Canup 2006). 

In the text 
Fig. 5 Future evolution of Jupiter’s spin axis projected in the plane of the obliquity and the precession constant α. Each panel corresponds to a value of the normalised polar moment of inertia of Jupiter given in title. The green bar shows the initial location of Jupiter’s spin axis according to our exploration interval of λ; the central mark is the value proposed by Ward & Canup (2006). The red curves show the centre of all firstorder secular spinorbit resonances (Cassini state 2) and the coloured areas represent their widths (same as Fig. 1). From bottom to top, the resonances are with ϕ_{6}, with ϕ_{4}, and with ϕ_{19} (see Table 2). The black dots show the numerical solutions obtained over a timespan of 5 Gyr from now; they evolve from bottom to top. According to the exact migration rate of the Galilean satellites, the timeline could somewhat contract or expand (see text). 

In the text 
Fig. 6 Period of small oscillations about the resonance centre for a resonance with ϕ_{4} or ϕ_{19}. Even though complete closedform solutions exist (see Haponiak et al. 2020), the smalloscillation limit leads to handier formulas, suitable for order of magnitude estimates. The resonant angles are σ_{4} = ψ + ϕ_{4} and σ_{19} = ψ + ϕ_{19}, respectively.Dashed curves are used for oscillations about Cassini state 2 before the separatrix appears. The appearance of the separatrix is marked by a blue dot. 

In the text 
Fig. 7 Maximum obliquity reached by Jupiter after 5 Gyr from now as a function of its normalised polar moment of inertia (right vertical axis) and the dissipation parameter of Jupiter at Io’s frequency (horizontal axis). The left vertical axis shows the current precession constant α of Jupiter. Some level curves are shown in red. 

In the text 
Fig. B.1 Level curves of the Hamiltonian function in the vicinity of resonance with ϕ_{4}. The resonant angle is σ_{4} = ψ + ϕ_{4}. Other terms are averaged and included up to the third order of their amplitudes (see Saillenfest et al. 2019). Each panel corresponds to a different value of the precession constant α. Equilibriumpoints (called Cassini states) are shown by black spots. The interior of the resonance is coloured red and the separatrix is shown with a thick red curve. The location and width of the resonance for continuous values of α can be seen in Fig. 5. In order to avoid being misled by coordinate singularities, Fig. B.2 shows the same level curves in a different set of coordinates. 

In the text 
Fig. B.2 Same as Fig. B.1, but using polar coordinates that are not singular for an obliquity ε = 0°. The outer black circle corresponds to an obliquity ε = 40°. 

In the text 
Fig. C.1 Poincaré surfaces of section showing the dynamics in the vicinity of resonances with ϕ_{4} and ϕ_{19}. Each graph corresponds to a different value of α (see titles). The separatrices of the two resonances taken separately are shown with coloured curves: red for ϕ_{4} and blue for ϕ_{19}. 

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.