A&A 412, 567-586 (2003)
DOI: 10.1051/0004-6361:20031539

Expressions for IAU 2000 precession quantities

N. Capitaine¹ - P. T. Wallace² - J. Chapront ¹

1 - Observatoire de Paris, SYRTE/UMR8630-CNRS, 61 avenue de l'Observatoire, 75014 Paris, France
2 - H.M. Nautical Almanac Office, Space Science and Technology Department, CLRC / Rutherford Appleton Laboratory, UK

Received 4 June 2003 / Accepted 19 September 2003

Abstract
A new precession-nutation model for the Celestial Intermediate Pole (CIP) was adopted by the IAU in 2000 (Resolution B1.6). The model, designated IAU 2000A, includes a nutation series for a non-rigid Earth and corrections for the precession rates in longitude and obliquity. The model also specifies numerical values for the pole offsets at J2000.0 between the mean equatorial frame and the Geocentric Celestial Reference System (GCRS). In this paper, we discuss precession models consistent with IAU 2000A precession-nutation (i.e. MHB 2000, provided by Mathews et al. 2002) and we provide a range of expressions that implement them. The final precession model, designated P03, is a possible replacement for the precession component of IAU 2000A, offering improved dynamical consistency and a better basis for future improvement. As a preliminary step, we present our expressions for the currently used precession quantities $\zeta_{\rm A}, \theta_{\rm A}, z_{\rm A}$, in agreement with the MHB corrections to the precession rates, that appear in the IERS Conventions 2000. We then discuss a more sophisticated method for improving the precession model of the equator in order that it be compliant with the IAU 2000A model. In contrast to the first method, which is based on corrections to the t terms of the developments for the precession quantities in longitude and obliquity, this method also uses corrections to their higher degree terms. It is essential that this be used in conjunction with an improved model for the ecliptic precession, which is expected, given the known discrepancies in the IAU 1976 expressions, to contribute in a significant way to these higher degree terms. With this aim in view, we have developed new expressions for the motion of the ecliptic with respect to the fixed ecliptic using the developments from Simon et al. (1994) and Williams (1994) and with improved constants fitted to the most recent numerical planetary ephemerides. We have then used these new expressions for the ecliptic together with the MHB corrections to precession rates to solve the precession equations for providing new solution for the precession of the equator that is dynamically consistent and compliant with IAU 2000. A number of perturbing effects have first been removed from the MHB estimates in order to get the physical quantities needed in the equations as integration constants. The equations have then been solved in a similar way to Lieske et al. (1977) and Williams (1994), based on similar theoretical expressions for the contributions to precession rates, revised by using MHB values. Once improved expressions have been obtained for the precession of the ecliptic and the equator, we discuss the most suitable precession quantities to be considered in order to be based on the minimum number of variables and to be the best adapted to the most recent models and observations. Finally we provide developments for these quantities, denoted the P03 solution, including a revised Sidereal Time expression.

Key words: astrometry - reference systems - ephemerides - time

1 Introduction

The IAU precession-nutation model in use until the implementation of the IAU 2000 Resolutions was composed of the IAU 1976 precession (Lieske et al. 1977) and IAU 1980 nutation (Wahr 1981; Seidelmann 1982). IAU Resolution B1.6 adopted in 2000 recommended that these models be replaced, beginning on 1 January 2003, by the IAU 2000 precession-nutation model: specifically the MHB 2000 model provided in Mathews et al. (2002) designated IAU 2000A or, for applications not requiring the utmost accuracy, its shorter version IAU 2000B.

The IAU 2000 precession-nutation model includes a new nutation series for a non-rigid Earth and corrections to the precession rates in longitude and obliquity. The revised precession-nutation model is oriented with respect to the International Celestial Reference System (ICRS) through a fixed 3D rotation between the mean equatorial frame at J2000.0 and the Geocentric Celestial Reference System (GCRS). This rotation, called the frame bias, includes the numerical values for the pole offset at J2000.0 that MHB 2000 specifies and a third number, the equinox offset at J2000.0, that MHB 2000 does not specify. The adopted equinox offset has only a second-order effect on the final transformation between celestial and terrestrial coordinates.

IAU Resolution B1.7 recommended that the motion of the Celestial Intermediate Pole (CIP) in the GCRS be realized "by the IAU 2000A model for precession and forced nutation for periods greater than two days plus additional time-dependent corrections provided by the International Earth Rotation and Reference Systems Service IERS (IERS) through appropriate astro-geodetic observations" (i.e. through VLBI observations).

It should be noted that the pre-2003 VLBI procedures used the IAU 1976 precession and the "total" nutations (i.e. the nutations themselves plus the contribution of the corrections to the precession rates plus the biases) and omitted the equinox offset. Consequently the MHB 2000 model left room for interpretation, in respect of how the frame bias was to be implemented, how the new precession rates were to be applied and, in particular, what was to be done about the unspecified equinox offset.

The implementations of the precession-nutation models IAU 2000A and B set out in the IERS Conventions 2000 follow the straightforward approach of updating the secular terms of precession only. Corresponding software implementations exist (IERS and Standards Of Fundamental Astronomy (SOFA)) and offer a variety of tools catering for a wide variety of applications, both classical and "CEO-based'' (i.e. based on the use of the Celestial Ephemeris Origin (CEO), cf. Sect. 2). These interpretations of the IAU 2000 resolution, although of practical utility for the next few years, are in fact dynamically inconsistent and suffer, except of course for the improvements in the precession rates, from the same limitations as the IAU 1976 precession in the precision of the coefficients and compliance with up to date models for the ecliptic motion.

An improved IAU 2000 precession model is therefore necessary. With the MHB 2000 precession rates as a starting point, it is possible to develop precession quantities consistent with IAU 2000 that are dynamically consistent. This requires solution of the dynamical equations for the precession motion of the celestial pole based on the MHB 2000 precession rates and on improved expressions for the motion of the ecliptic.

As well as more accurate models for the precession per se, additional products are possible, such as simple polynomial-based algorithms for generating the rotation matrix and mean-pole X,Ythat combine the frame bias and the precession.

In this paper we present expressions for the precession quantities consistent with the IAU 2000A model, as recommended in IAU 2000 Resolution B1.6, following the approach described above. Our approach makes a clear distinction between the precession of the ecliptic due to planetary perturbations and the precession of the equator due to the luni-solar and planetary torques on the oblate Earth, both motions being expressed with respect to inertial space.

We will use abbreviations for quoting papers to which we often refer in this work. These are respectively L77 for Lieske et al. (1977), MHB for Mathews et al. (2002), S94 for Simon et al. (1994) and W94 for Williams (1994). The present paper is designated P03.

Note that L77 uses a two-epoch formulation, allowing mean place at some starting epoch to be transformed to mean place at some final epoch without the requirement for either epoch to be J2000. In their formulation, t is time from the fundamental epoch J2000 to the final epoch, whereas T represents time from the starting epoch to J2000. Given the pre-eminence of J2000 as the fundamental epoch, it is now reasonable to reject the two-epoch approach as an unnecessary complication, especially as the same result can be achieved by two successive transformations, the first from the starting epoch to J2000 and the second from J2000 to the final epoch. For several years the IERS Conventions have provided developments for the single time argument t from the J2000.0 epoch, and so does W94 for the revised developments. It seems clear that the future developments for precession will use a single time argument, and we shall follow this approach here.

The final expressions provided in this paper are denoted "the P03 solutions''. They include (i) expressions (37) and (38) for the primary precession quantities relative to the ecliptic and equator, respectively, (ii) expressions (39) to (41) for the derived precession quantities for classical use, (iii) expressions (45), (49) and (50) for alternative quantities and (iv) expressions (42) and (43) for revised Sidereal Time. Note that the unit of time used in all the expressions of the paper is Julian century, denoted cy.

   
2 Transformation formulas

There are two equivalent bias-precession-nutation transformations from GCRS to ITRS, namely the new (CEO-based) transformation and the classical (equinox-based) transformation. These transformations are based on two different origins on the equator with quite different properties. The equinox is defined geometrically and has a complex and comparatively rapid motion along the instantaneous equator that is a consequence of the motion not only of the moving equator but of the moving ecliptic as well. The CEO, which is an implementation of the non-rotating origin (Guinot 1979) as recommended in IAU 2000 Resolution B1.8, in contrast is defined kinematically: from one moment to the next, it moves only at right-angles to the instantaneous equator, and no ecliptic is involved. This almost complete separation between the treatment of the precessing-nutating pole and the origin of "right ascension'' leads to a much simpler relationship between stellar hour angles and Universal Time (for more detail, see for example Capitaine et al. 2003a).

We use the symbol $\vec{R}$ for the GCRS-to-ITRS rotation matrix, omitting polar motion, with elements  $\vec{R}(i,j)$.

Using the usual notation, the CEO-based transformation is written as:

$$\vec{R}_{\rm new} = R_3(\theta-E-s) \cdot R_2(d)\cdot R_3(E),$$ (1)

where $\theta$ is the Earth Rotation Angle and

 

$$\begin{array}{rcl} &&E = \arctan(Y/X), \\ &&d = \arctan \left(\l... ...^{2}\right)\big/\left(1-X^{2}-Y^{2}\right)\right) ^{\!1/2} \right). \end{array}$$ (2)

X(t) and Y(t) are the components of the CIP unit vector in the GCRS, based on the IAU 2000A precession-nutation model and the corresponding biases $(\xi_0, \eta_0)$ and equinox offset at epoch $d\alpha_0$, and s provides the position of the CEO.

For practical reasons, X and Y are usually called "coordinates'' and their numerical expressions are multiplied by the factor $1~ 296 ~ 000''/2 \pi$ in order to provide in arcseconds the value of the corresponding "angles'' (strictly their sines) with respect to the z-axis of the GCRS.

The classical form of the transformation is written as

$$\vec{R}_{\rm class} = \vec{T} ~ \vec{N} ~ \vec{P} ~ \vec{B} ,$$ (3)

i.e. as the product of the individual rotation matrices $\vec{B}$ (bias) followed by $\vec{P}$ (precession) then $\vec{N}$(nutation) and finally $\vec{T}$ (Earth rotation):

 

$$\begin{array}{rcl} &&\vec{B} = R_1(-\eta_0) \cdot R_2(\xi_0) \cdo... ...ta\psi) \cdot R_1(\epsilon_{\rm A}) \\ &&\vec{T} = R_3({\rm GST}). \end{array}$$ (4)

The classical precession quantities $\psi _{\rm A}$, $\omega _{\rm A}$, $\epsilon_{\rm A}$ and  $\chi_{\rm A}$ are those defined by L77 and the nutation quantities  $\Delta\psi$ and $\Delta\epsilon$ are the luni-solar and planetary nutations. GST is Greenwich (apparent) Sidereal Time. Note that the precession matrix, $\vec{P}$, can be formed in several ways (W94), depending on which of the precession angles are used.

The four-angle formulation given above was chosen for the IERS/SOFA implementation of IAU 2000 (Sect. 4) because it enabled the specified precession-rate adjustments to be applied directly and unambiguously. In the case of the present work that is no longer a consideration and other choices are open. For example, the most common method is the three-angle formulation using  ${\zeta _{\rm A}}$, ${\theta _{\rm A}}$ and $z_{\rm A}$:

$$\vec{P} = R_3(-z_{\rm A}) \cdot R_2(+\theta_{\rm A}) \cdot R_3(-\zeta_{\rm A}).$$ (5)

