© ESO, 2011

1. Introduction

Europa and Ganymede will be the main targets of the future NASA-ESA Europa Jupiter System Mission (EJSM) (Blanc et al. 2009). One of the objectives of this mission is to determine the geophysical properties of both satellites and, in particular, to characterize water oceans beneath their icy shell. Such oceans might be crucial for the emergence of habitable worlds in the Solar System. Measuring the librational motion of the satellites can be a useful technique to prove the existence of a putative ocean beneath the ocean as suggested by Wu et al. (2001) for Europa.

The libration in longitude at the orbital period results from the variation in the orbital velocity of the satellite related to the orbital eccentricity (e.g. Comstock & Bills 2003). The librational motion in longitude considered here describes the principal axis oscillation of the body projected onto the equatorial plane of the satellite. Since the Keplerian orbital motion of the satellite is disturbed by gravitational interactions of the satellite with its neighborhood and the Sun, the elliptical elements will present a broad spectrum of oscillations. All these orbital oscillations modify the orientation of the satellite with respect to the central planet and result in libration motion. Henrard (2005a,b) has developed a Hamiltonian theory for the rigid rotation of Europa and shown that longer period librations are important, whereas direct effects of a torque due to Jupiter’s oblateness or Io’s point mass are negligible compared to the indirect effects on the orbital motion. Based on this approach, the rigid rotation of Galilean satellites has been investigated successively (Henrard & Schwanen 2004; Henrard 2005a; Noyelles 2009). A recent study of the librational motion of Enceladus, without a global ocean, further underlines the importance of long-period librations in the rotational motion of icy satellites (Rambaux et al. 2010).

The objective of this paper is to include the indirect effects of orbital perturbations on the rotational motion of Europa, Ganymede, and Callisto each of which contains a subsurface ocean. The presence of an ocean for the Galilean satellites is inferred from analysis of the geological structures at the surface of the satellites and magnetic measurements (Kivelson et al. 1999, 2000, 2002; Khurana 1998; Papalardo et al. 1999). Van Hoolst et al. (2008, 2009) and Baland & Van Hoolst (2010) have investigated the librational response of the Galilean satellites and Titan in the presence of a subsurface ocean and showed that it increases the librational amplitude. In contrast, the presence of a fluid core has a weak effect on the amplitude of librations (Baland & Van Hoolst 2010). In these studies the authors assumed that the orbit is purely Keplerian and focused on the libration in longitude at the orbital period. The amplitude of the libration essentially depends on (i) the magnitude of the external torque acting on the dynamical figure of the satellites; (ii) the thickness of the icy shell h; and (iii) torques between the different layers.

The paper is organized as follows: first, we present the outline of the librational model that follows the model developed in Van Hoolst et al. (2009). Secondly, we describe the orbital model of the Galilean satellites based on the ephemerides of Lainey et al. (2006). These numerical ephemerides ensure accuracy of at least a few tens of kilometers in the satellite position. The fourth section is dedicated to the interior structure models. The following sections are dedicated to describing the librational response for Europa, Ganyemede, and Callisto, and our conclusions are presented in the last section.

2. Librational model

Like our Moon in orbit around the Earth, the Galilean satellites are in synchronous spin-orbit resonance and present on average the same face towards their central body, i.e. Jupiter. Due to the finite eccentricity of the orbits, the orbital velocity of the moons is not constant. This variation leads to a gravitational torque from Jupiter on the moon’s dynamical figure that drives the librational response of the body. The physical libration represents the departure from the uniform rotational motion, and the libration in longitude describes the oscillatory motion in the equatorial plane of the body. The latter is described by the small libration angle γ = θ − M − θ₀, where M is the mean anomaly that evolves linearly with time, θ is the angle of rotation corresponding to the angle between the longest axis of the body and the line of the ascending node, and θ₀ is a constant representing the initial value of θ located at the ascending node of the moon’s orbit. As the obliquity of the Galilean satellites is expected to be small (Bills 2005) we focus our study on the libration in longitude.

