© ESO 2018

1. Introduction

The Jovian system holds a special place in planetary science at least since the time of Galileo, and Jupiter itself has been tracked by Babylonian and Chinese astronomers. Contemporary dynamical theories of the system have shown that at least for the Galilean moons, the major perturbations of the orbits arise from the oblateness of Jupiter and the mutual gravitational interactions among the satellites, thus making them part of a miniature solar system. Furthermore, in the last years of space exploration, several missions (Galileo, Juno, and Juice) have been sent or will be launched to refine our understanding of the system’s physical characteristics, and how they relate to the dynamics of the system.

The key dynamical aspect at the core of these studies is the so-called Laplace resonance, a three-body mean motion resonance coupling the dynamics of the moons Io, Europa, and Ganymede. This resonance, whose main geometrical effect is the prevention of a triple conjunction of the moons (see Fig. 1), is the product of two different 2:1 quasi-resonant interactions, which can be described in terms of the observed mean motions (n₁ for Io, n₂ for Europa, and n₃ for Ganymede) as (Brown 1977)$$\begin{array}{c}\hfill n_1-2n_2=0.{7395507361}^{\circ }/\mathrm{day},\\ \hfill n_2-2n_3=0.{7395507301}^{\circ }/\mathrm{day}.\end{array}$$(1)

Fig. 1. Jovicentric snapshots of the Jovian system, with T denoting Io’s mean orbital period. The moons are color-coded (Io in red, Europa in green, and Ganymede in blue).

In addition to these two, a weaker resonant interaction exists between Ganymede and Callisto, quantifiable as 3n₃ − 7n₄ = −0.049084320°/day (again, see Brown 1977).

The relations in Eq. (1) were first shown by Laplace (1805) to imply the existence of the librating angles$$\begin{array}{cc}& {\lambda }_1-2{\lambda }_2+{\varpi }_2\approx {180}^{\circ },\hfill \\ \hfill & {\lambda }_2-2{\lambda }_3+{\varpi }_2\approx 0^{\circ },\hfill \end{array}$$(2)

where λ_i denotes the mean longitudes of the moons and ϖ_i the longitudes of the perijoves. These equations combine in the Laplace argument φ_L ≡ λ₁ − 3λ₂ + 2λ₃ ≈ 180°. Furthermore, because ϖ₁ and ϖ₂ drift with similar rates (Yoder & Peale 1981), the system is affected by the additional librating relation λ₁ − 2λ₂ + ϖ₁ ≈ 0°.

Laplace’s theory was the first that could be considered relatively complete, but it was not specifically aimed to generate ephemerides for the system. Another model, similar in scope and important for our work, was developed some decades later by de Sitter (1931). It uses intermediary orbits, which as shown in Broer & Hanßmann (2016) are periodic solutions of the differential equations obtained by retaining exclusively resonant terms. Very recently, we have shown (Celletti et al. 2018a) that the stable configuration of this system differs from the actual configuration because the latter possesses a rotating angular combination that is fixed in the former.

At about the same time as de Sitter, Sampson (1921) developed his own theory on the motion of the Galilean moons, tailored to easily generate ephemerides tables. This was later corrected and expanded by Lieske (1977, 1997), whose E5 ephemerides are the reference today for precisely tracking the orbits of the Galilean moons, as reported, for example, by Musotto et al. (2002), Lainey et al. (2004), and Kosmodamianskii (2009).

In parallel to these “astronomical” works, the search for explicit approximate solutions of the system saw the introduction of additional theories, even Hamiltonian ones, starting from the 1970s. First in chronological order, Marsden (1966) applied von Zeipel’s method (Morrison 1966) in his PhD Thesis to average out the short-period terms, and he then solved the resulting differential equations for the long-period effects by successive substitutions. Later on, Sagnier (1975) and Ferraz-Mello (1975, 1979) worked with complex variables and Taylor expansions to find their solutions, particularly assuming fixed the resonant frequencies. Different sets of variables were used by Brown (1977), who exploited a Lie transform normalization (Kamel 1970; Giorgilli et al. 2017) to remove short-period terms in his non-Hamiltonian formulation. The latter was instead employed by Henrard (1984), who reduced his system to 4 d.o.f. (hereafter, short for degrees of freedom) and then explicitly computed amplitude and frequency of the Laplace resonance using a similar technique.

We here continue several of these works in the sense that our overall objective is the construction of approximate models for the study of the Laplace resonance over short timescales and not necessarily precise ephemerides of the Galilean moons (as done, e.g., in Celletti et al. 2018a). Similarly to Celletti et al. (2018b) and Lari (2018), the goal is extrapolating information on the dynamics of this particular resonance and provide a baseline reference for future analyses of its mechanisms.

By eschewing the search for precise ephemerides, we are able to obtain several advantages. In particular, we present models that use the fewest parameters possible to describe the Laplace resonance with a sufficient degree of precision, that is, the values for the main amplitude and frequency of the libration mostly agree with those reported in the literature (Lieske 1977; Musotto et al. 2002). Thanks to their simplicity, they are easy to reproduce and can be reused in conjunction with different objectives and approaches. The caveat, and this is an important point, is that they show a strong sensitivity to the initial conditions, especially with regard to the amplitude of φ_L. Still, we show that by exploiting the inherent dynamics of the system and normal form techniques, it is possible to construct a procedure to recover what we need, as long as a tradeoff on the initial conditions is taken into account.

The first step is to introduce a benchmark model for the amplitude and frequency of the resonant argument. While precise ephemerides are not a concern, we have found Lainey’s Cartesian formulation (Lainey 2002) to be particularly suitable for our purposes. Thus, after reviewing his model and trimming it to the fundamental dynamical contributions, we extend his work by evaluating the evolution of the Laplace argument over the period of 100 yr, comparing it to NASA’s ephemerides, and extracting the fundamental amplitude and frequency that we consider in the rest of the paper.

The main model of this work is Hamiltonian, which suits our analytical objectives because of the powerful results associated with this formulation. The basic, planar version of the model presented here rests on a chain of suitable series expansions and the retaining of terms corresponding to the 2:1 quasi resonant interactions described at the beginning (plus the secular effect of Jupiter’s J₂ harmonic). This degree of development is sufficient to analyze the validity of the model, and through numerical comparison with the Cartesian formulation, we show that under the same initial conditions, the Hamiltonian adequately approximates the Galilean moons’ dynamics in terms of eccentricities and semi-major axes. Amplitude and frequency of the Laplace argument, as well as the diagnostics n₁ − 2n₂ and n₂ − 2n₃, are instead more difficult to recover. In particular, under a suitable set of canonical coordinates, we apply the normalization procedure previously mentioned to focus on a certain equilibrium of the system, obtain a 1 d.o.f. normal form Hamiltonian, and exploit its characteristics to approximate the quantities we search for.

The paper is organized as follows. In Sect. 2 we present the Cartesian formulation employed as baseline for the history of the Laplace argument, and we compare it numerically with a set of ephemerides in order to justify the dynamical effects included in the model. The conclusions of this section are then exploited in Sect. 3 to construct the alternative, planar Hamiltonian formulation for the Laplace resonance. The validity of this is subsequently confirmed in Sect. 4, where a numerical comparison is drawn with the results associated with the Cartesian model and where we show how appropriate values for the desired quantities can be recovered. Finally, in Sect. 5 we summarize the main results of this work and present parallel and future directions of research.

2. Benchmark model for the Laplace argument

2.1. Introduction and reference frames