A revision of this 3-rotation formulation could bring some of the benefits of the MHB 2000 corrections to existing users without forcing them to make major changes to their procedures. However, this advantage would only be realized if the frame bias could be taken care of at the same time, and it turns out that this is not straightforward. The possibility of a concise formulation that does take into account both precession and frame bias will be considered in Sect. 7.3.1.

3 The IAU 1976 precession

3.1 Basis of the developments for precession

The development of the precession quantities depends upon models for the dynamical motion of (i) the ecliptic pole, relative to a fixed ecliptic, due to planetary perturbations and (ii) the celestial pole, due to luni-solar and planetary torques on the oblate Earth.

The basic quantities for the motion of the mean ecliptic pole are $\sin\pi_{\rm A} \sin\Pi_{\rm A}$ and $\sin\pi_{\rm A} \cos\Pi_{\rm A}$, which express the components of the Earth's orbital angular momentum upon a fixed ecliptic and which depend on the set of planetary masses adopted.

The basic quantities for the precession of the equator are the angles  $\psi _{\rm A}$ and  $\omega _{\rm A}$ of L77 that, although thought of as being the luni-solar precession in longitude and obliquity respectively, are actually the orientation parameters of the mean equator of date in the mean ecliptic frame at epoch. Note that the additional precession quantities considered in (4) are $\epsilon_{\rm A}$ (mean obliquity of date) for the inclination of the mean equator of date on the ecliptic of date and  $\chi_{\rm A}$ (the so-called "planetary precession'') for the contribution to the motion of the equinox that is due only to the precession of the ecliptic.

A model for the motion of the celestial pole can be derived from the dynamical equation (see Eq. (14) of L77) expressing the motion of the mean pole of date about the ecliptic pole, once given the values at the reference epoch for the mean obliquity of the ecliptic, $\epsilon_0$, for the speed of precession, and for the geodesic precession (de Sitter & Brouwer 1938). The expression for general precession (denoted $p_{\rm A}$ in L77) combines the precession in longitude of the equator and the precession of the ecliptic, the former being a function of the Earth's dynamical flattening and other constants related to orbital motion of the Moon and the Earth (see for example Kinoshita 1977; Dehant & Capitaine 1997). The geodesic precession is a general relativistic effect related to the rotation of the geocentric reference system with respect to the solar-system barycentric reference system (see for example Brumberg 1991).

3.2 The Lieske et al. (1977) expressions for the quantities

Expressions for the precession quantities, both for the equator and the ecliptic, were provided by L77 to be in agreement with the 1976 System of Astronomical Constants. These expressions are based on the 1976 values for the planetary masses, precession constant and J2000 obliquity, as well as on Newcomb's expression for the motion of the ecliptic.

The value for the precession constant was an observationally determined value of Newcomb's precessional constant, or rather speed of general precession in longitude, the precise interpretation of which is not obvious (Fricke 1971; Lieske 1985). The value for the geodesic precession is that of de Sitter & Brouwer (1938), $p_{\rm g}=1\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{92}$/cy.

The following L77 expressions provide the precession of the equator, based on the different angles mentioned in the previous section, the equatorial precession angles provided by (7) being the most often used in practice:

 

$$\psi_{\rm A} = 5038\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em... ...t^2-0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{001147}~t^3$$

$$\omega_{\rm A} = \epsilon_0 + 0\hspace{0.1em}'\hspace{-0.06em}'\h... ...t^2-0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{007726}~t^3$$

$$\epsilon_{\rm A} = \epsilon_0 - 46\hspace{0.1em}'\hspace{-0.06em}... ...2 + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{001813}~t^3$$

$$\chi_{\rm A} = 10\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.... ...2-0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{001125}~t^3 ,$$ (6)

with $\epsilon_0 = 84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{448}$ for the obliquity at J2000.0, or alternatively:

 

$$\zeta_{\rm A} = 2306\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3e... ...^2+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{017998}~t^3$$

$$\theta_{\rm A} = 2004\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3... ...^2- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{041833}~t^3$$

$$z_{\rm A} = 2306\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\... ...+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{018203}~t^3 .$$ (7)

The parameter t, used in the above expressions as well as in those below, is the elapsed time in Julian centuries since J2000 TT, defined by:

$$t=({\rm TT}-2000 {{\rm\ January\ 1d\ 12h}\ {\rm TT})}/36~525 ,$$ (8)

with TT in days.

   
3.3 Limitations

The IAU 1976 precession model has the following limitations:

(i)

The expressions are based on numerical values for the precession rates which have been shown from VLBI observations to be in error by about -3 mas/year in longitude and -0.25 mas/year in obliquity. The theoretical basis for the precession rate in obliquity, neglected in previous theories, has been developed in W94.

(ii)

The expressions are based on the 1976 value for the obliquity at J2000.0 which is known, based on LLR and planetary observations, to be in error by $0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{04}$.

(iii)

It uses the motion of the ecliptic which is based on Newcomb's solution and on the 1976 values for the planetary masses, whereas more recent analytical solutions for the ecliptic and improved numerical values for the planetary masses have been available for at least a decade.

(iv)

The numerical values of the coefficients in the expressions are given with a resolution of 0.1 milliarcsecond after a century, whereas the amplitudes of the MHB 2000 nutation are provided with 0.1 microarcsecond resolution.

(v)

The numerical expression for the geodesic precession is limited to the secular term which itself has a precision of the order of  $0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{01}$/cy.

(vi)

The developments of the expressions are limited to the third degree in t.

The error in the precession of the equator resulting from (i) is of the order of $0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{3}$/cy in longitude and $0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{025}$/cy in obliquity; the error in the precession of the ecliptic resulting from (ii) and (iii) is of the order of $0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{002}$/cy, 80% of the effect coming from the use of an improved theory and 20% from the use of improved values for the planetary masses.

Moreover, it has been shown that the 3-rotation and 4-rotation transformations using the L77 quantities (expressions (6) and (7)) show disagreement at 1 mas level over 2 centuries. See Capitaine et al. 2003a for more detail.

   
4 The IAU 2000 precession

4.1 The MHB 2000 model

The IAU 2000A Precession-Nutation model was adopted by IAU 2000 Resolution B1.6 to replace the IAU 1976 Precession (L77) and the IAU 1980 Theory of Nutation. The nutation series was generated by the convolution of the MHB 2000 transfer function with the rigid-Earth nutation series REN-2000 (Souchay et al. 1999), rescaled to account for the change in the dynamical ellipticity of the Earth implied by the observed correction to the lunisolar precession of the equator. It is based upon basic Earth parameters estimated from VLBI observations.

The resulting nutation series includes 678 lunisolar terms and 687 planetary terms and provides the direction of the celestial pole in the GCRS with an observed accuracy of 0.2 mas. The series includes the geodesic nutation (Fukushima 1991). On the other hand, the Free Core Nutation (FCN), which cannot be predicted rigorously, is not included in the IAU 2000A model, and sets a fundamental "noise level'' of a fraction of 1 mas if IAU 2000A is used as it is.

It should be noted that the IAU 2000A nutation series include nutations with very long periods ( $\rm i.e.> 250$ years) the contribution of which can be approximated, in $\mu$as, as:

$${\rm d}\psi = -1146 -159~t~$$

$${\rm d}\epsilon = +1141 -30~t .$$ (9)

In former precession-nutation models these long-period effects were included in the precession part. In models compliant with MHB 2000, the presence of these terms in the nutation part causes compensating changes in the precession part.

The IAU 2000 nutation series is associated with improved numerical values for the precession rate of the equator in longitude and obliquity:

 

$$\delta\psi_{\rm A} = (-0.29965 \pm 0.00040) ''/{\rm cy}$$

$$\delta\omega_{\rm A} = (-0.02524 \pm 0.00010) ''/{\rm cy}.$$ (10)

   
4.2 The frame bias at J2000

The IAU 2000 precession-nutation is associated with the constant offsets $\delta\psi_0$ and $\delta\epsilon_0$ of the direction of the CIP at J2000.0 from the direction of the pole of the GCRS, which have been estimated from VLBI data (Herring et al. 2002):

 

$$\delta\psi_0 = (-0.0417750 \pm 0.000025)~''$$

$$\delta\epsilon_0 = (-0.0068192 \pm 0.0000100)~''.$$ (11)

In contrast, $d\alpha_0$, the offset in right ascension of the mean equatorial frame with respect to the GCRS, cannot be derived directly from VLBI observations, which are insensitive at the first order to the position of the ecliptic. The determination of this offset thus requires the use of observations which are dependent on the position of the ecliptic. To take proper account of the equinox offset, it is necessary to use an equinox that (i) corresponds to an ecliptic dynamically consistent with the IAU's adopted precession-nutation model and (ii) can be provided by high accuracy observations.

The numerical value that has been used for the implementation of the IAU 2000 precession-nutation model is the GCRS right ascension of the mean dynamical equinox at J2000 as provided by Chapront et al. (2002) from a fit to LLR observations based jointly on the use of a dynamical theory for the Moon and of VLBI Earth orientation parameters:

 

$$d\alpha_0 = (- 0.0146 \pm 0.0005)''.$$ (12)

This value (12) for the equinox offset fulfils condition (i) as the MHB precession and nutation is based on the theory of the Earth's rotation for a rigid Earth (Souchay et al. 1999) which uses the analytical theories for the planets and the Moon and then refers to the same dynamical equinox (S94). It fulfils condition (ii) as the accuracy of the estimation is at least 20 times better than the previously available ones.

Note that the mean equinox of epoch derived in this way corresponds to the definition of the ecliptic in its "inertial'' sense, i.e. based on the motion of the orbital angular-momentum vector of the Earth-Moon barycenter. It differs by about 94 mas from the "rotational dynamical mean equinox of J2000.0'' (Standish 1981) as used in the past when referring to the FK5 equinox or to the origin of the JPL ephemerides DE200 and the ICRS position of which has been provided with an uncertainty of 10 mas by Folkner et al. (1994).

In the following, we will first deal with the part of the MHB 2000 model due to precession. We will then consider the polynomial part of the motion of the Celestial Intermediate Pole (CIP) with respect to the GCRS, including precession and bias.

4.3 Basic IAU 2000 expressions for precession

The most straightforward way to interpret the precession part of IAU 2000A is simply to add the longitude and obliquity rates to the existing L77 series for $\psi _{\rm A}$ and $\omega _{\rm A}$. The canonical formulation for the IAU 2000 precession matrix is thus naturally the 4-rotation one (see relations (4)), involving  $\epsilon_0$, $\psi _{\rm A}$, $\omega _{\rm A}$ and  $\chi_{\rm A}$.

The IAU 2000 expressions for the classical quantities $\psi _{\rm A}$ and  $\omega _{\rm A}$ have been provided in the IERS Conventions 2000 and should be regarded as the "defining model'', upon which the IERS and SOFA implementations are based:

 

$$\psi_{\rm A} = 5038\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em... ...t^2-0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{001147}~t^3$$

$$\omega_{\rm A} = \epsilon_0 - 0\hspace{0.1em}'\hspace{-0.06em}'\h... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{007726}~t^3 ,$$ (13)

with $\epsilon_0 = 84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{448}$.

Following such an approach, the above MHB 2000 precession corrections have also to be used to correct the linear terms (noted by index 1) of the quantities $\epsilon_{\rm A}$ (obliquity of the equator on the moving ecliptic) and $p_{\rm A}$ (general precession in longitude), whereas the expression for  $\chi_{\rm A}$ is unchanged. By adding the corrections ${\rm d}\omega_1$ to  $\epsilon_{\rm A1}$ and d$\psi_1$to $p_{\rm A1}$, respectively, one obtains:

 

$$p_{\rm A} = 5028\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\... ...2 - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{000006}~t^3$$

$$\epsilon_{\rm A} = \epsilon_0 - 46\hspace{0.1em}'\hspace{-0.06em}... ...2 + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{001813}~t^3$$

$$\chi_{\rm A} = 10\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{001125}~t^3 .$$ (14)

4.4 Expression for the IAU 2000A precession based on the X, Y coordinates of the CIP

The coordinates X and Y of the CIP in the GCRS have been developed (Capitaine et al. 2003a) to be consistent with the IAU 2000A nutation series and the basic expressions (13) for IAU 2000 precession. They used expression (12) for the equinox offset at J2000.0 and the celestial pole offsets at J2000.0, $\xi_0$ and $\eta_0$, which are derived from (11) as:

 

$$\xi_0 (-0.016617 \pm 0.000010)''$$

$$\eta_0 (-0.006819 \pm 0.000010)''.$$ (15)

The IAU 2000A X and Y expressions, provided in the IERS Conventions 2000, have the following form:

 

$$X_{\rm IAU2000} - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4272191}~t^2$$ =

