A&A 420, 1171-1183 (2004)
DOI: 10.1051/0004-6361:20034565

New accurate ephemerides for the Galilean satellites of Jupiter

I. Numerical integration of elaborated equations of motion

V. Lainey1 - L. Duriez1,2 - A. Vienne1,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

1 Introduction

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.

  
2 Equations of motion

2.1 Reference frame and notations

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  (Pi, mi). In particular, the index i will be zero in the case of the central body. The distance between two bodies Pi and Pj will be noted rij. When r will have only one index (let us say i) instead of two, it will denote the distance of body Pi to the central body P0. We will denote over the indixes 0,i,j,k... a bar or a hat when bodies P0,Pi,Pj,Pk,... 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 Pi=(xi,yi, zi), frequently requires to indicate by the variable $\gamma_i$, 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 $\Gamma_i$.

  
2.2 Equations of the ${\cal N}$-body problem with ponctual masses in a planetocentric frame with fixed axes

Let us start by modeling the Galilean system with an ${\cal N}+1$-body problem for which only the central body (Jupiter) is supposed to be oblate (or triaxial in a more general sense). The number ${\cal N}$ will refer to various values according to the number of bodies introduced into the model, in addition to the four Galilean satellites (Sun, planets, inner satellites...).

  
2.2.1 Expression of the additional oblateness forces

Let ${\cal N}$ bodies Pk of mass mk assumed ponctual and a central oblated body P0 of mass m0 attracting according to the Newtonian gravitation laws. We consider an orthonormal reference frame denoted (P0, x, y, z) centered on P0 with fixed axes.

We denote $\vec{r}_i=\overrightarrow{P_0P_i}$, $r_i=\vert\vec{r}_i\vert$, $r_{ij}=\vert\vec{r}_j-\vec{r}_i\vert$. Using the vectorial equality $\vec{\ddot{r}}_i=\ddot{\overrightarrow{OP_i}}-\ddot{\overrightarrow{OP_0}}$ where O is the origin of an arbitrary Galilean frame, we can easily write the equations of motion in the planetocentric frame  (P0,x,y,z). We have the following classical differential system

 
$\displaystyle \left\{
\begin{array}{lcl}
\ \vec{\ddot{r}}_1&=&\displaystyle{\fr...
...le{\frac{\vec{F}_{\cal N}}{m_{\cal N}}-\frac{\vec{F}_0}{m_0}}\end{array}\right.$     (1)

where $\vec{F}_k$ indicates the whole external forces exerted on the body Pk. $\vec{F}_{ki}$ is the exerted force on Pk by Pi. Introducing the gravitationnal potentials, we have $\vec{F}_{ki}=Gm_km_i{\bf\nabla}_kU_{ki}$, where $U_{ki}=U_{\bar{k}\bar{\imath}}=\frac{1}{r_{ki}}$.

For the forces induced by the central body, $\vec{F}_{k0}=\vec{F}_{\bar{k}\bar{0}}+
\vec{F}_{\bar{k}\hat{0}}$ with $\vec{F}_{\bar{k}\bar{0}}=Gm_km_0{\bf\nabla}_kU_{\bar{k}\bar{0}}$ and $\vec{F}_{\bar{k}\hat{0}}=Gm_km_0{\bf\nabla}_kU_{\bar{k}\hat{0}}$. This last potential represents the action of P0's triaxiality on Pk and is written, using the equatorial spherical coordinates $(r_k, \phi_k, \lambda_k)$ of Pk in the reference frame $(P_0,x^\prime ,y^\prime ,z^\prime )$ related to planet P0 (see Fig. 1), in the following form

 
$\displaystyle U_{\bar{k}\hat{0}}=U_{\bar{k}\hat{0}}^{(1)}+U_{\bar{k}\hat{0}}^{(2)}$     (2)

where
 
$\displaystyle U_{\bar{k}\hat{0}}^{(1)}=\sum_{n=2}^\infty - \frac{(E_r)^n}{r_k^{n+1}}J_nP_n(\sin \phi_k)$     (3)

and
 
$\displaystyle U_{\bar{k}\hat{0}}^{(2)} = \sum_{n=2}^\infty \frac{(E_r)^n}{r_k^{...
..._n^{(p)}(\sin \phi_k)
\left[c_{np}\cos p\lambda_k +s_{np}\sin p\lambda_k\right]$     (4)

with Er denoting the equatorial radius of P0. Note that only  $U_{\bar{k}\hat{0}}^{(1)}$ is a function of the latitude $\phi_k$.
  \begin{figure}
\par\includegraphics[width=6.9cm,clip]{0565fig1.eps}
\end{figure} Figure 1: Equatorial frame of the central body P0 precessing and rotating denoted $(P_0,x^\prime ,y^\prime ,z^\prime )$.
Open with DEXTER

Hence, the differential equation related to the motion of a body Pi can be rewritten as

 
$\displaystyle \vec{\ddot{r}}_i=-\frac{G(m_0+m_i)\vec{r}_i}{r_i^3}
+\sum_{j=1,j\...
...h}\hat{0}}+\sum_{j=1,j\neq i}^{\cal N}Gm_j{\bf\nabla}_jU_{\bar{\jmath}\hat{0}}.$     (5)

The last term of Eq. (5) as well as the mass mi in factor before the last term are generally forgotten in the equations of motion. Indeed, they represent terms of second order (order of the product mkJ2 at the most). Physically, these terms correspond to indirect forces resulting from the oblateness of the central body. We will preserve these terms in the continuation under the name of additional oblateness forces.

  
2.2.2 Expression of the disturbing potential $\textit{{U}}_{{\bar\textit{k}}\hat{0}}$ in a reference frame with fixed axes

For an explicit formulation of Eq. (5), we must write the gradient of the expression (2) relative to the Cartesian coordinates  (xk, yk, zk) in the planetocentric reference frame with fixed axes  (P0, x, y, z).

We have after calculation[*] for the first term

 
$\displaystyle \frac{\partial U_{\bar{k}\hat{0}}^{(1)}}{\partial\gamma_k}=
-\sum...
..._k)}{{\rm d}\sin \phi_k}
-\frac{(n+1)\gamma_kP_n(\sin \phi_k)}{r_k}\right]\cdot$     (6)