In this section we introduce a Cartesian model that precisely captures the evolution of the Laplace argument (over the period of 100 yr). Through straightforward numerical integrations, we compare the results of this model with the equivalent results associated with a given set of ephemerides, and we exploit these comparisons to single out the dynamical effects that are used in the successive Hamiltonian formulation. Key variables for these comparisons are not only amplitude and period of the Laplace argument, but because of the composite nature of the resonance, also the 2:1 quasi-resonant observed mean motion combinations n₁ − 2n₂ and n₂ − 2n₃.

The equations of motion are defined and integrated in a Jovicentric fixed frame. Here we used the “J2000” frame implemented in the NASA Spice toolkit ¹, from which we also extracted the set of ephemerides giving the history of the Laplace argument. As indicated in the Spice documentation, the previous frame is assumed to be individualized by Earth’s mean equatorial plane at the J2000 epoch (positive x-axis parallel to the vernal equinox direction) and the corresponding normal (with Earth’s spin axis indicating the positive z direction). In reality, it corresponds to the so-called International Celestial Reference Frame, but the two are separated by a rotation of just 0.1 arcsecond. For more detail, we refer to the Spice documentation².

The angular librating relations characterizing the Jovian system depend on the osculating orbital elements, which in turn are defined through Jovicentric equatorial coordinates. Thus, before visualization, the state vectors are transformed from the J2000 frame into a frame defined by Jupiter’s mean equatorial plane and the corresponding spin axis. The latter defines the z-axis of this frame, while the line of nodes obtained by intersecting the previous equatorial plane with the J2000 plane acts as the associated x-direction. We note that this new frame is non-inertial, since Jupiter’s spin axis is precessing (we do not take nutation into account). However, since this effect is very small, we considered it only for the representation of the ephemerides, while for the integrated flow, we froze the rotation at the initial integration epoch. A complete graphical visualization of the two frames (J2000 and mean equatorial) is given in Fig. 2.

Fig. 2. Precession and rotation of Jupiter’s mean equatorial frame (P0, x′, y′, z′) with respect to the J2000 frame (P0, x, y, z), represented through the Euler angles (ψ, I, χ).

This section is organized as follows. In Sect. 2.2 we introduce our notations and derive the general form of the equations of motion. Following this, in Sect. 2.3 we discuss the precession of Jupiter’s mean north pole, which, as shown in Sect. 2.4, plays a role in how the associated oblateness potential is shaped. In Sect. 2.5 we provide some technical details on the numerical simulations we performed, whose results are then described in Sect. 2.6.

To conclude, we observe that much of the content of this section, along with the notation we adopted, follows Lainey et al. (2004) and Lainey (2002). The focus of these works is the construction of an accurate set of ephemerides for the Galilean moons. Here instead we determine reference values and characteristics for the successive Hamiltonian formulation.

2.2. Equations of motion

The first step for constructing the Cartesian formulation of the Galilean moons’ dynamics is to consider the system they form with Jupiter as isolated. These bodies are denoted with the symbol P_i throughout, where the subscript i assumes higher values the greater the distance from Jupiter. Thus, Jupiter itself is denoted with P₀, Io with P₁, Europa with P₂, Ganymede with P₃, and Callisto with P₄.

Let r_i and r_ij denote the relative vectors $\overline{P_0{P}_i}$ and $\overline{P_i{P}_j}$, respectively, with r_i and r_ij indicating the associated Euclidean norms. Furthermore, let O denote the barycenter of the system. By considering the vectorial equality ${\mathit{r}}_i=\overline{O{P}_i}-\overline{O{P}_0}$, we can derive the equations describing the dynamics in a Jovicentric frame (P₀, x, y, z). These read$$\begin{array}{c}\hfill {\ddot{\mathit{r}}}_i=\frac{{\mathit{F}}_i}{m_i}-\frac{{\mathit{F}}_0}{m_0},i=1,\cdots ,4,\end{array}$$(3)

where m_i and m₀ denote the mass of the moon and the mass of Jupiter, respectively, while F_i and F₀ stand for the whole external forces acting on the two masses. With a notation similar to the one adopted for the relative vectors, we indicate with F_ij the force exerted by Pj over P_i. Thus, for example, if we introduce the simplified gravitational potential $U_{\mathit{ij}}=U_{\overline{i}\overline{j}}=\frac{1}{r_{\mathit{ij}}}$, the gravitational attraction exerted by P_j over P_i becomes F_ij = Gm_im_j∇_iU_ij.

In the notation for the gravitational potential we have placed a bar above the indices. This means that the body corresponding to that index is modeled as a point mass. In a similar manner, we model an eventual contribution due to an oblate shape by a hat above the corresponding subscript. For instance, the total force exerted by Jupiter over the ith moon can be written as$$\begin{array}{c}\hfill {\mathit{F}}_{i0}={\mathit{F}}_{\overline{i}\overline{0}}+{\mathit{F}}_{\overline{i}\widehat{0}},\end{array}$$(4)

where F_ī0̄ = Gm_im₀∇_iU_ī0̄ and F_ī0̂ = Gm_im₀∇_iU_ī0̂. This last term represents the action of the P₀ triaxiality over P_i, and it can be expressed using equatorial spherical coordinates (r_i, φ_i, λ _i), where the last two variables denote the latitude and longitude of P_i (the vector length is invariant under rotations), respectively. If we denote by R the equatorial radius of P₀, the potential U_ī0̂ describing the previous action can be written as$$\begin{array}{c}\hfill U_{\overline{i}\widehat{0}}=U_{\overline{i}\widehat{0}}^{(1)}+U_{\overline{i}\widehat{0}}^{(2)},\end{array}$$(5)

where$$U_{\overline{i}\widehat{0}}^{(1)}={\displaystyle \sum _{n = 2}^{\infty }\frac{R^n}{r_i^{n + 1}}}{J}_n{\tilde{P}}_n\left(\sin {\varphi }_i\right)$$(6)

and$$\begin{array}{c}\hfill U_{\overline{i}\widehat{0}}^{(2)}={\displaystyle \sum _{n=2}^{\infty }\frac{R^n}{r_i^{n+1}}{\displaystyle \sum _{m=1}^n{\tilde{P}}_n^m(sin{\phi }_i)}}[C_{\mathit{nm}}cosm{\lambda }_i+S_{\mathit{nm}}sinm{\lambda }_i].\end{array}$$(7)

The quantities J_n, C_nm, and S_nm are all constants depending on the particular primary, while ${\tilde{P}}_n^m$ denotes the associated Legendre polynomials (with P̃_n equal to ${\tilde{P}}_n^0$). For an in-depth treatment, see Kaula (2000) and Celletti & Gales (2018).

In conclusion, and this choice is justified in Sect. 2.6, we can construct a basic dynamical model for the Galilean moons by considering only their mutual gravitational interactions and the influence of Jupiter (modeled as an oblate planet). Consequently, the basic differential system for the ith moon becomes$$\begin{array}{c}\hfill {\ddot{\mathit{r}}}_i=-\frac{G(m_0+m_i)}{r_i^3}{\mathit{r}}_i+{\displaystyle \sum _{j=1,j\ne i}^{N}G{m}_j}(\frac{{\mathit{r}}_j-{\mathit{r}}_i}{r_{\mathit{ij}}^3}-\frac{{\mathit{r}}_j}{r_j^3})\\ \hfill +G(m_0+m_i){\mathrm{\nabla }}_i{U}_{\overline{i}\widehat{0}}+{\displaystyle \sum _{j=1,j\ne i}^{N}G{m}_j{\mathrm{\nabla }}_j{U}_{\overline{j}\widehat{0}}·}\end{array}$$(8)