$$-~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00004605}~t^4$$

$$+~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ...sum_{i} \sum_{j=0}^{2}\big[(a_{{\rm s},j})_i t^j \sin ({\rm\scriptstyle {ARG}})$$

$$+ (a_{{\rm c},j})_i t^j \cos ({\rm\scriptstyle {ARG}})\big] + \cdots,$$ (16)

 

$$Y_{\rm IAU2000} - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ...- 22\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4072510}~t^2$$ =

$$+~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ...+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00111306}~t^4$$

$$+~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ...\sum_{i} \sum_{j=0}^{2}\big[(b_{{\rm c},j})_i t^j\cos ({\rm\scriptstyle {ARG}})$$

$$+ (b_{{\rm s},j})_i t^j \sin ({\rm\scriptstyle {ARG}})\big] + \cdots$$ (17)

ARG stands for various combinations of the fundamental arguments of the nutation theory, including both luni-solar and planetary terms.

Taking into account precession only, these expressions, denoted P00, become:

 

$$X_{\rm P00} - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4272191}~t^2$$ =

$$-~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00004605}~t^4$$

$$+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00000598}~t^5,$$ (18)

$$Y_{\rm P00} -0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{... ...- 22\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4072510}~t^2$$ =

$$+~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ...+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00111306}~t^4$$

$$+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00000099}~t^5.$$ (19)

This shows that the polynomial part of the X and Y CIP coordinates originate from precession, except for the contribution from the frame bias and from cross nutation terms. The latter contribute as a constant term in Y of 132 $\mu$as and a secular term in X of 4 $\mu$as per century.

The contributions d $X_{\rm bias}$ and d $Y_{\rm bias}$ from the frame biases in X and Y are (in $\mu$as):

$${\rm d}X_{\rm bias} = -~16617\ - 2\ t^2 + 1\ \cos\Omega,$$

$${\rm d}Y_{\rm bias} = -~6819\ -142\ t + 1\ \sin\Omega,$$ (20)

the first term in each coordinate being the contribution from the celestial pole offset at J2000 and the following ones from the frame bias in right ascension. As these contributions are included in the Earth Orientation Parameters which are derived from VLBI observations, it appears preferable to consider the complete polynomial part of the expressions without separating precession, nutation and biases.

4.5 IAU 2000 expressions for the quantities z ${_{\rm A}}$, ${\zeta _{\rm A}}$, ${\theta _{\rm A}}$

The computation of improved expressions for the currently used precession quantities $z_{\rm A}$, ${\zeta _{\rm A}}$ and ${\theta _{\rm A}}$ can be derived in several ways from expressions (13) and (14) based on analytical or numerical solutions. The purpose is to provide expressions which are numerically equivalent to the 4-rotation expressions, based on the MHB corrections to the precession rates in $\psi$ and $\omega$, at one level of accuracy less than the MHB corrections themselves (i.e. $30~\mu$as/cy). This can be done in a semi-analytical way by solving the expressions for these quantities from the relations in the spherical triangles (Woolard & Clemence 1966) using, as basic developments, those for  $\psi _{\rm A}$ and  $\omega _{\rm A}$; this was carried out using the software GREGOIRE as a tool. The transformation can also be done numerically, by sampling the precession matrix throughout a chosen interval, decomposing the matrix into the three Euler angles concerned, and fitting polynomials in t. The results are quite consistent at the microarcsecond level over several centuries. Our approach was to use just enough terms, and coefficients with just enough precision, to match the resolution of the IAU 2000A model. The order of polynomial to be used was found to be t⁵ and the precision of the coefficients $0.1~\mu$as. The following series with a $0.1~\mu$as level of precision matches the canonical 4-rotation series to sub-microarcsecond accuracy over 4 centuries:

 

$$\zeta_{\rm A} = 2\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.... ... + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{3019015}~t^2$$

$$\qquad + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000002}~t^5$$

$$z_{\rm A} = - 2\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\h... ... + 1\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0947790}~t^2$$

$$\qquad + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000003}~t^5$$

$$\theta_{\rm A} = 2004\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0418251}~t^3$$

$$\qquad - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000001}~t^5.$$ (21)

These expressions have therefore been proposed as quantities of precession consistent with the IAU 2000A model, and were provided in the IERS Conventions 2000.

5 Precession-nutation models post IAU 2000

5.1 Areas of possible improvement

The only corrections that have been applied in the above calculations in order to be consistent with the IAU 2000A model are the MHB corrections to precession rates in longitude and obliquity. The value used for  $\epsilon_0$, as well as the expressions used for the motion of the ecliptic and for the other quantities of precession, were those of L77 (i.e. IAU 1976).

Such an approach is not satisfactory from a dynamical point of view because the precession rate corrections are understood only as a representation of linear terms in the observables and not as physical quantities. Moreover, such a precession model suffers, except for the precession rates, from the limitations of the IAU 1976 model mentioned in Sect. 3.3 and especially regarding the model for the precession of the ecliptic.

As noted by Woolard & Clemence (1966), the motion of the equator and the motion of the ecliptic are kinematically independent of each other, but to a small extent the secular variation in the motion of the equator depends dynamically upon the variations of the disturbing forces caused by the change in the average positions of the Sun and the Moon with the motion of the ecliptic. The improvement of the model for the precession of the equator therefore requires the use of an improved model for the motion of the ecliptic.

Based on the MHB corrections to the precession rates, the precession of the equator can be obtained in a more dynamically consistent way. The solution has to be based on:

(i)

improved expressions for the motion of the ecliptic,

(ii)

improved expressions for the contributions to precession and obliquity rates of the equator with respect to a fixed frame such as those provided in W94,

(iii)

the MHB estimates for the dynamical ellipticity of the Earth and for the effect of non-rigidity.

A more accurate expression for the geodesic precession is also required (Brumberg 1991) in order that it be consistent with the same level of accuracy as the other precession expressions. The following points should be noted:

1) Such further improvement, based on the MHB 2000 corrections to precession rates, complies with the recommendations of IAU 2000 Resolution B1.6, which encourages the development of new expressions for precession consistent with the IAU 2000A model.

2) An improved GMST expression would also be necessary as the expression directly depends on the precession in right ascension. On the other hand, neither the periodic part of the expression for the position of the CEO in the GCRS, nor that for the complementary terms in the equation of the equinoxes, will have to be modified, as the expected corrections to precession are below the sensitivity of the terms of these expressions to the precession quantities.

3) Regarding expressions (14), the formula for $p_{\rm A}$is no longer appropriate as the primary expression for precession, and the expression for  $\chi_{\rm A}$ is useful only as an intermediate in the coordinate transformation between TRS and CRS, as its effect is canceled out in the GST expression.

5.2 Method for improving the model for the precession of the equator

In contrast to the previous approach, a physical meaning can be given to the corrections to precession rates. The MHB correction d$\psi_1$ is then considered to be a correction to the speed of the mean celestial pole of date, P, about the ecliptic pole C (see Woolard & Clemence 1966 or L77) and d$\omega_1$ as the obliquity rate contribution with respect to inertial space (see W94). This provides improved constants for the precession model of the equator to be used in the kinematical relations.

This approach is more complicated, but more satisfactory than the previous one from a theoretical point of view. Its use does not mean that the observations are considered to be sensitive to an ecliptic. It means only that VLBI estimated corrections (d$\psi_1$and d$\omega_1$) to the precession rates with respect to a fixed ecliptic are used as constants of integration in solving the equations, thus providing corrections to the coefficients of the t² and t³ terms of the precession quantities.

There is however a problem of numerical consistency. The t²and t³ corrections that appear, arising from the MHB correction to the speed of the celestial pole, are indeed largely less than the corrections which would correspond to a better representation of the motion of the ecliptic. It would therefore be artificial to introduce corrections at the microarcsecond level when it is well known that corrections more than ten (or a hundred) times larger would be necessary. Moreover, the expressions of L77 are limited to terms of degree 3 or less, whereas there are non-negligible terms at degree 4, and the computations are therefore not consistent in principle. A development up to degree 4 is necessary.

5.3 Method for improving the model for the precession of the ecliptic

Since the adoption of the IAU 1976 precession, polynomial updates have been provided using improved ecliptic motion by Bretagnon & Chapront (1981), Laskar (1986) and S94. The latter paper provided expressions derived from improved precesssion rate, obliquity and also masses, and partial derivatives with respect to these quantities. These developments were used by W94, together with the newly introduced precession rate in obliquity and updated VLBI precession rate in longitude, to provide the most complete polynomials in time for the precession quantities. These latest expressions for the precession quantities are based on the analytical theory VSOP87 (Bretagnon & Francou 1988) for the motion of the ecliptic, which was fitted to the numerical ephemerides DE200/LE200, and on the IERS Standards 1992 system of planetary masses (McCarthy 1992). Improved expressions for the precession of the ecliptic can be provided by an appropriate fit of the semi-analytical solution of the motion of the ecliptic in the ICRS to the most accurate available numerical integration.

Table 1: Precessional quantities determined with different sources. Units: coefficients in arcsecond and time in century.

6 Improving the models for precession

6.1 Precession of the ecliptic

   
6.1.1 The precession quantities for the ecliptic

Hereafter we designate by P and Q the quantities $P= \sin\pi\sin\Pi$ and $Q=\sin\pi \cos\Pi$, where $\pi$ and $\Pi$are the osculating elements of the Earth-Moon barycenter orbit ($\pi$, the inclination on the ecliptic and $\Pi$, the longitude of the ascending node). The angles $\pi$ and $\Pi$ are referred to a fixed ecliptic for J2000.0. The planetary theory VSOP87 introduces slightly different quantities, $p=(\sin \pi/2) \sin \Pi$ and $q=(\sin \pi/2) \cos \Pi$. The corresponding precession quantities are denoted by $P_{\rm A}=\sin\pi_{\rm A}\sin\Pi_{\rm A}$ and $Q_{\rm A}=\sin\pi_{\rm A}\cos\Pi_{\rm A}$. They are time polynomials, easily derived from the secular developments of p and q in the Earth-Moon barycenter orbital motion. Such developments have been provided by S94 and W94.

The quantities p and q in VSOP87 are quasi-periodic functions of the time. They are expressed in the form of Poisson series whose arguments are linear combinations of the mean planetary longitudes; more precisely they involve a secular part, a periodic part (Fourier terms) and a mixed part (Poisson terms, i.e. Fourier terms with varying amplitudes). The secular parts provide $P_{\rm A}$ and $Q_{\rm A}$, which may be regarded as, respectively, the x and -y components of the secularly-moving ecliptic pole vector in a (right-handed) frame that has its x-axis through the J2000 (inertial) mean equinox and its z-axis through the J2000 ecliptic pole.

