A&A 420, 1171-1183 (2004)
DOI: 10.1051/0004-6361:20034565
V. Lainey^{1} - L. Duriez^{1,2} - A. Vienne^{1,2}
1 - Institut de Mécanique Céleste et de Calcul
des Éphémérides - Observatoire de Paris, UMR 8028 du CNRS,
77 avenue Denfert-Rochereau,
75014 Paris, France
2 - Université Lille1 IMCCE/LAL,
1 Impasse de l'Observatoire, 59000 Lille, France
Received 23 October 2003 / Accepted 15 March 2004
Abstract
We present a new theory of the four Galilean satellites Io, Europa, Ganymede and Callisto. This theory aims to deliver highly accurate ephemerides able to represent the Galilean satellites' motion over several centuries. It is based on the numerical integration of elaborated equations of motion. This first paper describes and tests many usually neglected perturbations. We are then able to retain some of them in the dynamical model for the Galilean system. A numerical method developed to adjust the model to the observations is given. We used a general formalism so it can be extended to systems other than the Galilean one. As an example of this method, we compare our model to the current E5 ephemerides of the four Galileans.
Key words: celestial mechanics - planets and satellites: individual: Jupiter - ephemerides
The Galilean satellites have been of special interest recently because of the Galileo spacecraft arrival. Besides the results involving the internal structure of the satellites, a lot of dynamical parameters are now well known. Hence the ephemerides of this system can be improved as a better modeling is now possible.
The ephemerides generally used for the Galilean system are based on Sampson-Lieske theory. The initial model was developed at the beginning of the last century by Sampson 1921, and improved by Lieske for the Voyager spacecraft needs (Lieske 1977). The last adjustment of this theory to the observations was performed in 1998 and is called E5 (Lieske 1998). Kaas et al. (1999) mentioned the possibility that E5 would need some more revision by comparison with the accuracy obtained with mutual event observations (in particular for Io). Indeed, the accuracy of such observations are estimated to be few tens of kilometers. Moreover, Mallama et al. (2000) also found important residuals using the eclipse observations of the Galilean satellites by Jupiter.
Highly accurate ephemerides are necessary to detect secular accelerations induced by tidal effects, after comparison with observations. Such accelerations can then be directly linked to the dissipation inside the satellites, and so to their internal structure.
In this first paper we give a numerical method able to provide ephemerides precise enough with respect to mutual event accuracy. We examine the equations of motion describing the Galilean system using an accurate modeling. We will introduce the triaxiality of the satellites, some additional Jovian oblateness forces and some other perturbations. In Sect. 3 we will test, by use of numerical integration, the perturbations we introduced to determine which perturbation should finally be taken into account. In Sect. 4 we will give the equations based on Moulton's method to perform an adjustment of the numerical model to the observations. In the final section, we will perform an adjustment to the E5 ephemerides.
The adjustment of our numerical model directly to the observations is still ongoing and will be presented in a second paper.
We develop the equations of motion in a planetocentric reference frame with fixed axes (that may be nonequatorial). Even though in such reference frame the expression of the disturbing potential induced by the central body's oblateness becomes rather complex, an equatorial reference frame would generate additional inertial forces in the system (equatorial precession and Jupiter rotation). Moreover, other complications arise in the formulation of the equations related to the oblateness or the triaxiality of the satellites. It would then be necessary to write the disturbing potential of the triaxialily of the satellites in a fixed reference frame at the first time (for instance J2000), and then in an equatorial reference frame of the central body, requirinq much more complicated equations.
In the following, the ponctual masses will be denoted (P_{i}, m_{i}). In particular, the index i will be zero in the case of the central body. The distance between two bodies P_{i} and P_{j} will be noted r_{ij}. When r will have only one index (let us say i) instead of two, it will denote the distance of body P_{i} to the central body P_{0}. We will denote over the indixes 0,i,j,k... a bar or a hat when bodies P_{0},P_{i},P_{j},P_{k},... related to each one of these indices will be respectively modeled as ponctual or oblated (see for example Eq. (8)).
The symmetry of the equations relative to the Cartesian coordinates P_{i}=(x_{i},y_{i}, z_{i}), frequently requires to indicate by the variable , any of the three Cartesian coordinates. In Sect. 4, we will even need to introduce this notation twice into the same equation, and we will note this second variable by .
Let bodies P_{k} of mass m_{k} assumed ponctual and a central oblated body P_{0} of mass m_{0} attracting according to the Newtonian gravitation laws. We consider an orthonormal reference frame denoted (P_{0}, x, y, z) centered on P_{0} with fixed axes.
We denote
,
,
.
Using the vectorial equality
where O is the origin of an arbitrary Galilean frame, we can easily write the equations of motion in the planetocentric frame
(P_{0},x,y,z). We have the following classical differential system
For the forces induced by the central body,
with
and
.
This last potential represents the action of P_{0}'s triaxiality on P_{k} and is written, using the equatorial spherical coordinates
of P_{k} in the reference frame
related to planet P_{0} (see Fig. 1), in the following form
Figure 1: Equatorial frame of the central body P_{0} precessing and rotating denoted . | |
Open with DEXTER |
Hence, the differential equation related to the motion of a body P_{i} can be rewritten as
For an explicit formulation of Eq. (5), we must write the gradient of the expression (2) relative to the Cartesian coordinates (x_{k}, y_{k}, z_{k}) in the planetocentric reference frame with fixed axes (P_{0}, x, y, z).
We have after calculation^{} for the first term
In this subsection, we will now present the terms to be introduced to take into account the triaxialily of all the bodies.
For convenience of notation, we consider the same reference frame centered on the body P_{0} with fixed axes (
P_{0}, x, y, z), although the results to follow remain valid for a reference frame centered on any other body. We still have the same differential system (1) but the force on P_{k} exerted by P_{l} is now
Taking into account Eq. (8), three terms should be finally added to
Eq. (5) to take into account the triaxiality of the bodies, which are
The explicit expression of accelerations of Eq. (9) can be done using the expressions found in the previous section. Hence, one can rewrite these three terms in the form
Figure 2: Equatorial reference frame of a satellite P_{l} related to the reference frame of fixed axes ( P_{0},x,y,z). | |
Open with DEXTER |
To obtain the expression of in the initial reference frame centered on body P_{0}, it is enough to substitute formally in the expressions of and the coordinates ( x_{k}, y_{k}, z_{k}) by ( x_{k}-x_{l}, y_{k}-y_{l}, z_{k}-z_{l}) and the angles by . In the case of and the change of coordinates simply amounts to replacing ( x_{k}, y_{k}, z_{k}) by ( -x_{k}, -y_{k}, -z_{k}). Moreover, the values of the equatorial radii and oblateness coefficients should be replaced by those relating to each satellite.
Finally, we obtain for the differential equation associated with a body P_{i} the expression
Other simplifications are also possible such as not considering the indirect perturbations (satellite-planet) for the coefficients J_{3}^{(k)}, J_{4}^{(k)}, J_{6}^{(k)} and even J_{2}^{(k)}. These perturbations are additional oblateness forces, but they are multiplied by the term with . These terms are thus even smaller than the additional oblateness forces related to the oblateness of the central body. However, we will preserve all these terms since they are necessary in the conservation of the first integrals of the system (see the Sect. 2.5).
Assuming that the satellites have their spin axis parallel with that of the central body (which is the case for most satellites close to their planet), a very simple method exists for taking into account the polar oblateness of the satellites (coefficients J_{n}^{(k)}), without using the previous equations. Indeed, let us recall that the disturbing potential related to the polar oblateness of the central body and of its satellites is
(15) |
In addition, we can deduce that during the adjustment to the observations, the coefficients J_{n} will be correlated to the values of the coefficients J_{n}^{(k)}, if the polar oblateness of the satellites is not taken into account. Thus, if one wishes to use the values J_{n}^{(k)} given by the space probes, it is necessary to carry out major corrections.
With the assumption of parallel spin axes, it is very easy to take into account in a numerical integration the ponctual satellite-oblate planet perturbation induced by a coefficient J_{2}^{(k)}. Indeed, one just needs to replace the value of the J_{2} coefficient by the new value during the integration of the satellite considered as oblate. For an analytical theory, it will be enough to modify the values of the coefficient J_{2} in the expression of the analytical series.
The coefficients c_{22}^{(k)} and s_{22}^{(k)} characterize the equatorial ellipticity of the satellites. Generally the reference frame associated with them has for the first direction the principal axes of inertia of the satellite. The equatorial deformation has a symmetry (until a certain precision) relative to one of the inertial axes. Hence, the coefficients s_{22}^{(k)} are equal to zero. And so, for the Galilean satellites, only the coefficients c_{22}^{(k)} will have to be considered.
The system (1) is incomplete if the Euler angles and , are changing with time (precession and rotation of bodies). Deriving the Euler equations that describe the orientation of the satellites is beyond the perspective of this paper (a non-rigid model should be used and initial conditions been required).
For Jupiter, it is convenient to take the results presented in Seidelmann et al. (2002), which explicitly give as a function of time the Jovian Euler angles.
The evolution of the Jovian rotation pole and the origin meridian line relative to the J2000 equinox at the time J2000.0 is given by the equalities (expressed in degrees)
The system (16) is thus directly usable and remains connected to the variables of our system by the relations , and where W_{0} is the longitude of the central meridian line refered to the axis (see Fig. 1).
We present the corrections that should be added to introduce relativistic effects. This section is strongly inspired by the work of Ferraz-Mello (1967). We start from the equation of the
-body problem (ponctual masses), developed in a Galilean reference frame using a relativistic formulation, given by de Sitter
= | |||
(17) |
In the case of the Jovian system, we can consider the relativistic terms depending only on the Jovian mass. In particular we neglect the relativistic effects induced by the Solar and satellites masses.
We thus obtain for a satellite P_{i} the equation
= | |||
(20) |
(21) |
We will use the conservation of energy of the system both to control the efficiency of our numerical integrator and to validate the computation of the equations of motion. The integral will not be preserved in a system for which some forces will be taken in an incomplete way.
Preserving the notations of the preceding sections and denoting M the sum of the masses of the system, we have for the energy E the following expression
(22) |
To obtain a complete explicit expression, it remains to express the gravitational potential of the system which will change with the perturbations included in the modeling. Let us consider the
-triaxial bodies problem relative to Eq. (13); we have after neglecting the terms
and
of Eqs. (11) and (12) the equality
(23) |
We present the method we used to test the influence of each perturbation, as well as the numerical tools that we have developed.
Figure 3: Variations of the mean longitudes multiplied by semi-major axes (in kilometers on the top) and of the inclinations multiplied by the semi-major axes (in kilometers on the bottom) for Io, Europa, Ganymede and Callisto under the effect of the Jovian precession. | |
Open with DEXTER |
Figure 4: Variations on positions (in kilometers) under the effect of the coefficients c_{22}^{(k)} with an exact spin-orbit resonance model for the satellites' rotation. | |
Open with DEXTER |
Figure 5: Variations on positions (kilometers) after substracting a straight line, under the action of the coefficients c_{22}^{(k)}. | |
Open with DEXTER |
We consider a first solution used as a reference over one century obtained by numerical integration of the four Galilean satellites and an oblated Jupiter using the coefficients J_{2}, J_{4}, J_{6}.
Then with the same initial conditions, we integrate a perturbed system by adding one of the perturbations to be tested to the main system. The differences between these two integrations thus provides the influence, in the sense of differences in the positions or variation of the elliptical elements of the studied perturbation. In practice we used a period of one century, which is very close to the time span of the astrometric observations. Thus, Fig. 3 to Fig. 5 present these differences (in km) as a function of time (expressed in years), over one century. The variations in longitude have been multiplied by the semi-major axes of reference to obtain differences in kilometers.
We should note that this method is not completely satisfactory as even small variation in mean motion could imply large differences in positions. We will discuss this at the end of the section.
The main system could include only the simplified equations which are generally used, namely
To test each perturbation, we programmed the corresponding equations of motion. The general formulation of the equations of Sect. 2.2 enabled us to apply our program to an unspecified dynamical system.
The integration subroutine implemented in the code is RA15 (Everhart's integrator, see Everhart 1985) based on the method of the Gauss-Radau polynomias. This routine of fifteenth order has the advantage of being both accurate and fast.
Finally, the initial conditions were taken from the theory G5 (Arlot 1982), available on the web server^{} of the IMCCE at the epoch of December 25, 1982, 0H00. Those are expressed in a J2000 reference frame in terrestrial time. The values of the masses and coefficients of triaxiality of the bodies were taken from Anderson et al. (1996) and Schubert (1994).
We used two complementary methods. The first method is based on the differences of positions and velocities given by the integrator after integrating back and forth over one century. This method enabled us to obtain an optimal step size equal to 0.08 day, giving errors of about a few tens of meters.
The second method consists of studying the variation of the integral of energy. As we already said, this method is not possible in all cases. The relative variation of this first integral was about 10^{-14} in all the cases we tested. The variation of energy on a back and forth integration then enables us to get metric differences (at least to obtain an order of magnitude from it). The result obtained, which appeared in agreement with the results of the first method, is a few tens of meters for one century of integration.
Table 1: Summary of the influence of each perturbation tested in decreasing order of magnitude. The differences refer to the most important ones for the four Galilean satellites. The name in parenthesis refers to the satellite the most perturbed.
The numerical integrations used in this section do not take into account the Sun, nor the other planets, in particular Saturn. Indeed, as we did not want to integrate the motion of all planets of the Solar system, the position of the Sun has been computed by an external program. Then, the control of the first integral is not possible anymore.
However the continuation of our work required us to include the solar perturbation in our numerical program. Thus, the Solar perturbation and that of Saturn were included using the numerical ephemerides DE406 of the^{} JPL. During this implementation in our program, we in particular wondered about the possible need to take for Jupiter mass the value resulting from DE406 and not that of Campbell & Synnott (1985). However the planetary theories represent the barycentre of the planetary systems (except for the Earth-Moon system). Admittedly it would be possible to cut off from the value of the Jupiter mass in DE406, the value of the masses of the Galilean satellites from the probes but this effect would remain artificial (it would still miss the mass of the non Galilean satellites). We thus did not use the mass of DE406 in our work.
A summary presenting all the perturbations tested with their differences on positions can be found in Table 1.
The study of some unusual perturbations was partially done in Lainey et al. (2001). We present here some graphs related to the perturbations which were not tested (or shown) then. We will also try to understand the real importance of the differences in positions obtained in Sect. 5.
The Jovian equatorial reference frame is not fixed in space. In particular, the angles and I for which the coefficients J_{n} are linked change with time. As this effect is partly a geometrical one and introduces only rotations, we can expect that it is the angular variables of the satellites that will change the most. An integration over one century has been done. Variations in mean longitudes and inclinations are shown in Fig. 3. One can see that Ganymede and Callisto are the most influenced by this perturbation. As the Jovian precession is basically a kinematic effect, the most influenced satellites are also the most distant. To opposite most perturbations that we tested, the mean longitude does not absorb the main effect of the perturbation, which reach^{} 110 kilometers. Hence the variation of the inclinations multiplied by the semi-major axes reaches 80 kilometers for Callisto. An analysis of the other elements reveals that the variations in semi-majors axes and eccentricity are negligible.
A, B and C in the equations of motion (see the page ). In Lainey et al. (2001), the study of J_{2}^{(k)} perturbation was estimated to 5000 kilometers over one century. We will here present the contribution of c_{22}^{(k)} coefficients. A good estimate of these coefficients for the Galilean satellites is known today (Schubert et al. 1994). They are estimated at approximately three tenths of the value of the respective J_{2}^{(k)} coefficients. To study their influence, it is necessary to introduce the rotation of the satellites. We modeled here only the exact spin-orbit resonance case ( ), the treatment of the satellites' rotation being deferred to a future work. The results for this perturbation are given in the Fig. 4. This perturbation is the most significant that we tested (approximately 9000 kilometers for Io), and affects essentially the longitudes. The oscillations on the straight lines are induced mainly by the differences of amplitude and frequency on the Laplacian libration. Indeed, we recall that Io, Europa and Ganymede are involved in a resonance of argument , where l denotes the mean longitude.
The variations in amplitude and frequency of the Laplacian libration can be seen directly as oscillations in Fig. 4 for the three resonant satellites. Figure 5 shows these variations after subtraction of the secular drifts. To determine the modifications in amplitude and frequency, we carried out a Fourier analysis on the two integrations concerned. The differences obtained over the amplitude and the period of the Laplacian libration are respectively only 0.01 degree and 0.003 day.
The majority of the perturbations induce mainly a linear term in time which is added to the mean longitudes. Taking into account such perturbations will rather poorly improve the development of the ephemerides, as it can be transfered in the determination of the initial conditions, and in particular that of the mean motion. However that does not remain completly without consequences since internal coherence of the theory is somewhat affected (see Sect. 5).
Moreover, if the ephemerides are only slightly improved, the physical model is considerably more complex. The rotation of Jupiter and spin-orbit resonances can be reserved for later work.
Table 1 presents the highest amplitude found for each perturbation we tested. We recall that to test the conservation of energy, the integrations reported in Table 1 did not take into account the Sun.
During the adjustment method, we will have to introduce the equations of variations, which are much more complex than the equations of motion. Thus we tried to limit modification to the most significant perturbations.
From Table 1, we decided to preserve in our dynamical model the additional oblateness forces and the triaxiality of the satellites using the coefficients J_{2}^{(k)} and c_{22}^{(k)}. These perturbations are the most significant by more than an order of magnitude than the others.
The perturbation by Amalthea was not retained, because we did not get a valid solution for this satellite over one century (period covered by the observations for the Galilean satellites). For the same reason, we did not included the Jovian precession.
In this section we will distinguish the physical parameters (such as the masses, the coefficients of oblateness, etc.) which we will denote by the vector , and the initial condition (which are the position-velocities of the bodies at an instant of reference t_{0}) of a theory. We denote these constant , or in short . We have to solve the following problem: How to determine optimal values for the parameters and initial conditions, so that the differences between the positions of the satellites delivered by our model and those resulting from the observations (O-C) are minimal in a least-squares sense.
Let us consider a set of observations
providing for discrete values of time t_{k} the observed positions
.
The theory (here the numerical integration), gives for these same values of time the computed positions
,
as solution of a differential system. Denoting
the dynamic flow of this system and
the projection in
of
such that
,
one has the equality
From the Taylor formula, we have
(26) |
(28) |
= | |||
= | (29) |
(30) |
(31) |
The differential system (32) requires us to write the vectors
and
as functions of the variables
.
First, we take for the expression of
the one given in Eq. (5) when modeling the satellites as ponctual
= | |||
= | |||
(33) |
Equation (32) can thus be rewritten, denoting
,
in the form
(34) |
(36) |
Table 2: Initial conditions found after fitting our numerical model to the E5 ephemerides. The positions are given in UA and the velocities in UA.j^{-1}.
Table 3: Parameters found after fitting our numerical model to the E5 ephemerides. The masses are given in Solar mass and the angles are in degrees.
Now we introduce the triaxiality of the satellites into the equations of variation, after neglecting the terms of the form
(term
of Eq. (9)). We replace in Eq. (35) the expression of
by the following one
= | |||
(38) |
In practice, Eq. (35) have been integrated simultaneously with the equations of motion.
We have adjusted our model to the Sampson-Lieske theory using the method presented in the preceding section.
The theory that we used to adjust our model is the last adjustment of Lieske's theory to the observations, called the E5 theory (Lieske 1998). To reduce the modeling differences, we did not include the additional oblateness forces and the satellites' triaxiality. Using one version of the theory instead of another (G5 for instance) does not change the results, as the theory remains the same and only the parameters change (Lieske 1977). Hence, just the internal differences with the two models will appear.
The adjustment was carried out in Cartesian equatorial coordinates J2000 using a sample of 3654 points per coordinate and with a step size of 10 days. The period covered by the adjustment is then approximately a century. This time span is close to the one we plan to use in the fit to the observations. To reduce the numerical errors, we integrated over fifty years back and forth starting at the epoch 01/01/1950 at 0H00.
The constants adjusted in our model were as follows:
The results are shown in Fig. 6 using the variables . We note that the discrepancies are about several hundred kilometers for each satellite, and of almost 800 kilometers for Europa. The root-mean-square (rms) position differences between our model and E5 are of 53.96 km, 127.50 km, 81.14 km and 91.19 km for Io, Europa, Ganymede and Callisto. These differences appear especially in the presence of long periods on the variable and in a less significant way the variable z_{i}. To understand such long period differences, we examine the variables used by Sampson.
Generally, the solution of a theory in celestial mechanics is given as a Fourier series with a polynomial amplitude (Poisson series), in which one can separate the short period terms and the long period terms. The long period terms are generally those having the strongest amplitudes (because of small dividers after integration), while the short period terms have small variable amplitudes, as a function of the long periods. This general characteristic with the theory expressed in elliptical elements is slightly different using the variables of Sampson because of the use of the coordinate ,
or equally to the coordinates z_{i}. Indeed, the nature of this coordinate (rather inappropriate to describe a disturbed Keplerian motion) leads to high oscillations in the amplitudes of z_{i} for the frequencies higher than the period of revolution of the satellite. In Lieske (1977), there are no long-period terms in the final development of the variables z_{i} at all.
Figure 6: Residuals after an adjustment over 100 years of our model to E5 ephemeride, for Io, Europa, Ganymede and Callisto ( from top to bottom) on the variables ( from left to right). The differences (in kilometers) are a function of time (in years) over one century. | |
Open with DEXTER |
Table 4: Variations introduced in the initial conditions of the satellites to minimize the influence of the additional oblateness forces and the triaxiality of the satellites. The values are given in kilometers for the positions and in kilometers per day for the velocities.
Table 5: Variations introduced in the parameters to minimize the influence of the additional oblateness forces and the triaxiality of the satellites.
One can wonder why the variations we found are more significant for the variable than the variable z_{i}. During our adjustment we had previously kept constant the values of the angles . The corresponding results then gave a stronger variation on z_{i} than that on . However, we ended up also adjusting the angles . The variations in z_{i} were then reflected on by a geometrical effect, thus giving the graphs of Fig. 6.
Another reason contributing to give these broad variations can be the use of different solar theories. To a lesser extent the influence of the Jovian precession (not taken into account in our model) could also slightly contribute to the residuals.
If the residuals shown in Fig. 6 represent the internal error of the Samspon-Lieske theory, one should expect to see such differences in the comparison to observations. Hence, in the paper of Mallama et al. 2000, some rms differences in distance of 62 km, 267 km, 142 km and 146 km for respectivly Io, Europa, Ganymede and Callisto were found. These residuals appear proportional to ours, but are higher. This is not surprising as the observations induce experimental errors, and as our residuals were obtained over one century instead of ten years for the observations used in Mallama et al. (2000).
On the other hand, the discrepancies for Io are the lowest ones. That differs from the interpretation by Kaas et al. (1999) that some fairly large (O-C)s in the relative positions of Io and Europa from mutual event observations in 1991 are likely due mainly to errors in Io's rather than in Europa's E5 positions, as found by us. Of course, our residuals are based on a spatial representation of E5. Hence some (O-C)s on right ascention and declination will be smaller, but that should not change the relative discrepancies. Nevertheless mutual events cover just one year and some local scattering in our graphs (see Fig. 6) could explain the residuals of the 1999 PHEMU campaign.
In Sect. 3, we had decided to preserve any small perturbations in our model without being able to define precisely the real influence of these latter. Indeed, although the difference between taking into account these perturbations or not was estimated, we noted that a small change in the constants of the theory would be enough to make these differences disappear (at least in a part).
Here, we determine precisely which constants are modified and how so as to minimize the lack of these perturbations in a model.
To complete this work, we consider again the adjustment method already used in the preceding section. We adjusted the main system (see Eq. (24)) with our final model (which takes into account the additional oblateness forces, and the triaxiality of the satellites).
The adjusted constants were the same as in the last subsection. Moreover, four iterations by the least-squares method were still carried out to optimize our adjustment. The differences obtained are ploted using the same variables , to understand the possible influence of the perturbations that we took into account and are neglected in the model of Sampson-Lieske.
We also chose a sample of 3654 points per coordinate with a step size of 10 days. Once again, to limit the numerical errors, we integrated over fifty years back and forth.
The differences that one cannot correct reach only a few kilometers (graphs not shown here, for more details see Lainey 2003). However the absorption of a perturbation is carried out at the price of a shift between the real physical quantities of the system and those of the model. Table 4 gives the variations induced to the initial conditions of the model to obtain our adjustment, and Table 5 gives the variations for the physical parameters.
It appears that if the variations in the initial conditions remain very weak, those in the parameters and especially the oblateness coefficients J_{2} and J_{4} are clearly greater than the error bars. As an indication, the error bars for these two coefficients (Campbell & Synnott 1985) are respectively of only 10^{-6} and whereas we obtain variations of and . This result is coherent with what was stated in Sect. 2.2. First of all the additional oblateness forces are directly related to the Jovian oblateness. Moreover, we had mentioned the possibility of absorbing the coefficients J_{2}^{(k)} of the satellites in the oblateness of the central body.
Finally, it may be surprising that the variations are more significant for J_{4} than for J_{2}. This effect comes directly from the correlation between these two coefficients. Actually, to minimize the variations in J_{4}, we have constrained this coefficient to remain constant during the first two adjustments with the least-squares procedure, and left it free in the two last.
We have determined which perturbations should be included to have an accurate modeling of the Galilean system. Hence, the additional oblateness forces and the satellites' triaxiality appeared non negligible. We gave explicitly the equations which should be used to adjust the numercial solution to the observations. We used this method to adjuste our model on the E5 ephemerides. The maximum residuals reached 250 km, 700 km, 500 km and 550 km for respectively Io, Europa, Ganymede and Callisto over one century. These residuals may be explained by the variables used in the Sampson-Lieske theory. Moreover they appear consistent with the results found by Mallama et al. (2000) as the authors had only a time span of comparison of ten years. As also discussed in Mallama et al. (2000) we did not find high residuals for Io as was hypothesized in Kaas et al. (1999) to explain mutual event (O-C)s, even if some local discrepancy induced by the internal error of E5 ephemerides could eventually explain their result.We will present an adjustment to the observations in a future paper. In particular the new ephemerides could be used to reduce the mutual phenomena, and should be accurate enough to detect the secular acceleration induced by tidal effects.
Finally, it is important to note that our formulation was done in a general way and can be applied to other satellite (planetary) systems.