We remark that unless indicated otherwise, the oblate potentials U_ī0̂ are restricted to the zonal terms in Eq. (6), and the corresponding series are truncated at J₂. Furthermore, in light of the masses considered, we take the sums up to N = 4. However, because of the general character of Eq. (8), we prefer to keep an equally general notation, which will facilitate understanding the next subsections. Of course, these equations can be extended by considering other masses or triaxiality contributions (as in Lainey et al. 2004). Finally, we remark that the last term in Eq. (8) and the factor associated with m_i in the preceding equation represent indirect forces resulting from the oblateness of the central body. Numerically, they are relatively small (O(m_iJ₂) at most), but they allow for a better preservation of the total energy of the system. Following Lainey et al. (2004), in the remaining paper we refer to them as additional oblateness forces.

2.3. Precession of Jupiter’s line of nodes

In this subsection we provide the formulas necessary to take into account the precession of Jupiter’s mean equatorial plane in the equations above (along with the rotation around the corresponding polar axis). This effect is taken into account solely to represent the evolution of the Laplace angle according to the ephemerides set, and not for the numerical integration of the benchmark model. We still describe it here in order to provide a general form for the oblateness terms introduced in the latter.

The evolution of the Jovian rotation pole and prime meridian line relative to the J2000 frame, without taking into account nutation terms, is constructed from the 2006 IAU report (rotation) and the corresponding 2009 version (as in the Spice toolkit). Specifically, the equations are³$$\begin{array}{cc}& \alpha =268.056595-0.006499T,\hfill \\ \hfill & \delta =64.495303+0.002413T,\hfill \\ \hfill & W=284.95+870.536d.\hfill \end{array}$$(9)

In the Eq. (9), all written in degrees α stands for the right ascension, δ for the declination, and W for the longitude of the prime meridian. Furthermore, T and d denote the times in Julian centuries (36 525 days) and Julian days (86 400 s), respectively, from the standard epoch J2000.

Through the angles above, we can define the rotation from the J2000 frame to the one individualized by the directions, at a given epoch of Jupiter’s mean north pole and equatorial node. This is done, as shown in Fig. 2, by considering the Euler angles$$\begin{array}{c}\hfill I=90-\delta ,\psi =\alpha +90,\chi =W.\end{array}$$(10)

In turn, these angles define the intrinsic rotation z − x̃ − z′, with z̃′ coinciding with z′ in the figure. This rotation is defined by the matrix (see for example Lainey 2002)$$M=\left(\begin{array}{ccc}\cos \chi \cos \psi -\sin \chi \sin \psi \cos I& \cos \chi \sin \psi +\sin \chi \cos I\cos \psi & \sin \chi \sin I\\ -\sin \chi \cos \psi -\cos \chi \sin \psi \cos I& -\sin \chi \sin \psi +\cos \chi \cos I\cos \psi & \cos \chi \sin I\\ \sin \psi \sin I& -\sin I\cos \psi & \cos I\end{array}\right)\cdot$$(11)

2.4. Gradient of the oblate potential

In order to implement Eq. (8), we need an explicit, general expression for the term U_ī0̂, which has to be valid also when the precession of Jupiter’s mean equatorial plane is taken into account. The formulas that we provide here can be found in Lainey (2002), along with a more extensive discussion of their derivation.

Limiting ourselves to zonal harmonics (which are the fundamental harmonics here), we first observe that in light of the previous section, we have$$\begin{array}{c}\hfill sin{\phi }_i=\frac{z_i^{\prime }}{r_i}=\frac{x_{i}sinIsin\psi -y_{i}cos\psi sinI+z_{i}cosI}{r_i}·\end{array}$$(12)

If we denote with γ_i alternatively each of the variables (x_i; y_i; z_i), then we have$$\begin{array}{c}\hfill \frac{\partial sin{\phi }_i}{\partial {\gamma }_i}=\frac{1}{r_i}(\frac{\partial z_i^{\prime }}{\partial {\gamma }_i}-\frac{{\gamma }_{i}sin{\phi }_i}{r_i})\end{array}$$(13)

and also$$\frac{\partial }{\partial {\gamma }_i}\left[\frac{{\tilde{P}}_n(\sin {\phi }_i)}{r_i^{n + 1}}\right]=\frac{1}{r_i^{n + 1}}\left[\frac{\partial \sin {\phi }_i}{\partial {\gamma }_i}\frac{\text{d}{\tilde{P}}_n(\sin {\phi }_i)}{\text{d}\sin {\phi }_{\text{i}}}-\frac{{\tilde{P}}_n(\sin {\phi }_i)(n+1){\gamma }_i}{r_i^2}\right]·$$(14)

By joining the previous equalities, we obtain for the derivatives of Eq. (6) the expression$$\frac{\partial U_{\overline{i}\widehat{0}}^{(1)}}{\partial {\gamma }_i}=-{\displaystyle \sum _{n\hspace{0.17em}=\hspace{0.17em}2}^{\infty }\frac{R^n{J}_n}{r_i^{n\hspace{0.17em}+\hspace{0.17em}2}}}\left[\left(\frac{\partial z_i^{\prime }}{\partial {\gamma }_i}-\frac{{\gamma }_i\sin {\phi }_i}{r_i}\right)\frac{\text{d}{\tilde{P}}_n(\sin {\phi }_i)}{\text{d}\sin {\phi }_{\text{i}}}-\frac{(n+1){\gamma }_i{\tilde{P}}_n(\sin {\phi }_i)}{r_i}\right].$$(15)

In conclusion, by limiting ourselves to consider only the harmonic J₂, we obtain explicitly$$\begin{array}{cc}\hfill \frac{\partial U_{\overline{i}\widehat{0}}^{(1)}}{\partial x_i}& =\frac{R^2{J}_2}{r_i^4}[\frac{x_i}{r_i}(\frac{15}{2}{sin}^2{\phi }_i-\frac{3}{2})-3sin{\phi }_{i}sinIsin\psi ],\hfill \\ \hfill \frac{\partial U_{\overline{i}\widehat{0}}^{(1)}}{\partial y_i}& =\frac{R^2{J}_2}{r_i^4}[\frac{y_i}{r_i}(\frac{15}{2}{sin}^2{\phi }_i-\frac{3}{2})+3sin{\phi }_{i}sinIcos\psi ],\hfill \\ \hfill \frac{\partial U_{\overline{i}\widehat{0}}^{(1)}}{\partial z_i}& =\frac{R^2{J}_2}{r_i^4}[\frac{z_i}{r_i}(\frac{15}{2}{sin}^2{\phi }_i-\frac{3}{2})-3sin{\phi }_{i}cos\psi ],\hfill \end{array}$$(16)

which can be further expanded by plugging them into Eq. (12).

2.5. Conservation of the energy

The equations of motion, Eq. (8), are propagated through the use of an adaptive Runge–Kutta algorithm of seventh order, with initial conditions given by the ephemerides of the moons at the date J2000. Since in Spice these are available only up to about 2100, we propagate the conditions only up to 100 yr (so that we can also compare our results to those of Lainey et al. 2004).