From the original expressions for the polynomial parts for p and q given in S94, which were computed using the IAU 1976 system of planetary masses, we have derived the polynomial developments for $P_{\rm A}$ and $Q_{\rm A}$. These expressions are very close to the developments used by W94 since the source is the same (VSOP87), the only change coming from the introduction of the IERS Standards 1992 masses (McCarthy 1992). The coefficients of these developments and of L77 are listed in Table 1 together with the improved values obtained from the procedures set out below.

When looking at the trends in the residuals between VSOP and various JPL source ephemerides (DE200, DE403 and DE405), the secular quantities ${\left(\frac{{\rm d}p}{{\rm d}t}\right)}_{t=0}$ and ${\left(\frac{{\rm d}q}{{\rm d}t}\right)}_{t=0}$ show systematic deviations that are probably due to the analytical solution VSOP87, independently of the reference frame, the constants of integrations and other physical parameters of the JPL reference ephemeris. The corresponding deviations in $P_{\rm A}$ and $Q_{\rm A}$ are, in the case of DE406: ${\left(\frac{{\rm d}P_{\rm A}}{{\rm d}t}\right)}_{t=0} \simeq -0.6$ mas/cy and ${\left(\frac{{\rm d}Q_{\rm A}}{{\rm d}t}\right)}_{t=0} \simeq -1.7$ mas/cy. The contributions due to the change of masses mentioned above are -0.1 mas/cy and -0.3 mas/cy respectively and are much smaller. An illustration of the periodic and secular deviations of P and Q between VSOP87 and DE406 is shown in Figs. 2 and 3, the time interval covering two millennia. The thickness of the curves brings out the residuals due to short-period terms whose amplitudes are smaller than 1 mas. The "noise'' produced by the short-period terms is the main limitation to the improvement of the secular variations described below.

6.1.2 Improving the developments

Our "improved'' ecliptic model comes from two independent descriptions of the path of the (inertial) ecliptic pole, one from VSOP87 and the other from DE406. VSOP87 provides the most detailed analytical model currently available of the periodic component of the ecliptic motion. DE406 contributes the best available accuracy, consistency with modern observations, and a long time span.

VSOP87 provides an ecliptic in the form of the osculating elements p,q, which are obtained from models that comprise a secular part (that we wish to improve upon) and a periodic part. The DE406 ecliptic is obtained from the Earth Moon Barycenter (EMB) position and velocity $(\mbox{\boldmath$r$ },\mbox{\boldmath$\dot{r}$ })$: the (inertial) ecliptic pole lies in the direction of r $\times$ $\dot{r}$ and its motion contains both secular and periodic effects. Combining these two ecliptics to obtain an improved secular model consists of finding the best rotation from the DE406 (=DE405) frame to ecliptic coordinates, then comparing the DE406 and VSOP87 p,q predictions, and finally minimizing the differences by making small adjustments to the VSOP87 secular models.

	Figure 1: Rotation angles $\phi$ and $\epsilon$ to rotate from equatorial to ecliptic frame.
Open with DEXTER

The procedure has one unavoidable weakness in that the final polynomial models for the secular part of p, q (and hence $P_{\rm A}$, $Q_{\rm A}$) are essentially empirical and have no regard for dynamical consistency. Although the VSOP87 model is theory-based, the DE406 predictions are, as far as this exercise is concerned, a source of observations, and the resulting polynomials empirically model those observations. However, they can be expected to be valid over the entire time span used (in our case 1000-3000 AD), with any unmodeled long-period terms in VSOP87 automatically corrected.

The details of the procedure we used to provide an improved ecliptic model were as follows:

1.

Adopt starting values for $\phi_0$, $\epsilon_0$ and $\psi_0$, the three Euler angles that relate the orientations of DE405/DE406 and VSOP87. We used $84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4}$ for $\epsilon_0$ and zero for the other two angles.

2.

Initialize to zero the current corrections to the VSOP87 polynomials. To match VSOP87, the two correction polynomials each consist of five coefficients, for powers of t from 1 to 5. (NB: The VSOP87 polynomials require t in millennia, whereas in this paper we use t in centuries.)

3.

(The procedure iterates from this point.) Using the current values for $\phi_0$ and  $\epsilon_0$, calculate the third Euler angle $\psi_0$ (see Fig. 1):

 

