2 Theory

 

2.1 Some words about the principal methods of computation of the nutation

 Let us recall some aspects of the general strategy of how to compute the orientation of the Earth in space. Basically nutation, together with precession, represents the motion of the Earth's axis with respect to an inertial frame. To be more precise, the Earth as a member of the solar system undergoes a translational motion usually described by the motion of its center of mass in space and a rotational motion around the center of mass (spin, nutation, precession). The precession and nutation have the same origin, namely the gravitational attraction between the Earth's equatorial bulge and the other bodies. For the rotational motion that force acts on the spinning Earth and makes it moving like a gyro under the influence of a torque. To describe its rotational motion theoretically we can find three possibilities.

1) First, employing the "classic'' computation method by using Lagrange theory and a force function. Woolard's classical treatment (1953) was formulated in frames of Lagrange theory. Taking the Euler angles $\theta, \psi$ and $\varphi$ as generalized coordinates q_i, i=1,2,3, the Lagrangian reads

$$\begin{eqnarray} {\cal L} = U + \frac{1}{2} (A \omega_1^2 + B \omega_2^2 + C \omega_3^2) ,\end{eqnarray}$$ (1)

where A, B and C are the principle moments of inertia and U is the so-called force function for the external forces. Relating the components $\omega_i$ of the rotation vector of the principal axes system of the Earth with the Euler angles using Euler's kinematical equations (shortly EKE), we find:

$$\begin{eqnarray} \omega_1 &=& {\dot{\psi}} \sin {\theta} \sin {\varphi} +{\dot... ...ega_3 &=& {\dot{\varphi}} + {\dot{\psi}} \cos {\theta} . \nonumber\end{eqnarray}$$ (2)

Next, we have Lagrange's equations for a conservative system:

$$\begin{eqnarray} \frac{\mathrm d}{{\mathrm d}t} \frac{\partial{\cal L}}{\partial \dot{q}_i} - \frac{\partial {\cal L}}{\partial q_i} = 0.\end{eqnarray}$$ (3)

Then the Euler's dynamical equations (shortly EDE) are:

$$\begin{eqnarray} A \dot{\omega}_1 + (C-B) \omega_2 \, \omega_3 &=& \frac{\sin \... ...A) \omega_1 \, \omega_2 &=& \frac{\partial U}{\partial \varphi} .\end{eqnarray}$$ (4)

Integrating these differential equations yields the expressions that are required for the theory of the motion of the Earth relative to its center of mass. This is - together with the calculation of the force function U for the gravitational action between the Earth and an external massive gravitating body - the main task for this computation method. Analytical ephemerides are (almost) exclusively used for the determination of the positions of the external body from which U can be developed. Depending on the coordinate systems involved (especially those used for the ephemerides), the computation of U can be - despite of all symbol manipulating programs - quite laborious. Also solving the differential Eqs. (4) analytically is not an easy task.

2) Another method is based on the Hamiltonian formalism. This approach was introduced by Andoyer (1923, 1926) and perfected in Kinoshita (1977), Kinoshita & Souchay (1990) and Souchay & Kinoshita (1996, 1997a). Kinoshita has developed the analytical theory describing the motion of the planes normal to the angular momentum vector, to the figure axis and to the rotation axis of the rigid Earth having the shape of a triaxial ellipsoid. The Hamiltonian formalism which was used is characterized by the use of Andoyer variables, a moving fundamental reference frame and Hori's perturbation method. The force function U is to be based on the analytical ephemerides, likewise in the case of using of the Lagrange theory.

Although the force function U can be considered as a small quantity thus allowing for a perturbative treatment, some special attention has to be given to the separation of the precession (long-periodic contribution) and the nutation (short-periodic contribution). Kinoshita (1972), K77 and Kinoshita & Souchay (1990) describe in detail how they use canonical transformations to tackle this problem.

3) The other approach was used to compute the nutation model H96NUT presented here. It differs from the classical method mainly in three aspects. First, instead of the force function U the external torque $\vec{T}$ is used directly and expressed entirely in terms of tidal quantities which can be taken from the tidal potential development HW95. Roosbeek & Dehant (1997) used also the external torque but they developed it from the ephemerides of the perturbing bodies. Second, the differential equations are solved directly (using also a perturbation method) in the time domain. Third, the coordinate system used for most of the calculations is the system of principal axes of the Earth. The transformation to the ecliptic of date is applied only in a final stage. All of this is described in the following sections. More detailed derivations of the formulas can be found in Hartmann (1996).

2.2 The tidal potential development HW95

