A&A 387, 700-709 (2002)
DOI: 10.1051/0004-6361:20020420

A new determination of lunar orbital parameters,
precession constant and tidal acceleration
from LLR measurements

J. Chapront - M. Chapront-Touzé - G. Francou

Observatoire de Paris - SYRTE - UMR 8630/CNRS, 61 avenue de l'Observatoire, 75014 Paris, France

Received 13 December 2001 / Accepted 5 March 2002

Abstract
An analysis of Lunar Laser Ranging (LLR) observations from January 1972 until April 2001 has been performed, and a new solution for the lunar orbital motion and librations has been constructed that has been named S2001. With respect to prior solutions, improvements in the statistical treatment of the data, new nutation and libration models and the addition of the positions of the observing stations to the list of fitted parameters have been introduced. Globally, for recent observations, our rms (root mean square error) is within 2 to 3 centimeters in the lunar distance. Special attention has been paid to the determination of the correction to the IAU76 luni-solar constant of precession, and the value of the secular acceleration of the Moon's longitude due to the tidal forces. The main results are:
-  correction to the constant of precession: $\Delta p$ = -0.302 $\pm$ $0.003 ''/{\rm cy}$,
-  tidal acceleration of the lunar longitude: $\Gamma$ = -25.858 $\pm$ $0.003 ''/{\rm cy}^{2}$.
The positions and velocities of the stations have also been determined. The results are consistent with the ITRF2000 determinations from SLR observations. The lunar theory ELP is referred to a dynamical system and introduces the inertial mean ecliptic of J2000.0. The positioning of the reference system of the theory with respect to ICRS is performed (and also with respect to some useful JPL numerical integrations). Finally the orientation of the celestial axes with respect to the ICRS reference system has been derived as well as the offsets of the Celestial Ephemeris Pole.

Key words: planets and satellites: individual: moon - astrometry - reference systems


   
1 The new solution S2001

Several analyses of LLR observations have been performed using the lunar theory ELP2000-96 (Chapront & Chapront-Touzé 1997) and an improved version of the lunar libration theory of Moons (1984) with numerical and analytical complements (Chapront et al. 1999a). A previous analysis, described in (Chapront et al. 1999b), covered the time span January 1972 until March 1998. It is referred to below as the solution S1998. We refer to this paper for the presentation of the principles of the analysis, discussion of the methods and results of the comparisons. More recently a new analysis was performed using LLR observations of McDONALD and CERGA before May 2000 (Chapront et al. 2000). Several improvements were introduced in the lunar ephemerides, mainly in the libration model, and also in the program of reduction (an up-to-date nutation model) and in the statistical treatment of the data (an adequate distribution of weights among the various observing stations and periods of observations). This solution is referred to as S2000. Now, on the same basis, we have enlarged the time span of observations until April 2001 and added a few parameters in the program of reduction (positions and velocities of the observing stations). This last solution is referred to as S2001.

In this paper, we shall also mention intermediate solutions in which the characteristics of the fit are the same as in S2001; the only change is the upper limit of the time span covered by the observations, which may vary within 5 years (1996-2001). This interval has been chosen considering that recent observations, more accurate than earlier ones, have a much larger weight in the determination of the fitted parameters. These intermediate solutions were done to study the evolution of the fitted values with the time interval of observations, in particular for tidal acceleration (Fig. 3), precession constant (Fig. 4), and obliquity (Fig. 5).

   
2 The residuals


 

 
Table 1: LLR residuals: time distribution of the rms (in centimeter). N is the number of normal points involved.
OBSERVATORY Time S1998 S2000 Time S2001 N
and instruments Interval rms rms Interval rms  
             
             
McDONALD 1972-1986 34.7 34.5 1972-1975 43.5 1487
Telescope 2.70 m       1976-1979 27.7 1035
and MLRS1       1980-1986 29.1  990
             
             
CERGA Rubis 1984-1986 18.2 18.8 1984-1986 18.7 1165
             
             
HALEAKALA 1987-1990 11.1  8.0 1987-1990  6.3  451
             
             
McDONALD 1987-1998  5.0   1987-1991  5.8  232
MLRS2 1987-2000    3.8 1991-1995  4.6  586
        1995-2001  3.3 1669
             
             
CERGA Yag 1987-1998  4.8   1987-1991  5.3 1574
  1987-2000    3.8 1991-1995  3.9 2044
        1995-2001  3.0 3273


Table 1 shows the residuals in distance and illustrates the global precision of the solution S2001. As it was already mentioned for S2000, an important gain of precision has been obtained in our new solutions compared to S1998. Consequently the unknowns in S2001 are determined with a better accuracy than in S1998.

In order to estimate the evolution in the quality of the observations we illustrate in Fig. 1 the time distribution of rms obtained with S2001, for the data provided by the 2 modern instruments: MLRS2 for McDONALD and Yag for the CERGA.

Since 1991 we have observed smaller residuals for the CERGA except around 1997. It is worthwhile to note that this period corresponds to an offset in the CERGA measurements (Mangin 1998), which is taken into account in our analysis by a global correction of 0.7 ns for the observations from January 13, 1997 to June 24, 1998. During the same period, we also observe in the determination of various parameters a kind of "accidental jump'' (see, for example, in Fig. 4, the "jump'' occurring during this period for the correction to the IAU76 constant of precession).


  \begin{figure} \par\includegraphics[angle=-90,width=8.8cm,clip]{MS2201f1.eps}\end{figure} Figure 1: Time evolution of rms for the 2 stations McDONALD and CERGA covering 14 years (solution S2001).
