Issue 
A&A
Volume 617, September 2018



Article Number  A35  
Number of page(s)  12  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201832856  
Published online  20 September 2018 
Element history of the Laplace resonance: a dynamical approach
^{1}
Department of Mathematics, University of Rome Tor Vergata, Via della Ricerca Scientifica 1, 00133
Rome, Italy
email: f.paita@live.it; celletti@mat.uniroma2.it
^{2}
Department of Physics, University of Rome Tor Vergata, Via della Ricerca Scientifica 1, 00133
Rome, Italy
email: pucacco@roma2.infn.it
Received:
19
February
2018
Accepted:
5
June
2018
Context. We consider the threebody mean motion resonance defined by the Jovian moons Io, Europa, and Ganymede, which is commonly known as the Laplace resonance. In terms of the moons’ mean longitudes λ_{1} (Io), λ_{2} (Europa), and λ_{3} (Ganymede), this resonance is described by the librating argument φ_{L} ≡ λ_{1} − 3λ_{2} + 2λ_{3} ≈ 180°, which is the sum of φ_{12} ≡ λ_{1} − 2λ_{2} + ϖ_{2} ≈ 180° and φ_{23} ≡ λ_{2} − 2λ_{3} + ϖ_{2} ≈ 0°, where ϖ_{2} denotes Europa’s longitude of perijove.
Aims. In particular, we construct approximate models for the evolution of the librating argument φ_{L} over the period of 100 yr, focusing on its principal amplitude and frequency, and on the observed mean motion combinations n_{1} − 2n_{2} and n_{2} − 2n_{3} associated with the quasiresonant interactions above.
Methods. First, we numerically propagated the Cartesian equations of motion of the Jovian system for the period under examination, and by comparing the results with a suitable set of ephemerides, we derived the main dynamical effects on the target quantities. Using these effects, we built an alternative Hamiltonian formulation and used the normal forms theory to precisely locate the resonance and to semianalytically compute its main amplitude and frequency.
Results. From the Cartesian model we observe that on the timescale considered and with ephemerides as initial conditions, both φ_{L} and the diagnostics n_{1} − 2n_{2} and n_{2} − 2n_{3} are well approximated by considering the mutual gravitational interactions of Jupiter and the Galilean moons (including Callisto), and the effect of Jupiter’s J_{2} harmonic. Under the same initial conditions, the Hamiltonian formulation in which Callisto and J_{2} are reduced to their secular contributions achieves larger errors for the quantities above, particularly for φ_{L}. By introducing appropriate resonant variables, we show that these errors can be reduced by moving in a certain actionangle phase plane, which in turn implies the necessity of a tradeoff in the selection of the initial conditions.
Conclusions. In addition to being a good starting point for a deeper understanding of the Laplace resonance, the models and methods described are easily generalizable to different types of multibody mean motion resonances. Thus, they are also prime tools for studying the dynamics of extrasolar systems.
Key words: celestial mechanics / planets and satellites: dynamical evolution and stability / methods: analytical / methods: numerical
© 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 socalled Laplace resonance, a threebody 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 quasiresonant interactions, which can be described in terms of the observed mean motions (n_{1} for Io, n_{2} for Europa, and n_{3} for Ganymede) as (Brown 1977)(1)
Fig. 1. Jovicentric snapshots of the Jovian system, with T denoting Io’s mean orbital period. The moons are colorcoded (Io in red, Europa in green, and Ganymede in blue). 