To verify the correctness of the numerical integration, we considered the temporal evolution of the total energy of the system. If we denote by M the sum of the system masses, then we can express the energy E as (Pucacco & Rosquist 2017)$$\begin{array}{c}\hfill E={\displaystyle \sum _{i=1}^N\frac{m_i{\dot{\mathit{r}}}_i^2}{2}}-\frac{1}{2M}{\left({\displaystyle \sum _{i=1}^N{m}_i{\dot{\mathit{r}}}_i}\right)}^2-U,\end{array}$$(17)

where U is the potential of the system. In light of the contributions considered, we can write this term as$$\begin{array}{c}\hfill U={\displaystyle \sum _{i=1}^{N}G{m}_i{m}_0}(\frac{1}{r_i}+U_{\overline{i}\widehat{0}})+{\displaystyle \sum _{i=1}^{N-1}{\displaystyle \sum _{j=i+1}^N\frac{G{m}_j{m}_i}{r_{\mathit{ij}}}}},\end{array}$$(18)

where the terms U_ī0̂ are defined in Eq. (5) and again N = 4. Figure 3 shows that the energy remains well preserved for the entire timescale we considered.

Fig. 3. Energy relative error (in normalized units) committed when propagating Eq. (8). Units are defined by Io’s mean orbital period and semi-major axis, and by Jupiter’s mass.

2.6. Results

We evaluate here the effectiveness of the Cartesian model we introduced in capturing the librating expressions associated with the Laplace resonance. We remark that since our interest is mainly analytical, our accuracy thresholds are fairly relaxed. Furthermore, in Table 1 we provide the values of the initial conditions used to generate the plots of this subsection.

As a first step, we plot in the top panel of Fig. 4 the history of the Laplace argument as derived from the ephemerides. We do not show it here, but taking the power spectrum of this signal reveals several frequencies with higher intensity than the frequency we are interested in (which occurs slightly later than 2000 days). A low-pass filter at 1000 days is enough to remove these higher frequencies, leaving a signal with a maximum amplitude of about 0.02°. Not only is this a lower value than reported by Lieske (1977) and Musotto et al. (2002), but when we compare our plot with Fig. 2 in Musotto et al. (2002), we still have higher frequencies (again, this can be confirmed by taking a fast Fourier transform). These discrepancies may be due to several factors, ranging from the different set of ephemerides to the function used for the filtering. For reference, our procedure is performed automatically in Mathematica 11, using the default options.

Fig. 4. Top: Laplace argument ephemerides history from the J2000 epoch, including the effect of Jupiter’s mean equatorial plane precession. Bottom: result after a low-pass filter at 1000 days (tails are a byproduct of the procedure). Ephemerides from the NASA Spice toolkit.

The top panel of Fig. 5 shows that the libration of the Laplace angle, complete with the desired period, appears already when we restrict our Cartesian model to only purely gravitational terms. The amplitude is extremely large, however: about 70°. As previously argued, to recover the (almost) correct history of the angle, it is sufficient to plug in the model Jupiter’s oblateness potential truncated at the zonal harmonic J₂. The second (integration) and third (ephemerides) panels of Fig. 5 show that the approximation is quite good for the timescale we considered (100 yr). We remark, however, that the variation increases from 0.05° to 0.15° with a modulation coherent with the libration main frequency, as is apparent in the bottom plot of the same figure. This suggests that additional effects acting on the mean longitudes may have to be taken into account for longer time spans.

Fig. 5. From top to bottom: integration history of the Laplace argument obtained without oblateness terms, including Jupiter’s J2 harmonic, the ephemerides history of the argument in the same time span, and the difference of the last two plots. The starting epoch is J2000, with initial conditions given by the corresponding ephemerides.

As an additional check, useful mainly for the next sections, we estimate the linear combinations of observed mean motions corresponding to the quasi-resonant interactions Io-Europa and Europa-Ganymede. These can be computed either via a fast Fourier transform or by numerical estimation of the periods. The first method is heavily influenced by the sampling frequency, particularly for the fast-traveling Io, thus here we rely on the second method.

We determined for each moon the orbital period T by considering the first return of the mean longitude to its initial value, and derived the observed mean motion as 2π/T. This computation was then repeated over 500 periods and the average was taken as reference value. The diagnostics corresponding to these values (for model (8)) are$$\begin{array}{cc}\hfill n_1-2n_2& =0.{739926}^{\circ }/\mathrm{day},\hfill \\ \hfill n_2-2n_3& =0.{740128}^{\circ }/\mathrm{day},\hfill \end{array}$$(19)

and the error on the nominal values (1) is O(10^−4°) per day.

3. Hamiltonian formulation

3.1. Introduction

In this section we introduce a Hamiltonian model for the problem at hand. The utility of doing so is well understood, and derives from several tools developed for this formulation that allow for a deeper understanding of the underlying dynamics. A complete discussion on this is beyond the scope of the paper and is reserved for parallel works (see Celletti et al. 2018a, b). Here we focus on presenting a basic version of this model and on validating it by comparing its results with the benchmark Cartesian model.

The section is organized as follows. In Sect. 3.2 we lay down our Hamiltonian in osculating orbital elements and introduce series expansions in the neighborhoods of the two near 2:1 mean motion resonances comprising the Laplace resonance, and retain only a finite number of resonant terms. The derivation is then completed in Sect. 3.3, where we introduce expressions for the perturbative effects and the canonical set of coordinates defining the shape of the function. Finally, in Sect. 3.4 we transform this Hamiltonian by using an appropriate set of resonant coordinates to focus our attention on the resonant variables of the system.

3.2. Hamiltonian in orbital elements

As stated before, we do not discuss the derivation of the Hamiltonian, since this is done in Celletti et al. (2018a), and more details are provided in Celletti et al. (2018b). Instead, we look at its components and provide the series expansions in orbital elements leading to its final form.

With the exclusion of oblateness effects, which we examine in the next subsection, we include in our Hamiltonian only the mutual gravitational contributions of Jupiter, Io, Europa, and Ganymede. As we show below, the choice of excluding the full influence of Callisto is a necessary step to facilitate some of the calculations, although it has consequences on some aspects of the dynamics.

In order to describe the Hamiltonian, we introduce the auxiliary variables$$\begin{array}{c}\hfill M_1=m_0+m_1,{\mu }_1=\frac{m_0{m}_1}{M_1},{\kappa }_1=\frac{m_1}{M_1},\\ \hfill M_2=M_1+m_2,{\mu }_2=\frac{M_1{m}_2}{M_2},{\kappa }_2=\frac{m_2}{M_2},\\ \hfill M_3=M_2+m_3,{\mu }_3=\frac{M_2{m}_3}{M_3},{\kappa }_3=\frac{m_3}{M_3}·\end{array}$$(20)

As described in Malhotra (1991), a possible way to introduce a Hamiltonian for the system is to use Jacobi coordinates. In addition to the motion of the barycenter, this leads to three main contributions: an unperturbed part generated by Jupiter’s pull, and a perturbation due to the mutual interactions of the satellites, which we can divide into a direct and an indirect part (see Murray & Dermott 1999, for example). As shown in Celletti et al. (2018b), these include terms in the variables κ_i, which can be expanded to first order into such variables to obtain a Jovicentric approximation of the Hamiltonian.

The latter can be written in terms of the osculating orbital elements as follows. If we discard the motion of the barycenter (which is constant), the unperturbed part can be separated into three unperturbed two-body energies as$$\begin{array}{c}\hfill H_K=-\frac{G{M}_1{\mu }_1}{2a_1}-\frac{G{M}_2{\mu }_2}{2a_2}-\frac{G{M}_3{\mu }_3}{2a_3},\end{array}$$(21)

