A&A 412, 567-586 (2003)
DOI: 10.1051/0004-6361:20031539
N. Capitaine^{1} - P. T. Wallace^{2} - J. Chapront ^{1}
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
,
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
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.
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 for the GCRS-to-ITRS rotation matrix, omitting polar motion, with elements .
Using the usual notation, the CEO-based transformation is written as:
For practical reasons, X and Y are usually called "coordinates'' and their numerical expressions are multiplied by the factor 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
(3) |
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
,
and :
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 and , 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 and 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 (mean obliquity of date) for the inclination of the mean equator of date on the ecliptic of date and (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, , for the speed of precession, and for the geodesic precession (de Sitter & Brouwer 1938). The expression for general precession (denoted 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).
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), /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:
The IAU 1976 precession model has the following limitations:
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.
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 (
years) the
contribution of which can be approximated, in as, as:
(9) |
The IAU 2000 nutation series is associated with improved numerical
values for the precession rate of the equator in longitude and
obliquity:
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:
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.
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 and . The canonical formulation for the IAU 2000 precession matrix is thus naturally the 4-rotation one (see relations (4)), involving , , and .
The IAU 2000 expressions for the classical quantities
and
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:
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
(obliquity of the equator on the
moving ecliptic) and
(general precession in longitude),
whereas the expression for
is unchanged. By
adding the corrections
to
and dto
,
respectively, one obtains:
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,
and ,
which are derived
from (11) as:
Taking into account precession only, these expressions, denoted P00, become:
The contributions d
and d
from the frame biases
in X and Y are (in as):
(20) |
The computation of improved expressions for the currently used
precession quantities ,
and
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
and ,
at one level of accuracy
less than the MHB corrections themselves (i.e. 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
and
;
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^{5} and the
precision of the coefficients as. The following series
with a as level of precision matches the canonical
4-rotation series to sub-microarcsecond accuracy over 4 centuries:
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 , 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:
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 is no longer appropriate as the primary expression for precession, and the expression for is useful only as an intermediate in the coordinate transformation between TRS and CRS, as its effect is canceled out in the GST expression.
In contrast to the previous approach, a physical meaning can be given to the corrections to precession rates. The MHB correction d 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 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 (dand d) 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^{2} and t^{3} terms of the precession quantities.
There is however a problem of numerical consistency. The t^{2}and t^{3} 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.
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.
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 and , 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 and . 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 and 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 and are, in the case of DE406: mas/cy and 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.
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
:
the
(inertial) ecliptic pole lies in the direction of
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 and 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 , ) 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:
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 and , the terms in t^{5} 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
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
,
and changed the obliquity by 3.329 mas from the inertial
value, to
.
Figure 2: Differences P(DE406) - P(VSOP) in . | |
Open with DEXTER |
Figure 3: Differences Q(DE406) - Q(VSOP) in . | |
Open with DEXTER |
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 and . Over the time interval [1900, 2100], we computed the linear regressions of the residuals: and . The tests consisted of varying the set of abscissas, the length of the time interval modified by 10 percent and a shift of the center by 10 years. Our conclusion from these tests is that our determination of the secular variations of and 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 , where A is the amplitude, the associated power of time, is the phase and the frequency, the time being reckoned from J2000.0. The period 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 , the constant term in q and its trend 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: . The planetary and lunar arguments are , , , , D = Delaunay's argument D (the difference of the mean longitudes of Moon and Sun), 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 was defined by the MHB 2000 precession-nutation model itself, but it was necessary to position the equinox by choosing a value for , thereby completing the 3D frame rotation.
The adopted value 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 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 would be justified. This would give a 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 change of 2.5 mas would be only one term of as/cy in Y, which is not of practical significance.
2. The value cannot be evaluated to better than about 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 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.
The precession of the equator with respect to a fixed frame can be derived by solving the differential equations for the precession quantities , and , 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 and 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 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 , 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:
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:
Using solutions (37) for and , 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 as for a century are (i) a bias in obliquity of as, that is in fact included in the adopted value for the mean obliquity at epoch, and (ii) secular terms of as/cy and as/cy in obliquity and longitude, respectively. It can be shown that the solution of Eqs. (24) is insensitive, within 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 and , which depend both on the motion of the equator and the motion of the ecliptic, as does the quantity , which is also required for computing .
The differential equations for
and
can be derived from the
expression for
and
respectively and the ecliptic
precession (see Laskar 1986; Simon et al. 1997; or W94 for the equation
in
), which gives:
A(p,q) | = | ||
B(p,q) | = | ||
C(p,q) | = | (27) |
The quantity
can be derived from geometrical
considerations in the spherical triangle N
,
where N is the node of the moving ecliptic in the ecliptic at
J2000 and
and
are the mean equinoxes of date
and of J2000 respectively. Using the definition of the quantities
and
as functions of elements
and
of
the spherical triangle, this gives:
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
and
,
that can be developed as polynomials in
t:
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_{4} effects in the direct luni-solar torque, (ii) direct planetary torque effect, (iii) J_{2} and planetary tilt effects, (iv) tides and 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, (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 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:
(30) |
The effect of geodesic precession provided in this table has been obtained by converting expression (31) for the effect in the quantities and into contributions (usually denoted ) and to the precession rates and 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_{0})_{1} and (u_{0})_{1} respectively. The (r_{0})_{1} term is the only one with a sufficiently large amplitude (of the order of ) to be sensitive to small changes in the value for the dynamical ellipticity of the Earth, which is one of the Basic Earth Parameters to be determined. The (r_{0})_{1} term is therefore of very special significance for computing the integration constant r_{0} to be used for solving the precession equations. Its value is derived from the adopted value for r_{0}, given the other contributions which are provided by the theory.
It should be noted that the quantity (r_{0})_{1} is related to the so-called "Newcomb's precessional constant'' at epoch, , by , whereas the quantity r_{0} 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 provided in Table 3 correspond to effects considered in the expression given by Kinoshita (1977) for the speed of luni-solar precession, f_{2000}, 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, and , of the equator used in the present paper.
This, in particular, introduces the following spurious contributions
to the estimated precession rates:
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 and on the theoretical values for the other contributions to r_{0} as well.
The change in (r_{0})_{1} that results from the differences between the values
used in this paper (see Tables 3 and 4) for providing our solutions
(P03
and P03) and the MHB ones is:
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 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 and before using them as the observed quantities for r_{0} and u_{0} respectively.
In contrast, the contribution (34) corresponding to the effect (iii) does not modify the value for r_{0} itself, but it modifies the value for (r_{0})_{1} (and consequently for ) that can be derived from the MHB observed value for r_{0} (see Table 4). This effect does not modify the secular term for the solution in , 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_{0}
and u_{0} taking into account the contributions that are listed on Table 3.
The r_{0} and u_{0} values to be used must be derived from MHB estimates, using
the following relations (in as/cy):
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.
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 and provided by Tables 3, 4 and 5. The differential equations for , and and were solved using the GREGOIRE software, starting from the IAU 1976 expressions, which provide numerical expressions for these quantities including corrections to the t^{2} and t^{3} terms with respect to the previous ones. Then, using relations in the spherical triangle mentioned earlier, the corresponding numerical expression for 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 as resolution for dates close to 2000, degrading to as resolution for values of t up to 10, corresponding to 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^{2} in and of 25 as/cy^{2} in . 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 and P03 solutions, are at a microarcsecond level (P03 being our final solution).
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 ,
,
,
,
,
,
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:
Table 7: Expressions for the coefficients of the precession quantities.
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.
The quantities
,
,
for the celestial pole and
,
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:
= | |||
(37) |
= | |||
(38) |
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 and ). 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
,
and
are given below.
The coefficients are in arcseconds and the time unit is TDB century.
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:
seconds, where is the UT1 and t is the TT, both expressed in Julian centuries after J2000.
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:
(44) |
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 , and . However, because of the non-zero pole displacement at epoch, the and 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
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:
(46) |
(47) |
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 as accuracy over a 400-year time span. A somewhat shorter version, consisting of 11 coefficients of as resolution, delivered 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 , and the product of this matrix with the GCRS vector gives the mean place of date.
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 , and and the MHB nutations in longitude and obliquity.
Taking into account the polynomial part only, the P03 expressions for X and Y are:
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 as to a
few 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 as:
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 . 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 ): 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 ): 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 ): 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 ): 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 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 with respect to F03 cannot be explained simply by the difference in the expressions for the ecliptic precession on which the solutions for 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^{2} 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.