The first point to be discussed is the computation of the tidal potential, which replaces the rather complicated calculation of the disturbing function U in K77 or KS90. Here, only a summary of the most important aspects is given. More details can be found in Hartmann & Wenzel (1995a,b), where a time-harmonic development of the Earth tide generating potential (tidal potential) called HW95 due to the Moon, the Sun and the planets is presented. It contains 12935 tidal waves and is the most accurate tidal potential series available today (see Wenzel 1996). The expansion is carried out to the degree 6 of the geopotential for the Moon, degree 3 for the Sun and degree 2 for the planets Mercury, Venus, Mars, Jupiter and Saturn. It is entirely based on the JPL standard numerical ephemerides DE200. The secular variation of the biggest amplitudes is computed for the period between 1850 and 2150. The indirect effects of the planets, mainly due to Venus and Jupiter, on the Earth's orbit and the Moon's orbit are taken into account implicitly by the use of the numerical ephemerides. This also applies to the so-called J₂-tilt effect, where the J₂ coefficient of the Earth causes the orbit of the Moon to be slightly tilted.

The tidal potential V(t) due to a specific body (labeled b) at the surface of the Earth is given by

$$\begin{eqnarray} V{(t)} &=& G M_{\rm b} \sum_{\ell=2}^{\ell_{\rm max}} \sum_{... ...ta_{\rm b}(t)) \cos[m (\lambda - \Lambda_{\rm b}(t))]. \nonumber\end{eqnarray}$$ (5)

Here, r, $\theta$, $\lambda$ denote the geocentric spherical coordinates of the station on (or near) the Earth's surface where r is the radial distance, $\theta$ is the co-latitude and $\lambda$ is the longitude and $r_{\rm b}$, $\theta_{\rm b}$, $\Lambda_{\rm b}$ are those of the center of mass of the specific planet, whose mass is denoted by $M_{\rm b}$. G is the gravitational constant. The fully normalized Legendre functions of degree $\ell$ and order m used for the computation of the catalogue of the tidal potential HW95 are denoted by $\overline{P}_{\ell m}$. Then, the catalogue HW95 has been computed by using an expansion of the form

$$\begin{eqnarray} V{(t)} &=& \sum_{\ell=2}^{\ell_{\rm max}} \sum_{m=0}^{\ell} \... ...alpha_i(T)) + S_i^{\ell m}(T) \sin(\alpha_i(T)) \right] \nonumber\end{eqnarray}$$ (6)

where a is the major semi axis of the Earth and T is the time reckoned from J2000 in Julian centuries. The time dependent tidal potential coefficients (shortly the tidal amplitudes) are given by

$$\begin{eqnarray} C_i^{\ell m}(T) &=& {\rm C0}_{i}^{\ell m} + T \cdot {\rm C1}_{... ...l m}(T) &=& {\rm S0}_{i}^{\ell m} + T \cdot {\rm S1}_{i}^{\ell m},\end{eqnarray}$$ (7) (8)

where the potential coefficients ${\rm C0}_{i}^{\ell m}$ and ${\rm S0}_{i}^{\ell m}$ have the dimension m² s^-2 and the linear drift coefficients ${\rm C1}_{i}^{\ell m}$ and ${\rm S1}_{i}^{\ell m}$ have the dimension m² s^-2 cy^-1. They are - together with the tidal arguments $\alpha_i(T)$ - the fundamental quantities in this paper from which everything else is deduced. The tidal arguments $\alpha_i(T)$ are computed from  

$$\alpha_i(T) = m \cdot \lambda + \sum_{j=1}^{11} k_{ij} \cdot \arg_j(T).$$ (9)

The integer coefficients k_ij are given in the HW95 catalogue, while the eleven astronomical arguments $\arg_j(T)$($j=1,\ldots,11:$ $\tau = t + h -s - \lambda$ is the mean local lunar time and t is mean solar time, s = mean lunar longitude, h = mean solar longitude, p = mean longitude of lunar perigee, $N^\prime$ = negative mean longitude of the lunar ascending node, $p_{\rm s}$ = mean longitude of solar perigee and the mean longitudes of Mercury, Venus, Mars, Jupiter and Saturn indicated by $l_{\rm Me}$, $l_{\rm V}$,$l_{\rm Ma}$, $l_{\rm J}$ and $l_{\rm S}$, respectively which are referred to the mean dynamical ecliptic and equinox of date) are computed from polynomials in time, that can be found in Simon et al. (1994). This set of arguments is also the basis for the nutation, but with $s, h, p, N^\prime$ and $p_{\rm s}$ replaced by $l_{\rm m}, l_{\rm s}, F, D$ and $\Omega$ using the well-known linear relations: $l_{\rm m}=s-p,\ l_{\rm s}=h-p_{\rm s},\ F=s+N^\prime,\ D=s-h$ and $\Omega=-N^\prime$. For the nutation the tidal arguments may be considered as linear functions of time. Corrections to this approximation are discussed in the Sect. 3.5.