Open with DEXTER

   
3 The fitted parameters

We list below the parameters that are fit in the solutions S1998 and S2000. All the angles and mean motions are referred to J2000.0.
- The geocentric lunar orbital parameters W1(0), W2(0), W3(0) (constants of the mean longitude and mean longitudes of perigee and node), $\nu=W_{1}^{(1)}$, $\Gamma$, E (sidereal mean motion, constants for inclination and eccentricity).
- The heliocentric orbital parameters of the Earth-Moon barycenter T(0), $\varpi ^{(0)}$ (constants of the mean longitude and mean longitude of perihelion), n', e' (sidereal mean motion and eccentricity).
- The bias parameters $\Delta W_{1}^{(2)}$, $\Delta W_{2}^{(1)}$, $\Delta W_{3}^{(1)}$ (observed corrections to the computed coefficient of the quadratic term of the lunar mean longitude, and the computed mean motions of perigee and node). $\Delta W_{1}^{(2)}$ yields an observed value of W1(2,T), the tidal part of the coefficient of the quadratic term of the mean longitude (half tidal secular acceleration).
- The 6 free libration parameters (parameters tied to the coefficients of the main free libration terms and values of the free libration arguments).
- The $3\times4$ reflector coordinates. The reflector coordinates are referred to lunar principal axes of inertia.
- The position angles $\phi$, $\epsilon$ and $\psi$ with respect to different systems of axes. Figure 2 illustrates the relative positions of various systems presented in Sect. 4.
- A correction to precession $\Delta p$: optional parameter.

In the solution S2001 we keep the same list as above and we add $5\times3\times 2$ optional parameters giving the positions and velocities of the 5 stations: McDONALD 2.70 m, MLRS1 and MLRS2, CERGA and Haleakala. The parameters are the equatorial rectangular coordinates X, Y, Z in the ITRF (position) and their derivatives $\dot{X}$, $\dot{Y}$, $\dot{Z}$ (velocity). Note that simultaneously fitting all the parameters has not been possible. The fits have been performed in several steps, but tests have been made in order to check the stability of the results. Indeed, strong correlations exist among some parameters that may weaken the accuracy of our determinations; in particular, it is the case of the variables related to the reference frame ($\phi$ and $\epsilon$) and the positions of the stations (X, Y, Z). $\Delta p$ (precession) and $\dot{\epsilon}$ (obliquity rate) are correlated with the velocities of the stations ($\dot{X}$, $\dot{Y}$, $\dot{Z}$); $\Delta p$ and the principal nutation term are also difficult to separate. We have adopted the following strategy. First, we determine the whole set of parameters mentioned above except the positions and velocities of the stations. Then fixing the value of $\phi$, we add the positions of the stations to the whole set and make a new improvement. Next we determine the velocities of the stations separately. Finally, fixing all the parameters, we perform a last analysis including $\Delta p$ and the principal term of $\Delta\psi$ (nutation in longitude) (see Sect. 10). At each step of the process we verify the coherence of the determinations; for example we verify that the introduction of the fitted values of X, Y, Z does not change significantly the value of $\phi$ if the first step is reiterated.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{MS2201f2.eps}\end{figure} Figure 2: Relative positions of the mean inertial ecliptic of J2000.0 with respect to ICRS, MCEP and JPL.
Open with DEXTER

   
4 Position angles of the inertial mean ecliptic J2000.0

We recall here the definition of the position angles of the inertial mean ecliptic J2000.0 with respect to various "equatorial'' reference systems (R). R stands either for ICRS (International Celestial Reference System), MCEP (Reference linked to the Mean Celestial Ephemeris Pole of J2000.0) or JPL (Reference system defined by a JPL numerical integration such as DE200, DE403 or DE405). We set:
- $\gamma_{2000}^{I}(R)$: ascending node of the inertial mean ecliptic J2000.0 on the equator of R;
- $\epsilon (R)$: Inclination of the inertial mean ecliptic on the equator of R;
- o(R): Origin of right ascensions in the equator of R;
- $\phi (R)$: Arc $o(R) \gamma_{2000}^{I}(R)$;
- $\psi (R)$: Arc $\gamma_{2000}^{I}({\rm ICRS})\gamma_{2000}^{I}(R)$.
Two solutions are investigated. They are denoted Sol. 1 (MCEP) and Sol. 2 (ICRS), corresponding to the reference systems MCEP or ICRS. In the reduction of LLR observations one has to transform the terrestrial coordinates of the station to celestial ones. Such a transformation depends of the daily values of the polar motion $x_{{\rm p}}, y_{{\rm p}}$, the difference UT1-UTC and a precession nutation matrix $P\times N$ which rotates the celestial instantaneous axes to a J2000.0 fixed celestial "equatorial'' system of axes.

In Sol. 1 (MCEP), the matrix $P\times N$ is provided by analytical solutions: polynomial expressions for the orientation of the Earth's equator (Williams 1994) and IERS 1996 theory of nutation (McCarthy 1996). The reference plane is the mean equator of the CEP for J2000.0. The corresponding system of axes is that of the MCEP.

In Sol. 2 (ICRS), $P\times N$ is computed via a system of corrections to a conventional set of values for the nutations in longitude and obliquity, i.e. $\delta\psi$ and $\delta\epsilon$ which are daily values provided by IERS (series C04). The corresponding system is the ICRS.

