The future large obliquity of Jupiter (original) (raw)

A&A 640, A11 (2020)

1IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Université de Lille, 75014 Paris, France
e-mail: melaine.saillenfest@obspm.fr
2Department of Mathematics, University of Pisa, Largo Bruno Pontecorvo 5, 56127 Pisa, Italy

Received: 18 May 2020
Accepted: 3 June 2020

Abstract

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 spin-axis of Jupiter is subject to the gravitational torques from its regular satellites and from the Sun. These torques evolve over time due to the long-term 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 spin-orbit 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 spin-orbit 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 non-zero 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 near-zero 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 near-zero 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 planetary-sized 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). Large-scale 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 well-known mechanism that can modify the obliquity of a planet is a so-called “secular spin-orbit 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 Earth-Moon 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 spin-orbit resonances can also take place very early in the history of a planet, that is, within the protoplanetary disc itself. More generally, secular spin-orbit 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 spin-orbit 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 spin-orbit coupling, now that the precession of Uranus’ spin axis is far from any first-order 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 spin-orbit 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 spin-orbit coupling, while tilting Saturn to its current orientation. In this regard, the spin-axis 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 long-term spin-axis 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 long-term 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 spin-axis dynamics of Jupiter. Due to their much smaller masses, the other satellites of Jupiter do not contribute noticeably to its spin-axis 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 spin-axis 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.

Table 1

Current obliquities of the planets of the Solar System.

2 Secular dynamics of the spin axis

2.1 Equations of motion

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

where the conjugate coordinates are X (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 semi-major 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, far-away 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, close-in satellites artificially increase the oblateness and the rotational angular momentum of the planet. In the close-in 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 semi-major axis, and the mean motion of the _k_th 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 close-in and far-away regimes. Using their article, we can verify that the Galilean satellites are in the close-in 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 long-term 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 quasi-periodic 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 (Lagrange-Laplace secular system), the z series only contains the frequencies g j, while the ζ series only contains the frequencies s j (see e.g. Murray & Dermott 1999 or Laskar et al. 2012). This is not the case in more realistic situations, as recalled for instance by Kreyche et al. (2020) in the context of obliquity dynamics. In planetary systems featuring mean-motion resonances, the spin-axis 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 spin-orbit resonant angles can only be of the form σ p = ψ + ϕ p, where p is a given index in the ζ series. Resonances featuring terms of the z series only appear at third order and beyond. For the 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 quasi-integrable nature of their dynamics. Therefore, the secular spin-orbit resonances are small and isolated from each other, and only first-order resonances play a substantial role.

Figure 1 shows the location and width of every first-order resonance for the spin-axis 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 billion-year timescale.

Table 2

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 spin-axis 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 sheet1 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 constant-density core and massless envelope to 0.221 for a constant-density envelope and point-mass 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 spin-axis of Jupiter is located very close to a strong secular spin-orbit 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 model-dependent 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 long-term spin-axis dynamics of a planet with migrating satellites is described by the Hamiltonian in Eq. (1), but where α is a slowly-varying function of time. In the Earth-Moon 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 far-satellite regime (see Boué & Laskar 2006). The Galilean satellites, on the contrary, are in the close-satellite regime, and their outward migration produces an increase of α, as shown by Eq. (5). This increase can be quantified using the long-term 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’ semi-major 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. 1Location and width of every first-order secular spin-orbit 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. 2Trajectory 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. 3Typical evolution of the semi-major 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. 4Evolution 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 short-period 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 quasi-periodic 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 allfirst-order secular spin-orbit 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 non-adiabatic 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 non-adiabatic crossings of secular spin-orbit 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 spin-axis 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 turn-off 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 adiabatic2. 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 long-term 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. 5Future 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 first-order secular spin-orbit 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. 6Period of small oscillations about the resonance centre for a resonance with _ϕ_4 or _ϕ_19. Even though complete closed-form solutions exist (see Haponiak et al. 2020), the small-oscillation 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. 7Maximum 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 long-term evolution of its spin-axis under the influence of its slowly migrating satellites.

Jupiter is located today near a strong secular spin-orbit 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 spin-orbit 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 close-satellite 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 2017-40-H.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.

Table A.1

Quasi-periodic decomposition of Jupiter’s eccentricity and longitude of perihelion (variable z).

Table A.2

Quasi-periodic 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 first-order secular spin-orbit 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.1Level 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.2Same 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 spin-orbit 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 long-term spin-axis 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 two-degree-of-freedom 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.1Poincaré 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 non-adiabatic 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 non-adiabatic 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., Nichols-Fleming, 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

Table 1

Current obliquities of the planets of the Solar System.

Table 2

First twenty terms of Jupiter’s inclination and longitude of ascending node in the J2000 equator and equinox reference frame.

Table A.1

Quasi-periodic decomposition of Jupiter’s eccentricity and longitude of perihelion (variable z).

Table A.2

Quasi-periodic decomposition of Jupiter’s inclination and longitude of ascending node (variable ζ).

All Figures

In the text

	Fig. 3Typical evolution of the semi-major 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

In the text

	Fig. 7Maximum 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

In the text

	Fig. B.2Same 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.1Poincaré 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