$$\begin{array}{rcl}\psi_0 & = & (\phi_0-\phi({\rm DE405})+\phi({\r... ...os{\epsilon_0} \\ & & +~\xi_0/ \sin{\epsilon_0}-\psi({\rm DE405}), \end{array}$$ (22)

where $\xi_0$ is 16.617 mas, the ICRS declination of the J2000 mean equinox, $\psi({\rm DE405})$ is 6.4 mas, the distance from the ICRS equinox to DE405 equinox, $\phi({\rm DE405})$ is 50.28 mas, the distance from the DE405 origin to the DE405 equinox, $\phi({\rm ICRS})$ is 55.42 mas, the distance from the ICRS origin to the ICRS equinox, and $d\alpha_0$ is 14.6 mas, minus the ICRS right ascension of the mean J2000 equinox. $\xi_0$is from the IAU 2000 frame bias (11) as presented in the IERS Conventions 2000, $\psi({\rm DE405})$, $\phi({\rm DE405})$, $\phi({\rm ICRS})$ and $d\alpha_0$ are from Chapront et al. (2002), the $d\alpha_0$ value being the adopted value (12). It should be noted that these numbers are not entirely consistent, because they come from observations. However, expression (22) defines how they have been combined to produce a working $\psi_0$ value.

4.

For a set of regularly-spaced times t (in Julian centuries of TDB) centered on J2000, interrogate DE406 to obtain the EMB heliocentric position and velocity and VSOP87 to obtain p and q. Our samples were spaced 0.7 years from J1000 to J3000, giving a t range of $\pm 10$.

5.

Transform the DE406 heliocentric position and velocity from the DE405 (i.e. equatorial) frame to the J2000 ecliptic frame by applying the rotation $R_3(\psi_0)R_1(\epsilon_0)R_3(-\phi_0)$ and, using the standard transformations, obtain the osculating elements $\pi$and $\Pi$.

6.

Compute $p=(\sin \pi/2) \sin \Pi$ and $q=(\sin \pi/2) \cos \Pi$ and subtract the VSOP87 values from Step 4 to give $\Delta p$ and $\Delta q$. The result of the last three steps is a table of t, $\Delta p$ and $\Delta q$.

7.

Using standard minimization techniques, fit polynomial coefficients in t_0-n to the sets of $t,\Delta p$ and $t,\Delta q$. As mentioned in Step 2, we chose n=5 to match the order of the VSOP87 p and q secular models.

8.

The t₀ terms of the resulting polynomials, p₀ and q₀, have a special status in that they describe the offsets at J2000 between the DE406 ecliptic and the one used to rotate the DE406 vectors into the ecliptic frame. They can thus be used to improve the $\epsilon_0$ and $\phi_0$ estimates, by adding 2q₀ to  $\epsilon_0$ and subtracting $2p_0/\sin{\epsilon_0}$ from $\phi_0$.

9.

The remaining polynomial terms, p_1-5 and q_1-5, are subtracted from the current corrections to the VSOP87 polynomials (see Step 2).

10.

Repeat from Step 3 until the $\Delta p$ and $\Delta q$polynomial coefficients have minimized. Convergence is rapid: no more than three iterations should be required.

11.

The final corrections are then subtracted from the VSOP87 polynomials, giving improved p,q models. These can be transformed into models for $P_{\rm A},Q_{\rm A}$ to give the final result.

The final DE405-to-ecliptic rotation angles were:

 

$$\phi_0 = 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{05132}$$

$$\epsilon_0 = 84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{40889}$$

$$\psi_0 = 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{03862}.$$ (23)

The improved $P_{\rm A}$ and $Q_{\rm A}$ series are as given as the last line (P03) of Table 1.

The number of quoted decimal places was chosen after observing the small variations seen when the time step was changed. Smaller steps than the adopted 0.7 years produced data sets so large that the accumulated rounding errors in the least-squares fitting procedure used began to influence the answers. Steps larger than 0.7 years degraded the accuracy through under-sampling. The acceptable margin was roughly a factor of two either side of the adopted time step. In the case of $P_{\rm A}$ and $Q_{\rm A}$, the terms in t⁵ were of marginal significance and have been omitted from Table 1.

Although the procedure set out above naturally produces an ecliptic according to the inertial definition it can readily be adapted to use the rotating definition instead (Standish 1981). The adaptation consists of expanding Step 5 to work in the rotating frame of the ecliptic of date rather than the fixed J2000 ecliptic frame. Because the required frame rotation and spin matrices are functions of the $P_{\rm A},Q_{\rm A}$series, these have to be determined from the currently estimated p,q series at each iteration, as well as at the end. We introduced the further refinement of editing the VSOP87 series for p,q to suppress annual terms, though this had only a small effect on the result. These adaptations produced an equinox displaced by -93.782 mas from the inertial equinox, to $\phi_{0r} = -0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{04246}$, and changed the obliquity by 3.329 mas from the inertial value, to $\epsilon_{0r} = 84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{41222}$.

	Figure 2: Differences P(DE406) - P(VSOP) in $P= \sin\pi\sin\Pi$.
Open with DEXTER

	Figure 3: Differences Q(DE406) - Q(VSOP) in $Q=\sin\pi \cos\Pi$.
Open with DEXTER

6.1.3 Accuracy estimates

It is useful to estimate the precision which can be assumed for the coefficients in the above determination. We have performed several tests to assess the influence of the fitting intervals on the determination of ${\left(\frac{{\rm d}P_{\rm A}}{{\rm d}t}\right)}_{t=0}$ and ${\left(\frac{{\rm d}Q_{\rm A}}{{\rm d}t}\right)}_{t=0}$. Over the time interval $\Delta t=$ [1900, 2100], we computed the linear regressions of the residuals: $P({\rm DE406})-P({\rm VSOP87})$ and $Q({\rm DE406})-Q({\rm VSOP87})$. The tests consisted of varying the set of abscissas, the length of the time interval $\Delta t$ modified by 10 percent and a shift of the center by $\pm$10 years. Our conclusion from these tests is that our determination of the secular variations of $P_{\rm A}$ and $Q_{\rm A}$ is better than 0.05 mas/cy close to J2000.0 although it is given with more digits in Table 1 for the purpose of internal computations.

The secular parts of the variables p and q are very sensitive to the long-period terms which have been retained in the analytical solution and to their accuracy. Table 2 shows the three terms of longest period which exist in VSOP87 in the variables pand q. Here, the general formulation for a Poisson term is $At^\alpha\cos (\omega t+\phi)$, where A is the amplitude, $\alpha$ the associated power of time, $\phi$ is the phase and $\omega$ the frequency, the time being reckoned from J2000.0. The period $P=2\pi/\omega$ is given in years. If we choose to develop with respect to the time the long-period term corresponding to the argument 4T-8Ma+3J and the very-long-period Poisson term corresponding to $T-D-L_{\rm Moon}$, the constant term in q and its trend ${\left(\frac{{\rm d}q}{{\rm d}t}\right)}_{t=0}$ will be slightly modified in the original solution (p is unchanged). The fit with DE406 compensates this discrepancy if it is performed on a sufficiently long time interval. On the other hand, a development of the term corresponding to 2J-5S will be not appropriate after few centuries.

Table 2: The long-period terms in VSOP87: $At^\alpha\cos (\omega t+\phi)$. The planetary and lunar arguments are $T =\rm Earth$, $Ma =\rm Mars$, $J = \rm Jupiter$, $S = \rm Saturn$, D = Delaunay's argument D (the difference of the mean longitudes of Moon and Sun), $L_{\rm Moon} =\rm mean$ longitude of the Moon.

In Sect. 4.2, we discussed the rotation between GCRS and the mean J2000 equator and equinox in the IAU 2000 model. The pole offset $(\delta\psi_0,\delta\epsilon_0)$ was defined by the MHB 2000 precession-nutation model itself, but it was necessary to position the equinox by choosing a value for $d\alpha_0$, thereby completing the 3D frame rotation.

The adopted value $d\alpha_0 = -14.6$ mas of the ICRS RA of the mean dynamical equinox at J2000 was that estimated by Chapront et al. (2002) from LLR data and VLBI EOP and therefore benefited from both LLR and VLBI data for placing the ecliptic in the ICRS. With an improved precession model in prospect, the question naturally arises or whether to retain this IAU 2000 frame bias or to consider introducing an improved frame bias.

As our goal is to remain compliant with the IAU 2000 resolutions as far as possible, it seems clear that the location of the J2000 mean pole must remain as defined by MHB 2000. However, if we look critically at $d\alpha_0$ we find that choices do in fact emerge, because of small inconsistencies between the LLR results and MHB 2000. In particular, in the Chapront et al. (2002) LLR determination, the mean equator at J2000 is based on the IERS 1996 precession-nutation model rather than MHB 2000 (which was not yet available at the time). This means that the LLR and MHB estimates for the frame bias of the J2000 mean equator with respect to the ICRS equator are slightly different, about 17.7 mas and 16.6 mas respectively. In view of this, it can be argued that an a posteriori adjustment to the published LLR $d\alpha_0$ would be justified. This would give a $d\alpha_0$ figure of -17.1 mas, a difference of -2.5 mas with respect to the value adopted for the IAU 2000 model. It should however be noted that:

1. The difference in the CIP motion resulting from a $d\alpha_0$change of 2.5 mas would be only one term of $24~\mu$as/cy in Y, which is not of practical significance.

2. The $d\alpha_0$ value cannot be evaluated to better than about $\pm$3 mas, given the uncertainties in the ICRS position both of the MHB J2000 mean equator and of the ecliptic.

Our conclusion is that although a revised $d\alpha_0$ could be introduced as a complement to an improved precession model, the inconvenience and opportunities for error that multiple frame bias models would lead to would not be repaid by worthwhile improvements in accuracy. We therefore recommend retaining the adopted value of -14.6 mas as a conventional value.

6.2 Precession of the equator

6.2.1 Equations for the precession quantities

 

The precession of the equator with respect to a fixed frame can be derived by solving the differential equations for the precession quantities $\psi _{\rm A}$, and $\omega _{\rm A}$, given the expressions for the various contributions to the precession rates (that are dependent on the orientation of the ecliptic with respect to the equator) and the precession of the ecliptic. This would theoretically require using expressions for these contributions referred to a fixed equator and ecliptic as was done by Woolard (1953) to provide the precession-nutation solution directly referred to the mean ecliptic frame of the epoch 1900. However, the most complete recent expressions for these contributions (for example in W94) are for the components $r_{\psi }$ and $r_{\epsilon }$ of the precession rates in longitude and obliquity respectively, expressed in an equatorial frame linked to the moving equinox. These components have therefore to be rotated using the quantity $\chi_{\rm A}$ for the planetary precession. The equations to be used are the first equations of formula (14) of L77 or (29) of W94, the last one including the additional component $r_{\epsilon }$, which was not present in the developments of L77.

The differential equations to be solved for providing the two basic quantities for the precession of the equator are then written as:

 

$$\sin\omega_{\rm A} {\rm d}\psi_{\rm A} /{\rm d}t = (r_\psi \sin\epsilon_{\rm A}) \cos\chi_{\rm A} - r_\epsilon \sin\chi_{\rm A} ,$$

$${\rm d}\omega_{\rm A} /{\rm d}t = r_\epsilon \cos\chi_{\rm A} + (r_\psi \sin\epsilon_{\rm A}) \sin\chi_{\rm A} .$$ (24)

These equations, which are relative to the angular momentum axis, ignore the second order terms appearing in the differential equations for body axes which give rise to additional terms, the so-called Oppolzer terms (see Woolard 1953). As the precession rates estimated from observations are relative to the CIP, these complementary terms in precession rates have first to be subtracted from the values estimated from observations before these values are used as the constants of integration $(r_\psi)_{t=0}$ and $(r_\epsilon)_{t=0}$ for solving Eqs. (24), and then have to be added to the solutions to provide the precession relative to the CIP.

The complementary precession terms for the CIP axis can be derived from classical expressions of Oppolzer terms (see for example Capitaine et al. 1985 or Hartmann et al. 1999) in the following form:

 

$$\delta \omega_{\rm A} = (A/(C\Omega)) ({\rm d}\psi_{\rm A} /{\rm d}t) \sin\epsilon_0$$

$$\delta \psi_{\rm A} = -\lbrack A/C\Omega\sin\epsilon_0\rbrack ({\rm d}\omega_{\rm A} /{\rm d}t).$$ (25)

A and C being the Earth's principal moments of inertia and $\Omega$ the mean Earth's angular velocity, ( $2\pi \times 36~525$)/cy.

Using solutions (37) for $\psi _{\rm A}$ and $\omega _{\rm A}$, and the 0.996726205 MHB value for the ratio A/C in expression (25), shows that the only corresponding Oppolzer terms with amplitudes larger than $1~\mu$as for a century are (i) a bias in obliquity of $8704~\mu$as, that is in fact included in the adopted value for the mean obliquity at epoch, and (ii) secular terms of $-4~\mu$as/cy and $+1~\mu$as/cy in obliquity and longitude, respectively. It can be shown that the solution of Eqs. (24) is insensitive, within $1~\mu$as precision, to contributions to the integration constants smaller than 1 mas, and consequently the secular components of Oppolzer terms can be ignored.

The differential Eqs. (24) require expressions for the precession quantities $\epsilon_{\rm A}$ and $\chi_{\rm A}$, which depend both on the motion of the equator and the motion of the ecliptic, as does the quantity $p_{\rm A}$, which is also required for computing $\chi_{\rm A}$.

The differential equations for $p_{\rm A}$ and $\epsilon_{\rm A}$ can be derived from the expression for $r_{\psi }$ and $r_{\epsilon }$ respectively and the ecliptic precession (see Laskar 1986; Simon et al. 1997; or W94 for the equation in $\epsilon_{\rm A}$), which gives:

 

$${\rm d}p_{\rm A} /{\rm d}t = r_\psi -\cot\epsilon_{\rm A} \lbrack A(p,q) \sin p_{\rm A}$$

$$~~~\qquad\quad+ B(p,q)\cos p_{\rm A} \rbrack- 2C(p,q)$$

$${\rm d}\epsilon_{\rm A} /{\rm d}t = r_\epsilon - B(p,q) \sin p_{\rm A} + A(p,q) \cos p_{\rm A} ,$$ (26)

p and q being the polynomial parts of the ecliptic variables defined in Sect. 6.1.1 and A(p,q), B(p,q) and C(p,q), functions of these quantities such that:

$$r \left[ \dot{q} + p(q \dot{p} - p \dot{q}) \right] ,$$               A(p,q) =

$$r \left[ \dot{p} -q(q \dot{p} - p \dot{q}), \right]$$ B(p,q) =

$$q \dot{p} - p \dot{q}$$ C(p,q) = (27)

with $r=2/ \sqrt{1-p^2- q^2}$.

The quantity $\chi_{\rm A}$ can be derived from geometrical considerations in the spherical triangle N $\gamma_1 \gamma$, where N is the node of the moving ecliptic in the ecliptic at J2000 and $\gamma_1$ and $\gamma$ are the mean equinoxes of date and of J2000 respectively. Using the definition of the quantities $P_{\rm A}$ and $Q_{\rm A}$ as functions of elements  $\Pi_{\rm A}$ and  $\pi_{\rm A}$ of the spherical triangle, this gives:

 

$$\sin\chi_{\rm A} \sin\omega_{\rm A} = P_{\rm A} \cos p_{\rm A} + Q_{\rm A} \sin p_{\rm A} .$$ (28)

A simultaneous solution of the four differential Eqs. (24) and (26) together with the geometrical relation (28) is necessary and an iterative procedure is required for providing a solution for the precession of the equator that is dynamically consistent. The solution includes the expressions for the basic quantities  $\psi _{\rm A}$, $\omega _{\rm A}$ and for the secondary quantities  $\epsilon_{\rm A}$, $\chi_{\rm A}$ and $p_{\rm A}$, once the precession of the ecliptic is given.

   
6.2.2 Theoretical contributions to the precession rates

The forms of the solution of the differential Eqs. (24) and (26) for the precession of the equator depend on the expressions for the precession rates $r_{\psi }$ and $r_{\epsilon }$, that can be developed as polynomials in t:

 

$$r_{\psi}=r_0+r_1t + r_2 t^2 + r_3 t^3 ,$$

$$r_{\epsilon}= u_0 + u_1t + u_2 t^2 + u_3 t^3 .$$ (29)

The constants of integration for solving Eqs. (24) and (26) are the constant terms r₀ and u₀ in (29).

The theoretical contributions to the precession rate in longitude have been described in detail for a rigid Earth for example by Kinoshita (1977), Laskar (1986), W94, Souchay & Kinoshita (1996) and Roosbeek & Dehant (1998).

A complete evaluation of the theoretical contributions to the precession rates for a non-rigid Earth has been provided in W94 including, for the first time, the obliquity rate with respect to space and the tilt effects. This provided constant and time-dependent components in longitude and obliquity, corresponding to (i) first order, second order and J₄ effects in the direct luni-solar torque, (ii) direct planetary torque effect, (iii) J₂ and planetary tilt effects, (iv) tides and $\dot{J_2}$ effects and (v) geodesic precession.

The contributions (i) to (iv) result from the physical effect of the luni-solar and planetary torques on the oblate Earth, for which the dynamical ellipticity of the Earth, $H_{\rm d}=\lbrack C-(A+B)/2 \rbrack / C$ (A, B and C being the Earth's principal moments of inertia), is an essential parameter.

In contrast, the geodesic precession (v) originates from the relativistic rotation of the "dynamically non-rotating'' geocentric frame, in which the precession equations are solved for the dynamical effects (i) to (iv), with respect to the "kinematically non-rotating'' GCRS in which precession-nutation is actually observed. Consequently $H_{\rm d}$ is a scaling factor for all the effects, except for the geodesic precession. The initial numerical values for the components (except for the geodesic precession) are also dependent on the value for the mean obliquity of the ecliptic at J2000.

The computations performed in the present work for solving the precession equations make use of:

1.

MHB values (10) or equivalently (13) for the precession for computing the constants of integration r₀ and u₀ (see Sect. 6.2.4 for the details of this computation);

2.

the LLR estimated value, $\epsilon_0 = 84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{406}$ (see Fig. 1 and Sect. 6.1.1) for the mean obliquity at J2000, that is also the IERS 2000 value;

3.

W94 values described above for the various contributions (i) to (iv) to precession rates (Table 4 of W94 and Williams 1995), rescaled by the MHB values for the dynamical ellipticity, H$_{\rm d}$ of the Earth and for the mean obliquity of the ecliptic $\epsilon_0$, adopted here;

4.

the $\epsilon$-dependence of each component as provided in Table 4 of W94;

5.

MHB "non-rigidity contribution'' to the precession rate in longitude due to non-linear terms in the torque equations, which is, in microarcseconds (Mathews 2002):

$${\rm d}_{\rm nr}\psi=-21050~t;$$ (30)

6.

a complete expression for the geodesic precession from Brumberg (1991), Bretagnon et al. (1997) and Brumberg (2003), such as, in microarcseconds:

 

$$\psi_{\rm g} = 1919882.7~t - 503.9~t^2 - 0.7~t^3$$

$$\omega_{\rm g} = 1.0~t + 19.5~t^2 -4.7~t^3.$$ (31)

Table 3 provides the theoretical values used in our work for the different precession rate contributions, the origin of which is indicated in the first column. The form of the dependence of each effect on the obliquity of the ecliptic, $\epsilon$, that has to be considered when integrating the equations, is provided in the second column.

The effect of geodesic precession provided in this table has been obtained by converting expression (31) for the effect in the quantities  $\psi _{\rm A}$ and $\omega _{\rm A}$ into contributions $r_{\psi_{\rm g}}$ (usually denoted $P_{\rm g}$) and $r_{\epsilon{_{\rm g}}}$ to the precession rates $r_{\psi }$ and  $r_{\epsilon }$ respectively, these quantities being related through expression (24).

The main contributions to the precession rates that originate from the first order terms of the luni-solar torque are denoted (r₀)₁ and (u₀)₁ respectively. The (r₀)₁ term is the only one with a sufficiently large amplitude (of the order of $5000''/\rm cy$) to be sensitive to small changes in the value for the dynamical ellipticity $H_{\rm d}$of the Earth, which is one of the Basic Earth Parameters to be determined. The (r₀)₁ term is therefore of very special significance for computing the integration constant r₀ to be used for solving the precession equations. Its value is derived from the adopted value for r₀, given the other contributions which are provided by the theory.

It should be noted that the quantity (r₀)₁ is related to the so-called "Newcomb's precessional constant'' at epoch, ${\cal P}_0$, by $(r_0)_1={\cal P}_0\cos\epsilon_0$, whereas the quantity r₀ corresponds to the value at epoch of the so-called "speed of luni-solar precession'' in previous IAU models in which the planetary contribution to the precession of the equator was ignored. More generally, the largest contributions to $r_{\psi }$ provided in Table 3 correspond to effects considered in the expression given by Kinoshita (1977) for the speed of luni-solar precession, f₂₀₀₀, whereas the other contributions correspond to the additional effects mentioned above.

Table 3: Theoretical contributions (from W94, Williams (1995) and MHB) to the precession rates, $r_{\psi }$ and $r_{\epsilon }$, of the equator used in the present paper.

   
6.2.3 Spurious contributions to the estimated rates

The MHB precession rate corrections were estimated in a fit of Basic Earth Parameters (BEP), that are part of the precession-nutation theory, to VLBI series of celestial pole offsets. The estimated precession rates are therefore dependent (i) on the parameters that were used as observations, (ii) on the procedure that was used for estimating these observed parameters and (iii) on the MHB fit. The contributions of these different effects have been evaluated in order to be subtracted from the MHB estimates before using these values as constants of integration.

(i) Effect of the "observed'' parameters

The nutation parameters determined from the analysis of the VLBI data are corrections to the precession-nutation model expressed as celestial pole offsets in the form d$\psi$ and d$\epsilon$. It should be noted that these celestial pole offsets were determined using the conventional L77 value for the obliquity of the ecliptic at J2000 to which VLBI observations are actually not sensitive. This means that, whereas the VLBI estimate for the precession rate is provided for the correction d$\psi_1$ to the t term $\psi_1$ of the precession in longitude, the estimate in fact includes the effect of the conventional scaling factor $\sin\epsilon_{0}$through the precession nutation matrix product that was used in the VLBI computation. This is equivalent to saying that the parameter to which VLBI observations are actually sensitive is not d $\psi _{\rm A}$ itself, but is d $\psi_{\rm A} \sin\epsilon_{A}$ (i.e. the GCRS X-coordinate of the CIP). Consequently, a change in the conventional value for the obliquity of the ecliptic at J2000 would result in a spurious change in the corresponding estimate for $\psi_1$. The spurious effect in the precession rate resulting from a change d $\epsilon_0$ from the L77 value $(84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{448})$ to the IERS 2000 value $(84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{406})$for  $\epsilon_0$, is:

 

$${\rm d}_1\psi_1 = ~ \psi_1\ {\rm d}\epsilon_0 \cot \epsilon_{0} = ~ -2366~\mu{\rm as}/\rm cy .$$ (32)

(ii) Effect of the pre-2003 VLBI procedures

As mentioned in the introduction, the pre-2003 VLBI procedures did not use the rigorous transformation as described in Sect. 2, but used a procedure which a) considered the precession and frame biases corrections as if they were nutations and b) omitted the equinox offset.

This, in particular, introduces the following spurious contributions to the estimated precession rates:

 

$${\rm d}_2\psi_1 = \eta_0 \psi_1 \cot\epsilon_0 = - 384~\mu{\rm as}/\rm cy$$

$${\rm d}_2\omega_1 = - \xi_0 \psi_1 \cos\epsilon_0 - {\rm d}\alpha_0 \psi_1 \sin\epsilon_0= +514~\mu{\rm as}/\rm cy.$$ (33)

(iii) Effect of the MHB fit

The MHB procedure to estimate the precession rates $\lbrack (r_{\psi})_{t=0}\rbrack _{\rm MHB}$ in longitude and $\lbrack(r_{\epsilon})_{t=0} \rbrack _{\rm MHB}$ in obliquity included three steps (Mathews 2002). The first step was to estimate VLBI corrections to the luni-solar precession rate in longitude, P, to a number of prograde and retrograde nutation amplitudes and to the precession rate in obliquity. The second step was to fit the theoretical expressions for P and the nutation amplitudes to their observational estimates in order to get the dynamical ellipticity $H_{\rm d}$ as one of the BEPs. The last step was to use the estimated value for $H_{\rm d}$ to derive the MHB precession rate in longitude as being the sum of the value, $P_{\rm r}(H_{\rm d})$, corresponding to a rigid Earth and the other contributions to the precession rate (see Sect. 6.2.2), including that from non-rigidity. The value used for $({\rm d}P_R/{\rm d}H_{{\rm d}R})$ was 1 539 706 mas/yr from SMART97 (Bretagnon et al. 1998), based itself on the W94 precession rate correction and $\epsilon_0$value, $84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{410}$.

Table 4: Initial values at J2000.0 used for computing the precession of the equator.

This shows that MHB estimate for the precession rate in longitude is dependent on the value that is used for the factor $({\rm d}P_R/{\rm d}H_{{\rm d}R})$and on the theoretical values for the other contributions to r₀ as well.

The change in (r₀)₁ that results from the differences between the values used in this paper (see Tables 3 and 4) for providing our solutions (P03 $_{\rm prel}$ and P03) and the MHB ones is:

 

$${\rm d}(r_0)_1=3914~\mu{\rm as}/\rm cy .$$ (34)

   
6.2.4 Evaluation of the constants of integration

The constants of integration to be used for solving the precession equations must take proper account of all the perturbing effects that have been evaluated in Sect. 6.2.3.

First, the contribution (32) corresponding to the effect (i) has to be subtracted from the MHB estimate for $\psi_1$ before it is considered to be the physical precession rate along the ecliptic of epoch.

Second, the contribution (33) corresponding to the effect (ii) has to be subtracted from the observed values for $\psi_1$ and $\omega_1$ before using them as the observed quantities for r₀ and u₀ respectively.

In contrast, the contribution (34) corresponding to the effect (iii) does not modify the value for r₀ itself, but it modifies the value for (r₀)₁ (and consequently for $H_{\rm d}$) that can be derived from the MHB observed value for r₀ (see Table 4). This effect does not modify the secular term for the solution in $\psi _{\rm A}$, and only very slightly that for its quadratic term.

The integration constants for the precession Eqs. (24) and (26) are computed from the MHB estimates for r₀ and u₀ taking into account the contributions that are listed on Table 3. The r₀ and u₀ values to be used must be derived from MHB estimates, using the following relations (in $\mu$as/cy):

 

$$r_0 = (r_0)_{\rm MHB} + 2757$$

$$u_0 = (u_0)_{\rm MHB} - 514 .$$ (35)

The contributions (32) and (33) to the MHB estimates and that from nutation (see Sect. 4.4) have not been considered in other existing solutions for precession. For comparison purposes with these solutions, we will therefore keep the preliminary solution not corrected for these spurious contributions, denoting it P03 $_{\rm prel}$. The solution corresponding to relation (35), that is considered as being the final solution, is denoted P03.

Table 4 provides a comparison between numerical values for the fundamental parameters corresponding to different solutions and Table 5 provides the corresponding contributions to precession rates that have been taken into account. The first column in this Table provide the first order terms in precession rates in longitude and obliquity respectively, that can be derived from the observed values, given the other theoretical contributions.

Table 5: Contributions in arcseconds to the precession rates used in different models.

6.2.5 Solution for precession quantities compliant with IAU 2000A

Following the approach described in Sect. 6.2.1, we used analytical and semi-analytical tools for computing corrected expressions for the precession quantities, given the motion of the ecliptic as obtained in Sect. 6.1.1 and the expressions for the quantities $r_{\psi }$ and $r_{\epsilon }$ provided by Tables 3, 4 and 5. The differential equations for $\psi _{\rm A}$, $\omega _{\rm A}$ and $\epsilon_{\rm A}$ and $p_{\rm A}$ were solved using the GREGOIRE software, starting from the IAU 1976 expressions, which provide numerical expressions for these quantities including corrections to the t² and t³ terms with respect to the previous ones. Then, using relations in the spherical triangle mentioned earlier, the corresponding numerical expression for  $\chi_{\rm A}$ was derived. Two iterations of this process achieved convergence at a sub-microarcsecond level, providing the final solutions.

The program was tested against IAU 1976 expressions and against W94 solutions. With the input of the ecliptic expressions and the numerical values for the precession rate contributions, the process was able to reproduce the polynomials in the L77 and W94 papers up to the last digits of these expressions.

Table 6 compares the developments, obtained from different solutions, for the two basic quantities for the precession of the equator. (Here, and elsewhere in the paper, each coefficient is quoted to a number of digits that delivers $1~\mu$as resolution for dates close to 2000, degrading to $10~\mu$as resolution for values of t up to $\pm$10, corresponding to $\pm$1 millennium. It should be understood that this is a numerical convention and does not necessarily imply that any given coefficient is known to the quoted accuracy.)

Table 6: Comparisons between solutions for the precession of the equator and expressions for the precession rates.

It should be noted that the improvement from the IAU 2000 to the P03 solutions is of the order of 5 mas/cy² in $\psi _{\rm A}$and of 25 $\mu$as/cy² in $\omega _{\rm A}$. It should also be noted that, except for the secular term due to the correction applied in removing spurious effects, the differences between the coefficients of the expressions of our P03 $_{\rm prel}$ and P03 solutions, are at a microarcsecond level (P03 being our final solution).

6.2.6 Analytical expressions for the coefficients

Analytical expressions for the precession quantities are necessary for understanding each contribution of the precession of the equator and the precession of the ecliptic, and the coupling effects between them. Hilton (2002) proposed an extension of the L77 analytical expressions to IAU 2000 as functions of the L77 coefficients $\psi_1$, $\psi'_1$, $\psi''_1$, $\omega_1$, $\omega'_1$, $\omega''_1$, etc. As the computations performed in this paper include various precession rate contributions that have different effects in the solutions according to their epsilon-dependence, we have instead chosen the precession rates in longitude and obliquity as provided by (29) to be the basic quantities for the precession of the equator. The basic expressions for the precession of the ecliptic are the polynomial developments of the coordinates of the ecliptic pole, which are expressed as:

 

$$P_{\rm A}=s_1t+s_2t^2+s_3t^3+ s_4 t^4 + s_5 t^5$$

$$Q_{\rm A}=c_1t+c_2t^2+c_3t^3+ c_4 t^4 + c_5 t^5.$$ (36)

We used the Maple software to derive analytically the expressions for the coefficients of the solutions of Eqs. (24) and (26). These expressions include the additional $\omega_1$ term (=u₀), which was absent in L77. Table 7 provides the resulting expressions, up to the third degree, for the coefficients of the precession quantities as functions of the coefficients of the expressions (29) and (36).

Table 7: Expressions for the coefficients of the precession quantities.

7 An improved precession model

The only independent quantities for precession are the two direction cosines of the pole of the equator and those of the pole of the ecliptic. In this paper, we have clearly separated precession of the equator and precession of the ecliptic and we have obtained the developments of the quantities as functions of time through two independent approaches.

7.1 Final expressions for the primary angles

The quantities $\psi _{\rm A}$, $\omega _{\rm A}$, for the celestial pole and $P_{\rm A}$, $Q_{\rm A}$ for the ecliptic will be considered as being the primary precession quantities from which the others can be derived. Their developments are given below. The coefficients are in arcseconds and the time unit is 1 century of TDB (or TT in practice).


Precession of the equator, P03 solution:

 

$$\psi_{\rm A} 5038\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00114045}~t^3$$ =

$$+\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000951}~t^5$$