The numerical standard ephemerides DE200 (Standish & Williams 1981) supplied by Jet Propulsion Laboratory, Pasadena, have been used to compute the celestial geocentric vector $\vec{r}_{\mathrm{cel.}}$for the time interval 1850 to 2150. The conversion to terrestrial geocentric vector $\vec{r}_{\mathrm{terr.}}$uses the well-known relation  

$$\vec{r}_{\rm terr.} = {\cal S} \cdot {\cal N} \cdot {\cal P} \cdot \it \vec{r}_{\rm cel.}.$$ (10)

The mean obliquity, the precession constant and the precession formulas of Simon et al. (1994) have been used for ${\cal P}$including corrections for the new values of planetary masses. The IAU 1980 Nutation theory (Seidelmann 1982) has been used for ${\cal N}$ with the Delaunay arguments replaced by those of Simon et al. (1994). To assure the precision of the tidal potential of 10^-8 m² s^-2 we need the coordinates of the attracting bodies with an accuracy better than 30 mas. Taking into account the known difference of the IAU 1980 Nutation theory and recent observations which are less then 10 mas we are justified to use the IAU 1980 Nutation theory. The model of Earth used in the computation for ${\cal S}$ rotates according to the expression of Aoki et al. (1982, Eq. 14). The equation of equinoxes was added to obtain the Greenwich Apparent Sidereal Time (GAST)  

$${\rm GAST} = {\rm GMST} + \Delta\Psi \cdot \cos \varepsilon .$$ (11)

All numerical values for the astronomical constants were taken from the IERS 1992 Standards (McCarthy 1992).

It is important to note that the transformation from (5) to (6) requires the computation of a time series of the tidal potential (about 300 years were used) and a spectral analysis to find the individual terms in the expansion (6) see Hartmann (1996) and Hartmann & Wenzel (1995a,b) for details. In addition, the integer numbers k_ij for the tidal and nutational argument had to be found from the particular value of the interpolated Fourier-frequency. The numerical separation of waves with very small frequency differences (e.g. $2 p_{\rm s}$ corresponding to a period of 10 468 years which also occurs as a difference of two nutation arguments) as well as the final total fit of the tidal amplitudes caused some problems that will appear later also in the corresponding nutation series (see also Roosbeek 1996).

2.3 The relationship between the torque and the tidal potential

Next, the relationship between the external torque $\vec{T}$and the HW95 tidal potential quantities will be calculated. In Hartmann et al. (1994) it is shown that this formula can be derived without integration via symmetric, trace-free tidal moments. Using again the system based on the principal axes of the Earth we find for the main contribution of the tidal potential of degree $\ell = 2$(omitting the sum i over all tides whenever possible) the exact expression  

$$\vec{T} = \frac{\sqrt{15}}{a^2} \pmatrix{ (C-B) \; [ S_i^{... ...i^{22} \cos(\alpha_i(T)) + C_i^{22} \sin(\alpha_i(T)) ] \cr},$$ (12)

where a is taken to be the semi-major axis of the Earth from IERS 1992 Standards (McCarthy 1992). For an axi-symmetric Earth (A = B) only the tesseral part of the tidal potential (order m=1) is involved, but this is not assumed in the following in order to obtain the triaxiality terms. The sectorial part (m=2) slightly influences the length of day as shown in e.g. Wünsch (1991) but this can be ignored.

Similar relations can be found for the higher parts ($\ell=3,4,~\dots$) of the tidal potential (5). For the nutation coming from these parts the higher zonal mass multipole moments (J₃, J₄) are the dominant contributions to the torque. Then, the torque can again be expressed (now entirely) in terms of the tesseral part (m=1) of the tidal potential. Due to the advantageous representation (6) of the HW95 tidal potential, these contributions to the nutation can be obtained by replacing the scale factor and the tidal amplitudes (see Hartmann et al. 1996; for more details). To simplify the following calculations, the factor $\sqrt{15}/a^2$is absorbed into the tidal amplitudes.

2.4 The free rotational motion of the Earth

The rotational equation of motion in the figure axes system, in which the torque is given by (12), reads  

$$\frac{\mathrm d}{{\mathrm d}t} \vec{S} + \vec{\omega} \times \vec{S} = \vec{T} ,$$ (13)

where the torque $\vec{T}$ is included for later purposes; it vanishes for free motion. The spin (or angular momentum) vector $\vec{S}$ of the Earth is related to the rotation vector $\vec{\omega}$by the moment of inertia tensor $\tens{C}$, which - again in the figure axis system - reads  

$$\vec{S} = \tens{C} \cdot \vec{\omega} = \pmatrix{ A & 0 & 0 \cr 0 & B & 0 \cr 0 & 0 & C \cr } \cdot \vec{\omega} .$$ (14)

Written out in components this leads again to EDE given in (4) and which therefore relate the components of the rotation vector $\vec{\omega}$in the figure axes system to those of the torque, or with the help of the result of the last section, to the tidal quantities. The exact solution of EDE with vanishing torque can be deduced by means of elliptical functions. However, their modulus m

$$m= \frac{\omega_{10}^2}{\omega_{30}^2}\frac{A(B-A)}{C(C-B)}$$ (15)

is so small for the Earth ($\approx\!10^{-14}$)that they can be replaced by their trigonometrical counterparts. Therefore we proceed with

$$\begin{eqnarray} \omega_1(t) &\approx& \omega_{10} \cos( \mu \, (t-t_0) ) \nonu... ...\mu \, (t-t_0) ) \\ \omega_3(t) &\approx& \omega_{30} \nonumber .\end{eqnarray}$$ (16)