where a_i denotes the osculating semi-major axis of the ith moon. For the direct and indirect terms of the perturbation, we can work in a neighborhood of the exact 2:1 mean motion commensurabilities corresponding to the quasi-resonant interactions Io-Europa and Europa-Ganymede and, after some transformations, obtain an expansion up to second order in the eccentricities. We remark that in this expansion we chose to retain only low-frequency terms associated with the resonant angles. Furthermore, it is important to note that we are working with two separate quasiresonances constituting the Laplace resonance, and not with the latter itself. As we show in Sect. 3.4, it is possible to introduce this argument with a suitable canonical change of variables.

We obtain for the perturbing terms the formulae (as mentioned, compare with Malhotra 1991 and Murray & Dermott 1999)$$\begin{array}{cc}& H_P^{1,2}=-\frac{G{m}_1{m}_2}{a_2}(B_0({\alpha }_{1,2})+f_1^{1,2}{e}_{1}cos(2{\lambda }_2-{\lambda }_1-{\varpi }_1)\hfill \\ \hfill & +f_2^{1,2}{e}_{2}cos(2{\lambda }_2-{\lambda }_1-{\varpi }_2)+f_3^{1,2}(e_1^2+e_2^2)\hfill \\ \hfill & +f_4^{1,2}{e}_1{e}_{2}cos({\varpi }_1-{\varpi }_2)+f_5^{1,2}{e}_1{e}_{2}cos(4{\lambda }_2-2{\lambda }_1-{\varpi }_1-{\varpi }_2)\hfill \\ \hfill & +f_6^{1,2}{e}_1^{2}cos(4{\lambda }_2-2{\lambda }_1-2{\varpi }_1)+f_7^{1,2}{e}_2^{2}cos(4{\lambda }_2-2{\lambda }_1-2{\varpi }_2)),\hfill \\ \hfill & H_P^{2,3}=-\frac{G{m}_2{m}_3}{a_3}(B_0({\alpha }_{2,3})+f_1^{2,3}{e}_{2}cos(2{\lambda }_3-{\lambda }_2-{\varpi }_2)\hfill \\ \hfill & +f_2^{2,3}{e}_{3}cos(2{\lambda }_3-{\lambda }_2-{\varpi }_3)+f_3^{2,3}(e_2^2+e_3^2)\hfill \\ \hfill & +f_4^{2,3}{e}_2{e}_{3}cos({\varpi }_2-{\varpi }_3)+f_5^{2,3}{e}_2{e}_{3}cos(4{\lambda }_3-2{\lambda }_2-{\varpi }_2-{\varpi }_3)\hfill \\ \hfill & +f_6^{2,3}{e}_2^{2}cos(4{\lambda }_3-2{\lambda }_2-2{\varpi }_2)+f_7^{2,3}{e}_3^{2}cos(4{\lambda }_3-2{\lambda }_2-2{\varpi }_3)),\hfill \\ \hfill & H_P^{1,3}=-\frac{G{m}_1{m}_3}{a_3}(B_0({\alpha }_{1,3})+f_3^{1,3}(e_1^2+e_3^2)\hfill \\ \hfill & +f_4^{1,3}{e}_1{e}_{3}cos({\varpi }_1-{\varpi }_3))·\hfill \end{array}$$(22)

The functions $f_k^{i,j}$ are linear combinations of the Laplace coefficients $b_s^{\left(n\right)}$ and the associated derivatives (see, e.g., Murray & Dermott 1999 or Showman & Malhotra 1997), which in turn are functions of the semi-major axis ratios ${\alpha }_{i,j}=\frac{a_i}{a_j}$. The terms B₀(α_i,j) equal $\frac{1}{2}{b}_{1/2}^{(0)}({\alpha }_{i,j})-1$. We also note that all the inclinations have been set to zero, and therefore the model we consider is planar. This assumption is reasonable (although it has some consequences) since the inclinations are on the order of 10^−1° at most.

3.3. Canonical elements and the perturbative effects

To conclude our construction, we need to select an appropriate set of canonical variables in which to express the Hamiltonian above. For this work we consider the modified Delaunay elements defined by the equations (Malhotra 1991)$$\begin{array}{cc}& {\lambda }_i,L_i={\mu }_i\sqrt{G{M}_i{a}_i}\simeq m_i\sqrt{G{m}_0{a}_i},\hfill \\ \hfill & -{\varpi }_i,G_i\equiv L_i(1-\sqrt{1-e_i^2})\simeq \frac{1}{2}{L}_i{e}_i^2.\hfill \end{array}$$(23)

Because we consider only three moons, the final Hamiltonian has six degrees of freedom. Thanks to our choice in the terms to retain, and with the change of coordinates in the next section, it is possible to reduce the complexity of this Hamiltonian even more.

The same type of steps that we have taken for the gravitational part of the Hamiltonian has to be considered for the perturbations taken into account. For example, we showed in the previous section that Jupiter’s oblateness potential, truncated at the J₂ harmonic, is the single most important factor for recovering the amplitude of the Laplace argument. Since we are interested in recovering only the latter’s main frequency, for the Hamiltonian model we only consider the secular part of the perturbation provided by the J₂ harmonic. This can be expressed in terms of the osculating orbital elements as$$\begin{array}{c}\hfill H_J=-{\displaystyle \sum _{i=1}^3\frac{G{M}_i{\mu }_i}{2a_i}}[J_2{\left(\frac{R}{a_i}\right)}^2(1+\frac{3}{2}{e}_i^2)-\frac{3}{4}],\end{array}$$(24)

and is trivially converted via the canonical elements introduced before. Similarly, we can reintroduce Callisto as a perturbation. If we only consider its secular part, we have$$\begin{array}{c}\hfill H_{\mathrm{C}}=-{\displaystyle \sum _{i=1}^3\frac{G{m}_i{m}_4}{a_4}}[B_0\left({\alpha }_{i,4}\right)+\frac{{\alpha }_{i,4}}{8}{b}_{3/2}^{(1)}\left({\alpha }_{i,4}\right)(e_i^2+e_{\mathrm{C}}^2)].\end{array}$$(25)

Unless explicitly mentioned, this last perturbation is not considered in the Hamiltonian model, since the Laplace resonance does not directly involve Callisto.

3.4. Resonant coordinates

As described above, we can introduce a canonical change of coordinates to highlight the Laplace argument as a variable of the system. This is defined in terms of the modified Delaunay elements (23) as$$\begin{array}{cc}\hfill q_1& =2{\lambda }_2-{\lambda }_1-{\varpi }_1,P_1=G_1,\hfill \\ \hfill q_2& =2{\lambda }_2-{\lambda }_1-{\varpi }_2,P_2=G_2,\hfill \\ \hfill q_3& =2{\lambda }_3-{\lambda }_2-{\varpi }_3,P_3=G_3,\hfill \\ \hfill q_4& =3{\lambda }_2-2{\lambda }_3-{\lambda }_1,P_4=\frac{1}{3}(L_2-2(G_1+G_2)+G_3),\hfill \\ \hfill q_5& ={\lambda }_1-{\lambda }_3,P_5=\frac{1}{3}(3L_1+L_2+G_1+G_2+G_3),\hfill \\ \hfill q_6& ={\lambda }_3,P_6=L_1+L_2+L_3-G_1-G_2-G_3.\hfill \end{array}$$(26)