In both solutions, corrections to the precession constant and to the obliquity are fit. Sol. 1 (MCEP) involves the theoretical value of the obliquity rate due to (Williams 1994), -46.8340''/cy. Corrections to the principal terms of the nutations in longitude and obliquity are also fit in Sol. 1 (MCEP) (see Sect. 10).

   
5 Orbital motion


 

 
Table 2: Solution S2001: corrections to the nominal values of the orbital parameters formerly fit to DE200 (in arcsecond, except for $\nu $ and n', in arcsecond/cy).
     
Variable Sol. 1 (MCEP) Sol. 2 (ICRS)
     
     
W1(0) -0.1218 $\pm$ 0.0002 -0.0775 $\pm$ 0.0002
W2(0) -0.0673 $\pm$ 0.0002 -0.0229 $\pm$ 0.0002
W3(0) -0.1155 $\pm$ 0.0005 -0.0717 $\pm$ 0.0005
$\nu $ -0.3978 $\pm$ 0.0008 -0.4033 $\pm$ 0.0008
$\Gamma$   0.0008 $\pm$ 0.0000   0.0009 $\pm$ 0.0000
E   0.0002 $\pm$ 0.0000   0.0002 $\pm$ 0.0000
T(0) -0.0770 $\pm$ 0.0002 -0.0326 $\pm$ 0.0002
$\varpi'^{(0)}$ -0.0587 $\pm$ 0.0003 -0.0140 $\pm$ 0.0003
n'   0.0304 $\pm$ 0.0008   0.0258 $\pm$ 0.0008
e'   0.0000 $\pm$ 0.0000   0.0000 $\pm$ 0.0000



 

 
Table 3: Solution S2001: fitted value of tidal part of the quadratic term of the mean longitude (in arcsecond/cy2) and observed corrections to the mean motion of perigee and node (in arcsecond/cy).
Variable Sol. 1 (MCEP) Sol. 2 (ICRS)
     
     
W1(2,T) -12.9257 $\pm$ 0.0016 -12.9290 $\pm$ 0.0016
$\Delta W_{2}^{(1)}$    0.0334 $\pm$ 0.0008     0.0317 $\pm$ 0.0008
$\Delta W_{3}^{(1)}$  -0.3806 $\pm$ 0.0097   -0.3611 $\pm$ 0.0097


The main results are summarized in Tables 2 and 3. They can be compared to analogous tables in Chapront et al. (1999b, 2000). The uncertainties reported in these tables are formal errors.

Table 2 shows the corrections to the nominal values of the orbital parameters of the Moon which were formerly fit to JPL ephemerides DE200. These basic parameters are regarded as nominal values in our lunar ephemeris ELP2000. They have been used in our work since 1982 (Chapront-Touzé & Chapront 1983, 1988) and maintained to facilitate further comparison. The origin of angles is $\gamma_{2000}^{I}({\rm MCEP})$ in Sol. 1 (MCEP) and $\gamma _{2000}^{I}({\rm ICRS})$ in Sol. 2 (ICRS). The differences between the angles W1(0), W2(0), W3(0) in Sol. 1 (MCEP) and in Sol. 2 (ICRS) are close to the value of :

\begin{eqnarray*}\psi ({\rm MCEP}) = \gamma_{2000}^{I}({\rm ICRS}) \gamma_{2000}^{I}({\rm MCEP}) = 0.0445'' \end{eqnarray*}


given below (see Table 9). For all the quantities, the differences, corrected by $\psi$ for angles, are less than 3$\sigma$ except for $\nu $ and n'. We note that the differences $\nu _{{\rm MCEP}}-\nu _{{\rm ICRS}}$ and $n'_{{\rm MCEP}}-n'_{{\rm ICRS}}$ are almost equal, around 0.0050''/cy, and do not vary significantly with the time interval of observations in the intermediate solutions (1996-2001) (see Sect. 1). This suggests a slight drift between the two frames that is not yet completely clarified.


 

 
Table 4: Angular mean elements of the Moon and the Earth-Moon barycenter and Delaunay arguments deduced from S2001. The origin of the angles is $\gamma _{2000}^{I}({\rm ICRS})$ (see Fig. 2); t is the time in Julian centuries reckoned from J2000.0.
W1 = $218^{\circ}18'59\farcs878~2 ~+1~732~559~343\farcs332~8~t ~-\hphantom{1}6\farcs870~0~t^2 ~+0\farcs006~604~t^3 ~-0\farcs000~031~69~t^4$
W2 = $\hphantom{1}83^{\circ}21'11\farcs651~8 ~+\hphantom{1}\hphantom{1}~14~643~420\... ...s330~4~t ~-38\farcs263~9~t^2 ~-0\farcs045~047~t^3 ~+0\farcs000~213~01~t^4$
W3 = $125^{\circ}02'40\farcs326~5 ~-\hphantom{1}\hphantom{1}\hphantom{1}~6~967~919\... ...~+\hphantom{1}6\farcs359~3~t^2 ~+0\farcs007~625~t^3 ~-0\farcs000~035~86~t^4$
T = $100^{\circ}27'59\farcs188~0 ~+\hphantom{1}~129~597~742\farcs301~6~t ~-\hphantom{1}0\farcs020~2~t^2 ~+0\farcs000~009~t^3 ~+0\farcs000~000~15~t^4$
$\varpi'$ = $102^{\circ}56'14\farcs413~6 ~+\hphantom{1}\hphantom{1}\hphantom{1}\hphantom{1... ...tom{1}0\farcs532~7~t^2 ~-0\farcs000~138~t^3 \hphantom{~+0.000~000~00''~t^0}$
l = $134^{\circ}57'48\farcs226~4 ~+1~717~915~923\farcs002~4~t ~+31\farcs393~9~t^2 ~+0\farcs051~651~t^3 ~-0\farcs000~244~70~t^4$
l' = $357^{\circ}31'44\farcs774~4 ~+\hphantom{1}~129~596~581\farcs073~3~t ~-\hphantom{1}0\farcs552~9~t^2 ~+0\farcs000~147~t^3 ~+0\farcs000~000~15~t^4$
F = $\hphantom{1}93^{\circ}16'19\farcs551~7 ~+1~739~527~263\farcs217~9~t ~-13\farcs229~3~t^2 ~-0\farcs001~021~t^3 ~+0\farcs000~004~17~t^4$
D = $297^{\circ}51'00\farcs690~2 ~+1~602~961~601\farcs031~2~t ~-\hphantom{1}6\farcs849~8~t^2 ~+0\farcs006~595~t^3 ~-0\farcs000~031~84~t^4$


Table 3 shows the observed value W1(2,T), and the bias in the mean motions W2(1), W3(1). It is worthwhile to note the agreement of the three quantities determined independently in the two solutions; it was not the case in S1998 for W1(2,T) and W3(1).

Table 4 gives the expressions of the angular mean elements of the Moon and of the Earth-Moon barycenter, and the corresponding Delaunay arguments deduced from the S2001 Sol. 2 (ICRF):
- the mean longitude and the mean longitudes of perigee and node of the Moon (W1, W2, W3);
- the heliocentric orbital parameters of the Earth-Moon barycenter (T, $\varpi'$);
- the Delaunay arguments: l=W1-W2, $l'=T-\varpi '$, F= W1- W3 and $D=W_1-T+180^{\circ}$.
The mean motions are referred to J2000.0 and the origin of the angles is $\gamma _{2000}^{I}({\rm ICRS})$ (see Fig. 2).


 

 
Table 5: Fit of the positions and velocities of the two observing LLR stations CERGA (Yag) and McDONALD (MLRS2). Units: positions X, Y, Z in meters; velocities $\dot{X},~\dot{Y},~\dot{Z}$ in meter/year.
CERGA (Yag) X Y Z $\dot{X}$ $\dot{Y}$ $\dot{Z}$
             
             
References            
ITRF1994 Epoch 1993.0  4581692.254   556195.961  4389355.019 -0.0122  0.0194  0.0084
ITRF1996 Epoch 1997.0  4581692.217   556196.039  4389355.075 -0.0120  0.0189  0.0106
ITRF2000 Epoch 1997.0  4581692.181   556196.024  4389355.072 -0.0131  0.0189  0.0101
Corrections            
S2001 Sol. 2 (ICRS) -0.012 -0.006 -0.001 -0.0008  0.0000    0.0036
Comparisons            
ITRF96 - ITRF2000  0.036  0.015  0.003      
S2001*  - ITRF96 -0.027 -0.008 -0.017      
S2001*  - ITRF2000  0.009  0.007 -0.014      
S2001** - ITRF2000  0.005  0.007  0.001      
             
             
McDONALD (MLRS2) X Y Z $\dot{X}$ $\dot{Y}$ $\dot{Z}$
             
             
References            
ITRF1994 Epoch 1993.0 -1330021.390 -5328403.336  3236481.726 -0.0116 -0.0036 -0.0073
ITRF1996 Epoch 1997.0 -1330021.431 -5328403.340  3236481.672 -0.0117 -0.0034 -0.0064
ITRF2000 Epoch 1997.0 -1330021.440 -5328403.330  3236481.675 -0.0125 -0.0001 -0.0065
Corrections            
S2001 Sol. 2 (ICRS) -0.006  0.005 -0.031 -0.0004 -0.0007   -0.0002
Comparisons            
ITRF96 - ITRF2000  0.009 -0.010 -0.003      
S2001*  - ITRF96 -0.015  0.009 -0.003      
S2001*  - ITRF2000 -0.006 -0.001 -0.006      
S2001** - ITRF2000 -0.008 -0.004 -0.007      


   
6 Positions and motions of the stations

In the solutions S2000 and earlier we set the positions and motions of the stations to constant values deduced from the ITRF94. Increasing the interval of observations, we noticed slight trends of quantities related to the reference frame, in particular in the obliquity of the ecliptic $\epsilon$. Furthermore, with other sets of positions and velocities for the stations, deduced from ITRF96 or ITRF2000, we observed significant changes in the obliquity. Hence, we modified our list of fitted parameters in the analysis, adding the positions and velocities of the stations as new unknowns. As initial values, we adopted the positions deduced from ITRF94 and the velocities of ITRF2000.
The main results are the following:
- The trend in $\epsilon$, when the time interval of observations increases, is significantly reduced.
- The residuals for Haleakala are significantly reduced from 8 cm in S2000 to 6 cm in S2001 (see Table 1);
- The corrections to the positions of the stations are in good agreement with the positions given in ITRF2000.
The last point calls for some comments. The corrections to the positions determined either with a large set of parameters (see Sect. 3) or independently are very close, which supports the relevance of our present determinations. The corrections obtained in Sol. 1 (MCEP) and Sol. 2 (ICRS), are also very close.

In Table 5 we present our fit in Sol. 2 (ICRS) for the positions of the 2 operating stations: CERGA (Yag) and McDONALD (MLRS2). The "references'' are deduced from ITRF publications (Boucher et al. 1996, 1998); in the case of MLRS2 we used corrections of eccentricity to the nearby SLR station (Ries 1999). The "corrections'' S2001 Sol. 2 (ICRS) must be applied to ITRF94 positions (Epoch 1993.0) and to ITRF2000 velocities, and have been fit independently. For the "comparisons'', our corrected positions at the epoch 1997.0 were computed in two different ways: using the ITRF2000 velocities (positions S2001*) or using the corrected velocities (positions S2001**). The positions S2001* show a better agreement with ITRF2000 than with ITRF96.

For the CERGA the results S2001** are even closer. Nevertheless, the correction to $\dot{Z}$ (0.0036 m/year) seems to be too large compared to the reference value (0.0101 m/year in ITRF2000). New determinations of the positions of LLR and SLR stations at CERGA have been obtained recently by (Nicolas 2000), using the observations of the satellites LAGEOS 1 and 2 from 1997 to 1999. They show local and seasonal displacements of the LLR station, reaching a few centimeters with a probable tendency of a few mm/year. These small displacements are ignored in ITRF2000, but the recent LLR observations are sensitive to this level of accuracy. The order of magnitude of such displacements is in accordance with our corrections.

   
7 The tidal acceleration of the Moon

The tidal component of the secular acceleration of the Moon's longitude is a fundamental parameter of the evolution of the Earth-Moon system. The tidal dissipation is due to a misalignment of the tidal bulge of the Earth relative to the Earth-Moon direction. This bulge exerts a secular torque and most of the effect comes from the ocean tides. It produces a secular negative acceleration of the Moon and a decrease in the Earth's rotation rate (increase of the length of day). A consequence of the negative acceleration in the lunar longitude of approximately $-25.8''/{\rm cy}^{2}$ is the well-known displacement of the Moon's barycenter that corresponds to an increase of the Earth-Moon distance of 3.8 cm/year.

Let us examine the quadratic terms of the lunar mean longitude W1(2) as it appears in the analytical lunar theory ELP2000-96 (Chapront et al. 1997). Here the constants and parameters are those from DE245, in particular Love numbers and time delays (see Sect. 8). W1(2) contains mainly 3 contributions (in $''/{\rm cy}^{2}$):
[1]       5.8665    Planetary perturbations;
[2]       0.1925    Earth's figure perturbations;
[3]  -12.8125    Tides.
These contributions have to be multiplied by 2 if we speak of "secular acceleration'' components.

In [1] the main effect is due to the secular variation of the solar eccentricity. It represents in fact a Taylor expansion of long periodic perturbations, mentioned also as "secular perturbations'' in classical celestial mechanics, with periods of the solar perigee and node of several ten thousand years.

In [2] the contribution is of the same nature. Hence, [1] and [2] are due to very long period effects. They have to be distinguished from the tides that induce dissipative forces. An analytical solution for the lunar motion also provides secular terms of higher degree in time (t3-terms, t4-terms, etc.). With the knowledge of all the secular components in the longitude, we are able to isolate the tidal acceleration from other perturbations.

We gather in Table 6 a non-exhaustive list of determinations of the tidal secular acceleration of the lunar longitude. The most recent values have been obtained with LLR observations. We note, for this type of determination, a significant improvement of the precision with more observations. It is also worth noticing that the most recent determination around $-25.86''/{\rm cy}^{2}$ comes closer to the value of Morrison & Ward (1975), $-26''/{\rm cy}^{2}$, obtained with an analysis of optical observations including occultations and planetary transits covering a time interval of approximately 2 centuries.


   
Table 6: Tidal acceleration of the lunar mean longitude (in arcsecond/cy2).
Authors  Value  Publication
Spencer Jonesa -22 1939
Oesterwinter & Cohena -38 1975
Morrison & Warda -26 1975
Mullerb -30 1976
Calame & Mulhollandc -24.6 1978
Ferrari et al.d -23.8 1980
Dickey et al.c -23.8 1982
Dickey and Williamsc -25.10 1982
Newhall et al.c -24.90 1988
Chapront Touzé et al.c -25.62 1997
(Solution S1998)c -25.78 1999b
(Solution S2000)c -25.836 2000
(this paper, solution 2001)c -25.858  


Type of observations: a Occultations, b eclipses, c LLR, d LLR and Lunar orbiter.


It is interesting to examine the intrinsic values of the tidal acceleration in various JPL numerical integrations (see Table 7). These values do not appear explicitly in the lists of parameters that are provided with each lunar ephemeris, but we have computed them from those parameters and from the models described in Sect. 8. It is worth noticing that the difference between the tidal secular acceleration in DE405 and our determination in S2001 is about $0.03''/{\rm cy}^{2}$ which gives an idea of the present uncertainty for this fundamental lunar parameter.


 

 
Table 7: Tidal acceleration of the lunar mean longitude in various JPL ephemerides (in arcsecond/cy2) (Standish 1982, 1995, 1998).
JPL ephemeris Value Publication
DE200    -23.895    1982
DE245    -25.625    1990
DE403    -25.580    1995
DE405    -25.826    1998


The expression of the lunar mean longitude of the Moon W1 has the following secular expansion: W1 = W1(0) + W1(1) t + W1(2) t2 +..., where t is the time in century reckoned from J2000.0; W1(0) is the constant term; W1(1) is the sidereal mean motion for J2000.0; W1(2) is the total half secular acceleration of the Moon. We have investigated the convergence of these quantities to the S2001 values given in Table 4, when the upper limit of the time interval of observations increases. The results are given in Fig. 3 where $\Delta W_1^{(0)}$, $\Delta W_1^{(1)}$ and $\Delta W_{1}^{(2)}$ are the corrections to S2001 values obtained in the intermediate solutions (1996-2001). In particular the variation of $\Delta W_{1}^{(2)}$ shows the evolution of the fitted value of the tidal acceleration when using more and more recent LLR observations: as mentioned above W1(2) is the sum of several contributions listed [1], [2] and [3] at the beginning of the present section, and we may assume that [1] and [2] have been computed with sufficient accuracy through the secular terms in ELP theory; only [3] (half tidal acceleration) needs to be fit. This graph allows us to ensure nowadays a realistic precision in the knowledge of the secular acceleration of better than $0.03''/{\rm cy}^{2}$, in agreement with the conclusion resulting from the comparison of the S2001 value to the DE405 one.


  \begin{figure} \par\includegraphics[angle=-90,width=8.8cm,clip]{MS2201f3.eps}\end{figure} Figure 3: Time evolution of the corrections $\Delta $ to the secular components of the mean longitude of the Moon, W1 = W1(0) + W1(1) t + W1(2) t2, when increasing the set of LLR measurements.
Open with DEXTER

   
8 An analytical approach to the tidal perturbations

The tidal perturbations of the Moon have various origins:
[1] the Earth deformations due to the Moon;
[2] the Earth deformations due to the Sun;
[3] the deformations of the Moon itself by the Earth;
[4] the deformations of the Moon itself by the Sun.

The main effect arises from [1]. Below, we restrict ourselves to this tidal component. The additional potential of a non rigid Earth acting on the Moon has the well-known classical form (Lambeck 1980):

\begin{eqnarray*}&& \Delta U=\frac{GM}{r^*}\sum_{l=2}^{\infty}\left(\frac{R}{r}\... ...P_{lm}(\sin \phi)P_{lm}(\sin \phi^*) \cos m(\lambda-\lambda^*), \end{eqnarray*}


where G is the gravitational constant and M the mass of the Moon, R the equatorial radius of the Earth; l and m are integers. kl,m are Love numbers and Pl,m Legendre functions. $\lambda$, $\phi$ and r are the spherical coordinates of a point outside of the Earth at time t, in terrestrial axes; $\lambda^{*}$, $\phi^{*}$ and r* are the same quantities for the Moon's barycenter, the symbol star (*) meaning that the coordinates are evaluated at the time $t^{*} = t - \tau_{l,m}$ where $\tau_{l,m}$ is a delay in the deformation related to the harmonic of index (l,m) attached to the Love number kl,m. For an elastic Earth: $\tau_{l,m}=0$; in case of an anelastic Earth: $\tau_{l,m}\not=0$. If now we limit ourselves to l=2 and express the above formula in terms of right ascension $\alpha$ and declination $\delta$, we derive for the disturbing function acting on the Earth-Moon vector:

\begin{eqnarray*}&& \Delta U=GM\left( 1+\frac{M}{E}\right)\frac{R^5}{{r^*}^3r^3}... ...{22}(\sin\delta^*) \cos(2\alpha-2\alpha^*-2\omega\tau_{22})~], \end{eqnarray*}


where $\omega$ is the angular velocity of the Earth, and E is the mass of the Earth. The index (2, 1) induces the diurnal tides; the index (2, 2) induces the semi-diurnal ones.

In the simplified model of the JPL numerical integration DE200 the following approximations where done: k20=k21=k22=k; $\tau_{21}=\tau_{22}=\tau$; $(\omega \tau)^{2}$ is neglected.

In DE245, and further JPL integrations (DE403 and DE405), one simply puts $\tau_{20}=0$, and k21=k22. In the analytical solution ELP we follow the same way, and substitute in $\Delta U$ analytical series for the lunar coordinates. After integration of the differential equations, the main effect consists of two contributions in the mean longitude:

\begin{eqnarray*}\Delta W_{1} = a t^{2} + b \cos \Omega. \end{eqnarray*}



 

 
Table 8: Contributions of harmonics in the evaluation of a and b coefficients in the lunar mean longitude $\Delta W_{1} = a t^{2} + b \cos \Omega $. k: Love number; $\tau $: time delay.
Parameter (2, 0) (2, 1) (2, 2) Total
         
         
DE200        
k    0.30   0.30 0.30  
$\tau $ (day)    0   0.006460 0.006460  
a ( $''/{\rm cy}^{2}$)   -0.90 -10.55   -11.45  
b (10-5'')    0   58 -138   -80  
         
         
DE403        
k    0.34   0.30 0.30  
$\tau $ (day)    0   0.014350 0.006772  
a ( $''/{\rm cy}^{2}$)   -2.0 -11.06   -13.06  
b (10-5'')    0   127 -144   -17  
         
         
DE405        
k    0.34   0.30 0.30  
$\tau $ (day)    0   0.012909 0.006942  
a ( $''/{\rm cy}^{2}$)   -1.8 -11.34   -13.14  
b (10-5'')    0   114 -148   -34  



 

 
Table 9: Position angles of the inertial mean ecliptic of J2000.0 with respect to equatorial celestial system (R) in S2001. The uncertainties are formal errors. Units: arcsecond.
R $\epsilon - 23^{\circ}26'21\arcsec$ $\phi$ $\psi$ Mean Epoch
         
         
ICRS 0.41100 $\pm$ 0.00005 -0.05542 $\pm$ 0.00011   Dec. 1994
MCEP 0.40564 $\pm$ 0.00009 -0.01460 $\pm$ 0.00015 0.0445  $\pm$  0.0003 Dec. 1994
DE403 0.40928 $\pm$ 0.00000 -0.05294 $\pm$ 0.00001 0.0048  $\pm$  0.0004 Jan. 1985
DE405 0.40960 $\pm$ 0.00001 -0.05028 $\pm$ 0.00001 0.0064  $\pm$  0.0003 Jan. 1990


The first term corresponds to the secular acceleration $2\times a$. The second term is a periodic term in $\Omega$ with the period of the ascending node of the Moon on the ecliptic (18.6 years).

We gather in Table 8 various evaluations of a and b depending on the model, and the corresponding values of Love numbers and delay. We notice that the b coefficient is smaller in the case of DE403 and DE405 than for DE200. It is mainly due to a partial cancellation of the coefficients arising from (2, 1) and (2, 2), which is not the case for DE200. It should be noted that the value of the secular acceleration ($2\times a$) from Table 8, Col. 5, is not complete as mentioned above. It represents nevertheless the main tidal contributions that are slightly different from the total contributions given in Table 7.

In JPL numerical integrations, Love numbers and time delays are fitted parameters, while in our solution only a is fit. In both cases, because of the long period of the argument $\Omega$, a long and accurate set of data is necessary. Hence, we understand why the value ($2\times a$) of the secular acceleration has been significantly improved since DE200, which was fit on less that 15 years of observations (Williams et al. 1978). A better knowledge of this parameter has benefited from an increasing range of data and an improvement in the quality of observations.

   
9 Orientation of the celestial axes

We gather in Table 9 our new determinations of the position angles $\phi$, $\epsilon$ and $\psi$ (see Sect. 4). The angle $\psi$ is obtained through the 2 different versions of the mean longitude of the Moon W1, which is evaluated in the ICRS (with Sol. 2), and in the reference system R, R being the MCEP with Sol. 1, or the reference system of any JPL numerical integration DEn with the solution fit to the lunar ephemeris of DEn. Such a solution is deduced from an analysis of the same nature as those using the LLR observations themselves, by considering the JPL lunar ephemeris as an "observational model''. Hence, we compute the difference:

\begin{eqnarray*}\psi(R)&=& W_{1}(ICRS)- W{_1}(R) \\ \hphantom{\psi(R)}&=& W_{1... ...\hphantom{\psi(R)}&&+[W_{1}^{(2)}(ICRS) - W_{1}^{(2)}(R)]t^{2}, \end{eqnarray*}


where t is the time at a mean epoch, reckoned in centuries from J2000.0.

The mean epoch for MCEP arises directly from the least-squares fit; it is the weighted time related to the distribution of weights of the sub-groups. For the JPL ephemerides, the mean epochs are mentioned in the literature and correspond to JPL's fits. In $\psi (R)$ the linear term corresponds to the difference of sidereal mean motions in the 2 systems: $\nu(ICRS)-\nu(R)$; the quadratic term corresponds to half the difference between the tidal parts of the acceleration of the mean longitude in the 2 systems: W1(2,T)(ICRS)-W1(2,T)(R), the non tidal parts being the same.

The $\psi$ function that we obtain for DE405 is:
$\psi(DE405) = 0.00832 + 0.01793 t - 0.01555 t^{2}$ (arcsecond).
From (Standish 2000), the reference system tied to the ephemeris is based on VLBI observations (Magellan spacecraft to Venus and Phobos approach to Mars) made between 1989 and 1994. Hence, we have arbitrarily chosen 1990 Jan. 1 for the mean epoch of DE405.

Using the quantities $\phi$ and $\psi$ of Table 9, we make the projection on the ICRS "equator'' of the origin of right ascension o(DE405) which is distant from o(ICRS) by less than one mas. We find:

o(ICRS)o(DE405) = 0.7 mas (Epoch, 1990 Jan. 1).

For DE403 we obtain:

o(ICRS)o(DE403) = 1.9 mas (Epoch, 1985 Jan. 1).

These results concerning DE405 agree with the fact that the numerical integration is oriented onto the ICRS (Standish 1998).

   
10 The precession constant

The correction to the IAU76 constant of precession in (Chapront et al. 2000) was:

$\displaystyle % \Delta_{1}p = -0.3164 \pm 0.0027 ''/{\rm cy}.$     (1)

This value was presented at the IAU meeting in Manchester (Fukushima 2000a). It was obtained via the solution S2000, using Sol. 1 (MCEP), as mentioned above. With Sol. 1 (MCEP) of S2001, we get now:

\begin{eqnarray*}\Delta_{1}p = -0.3364 \pm 0.0027 ''/{\rm cy}. \end{eqnarray*}


The noticeable divergence of this new determination with respect to the previous one in (1) caught our attention and we noticed that in Sol. 2 (ICRS) a residual $\Delta_{2}p$ arises:

\begin{eqnarray*}\Delta_{2}p = -0.0316 \pm 0.0027 ''/{\rm cy}. \end{eqnarray*}


Sol. 1 (MCEP) and Sol. 2 (ICRS) use the same observations (the main source of errors) and the same models, except for the motion of the reference frame due to precession and nutation. Hence, the differences between the corrections $\Delta_{1}p$ and $\Delta_{2}p$ are mainly due to the precession-nutation models. If we assume that C04 $\delta\psi$ and $\delta\epsilon$ series used in Sol. 2 (ICRS) based on VLBI observations, contributes ideally to the precession-nutation matrix $P\times N$, the difference $\Delta_{1}p - \Delta_{2}p$ gives an estimate of the corrections that should be applied to the IAU76 precession constant. In particular this difference eliminates the effects of an improper motion of the stations, errors in the EOP series for Universal Time and polar motion, and a local bias produced by the observations themselves. Contrariwise all systematic effects or local errors in C04 $\delta\psi$ and $\delta\epsilon$ series are still in the difference.


 

 
Table 10: Correction to the IAU76 constant of precession $\Delta p$ (in arcsecond/cy) and offsets of Celestial Ephemeris Pole at J2000.0 $-\psi \sin\epsilon$ and $\Delta \epsilon $ (in arcsecond). The uncertainties are formal errors.
Method Source $\Delta p$ $-\psi \sin\epsilon$ $\Delta \epsilon $
         
         
VLBI Fukushima (2000a)  -0.297  $\pm$  0.004  -0.0167 $\pm$ 0.0005 -0.0049 $\pm$ 0.0003
LLR S1998  -0.344  $\pm$  0.004  -0.0183 $\pm$ 0.0004 -0.0056 $\pm$ 0.0002
LLR S2000  -0.316  $\pm$  0.003  -0.0173 $\pm$ 0.0004 -0.0054 $\pm$ 0.0002
LLR S2001 (this paper)  -0.302  $\pm$  0.003  -0.0177 $\pm$ 0.0004 -0.0054 $\pm$ 0.0002
  IAU 2000A (Mathews et al. 2002)  -0.29965    


If we suppose that C04 $\delta\psi$ and $\delta\epsilon$ series do not contain any secular trends or bias, a more appropriate correction to the IAU76 constant of precession should then be done with:

$\displaystyle \Delta p = \Delta_{1}p - \Delta_{2}p = -0.3048 ''/{\rm cy}.$     (2)

Figure 4 illustrates the evolution of $\Delta_{1}p$ and $\Delta p = \Delta_{1}p - \Delta_{2}p$ when the upper limit of the time span covered by the fit varies. We see that, though $\Delta_{1}p$ varies, $\Delta p$ remains constant around the value -0.302''/cy, and we propose now this value as the LLR correction value.


  \begin{figure} \par\includegraphics[angle=-90,width=8.8cm,clip]{MS2201f4.eps}\end{figure} Figure 4: Evolution of the correction to the IAU76 constant of precession with the upper limit of the time span covered by the fit.
Open with DEXTER

The obliquity shows a similar phenomenon illustrated by Fig. 5. We note $\Delta_{1}\epsilon$ the correction obtained with Sol. 1 (MCEP) to a reference value of obliquity, $\Delta_{2}\epsilon$ the similar quantity obtained with Sol. 2 (ICRS), and $\Delta \epsilon $ the difference: $\Delta\epsilon = \Delta_{1}\epsilon - \Delta_{2}\epsilon$. A trend in $\Delta_{1}\epsilon$ is apparent in Fig. 5 (about $0.008\arcsec/$cy) but the difference $\Delta \epsilon $ is almost a constant. The trend in $\Delta_{1}\epsilon$ could make one believe in a correction to the obliquity rate -46''.8340/cy adopted in Sol. 1 (MCEP) from Williams (1994), but the similar trend in $\Delta_{2}\epsilon$ makes this hypothesis vanish. As for precession, the trends in $\Delta_{1}\epsilon$ and $\Delta_{2}\epsilon$ are rather due to an improper motion of the stations or to a local bias produced by the observations themselves. Note that trends are larger when positions and velocities of the stations are not improved (see Sect. 6).


  \begin{figure} \par\includegraphics[angle=-90,width=8.8cm,clip]{MS2201f5.eps}\end{figure} Figure 5: Evolution of the correction to the obliquity with the upper limit of the time span covered by the fit.
Open with DEXTER

In Table 10, we bring together our LLR determinations for the correction to the IAU76 constant of precession with the last values obtained by VLBI (Fukushima 2000b and Mathews et al. 2002), and the best estimates for the offsets of Celestial Ephemeris Pole at J2000.0, $-\psi \sin\epsilon(MCEP)$ and $\Delta \epsilon=\epsilon (MCEP)-\epsilon(ICRS)$. The last two quantities are denoted as $\theta_2$ and $-\theta_1$ in Chapront et al. (1999b), and $\Delta\psi \sin\epsilon_0$ and $\Delta \epsilon_0$ by Fukushima. We note that our values for $\Delta p$ are significantly different in S2001 and S1998. The nutation model and the weight distribution are deciding factors for the improvement of the solution. Now the value for $\Delta p$ obtained by LLR and VLBI converge nicely with a separation smaller than 0.03 mas/year.

We have also performed an analysis including the precession and the principal terms of nutations in longitude and obliquity. Although there is a strong correlation between precession and nutation, the final correction to the above $\Delta p$ is small ( +0.0082 ''/cy) and the amplitudes of the principal terms in nutations are not sensibly modified (0.5 mas for longitude and 0.02 mas for obliquity) within the formal errors.

11 Conclusion

The complete set of LLR observations now covers a time interval longer than 30 years. During the last ten years the precision in the measurements has been improved noticeably. Presently an individual measurement shows at CERGA an error around 30-60 ps in time, which corresponds to 5 to 10 mm in the one-way distance (Chapront & Mignard 2000). Hence the quality in the determination of several parameters of the Earth-Moon system has been improved correspondingly. This is the case in particular for the precession constant and the secular acceleration in the Moon's longitude. Increasing the precision and the length of the observing time, the models have to be refined. As mentioned above, small trends in the obliquity and bias in the mean motions of the Sun and Moon have not been completely eliminated from our analysis and that should be improved in the near future.

Acknowledgements
The authors are grateful to the LLR staff at CERGA and McDONALD for providing the observations used. They thank D. D. McCarthy for helpful comments on the manuscript.

References

 

Copyright ESO 2002