And for the second term, we have
 
                            $\displaystyle \frac{\partial U_{\bar{k}\hat{0}}^{(2)}}{\partial\gamma_k}$ = $\displaystyle \sum_{n=2}^\infty
\sum_{p=1}^n\frac{(E_r)^n}{r_k^{n+1}}\bigg\{\fr...
...al z_k^\prime}
{\partial\gamma_k}-\frac{\gamma_k\sin \phi_k}{r_k}\right)\right.$  
    $\displaystyle \times \left.\left.\frac{{\rm d}P_n^{(p)}(\sin \phi_k)}{{\rm d}\s...
...hi_k}\right.-\frac{(n+1)\gamma_kP_n^{(p)}(\sin \phi_k)}{r_k}_{\frac{}{}}\right]$  
    $\displaystyle \times \left[{c_{np}\cos p\lambda_k+s_{np}\sin p\lambda_k}_{\frac{}{}}\right]$  
    $\displaystyle +P_n^{(p)}(\sin\phi_k)p\frac{\frac{\partial y_k^{\prime}}{\partia...
...rtial\gamma_k}}{x^{\prime2}_k+y^{\prime2}_k}\left[-c_{np}\sin p\lambda_k\right.$  
    $\displaystyle +\left.s_{np}\cos p\lambda_k\right]_\frac{}{}\bigg\}\cdot$ (7)

In these two last formulae, $x_k^{\prime}$, $y_k^{\prime}$ and $z_k^{\prime}$ are the Cartesian coordinates of Pk in the planetocentric reference frame $(P_0,x^\prime ,y^\prime ,z^\prime )$ (see Fig. 1). In the case of Jupiter for which only the coefficients J2,J3,J4,J6,c22 and s22 are known (Campbell & Synnott 1985), Eq. (6) will be used with n=2,3,4,6 and Eq. (7) with n=p=2.

2.3 Equations of the ${\cal N}$ triaxial body problem in a planetocentric reference frame with fixed axes

In this subsection, we will now present the terms to be introduced to take into account the triaxialily of all the bodies.

  
2.3.1 Equations of motion for the centres of mass