Equation (22) shows that with the terms that we decided to retain in our expansions, the variables q₅ and q₆ do not appear in the Hamiltonian (P₆ is the total angular momentum of the system; Laskar & Petit 2017). Under the previous change of coordinates, the Hamiltonian H(q_i, P_i) therefore does possess four d.o.f.

Furthermore, and this is explained in greater detail in Celletti et al. (2018a) or Broer & Hanßmann (2016), it can be proved that the Hamiltonian, with the gravitational terms restricted to first order in the eccentricity, has exactly eight equilibrium configurations (excluding the symmetries of the system). The only stable configuration, whose discovery is generally attributed to de Sitter (1931), is given by$$\begin{array}{c}\hfill \{q_1,q_2,q_3,q_4\}=\{0,\pi ,\pi ,\pi \}.\end{array}$$(27)

This equilibrium does not correspond to the actual state of the system, since we know from observations that the angle q₃ rotates. This is an important detail and is exploited in the next section to facilitate the search for the correct amplitude of the resonant argument.

4. Numerical simulations

4.1. Comparison of the models

Here we numerically compare the Cartesian and Hamiltonian formulations in order to measure how well the latter approximates the main amplitude and frequency of the Laplace argument, and the Jovian moons’ dynamics (except for Callisto, of course). The initial conditions are the same for the two models and correspond to the ephemerides at the date J2000, which we also used in Sect. 2.6. For the Hamiltonian model, they are expressed in terms of the osculating orbital elements, whose values are listed in Table 3.

As a first step, in Fig. 6 we superimpose the Hamiltonian evolution of the Laplace argument on the Cartesian evolution presented in Fig. 5. In the top panel, where we exclude the oblateness contribution for both models, the difference is relatively small, both in period and amplitude. However, if we retain the J₂ terms (bottom panel), the corresponding amplitude is out of scale of at least an order of magnitude, even though the main period is still well determined by the Hamiltonian formulation. Of course, it is to be expected that the same initial conditions do not lead to the same evolution of the argument, and it is still possible to recover the correct amplitude, if with some compromise. We show this in the next subsection.

Fig. 6. Top: Laplace argument integration history obtained without oblateness terms. Bottom: same evolution with the addition of the J2 terms. In red we show the results from integration of the Cartesian model, in green equivalent from the Hamiltonian.

We examine the Jovian moons’ dynamics stemming from the Hamiltonian’s approximation in more detail. In particular, we consider the “shape” of the orbit, that is, the evolution of the moons’ eccentricities and semi-major axes. In Fig. 7 we plot the history of the former for the Hamiltonian formulation (top panel, H_C not included) and the Cartesian one (bottom panel). In all cases, the absolute error achieves a maximum value of O(10⁻³). In turn, this implies different values for the maximum relative error, with peaks of 40% for Io, a steady increase of up to 15% for Europa, and occasional ventures beyond 100% for Ganymede (although the average is about 50−60%). Furthermore, it is evident that an additional frequency in Ganymede’s eccentricity evolution is present in the Cartesian case. Particularly for this last effect, the primary suspect might be Callisto, whose gravitational contribution is not taken into account in Fig. 7. Introducing its secular effect H_C in the Hamiltonian formulation does not significantly change the results, however. Thus, it seems that this lack of accuracy is a price to pay if the full effect of Callisto is not taken into account (and therefore the number of d.o.f. is kept small).

Fig. 7. Top: Jovian moons’ eccentricity evolution as computed with the Hamiltonian model. Bottom: corresponding histories derived from the Cartesian model The moons are color-coded (Io in red, Europa in green, and Ganymede in blue).

In Fig. 8 we consider the absolute error of the history of the semi-major axes as obtained from the two formulations. By dividing them for the semi-major axis values obtained with the Cartesian integration, it is easy to check that the relative error committed has an upper bound of about 0.001%. Naturally, since we approximate the real system with a planar system and the inclinations are more or less on the same order, the largest absolute error is committed for Ganymede, which is the moon most distant from Jupiter.

Fig. 8. History of the absolute errors on the Jovian moons’ semi-major axes as computed in the Cartesian and Hamiltonian formulations. The moons are color-coded (Io in red, Europa in green, and Ganymede in blue).

Finally, in Fig. 9 we analyze the different model histories for the angle q₃, whose character, as stated in Sect. 3.4 and as shown in the next subsection, determines whether we are in the Laplace resonance. In the first case, this angle rotates with a period of about 485 days, as apparent in the top plot from its Cartesian evolution (red). Small period oscillations are present, but they are averaged out when we consider the same angle in the Hamiltonian case (green). Furthermore, as shown in the middle plot, the agreement of the two models on this angle deteriorates with time, which implies a small error in the related drift rate for the two formulations. We remark that as the bottom plot suggests, a similar, less severe displacement exists between the Cartesian model and the ephemerides set. Moreover, these accuracy losses are comparable in terms of magnitude of the relative error with those exhibited by the Laplace argument, as it is apparent by matching the bottom panels of Figs. 5 and 6 with the equivalent panels in Fig. 9.

Fig. 9. From top to bottom: Cartesian and Hamiltonian histories for the angle q3, corresponding error on the evolution, and similar error between the Cartesian formulation and the ephemerides set.

For reference, points in the plots are taken each day, with the exception of the pictures for the semi-major axes, where they are taken every 10 days, and the top plot of Fig. 9, where they are taken every hour. This convention is also maintained in the next subsection, where we show why the main resonant behavior is not recovered and how this can be remedied.

4.2. 4.2. Conditions for a “correct” Hamiltonian resonance

We showed in the previous section that taking the ephemerides as initial conditions for the Hamiltonian model does not lead to a correct approximation of the Laplace argument amplitude. Similarly, the observed quasi-resonant mean motion combinations slightly disagree with the nominal values (1). In particular, the average over a period of 100 yr gives$$\begin{array}{cc}\hfill n_1-2n_2& =0.{734517}^{\circ }/\mathrm{day},\hfill \\ \hfill n_2-2n_3& =0.{734591}^{\circ }/\mathrm{day},\hfill \end{array}$$(28)

with an error of about 0.005°/day. This is reasonable, since there are several approximations in the Hamiltonian model, ranging from its planarity to the more important restrictions on the selected resonant terms.

Fortunately, we can exploit its dynamics, as seen from the lens of the resonant variables in Eq. (26), to revert the problem and determine the initial conditions that generate values closer to the nominal ones. The first step, as shown in Celletti et al. (2018a), is to restrict our resonant Hamiltonian to first order in the eccentricity in the gravitational terms, and retain the full secular contribution of Jupiter’s oblateness in Eq. (24), so that we can employ normal form techniques to approximate the dynamics on the phase plane (q₃, P₃; similarly to Henrard 1984 and Bucciarelli et al. 2016).

As mentioned in Sect. 3.4, the first-order Hamiltonian has exactly one stable equilibrium, and the difference between this (“de Sitter”) equilibrium and the actual state of the system lies in the rotating character of the angle q₃. Correspondingly, in Fig. 10, where the phase space of the normal form Hamiltonian is plotted, the first appears as an elliptic fixed point, while the second lies in the vicinity of one of the rotational tori (blue curves).

Fig. 10. Phase plane for the normal form Hamiltonian obtained by restricting Eq. (22) to first-order terms in the eccentricity and retaining the full secular effect of Jupiter’s J2 harmonic in Eq. (24) (along with Eq. (21), of course). The full blue curve is obtained by integrating Hamilton’s equations with Eq. (30) as initial conditions, except for q3 = 0 and P3 as in Eq. (31). These are used also for the dashed blue curve, where Eq. (25) and the second-order terms of Eq. (22) are also taken into account.