The integration constants are $\omega_{10}, \omega_{30}$ and t₀. The two other constants are defined by  

$$\omega_{20}^2 \equiv \omega_{10}^2 \frac{A}{B} \frac{C-B}{C... ...ad \mu^2 \equiv \frac{C-A}{A} \frac{C-B}{B} \, \omega_{30}^2 .$$ (17)

For the Earth the ratios $\omega_{10}/\omega_{30} \approx \omega_{20}/\omega_{30}$ describing the angle between the rotation axis and the figure axis are a few $0\hbox{$.\!\!^{\prime\prime}$}1 \approx 0.5 ~ 10^{-6}$ rad. Finally EKE in Eq. (2), must be solved with the left hand side being given by (16). The solution of the rearranged EKE

$$\begin{eqnarray} {\dot{\psi}} \sin {\theta} &=& \sin {\varphi} \, \omega_1(t) ... ...{\varphi}} &=& \omega_3(t) - {\dot{\psi}} \cos {\theta} \nonumber\end{eqnarray}$$ (18)

can be obtained by integration and gives with good approximation:

$$\begin{eqnarray} {\psi}(t) &=& \psi_0 - \frac{A_+}{\sin \theta_0} \cos \varphi_... ...eta_0 \{ A_+ \cos \varphi_{+} + A_- \cos \varphi_{-} \} \nonumber\end{eqnarray}$$ (19)

with $\varphi_\pm \equiv \varphi_0+[\omega_{30} \pm \mu] \, t$and constants

$$\begin{eqnarray} A_+ &\equiv& \frac{\omega_{10} + \omega_{20}}{2 \, (\omega_{30... ...}}{2 \, (\omega_{30} - \mu)} \approx -2.4~ 10^{-9} \, {\rm rad} .\end{eqnarray}$$ (20) (21)

$\psi_0, \theta_0$ and $\varphi_0$ are three integration constants. Neglected terms contribute with less than 1 $\mu$as or 5  10^-12 rad. Note that the result (19) does not give the usual nutation angles but the Euler angles between the figure axis system and an inertial system. How to extract the nutation angles from them will be explained in the Sect. 2.7.

2.5 Perturbation theory for the forced motion of the Earth

2.5.1 Solution of EDE with external torque

Now, the external luni-solar and planetary torques are included. It turns out that in order to satisfy the first two EDE's, (4), to our level of accuracy $\omega_3(t) = \omega_{30} = \rm const$ is sufficient. Then, by using the substitutions

$$\begin{eqnarray} \omega_1(t) &=& \omega_1^{\rm free} + (\omega_{\rm 1c}+t \, \... ...t}) \sin(\phi_i + \omega_i \,t) \\ \omega_3(t) &=& \omega_{30} ,\end{eqnarray}$$ (22) (23) (24)

where $\omega_1^{\rm free}, \omega_2^{\rm free}$are the free solutions according to (16) describing the Eulerian (=torque-free) motion of the Earth, one can satisfy the EDE by an appropriate choice of the eight constant coefficients $\omega_{\rm 1c}, \omega_{\rm 1ct},~\ldots$They depend on the tidal amplitudes C²¹_i, S²¹_i, tidal phases $\phi_i$ and frequencies $\omega_i$ and on the constants $\omega_{10}, \omega_{20}, \omega_{30}, \mu$ introduced in the last section. Further corrections on $\omega_3$ are not necessary for the computation of the nutations since they are smaller than $\omega_{30}$by a relative factor of at least 10^-9.