$$\omega_{\rm A} \epsilon_0 - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hs... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00772503}~t^3$$ =

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000003337}~t^5$$ (37)

with $\epsilon_0 = 84381\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{406}$.

Precession of the ecliptic, P03 solution:

 

$$P_{\rm A} +\ 4\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00022466}~t^3$$ =

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000120}~t^5$$

$$Q_{\rm A} -\ 46\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1e... ...+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00052413}~t^3$$ =

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000172}~t^5.$$ (38)

7.2 Final expressions for the derived angles and GMST

The classical "general precession,'' which mixes the motion of the equator in the GCRS and the motion of the ecliptic in the ICRS (and moreover may not be defined in the framework of General Relativity without fundamental problems) should no longer be regarded as a primary precession quantity. It is considered here as a derived quantity, along with the other precession quantities that can be obtained from the primary ones.

The expression for GMST(UT1,TT) provided by Capitaine et al. (2003b) must be revised in order to take into account the improvement in the expressions for the precession quantities (mainly $\chi_{\rm A}$ and $\psi _{\rm A}$). In contrast, the expressions for the periodic part of the quantity s and for the complementary terms in the equation of the equinoxes (see Capitaine et al. 2003a,b) are unchanged.

The P03 developments for the quantities $\epsilon_{\rm A}$, $p_{\rm A}$ and $\chi_{\rm A}$ are given below. The coefficients are in arcseconds and the time unit is TDB century.

 