The rotation of a body composed of a rigid icy shell, liquid ocean, and rigid solid interior can be described by the angular momentum equations for each layer l: $\frac{\mathrm{d}{{H}}_{{l}}}{\mathrm{d}\mathit{t}}\mathrm{=}{{\Gamma }}_{{l}}$(1)where H_l is the angular momentum and Γ_l is the sum of Jupiter’s gravitational torque and internal gravitational and pressure torques. As we neglect the obliquity this vectorial set of equations reduces to three equations projected along the rotation axis. In addition, the total torque on the internal ocean reduces to zero for small differential rotations of layers with respect to the synchronous rotation, and the dynamical equations can be expressed as (Van Hoolst et al. 2009): $\begin{array}{ccc}{\mathit{C}}_{\mathit{s}}\begin{array}{}\mathrm{¨}\\ {\mathit{\gamma }}_{\mathit{s}}\end{array}\mathrm{+}{\mathit{C}}_{\mathrm{o}}\begin{array}{}\mathrm{¨}\\ {\mathit{\gamma }}_{\mathrm{o}}\end{array}\mathrm{+}{\mathit{C}}_{\mathit{i}}\begin{array}{}\mathrm{¨}\\ {\mathit{\gamma }}_{\mathit{i}}\end{array}& \mathrm{=}& \frac{\mathrm{3}}{\mathrm{2}}\frac{\mathit{Gm}}{{\mathit{d}}^{\mathrm{3}}}\mathrm{\left[}\mathrm{\right(}{\mathit{B}}_{\mathit{s}}\mathrm{-}{\mathit{A}}_{\mathit{s}}\mathrm{)}\sin \mathrm{2}{\mathit{\psi }}_{\mathit{s}}\\ & & \mathrm{+}\mathrm{(}{\mathit{B}}_{\mathrm{o}}\mathrm{-}{\mathit{A}}_{\mathrm{o}}\mathrm{)}\sin \mathrm{2}{\mathit{\psi }}_{\mathrm{o}}\\ & & \mathrm{+}\mathrm{(}{\mathit{B}}_{\mathit{i}}\mathrm{-}{\mathit{A}}_{\mathit{i}}\mathrm{)}\sin \mathrm{2}{\mathit{\psi }}_{\mathit{i}}\mathrm{]}\\ {\mathit{C}}_{\mathit{i}}\begin{array}{}\mathrm{¨}\\ {\mathit{\gamma }}_{\mathit{i}}\end{array}& \mathrm{=}& \frac{\mathrm{3}}{\mathrm{2}}\frac{\mathit{Gm}}{{\mathit{d}}^{\mathrm{3}}}\mathrm{\left[}\mathrm{\right(}{\mathit{B}}_{\mathit{i}}\mathrm{-}{\mathit{A}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathrm{(}{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\sin \mathrm{2}{\mathit{\psi }}_{\mathit{i}}\\ & & \end{array}$where the first equation represents the librational motion of the whole satellite and the second equation describes the librational motion of the interior. The angle γ_l is the librational angle of layer l where l stands for s = shell, o = ocean, or i = interior. G is the gravitational constant, m the mass of Jupiter and d the relative distance between Jupiter and the moon. The angle ψ_l is the angle between the longest axis and the direction of the satellite to the Jupiter. It is equal to ψ_l = v − θ_l where v is the true longitude of the satellite from the line of the ascending node of the orbit. The right-hand-side terms depending on ψ_l describe Jupiter’s gravitational torque on the different internal layers. The constant K_int represents the magnitude of the internal gravitational and ocean pressure coupling between the shell and the interior due to the misalignment of the two layers. Details of these torques are presented in Van Hoolst et al. (2009) and references therein. The principal moments of inertia for each layer A_l < B_l < C_l are assumed to be constant. We consider that the solid layers have infinite rigidity. Hence elastic deformations and the phase lag due to anelastic response were neglected (Baland & Van Hoolst 2010). The polar moment of inertia C_l is defined by: ${\mathit{C}}_{\mathit{l}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }}{\mathrm{15}}{\mathit{\rho }}_{\mathit{l}}\left[{\mathit{r}}_{\mathrm{0}\mathit{,l}}^{\mathrm{5}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\alpha }}_{\mathit{l}}\right)\mathrm{-}{\mathit{r}}_{\mathrm{0}\mathit{,l}\mathrm{-}\mathrm{1}}^{\mathrm{5}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\alpha }}_{\mathit{l}\mathrm{-}\mathrm{1}}\right)\right]$(4)where ρ_l is the density, r_0,l is the mean radial coordinate of the outer surface of layer l from the center of the satellite, and α_l =  [(a_l + b_l)/2 − c_l]/[(a_l + b_l)/2]  is the polar flattening due to the centrifugal and static tidal potentials with a_l,b_l,c_l the radii of the three principal axes of the ellipsoids layers a_l > b_l > c_l. The equatorial moment of inertia difference B_l − A_l resulting from the static tidal potential is: ${\mathit{B}}_{\mathit{l}}\mathrm{-}{\mathit{A}}_{\mathit{l}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }}{\mathrm{15}}{\mathit{\rho }}_{\mathit{l}}\mathrm{[}{\mathit{r}}_{\mathrm{0}\mathit{,l}}^{\mathrm{5}}{\mathit{\beta }}_{\mathit{l}}\mathrm{-}{\mathit{r}}_{\mathrm{0}\mathit{,l}\mathrm{-}\mathrm{1}}^{\mathrm{5}}{\mathit{\beta }}_{\mathit{l}\mathrm{-}\mathrm{1}}\mathrm{]}$(5)where β_l is the equatorial flattening equal to (a_l − b_l)/a_l. The contribution $\mathrm{(}{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{)}$ is the moment of inertia difference for the volume of the interior with a constant density ρ_o$\mathrm{(}{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{)}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }}{\mathrm{15}}{\mathit{\rho }}_{\mathrm{o}}{\mathit{\beta }}_{\mathit{i}}{\mathit{r}}_{\mathrm{0}\mathit{,i}}^{\mathrm{5}}\mathit{,}$(6)and the corresponding term in Eq. (3) expresses the pressure torque exerted at the interior-ocean interface.