Open with DEXTER 
In addition to these two, a weaker resonant interaction exists between Ganymede and Callisto, quantifiable as 3n_{3} − 7n_{4} = −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(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} ≡ λ_{1} − 3λ_{2} + 2λ_{3} ≈ 180°. Furthermore, because ϖ_{1} and ϖ_{2} drift with similar rates (Yoder & Peale 1981), the system is affected by the additional librating relation λ_{1} − 2λ_{2} + ϖ_{1} ≈ 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 shortperiod terms, and he then solved the resulting differential equations for the longperiod effects by successive substitutions. Later on, Sagnier (1975) and FerrazMello (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 shortperiod terms in his nonHamiltonian 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_{2} 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 semimajor axes. Amplitude and frequency of the Laplace argument, as well as the diagnostics n_{1} − 2n_{2} and n_{2} − 2n_{3}, 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 quasiresonant observed mean motion combinations n_{1} − 2n_{2} and n_{2} − 2n_{3}.
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 ^{1}, 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 xaxis 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 socalled 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^{2}.
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 zaxis of this frame, while the line of nodes obtained by intersecting the previous equatorial plane with the J2000 plane acts as the associated xdirection. We note that this new frame is noninertial, 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 (P_{0}, x′, y′, z′) with respect to the J2000 frame (P_{0}, x, y, z), represented through the Euler angles (ψ, I, χ). 

Open with DEXTER 
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_{0}, Io with P_{1}, Europa with P_{2}, Ganymede with P_{3}, and Callisto with P_{4}.
Let r_{i} and r_{ij} denote the relative vectors and , 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 , we can derive the equations describing the dynamics in a Jovicentric frame (P_{0}, x, y, z). These read(3)
where m_{i} and m_{0} denote the mass of the moon and the mass of Jupiter, respectively, while F_{i} and F_{0} 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 , the gravitational attraction exerted by P_{j} over P_{i} becomes F_{ij} = Gm_{i}m_{j}∇_{i}U_{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(4)
where F_{ī0̄} = Gm_{i}m_{0}∇_{i}U_{ī0̄} and F_{ī0̂} = Gm_{i}m_{0}∇_{i}U_{ī0̂}. This last term represents the action of the P_{0} 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_{0}, the potential U_{ī0̂} describing the previous action can be written as(5)
The quantities J_{n}, C_{nm}, and S_{nm} are all constants depending on the particular primary, while denotes the associated Legendre polynomials (with P̃_{n} equal to ). For an indepth 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(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_{2}. 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_{i}J_{2}) 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^{3}(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(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)(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(12)
If we denote with γ_{i} alternatively each of the variables (x_{i}; y_{i}; z_{i}), then we have(13)
By joining the previous equalities, we obtain for the derivatives of Eq. (6) the expression(15)
In conclusion, by limiting ourselves to consider only the harmonic J_{2}, we obtain explicitly(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)(17)
where U is the potential of the system. In light of the contributions considered, we can write this term as(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 semimajor axis, and by Jupiter’s mass. 

Open with DEXTER 
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.
SPICE Cartesian elements for the Galilean moons in the Jovicentric equatorial frame.
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 lowpass 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 lowpass filter at 1000 days (tails are a byproduct of the procedure). Ephemerides from the NASA Spice toolkit. 

Open with DEXTER 
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_{2}. 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 J_{2} 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. 

Open with DEXTER 
As an additional check, useful mainly for the next sections, we estimate the linear combinations of observed mean motions corresponding to the quasiresonant interactions IoEuropa and EuropaGanymede. 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 fasttraveling 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(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(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 twobody energies as(21)
where a_{i} denotes the osculating semimajor 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 quasiresonant interactions IoEuropa and EuropaGanymede 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 lowfrequency 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)(22)
The functions are linear combinations of the Laplace coefficients and the associated derivatives (see, e.g., Murray & Dermott 1999 or Showman & Malhotra 1997), which in turn are functions of the semimajor axis ratios . The terms B_{0}(α_{i,j}) equal . 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)(23)
Parameters for the simulations.
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_{2} 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_{2} harmonic. This can be expressed in terms of the osculating orbital elements as(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(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(26)
Equation (22) shows that with the terms that we decided to retain in our expansions, the variables q_{5} and q_{6} do not appear in the Hamiltonian (P_{6} 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(27)
This equilibrium does not correspond to the actual state of the system, since we know from observations that the angle q_{3} 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_{2} 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 J_{2} terms. In red we show the results from integration of the Cartesian model, in green equivalent from the Hamiltonian. 

Open with DEXTER 
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 semimajor 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^{−3}). 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 colorcoded (Io in red, Europa in green, and Ganymede in blue). 

Open with DEXTER 
In Fig. 8 we consider the absolute error of the history of the semimajor axes as obtained from the two formulations. By dividing them for the semimajor 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’ semimajor axes as computed in the Cartesian and Hamiltonian formulations. The moons are colorcoded (Io in red, Europa in green, and Ganymede in blue). 

Open with DEXTER 
Finally, in Fig. 9 we analyze the different model histories for the angle q_{3}, 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 q_{3}, corresponding error on the evolution, and similar error between the Cartesian formulation and the ephemerides set. 

Open with DEXTER 
For reference, points in the plots are taken each day, with the exception of the pictures for the semimajor 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 quasiresonant mean motion combinations slightly disagree with the nominal values (1). In particular, the average over a period of 100 yr gives(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_{3}, P_{3}; similarly to Henrard 1984 and Bucciarelli et al. 2016).
As mentioned in Sect. 3.4, the firstorder 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_{3}. 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 firstorder terms in the eccentricity and retaining the full secular effect of Jupiter’s J_{2} 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 q_{3} = 0 and P_{3} as in Eq. (31). These are used also for the dashed blue curve, where Eq. (25) and the secondorder terms of Eq. (22) are also taken into account. 

Open with DEXTER 
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(29)
where the units of measure are determined by Io’s semimajor 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_{1}, q_{2}, q_{3}, q_{4}} = {0, π, π, π} by employing a rootfinding 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,(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_{3} argument. It is therefore sufficient to vary the action P_{3} and therefore the eccentricity of Ganymede (if a fixed semimajor 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(31)
where the unit of length is assumed normalized with respect to Io’s semimajor 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 semimajor axis). Furthermore, to avoid being captured in the librational regime, we need to chose the initial value of q_{3} properly. We here selected q_{3} = 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_{3} and q_{3}, 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 firstorder 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 q_{3} = 0 and P_{3} as in Eq. (31). Bottom: corresponding filtering at 1000 days. 

Open with DEXTER 
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(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., semimajor 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 quasiresonant mean motion relations of the pairs IoEuropa and EuropaGanymede. 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_{2} (and, of course, including the pointmass 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 actionangle 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 nonoblateness 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 2013056RO.1, and we acknowledge the GNFM/INdAM. We are also grateful to S. FerrazMello 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
 Broer, H. W., & Hanßmann, H. 2016, Indagationes Mathematicae, 27, 1305 [CrossRef] [Google Scholar]
 Brown, B. 1977, Celest. Mech. Dyn. Astron., 16, 229 [CrossRef] [Google Scholar]
 Bucciarelli, S., Ceccaroni, M., Celletti, A., & Pucacco, G. 2016, Annali di Matematica Pura e Applicata, 195, 489 [CrossRef] [MathSciNet] [Google Scholar]
 Celletti, A., & Gales, C. 2018, SIAM J. Appl. Dyn. Sys., 17, 203 [CrossRef] [Google Scholar]
 Celletti, A., Paita, F., & Pucacco, G. 2018a, Celest. Mech. Dyn. Astron., 130, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Celletti, A., Paita, F., & Pucacco, G. 2018b, The Dynamics of LaplaceLike Resonances, in press [Google Scholar]
 de Sitter, W. 1931, MNRAS, 91, 706 [NASA ADS] [CrossRef] [Google Scholar]
 FerrazMello, S. 1975, Celest. Mech. Dyn. Astron., 12, 27 [CrossRef] [Google Scholar]
 FerrazMello, S. 1979, Mathematical and Dynamical Astronomy Series (Sāo Paulo: Universidad, Instituto Astronômico e Geofísico) [Google Scholar]
 Giorgilli, A., Locatelli, U., & Sansottera, M. 2017, Reg. Chao. Dyn., 22, 54 [CrossRef] [Google Scholar]
 Henrard, J. 1984, Celest. Mech. Dyn. Astron., 34, 255 [CrossRef] [Google Scholar]
 Kamel, A. 1970, Celest. Mech. Dyn. Astron., 3, 90 [CrossRef] [Google Scholar]
 Kaula, W. 2000, Theory of Satellite Geodesy, Dover Books on Earth Sciences (Mineola, New York: Dover Publications, Inc.) [Google Scholar]
 Kosmodamianskii, G. 2009, Sol. Syst. Res., 43, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Lainey, V. 2002, PhD Thesis, Observatoire de Paris, Paris, France [Google Scholar]
 Lainey, V., Duriez, L., & Vienne, A. 2004, A&A, 420, 1171 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laplace, P.S. 1805, Traité de Mécanique Céleste (Paris, France: Chez Courcier), IV [Google Scholar]
 Lari, G. 2018, Celest. Mech. Dyn. Astron., 130, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Petit, A. 2017, A&A, 605, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lieske, J. 1977, A&A, 56, 333 [NASA ADS] [Google Scholar]
 Lieske, J. 1997, A&AS, 129, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malhotra, R. 1991, Icarus, 94, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Marsden, B. 1966, PhD Thesis, Yale University, USA [Google Scholar]
 Morrison, J. 1966, Prog. Astronaut. Rocketry, 17, 117 [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Musotto, S., Varadi, F., Moore, W., & Schubert, G. 2002, Icarus, 159, 500 [NASA ADS] [CrossRef] [Google Scholar]
 Pucacco, G., & Rosquist, K. 2017, J. Geom. Phys., 115, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Sagnier, J. 1975, Celest. Mech. Dyn. Astron., 12, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Sampson, R. 1921, Mem. R. Astron. Soc., 63, 1 [NASA ADS] [Google Scholar]
 Showman, A., & Malhotra, R. 1997, Icarus, 127, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Yoder, C., & Peale, S. 1981, Icarus, 47, 1 [NASA ADS] [CrossRef] [Google Scholar]
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(A.1)
and Callisto’s contribution, Eq. (25), can be split into three separate terms (one each for Io, Europa, and Ganymede) as(A.2)
All Tables
SPICE Cartesian elements for the Galilean moons in the Jovicentric equatorial frame.
All Figures
Fig. 1. Jovicentric snapshots of the Jovian system, with T denoting Io’s mean orbital period. The moons are colorcoded (Io in red, Europa in green, and Ganymede in blue). 

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

Open with DEXTER  
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 semimajor axis, and by Jupiter’s mass. 

Open with DEXTER  
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 lowpass filter at 1000 days (tails are a byproduct of the procedure). Ephemerides from the NASA Spice toolkit. 

Open with DEXTER  
In the text 
Fig. 5. From top to bottom: integration history of the Laplace argument obtained without oblateness terms, including Jupiter’s J_{2} 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. 

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

Open with DEXTER  
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 colorcoded (Io in red, Europa in green, and Ganymede in blue). 

Open with DEXTER  
In the text 
Fig. 8. History of the absolute errors on the Jovian moons’ semimajor axes as computed in the Cartesian and Hamiltonian formulations. The moons are colorcoded (Io in red, Europa in green, and Ganymede in blue). 

Open with DEXTER  
In the text 
Fig. 9. From top to bottom: Cartesian and Hamiltonian histories for the angle q_{3}, corresponding error on the evolution, and similar error between the Cartesian formulation and the ephemerides set. 

Open with DEXTER  
In the text 
Fig. 10. Phase plane for the normal form Hamiltonian obtained by restricting Eq. (22) to firstorder terms in the eccentricity and retaining the full secular effect of Jupiter’s J_{2} 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 q_{3} = 0 and P_{3} as in Eq. (31). These are used also for the dashed blue curve, where Eq. (25) and the secondorder terms of Eq. (22) are also taken into account. 

Open with DEXTER  
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 q_{3} = 0 and P_{3} as in Eq. (31). Bottom: corresponding filtering at 1000 days. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.