The future large obliquity of Jupiter

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{\deg} to 37{\deg} after 5 Gyrs 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.


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. 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 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 plan- Notes. Mercury's obliquity is taken from Konopliv et al. (2020). Other values are taken from Murray & Dermott (1999) who cite the compilation made by Yoder (1995).
etary systems. The overlap of such resonances produces a large chaotic region for the spin axis of the terrestrial planets of the Solar System (see . 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;. The Moon currently protects the Earth from large chaotic variations in its obliquity 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 et al. 2014;Deitrick et al. 2018b,a;Shan & Li 2018;Millholland & Laughlin 2018Quarles 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 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. 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 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.

Equations of motion
The spin-axis dynamics of an oblate planet subject to the lowestorder term of the torque from the Sun is given for instance by  or Néron de Surgy & Laskar (1997). Far from spin-orbit resonances, and due to the weakness of the torque, the long-term evolution of the spin axis is accurately described by the secular Hamiltonian function (i.e. averaged over rotational and orbital motions). This Hamiltonian can be written where the conjugate coordinates are X (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 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 αX/(1 − e 2 ) 3/2 . The quantity α is called the precession constant. It depends on the spin rate of the planet and of its mass distribution, through the formula: where G is the gravitational constant, m is the mass of the Sun, ω is the spin rate of the planet, a is its semi-major axis, J 2 is its second zonal gravity coefficient, and λ is its normalised polar moment of inertia. We retrieve the expression given for instance by Néron de Surgy & Laskar (1997) by noting that 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: where m k , a k , and n k are the mass, the semi-major 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 close-in and faraway 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 J 2 and λ with terms of order e 2 k and η 2 k , 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 J 2 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).

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: 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 over time t, with frequencies µ k and ν k : 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 1999or 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(Quillen et al. , 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 ; Néron de Surgy & Table 2. First twenty terms of Jupiter's inclination and longitude of ascending node in the J2000 equator and equinox reference frame. Notes. Due to the secular resonance (g 1 − g 5 ) − (s 1 − s 2 ), an additional fundamental frequency γ appears in terms 18 and 20 (see Laskar 1990). (*) There is a typographical error in Laskar (1990) in the identification of the 8th term. Laskar 1997;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 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 ·yr −1 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.

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 Article number, page 3 of 12 A&A proofs: manuscript no. satgal2Jupiter 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).  (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.  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 ·yr −1 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.

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 nominal orientation. 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 ψ.

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 Gyrs. 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 all first-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. 6 shows 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 Each panel corresponds to a value of the normalised polar moment of inertia of Jupiter λ = C/(MR 2 eq ) 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 Gyrs 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 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.
λ 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 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 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 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 2 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).  If the Galilean satellites migrate faster than shown in Fig. 3, the obliquity reached in 5 Gyrs 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 (k 2 /Q) 0,1 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 Gyrs for λ sampled in our exploration interval and (k 2 /Q) 0,1 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 Gyrs 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 Gyrs in the most unfavourable case of Fig. 7.

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 λ = C/(MR 2 eq ) 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 values of λ, 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 Gyrs. For the value λ ≈ 0.252 obtained by Helled 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. 2009Lainey et al. , 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. 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. Notes. This solution has been directly obtained from Laskar (1990) as explained in the text. The phases θ (0) k are given at time J2000.