2.5.2 Solution of EKE with external torque

To solve EKE in the form (18), we have to substitute the solution (22-24) into their right hand sides. Rewriting the products of the trigonometric functions as sums and subtractions of their arguments, we obtain

$$\begin{eqnarray} \sin \theta \, \dot \psi =&& \bar a_3 + \bar a_4 \, t + \fra... ... \\ [1mm] \dot \varphi =&& \omega_{30} - \cos \theta \, \dot \psi\end{eqnarray}$$ (25) (26) (27)

with new constants $a_{1,~\ldots,~8}$and $b_{1,~\ldots,~8}$ simply given as sums and differences of the former $\omega_{\rm 1c}$, $\omega_{\rm 1ct}$, $\ldots$It turns out that the constants $a_3 = -b_1 =(\omega_{\rm 1s} + \omega_{\rm 2c})/2$ are the largest. The constant terms $\bar a_3, \bar a_4, \bar b_3$ and $\bar b_4$ denote the contributions from the K₁-tide for which $\varphi - \phi_i - \omega_i t$ vanishes and have to be separated from the other tides in the fourth lines since they determine the precession. The terms including $\omega_{10} \pm \omega_{20}$ result from the torque-free motion. Terms including arguments $(\varphi-\phi_i-\omega_i\,t)$will give the "real'' in- and out-of-phase nutation terms, while those including $(\varphi+\phi_i+\omega_i\,t)$will give the terms coming from the triaxiality of the Earth (also in- and out-of-phase, but the latter are negligible). The task is now to solve these differential equations by means of successive approximations.

2.5.3 Solution of EKE with external torque at the first order

In a first approximation, we neglect the second term in (27) and put $\sin \theta$ equal to $\sin \theta_0$ in (25), where $\theta_0$ is the constant mean obliquity at J2000. Then, by integration we find immediately:

$$\begin{eqnarray} \varphi_1 &=& \varphi_0 + \omega_{30} \, t \\ [1mm] \sin \the... ...omega_i^2} - \frac{b_6}{\Sigma\omega_i} \, t \right] \nonumber ,\end{eqnarray}$$ (28) (29) (30)

with the abbreviations  

$$\begin{eqnarray} \Delta \phi_i =& \phi_i - \varphi_0 , \qquad& \Sigma \phi_i = ... ...{30}, \qquad& \Sigma\omega_i = \omega_i + \omega_{30}. \nonumber \end{eqnarray}$$ (31)

The interpretation of the various terms is quite obvious. $\varphi_1$ describes (at the first order) the Earth's rotation around the figure axis. The first lines of (29) and (30) are secular contributions describing the precession. Next comes the contribution from the torque-free motion (the terms containing A₊ and A_-). Finally, the in- and out-of-phase short-periodic terms with periods approximately between 4 days and 18.6 years and the in- and out-of-phase terms coming from the triaxiality of the Earth with periods of about half a day occur. It should be mentioned that this solution also gives the Euler angles and not the nutation angles themselves (see below).

2.5.4 Solution of EKE with external torque at the second order

Now second order corrections necessary to obtain the desired accuracy will be computed. They need to be applied only to the first order terms being large enough, say 10 mas, i.e. to the precession and the in-phase nutation, but not to the free motion, the triaxiality terms and the out-of-phase nutation terms.

First, the correction for the angle $\varphi$ due to the neglected second term in (27) can be computed from  

$$\dot{\varphi}_2 = - \cos\theta \dot{\psi}$$ (32)

and leads to (taking into account for $\psi$ only the precession $\bar a_3$ and the nutation a_3i in longitude)

$$\begin{eqnarray} \varphi_2 &\approx &-\cos \theta_0 \, \psi_1 \approx \nonumbe... ...}}{\Delta\omega_i} \sin(\Delta \phi_i+\Delta\omega_i\,t) \right]. \end{eqnarray}$$ (33)

From now a subscript i on the constants $a_3, b_1,~\ldots$is written explicitly to show that they depend on the i-th tide. At the first degree of approximation we get: $a_{3i} = -b_{1i} \approx ( S^{21}_i / \omega_i)H_{\mathrm{dyn}}$,where S²¹_i is the tidal amplitude and $\omega_i$ the tidal frequency. An interesting feature of (33) is that the first contribution to $\varphi_2$,which is nothing else but the precession in right ascension, can be included into a modified Earth rotation speed $\tilde{\omega}_{30} = \omega_{30} - \cot \theta_0 \bar a_3$.In other words, $\omega_{30}$ and $\tilde{\omega}_{30}$ stand for the Earth's rotation speed in a precessing and inertial coordinate system, respectively. The other contribution in (33) is simply the first order difference between GMST and GAST which is the "equation of equinoxes''. Since it is a periodic contribution with a small amplitude an expansion only at the first order is possible.