$$\epsilon_{\rm A} = \epsilon_0 - 46\hspace{0.1em}'\hspace{-0.06em}... ...+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00200340}~t^3$$

$$\qquad\quad-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\h... ...0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000434}~t^5$$

$$p_{\rm A} = 5028\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\... ...+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00007964}~t^3$$

$$\qquad\quad-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\h... ...0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000383}~t^5$$

$$\chi_{\rm A} = 10\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00121197}~t^3$$

$$\qquad\quad+\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\h... ...\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000560}~t^5.$$ (39)

P03 developments for the equatorial precession angles ${\zeta _{\rm A}}$, $z_{\rm A}$and ${\theta _{\rm A}}$ can be derived from the above developments for $\psi _{\rm A}$, $\omega _{\rm A}$, $\epsilon_{\rm A}$and $\chi_{\rm A}$:

 

$$\zeta_{\rm A} 2\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{6... ... + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{2988499}~t^2$$ =

$$+\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000003173}~t^5$$

$$z_{\rm A} -\ 2\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ... + 1\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0927348}~t^2$$ =

$$+\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000002904}~t^5$$

$$\theta_{\rm A} 2004\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...- 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{04182264}~t^3$$ =

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000001274}~t^5 .$$ (40)

P03 developments for the ecliptic precession angles $\pi_{\rm A}$ and $\Pi_{\rm A}$ can be derived from the developments for the basic ecliptic quantities $P_{\rm A}$ and $Q_{\rm A}$.

 

$$\pi_{\rm A}=\ +~ 46\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em... ... 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{00012559}\ t^3$$

$$\qquad ~ ~~+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hs... ...\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000022}\ t^5$$

$$\Pi_{\rm A} = 629546\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3e... ... + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{157992}\ t^2$$

$$~ ~~\qquad-~ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hs... ...\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{000000072}\ t^5.$$ (41)

The Greenwich Mean Sidereal Time, expressed in terms of the Earth Rotation Angle $\theta$, becomes:

 

$${\rm GMST_{P03}(UT1,}~t) \theta({\rm UT1}) + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{014506}$$ =

$$+\ 4612\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.... ... + 1\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{3915817}~t^2$$

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ... 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{000029956}~t^4$$

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000000368}~t^5.$$ (42)

Note that the secular term in this expression includes a contribution of $- 3872~\mu$as from nutation (for more details, see Capitaine et al. 2003b).

This can be converted from angle units to time units by writing out the ERA and expressing the coefficients in seconds of Sidereal Time. To a resolution of 0.1 microsecond this gives:

 

$${\rm GMST}_{P03}(t_{\rm u},t) UT1 +\ 24110.5493771$$ =

$$+\ 8640184.79447825~t_u$$

$$+\ 307.4771013~(t-t_u)$$

$$+\ 0.092772110~t^2 -\ 0.0000002926~t^3$$

$$-\ 0.00000199708~t^4$$

$$-\ 0.000000002454~t^5$$ (43)

seconds, where $t_{\rm u}$ is the UT1 and t is the TT, both expressed in Julian centuries after J2000.

7.3 Expressions for combined frame bias and precession

   
7.3.1 The rotation vector approach

Although the most general, and potentially most compact, way of formulating a model for predicting the orientation of the Earth's axis is as a combined bias-precession-nutation effect, most existing applications treat precession and nutation separately as well as having no concept of frame bias. In such cases it would be convenient to regard bias as a component of precession and to express the bias-precession combination in a single model. Such a model could then replace the existing precession model, enabling GCRS coordinates to be transformed directly into mean place, leaving the nutation (luni-solar plus planetary) to be applied as usual in order to obtain true place.

Therefore, our goal is a compact formulation for the matrix:

$$\vec{R}_{PB} = \vec{P} ~ \vec{B} ,$$ (44)

where $\vec{R}_{PB}$ combines the individual rotation matrices for frame bias ($\vec{B}$) followed by precession ($\vec{P}$). The bias matrix $\vec{B}$ is given by the first of expressions (4). The precession matrix $\vec{P}$ is given by the second of expressions (4), using the improved $\epsilon_0$, $\psi _{\rm A}$, $\omega _{\rm A}$ and $\chi_{\rm A}$ from expressions (37) and (39).

The obvious formulation for such a combined bias-precession model is the traditional 3-rotation approach (5) used in the IAU 1976 precession model. At first sight, this offers the prospect of being able to replace the IAU 1976 precession matrix with one for IAU 2000 combined bias-precession simply by introducing revised polynomial expressions for ${\zeta _{\rm A}}$, $z_{\rm A}$and ${\theta _{\rm A}}$. However, because of the non-zero pole displacement at epoch, the $z_{\rm A}$ and ${\zeta _{\rm A}}$ angles representing such a combined bias-precession rotation undergo rapid changes around J2000 that make polynomials in t impractical. Thus a simple revision of IAU 1976 is unfortunately not feasible if the frame bias is to be included. This effect has also been noted by Fukushima (2003).

An alternative approach is to model the Cartesian components of what we may call the rotation vector. It is well known that any finite rotation of the coordinate frame can be expressed as the "Euler axis and angle'', which are, respectively, the unit vector along the axis of rotation and the amount of rotation. These can be combined in various ways, including the four "Euler symmetrical parameters'', often expressed in quaternion form, and the three components of the "Gibbs vector''; see Wertz (1986). A particularly straightforward three-component option is simply to scale the Euler-axis unit vector by the amount of rotation in radians. This "rotation vector'' approach proves more efficient for representing the bias-precession than either the quaternion or the Gibbs vector, because the precession approximates a constant rotation about a fixed point (namely the ecliptic pole).

To obtain expressions for the components of the rotation vector, we generated $\vec{R}_{PB}$ for a series of dates between 1800 and 2200, used standard transformations to express the matrix as the Euler axis and angle, and fitted polynomials in t to the components of the resulting rotation vector:

 

$$x_{\rm r} = + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\h... ... + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000236}~t^2$$

$$\qquad - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000004}~t^4$$

$$y_{\rm r} = - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\h... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4294924}~t^2$$

$$\qquad - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace... ... + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000092}~t^4$$

$$z_{\rm r} = - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\h... ... - 1\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{3915844}~t^2$$

$$\qquad + 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace... ...+ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000300}~t^4.$$ (45)

To express the vector $(x_{\rm r},y_{\rm r},z_{\rm r})$ as the rotation matrix $\vec{R}_{PB}$, we first decompose it into the amount of rotation in radians:

$$\phi = \left(x_{\rm r}^2 + y_{\rm r}^2 + z_{\rm r}^2\right)^{1/2}$$ (46)

and the rotation-axis unit vector:

$$x = x_{\rm r} / \phi,~~ y = y_{\rm r} / \phi,~~ z = z_{\rm r} / \phi.$$ (47)

Then, writing $s = \sin{\phi}$,  $c = \cos{\phi}$  and f = 1 - c,  we form the matrix elements as follows:

$$\vec{R}_{PB} \simeq \left( \begin{array}{lll} ~xxf~+~c~~~ & ... ...\ ~zxf~+~sy~~~ & zyf~-~sx~~~ & zzf~+~c~~ \end{array} \right).$$ (48)

Note that the algorithm is computationally very efficient: only one square root and two trigonometrical functions (of the same angle) are required, in addition to arithmetic operations.

In expressions (45), the choice of polynomial order and coefficient resolution was made on the basis that the formulation should reproduce the rigorous method to $1~\mu$as accuracy over a 400-year time span. A somewhat shorter version, consisting of 11 coefficients of $1~\mu$as resolution, delivered $50~\mu$as accuracy, and other compromises are of course possible.

To summarize, given the date t in Julian centuries after J2000, expressions (45) to (48) can be evaluated to generate the matrix $\vec{R}_{PB}$, and the product of this matrix with the GCRS vector gives the mean place of date.

7.3.2 The mean pole X and Y approach

The position of the CIP in the GCRS, as defined by IAU Resolution B1.7, is a function of frame bias, precession and nutation. The CIP coordinates can be obtained by evaluating the corresponding sequence of nine or ten coordinate rotations, or to use the simple expressions for X and Y derived from the improved expressions (37) and (39) for the classical precession quantities $\psi _{\rm A}$, $\omega _{\rm A}$ and $\chi_{\rm A}$ and the MHB nutations in longitude and obliquity.

Taking into account the polynomial part only, the P03 expressions for X and Y are:

 

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ... - 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4297829}~t^2$$                           X_P03 =

$$-~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ... 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{000007578}~t^4$$

$$+\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000059285}~t^5$$ (49)

 