The two blue curves are both computed from direct integration of Hamilton’s equations; the full curve corresponds to the sum of the terms in Eqs. (21) and (22) (restricted to first order in the eccentricity), and (24), all translated into the coordinates of Eq. (23). The dashed curve instead adds to the previous addenda the terms at second order in the eccentricity for Eq. (22) and the secular contribution in Eq. (25) of Callisto. The red curve is outlined only to remark its boundary status between the rotational and librational regimes (it is not a formal separatrix). The normal form Hamiltonian associated with the phase plane is given by$$\begin{array}{cc}& H_{\mathrm{d}}S=-0.003638P_3-1.6497P_3^2\hfill \\ \hfill & -6.1362\times {10}^{-6}\sqrt{P_3}cos{q}_3+1.1155\times {10}^{-5}{P}_3^{\frac{3}{2}}cos{q}_3,\hfill \end{array}$$(29)

where the units of measure are determined by Io’s semi-major axis, its mean orbital period in days, and by Jupiter’s mass.

The blue curves were determined by starting from the numerical location of the de Sitter equilibrium, which is obtained from Hamilton’s equations evaluated at {q₁, q₂, q₃, q₄} = {0, π, π, π} by employing a root-finding algorithm for the momenta P_i. Differently from Celletti et al. (2018a), for this computation we derived the de Sitter equilibrium by also including the terms at second order in the eccentricity in the gravitational terms (those pertaining to the oblateness remained unchanged). The values obtained in this way were subsequently converted into Delaunay elements L_i and G_i. From these, we can derive the approximate location of the equilibrium in terms of the orbital elements, that is,$$\begin{array}{cc}\hfill a_1^{\star }& =422043.201\mathrm{km},e_1^{\star }=0.00422405,\hfill \\ \hfill a_2^{\star }& =671235.280\mathrm{km},e_2^{\star }=0.00948442,\hfill \\ \hfill a_3^{\star }& =1070501.267\mathrm{km},e_3^{\star }=0.00063468.\hfill \end{array}$$(30)

Furthermore, we point out that in addition to having been computed with additional terms in the Hamiltonian, these values differ from those in Celletti et al. (2018a) because the initial conditions there are taken from Lainey et al. (2004).

As mentioned before, the “real” Laplace resonance is separated from the de Sitter equilibrium by the rotation and not libration of the q₃ argument. It is therefore sufficient to vary the action P₃ and therefore the eccentricity of Ganymede (if a fixed semi-major axis is assumed) to distinguish these two. We found that a good value for recovering the main amplitude of the Laplace argument is given by$$\begin{array}{c}\hfill P_3={10}^{-8},\end{array}$$(31)

where the unit of length is assumed normalized with respect to Io’s semi-major axis. This value is about two orders of magnitude lower than the value in Celletti et al. (2018a; even accounting for the difference in Io’s semi-major axis). Furthermore, to avoid being captured in the librational regime, we need to chose the initial value of q₃ properly. We here selected q₃ = 0.

In Fig. 11 we plot the time evolution of the Hamiltonian Laplace argument, with initial conditions corresponding to the de Sitter equilibrium at second order in the eccentricity, except for P₃ and q₃, which are taken as above. The results are clearly greatly improved with respect to Fig. 6, with the amplitudes of both the normal and filtered angle very close to those in Fig. 4. The main difference lies in the lack of additional frequencies after filtering, since they have been smoothed out in the Hamiltonian construction. We remark that key to this accurate result is having taken as initial conditions those corresponding to the de Sitter equilibrium at second order in the eccentricity. Using the first-order equilibrium leads to an amplitude of one order of magnitude more for both the normal and the filtered angle.

Fig. 11. Top: Hamiltonian history of the Laplace argument obtained from the initial conditions of the de Sitter equilibrium, Eq. (30), except for q3 = 0 and P3 as in Eq. (31). Bottom: corresponding filtering at 1000 days.

Interestingly enough, while the diagnostics do vary, the error with respect to the nominal values in Eq. (1) remains similar to Eq. (28). With the initial conditions above, the average over a period of 100 yr becomes$$\begin{array}{cc}\hfill n_1-2n_2& =0.{733692}^{\circ }/\mathrm{day},\hfill \\ \hfill n_2-2n_3& =0.{733676}^{\circ }/\mathrm{day}.\hfill \end{array}$$(32)

Several factors may contribute to this discrepancy, and it is difficult to pinpoint an exact cause.

To summarize, by exploiting the Hamiltonian dynamics of the system, we have highlighted a way to recover the correct values for the quantities we searched for. However, a quick comparison with the different conditions and results of Celletti et al. (2018a) shows that this procedure (and therefore the Hamiltonian model) is highly sensitive to the initial conditions. It clearly shows, however, that the approximate nature of the model requires a tradeoff when certain values are searched for. In particular, if the nominal evolution of the Laplace resonance is searched for, the Hamiltonian “shape” of the orbits (e.g., semi-major axes and eccentricities) cannot be very close to the real shape.

5. Conclusion

We have considered the problem of constructing approximate models for the dynamics of the Galilean moons in the Jovian system, with a particular focus on recovering the time evolution of the Laplace resonance. To do this, we considered as benchmark a set of ephemerides extracted from NASA’s Spice toolkit, that we defined as elements to recover the corresponding amplitude and period of the main frequency of the Laplace argument, and the related quasi-resonant mean motion relations of the pairs Io-Europa and Europa-Ganymede. Subsequently, we introduced a Jovicentric Cartesian model of the system following Lainey (2002), and exploited it to numerically show that within the scale of 100 yr, a sufficient approximation of the Laplace argument can be obtained by considering only the mutual gravitational interactions among the moons, along with Jupiter’s oblateness potential truncated at the harmonic J₂ (and, of course, including the point-mass gravitational term).

Starting from the information obtained from the Cartesian formulation, we introduced an alternative Hamiltonian model that is more suitable to analytically delve into the mechanisms of the resonance. In order to contain its complexity, we considered the Jovian system as planar, concentrated on the dynamics of the moons involved in the resonance, and restricted the oblateness perturbation to its secular part.

Numerical comparisons with the Cartesian model showed that if the same initial conditions are used, the Hamiltonian formulation determines the main period of the Laplace argument, but it fails to recover the associated amplitude by an order of magnitude. This is reflected in the mean motion diagnostics, which presents an absolute error on the order of 10^−3°/day with respect to the corresponding Cartesian formulation.

By introducing appropriate resonant coordinates, it is apparent that the actual Laplace configuration is one rotating angle away from the only stable equilibrium of the Hamiltonian. Thus, by exploiting normal form techniques to project the Hamiltonian flow onto the appropriate action-angle phase plane, we first numerically located the stable equilibrium, then moved the action variable up to the rotational regime in search of a value of this variable that would generate a good approximation of the amplitude of the Laplace argument. We used a much lower value than Celletti et al. (2018a). Additionally, the equilibrium was computed using the full Hamiltonian in place of restricting non-oblateness terms to first order in the eccentricity. These two factors combined led to an excellent approximation of the Laplace angle, with an amplitude comparable to the values reported in the literature. Following the chain of canonical transformations backward shows that a good tradeoff to obtain the desired values with this procedure lies in Ganymede’s eccentricity.