The second order correction $\theta_2$ in obliquity is then computed. Again EKE have to be solved but now with $\varphi = \varphi_0 + \tilde{\omega}_{30} \, t + ({a_{3i}}/{\Delta\omega_i}) \allowbreak \sin(\Delta \phi_i+\Delta\omega_i\,t)$ instead of $\varphi = \varphi_0 + \omega_{30} \, t$in the arguments on the right hand side of (26). First, we get the solution (30) with $\omega_{30}$ being replaced by $\tilde{\omega}_{30}$(which also transforms $\Delta\omega_i = \omega_i - \omega_{30}$into $\Delta\tilde{\omega}_i = \omega_i - \tilde{\omega}_{30}$). Although this modified solution may be considered to contain "second order'' contributions, we still call these terms first order terms here. Second, we get a true second order term from the expansion with respect to the periodic correction in $\varphi_2$.By integrating

$$\begin{eqnarray} \dot\theta_2 &=& \sum_{i,j} - \cot \theta_0 \frac{a_{3j}}{\Del... ...mega_i} \, t) \sin (\Delta \phi_j + \Delta \tilde{\omega}_j \, t)\end{eqnarray}$$ (34)

and after some rearranging and using b_1i = - a_3i, we obtain

$$\begin{eqnarray} \theta_2 &=& \sum_{i,j} \frac{\cot \theta_0}{4} \frac{a_{3i}}... ...tilde\omega_i}+\Delta \tilde{\omega}_j]\,t) \biggr\} \nonumber .\end{eqnarray}$$ (35)

The second order correction $\psi_2$ in longitude is computed in a similar way. The corrections on the l.h.s. of (25) come from the first order nutation in obliquity denoted by $\bar \theta_1$(the precession in obliquity is neglected here) and a correction to the first order nutation $\psi_1$

$$\begin{eqnarray} \sin \theta \, \dot{\psi} &\approx& \sin(\theta_0 + \bar \the... ...theta_1 \cos \theta_0 \dot{\psi}_1 + \sin \theta_0 \dot{\psi}_2 .\end{eqnarray}$$ (36)

The underlined term is of first order, while the other two terms describe second order contributions. There is also one term from the expansion of the right hand side of (25) with respect to the periodic correction in $\varphi_2$ (see Eq. (33)). Omitting the detailed derivation, we find in total three contributions. First, the first order solution (29) with all $\omega_{30} [\Delta\omega_i]$ being replaced by $\tilde{\omega}_{30} [\Delta\tilde{\omega}_i]$.Second, what we can call the "precession on nutation'' effect in longitude (coming from the $\bar \theta_1$term in (36))  

$$\sin \theta_0 \, \psi_2^{(1)} = \sum_{i,j} - \cot \theta_0... ...ega_i}^2} \sin (\Delta \phi_i + {\Delta\tilde\omega_i} \, t) .$$ (37)

Third, the nutation term at the second order in longitude which comes from the remaining second order terms:

$$\begin{eqnarray} \sin \theta_0 \, \psi_2^{(2)} &=& \sum_{i,j} \frac{\cot \thet... ...lta \phi_j+[{\Delta\tilde\omega_i}+\Delta \tilde{\omega}_j]\,t) . \end{eqnarray}$$ (38)

Summing up the results from the previous subsections, the complete analytical solution for the Euler angles $\psi, \theta, \varphi$ between the figure axis system and an inertial system for the motion of the Earth according to a (mainly luni-solar) torque of the form (12) has been computed for our level of accuracy. All of the second order terms are proportional to the product of two tidal amplitudes by definition.

2.6 Coordinate systems for precession and nutation

The definition of the precessional and nutational quantities can be seen in Fig. 1.

  

Figure 1: Coordinate systems and variables for precession and nutation

The subscripts D and F stand for date and fixed in the following. Important coordinate systems are:

1.

the mean equator of date, EQ$_{\rm D}$,

2.

the mean fixed equator J2000, EQ$_{\rm F}$,

3.

the mean ecliptic of date, EC$_{\rm D}$,

4.

the mean fixed ecliptic J2000, EC$_{\rm F}$

5.

the system related to the figure axis of the Earth, denoted by $\sum_{\displaystyle \pi}$ here, which is identical to the true equator of date counting the origins from the axis of the smallest moment of inertia A or - by using the axisymmetric approximation for the Earth and a constant rotation which can be included in $\varphi_0$ (see Eq. (19)) - from Greenwich Meridian.