$$-\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em... ...- 22\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{4072747}~t^2$$                           Y_P03 =

$$+~0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}... ... 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{001112526}~t^4$$

$$+\ 0\hspace{0.1em}'\hspace{-0.06em}'\hspace{-0.3em}.\hspace{0.1em}{0000001358}~t^5.$$ (50)

The purpose of this paper is to provide positions of the CIP computed to an accuracy of a few microarcseconds over a time span of a few hundred years, meeting the requirements of high-accuracy applications. For such a time span, the X, Y formulation achieves a similar accuracy to the classical one based on the basic precession quantities. The GCRS position of the CIP provided by expressions (49) and (50) is thus in agreement with that provided by expressions (37) and (38) to microarcsecond accuracy for a few centuries. The use of the X,Y coordinates for long-term studies would require a different development in trigonometric functions of the precession angles, which will not be discussed here.

It should be noted that changing from the P00 precession solution to the P03 solution gives rise to changes of the order of 0.1 $\mu$as to a few $\mu$as for a century in a few terms of the periodic part of the expressions for the CIP X, Y and in the polynomial part of the quantity s+XY/2 that provides the GCRS position of the CEO. The P03 expression for this latter quantity is, in $\mu$as:

 

      (s + XY/2)_P03 = 94.0 + 3808.65 t - 122.68 t²

$$-\ 72574.11 t^3 + 27.98 t^4 + 15.62 t^5.$$ (51)

The only change larger than $0.5~\mu$as with respect to the IAU 2000 expression for the GCRS position of the CEO is of $2.7~\mu$as in the quadratic term, the other changes all being less than $0.5~\mu$as.

8 Summary and concluding remarks

In the work described in this paper we have computed new expressions for precession consistent with the IAU 2000A precession-nutation model. The precession of the ecliptic has been derived from the analytical theory VSOP87 fitted to the JPL ephemerides DE406 for improving the polynomial terms in the expression for the component of the EMB orbital angular momentum with respect to a fixed ecliptic. It uses the value for the mean obliquity of the ecliptic at J2000 as derived from a fit of the dynamical theory for the Moon to LLR observations. The equinox offset in the GCRS has been derived from this fit based on VLBI Earth Orientation Parameters. The model for the precession of the equator has been obtained by solving the dynamical precession equations based both on the most recent expressions for the theoretical contributions to precession (W94) and on the MHB estimates of the precession rates. In this computation, we took proper account of all the perturbing effects on the observed quantities.

We have moreover discussed the most suitable precession quantities to be considered in order to be based on the minimum number of variables and to be the best adapted to the most recent models and observations.

The "rotation vector'' method, that is able to express bias plus precession, has some intuitive appeal, as it specifies the point near the ecliptic pole about which the rotation occurs, scaled by the accumulated precession. The rotation vector expressions (45) are ideal for rapid generation of the precession matrix in practical software for such purposes as pointing telescopes, and can be extended in order and resolution as required to meet more demanding needs.

We also provide the solutions for the X, Y coordinates of the CIP in the GCRS which include precession, nutation and frame biases and therefore tie the precession nutation directly to the ICRS by simply providing where the pole is in the sky. Their advantage is of being close to the parameters that are actually observed by VLBI, which is the best way of determining the precession-nutation motion and is not sensitive to an ecliptic. Note for example that the precession in longitude that is derived from VLBI is in fact the projection of the X-coordinate of the CIP on a conventional ecliptic and changing the obliquity of this conventional ecliptic will change the value in longitude, whereas the corresponding value for X is independent of the ecliptic (See Sect. 6.2.3).

Recent papers by other authors, using different methods, furnish useful comparisons with our results. For example, Bretagnon et al. (2003) have provided precession expressions derived from the theory SMART97 of the rotation of a rigid Earth using the MHB observed precession rate in $\psi _{\rm A}$. The precession of the ecliptic is provided by VSOP87 plus improved values for planetary masses (IERS 1992). Fukushima (2003) has used a fit to numerical ephemerides for improving the precession of the ecliptic and a fit to VLBI data for improving the precession of the equator, given his solution for nutation. Table 8 provides a comparison of the different expressions (denoted B03 and F03 respectively) and Figs. 4 to 7 compare the final solution of this paper (denoted P03) for the precession of the ecliptic and the precession of the equator.

Table 8: Comparisons between expressions for the basic precession quantities.

	Figure 4: Comparisons of models for the precession of the ecliptic (quantity $P_{\rm A}$): differences with respect to the solution P03 provided in this paper.
Open with DEXTER

	Figure 5: Comparisons of models for the precession of the ecliptic (quantity $Q_{\rm A}$): differences with respect to the solution P03 provided in this paper.
Open with DEXTER

	Figure 6: Comparisons of models for the precession of the equator (quantity $\psi _{\rm A}$): differences with respect to the solution P03 provided in this paper.
Open with DEXTER

	Figure 7: Comparisons of models for the precession of the equator (quantity  $\omega _{\rm A}$): differences with respect to the solution P03 provided in this paper.
Open with DEXTER

Regarding the precession of the ecliptic, we note the good agreement with the B03 solution (in the graph the B03 and W94 plots are practically indistinguishable), differing only, at the level of precision provided by the figures, by a secular trend, whereas there are very large discrepancies with respect to L77 (IAU 1976) and F03.

Regarding the precession of the equator, we note the quadratic difference in $\psi _{\rm A}$ with respect to the IAU 2000 and B03 solutions which can be explained by the fact that both solutions (except for their secular term) are relative to a rigid Earth. The large discrepancy in $\omega _{\rm A}$ with respect to F03 cannot be explained simply by the difference in the expressions for the ecliptic precession on which the solutions for $\omega _{\rm A}$ are based.

The work described here has taken advantage of the best available observations (VLBI for the equator, LLR for the ecliptic), of the most recent theories for the Earth (VSOP87) and the Moon (ELP2000), and the most precise numerical ephemerides (DE406), to develop expressions of the precession of the equator and the ecliptic that are compliant with the IAU 2000 resolutions and that are dynamically consistent.

Future VLBI observations covering a longer period of time will allow improved separation between the estimates of precession rates and amplitudes of long period nutation. They will also allow the t² terms in the developments to be estimated, providing an indirect access to the motion of the ecliptic.

Acknowledgements

We are grateful to the referee of the paper, J. Vondrák, for his valuable suggestions to improve the text. We also thank J. Hilton and G. Kaplan for their useful tests of the solution provided in this paper and P.M. Mathews for valuable discussion.

References

 
• Bretagnon, P., & Chapront, J. 1981, A&A, 103, 103 In the text NASA ADS  
• Bretagnon, P., & Francou, G. 1988, A&A, 202, 309 In the text NASA ADS  
• Bretagnon, P., Rocher, P., & Simon, J.-L. 1997, A&A, 319, 305 In the text NASA ADS  
• Bretagnon, P., Fienga, A., & Simon, J.-L. 2003, A&A, 400, 785 In the text NASA ADS  
• Brumberg, V. 1991, Essential Relativistic Celestial Mechanics (Adam Hilger Publ.) In the text  
• Brumberg, V. 2003, Private communication In the text  
• Capitaine, N., Williams, J. G., & Seidelmann, P. K. 1985, A&A, 146, 381 In the text NASA ADS  
• Capitaine, N., Chapront, J., Lambert, S., & Wallace, P. T. 2003a, A&A, 400, 1145 In the text NASA ADS  
• Capitaine, N., Wallace, P. T., & McCarthy, D. D. 2003b, A&A, 406, 1135 NASA ADS  
• Chapront, J., Chapront-Touzé, M., & Francou, G. 2002, A&A, 387, 700 In the text NASA ADS  
• Dehant, V., & Capitaine, N. 1997, Celest. Mech. Dyn. Astr., 65, 439 In the text NASA ADS  
• de Sitter, & Brouwer 1938, Bull. Astron. Inst. Neth., 8, 213 In the text NASA ADS  
• Folkner, W. M., Charlot, P., Finger, et al. 1994, A&A, 287, 279 In the text NASA ADS  
• Fricke, W. 1971, A&A, 13, 298 In the text NASA ADS  
• Fukushima, T. 1991, A&A, 244, L11 In the text NASA ADS  
• Fukushima, T. 2003, AJ, 126, 1 In the text NASA ADS  
• Guinot, B. 1979, in Time and the Earth's Rotation, ed. D. D. McCarthy, & J. D. Pilkington (D. Reidel Publishing Company), 7 In the text  
• Hartmann, T., Soffel, M., & Ron, C. 1999, A&AS, 134, 271 In the text NASA ADS  
• Herring, T., Mathews, P. M., & Buffett, B. A. 2002, J. Geophys. Res., 107, B4, 10.1029/ 2001JB000390 In the text  
• Hilton, J. 2002, IERS Technical Note 29, ed. N. Capitaine, D. Gambis, D. D. McCarthy, et al. (Frankfurt am Main: Verlag des Bundesamts für Kartographie und Geodäsie), 101 In the text  
• IERS Conventions 2000, http://www.usno.mil/Conventions2000, draft  
• Kinoshita, H. 1977, Celest. Mech., 15, 277 In the text NASA ADS  
• Laskar, J. 1986, A&A, 157, 59 In the text NASA ADS  
• Lieske, J. H. 1985, Celest. Mech., 37, 209 In the text NASA ADS  
• Lieske, J. H., Lederle, T., Fricke, W., & Morando, B. 1977, A&A, 58, 1 In the text NASA ADS  
• Mathews, P. M. 2002, private communication In the text  
• Mathews, P. M., Herring, T. A., & Buffett, B. A. 2002, J. Geophys. Res., 107, B4, 10.1029/2001JB000165  
• McCarthy, D. D. (ed.) 1992 , IERS Standards, IERS Technical Note 13 (Paris: Observatoire de Paris) In the text  
• Newcomb, S. 1894, Astron. Papers Am. Ephemeris 5, Part 4, 301  
• Roosbeck, F., & Dehant, V. 1998, Celest. Mech., 70, 215  
• Seidelmann, P. K. 1982, Celest. Mech., 27, 79 In the text NASA ADS  
• Simon, J. L., Bretagnon, P., Chapront, J., et al. 1994, A&A, 282, 663 In the text NASA ADS  
• Simon, J.-L., Chapront-Touzé, M., Morando, B., & Thuillot, W. 1997, Introduction aux Ephémérides astronomiques, Supplément à la Connaissance des Temps (Éditions de Physique), 127 In the text  
• Souchay, J., Loysel, B., Kinoshita, H., & Folgueira, M. 1999, A&AS, 135, 111 In the text NASA ADS  
• Souchay, J., & Kinoshita, H. 1996, A&A, 312, 1017 In the text NASA ADS  
• Standish, E. M. 1981, A&A, 101, L17 In the text NASA ADS  
• Standish, M. 2000, JPL Planetary and Lunar Ephemerides DE406  
• Wahr, J. H., 1981, Geophys. J. R. Astr. Soc., 64, 705  
• Wertz, J. M. 1986, Spacecraft Attitude Determination and Control, ed. R. James, & Wertz, Astrophysics and Space Science Library, vol. 73 (D. Reidel Publishing Company) In the text  
• Williams, J. G. 1994, AJ, 108, 711 In the text NASA ADS  
• Williams, J. G. 1995, AJ, 110, 1420 In the text NASA ADS  
• Woolard, E. W. 1953, Astr. Pap. Amer. Ephem. Naut. Almanach XV, I, 1 In the text  
• Woolard, E. W., & Clemence, G. M. 1966, Spherical Astronomy (New-York, London: Academic Press) In the text