Parallel and future work involve a thorough use of the Hamiltonian introduced here to determine the mechanisms of multibody resonances. In particular, Celletti et al. (2018b) considered pairs of commensurabilities different from the two 2:1 studied here, and we classify the dynamics on the base of the eccentricity order associated to the relevant resonant terms. This is a necessary step to study systems different from the Jovian one, which of course means extending the study beyond our solar system. Delving deeper into the mechanisms of these resonances implies extending the time scales considered here, which in turn means considering additional conservative and dissipative perturbations for the models.

Acknowledgments

We thank the Italian Space Agency for partially supporting this work through the grant 2013-056-RO.1, and we acknowledge the GNFM/INdAM. We are also grateful to S. Ferraz-Mello for the many fruitful discussions, to L. Iess for his support, and the referee for their useful comments. Finally, F.P. thanks M.J. Mariani, F. De Marchi, and D. Durante for their help with the SPICE toolkit, while A.C. acknowledges the MIUR Excellence Department Project awarded to the Department of Mathematics of the University of Rome “Tor Vergata” (CUP E83C18000100006).

References

Appendix A: Explicit Laplace coefficients

For comparison reasons, in this short appendix we provide an explicit expression for the Hamiltonian terms containing any type of Laplace coefficient (computed through series expansions), that is, the perturbative functions in Eq. (22) and Callisto’s contribution, Eq. (25). The formulae used for the coefficients can be found in the references mentioned in the appropriate sections, while the values used for the computations are given in the tables of this paper. Of course, the numbers are truncated at a significant enough digit.

The perturbative terms in Eq. (22) are$$\begin{array}{cc}\hfill H_{\mathrm{P}}^{1,2}& =-\frac{G{m}_1{m}_2}{a_2}(0.12954+0.384233e_{1}cos(2{\lambda }_2-{\lambda }_1-{\varpi }_1)\hfill \\ \hfill & -1.18409e_{2}cos(2{\lambda }_2-{\lambda }_1-{\varpi }_2)+0.424985(e_1^2+e_2^2)\hfill \\ \hfill & +1.67945e_1{e}_{2}cos({\varpi }_1-{\varpi }_2)\hfill \\ \hfill & +3.57298e_1{e}_{2}cos(4{\lambda }_2-2{\lambda }_1-{\varpi }_1-{\varpi }_2)\hfill \\ \hfill & -4.92864cos(4{\lambda }_2-2{\lambda }_1-2{\varpi }_1)\hfill \\ \hfill & -0.569703e_2^{2}cos(4{\lambda }_2-2{\lambda }_1-2{\varpi }_2)),\hfill \\ \hfill H_{\mathrm{P}}^{2,3}& =-\frac{G{m}_2{m}_3}{a_3}(0.12861+0.379622e_{2}cos(2{\lambda }_3-{\lambda }_2-{\varpi }_2)\hfill \\ \hfill & -1.17533e_{3}cos(2{\lambda }_3-{\lambda }_2-{\varpi }_3)+0.420341(e_2^2+e_3^2)\hfill \\ \hfill & +1.65729e_2{e}_{3}cos({\varpi }_2-{\varpi }_3)\hfill \\ \hfill & +3.54454e_2{e}_{3}cos(4{\lambda }_3-2{\lambda }_2-{\varpi }_2-{\varpi }_3)\hfill \\ \hfill & -4.87651e_2^{2}cos(4{\lambda }_3-2{\lambda }_2-2{\varpi }_2)\hfill \\ \hfill & -0.561562e_3^{2}cos(4{\lambda }_3-2{\lambda }_2-2{\varpi }_3)),\hfill \\ \hfill H_{\mathrm{P}}^{1,3}& =-\frac{G{m}_1{m}_3}{a_3}(0.04267+0.0793481(e_1^2+e_3^2)\hfill \\ \hfill & +0.206309e_1{e}_{3}cos({\varpi }_1-{\varpi }_3)),\hfill \end{array}$$(A.1)

and Callisto’s contribution, Eq. (25), can be split into three separate terms (one each for Io, Europa, and Ganymede) as$$\begin{array}{cc}& H_{\mathrm{C}}^1=-\frac{G{m}_1{m}_4}{a_4}[0.01293+0.0207567(e_i^2+e_{\mathrm{C}}^2)],\hfill \\ \hfill & H_{\mathrm{C}}^2=-\frac{G{m}_2{m}_4}{a_4}[0.03427+0.0615473(e_i^2+e_{\mathrm{C}}^2)],\hfill \\ \hfill & H_{\mathrm{C}}^3=-\frac{G{m}_3{m}_4}{a_4}[0.09993+0.253539(e_i^2+e_{\mathrm{C}}^2)].\hfill \end{array}$$(A.2)

All Tables

All Figures

	Fig. 1. Jovicentric snapshots of the Jovian system, with T denoting Io’s mean orbital period. The moons are color-coded (Io in red, Europa in green, and Ganymede in blue).
In the text

	Fig. 2. Precession and rotation of Jupiter’s mean equatorial frame (P0, x′, y′, z′) with respect to the J2000 frame (P0, x, y, z), represented through the Euler angles (ψ, I, χ).
In the text

	Fig. 3. Energy relative error (in normalized units) committed when propagating Eq. (8). Units are defined by Io’s mean orbital period and semi-major axis, and by Jupiter’s mass.
In the text

	Fig. 4. Top: Laplace argument ephemerides history from the J2000 epoch, including the effect of Jupiter’s mean equatorial plane precession. Bottom: result after a low-pass filter at 1000 days (tails are a byproduct of the procedure). Ephemerides from the NASA Spice toolkit.
In the text

	Fig. 5. From top to bottom: integration history of the Laplace argument obtained without oblateness terms, including Jupiter’s J2 harmonic, the ephemerides history of the argument in the same time span, and the difference of the last two plots. The starting epoch is J2000, with initial conditions given by the corresponding ephemerides.
In the text

	Fig. 6. Top: Laplace argument integration history obtained without oblateness terms. Bottom: same evolution with the addition of the J2 terms. In red we show the results from integration of the Cartesian model, in green equivalent from the Hamiltonian.
In the text

	Fig. 7. Top: Jovian moons’ eccentricity evolution as computed with the Hamiltonian model. Bottom: corresponding histories derived from the Cartesian model The moons are color-coded (Io in red, Europa in green, and Ganymede in blue).
In the text

	Fig. 8. History of the absolute errors on the Jovian moons’ semi-major axes as computed in the Cartesian and Hamiltonian formulations. The moons are color-coded (Io in red, Europa in green, and Ganymede in blue).
In the text

	Fig. 9. From top to bottom: Cartesian and Hamiltonian histories for the angle q3, corresponding error on the evolution, and similar error between the Cartesian formulation and the ephemerides set.
In the text

Fig. 10. Phase plane for the normal form Hamiltonian obtained by restricting Eq. (22) to first-order terms in the eccentricity and retaining the full secular effect of Jupiter’s J2 harmonic in Eq. (24) (along with Eq. (21), of course). The full blue curve is obtained by integrating Hamilton’s equations with Eq. (30) as initial conditions, except for q3 = 0 and P3 as in Eq. (31). These are used also for the dashed blue curve, where Eq. (25) and the second-order terms of Eq. (22) are also taken into account.

In the text

	Fig. 11. Top: Hamiltonian history of the Laplace argument obtained from the initial conditions of the de Sitter equilibrium, Eq. (30), except for q3 = 0 and P3 as in Eq. (31). Bottom: corresponding filtering at 1000 days.
In the text