Williams (1994) presents four methods to make the transformation from the mean fixed equator and equinox J2000 (e.g. used in JPL ephemerides DE200) via the mean equinox of date $\rm \gamma_D^{\mathrm{mean}}$to the true equator and equinox of date. $R_i(\alpha)$ denotes a rotation with angle $\alpha$ around the i-axis. For the meaning of the various variables see also Williams (1994) and Fig. 1.

1$^\mathrm{st}$method: 

($\gamma_{\rm F}^{\mathrm{mean}} \to Q \to \gamma_{\rm D}^{\mathrm{mean}} \to \gamma_{\rm D}^{\mathrm{true}}$)

$$\begin{eqnarray} \underbrace{R_1(-\varepsilon_A-\Delta\varepsilon_{\rm D}) \; R... ...(90^\circ\!-\zeta_{\rm A})} _{\displaystyle\cal{P}} , \nonumber \end{eqnarray}$$ (39)

where ${\cal N}$ and ${\cal P}$ are the usual nutation and precession matrix, respectively.

2$^\mathrm{nd}$method: 

($\gamma_{\rm F}^{\mathrm{mean}} \to \gamma_1 \to \gamma_{\rm D}^{\mathrm{mean}} \to \gamma_{\rm D}^{\mathrm{true}}$)

$$\begin{eqnarray} R_1(-\varepsilon_{\rm A}-\Delta\varepsilon_{\rm D}) R_3(-\Delt... ...ga_{\rm A}) R_3(-\psi_{\rm A}) R_1(\varepsilon_0) . \nonumber \end{eqnarray}$$ (40)

3$^\mathrm{rd}$method: 

($\gamma_{\rm F}^{\mathrm{mean}} \to N \to \gamma_{\rm D}^{\mathrm{mean}} \to \gamma_{\rm D}^{\mathrm{true}}$)

$$\begin{eqnarray} R_1(-\varepsilon_{\rm A}-\Delta\varepsilon_{\rm D}) R_3(-\Pi_... ...\pi_{\rm A}) \; R_3(\Pi_{\rm A}) \; R_1(\varepsilon_0) , \nonumber\end{eqnarray}$$ (41)

where $\Pi_{\rm A}+p_{\rm A}=\sigma_{\rm A}+\eta_{\rm A}$, see Williams (1994).

4$^\mathrm{th}$method: 

($\gamma_{\rm F}^{\mathrm{mean}} \to \gamma_2 \to \gamma_{\rm D}^{\mathrm{mean}} \to \gamma_{\rm D}^{\mathrm{true}}$)  

$$R_1(-\varepsilon_{\rm A}-\Delta\varepsilon_{\rm D}) R_3(-... ...{\rm D}) R_1(\varepsilon_{\rm A}^\prime) R_3(\xi_{\rm A}) .$$ (42)

And we have another method for later purpose in Sect. 2.8:

5$^\mathrm{th}$method: 

( $\gamma_{\rm F}^{\mathrm{mean}} \to \gamma_3 \to \gamma_{\rm D}^{\mathrm{true}}$)

$$\begin{eqnarray} R_3(\chi_{\rm A}+\Delta\chi_{\rm F}) R_1(-\varepsilon_{\rm A}... ...3(-\psi_{\rm A}-\Delta\psi_{\rm F}) R_1(\varepsilon_0) . \nonumber\end{eqnarray}$$ (43)

These transformations are needed to relate the Euler angles computed in the last subsection to the nutation quantities used in the nutation theories.

2.7 The relation between the Euler angles and the nutation and precession variables with respect to the mean fixed ecliptic

 From the previous subsection and with the help of Fig. 1 we can extract the nutation quantities from the Euler angles shown above. Since the Euler angles of the figure axes system directly transform from $\sum_{\displaystyle \pi}$ to the inertial system represented by the mean fixed ecliptic J2000 (and vice versa) one finds by comparing  

$$R_3({\phi}) R_1(-{\theta}) R_3(-{\psi}) R_1(\varepsilon_0)$$ (44)

with (43), or regarding Fig. 1, that the following relations hold

$$\begin{eqnarray} {\psi}(T)&\equiv\;\psi_{\rm F}^{\mathrm{fig}}=&\psi_{\rm A}+ \... ...fig}}=&\chi_{\rm A}+ \Delta\chi_{\rm F}+{\rm GAST}(T) \nonumber . \end{eqnarray}$$ (45)