For convenience of notation, we consider the same reference frame centered on the body P0 with fixed axes ( P0, 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 Pk exerted by Pl is now

 
$\displaystyle \vec{F}_{kl}=\vec{F}_{\bar{k}\bar{l}}+\vec{F}_{\hat{k}\bar{l}}+\vec{F}_{\bar{k}\hat{l}}+
\vec{F}_{\hat{k}\hat{l}}.$     (8)

More detailed study (see Krivov 1993 or Borderies 1978) of the term ${\vec F}_{\hat{k}\hat{l}}$ would show that it is a very small term, of the order of the product J2(k)J2(l) which we will neglect subsequently[*]. The terms of the form $\vec{F}_{\hat{0}\hat{l}}$ will also be neglected, for the same reason.

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

 
$\displaystyle \vec{A}+\vec{B}+\vec{C}=(m_0+m_i)\frac{\vec{F}_{\hat{\imath}\bar{...
...-\sum_{j\neq i, j\neq 0}^{\cal N}\frac{\vec{F}_{\bar{0}\hat{\jmath}}}{m_0}\cdot$     (9)

The first term (noted $\vec{A}$) is the direct perturbation ponctual planet-oblate satellite and a component of the indirect part of the same perturbation (for the mass mi). The second term (noted $\vec{B}$) is the direct mutual perturbations oblate satellite-ponctual satellite. The third term (noted $\vec{C}$) comes from the remainder of the indirect perturbations oblate satellite-ponctual planet.

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

  
                                 $\displaystyle \mbox{\vec{A}}$ = $\displaystyle -(m_0+m_i)\frac{\vec{F}_{\bar{0}\hat{\imath}}}{m_im_0}=-G(m_0+m_i)\nabla_0U_{\bar{0}\hat{\imath}}$ (10)
$\displaystyle \mbox{\vec{B}}$ = $\displaystyle \sum_{j\neq i,j\neq 0}^{\cal N}G\left(m_j\nabla_iU_{\bar{\imath}\hat{\jmath}}-m_j\nabla_jU_{\bar{\jmath}\hat{\imath}}\right)$ (11)
$\displaystyle \mbox{\vec{C}}$ = $\displaystyle -\sum_{j\neq i,j\neq 0}^{\cal N}Gm_j\nabla_0U_{\bar{0}\hat{\jmath}}$ (12)

where the notation $U_{\bar{k}\hat{l}}$ indicates the disturbing potential induced on the ponctual body Pk by the triaxiality of body Pl. However, the equations of Sect. 2.2 provide, more generally, the expression of ${\bf\nabla}_kU_{\bar{k}\hat{l}}$ in Cartesian coordinates centered on an unspecified body whose Euler angles ( $\psi_l, I_l, \chi_l$) are reported to the reference frame with fixed axes of directions (Pl,x,y,z) (see Fig. 2).
  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{0565fig2.eps}
\end{figure} Figure 2: Equatorial reference frame of a satellite Pl related to the reference frame of fixed axes ( P0,x,y,z).
Open with DEXTER

To obtain the expression of ${\bf\nabla}_kU_{\bar{k}\hat{l}}$ in the initial reference frame centered on body P0, it is enough to substitute formally in the expressions of $\frac{\partial U_{\bar{k}\hat{0}}^{(1)}}{\partial\gamma_k}$ and $\frac{\partial U_{\bar{k}\hat{0}}^{(2)}}{\partial\gamma_k}$ the coordinates ( xk, yk, zk) by ( xk-xl, yk-yl, zk-zl) and the angles $(\psi, I, \chi)$ by $(\psi_l, I_l, \chi_l)$. In the case of $\vec{A}$ and $\vec{C}$ the change of coordinates simply amounts to replacing ( xk, yk, zk) by ( -xk, -yk, -zk). 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 Pi the expression

 
                       $\displaystyle \vec{\ddot{r}}_i$ = $\displaystyle -G(m_0+m_i)\left(\frac{\vec{r}_i}{r_i^3}-{\bf\nabla}_iU_{\bar{\imath}\hat{0}}+\nabla_0U_{\bar{0}\hat{\imath}}\right)$  
    $\displaystyle +\sum_{j=1,j\neq i}^{\cal N}Gm_j\left(\frac{\vec{r}_j-\vec{r}_i}{...
...t{\imath}}
+\nabla_iU_{\bar{\imath}\hat{\jmath}}-\frac{\vec{r}_j}{r_j^3}\right.$  
    $\displaystyle +\left._\frac{}{}\nabla_jU_{\bar{\jmath}\hat{0}}-\nabla_0U_{\bar{0}\hat{\jmath}_\frac{}{}}\right).$ (13)

The explicit writing of Eq. (13) is rather complex; in practice it should be enough to use the first values of n and p for the coefficients Jn(k) and cnp(k) representing the triaxiality of the satellites. Furthermore very few of them are still known. Indeed, only the values of the coefficients J2(k) and c22(k) (which are the most significant) of the Galilean satellites were estimated from the trajectories of the space probes.

Other simplifications are also possible such as not considering the indirect perturbations (satellite-planet) for the coefficients J3(k), J4(k), J6(k) and even J2(k). These perturbations are additional oblateness forces, but they are multiplied by the term  $\frac{(E_r^{(k)})^n}{r_{kl}^{n+1}}$ with $n\geq2$. 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 Jn(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

 
$\displaystyle U_{\bar{k}\hat{0}}+U_{\bar{0}\hat{k}}=-\sum_{n=2}^\infty\frac{(E_...
...m_{m=2}^\infty\frac{(E_r^k)^m}{r_k^{m+1}}J_m^{(k)}P_m\left(\sin\phi_0^k\right).$     (14)

In the case where the north poles stay parallel, we have the equality $\phi_k=-\phi_0^k$. Then we have by oddness of the sinus function and the properties of Legendre polynomials
$\displaystyle U_{\bar{k}\hat{0}}+U_{\bar{0}\hat{k}} = -\sum_{n=2}^\infty\frac{(...
..._n(\sin\phi_k)\left(J_n+(-1)^n\left(\frac{E_r^k}{E_r}\right)^nJ_n^{(k)}\right).$     (15)

If the mutual perturbations oblate satellite-oblate satellite remain negligible, one can modify the values of the coefficients Jn of the central body consequently.

In addition, we can deduce that during the adjustment to the observations, the coefficients Jn will be correlated to the values of the coefficients Jn(k), if the polar oblateness of the satellites is not taken into account. Thus, if one wishes to use the values Jn(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 J2(k). Indeed, one just needs to replace the value of the J2 coefficient by the new value $(J_2+(\frac{E_r^k}{E_r})^2J_2^{(k)})$ during the integration of the satellite considered as oblate. For an analytical theory, it will be enough to modify the values of the coefficient J2 in the expression of the analytical series.

The coefficients c22(k) and s22(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 s22(k) are equal to zero. And so, for the Galilean satellites, only the coefficients c22(k) will have to be considered.

2.3.2 Precession and rotation of the central body

The system (1) is incomplete if the Euler angles $(\psi, I, \chi)$ and $(\psi_k,I_k,\chi_k)_{1\leq k\leq {\cal N}}$, 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)

 
$\displaystyle \left\{
\begin{array}{rcl}
\alpha_0&=&268.05-0.009T\\  [2mm]
\delta_0&=&64.49+0.003T\\  [2mm]
W&=&284.95+870.536642d
\end{array}\right.$     (16)

where T denotes the time in Julian century of 36 525 days and d the time in Julian days from the standard epoch J2000.

The system (16) is thus directly usable and remains connected to the variables of our system by the relations $I=90^\circ-\delta_0$, $\psi=\alpha_0+90^\circ$ and $\chi=W+W_0$ where W0 is the longitude of the central meridian line refered to the $(P_0,\tilde{x}^\prime)$ axis (see Fig. 1).

2.4 Post-Newtonian correction

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 ${\cal N}$-body problem (ponctual masses), developed in a Galilean reference frame using a relativistic formulation, given by de Sitter

                                   $\displaystyle \vec{\ddot{r}}_i$ = $\displaystyle -\sum_{j\neq i}Gm_j\frac{\vec{r}_i-\vec{r}_j}{r_{ij}^3}+\sum_{j\neq i}\left\{\frac{G^2}{c^2}m_im_j\frac{\vec{r}_i-\vec{r}_j}{r_{ij}^4}\right.$  
    $\displaystyle +4\sum_{k\neq i}\frac{G^2}{c^2}m_jm_k\frac{\vec{r}_i-\vec{r}_j}{r...
...\neq j, k\neq i}\frac{G^2}{c^2}m_jm_k\frac{\vec{r}_i-\vec{r}_j}{r_{ij}^3r_{jk}}$  
    $\displaystyle +\frac{4G}{c^2}m_j[\vec{\dot{r}}_i\cdot (\vec{r}_i-\vec{r}_j)]\fr...
...i-
\vec{\dot{r}}_j}{r_{ij}^3}+\frac{7G}{2c^2}m_j\frac{\vec{\ddot{r}}_j}{r_{ij}}$  
    $\displaystyle -\frac{3G}{c^2}m_j[\vec{\dot{r}}_j\cdot(\vec{r}_i-\vec{r}_j)]\frac{\vec{\dot{r}}_i
-\vec{\dot{r}}_j}{r_{ij}^3}$  
    $\displaystyle +\frac{G}{2c^2}m_j[2\vec{\dot{r}}_i+3(\vec{\dot{r}}_j\cdot \vec{n}_{ij})^2-4
(\vec{\dot{r}}_i-\vec{\dot{r}}_j)$  
    $\displaystyle -\left.\vec{\ddot{r}}_j\cdot(\vec{r}_j-\vec{r}_i)]\frac{\vec{r}_i-\vec{r}_j}{r_{ij}^3}\right\}$ (17)

where c is the speed of light, $\vec{n}_{ij}$ the unit vector of direction  $\overrightarrow{P_jP_i}$, so $\vec{n}_{ij}=\frac{{\bf r}_i-{\bf r}_j}{r_{ij}}$.

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 Pi the equation

 
                                $\displaystyle \vec{\ddot{r}}_i$ = $\displaystyle -\sum_{j\neq i}Gm_j\frac{\vec{r}_i-\vec{r}_j}{r_{ij}^3}+4\frac{G^2}{c^2}m_0^2\frac{\vec{r}_i-
\vec{r}_0}{r_{i0}^4}$  
    $\displaystyle +\frac{4G}{c^2}m_0[\vec{\dot{r}}_i\cdot (\vec{r}_i-\vec{r}_0)]\fr...
...i
-\vec{\dot{r}}_0}{r_{i0}^3}+\frac{7G}{2c^2}m_0\frac{\vec{\ddot{r}}_0}{r_{i0}}$  
    $\displaystyle -\frac{3G}{c^2}m_0\left[\vec{\dot{r}}_0\cdot(\vec{r}_i-\vec{r}_0)\right]\frac{\vec{\dot{r}}_i
-\vec{\dot{r}}_0}{r_{i0}^3}+\frac{G}{2c^2}m_0$  
    $\displaystyle \times \left[2\vec{\dot{r}}_i+3(\vec{\dot{r}}_0\cdot
\vec{n}_{i0}...
...c{\dot{r}}_i-\vec{\dot{r}}_0)-\vec{\ddot{r}}_0\cdot(\vec{r}_0-\vec{r}_i)\right]$  
    $\displaystyle \times \frac{\vec{r}_i-\vec{r}_0}{r_{i0}^3}$ (18)

and for the central body
 
$\displaystyle \vec{\ddot{r}}_0+\sum_{j\neq 0}Gm_j\frac{\vec{r}_0-\vec{r}_j}{r_{0j}^3}=0.$     (19)

Introducing Eq. (19) into Eq. (18) gives
                             $\displaystyle \vec{\ddot{r}}_i$ = $\displaystyle -\sum_{j\neq i}Gm_j\frac{\vec{r}_i-\vec{r}_j}{r_{ij}^3}+4\frac{G^2}{c^2}m_0^2\frac{\vec{r}_i-
\vec{r}_0}{r_{i0}^4}$  
    $\displaystyle +\frac{Gm_0}{c^2r_{i0}^2}\left\{[4\vec{\dot{r}}_i\cdot\vec{n}_{i0...
...\dot{r}}_0+\frac{3}{2}(\vec{\dot{r}}_0\cdot\vec{n}_{i0})^2
\vec{n}_{i0}\right\}$  
    $\displaystyle -\frac{Gm_0}{c^2r_{i0}^2}
\left[\vec{\dot{r}}_i^2+2\vec{\dot{r}}_0^2-4\vec{\dot{r}}_i\cdot \vec{\dot{r}}_0\right]\vec{n}_{i0}.$ (20)

Hence, taking into account the relativistic effects, in a Jovicentric reference frame, implies adding the final expression
$\displaystyle \ddot{\overrightarrow{OP_i}}-\ddot{\overrightarrow{OP_0}}-\ddot{\...
..._i\ {\bf\cdot}\
\vec{n}_i)\ \vec{\dot{r}}_i - \vec{\dot{r}}_i^2\ {\bf n}_i\big]$     (21)

to Eq. (13).

  
2.5 Integral of energy

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

$\displaystyle E=\sum_{i=1}^{\cal N}\frac{1}{2}m_i\vec{\dot{r}}_i^2-\frac{1}{2M}\left(\sum_{i=1}^{\cal N}m_i
\vec{\dot{r}}_i\right)^2-U$     (22)

where U is the gravitational potential of the system.

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 ${\cal N}+1$-triaxial bodies problem relative to Eq. (13); we have after neglecting the terms $\vec{B}$ and $\vec{C}$ of Eqs. (11) and (12) the equality

$\displaystyle U = \sum_{i=1}^{\cal N}Gm_im_0\left(\frac{1}{r_i}+U_{\bar{\imath}...
...ght)
+\sum_{i=1}^{{\cal N}-1}\sum_{k=i+1}^{{\cal N}}\frac{Gm_km_i}{r_{ik}}\cdot$     (23)

In this last expression we have to express $U_{\bar{\imath}\hat{0}}$, that we can obtain easily expressing the Legendre polynomias. In the same way one can write the potential $U_{\bar{0}\hat{\imath}}$, using Sect. 2.3.1.

  
3 Magnitude of the perturbations

We present the method we used to test the influence of each perturbation, as well as the numerical tools that we have developed.

  \begin{figure}
\par\includegraphics[angle=-90,width=8cm,clip]{0565fig3.ps}\vspace*{3mm}
\includegraphics[angle=-90,width=8cm,clip]{0565fig4.ps}
\end{figure} 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


  \begin{figure}
\par\includegraphics[angle=-90,width=8.5cm,clip]{0565fig5.ps}
\end{figure} Figure 4: Variations on positions (in kilometers) under the effect of the coefficients c22(k) with an exact spin-orbit resonance model for the satellites' rotation.
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=-90,width=8.5cm,clip]{0565fig6.ps}
\end{figure} Figure 5: Variations on positions (kilometers) after substracting a straight line, under the action of the coefficients c22(k).
Open with DEXTER

3.1 The test method

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 J2, J4, J6.

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

 
$\displaystyle \vec{\ddot{r}}_i=-\frac{G(m_0+m_i)\vec{r}_i}{r_{i}^3}
+\sum_{j=1,...
...^3}-\frac{\vec{r}_j}{r_{0j}^3}\right)
+Gm_0{\bf\nabla}_iU_{\bar{\imath}\hat{0}}$     (24)

but we considered instead (except for the additional oblateness forces) a more complete expression which preserves the integral of energy by adding the terms $\sum_{k=1}^{\cal N}\frac{\vec{F}_{\hat{0}\bar{k}}}{m_0}$, to find the complete Eq. (5).

3.2 The code

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).

3.2.1 Control of numerical error

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.

3.2.2 The Solar and planetary perturbations

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.

3.3 Analysing the unusual perturbations

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.

3.3.1 Study of the Jovian precession

The Jovian equatorial reference frame is not fixed in space. In particular, the angles $\psi$ and I for which the coefficients Jn 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.

  
3.3.2 Effect of the satellites' triaxiality

A, B and C in the equations of motion (see the page [*]). In Lainey et al. (2001), the study of J2(k) perturbation was estimated to 5000 kilometers over one century. We will here present the contribution of c22(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 J2(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 ( $\lambda_0^{(k)}=0$), 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 $l_1-3l_2+2l_3\simeq\pi$, where l denotes the mean longitude.

3.3.3 Variations on the amplitude and the frequency of Laplacian libration

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.

3.4 Summary of the perturbations

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.

3.4.1 Classification of the perturbations

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.

3.4.2 The perturbations selected

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 J2(k) and c22(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.

  
4 Partial derivatives

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 $\vec{p}=(p_1, \cdots, p_M)$, and the initial condition (which are the $6{\cal N}$ position-velocities of the ${\cal N}$ bodies at an instant of reference t0) of a theory. We denote these constant $\vec{c}=(\vec{r}_1(t_0),\cdots,
\vec{r}_{\cal N}(t_0);\vec{\dot{r}}_1(t_0),\cdots,$ $\vec{\dot{r}}_{\cal N}(t_0);\vec{p})$, or in short $\vec{c}=(\vec{r}^{(0)},\vec{\dot{r}}^{(0)},\vec{p})$. 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.

4.1 Equations of variations

Let us consider a set of observations $(t_k, \vec{r}_i^{\prime(k)})_{1\leq k\leq K}$ providing for discrete values of time tk the observed positions $\vec{r}_i^{\prime(k)}=\vec{r}_i^\prime(t_k)$. The theory (here the numerical integration), gives for these same values of time the computed positions $(t_k, \vec{r}_i^{(k)})$, as solution of a differential system. Denoting $\varphi_t
(\vec{r}^{(0)}, \vec{\dot{r}}^{(0)}, \vec{p})$ the dynamic flow of this system and $\tilde{\varphi}_t^i$ the projection in  ${{\mathbb R}^3}$ of $\varphi$ such that $\tilde{\varphi}_t^i=\vec{r}_i$, one has the equality

 
                            $\displaystyle \Delta \vec{r}_i^{(k)}$ = $\displaystyle \left(\vec{r}_i^{\prime(k)}-\vec{r}_i^{(k)}\right)$  
  = $\displaystyle \left(\tilde{\varphi}_{t_k}^i
\left(\vec{r}^{\prime(0)},\vec{\dot...
..._{t_k}^i\left(\vec{r}^{(0)},\vec{\dot{r}}^{(0)},
\vec{p}\right)\frac{}{}\right)$ (25)

where $\vec{r}^{\prime(0)},\vec{\dot{r}}^{\prime(0)},$ and $\vec{p}^\prime$ correspond to the constants of the true system.

From the Taylor formula, we have

$\displaystyle \left(\vec{r}_i^{\prime(k)}-\vec{r}_i^{(k)}\right)=\sum_{l=1}^{6{...
...hi}_{t_k}^i}{\partial c_l}(\vec{c})\cdot\Delta c_l+O\left((\Delta c_l)^2\right)$     (26)

where $\delta c_l=c_l^\prime-c_l$. Considering the $\delta c_l$ sufficiently small, we thus have
 
$\displaystyle \Delta \vec{r}_i^{(k)}\approx\sum_{l=1}^{6{\cal N}+M}\frac{\partial\tilde{\varphi}_{t_k}^i}{\partial c_l}(\vec{c})\cdot\Delta c_l$     (27)

which is a system of three equations with $6{\cal N}+M$ unknown for each observation. This system will be solved by the least-squares method. However let us notice that the last expression requires us to know for the times tk the values of $\frac{\partial\tilde{\varphi}_{t_k}^i}{\partial c_l}({\vec c})$, or more explicitly $\frac{\partial \tilde{\varphi}_{t_k}}{\partial {\vec r}^{(0)}}$, $\frac{\partial \tilde{\varphi}_{t_k}}{\partial \dot{{\vec r}}^{(0)}}$ and $\frac{\partial \tilde{\varphi}_{t_k}}{\partial {\vec p}}$. The calculation of these quantities can be done solving a differential system thanks to Moulton's method. We have
$\displaystyle \frac{{\rm d}^2\vec{r}_i}{{\rm d}t^2}=\frac{\vec{F}_i}{m_i}(\vec{r},\vec{\dot{r}},
\vec{p}) .$     (28)

Denoting cl an unspecified constant of the theory ( $\vec{r}^{(0)}, \vec{\dot{r}}^{(0)}$, or $\vec{p})$, we obtain
                               $\displaystyle \frac{\partial}{\partial c_l}\left(\frac{{\rm d}^2\vec{r}_i}{{\rm d}t^2}\right)$ = $\displaystyle \frac{\partial}{\partial c_l}\frac{\vec{F}_i}{m_i}(\vec{r},\vec{\dot{r}},
\vec{p})$  
  = $\displaystyle \frac{1}{m_i}\left(\frac{\partial \vec{F}_i}{\partial\vec{r}}\fra...
...{\partial\vec{F}_i}{\partial\vec{p}}\frac{\partial\vec{p}}{\partial c_l}\right)$ (29)

$\vec{F}_i$ is generally independent of $\vec{\dot{r}}$, so that denoting $\gamma_j$ any of the coordinates of a body Pj, we have
$\displaystyle \frac{\partial}{\partial c_l}\left(\frac{{\rm d}^2
\vec{r}_i}{{\r...
...\partial \vec{F}_i}{\partial p_m}\frac{\partial p_m}{\partial c_l}\right) \cdot$     (30)

As cl is independent of time, we have the equality
$\displaystyle \frac{{\rm d}^2}{{\rm d}t^2}\left(\frac{\partial \vec{r}_i}{\part...
...}{{\rm d}t^2}\right)=\frac{1}{m_i}\frac{\partial \vec{F}_i}{\partial c_l} \cdot$     (31)

That leads to the following differential system
 
$\displaystyle \frac{{\rm d}^2}{{\rm d}t^2}\left(\frac{\partial \tilde{\varphi}_...
...\frac{\partial
\vec{F}_i}{\partial p_m}\frac{\partial p_m}{\partial c_l}\right)$     (32)

of the form $\frac{{\rm d}^2\vec{X}^i}{{\rm d}t^2}=\sum_{j=1}^{\cal N}\vec{g}_i^j(t)
\vec{X}^j+\vec{h}_i(t)$ with initial conditions $(0, \cdots 0,1,0,\cdots 0)$ (the 1 being located at the row l) for each component of $\vec{X}^i$. We finally need to solve $6{\cal N}+M$ differential systems (for each constant of the theory), of $3{\cal N}$ equations each.

  
4.2 Form of the equations of variation

The differential system (32) requires us to write the vectors $\frac{\partial \vec{F}_i}{\partial \vec{r}}$ and $\frac{\partial \vec{F}_i}{\partial \vec{p}}$ as functions of the variables $(\vec{r}, \vec{\dot{r}}, \vec{p})$. First, we take for the expression of $\vec{F}_i$ the one given in Eq. (5) when modeling the satellites as ponctual

                              $\displaystyle \frac{\vec{F}_i}{m_i}$ = $\displaystyle -\frac{G(m_0+m_i)\vec{r}_i}{r_i^3}
+\sum_{j=1,j\neq i}^{\cal {\ca...
...}Gm_j\left(\frac{\vec{r}_j-\vec{r}_i}{r_{ij}^3}-
\frac{\vec{r}_j}{r_j^3}\right)$  
    $\displaystyle +G(m_0+m_i){\vec \nabla}_iU_{\bar{\imath}\hat{0}}+\sum_{j=1,j\neq i}^{\cal {\cal N}}Gm_j
{\vec \nabla}_jU_{\bar{\jmath}\hat{0}}.$  

Let us write now $\frac{\partial \vec{F}_i}{\partial \vec{r}}$. Denoting $\gamma_n$ any of the Cartesian coordinates of body Pn, after calculations
                            $\displaystyle \frac{1}{m_i}\frac{\partial\vec{F}_i}{\partial\gamma_n}$ = $\displaystyle \frac{G(m_0+m_i)}{r_i^3}\left[\frac{3\gamma_n\vec{r}_i}{r_i^2}-\zeta\right]\delta_{in}$  
    $\displaystyle +\frac{Gm_n}{r_n^3}\left(\frac{3\gamma_n\vec{r}_n}{r_n^2}-\zeta\right)(1-\delta_{ni})$  
    $\displaystyle +\frac{Gm_n}{r_{in}^3}\left(\zeta-\frac{3(\gamma_n-\gamma_i)(\vec{r}_n-
\vec{r}_i)}{r_{in}^2}\right)(1-\delta_{ni})$  
    $\displaystyle -\sum_{j=1,j\neq i}^{\cal N}\frac{Gm_j}{r_{ij}^3}\left(\zeta+\frac{3(\gamma_i-\gamma_j)(
\vec{r}_j-\vec{r}_i)}{r_{ij}^2}\right)\delta_{in}$  
    $\displaystyle +\left(G(m_0+m_i)\frac{\partial}{\partial \gamma_n}\nabla_iU_{\bar{\imath}\hat{0}}\right)\delta_{in}$  
    $\displaystyle +\left(Gm_n\frac{\partial}{\partial\gamma_n}\nabla_nU_{\bar{n}\hat{0}}\right)(1-\delta_{in})$ (33)

where $\delta_{in}$ is the Kronecker symbol, and $\zeta=1$ if $\gamma_n$ is of the type of $\vec{F}_i$, $\zeta=0$ if not[*].

Equation (32) can thus be rewritten, denoting $\vec{f}_{il}=\frac{\partial\tilde{\varphi}_t^i}{\partial c_l}$, in the form

$\displaystyle \frac{{\rm d}^2\vec{f}_{il}}{{\rm d}t^2} = \frac{1}{m_i}\sum_{n=1...
...c{F}_i}{\partial z_n}f_{nl}^{(3)}+\frac{\partial\vec{F}_i}{\partial c_l}\right)$     (34)

which is under an explicit form
 
                             $\displaystyle \frac{{\rm d}^2\vec{f}_{il}}{{\rm d}t^2}$ = $\displaystyle \frac{G(m_0+m_i)}{r_i^3}\left[\frac{3(\vec{r}_i\cdot
\vec{f}_{il})\ \vec{r}_i}{r_i^2}- \vec{f}_{il}\right]$  
    $\displaystyle +G(m_0+m_i)\ \vec{J}_{il}+\sum_{n=1,n\neq i}^{\cal N}Gm_n\ \vec{J}_{nl}$  
    $\displaystyle +\sum_{n=1,n\neq i}^{\cal N}\Bigg\{\frac{Gm_n}{r_n^3}\left(\frac{3(\vec{r}_n\cdot
\vec{f}_{nl}) \vec{r}_n}{r_n^2}-\vec{f}_{nl}\right)$  
    $\displaystyle +\frac{Gm_n}{r_{in}^3}\left[_{\frac{}{}}(\vec{f}_{nl}-\vec{f}_{il})_{\frac{}{}}\right.$  
    $\displaystyle +\left.\frac{3(\vec{r}_n-\vec{r}_i)\cdot(\vec{f}_{il}-\vec{f}_{nl})}{r_{in}^2}
(\vec{r}_n-\vec{r}_i)\right]_{\frac{}{}}\Bigg\}$  
     $\displaystyle +\frac{1}{m_i}\frac{\partial \vec{F}_i}{\partial c_l}$ (35)

where
$\displaystyle \vec{J}_{kl}=\left(\frac{\partial}{\partial x_k}\nabla_kU_{\bar{k...
...ft(\frac{\partial}{\partial z_k}\nabla_kU_{\bar{k}\hat{0}}\right)f_{kl}^{(3)} .$     (36)

It is clear that the term $\frac{\partial\vec{F}_i}{\partial c_l}$ is equal to zero when the adjusted constant is an initial condition. In the opposite case, its expression results easily from the developments written in Sect. 2.2. Moreover, the explicit calculation of the terms in  $\vec{J}_{kl}$ has the form
 
                         $\displaystyle \frac{\partial}{\partial\Gamma_k}
\left(\frac{\partial U_{\bar{k}\hat{0}}}{\partial\gamma_k}\right)$ = $\displaystyle \sum_{n=2}^\infty\frac{(E_r)^nJ_n}{r_k^{n+2}}\Bigg\{\bigg(
\frac{(n+1)P_n(\sin\phi_k)}{r_k}$  
    $\displaystyle +\frac{\sin\phi_k\frac{{\rm d}P_n(\sin\phi_k)}{{\rm d}\sin\phi_k}...
...\bigg)
\cdot\left.\left(\zeta-\frac{(n+3)\Gamma_k\gamma_k}{r_k^2}\right)\right.$  
    $\displaystyle +\frac{\partial\sin\phi_k}{\partial\Gamma_k}\bigg(\frac{(n+2)\gamma_k
\frac{{\rm d}P_n(\sin\phi_k)}{{\rm d}\sin\phi_k}}{r_k}$  
    $\displaystyle +\frac{\gamma_k\sin\phi_k\frac{{\rm d}^2P_n(\sin\phi_k)}{{\rm d}\sin^2\phi_k}}{r_k}$  
    $\displaystyle -\frac{\partial z_k^\prime}{\partial\gamma_k}
\frac{{\rm d}^2P_n(\sin\phi_k)}{{\rm d}\sin^2\phi_k}\bigg)$  
    $\displaystyle +\frac{(n+2)\Gamma_k\left(\frac{\partial z_k^\prime}{\partial\gamma_k}
\frac{{\rm d}P_n(\sin\phi_k)}{{\rm d}\sin\phi_k}\right)}{r_k^2}\Bigg\}\cdot$ (37)

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 $U_{\hat{\jmath}\bar{\imath}}$ (term $\vec{B}$ of Eq. (9)). We replace in Eq. (35) the expression of  $\vec{J}_{kl}$ by the following one

                                $\displaystyle \vec{J}_{kl}$ = $\displaystyle \left(\frac{\partial}{\partial x_k}\nabla_kU_{\bar{k}\hat{0}}+\frac{\partial}{\partial x_0^k}\nabla_0U_{\bar{0}\hat{k}}\right)f_{kl}^{(1)}$  
    $\displaystyle +\left(\frac{\partial}{\partial y_k}\nabla_kU_{\bar{k}\hat{0}}+\frac{\partial}{\partial y_0^k}\nabla_0U_{\bar{0}\hat{k}}\right)f_{kl}^{(2)}$  
    $\displaystyle +\left(\frac{\partial}{\partial z_k}\nabla_kU_{\bar{k}\hat{0}}+\frac{\partial}{\partial z_0^k}\nabla_0U_{\bar{0}\hat{k}}\right)f_{kl}^{(3)}$ (38)

where (x0k,y0k,z0k) are the coordinates of body P0 in the reference frame centered on body Pk. Then one can refer to Eq. (37) for an explicit writing.

In practice, Eq. (35) have been integrated simultaneously with the equations of motion.

  
5 Adjustment to E5 ephemerides

We have adjusted our model to the Sampson-Lieske theory using the method presented in the preceding section.

  
5.1 Solution E5

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 $(\epsilon_k, \beta_k)$ 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 $6{\cal N}$ initial conditions (positions and velocities) denoted  $(\vec{r}_k,\vec{\dot{r}}_k)$ of the four satellites;
-
the masses mk of the four satellites;
-
the mass m0 of Jupiter;
-
the Jovian oblateness coefficients J2 and J4;
-
the angles $(\psi, I)$ relating to the orientation of the Jovian pole.
Four iterations using the least-squares method were carried out to optimize our adjustment. To interpret the differences obtained, we choose to represent them in a Jovian equatorial frame with cylindrical coordinates $(\rho ,\nu ,z)$ and thus very close to the variables of Sampson (1921).

5.1.1 Results

The results are shown in Fig. 6 using the variables $(\rho_i,\nu_i,z_i)$. 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 $\nu_i$ and in a less significant way the variable zi. To understand such long period differences, we examine the variables used by Sampson.

5.1.2 The problem induced by cylindrical coordinates

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 $\zeta_i$, or equally to the coordinates zi. Indeed, the nature of this coordinate (rather inappropriate to describe a disturbed Keplerian motion) leads to high oscillations in the amplitudes of zi 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 zi at all.

  \begin{figure}
\par\includegraphics[angle=-90,width=5.8cm,clip]{0565fig7.ps}\hsp...
...mm}
\includegraphics[angle=-90,width=5.8cm,clip]{0565fig18.ps}\par
\end{figure} 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 $(\rho ,\nu ,z)$ ( 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 $\nu_i$ than the variable zi. During our adjustment we had previously kept constant the values of the angles $(\psi, I)$. The corresponding results then gave a stronger variation on zi than that on $\nu_i$. However, we ended up also adjusting the angles $(\psi, I)$. The variations in zi were then reflected on $\nu_i$ 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.

5.2 Comparison with the observation accuracy

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.

5.3 Contribution of the small perturbations

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.

5.3.1 Adjustment of the main system to our final 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 $(\rho_i,\nu_i,z_i)$, 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.

5.3.2 Results

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 J2 and J4 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 $5 \times 10^{-6}$ whereas we obtain variations of $6,14 \times 10^{-6}$ and $3,58 \times 10^{-5}$. 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 J2(k) of the satellites in the oblateness of the central body.

Finally, it may be surprising that the variations are more significant for J4 than for J2. This effect comes directly from the correlation between these two coefficients. Actually, to minimize the variations in J4, 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.

6 Conclusion

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.

References



Copyright ESO 2004