We express the torque of Jupiter on the ocean as the sum of the torques on the top and bottom parts of the ocean. Therefore, the second term inside the brackets in the right-hand side of Eq. (2) can be rewritten as $\mathrm{(}{\mathit{B}}_{\mathrm{o}}\mathrm{-}{\mathit{A}}_{\mathrm{o}}\mathrm{)}\sin \mathrm{2}{\mathit{\psi }}_{\mathrm{o}}\mathrm{=}\mathrm{(}{\mathit{B}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}\sin \mathrm{2}{\mathit{\psi }}_{\mathit{s}}\mathrm{-}\mathrm{(}{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{)}\sin \mathrm{2}{\mathit{\psi }}_{\mathit{i}}$(7)where $\mathrm{(}{\mathit{B}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }}{\mathrm{15}}{\mathit{\rho }}_{\mathrm{o}}{\mathit{\beta }}_{\mathrm{o}}{\mathit{r}}_{\mathrm{0}\mathit{,}\mathrm{o}}^{\mathrm{5}}$(8)is the equatorial moment of inertia difference of the top of the ocean.

Because the eccentricity and the physical libration are small, we retain only first-order terms in these quantities. We then have $\begin{array}{ccc}{\mathit{C}}_{\mathit{s}}\begin{array}{}\mathrm{¨}\\ {\mathit{\gamma }}_{\mathit{s}}\end{array}\mathrm{+}\mathrm{\left\{}\mathrm{3}{\mathit{n}}^{\mathrm{2}}\mathrm{\right[}\mathrm{(}{\mathit{B}}_{\mathit{s}}\mathrm{-}{\mathit{A}}_{\mathit{s}}\mathrm{)}\mathrm{+}\mathrm{(}{\mathit{B}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}\mathrm{]}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{int}}\mathrm{\}}{\mathit{\gamma }}_{\mathit{s}}\mathrm{-}\mathrm{2}{\mathit{K}}_{\mathrm{int}}{\mathit{\gamma }}_{\mathit{i}}\mathrm{=}& & \\ \mathrm{3}{\mathit{n}}^{\mathrm{2}}\mathrm{\left[}\mathrm{\right(}{\mathit{B}}_{\mathit{s}}\mathrm{-}{\mathit{A}}_{\mathit{s}}\mathrm{)}\mathrm{+}\mathrm{(}{\mathit{B}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\mathrm{(}\mathit{v}\mathrm{-}\mathit{M}\mathrm{-}{{\mathit{\theta }}_{\mathrm{0}}}_{\mathit{s}}\mathrm{)}& & \\ {\mathit{C}}_{\mathit{i}}\begin{array}{}\mathrm{¨}\\ {\mathit{\gamma }}_{\mathit{i}}\end{array}\mathrm{+}\mathrm{\left\{}\mathrm{3}{\mathit{n}}^{\mathrm{2}}\mathrm{\right[}\mathrm{(}{\mathit{B}}_{\mathit{i}}\mathrm{-}{\mathit{A}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathrm{(}{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{)}\mathrm{]}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{int}}\mathrm{\}}{\mathit{\gamma }}_{\mathit{i}}\mathrm{-}\mathrm{2}{\mathit{K}}_{\mathrm{int}}{\mathit{\gamma }}_{\mathit{s}}\mathrm{=}& & \\ \mathrm{3}{\mathit{n}}^{\mathrm{2}}\mathrm{\left[}\mathrm{\right(}{\mathit{B}}_{\mathit{i}}\mathrm{-}{\mathit{A}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathrm{(}{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\mathrm{(}\mathit{v}\mathrm{-}\mathit{M}\mathrm{-}{{\mathit{\theta }}_{\mathrm{0}}}_{\mathit{i}}\mathrm{)}& & \end{array}$where we used θ_l = M + γ_l + θ₀ in the dynamical equation. In this linearized equation in eccentricity the quantity (Gm/d)³ can be taken equal to n² by using Kepler’s third law. If the orbit of the satellite is purely Keplerian v − M can be developed as a Fourier series of the mean anomaly as a function of the eccentricity. However, due to gravitational interactions with the other satellites and the Sun, the difference v − M − θ₀ is densified by Fourier coefficients of the frequencies corresponding to the longitude and nodes of the other Galilean satellites and the Sun (see Lainey et al. 2006, and the following section). The equation of the center v − M − θ₀ is thus expressed as a Fourier series $\mathit{v}\mathrm{-}\mathit{M}\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{=}\sum _{\mathit{j}}{\mathit{H}}_{\mathit{j}}\sin \mathrm{(}{\mathit{\omega }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\alpha }}_{\mathit{j}}\mathrm{)}$(11)where H_j,ω_j and φ_j are the magnitude, frequency and phase of the orbital perturbations described in the next section.

By substituting Eq. (11) into the libration Eqs. (9) and (10), we obtain the forced solution for the libration angle of the shell γ_s and the interior γ_i as ${\mathit{\gamma }}_{\mathit{s}}\mathrm{=}\sum _{\mathit{j}}{\mathit{\gamma }}_{\mathit{s}}^{\mathit{j}}\sin \mathrm{(}{\mathit{\omega }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\alpha }}_{\mathit{j}}\mathrm{)}\mathrm{;}{\mathit{\gamma }}_{\mathit{i}}\mathrm{=}\sum _{\mathit{j}}{\mathit{\gamma }}_{\mathit{i}}^{\mathit{j}}\sin \mathrm{(}{\mathit{\omega }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\alpha }}_{\mathit{j}}\mathrm{)}$(12)with the amplitude of sine terms expressed as $\begin{array}{ccc}{\mathit{\gamma }}_{\mathit{s}}^{\mathit{j}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{C}}_{\mathit{i}}{\mathit{C}}_{\mathit{s}}}\frac{\mathrm{1}}{{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}\left[\frac{{\mathit{\phi }}_{\mathrm{1}}^{\mathit{s}}}{{\mathit{\omega }}_{\mathit{j}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\phi }}_{\mathrm{2}}^{\mathit{s}}}{{\mathit{\omega }}_{\mathit{j}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}\right]\\ {\mathit{\gamma }}_{\mathit{i}}^{\mathit{j}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{C}}_{\mathit{i}}{\mathit{C}}_{\mathit{s}}}\frac{\mathrm{1}}{{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}\left[\frac{{\mathit{\phi }}_{\mathrm{1}}^{\mathit{i}}}{{\mathit{\omega }}_{\mathit{j}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\phi }}_{\mathrm{2}}^{\mathit{i}}}{{\mathit{\omega }}_{\mathit{j}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}\right]\mathrm{·}\end{array}$ The expressions for the resonant frequencies ω₁ and ω₂ and the resonance strengths ${\mathit{\phi }}_{\mathrm{1}}^{\mathit{k}}$ and ${\mathit{\phi }}_{\mathrm{2}}^{\mathit{k}}$ are given in Appendix A. The proper modes have been described in Van Hoolst et al. (2008, 2009) and correspond to the oscillations at which the body will librate if it is slightly shifted from the dynamical equilibrium. Usually, dissipation acting on long time span will damp the proper mode oscillations to small amplitudes. Therefore we neglect these terms in the libration solution Eq. (12).

3. Orbital model

We use the numerical ephemerides of Lainey et al. (2006) to obtain an accurate representation of the orbit of the Galilean satellites. The accuracy of these ephemerides is 20 km for Europa and Ganymede, and 35 km for Callisto. The reference frame is centered on Jupiter and defined by the Jupiter equatorial plane at J1950. The orbital motion of the satellites is given in either Cartesian coordinates or classical geometrical elements. The Fourier series of the true longitude v are determined using the TRIP software (Gastineau & Laskar 2008) based on the work of Laskar (1988, 2005).Tables 1–3 show the main frequencies for each body. The minimum magnitude¹ is chosen as 1″ for Europa, Ganymede, and Callisto.

The different terms in Tables 1–3 can be identified by using the frequencies listed in the tables of Lainey et al. (2006). We follow the same notation in Lainey et al. (2006), m = 1,2,3,4 represents Io, Europa, Ganymede, and Callisto. The angular variables L_m are the linear part of the mean longitudes of the four Galilean satellites, ϖ_m the longitudes of their pericenters, Ω_m the longitudes of their nodes, ν the great inequality L₁ − 2L₂ or L₂ − 3L₃ + π, L_s the linear part of the mean longitude of the Sun (corresponding to Jupiter’s frequency around the Sun), and Ψ the argument of the Laplacian libration. The variable ρ is the De Haerdtl inequality 3L₃ − 7L₄.

The frequency analysis highlights two different timescales: (i) short-period oscillations, with periods close to the orbital period of the satellites, which stem from the orbital eccentricity of the bodies and resonant interactions among satellites; and (ii) long-period oscillations resulting from the Laplace resonance (for Europa and Ganymede), interaction with the Sun, and secular perturbations between satellites. The main forcing in Tables 1–3 is related to the equation of the center (difference between mean and true anomaly) and has an amplitude equal to two times the eccentricity. The associated frequency is the mean anomaly except for Europa for which it is given by the combination L₁ − L₂. Indeed for bodies in resonance, the eccentricity is the combination of a forced component due to the resonance and a free component dependent on the initial conditions (e.g. Ferraz-Mello 1979; Greenberg 2005). For Europa, the forced eccentricity is larger than the free eccentricity and so it is dominant in the series of v associated with the period L₁ − L₂ as shown in Table 1.

We conclude this section by noting that the frequencies of the orbital motions are not fixed on long time-scale but slowly vary in time at an observable level (Lainey et al. 2009). The variations are very slow and the orbital perturbations may be developed as Poisson series as for the Earth-Moon system. Although the introduction of Poisson terms in the orbital series will imply Poisson term and out-of-phase terms in the quasi-periodic librational development, we do not investigate such terms because their amplitudes are expected to be small.

4. Interior models

The amplitudes of the torques which drive the librations in longitude depend on the equatorial flattening of the satellite. The gravitational quadrupole moments, determined from the Galileo spacecrafts radio tracking data, are consistent with the assumption that the satellites are in tidal and rotational equilibrium (Schubert et al. 2004). The equilibrium figure and the flattening of the internal layers due to the centrifugal and tidal potentials, are calculated by using the Clairaut equation (Jeffreys 1952) extended to satellites in spin-orbit resonance (Van Hoolst et al. 2008), for a given internal density distribution. The interior models for the density are required to be consistent with the satellite mass and radius, and the mean moment of inertia I = (A + B + C)/3 determined from the quadrupole moments (see Table 5). These data by themselves can not be inverted to a unique configuration. With such limited constraints, we examine simple layered models with homogenous densities. All models are composed of an ice shell including a subsurface ocean and a rocky mantle. Oceans are predicted for all three satellites (Kivelson et al. 1999, 2000, 2002; Khurana 1998; Papalardo et al. 1999) and theoretical models shows that internal oceans are conceivable for Europa, Ganymede and for a partly differentiated Callisto with thicknesses as large as few hundred km (Schubert et al. 2004). The Europa ocean is likely situated above the rocky mantle. For Ganymede and Callisto, whose mean densities and moments of inertia suggest much thicker ice/ocean shells, the ocean may have a larger thickness and be between two ice layers. The thickness of the ice shell can be calculated theoretically by modeling internal radiogenic heating and tidal dissipation. The resulting thicknesses however are not always in good agreement with those estimated from geological evidence, especially for Europa (Billings & Kattenhorn 2005). A metallic core is considered for Europa and Ganymede (additional evidence comes from Ganymede’s intrinsic magnetic field, Kivelson et al. 1996), but not for the partially differentiated Callisto. Interior models constrained only by the mean density and the mean moment of inertia requires that choices have to be made for the values of some parameters such as the density and thickness of the ice shell and/or the subsurface ocean. In the present study, the densities of the shell and ocean are assumed to vary between 800 and 1200 kg/m³. We considered the shell thickness variations between 40 km and 120 km except for Europa where it varies between 5 km and 45 km. The simple interior models developed with homogenous density layers are consistent with previous studies (Schubert et al. 2004; Spohn & Schubert 2003; Sohl et al. 2002; Van Hoolst et al. 2008; Hussmann et al. 2006). A large number (more than 20 for each satellite) of models is constructed to obtain plausible ranges of moments of inertia, equatorial flattenings and associated torques. Table 4 summarizes the density profile for selected models used to compute the libration amplitude of Tables 6–8.

5. Europa’s librations

5.1. Description of main librations

Fig. 1 Librational response times the mean radius of the icy shell of Europa by taking into account the orbital perturbations over 4500 days a) and a zoom over 800 days b). The amplitude increase in panel a) is due to the contribution of long-periodic Fourier terms with different phases. The initial date for the plot is January 1st 2026 corresponding to the expected date of arrival of EJSM at Jupiter. The interior model comes from the Table 4.

Fig. 2 Variation, with shell thickness (h), of the amplitude multiplied by the mean radius for the 4 main librations, 3.52 days a), 7.05 days b), 482 days c), and 485 days d).

The librational response driven by the external torque exerted by Jupiter on Europa is obtained from Eqs. (12)–(13) for the orbital perturbations (Table 1). We first calculate the librations for the interior structure of Europa composed of a solid interior (core plus mantle), a water ocean, and a solid icy shell listed in Table 4. The temporal behavior of the libration from January 1st 2026, date of expected arrival EJSM, is shown in Figs. 1a,b over 4500 and 800 days, respectively. The main libration is seen to have a short period but its amplitude appears to be strongly modulated on long periods of about 480 days. The short-period libration has an amplitude of about 200″ corresponding to a surface displacement at the equator of 1.4 km and the long period of 480 days has almost the same amplitude. Here, the amplitude is given as a surface displacement at the equator by multiplying the libration angle by the mean radius of the satellite.

A deeper insight into the motion is obtained by analyzing directly the librational motion in frequency domain. Table 6 lists the physical librations in order of decreasing amplitude. The main libration is a 3.52 days oscillation, i.e. the mean mean longitude of Europa perturbed by Io (L₁ − L₂). The forcing for this term is mainly related to the eccentricity of the orbit (see Sect. 3). Next, librations associated with long periods of 482.06 days, 485.60 days, and 462.51 days, which are related to the precession of the nodes of the orbits of satellites (ν + ϖ_k where k stands for each Galilean satellites), have amplitudes around 70 and 10% smaller than the main term. The influence of the resonant interactions and the Sun are the 5th and 6th most important librations. We also note the presence of the L₂ − L₃ oscillation, with a period approximatively twice the orbital period, although its magnitude is relatively small.

The comparison of the magnitude of the forcing terms in Table 1 and the amplitude of librations in Table 6 shows that for the low frequency (long period) terms the forcing magnitude H_j and the libration amplitude ${\mathit{\gamma }}_{\mathit{s}}^{\mathit{j}}$ are similar whereas the amplitudes of the high frequency (short period) terms are strongly diminished with respect to the forcing. The separation between high and low frequencies is with respect to the proper frequencies and is due to the presence of the forcing frequency in the denominator (Eq. (13)), which decreases the libration amplitude for short periods. This dynamical behavior corresponds to the dynamics of a coupled forced oscillator where the restoring force is Jupiter’s gravitational torque, the coupling force is K_int(γ_s − γ_i), and the moons inertia is ${\mathit{C}}_{\mathit{s}}\begin{array}{}{}_{\mathrm{¨}}\\ {\mathit{\gamma }}_{\mathit{s}}\end{array}$ approximated as ${\mathit{C}}_{\mathit{s}}{\mathit{\omega }}_{\mathit{j}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathit{s}}^{\mathit{j}}$ for quasi-periodic librational response (Eqs. (9) and (10)). For small frequencies, the gravitational torque is dominant and the orientation of the moon follows the direction of the force by keeping the same face toward Jupiter, i.e. the amplitude of the libration (γ_l for l = s,i) equals to the magnitude of the longitude perturbation (v) and hence ψ_s = 0 (ψ_s = v − θ_s defined in Sect. 2). We note that in this case the coupling torque vanished since the three angles (γ_s,γ_i,v) are in phase. On contrast, for high frequency oscillations the inertia of the moon has a strong impact and the moon does not orient its long axis toward Jupiter anymore. This behavior is similar to the case of Enceladus (Rambaux et al. 2010).

5.2. Influence of geophysical parameters

Figure 2a shows that the main geophysical parameter that drives the libration at the orbital frequency is the thickness of the icy shell h (Van Hoolst et al. 2008). Here a similar hyperbolic behavior dependence of libration amplitude on shell thickness is obtained for the libration at 7.05 days (Fig. 2b). On contrast, the amplitudes of the librations at 482 days and 485 days depend only very weakly on shell thickness (see Figs. 2c,d), in accordance with the observation above that the libration amplitude at low frequency is almost equal to the restoring force, i.e. Jupiter’s gravitational interaction, magnitude. By introducing the condition that the forcing frequency ω_j ≪ ω₁,ω₂ in the librational solution (Eq. (12)), the libration amplitude of the shell can be approximated correct up to the third order in ω_j as ${\mathit{\gamma }}_{\mathit{s}}^{\mathit{j}}\mathrm{~}{\mathit{H}}_{\mathit{j}}\mathrm{+}{\mathit{H}}_{\mathit{j}}\frac{\mathrm{\left(}\mathrm{2}{\mathit{K}}_{\mathrm{int}}\mathrm{\right(}{\mathit{C}}_{\mathit{s}}\mathrm{+}{\mathit{C}}_{\mathit{i}}\mathrm{)}\mathrm{+}{\mathit{C}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{)}}{\mathrm{\left(}\mathrm{2}{\mathit{K}}_{\mathrm{int}}\mathrm{\right(}{\mathit{K}}_{\mathit{s}}\mathrm{+}{\mathit{K}}_{\mathit{i}}\mathrm{)}\mathrm{+}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{)}}{\mathit{\omega }}_{\mathit{j}}^{\mathrm{2}}\mathrm{+}\mathit{O}\mathrm{(}{{\mathit{\omega }}_{\mathit{j}}^{\mathrm{4}}}^{\mathrm{)}}\mathit{.}$(15)where the K_j (j = i,s) are defined in the Appendix A and depend on interior parameters. At first order in ω_j, the amplitude ${\mathit{\gamma }}_{\mathit{s}}^{\mathit{j}}$ is equal to the forcing magnitude H_j. The small dependence on the interior is given by the coefficient of the term in the square of the forcing frequency. The dependence on h is very weak, with differences in amplitude of the order of a meter (Figs. 2c, d) that is too small to be observed. Therefore, it will be not possible to obtain information on the interior by measuring these long-period librations.

The amplitude of high frequency librations (Figs. 2a,b) depends not only on the magnitude of the coupling (H_j in Eq. (11)) but also on the proximity of the forcing frequency (ω_j) to the eigenvalues of the proper modes ω₁ and ω₂ as shown in Eq. (13). Whereas ω₂ (between 2π/52 and 2π/60 days^-1 for our models) is far from the forcing frequencies, ω₁ can be close to the largest frequencies for small shell thickness (Fig. 3). For small thickness of the icy shell the proper period (2π/ω₁) tends toward a few days (Fig. 3), i.e. very close to the forcing periods of 7.05 days and 3.52 days, and leads to the resonance shape curves shown in Figs. 2a,b. For the 7.05 days libration the increase is significant because the resonance is crossed for h ~ 10 km. In the linear model developed in Sect. 2 the amplitude tends toward infinity that is contrary to the hypothesis of small libration amplitude. In order to compare with a more realistic libration model at the resonance, we numerically integrated the non-linear Eq. (2). The resulting libration is shown with red triangles. Outside the resonance both results coincide and the difference appears only at the resonance as expected. The two lines of resonance (blue and violet) have been computed by assuming that only B_s − A_s and C_s depend on h and keeping all other geophysical parameters constants. In this case, the amplitude can be expressed as ${\mathit{\gamma }}_{\mathit{s}}^{{\mathit{L}}_{\mathrm{2}}\mathrm{-}{\mathit{L}}_{\mathrm{3}}}\mathrm{=}\frac{\mathrm{35}\hspace{0.17em}\mathrm{009.990}\mathrm{+}\mathrm{0.0145}{\mathit{\rho }}_{\mathit{s}}\mathit{h}}{\mathrm{65}\hspace{0.17em}\mathrm{510.895}\mathrm{-}\mathrm{7.793}{\mathit{\rho }}_{\mathit{s}}\mathit{h}}$(16)where ρ_s is the density of the shell expressed in kg/m³ and the icy thickness h is expressed in km. The blue dashed curve is plotted for a density equal to 800 kg/m³ while the violet line is plotted for a density of 1000 kg/m³. These curves highlight the importance of the shell density on the position of the resonance (in the frequency space). Therefore, the amplitude of the libration is sensitive to the icy thickness h and to the density of the ice shell because the last parameter controls the position of the resonance.

The observational determination of the two librations at 3.52 and 7.05 days might separate the contribution of the icy thickness h and icy density ρ_s. These two geophysical parameters are crucial for geophysical, geological, and astrobiological studies of Europa. The presence of a resonance in the librational motion of Europa may increase significantly the oscillatory motion of the satellite and thus tidal heating and surface faulting of the moon. However, to determine the amplitude of librations in resonance, the rotational model needs to take into account dissipative torques that will bound the highest amplitude and modify slightly the location of the resonance. In the available observations no signature of such dramatic events have been found but it cannot be ruled out that today or during its history Europa has reached such resonance. Hussmann & Spohn (2004) showed that the icy thickness of Europa may vary during time due to the variation of the eccentricity and may have reached values as small as 3 km.

Fig. 3 Variation of the period of the eigenvalue of the first proper mode for Europa.

6. Ganymede’s librations

Fig. 4 Librational response times the mean radius of Ganymede from orbital perturbations. The main libration is an oscillation at the longitude of the Sun period. The initial date for the plot is January 1st 2026.

In contrast to Europa, Ganymede’s largest libration is due to solar perturbations L_s and is not of orbital period (Fig. 4 and Table 7) because Ganymede is further away from Jupiter. As can be seen in Table 7, the second to fourth largest libration terms are due to the precession of the nodes ν + ϖ_m with m = 2,3,4. The libration in longitude at the orbital period 7.15 days is only the 6th largest term and is almost 3.8 times smaller than the main term. We recall that this mode is essentially due to the eccentricity of the orbit (see Sect. 3). In the case of Ganymede the eccentricity is around 0.0013, about a factor 7 smaller than the eccentricity of Europa (0.0094). Since the long-period librations are dominant, a complete study of the librational motion of Ganymede including the long-period librations is clearly required for any rotational data analysis of Ganymede.

Fig. 5 Variation of the period of the first eigenvalue of the proper mode for Ganymede. For an icy shell of thickness equal to 60 kilometers, the proper mode is of the order of 50 days that is very close to the forcing mode of L3 − 2L4 + ϖ4 and leads to a large resonance effects.

As explained for Europa, the amplitude of the long-period librations almost does not depend on the existence of a subsurface ocean and can therefore not be used to study the interior structure of Ganymede. The dependence of the libration at orbital period on the thickness h of the icy shell for Ganymede is similar to the case of Europa (see also Baland & Van Hoolst 2010). Of particular interest here is that the libration with a period of 50 days associated with the forced libration L₃ − 2L₄ + ϖ₄ related to the interaction with Callisto may be resonant with the first proper mode, which has a period tending toward 50 days for an icy shell thickness of h ~ 60 kilometers (Fig. 5). Such a resonance might lead to a large libration amplitude at that period if the forced and proper frequencies coincide.

7. Callisto’s librations

The most distant satellite of the Galilean moons is Callisto and contrary to Io, Europa, and Ganymede, it is not involved in the Laplace resonance. Therefore, its libration series (Table 8) is dominated by the Sun perturbations (at 4332 days) and there are no components containing the great inequality. The behavior of the librational response of Callisto is shown in Fig. 6 with the geophysical parameters listed in Table 4. The second main libration is the long period of  − ϖ₃ + ϖ₄ − Ω₃ + Ω₀. The third oscillation is the libration in longitude at the orbital period related to the eccentricity. Its amplitude is almost 10 times smaller than the main term. The long-period librations are even more dominant than for Ganymede and Europa due to Callisto’s larger distance from Jupiter.

The long period behavior is similar to that for Europa and Ganymede in the sense that the amplitudes of long-period librations are close to the magnitudes of the forcing terms. However, in the short period regime a resonance with a proper mode does not seem to be possible. The proper periods are in the range of 250–550 days for P₁ and 775–1175 days for P₂, which are far from any forcing frequency.

Fig. 6 Librational response multiplied by the mean radius of Callisto from orbital perturbations. The main libration is a long period at 4332 days due to the Sun orbital perturbation. The initial date for the plot is January 1st 2026.

8. Discussion and conclusion

In this paper, we have studied the effect of orbital perturbations on the librational response of Europa, Ganymede, and Callisto. The orbital perturbations are induced by mutual gravitational interactions among satellites, especially the resonant interactions related to the Laplace resonance for Europa and Ganymede, and the solar gravitational interaction. These perturbations, which have been determined from the most recent numerical ephemerides of the Galilean moons (Lainey et al. 2006), introduce a wide range of frequencies in the librations. The resulting librations can be classified in two groups depending on the value of the frequencies with respect to the proper frequencies. For small frequencies, the restoring torque is dominant and the moons librate by keeping the same face toward Jupiter in phase with the perturbations. The resulting amplitude is almost equal to the magnitude of the orbital perturbations (at first order in frequencies). For the second group, the frequencies are long and around the orbital frequency. In this case, the inertia of the moon plays an important role and the amplitude is reduced with respect to the magnitude of the perturbations. Nevertheless, they can be strongly increased by a resonance between a proper frequency and forced frequencies. The limit between the short and long periods is essentially controlled by the highest proper frequency (ω₁).

The short-period libration amplitudes may increase in the presence of a subsurface oceans as shown in Table 9. The observed libration amplitude depends on the strength of the external torque and are modified by the internal structure. They would provide information on the exact coupling mechanisms between the shell and the interior and can be used to better model satellite’s interior. In contrast to previous studies, which only considered libration at orbital period, we used the whole spectrum of the Galilean satellites to study the dependence of librations on the internal structure. In particular for Europa, we identified a second short-period libration at twice the orbital period, whose amplitude, although an order of magnitude smaller than the amplitude of the libration at 3.52 days, could also be used to infer information on the thickness and density of the icy shell relevant for geophysical, geological, and astrobiological studies (Eq. (16)).

The long-period librations are related to the orbital motion of Jupiter around the Sun with a period of 4332 days and to the motion of the orbital nodes and great inequality ν, leading to periods of  ~480 days. In contrast to the high frequency librations, the amplitude of the small frequency librations does almost not depend on the distribution of mass inside the satellites. As a consequence, this result is robust and does not depend on the coupling mechanisms. In addition, it is important to underline that an adequate description of the rotational motion of the satellites requires to take into account the long-period librations because their amplitudes are not negligible and even dominant for Ganymede and Callisto. If such librations are neglected in the reduction process, the fitted rotational speed will be increased (or decreased depending on the phase) leading to a misstatement of non-synchronous rotation of satellites.

The oscillatory motions of the icy Galilean moons present a wide spectrum of frequencies due to orbital perturbations. Some of the librations may present large amplitudes due to a resonance with one of the proper modes. The exact position of the resonance in frequency space, i.e. the value of the proper frequencies, is sensitive to the values of the geophysical parameters introduced in our model. By exploring a large range of possible values of structure parameters, we have shown that such resonances might exist for Europa and Ganymede for realistic ranges of parameters. Although no observational evidence is available for a resonant libration for Europa or Ganymede, it cannot be ruled out that these bodies have crossed a libration resonance during their histories, leading to potential surface faulting and increased tidal heating.

In the present libration model, we have taken into account the gravitational torque exerted by Jupiter and both pressure and internal gravitational couplings, but neglected (visco-)elastic effects. The largest elastic effect can be expected to be due to the tidal forcing of Jupiter and, by analogy to solid rotation, it could contribute up to 10% to the libration amplitude at orbital frequency (Baland & Van Hoolst 2010). The effect of inelastic tides on libration is most likely very small, as shown e.g. for Enceladus (Rambaux et al. 2010) without liquid ocean. Viscous relaxation of the ice shell is also thought not to be efficient in changing the libration because it is too slow for short periodic librations and not effective for the long-period librations for which the satellites always keep the same face toward Jupiter (see Sect. 5). In addition, friction at the shell/ocean interface and flow in the ocean (Noir et al. 2009; Tyler 2008) could modify the amplitude of the librations.

The librational motion of the Galilean satellites may be inferred by spacecraft measurements by tracking surface landmarks, by fitting the shape of the satellites from surface images (Tiscareno et al. 2010), and/or by altimeter measurements. The accuracy obtained in the pole location of Titan from the different flybys by the radar Cassini instrument is below 1 km. However, Wu et al. (2001) predicted for Europa an accuracy of 0.02 km with 15 days of tracking from 2 Earth stations of an Europa orbiter. At such level of precision (2.8″) all librations presented in Table 6 may be detected. Let us mention that there is not yet enough powerful techniques to reach the Galilean moons directly from the Earth as one for Mercury by radar interferometry (Margot et al. 2009) or the Moon by Lunar Laser Ranging (Williams et al. 2001) but we could expect that some techniques will be improved in the near future.

Acknowledgments

N.R. thanks V. Lainey for fruitful discussions on the orbital ephemerides of Galilean satellites.

References

Appendix A: Expression of the librational angles

We introduce the following notation: ${\mathit{K}}_{\mathit{s}}\mathrm{=}\mathrm{(}{\mathit{B}}_{\mathit{s}}\mathrm{-}{\mathit{A}}_{\mathit{s}}\mathrm{)}\mathrm{+}\mathrm{(}{\mathit{B}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}$ and ${\mathit{K}}_{\mathit{i}}\mathrm{=}\mathrm{(}{\mathit{B}}_{\mathit{i}}\mathrm{-}{\mathit{A}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathrm{(}{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{)}$. ${\mathit{\phi }}_{\mathrm{1}}^{\mathit{s}}\mathrm{=}{\mathit{H}}_{\mathit{j}}\mathrm{(}\mathrm{-}{\mathit{K}}_{\mathit{s}}{\mathit{C}}_{\mathit{i}}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathrm{int}}\mathrm{+}{\mathit{K}}_{\mathit{s}}\mathit{Ki}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{int}}{\mathit{K}}_{\mathit{i}}\mathrm{)}$(A.1)${\mathit{\phi }}_{\mathrm{2}}^{\mathit{s}}\mathrm{=}\mathrm{-}{\mathit{H}}_{\mathit{j}}\mathrm{(}\mathrm{-}{\mathit{K}}_{\mathit{s}}{\mathit{C}}_{\mathit{i}}{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathrm{int}}\mathrm{+}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{int}}{\mathit{K}}_{\mathit{i}}\mathrm{)}$(A.2)${\mathit{\phi }}_{\mathrm{1}}^{\mathit{i}}\mathrm{=}{\mathit{H}}_{\mathit{j}}\mathrm{(}\mathrm{-}{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{s}}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{int}}{\mathit{K}}_{\mathit{i}}\mathrm{+}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathrm{int}}\mathrm{)}$(A.3)${\mathit{\phi }}_{\mathrm{2}}^{\mathit{i}}\mathrm{=}\mathrm{-}{\mathit{H}}_{\mathit{j}}\mathrm{(}\mathrm{-}{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{s}}{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{int}}{\mathit{K}}_{\mathit{i}}\mathrm{+}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathit{s}}{\mathit{K}}_{\mathrm{int}}\mathrm{)}$(A.4)$\begin{array}{ccc}\mathrm{\Delta }& \mathrm{=}& \mathrm{4}\mathrm{(}{\mathit{C}}_{\mathit{i}}\mathrm{+}{\mathit{C}}_{\mathit{s}}{\mathrm{)}}^{\mathrm{2}}{\mathit{K}}_{\mathrm{int}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{C}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{-}{\mathit{C}}_{\mathit{i}}{\mathit{K}}_{\mathit{s}}\mathrm{\left)}\mathrm{\right(}{\mathit{C}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{-}{\mathit{C}}_{\mathit{i}}{\mathit{K}}_{\mathit{s}}\\ & & \mathrm{+}\mathrm{4}{\mathit{K}}_{\mathrm{int}}\mathrm{(}{\mathit{C}}_{\mathit{s}}\mathrm{-}{\mathit{C}}_{\mathit{i}}\mathrm{)}\mathrm{)}\end{array}$(A.5)${\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}{\mathit{C}}_{\mathit{i}}{\mathit{K}}_{\mathrm{int}}\mathrm{+}{\mathit{C}}_{\mathit{i}}{\mathit{K}}_{\mathit{s}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathit{s}}{\mathit{K}}_{\mathrm{int}}\mathrm{+}{\mathit{C}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{+}\sqrt{\mathrm{\Delta }}}{\mathrm{2}{\mathit{C}}_{\mathit{s}}{\mathit{C}}_{\mathit{i}}}$(A.6)${\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}{\mathit{C}}_{\mathit{i}}{\mathit{K}}_{\mathrm{int}}\mathrm{+}{\mathit{C}}_{\mathit{i}}{\mathit{K}}_{\mathit{s}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathit{s}}{\mathit{K}}_{\mathrm{int}}\mathrm{+}{\mathit{C}}_{\mathit{s}}{\mathit{K}}_{\mathit{i}}\mathrm{-}\sqrt{\mathrm{\Delta }}}{\mathrm{2}{\mathit{C}}_{\mathit{s}}{\mathit{C}}_{\mathit{i}}}\mathrm{·}$(A.7)

All Tables

All Figures

Fig. 1 Librational response times the mean radius of the icy shell of Europa by taking into account the orbital perturbations over 4500 days a) and a zoom over 800 days b). The amplitude increase in panel a) is due to the contribution of long-periodic Fourier terms with different phases. The initial date for the plot is January 1st 2026 corresponding to the expected date of arrival of EJSM at Jupiter. The interior model comes from the Table 4.

In the text

	Fig. 2 Variation, with shell thickness (h), of the amplitude multiplied by the mean radius for the 4 main librations, 3.52 days a), 7.05 days b), 482 days c), and 485 days d).
In the text

	Fig. 3 Variation of the period of the eigenvalue of the first proper mode for Europa.
In the text

	Fig. 4 Librational response times the mean radius of Ganymede from orbital perturbations. The main libration is an oscillation at the longitude of the Sun period. The initial date for the plot is January 1st 2026.
In the text

	Fig. 5 Variation of the period of the first eigenvalue of the proper mode for Ganymede. For an icy shell of thickness equal to 60 kilometers, the proper mode is of the order of 50 days that is very close to the forcing mode of L3 − 2L4 + ϖ4 and leads to a large resonance effects.
In the text

	Fig. 6 Librational response multiplied by the mean radius of Callisto from orbital perturbations. The main libration is a long period at 4332 days due to the Sun orbital perturbation. The initial date for the plot is January 1st 2026.
In the text