$\psi_{\rm A}, \omega_{\rm A}, \chi_{\rm A}$ are respectively the precession in longitude, in obliquity and the planetary precession. T is the time in Julian centuries measured from J2000. Thus there is a very simple relationship between the Euler angles and the nutation angles $\Delta\psi_{\rm F}, \Delta\varepsilon_{\rm F}$referred to the mean fixed ecliptic J2000. However, traditionally the nutation series are referred to the mean ecliptic of date so the transformation described in the next section is needed.

2.8 The relation between the nutation angles referred to the mean fixed ecliptic and those referred to the mean ecliptic of date

 The definition of the nutation angles with respect to the mean ecliptic of date can also be seen in Fig. 1. Thus, $\Delta\psi_{\rm D}$ and $\Delta\varepsilon_{\rm D}$ transform from mean EC$_{\rm D}$ to $\sum_{\displaystyle \pi}$.

Comparing the transformation (43) to anyone of the four methods (39-42) and solving for $\Delta\psi_{\rm D}$ and $\Delta\varepsilon_{\rm D}$we obtain (using the numerical values for the other precessional quantities given in Simon et al. 1994) 0pt

$$\begin{eqnarray} \Delta\psi_{\rm D}= &\phantom{+} \Delta\psi_{\rm F} \, (1 & + 0... ... (1& - 0.000\,000\,003 \, T^2 - 0.000\,000\,001\, T^3). \nonumber\end{eqnarray}$$ (46) (47)

Therefore, only the small time-dependent part of the nutation angles is affected by this transformation while the constant part remains unchanged. Whether the higher powers of T produce any significant terms is discussed in Sect 3.5.

2.9 The geodesic precession and nutation

For completeness the relativistic contribution to precession and nutation is also given here. Using the results of Fukushima (1991) or of Damour et al. (1993) and the numerical expressions of Simon et al. (1994) we obtain for the additional terms in longitude (there are no terms in obliquity)

$$\begin{eqnarray} \Psi &=& \frac{3 G M_{\odot}}{2 \, c^2 a} \left[ (1 + e^2) l_... ... s} + \frac{9}{4} e^2 \sin 2l_{\rm s} \ldots \right], \nonumber \end{eqnarray}$$

where a, e, and $l_{\rm s}$ are respectively the semi-major axis, eccentricity, and mean anomaly of the Earth's orbit.

$$\begin{eqnarray} \Psi &=& 0\hbox{$.\!\!^{\prime\prime}$}019\,062 + 1\hbox{$.\... ...} + 0\hbox{$.\!\!^{\prime\prime}$}000\,001\,9 \, \sin 2l_{\rm s}.\end{eqnarray}$$ (48)

The secular term is the geodesic precession $p_{\rm g}$ and the periodic terms are called the geodesic nutation $\Delta \psi_{\rm g}$. Considering that the kinematically non-rotating reference system is to be adopted these terms should be subtracted from the precession and from the nutation series, respectively.

2.10 The other axes - the Oppolzer terms

The theories of nutation mentioned in Sect. 2.1 compute the nutation of the angular momentum axis and give the transformation to obtain the nutation of the other two axes (figure axis and rotation axis) by adding the so called Oppolzer terms. The nutation theory presented here which starts from the tidal potential gives the nutation of the figure axis. To get the nutation for the other two axes we have therefore to use some combination of the classical Oppolzer terms

$$\begin{eqnarray} {\Delta\psi}_{\mathrm{fig} \to \mathrm{rot}} &=& {\Delta\psi}_... ...\mathrm{Opp}} - {\Delta\varepsilon}_{\mathrm{fig}}^{\mathrm{Opp}} \end{eqnarray}$$ (49)

(see K77 and ZG89). Using some geometrical relations, and vanishing the terms under our threshold, the result derived in Hartmann (1996) is

$$\begin{eqnarray} {\Delta\psi}_{\mathrm{fig} \to \mathrm{rot}} &=& - \frac{{\De... ...er \\ & & - \frac{p_A}{\tilde{\omega}_{30}} \sin \varepsilon_0 ,\end{eqnarray}$$ (50) (51)

where $\psi_{\mathrm{in}}$ and $\varepsilon_{\mathrm{in}}$ denote the in-phase terms only.

To get the nutation for the angular momentum axis we have to replace $\tilde{\omega}_{30}$ in Eqs. (50-51) by $C/A \cdot \tilde{\omega}_{30}$.The constant offset in obliquity due to the precession in longitude $p_{\rm A}$ (-8711.90 $\mathrm \mu$as and -8683.38 $\mathrm \mu$as for the rotation axis and angular momentum axis, respectively) is well-known, see e.g. Woolard (1953) or Capitaine et al. (1985). For another way to derive the Oppolzer terms, except for the constant offset, see ZG89.