GENERAL THEORY OF THE ROTATION OF THE NON-RIGID EARTH AT THE SECOND ORDER. I. THE RIGID MODEL IN ANDOYER VARIABLES

Present astrometric techniques allow for determining the transformation between the terrestrial and celestial reference frames with great accuracies, demanding a constant improvement of Earth rotation modeling. From the perspective of theoretical investigations, it reflects the need for developing complete models that necessarily must incorporate second-order effects, which are of two different types accordingly they come from direct new small contributions, or from crossed combinations of first-order terms, that is to say, non-linear second-order terms in the sense of perturbation theories. Some of the current studies (Dehant et al. 2008) provide contributions at the microarcsecond (μas) level. At this truncation level, it is also indispensable to analyze the internal consistency of the existing theories, examining the way in which the second-order terms are mathematically tackled. In this sense, let us point out that the official precession–nutation model IAU 2000A (adopted at the 24th General Assembly of the International Astronomical Union), based on the transfer functions of Mathews et al. (2002), is strictly speaking only valid at the first order. It is due to the fact that only at the first order the approximation of the external gravitational torque perturbing the Earth rotation is independent of the Earth model. Therefore, with this formulation, it is not possible to incorporate in a mathematical consistent and systematic way the above-mentioned non-linear second-order contributions.

A clear evidence of the importance of these kinds of effects was reported in Ferrándiz et al. (2004), where for the first time it was shown that the internal structure of the Earth does affect the precession rate in longitude in a non-negligible quantity. That discovery was the consequence of applying the Hamiltonian methods to study the second-order effects in the precessional motion of a two-layer Earth model composed of a rigid mantle and a fluid core, the so-called Poincaré model. Bearing in mind these considerations, we are developing a second-order analytical theory of the rotational motion of the non-rigid Earth. The construction of a theory of these characteristics is a formidable task, since to account fully for the second-order effects it would be necessary to consider, besides the spin–spin coupling, the interactions between the orbital motion of the Moon and the rotational motion of the non-rigid Earth. In addition, it also would be necessary to discuss the possible influence of the non-sphericity of the Moon, i.e., the coupling between the Earth and the Moon rotations, and of the interaction with the motions of the Sun.

To tackle all these inherent complexities in a way that allows obtaining a solution, we have divided the problem into different parts. In the first one, we will analyze the effects of the non-rigidity of Earth in the spin–spin coupling, since it is expected that it plays a more important role than in the case of the spin–orbit coupling. Then, in a second stage, we will examine the interesting and almost unexplored question about how the non-rigidity of the Earth affects the interaction between the orbital motion of the Moon and the rotational motion of the Earth.

Here, we start the first stage of the above described investigation by considering a general and complete second-order theory of the rotation of the non-rigid Earth concerned with spin–spin coupling. From now on, and for the sake of conciseness, we will refer simply to the rotation of the non-rigid Earth at the second order but, as we have mentioned, it must be understood that here we are only focused on the spin–spin coupling aspect. The theory has a general character since, as opposed to other recent investigations, it comprises not only the precession and the nutation but also the whole rotational motion. In this way, it extends previous first-order Hamiltonian studies of the non-rigid Earth (e.g., Getino & Ferrándiz 1995, 2001; Escapa et al. 2001), providing a unified way to work out analytical formulae that model the nutation–precession, the polar motion, and the length of day (LOD) of a non-rigid Earth model. On the other hand, the theory is complete in the sense that the whole expression of the considered rotational Hamiltonian of the system is worked out, contrary to other treatments that assume a partial expression of the same Hamiltonian.

In this regard, let us point out that although other approaches are probably possible, we keep working within the Hamiltonian framework since in our opinion this approach is specially suitable to tackle analytically the problem due to, among other advantages, the existence of canonical sets particularly well suited to describe the orbital and rotational motions (Delaunay and Andoyer variables), and the possibility of applying canonical perturbation theories.

Let us recall that, even in the case of a rigid Earth model, it is not a direct task to develop a second-order analytical theory of the rotational motion of the Earth. It is due to the intrinsic complexities of second-order solutions and to the six-dimensional structure of the phase space of this dynamical system. When generalizing this theory to non-rigid models the difficulties grow up exponentially, since we increase the dimension of the phase space and, what is more important, we have to incorporate more realistic mechanisms into the models. Furthermore, the inclusion of dissipation processes in the framework of a second-order theory makes its formulation more cumbersome, since one should consider general, non-canonical, perturbation methods or transforming the system of differential equations describing the motion into a canonical one by doubling the dimension of the original phase space.

Under these circumstances, it has been necessary to partially leave the first-order canonical approaches, since they are not suitable for second-order theories, and to develop a new formalism that made the treatment of the second-order contributions easier. This formalism takes advantage of the use of a matrix formulation that reduces the difficulties in the analytical computation of the new contributions, as well as enables the inclusion of additional effects not yet modeled in a systematic way. In order to underline the essential points of the new formalism and to make them intelligible, we will present the results of the investigation concerning the spin–spin coupling on the non-rigid Earth in different papers. In subsequent articles, we will report the fundamentals of the new analytical formalism, its application to studying the rotational motion of the non-rigid Earth and a comparison with the methods used in the precession–nutation model IAU 2000A, showing that the application of the transfer function to the rigid second-order contributions is completely misleading.

In this paper, we present a complete second-order analytical treatment of the rotation of the rigid Earth relative to rotation–rotation interaction. The reason to begin studying the rotation of the non-rigid Earth at the second order with a rigid Earth model is due to the fact that the rigid Earth is a kind of benchmark for non-rigid models, in the sense that it is expected that when reducing the non-rigid Earth treatments to the rigid ones we recover the rigid Earth results. Obviously, it does not guarantee the complete correctness of the non-rigid formulation, but it is an indication in its favor. It has been the case, for example, of the first-order Hamiltonian non-rigid analytical theories, which allow recovering the classical analytical results of the rigid Earth given in Kinoshita (1977).

With respect to the second-order Hamiltonian theories in Souchay et al. (1999; see also Kinoshita & Souchay 1990), a complete reconstruction of the nutation of the rigid Earth with a truncation level of 0.1 μas was carried out. This Hamiltonian investigation incorporated the second-order terms as one fundamental aspect of the theory. In particular, these authors considered two kinds of second-order effects: one due to the interaction of the Earth rotation with itself, the spin–spin coupling, and the other one due to the interaction of the Earth rotation with the orbital motion of the Moon, the spin–orbit coupling. However, the treatment of the spin–spin coupling presents some shortcomings. In contrast to the treatment given in Kinoshita (1977), analytical expressions for the different aspects of the rotational motion in their final, simplified form which had allowed users to form their theoretical interpretation or an independent numerical evaluation in a direct manner were not provided. Instead, they gave some intermediate formula from which the reported numerical values could be derived with the use of proper software procedures. In addition, this theory was not general since the whole rotational motion of the Earth was not discussed, but only the motion of the equatorial plane by assuming that when computing the second-order effects, the rotational angular momentum and the axis of the maximum principal moment of inertia can be considered coincident. There are also other Hamiltonian treatments of the rotation of the rigid Earth based on a specific symbolic processor (Navarro & Ferrándiz 2002). Although this approach has provided important corrections to the Souchay et al. (1999) tables, as in the case of planetary nutations, the second-order effects are computed partially and without providing any final analytical formulae.

Hence, due to this lack of a complete analytical treatment, we had developed a second-order theory of the rotation of the rigid Earth which provides explicit analytical expressions for the precession–nutation, the polar motion, and the LOD, arising from the spin–spin coupling. As we have remarked, these final analytical formulae will be very valuable in the non-rigid case and will allow determining the contribution of the non-rigidity in the second-order terms. We have performed this investigation in the same Hamiltonian framework as that introduced by Kinoshita (1977), aiming to comprehensively extend this cornerstone investigation at the second order. However, there are some significative differences with respect to this work and, specially, with respect to the subsequent treatments given to the second-order terms in Kinoshita & Souchay (1990) and Souchay et al. (1999). For example, one is concerned with the splitting of the Hamiltonian to apply the perturbation methods: our approach allows a more direct computation of the generating function, since only the kinetic energy of the model is considered in the unperturbed Hamiltonian. Another important difference refers to distinguishing, also in the second order, the rotational angular momentum from the axis of the maximum principal moment of inertia as it is done in the first-order case, that is to say, to consider a complete treatment. In this way, we remove the simplification assumed in Souchay et al. (1999), a consequence being the appearance of some new contributions in the motion of the equatorial plane whose numerical value is greater than the threshold of 0.1 μas established in Souchay et al. (1999). In addition, these terms play a relevant role in the case of the Poincaré model due to the fluid resonance, as was shown in Ferrándiz et al. (2004).

The structure of the paper is as follows. In Section 2, after introducing the different nomenclature for the Earth rotation problem, we review the fundamentals of the Hamiltonian formalism to mathematically model the rotation of the rigid Earth, providing the equations of motion of the system and giving a sketch of the second-order perturbation method (Hori 1973). In particular, we consider that the Earth rotation motion is perturbed by the main term V in the expansion of the Earth geopotential, J₂ term, without performing any simplification in its expression as a function of the Andoyer canonical variables. Namely, as in the first-order case, we assume that the rotational angular momentum and the axis of the maximum principal moment of inertia are not coincident, but form a definite angle σ different from 0. In addition, the Hamiltonian of the system incorporates a complementary term E due to the fact that the Earth rotational motion is referred to the non-inertial moving ecliptic. In Section 3, we apply the canonical theory previously presented to study the rotation of the rigid Earth at the second order. In contrast to other Hamiltonian investigations, we treat both V and E as first-order perturbations over the unperturbed Hamiltonian of the system, given by its kinetic energy T. In this way, the second-order contributions are a consequence of the interaction between the J₂ term and itself, and the interaction between E and the J₂ term. Once the Hamiltonian of the system is divided into those two parts according to their order of magnitude, we implement the perturbation algorithm determining the zeroth-, first-, and second-order terms of the transformed Hamiltonian, as well as its associated generating function of the first and second order. This computation is mainly based on the evaluation of different Poisson brackets and their splitting in periodic and secular parts. From these expressions, it is possible to determine the evolution of any function of the canonical variables of the system at the second order. In this sense, the theory presented here is a general one, since it is applicable to studying any aspect of the rotational motion.

In particular, in Section 4, we determine the motion of the equatorial plane, the time evolution of the angle ϕ, essential in determining LOD variations, and the polar motion. Specifically, we provide the final, simplified analytical formulae of functions that allow us to obtain the amplitudes of the short-period part of the motion, as well as of the rates characterizing the secular part of the motion. Some of these functions emerge as a consequence of considering σ ≠ 0 in the second-order developments. In Section 5, we summarize the origin of the second-order contributions according to the different parts of the first-order interactions, giving a comprehensive comparison with the results obtained by Souchay et al. (1999) for the spin–spin coupling. This comparison focuses on the differences originated from removing the simplification σ = 0 assumed in Souchay et al. (1999), as well as from the different splitting of the Hamiltonian of the system taken in both approaches. We have restricted ourselves to comparing the motion of the equatorial plane, since it is the only aspect of the Earth rotational motion considered in Souchay et al. (1999). We find a good agreement among the terms coming from the simplification σ = 0, although this agreement has had to be established numerically, since in Souchay et al. (1999) equivalent analytical expressions to those ones obtained in this work are not given. However, we also find some additional terms, which come from taking σ ≠ 0, whose numerical value is greater than the truncation level of 0.1 μas established in Souchay et al. (1999). Specifically, we have found three terms in longitude greater or equal than 0.1 μas and nine terms in obliquity, as well as a new contribution to the precession rate in longitude of 0.692 mas cy⁻¹. These terms should be considered to reach the above-mentioned truncation level. Finally, in order to facilitate the reading of the article, we have included five appendices and different tables where the explicit analytical formulae of the second-order contributions, some intermediate results, and the different numerical parameters used in this work to compute some second-order contributions are presented.

2.1. Description of the Rotational Motion

The Earth's rotational motion around its center of mass is defined by giving the temporal dependence of a rotational operator that connects two reference frames with a common origin: a celestial reference frame, inertial or non-inertial, with a terrestrial reference frame, or body-fixed frame, which is attached to the Earth in some prescribed way. From a mathematical point of view, these two reference frames are described, respectively, with the help of two orthogonal coordinate triads $\left(\vec{e}_{X},\vec{e}_{Y},\vec{e} _{Z}\right)$ and $\left(\vec{e}_{x},\vec{e}_{y},\vec{e}_{z}\right)$. The rotational operator is represented by a rotation matrix whose entrances are regular functions of time.

There are many possibilities to characterize the rotation matrix (see Fukushima 2008 for a recent review), the Euler angles in the sequence 1–3–1 being the most usual choice in Earth rotation studies; we denote them as ψ (longitude), θ (obliquity), and ϕ (see Figure 1 in Appendix A). The time evolution of these angles is a consequence of the different interactions acting on the Earth, as well as of its internal mechanical properties, and can be determined by solving the differential equations of the Earth's rotational dynamics. The solution of these equations provides the values of the Euler angles and their first derivatives $\dot{\psi }$, $\dot{\theta }$, and $\dot{\phi }$ for each instant, therefore giving a complete knowledge of the Earth rotational motion.

Zoom In Zoom Out Reset image size

Alternatively, the rotational motion is also determined by giving the values of the Euler angles and the components ω_x, ω_y, ω_z of the angular velocity vector with respect to the terrestrial reference frame at each instant, since those components are functions of the Euler angles and their first derivatives through the Euler kinematical equations (Kinoshita 1977). This description characterizes all the different classical aspects of the Earth rotational motion, namely:

• 1.  
  The angles ψ and θ determine the components of the z-axis of the terrestrial reference frame in the celestial reference frame, since

  $$\begin{eqnarray} \vec{e}_{z}(t)&=&(\sin \theta (t)\sin \psi (t),\nonumber\\ &&-\sin \theta (t)\cos \psi (t), \cos \theta (t)) ^{T}, \end{eqnarray} \tag{ 1 }$$

  the superscript T denoting the transpose. Usually, the z-axis of the terrestrial reference frame has the same direction and sense as the axis of the maximum principal moment of inertia; in this case, it is named figure axis, the equator of the Earth being the perpendicular plane to this axis. Hence, the angles ψ and θ can also be viewed as the longitude of the node and the inclination of the equatorial plane. The functions of time θ(t) and ψ(t) are expressed as quasi-polynomials in t. The secular (non-bounded) and long-period part of these functions is known as precessional motion or simply precession, whereas the short-period part is known as nutational motion, or nutation, and represents a motion of quasi-periodic nature.
• 2.  
  The components ω_x and ω_y provide the polar motion, which describes the position of the angular velocity vector in the xy-plane of the terrestrial reference frame or, in other words, the difference between the figure axis and the angular velocity vector.
• 3.  
  The time evolution of ϕ as well as of the component ω_z are the basis of measurement of time by the Earth rotational motion, that is to say, of the sidereal time and LOD variations.

2.2. Hamiltonian Formalism

Although it is possible to develop a dynamical theory of the rotational motion of the Earth in terms of the Euler angles by means of the Euler or Lagrangian equations (e.g., Woolard 1953; Bretagnon et al. 1998), the treatment of this problem within the framework of Hamiltonian mechanics has been extremely productive. As a matter of fact, some of the most accurate theories for the rotation of rigid and non-rigid Earth models (Souchay et al. 1999; Getino & Ferrándiz 1995, 2001; Escapa et al. 2001) have employed this approach, following the pioneering work of Kinoshita (1977). Next, we sketch the basic guidelines of Hamiltonian theory for the rotation of a rigid Earth, referring the reader to the investigation of Kinoshita (1977) for a comprehensive treatment.

Let us consider a rigid Earth model whose principal moments of inertia are denoted as A, B, and C, the axis of each moment being coincident with the x, y, and z axes of the terrestrial frame, whose origin is located at the center of mass of the Earth. Since the Earth is practically an axial symmetric body, we will consider that the equatorial moments of inertia are equal, that is to say, we will take A = B < C. From the perspective of the Hamiltonian methods, the rotational behavior of this system is described through a canonical set of variables, being the most expedient choice the Andoyer canonical set (Andoyer 1923) denoted as Λ, M, N, λ, μ, ν. The variables Λ, M, N are the canonical momenta, whereas λ, μ, ν are the canonical coordinates. The canonical momenta have a clear dynamical meaning (Kinoshita 1977): M is the modulus of the angular momentum vector $\vec{L}$, and Λ and N are, respectively, its projections on the $\vec{e} _{Z}$ and $\vec{e}_{z}$ axes. Hence, it is also useful to introduce the auxiliary variables σ and I

$$\begin{equation} \cos I=\frac{\Lambda }{M},\quad \cos \sigma =\frac{N}{M}, \end{equation} \tag{ 2 }$$

in such a way that I is the angle between the angular momentum vector and the Z-axis of the celestial reference frame and σ is the angle between the angular momentum vector and the figure axis (see Figure 1 in Appendix A). Further details of the geometrical and dynamical interpretation of the Andoyer set can be found in Kinoshita (1977) and in Efroimsky & Escapa (2007).

2.3. The Equations of Motion

The equations of motion of the rotation of the rigid Earth are derived from the Hamiltonian of the system, which for our model has the form

$$\begin{equation} \mathcal {H}=T+V+E. \end{equation} \tag{ 3 }$$

The first term is the rotational kinetic energy, the second one is the potential energy, and the third one is introduced when considering a non-inertial celestial reference frame (Kinoshita 1977). The kinetic energy is given by

$$\begin{equation} T=\frac{1}{2A}(M^{2}-N^{2})+\frac{1}{2C}N^{2}. \end{equation} \tag{ 4 }$$

The potential energy is due to the gravitational perturbation exerted by the Moon and the Sun. Since in this investigation we are working out a second-order theory, we will consider the crossed interactions coming from the complete main term in the Earth geopotential expansion (J₂ term). From Kinoshita (1977), this terms turns out to be of the form

$$\begin{eqnarray} V &=&{\sum _{p=S,M} }k_{p}^{\prime }\sum _{i}\Bigg[ \frac{1}{2}(3\cos ^{2}\sigma -1)B_{i}\cos \Theta _{i}\nonumber\\ &&- \frac{1}{2}\sin 2\sigma \sum _{\tau =\pm 1}C_{i,\tau }\cos (\mu -\tau \Theta _{i})\nonumber \\ &&+\frac{1}{4}\sin ^{2}\sigma \sum _{\tau =\pm 1}D_{i,\tau }\cos (2\mu -\tau \Theta _{i})\Bigg], \end{eqnarray} \tag{ 5 }$$

where the subindex p stands for the perturbing body (S = Sun, M = Moon). The parameter k'_p = 3κ²m_p(C − A)/a³_p characterizes the gravitational intensity generated by the external body. The coefficients B_i, C_i,τ, and D_i,τ, where τ can take the values +1 or −1, depend on the auxiliary variable I and are of the form:

$$\begin{eqnarray} B_{i} &=&-\!\frac{1}{6}(3\cos ^{2}I-1)A_{i}^{0)}-\frac{1}{2}\sin 2IA_{i}^{1)}- \frac{1}{4}\sin ^{2}IA_{i}^{2)}, \nonumber \\ C_{i,\tau } &=&-\!\frac{1}{4}\sin 2IA_{i}^{0)}+\frac{\tau }{4}\sin I(1+\tau \cos I)A_{i}^{1)}\nonumber\\ && +\frac{1}{2}(1+\tau \cos I)(-1+2\tau \cos I)A_{i}^{2)}, \nonumber \\ D_{i,\tau } &=&-\!\frac{1}{2}\sin ^{2}IA_{i}^{0)}+\tau \sin I(1+\tau \cos I)A_{i}^{1)}\nonumber\\ && -\frac{1}{4}(1+\tau \cos I)^{2}A_{i}^{2)}. \end{eqnarray} \tag{ 6 }$$

As for the argument Θ_i, we have that

$$\begin{equation} \Theta _{i}=m_{1i}l_{M}+m_{2i}l_{S}+m_{3i}F+m_{4i}D+m_{5i}\Omega, \end{equation} \tag{ 7 }$$

where l_M, l_S, F, D, and Ω are functions of the Delaunay variables of the Moon and the Sun. Within our level of approximation, we can assume that

$$\begin{equation} \frac{d\Theta _{i}}{dt}=n_{i}, \end{equation} \tag{ 8 }$$

with n_i being constant. Let us stress that the canonical variable λ is implicitly contained in Ω through

$$\begin{equation} \Omega =\Omega _{0}-\lambda \Rightarrow \frac{\partial \Theta _{i}}{\partial \lambda }=-m_{5i,} \end{equation} \tag{ 9 }$$

where Ω₀ is the mean longitude of the Moon referred to the origin of longitude on the ecliptic of date. The numerical values of the coefficients A^j)_i, Θ_i, and n_i, as well as of the list of the five integer numbers (m_1i, m_2i, m_3i, m_4i, m_5i) associated to each value of the index i, depend on the orbital theories of the Moon and the Sun, and are given in Kinoshita (1977) and updated in Kinoshita & Souchay (1990; also see Tables 8 and 9).

The term E is due to the non-inertial character of the celestial reference frame. This frame is coincident with the ecliptic of date, since in this way the development of the potential energy is considerably simplified (Kinoshita 1977). This reference is moving with respect to the fixed ecliptic of epoch in a known manner; in particular, this motion is described by giving the temporal dependence of two angles (see Kinoshita 1977)

$$\begin{equation} \begin{array}{c} \sin \pi _{1}\sin \Pi _{1}=pt+p^{\prime }t^{2}+p^{\prime \prime }t^{3} \nonumber \\ \sin \pi _{1}\cos \Pi _{1}=qt+q^{\prime }t^{2}+q^{\prime \prime }t^{3} \end{array}\!. \end{equation} \tag{ 10 }$$

The expression of E (Kinoshita 1977; Souchay et al. 1999) can be cast as

$$\begin{equation} E=\Lambda e_{1}+M\sin I (e_{2}\cos \lambda +e_{3}\sin \lambda), \ \end{equation} \tag{ 11 }$$

 e_i being functions of time of the form

$$\begin{equation} \left. \begin{array}{l} e_{1}(t)=\left(1-\cos \pi _{1}\right) \displaystyle \ \frac{d\Pi _{1}}{dt}, \nonumber \\ \nonumber \\ e_{2}(t)=\sin \pi _{1}\cos \Pi _{1}\displaystyle \ \frac{d\Pi _{1}}{dt}+\sin \Pi _{1}\displaystyle \ \frac{d\pi _{1}}{dt}, \nonumber \\ \nonumber \\ e_{3}(t)=\sin \pi _{1}\sin \Pi _{1}\displaystyle \ \frac{d\Pi _{1}}{dt}-\cos \Pi _{1}\displaystyle \ \frac{d\pi _{1}}{dt}. \end{array} \right. \end{equation} \tag{ 12 }$$

Since π₁ ≪ 1, we can approximate the former expressions obtaining

$$\begin{eqnarray} e_{1}(t) &=&0+O\big(\pi _{1}^{2}\big), \nonumber \\ e_{2}(t) &=&p+2p^{\prime }t+3p^{\prime \prime }t^{2}+O\big(\pi _{1}^{2}\big), \nonumber \\ e_{3}(t) &=&-\!q-2q^{\prime }t-3q^{\prime \prime }t^{2}+O\big(\pi _{1}^{2}\big). \end{eqnarray} \tag{ 13 }$$

In view of the smallness of magnitudes, the p', p'',  q', and q'' (Kinoshita 1977), and considering that in this investigation we are focusing on second-order effects, we will take e₁(t), e₂(t), and e₃(t) as constant parameters.

Once the Hamiltonian of the system has been constructed, we can obtain the time evolution of the Andoyer set by solving the Hamilton equations

$$\begin{equation} \frac{dp}{dt}=-\frac{\partial \mathcal {H}}{\partial q},\quad \frac{dq}{dt}=\frac{ \partial \mathcal {H}}{\partial p}, \end{equation} \tag{ 14 }$$

where p and q stand for each pair of the conjugated canonical variables (M, μ), (N, ν), and (Λ, λ). In this way, the behavior of the Earth's rotational motion can be completely modeled, since any aspect of the Earth's rotation can be mathematically expressed with the help of a function f(p, q) depending on the Andoyer canonical variables. In turn, the solution of Equation (14) provides the time evolution of the function f through

$$\begin{equation} \frac{df}{dt}=\left\lbrace f;\mathcal {H}\right\rbrace, \end{equation} \tag{ 15 }$$

where {−; −} denotes the Poisson bracket in the Andoyer canonical set, explicitly

$$\begin{eqnarray} \left\lbrace f;\mathcal {H}\right\rbrace &=&\frac{\partial f}{\partial \mu }\frac{ \partial \mathcal {H}}{\partial M}-\frac{\partial \mathcal {H}}{\partial \mu } \frac{\partial f}{\partial M}+\frac{\partial f}{\partial \nu }\frac{\partial \mathcal {H}}{\partial N}-\frac{\partial \mathcal {H}}{\partial \nu }\frac{ \partial f}{\partial N} \nonumber \\ &&+\frac{\partial f}{\partial \lambda }\frac{\partial \mathcal {H}}{\partial \Lambda }-\frac{\partial \mathcal {H}}{\partial \lambda }\frac{\partial f}{ \partial \Lambda }. \end{eqnarray} \tag{ 16 }$$

Because both the function f and the Hamiltonian $\mathcal {H}$ may contain the auxiliary variables I and σ, it is useful to consider the relationships

$$\begin{eqnarray} \frac{\partial }{\partial M} &=&\left(\frac{\partial }{\partial M}\right) + \frac{\cos \sigma }{M\sin \sigma }\frac{\partial }{\partial \sigma }+\frac{ \cos I}{M\sin I}\frac{\partial }{\partial I}, \nonumber \\ \frac{\partial }{\partial N} &=&\left(\frac{\partial }{\partial N}\right) - \frac{1}{M\sin \sigma }\frac{\partial }{\partial \sigma }, \nonumber \\ \frac{\partial }{\partial \Lambda } &=&\left(\frac{\partial }{\partial \Lambda }\right) -\frac{1}{M\sin I}\frac{\partial }{\partial I}, \end{eqnarray} \tag{ 17 }$$

where the derivatives inside the parentheses refer to the partial derivatives with respect to M,  N, and Λ when these variables explicitly appear in the function that is being derived.

By doing so, we can derive the time evolution of the Euler angles if we know their expression in terms of the Andoyer variables. They turn out to be

$$\begin{eqnarray} \theta &=&I+\sigma \cos \mu +\sigma ^{2}\frac{\cos I}{2\sin I}\sin ^{2}\mu +O(\sigma ^{3}), \nonumber \\ \psi &=&\lambda +\sigma \frac{1}{\sin I}\sin \mu -\sigma ^{2}\frac{\cos I}{ \sin ^{2}I}\sin \mu \cos \mu +O(\sigma ^{3}), \nonumber \\ \phi &=&\mu +\nu -\sigma \frac{\cos I}{\sin I}\sin \mu \nonumber\\ &&+\sigma ^{2}\frac{ (1+\cos ^{2}I)}{2\sin ^{2}I}\sin \mu \cos \mu +O(\sigma ^{3}). \end{eqnarray} \tag{ 18 }$$

A new way to obtain these formulae is explained in the Appendix A. Let us note that up to the first order in σ these expressions agree with those obtained by Kinoshita (1977); however, the second-order contributions in σ are not provided in Kinoshita (1977) and must be considered when constructing a second-order theory. The expansions in the auxiliary variable σ of the different functions of the rotational motion are common in Earth rotation theories, since the value of this angle is about 10⁻⁶ rad (Kinoshita 1977). In addition, these expansions usually allow one to express the explicit dependence of a function on Andoyer variables in a simple form. In a similar way, the evolution of the components of the angular velocity in the terrestrial frame are obtained (Kinoshita 1977) with the help of the relations:

$$\begin{equation} \omega _{x}=\frac{M}{A}\sin \sigma \sin \nu,\ \ \omega _{y}=\frac{M}{A}\sin \sigma \cos \nu,\,\,\,\omega _{z}=\frac{N}{C}. \end{equation} \tag{ 19 }$$

2.4. Integration of the Equations of Motion

The direct integration of the Hamiltonian equations in the form of Equation (14) is completely unfeasible. However, one of the advantages of tackling this problem within the Hamiltonian framework is the possibility to apply the theory of the general perturbations based on Lie series (Hori 1966), since the original Hamiltonian of the system can be put into the form

$$\begin{equation} \mathcal {H(}p,q)=\mathcal {H}_{0}\mathcal {(}p,q)+\mathcal {H}_{1}\mathcal {(} p,q)+\mathcal {H}_{2}\mathcal {(}p,q), \end{equation} \tag{ 20 }$$

$\mathcal {H}_{i}\mathcal {\ }$being of order O(εⁱ), where i is a non-negative integer, and ε is a small parameter measuring the perturbation intensity. Let us sketch this procedure at the second order in ε following a similar method to that of Kinoshita (1977).

The algorithm consists of performing a canonical transformation from the actual canonical set (p, q) to a new one (p*, q*). This transformation is given at the second order by the generating function

$$\begin{equation} \mathcal {W}=\mathcal {W}_{1}+\mathcal {W}_{2}, \end{equation} \tag{ 21 }$$

with $\mathcal {W}_{i}\,{=}\,O(\varepsilon ^{i})$, which depends on the transformed set (p*, q*) of canonical variables. The transformed Hamiltonian at the second order is also of the form

$$\begin{equation} \mathcal {H}^{\ast }(p^{\ast },q^{\ast })=\mathcal {H}_{0}^{\ast }(p^{\ast },q^{\ast })+\mathcal {H}_{1}^{\ast }(p^{\ast },q^{\ast })+\mathcal {H} _{2}^{\ast }(p^{\ast },q^{\ast }), \end{equation} \tag{ 22 }$$

being $\mathcal {H}_{i}^{\ast }(p^{\ast },q^{\ast })=O(\varepsilon ^{i})$. In addition, some extra conditions are imposed on $\mathcal {H}_{i}^{\ast }$ in order to ensure that $\mathcal {H}^{\ast }$ is easier to integrate than $\mathcal {H}$. In particular, we force $\mathcal {H}^{\ast }$ to be free from periodic terms, that is to say, we combine the Lie transformation with an averaging method. By doing so, the transformed Hamiltonian $\mathcal {H} ^{\ast }$ and the generating function $\mathcal {W}$ are determined by the so-called equations of the method (Hori 1966, 1973), which can be written up to the second order as

$$\begin{eqnarray} \mathcal {H}_{0}^{\ast } &=&\mathcal {H}_{0},\ \ \mathcal {H}_{1}^{\ast }= \mathcal {H}_{1\rm {sec}}, \nonumber\\ \mathcal {H}_{2}^{\ast }&=&\mathcal {H}_{2\rm {sec} }+\displaystyle \ \frac{1}{2}\lbrace \mathcal {H}_{1}+\mathcal {H}_{1\rm {sec} };\mathcal {W}_{1}\rbrace _{\rm {sec}}, \nonumber \\ && \nonumber \\ \mathcal {W}_{1} &=&\displaystyle \ \int _{UP}\mathcal {H}_{1\rm {per}}\ dt, \nonumber\\ \mathcal {W}_{2}&=&\displaystyle \ \int _{UP}\mathcal {H}_{2\rm {per}}\ dt+ \displaystyle \ \frac{1}{2}\lbrace \mathcal {H}_{1}+\mathcal {H}_{1\rm {sec}}; \mathcal {W}_{1}\rbrace _{\rm {per}}, \end{eqnarray} \tag{ 23 }$$

where the Poisson brackets are computed in the (p*, q*) canonical set, and the subscripts per and sec denote the periodic or secular part of the corresponding function. Following Hori (1966), these ones are defined through the expressions

$$\begin{equation} \mathcal {G}_{\rm {sec}}=\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0\ UP}^{T}\mathcal {G} dt,\quad \mathcal {G}_{\rm {per}}=\mathcal {G}-\mathcal {G}_{\rm {sec}}. \end{equation} \tag{ 24 }$$

The integrals appearing in Equations (23) and (24) are evaluated along the solutions of the unperturbed problem generated by the Hamiltonian $\mathcal {H }_{0}^{\ast }$, which is obtained by literally substituting the variables (p, q) by the variables (p*, q*) in $\mathcal {H}_{0}$. The time evolution of the transformed canonical variables (p*, q*) are determined by solving the Hamiltonian equations

$$\begin{equation} \frac{dp^{\ast }}{dt}=-\frac{\partial \mathcal {H}^{\ast }}{\partial q^{\ast } },\quad \frac{dq^{\ast }}{dt}=\frac{\partial \mathcal {H}^{\ast }}{\partial p^{\ast }}. \end{equation} \tag{ 25 }$$

In this way, the variation of any function of the canonical variables f(p, q) can be known at the second order by the formula (Hori 1973)

$$\begin{equation} f(p,q)=f^{\ast }(p^{\ast },q^{\ast })+\Delta f(p^{\ast },q^{\ast }), \end{equation} \tag{ 26 }$$

where

$$\begin{eqnarray} f^{\ast }(p^{\ast },q^{\ast })&=& f(p^{\ast },q^{\ast }); \nonumber\\ \Delta f &=&\Delta _{1}f+\Delta _{2}f+\Delta _{3}f\nonumber\\ &\rightarrow &\left\lbrace \begin{array}{@{}l} \Delta _{1}f=\left\lbrace f^{\ast };\mathcal {W}_{1}\right\rbrace \\ \Delta _{2}f=\left\lbrace f^{\ast };\mathcal {W}_{2}\right\rbrace \\ \Delta _{3}f=\displaystyle \ \frac{1}{2}\lbrace \lbrace f^{\ast };\mathcal {W} _{1}\rbrace ;\mathcal {W}_{1}\rbrace \end{array} \right.. \end{eqnarray} \tag{ 27 }$$

By doing so, the determination of the transformed Hamiltonian $\mathcal {H} ^{\ast }$ and the generating function $\mathcal {W}$ allows describing the time evolution, at the second order in the small parameter ε, of any aspect of the Earth rotational motion. In this sense, the Hamiltonian model that we are developing provides a general theory of the rotation of the Earth.

Below, we will apply the former procedure to obtain a second-order analytical integration of the rotational motion of the rigid Earth. In particular, we split the Hamiltonian into the form

$$\begin{equation} \mathcal {H=H}_{0}+\mathcal {H}_{1}\rightarrow \left\lbrace \begin{array}{@{}l@{}} \mathcal {H}_{0}=T \nonumber \\ \mathcal {H}_{1}=E+V \end{array} \right.\!. \end{equation} \tag{ 28 }$$

This decomposition is justified since we have the following orders of magnitude (Getino & Ferrándiz 1995; Kinoshita 1977) for the perturbing terms E and V

$$\begin{equation} \frac{E}{T}\sim 10^{-7}, \frac{V_{\rm{Moon}}}{T}\sim 6\times 10^{-8}, \frac{V_{\rm{Sun}}}{T}\sim 3\times 10^{-8}. \end{equation} \tag{ 29 }$$

It is important to point out that there is one remarkable difference with respect to the procedure followed by Kinoshita (1977): in our approach, both the terms E and V are considered first-order perturbations. This choice simplifies a precise computation of the generating function $\mathcal {W}$, since in the unperturbed problem all the momenta are constant of motion, our separation being valid from the perspective of constructing an ephemeris second-order analytical theory.

This way of implementing the Hori method allows us to identify the different kinds of contributions to any function f. In particular, with the help of Equation (26), we can distinguish a part in the time evolution of f that depends on the arguments Θ_i, that is to say, a short-period part which comes from Δf. On the contrary, the term f(p*, q*) is free from those arguments and can be computed through the transformed Hamiltonian $\mathcal {H}^{\ast }$

$$\begin{equation} \frac{df^{\ast }}{dt}=\lbrace f^{\ast };\mathcal {H}^{\ast }\rbrace =\lbrace f^{\ast };\mathcal {H}_{0}^{\ast }\rbrace +\lbrace f^{\ast };\mathcal {H} _{1}^{\ast }+\mathcal {H}_{2}^{\ast }\rbrace. \end{equation} \tag{ 30 }$$

So, it is the sum of the unperturbed part given by the Hamiltonian $\mathcal { H}_{0}^{\ast }$, plus first- and second-order contributions caused by the perturbation. The secular, or long period, part of f is included in these terms.

Next, we compute the zero-order terms, with the solution of the unperturbed problem, and the first- and second-order terms with the aim of implementing the perturbation method. From Equation (28), it is followed that the second-order terms are due to the interaction of the ecliptic term E with the gravitational potential energy (J₂ term), and to the interaction of the gravitational potential energy with itself; this is the spin–spin coupling. In this way, the second-order contributions are proportional to k'_pk'_q or k'_pe_i, where the parameters k'_p and e_i are defined in Equations (5) and (11).

3.1. Zeroth-order Terms

In order to perform the analytical perturbation method, it is necessary to obtain the solutions of the unperturbed Hamiltonian

$$\begin{equation} \mathcal {H}_{0}^{\ast }=\mathcal {H}_{0}(q^{\ast },p^{\ast })=\frac{1}{2A} (M^{\ast 2}-N^{\ast 2})+\frac{1}{2C}N^{\ast 2}, \end{equation} \tag{ 31 }$$

which are given by the canonical equations

$$\begin{equation} \frac{dq^{\ast }}{dt}=\frac{\partial \mathcal {H}_{0}^{\ast }}{\partial p^{\ast }},\,\quad \frac{dp^{\ast }}{dt}=-\frac{\partial \mathcal {H}_{0}^{\ast } }{\partial q^{\ast }}, \end{equation} \tag{ 32 }$$

obtaining

$$\begin{eqnarray} \Lambda ^{\ast } &=&\Lambda _{0}^{\ast },\quad M^{\ast }=M_{0}^{\ast },\quad N^{\ast }=N_{0}^{\ast }, \nonumber \\ \lambda ^{\ast } &=&\lambda _{0}^{\ast },\quad \mu ^{\ast }=n_{\mu }^{\ast }t+\mu _{0}^{\ast },\quad \nu ^{\ast }=n_{\nu }^{\ast }t+\nu _{0}^{\ast }, \end{eqnarray} \tag{ 33 }$$

where the subscript 0 denotes the constant values, and n*_μ and n*_ν are the corresponding mean motions given by

$$\begin{equation} n_{\mu }^{\ast }=\frac{1}{A}M^{\ast }=\frac{1}{A}M_{0}^{\ast },\ \ n_{\nu }^{\ast }=\left(\frac{1}{C}-\frac{1}{A}\right) N^{\ast }=\left(\frac{1}{C}- \frac{1}{A}\right) N_{0}^{\ast }. \end{equation} \tag{ 34 }$$

As for the auxiliary variables, and considering their definition of Equation (2), we get

$$\begin{equation} \sigma ^{\ast }=\sigma _{0}^{\ast },\,\quad I^{\ast }=I_{0}^{\ast }. \end{equation} \tag{ 35 }$$

3.2. First-order Terms

According to Equations (23) and (24), we should take $\mathcal {H} _{1}^{\ast }$ as the secular or long-period part of $\mathcal {H}_{1}$. From Equations (3), (5), (11), and (33), we will have

$$\begin{equation} \mathcal {H}_{1}^{\ast }=\mathcal {H}_{1\rm{sec}}=E+V_{\rm{sec}}. \end{equation} \tag{ 36 }$$

The secular part of V, $V_{\rm {sec}}$, comes from a term corresponding to Θ_i = 0

$$\begin{equation} V_{\rm{sec}}={\sum _{p=S,M} }k_{p}^{\prime }\frac{1}{2}(3\cos ^{2}\sigma ^{\ast }-1)B_{0}^{\ast }, \end{equation} \tag{ 37 }$$

with B*₀ referring to the value obtained for (m_1i, m_2i, m_3i, m_4i, m_5i) = (0, 0, 0, 0, 0), a term that we will denote as i = 0. Then, $\mathcal {H}_{1\rm {per}}$ will be given by

$$\begin{eqnarray} \mathcal {H}_{1 \rm{per}} &=&\mathcal {H}_{1}-\mathcal {H}_{1 \rm{sec}}= {\sum _{p=S,M} }k_{p}^{\prime }\Bigg[ \frac{1}{2}(3\cos ^{2}\sigma ^{\ast }-1) \nonumber \\ &\times & \sum _{i\ne 0}\! B_{i}^{\ast }\cos \Theta _{i}^{\ast }\,{-}\,\frac{1}{2}\sin 2\sigma ^{\ast }\! \sum _{i}\! \sum _{\tau =\pm 1}C_{i,\tau }^{\ast }\cos (\mu ^{\ast }\,{-}\,\tau \Theta _{i}^{\ast }) \nonumber \\ & + &\frac{1}{4}\sin ^{2}\sigma ^{\ast }\sum _{i}\sum _{\tau =\pm 1}D_{i,\tau }^{\ast }\cos (2\mu ^{\ast }-\tau \Theta _{i}^{\ast })\Bigg ], \end{eqnarray} \tag{ 38 }$$

and with the help of Equations (8), (23), and (33), we can compute the generating function⁴ $\mathcal {W}_{1}$:

$$\begin{eqnarray} \mathcal {W}_{1} &=&{\sum _{p=S,M} }k_{p}^{\prime }\left[ \frac{1}{2} (3\cos ^{2}\sigma ^{\ast }-1)\sum _{i\ne 0}\frac{B_{i}^{\ast }}{n_{i} }\sin \Theta _{i}^{\ast }\right. \nonumber \\ &&-\frac{1}{2}\sin 2\sigma ^{\ast }\sum _{i}\sum _{\tau =\pm 1}\frac{C_{i,\tau }^{\ast }}{n_{\mu }^{\ast }-\tau n_{i}}\sin (\mu ^{\ast }-\tau \Theta _{i}^{\ast }) \nonumber \\ &&\left. +\,\frac{1}{4}\sin ^{2}\sigma ^{\ast }\sum _{i}\sum _{\tau =\pm 1}\frac{ D_{i,\tau }^{\ast }}{2n_{\mu }^{\ast }-\tau n_{i}}\sin (2\mu ^{\ast }-\tau \Theta _{i}^{\ast })\!\right]\!.\quad \ \ \ \ \end{eqnarray} \tag{ 39 }$$

It is important to note that, since in the unperturbed problem λ* = λ*₀, the frequency n_i does not depend on any canonical momenta of the transformed variables.⁵ It represents a difference with respect to the investigation by Kinoshita (1977). For the sake of simplicity, from now on we will omit the asterisk used in the expressions of functions that depend on the transformed canonical variables.

3.3. Second-order Terms

Considering that $\mathcal {H}_{1}=\mathcal {H}_{1\rm {sec}}+\mathcal {H}_{1 \rm {per}}$, and since $\left\lbrace \mathcal {H}_{1\rm {sec}};\mathcal {W} _{1}\right\rbrace _{\rm {sec}}=0$, we can write

$$\begin{eqnarray} \mathcal {H}_{2}^{\ast } &=&\frac{1}{2}\lbrace \mathcal {H}_{1\rm{per}}; \mathcal {W}_{1}\rbrace _{\rm {sec}}, \nonumber \\ \mathcal {W}_{2} &=&\mathcal {W}_{2s}+\mathcal {W}_{2p}\rightarrow \left\lbrace \begin{array}{@{}l@{}} \mathcal {W}_{2s}=\displaystyle \ \int _{UP}\lbrace \mathcal {H}_{1\rm{sec}}; \mathcal {W}_{1}\rbrace dt \\ \ms \mathcal {W}_{2p}=\displaystyle \ \frac{1}{2}\int _{UP}\lbrace \mathcal {H}_{1 \rm {per}};\mathcal {W}_{1}\rbrace _{\rm {per}} dt. \end{array} \right. \end{eqnarray} \tag{ 40 }$$

Therefore, it is necessary to compute the following Poisson brackets

$$\begin{equation} \mathcal {C}_{S}=\left\lbrace \mathcal {H}_{1\rm {sec}};\mathcal {W}_{1}\right\rbrace \,\quad \mathcal {C}_{P}=\frac{1}{2}\lbrace \mathcal {H}_{1\rm {per}};\mathcal {W} _{1}\rbrace. \end{equation} \tag{ 41 }$$

To facilitate these computations, we will make an expansion in the angle σ keeping all the terms in σ⁰, σ¹, and σ². After some algebra, we get

$$\begin{eqnarray} \mathcal {C}_{S} &=&{\sum _{p=S,M}}k_{p}^{\prime }\Bigg[ \left(1-\frac{3}{2}\sigma ^{2}\right)\displaystyle \sum _{i\ne 0}\left(\mathcal {C}_{s0}^{\rm in}\cos \Theta _{i}+ \mathcal {C}_{s0}^{\rm out}\sin \Theta _{i}\right) \nonumber \\ &&+\sigma \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\big(\mathcal {C} _{s1}^{\rm in}\cos (\mu -\tau \Theta _{i})+\mathcal {C}_{s1}^{\rm out}\sin (\mu -\tau \Theta _{i})\big)\nonumber\\ && +\sigma ^{2}\displaystyle \sum _{i} \displaystyle \sum _{\tau =\pm 1}\left(\mathcal {C}_{s2}^{\rm in}\cos (2\mu -\tau \Theta _{i})+\mathcal {C} _{s2}^{\rm out}\sin (2\mu -\tau \Theta _{i})\right)\!\! \Bigg] \nonumber \\ &&+\displaystyle \sum _{p,q=S,M}\frac{k_{p}^{\prime }k_{q}^{\prime }}{M\sin I} \Bigg[ (1-3\sigma ^{2})\displaystyle \sum _{i\ne 0}\mathcal {C}_{s0}^{a}\cos \Theta _{i} \nonumber \\ &&+\sigma \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\mathcal {C} _{s1}^{a}\cos (\mu -\tau \Theta _{i}) \nonumber \\ &&+ \sigma ^{2}\displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1} \mathcal {C}_{s2}^{a}\cos (2\mu -\tau \Theta _{i})\Bigg] +O(\sigma ^{3}) \end{eqnarray} \tag{ 42 }$$

and

$$\begin{eqnarray} \mathcal {C}_{P} &=&\displaystyle \sum _{p,q=S,M}\frac{k_{p}^{\prime }k_{q}^{\prime }}{M\sin I}\displaystyle \sum _{\tau,\rho =\pm 1}\Bigg\lbrace (1-3\sigma ^{2})\nonumber\\ &&\times\displaystyle \sum _{i\ne 0}\displaystyle \sum _{j\ne 0}\mathcal {C} _{p0}^{a}\cos (\tau \Theta _{i}-\rho \Theta _{j}) \nonumber \\ &&+\displaystyle \sum _{i}\displaystyle \sum _{j\ne 0}\left[ \sigma \mathcal {C} _{p1}^{b}\cos (\mu -\tau \Theta _{i}+\rho \Theta _{j})\right.\nonumber\\ &&\left. +\, \sigma _{\ }^{2} \mathcal {C}_{p2}^{b}\cos (2\mu -\tau \Theta _{i}+\rho \Theta _{j})\right] \nonumber \\ &&+\displaystyle \sum _{i}\displaystyle \sum _{j}\left[ \left(\mathcal {C} _{p0}^{c}+\sigma ^{2}\mathcal {C}_{p2}^{c1}\right) \cos (\tau \Theta _{i}-\rho \Theta _{j})\right.\nonumber\\ &&+\sigma \mathcal {C}_{p1}^{c}\cos (\mu -\tau \Theta _{i}+\rho \Theta _{j}) \nonumber \\ && \left. +\,\sigma ^{2}\mathcal {C}_{p2}^{c2}\cos (2\mu -\tau \Theta _{i}-\rho \Theta _{j})\right] \Bigg\rbrace +O(\sigma ^{3}). \end{eqnarray} \tag{ 43 }$$

The expression of the different functions C^β_α in Equations (42) and (43) is explicitly given in Appendix B.

According to Equations (41) and (43), the second-order transformed Hamiltonian $\mathcal {H}_{2}^{\ast }$ is the secular part of $\mathcal {C}_{P}$, therefore,

$$\begin{eqnarray} \mathcal {H}_{2}^{\ast }&=&\displaystyle \sum _{p,q=S,M}\frac{k_{p}^{\prime }k_{q}^{\prime }}{M\sin I}\Bigg\lbrace \displaystyle \sum _{\tau,\rho =\pm 1}\Bigg[ (1-3\sigma ^{2})\displaystyle \sum _{i\ne 0}\displaystyle \sum _{j\ne 0}\mathcal {C }_{p0}^{a}\nonumber\\ &&+\displaystyle \sum _{i}\displaystyle \sum _{j}\left(\mathcal {C} _{p0}^{c}+\sigma ^{2}\mathcal {C}_{p2}^{c1}\right) \Bigg]\Bigg\rbrace _{\tau \Theta _{i}=\rho \Theta _{j}}. \end{eqnarray} \tag{ 44 }$$

Taking into account the Equations (41), (42), and (43), the second-order generating functions $\mathcal {W}_{2s}$ and $\mathcal {W} _{2p}$ are computed by integrating the periodic part of $\mathcal {C}_{s}$ and $\mathcal {C}_{p}$ over the unperturbed problem. In this way, we obtain

$$\begin{eqnarray} &&\mathcal {W}_{2s} ={\sum _{p=S,M}}k_{p}^{\prime }\Bigg[ \left(1-\frac{3}{2}\sigma ^{2}\right)\displaystyle \sum _{i\ne 0}\frac{\left(\mathcal {C}_{s0}^{\rm in}\sin \Theta _{i}-\mathcal {C}_{s0}^{\rm out}\cos \Theta _{i}\right) }{n_{i}} \nonumber \\ &&\quad+\sigma \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\frac{\left(\mathcal {C}_{s1}^{\rm in}\sin (\mu -\tau \Theta _{i})-\mathcal {C}_{s1}^{\rm out}\cos (\mu -\tau \Theta _{i})\right) }{n_{\mu }-\tau n_{i}}\nonumber \\ &&\quad+ \sigma ^{2}\displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\frac{ \left(\mathcal {C}_{s2}^{\rm in}\sin (2\mu -\tau \Theta _{i})-\mathcal {C} _{s2}^{\rm out}\cos (2\mu -\tau \Theta _{i})\right) }{2n_{\mu }-\tau n_{i}} \Bigg] \nonumber \\ &&\quad+\displaystyle \sum _{p,q=S,M}\frac{k_{p}^{\prime }k_{q}^{\prime }}{M\sin I} \Bigg[ (1-3\sigma ^{2})\displaystyle \sum _{i\ne 0}\frac{\mathcal {C} _{s0}^{a}\sin \Theta _{i}}{n_{i}} \nonumber \\ &&\quad+\sigma \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\frac{\mathcal {C} _{s1}^{a}\sin (\mu -\tau \Theta _{i})}{n_{\mu }-\tau n_{i}} \nonumber \\ &&\quad+ \sigma ^{2}\displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\frac{ \mathcal {C}_{s2}^{a}\sin (2\mu -\tau \Theta _{i})}{2n_{\mu }-\tau n_{i}} \Bigg] +O(\sigma ^{3}) \end{eqnarray} \tag{ 45 }$$

and

$$\begin{eqnarray} \mathcal {W}_{2p} &=&\displaystyle \sum _{p,q=S,M}\frac{k_{p}^{\prime }k_{q}^{\prime }}{M\sin I}\displaystyle \sum _{\tau,\rho =\pm 1}\Bigg\lbrace (1-3\sigma ^{2})\nonumber\\ &&\times\displaystyle \sum _{i\ne 0}\displaystyle \sum _{j\ne 0}\displaystyle \underbrace{\frac{\mathcal {C} _{p0}^{a}\sin (\tau \Theta _{i}-\rho \Theta _{j})}{\tau n_{i}-\rho n_{j}}}_{ \tau \Theta _{i}\ne \rho \Theta _{j}} \nonumber \\ &&+\displaystyle \sum _{i}\displaystyle \sum _{j\ne 0}\Bigg[ \sigma \frac{\mathcal {C }_{p1}^{b}\sin (\mu -\tau \Theta _{i}+\rho \Theta _{j})}{n_{\mu }-\tau n_{i}+\rho n_{j}}\nonumber\\ &&+\sigma ^{2}\frac{\mathcal {C}_{p2}^{b}\sin (2\mu -\tau \Theta _{i}+\rho \Theta _{j})}{2n_{\mu }-\tau n_{i}+\rho n_{j}}\Bigg]\nonumber \\ &&+\displaystyle \sum _{i}\displaystyle \sum _{j}\Bigg[ \displaystyle \underbrace{\frac{\big(\mathcal {C} _{p0}^{c}+\sigma ^{2}\mathcal {C}_{p2}^{c1}\big) \sin (\tau \Theta _{i}-\rho \Theta _{j})}{\tau n_{i}-\rho n_{j}}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}\nonumber\\ &&+\sigma \frac{\mathcal {C} _{p1}^{c}\sin (\mu -\tau \Theta _{i}+\rho \Theta _{j})}{n_{\mu }-\tau n_{i}+\rho n_{j}} \nonumber \\ && +\sigma ^{2}\frac{\mathcal {C} _{p2}^{c2}\sin (2\mu -\tau \Theta _{i}-\rho \Theta _{j})}{2n_{\mu }-\tau n_{i}-\rho n_{j}}\Bigg] \Bigg\rbrace +O(\sigma ^{3}). \end{eqnarray} \tag{ 46 }$$

Let us underline that, as we have explicitly indicated, in the expression $\mathcal {W}_{2p}$, the sums for some terms exclude the combinations that come from secular contributions.

4.1. Motion of the Equatorial Plane

As we have pointed out in Section 2, the Euler angles ψ (longitude) and θ (obliquity) determine the longitude of the node and the inclination of the equatorial plane. These angles are related to the Andoyer variables by means of Equation (18). Following Kinoshita (1977), the study of the evolution of those angles is usually broken down into two parts:

$$\begin{equation} \psi =\lambda +(\psi -\lambda),\ \quad \theta =I+(\theta -I), \end{equation} \tag{ 47 }$$

with

$$\begin{eqnarray} \psi -\lambda &=&\sigma \frac{1}{\sin I}\sin \mu -\sigma ^{2}\frac{\cos I}{ \sin ^{2}I}\sin \mu \cos \mu +O(\sigma ^{3}), \nonumber \\ && \nonumber \\ \theta -I &=&\sigma \cos \mu +\sigma ^{2}\frac{\cos I}{2\sin I}\sin ^{2}\mu +O(\sigma ^{3}).\ \end{eqnarray} \tag{ 48 }$$

In this way, the motion of the equatorial plane is described through the motion of a plane perpendicular to the angular momentum vector, referred to as the Andoyer plane, determined by the angles λ and I, plus some additional terms depending on σ, which represent the motion of the equatorial plane with respect to the Andoyer plane. At the first order (Kinoshita 1977), the short-period terms of the Andoyer plane are known as Poisson terms since they satisfy the Poisson equations in first-order periodic perturbation. The short-period part in the motion of the equatorial plane with respect to the Andoyer plane is known as the Oppolzer terms. For the sake of simplicity, we will extend this nomenclature to the second-order case.

First, let us study the motion of the Andoyer plane. According to Equation (26), the evolution of these terms is given by

$$\begin{equation} \lambda =\lambda ^{\ast }+\Delta \lambda, \,\quad I=I^{\ast }+\Delta I \end{equation} \tag{ 49 }$$

The solution of λ* and I* is obtained from Equation (30)

$$\begin{equation} \frac{d\lambda ^{\ast }}{dt}=\lbrace \lambda ^{\ast };\mathcal {H}^{\ast }\rbrace, \,\quad \frac{dI^{\ast }}{dt}=\lbrace I^{\ast };\mathcal {H}^{\ast }\rbrace, \end{equation} \tag{ 50 }$$

providing the secular evolution of these two angles, since the solutions Δλ and ΔI are short periodic. Taking into account that the transformed Hamiltonian $\mathcal {H}^{\ast }=\mathcal {H}_{0}^{\ast }+ \mathcal {H}_{1}^{\ast }+\mathcal {H}_{2}^{\ast }$ is given by Equations (31), (36), (37), and (44), and performing the corresponding computations, we obtain

$$\begin{eqnarray} \frac{d\lambda ^{\ast }}{dt} &=&S_{E}^{L}+{\sum _{p=S,M} }\frac{k_{p} }{\sin I^{\ast }}S_{1}^{L}+\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin ^{2}I^{\ast }}\nonumber\\ &&\times \left\lbrace \displaystyle \sum _{\tau,\rho =\pm 1}\left[ \displaystyle \sum _{i,j\ne 0} S_{2a}^{L}+\displaystyle \sum _{i,j} S_{2b}^{L}\right]\right\rbrace _{\tau \Theta _{i}=\rho \Theta _{j}}+O(\sigma ^{\ast }), \nonumber \\ \frac{dI^{\ast }}{dt} &=&S_{E}^{O}, \end{eqnarray} \tag{ 51 }$$

with k_p, q = k'_p, q/M. In these formulae, the superscripts L and O stand, respectively, for longitude and obliquity, while the subscripts E, 1 and 2a,  2b  indicate the origin of each term: S^L,O_E are due to the effect of the ecliptic (Equation (11)), S^L₁ comes from $\mathcal {W}_{1}$ (first-order terms), and S^L_2a,  S^L_2b come from $\mathcal {W}_{2}$ (second-order terms). Explicit expressions for each function S^β_α are listed in Appendix C.

On the other hand, the terms Δλ and ΔI give the Poisson terms, or the nutations of the angles λ and I, which taking into account Equation (27) are

$$\begin{eqnarray} \Delta \lambda &=& \lbrace \lambda ^{\ast };\mathcal {W}_{1}\rbrace + \lbrace \lambda ^{\ast };\mathcal {W}_{2}\rbrace +\frac{1}{2}\lbrace \lbrace \lambda ^{\ast };\mathcal {W}_{1}\rbrace ;\mathcal {W}_{1}\rbrace, \nonumber \\ && \nonumber \\ \Delta I &=&\lbrace I^{\ast };\mathcal {W}_{1}\rbrace +\lbrace I^{\ast }; \mathcal {W}_{2}\rbrace +\frac{1}{2}\lbrace \lbrace I^{\ast };\mathcal {W} _{1}\rbrace ;\mathcal {W}_{1}\rbrace. \end{eqnarray} \tag{ 52 }$$

With the help of Equations (39), (41), (45), and (46), the corresponding nutations can be finally arranged as

$$\begin{eqnarray} \Delta \lambda &=&{\sum _{p=S,M} }\frac{k_{p}}{\sin I}\displaystyle \sum _{i\ne 0}\big[ \big(\mathcal {L}_{1}+\mathcal {L}_{E}^{\rm in}\big)\sin \Theta _{i}+ \mathcal {L}_{E}^{\rm out}\cos \Theta _{i}\big] \nonumber \\ &&+\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin ^{2}I}\Bigg[ \displaystyle \sum _{ i\ne 0}\mathcal {L}_{2}^{s}\sin \Theta _{i}\nonumber \\ &&+\displaystyle \sum _{\tau,\rho =\pm 1}\Bigg(\displaystyle \sum _{i,j\ne 0}\big(\displaystyle \underbrace{\mathcal { L}_{2}^{p1}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+\mathcal {L}_{3}^{1}\big)\nonumber\\ &&+\displaystyle \sum _{i,j}\big(\displaystyle \underbrace{\mathcal {L}_{2}^{p2}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+ \mathcal {L}_{3}^{2}\big)\Bigg) \sin (\tau \Theta _{i}-\rho \Theta _{j})\Bigg] +O(\sigma). \nonumber \\ \end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray} \Delta I &=&{\sum _{p=S,M} }\frac{k_{p}}{\sin I}\displaystyle \sum _{i\ne 0}\big[ \big(\mathcal {O}_{1}+\mathcal {O}_{E}^{\rm in}\big)\cos \Theta _{i}+\mathcal { O}_{E}^{\rm out}\sin \Theta _{i}\big] \nonumber \\ &&+\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin ^{2}I}\Bigg[ \displaystyle \sum _{ i\ne 0}\mathcal {O}_{2}^{s}\cos \Theta _{i} \nonumber \\ &&+ \displaystyle \sum _{\tau,\rho =\pm 1}\Bigg(\displaystyle \sum _{i,j\ne 0}\big(\displaystyle \underbrace{\mathcal { O}_{2}^{p1}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+\mathcal {O}_{3}^{1}\big)\nonumber\\ &&+\displaystyle \sum _{i,j}\big(\displaystyle \underbrace{\mathcal {O}_{2}^{p2}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+ \mathcal {O}_{3}^{2}\big)\Bigg) \cos (\tau \Theta _{i}-\rho \Theta _{j})\Bigg] +O(\sigma). \nonumber \\ \end{eqnarray} \tag{ 54 }$$

From these expressions, we can observe the existence of out-of-phase terms $\mathcal {L}_{E}^{\rm out}$ and $\mathcal {O}_{E}^{\rm out}$ that are originated by the moving ecliptic. In addition, the term E also gives raise to in-phase terms $\mathcal {L}_{E}^{\rm in}$ and $\mathcal {O}_{E}^{\rm in}$. In the formulae (53) and (54), the first-order effects are given by $\mathcal {L}_{1}$ and $\mathcal {O}_{1}$ and the subscript 2 stands for second terms arising from $\mathcal {W}_{2}$, where $\mathcal { L}_{2}^{s}$ and $\mathcal {O}_{2}^{s}$ come from $\mathcal {W}_{2s}$, and $\mathcal { L}_{2}^{p1,2}$ and $\mathcal {O}_{2}^{p1,2}$ correspond to $\mathcal {W}_{2p}$. Finally, terms in $\mathcal {L}_{3}^{1,2}$ and $\mathcal {O}_{3}^{1,2}$ are due to Δ₃f in Equation (27). These coefficients are listed in Appendix D.

In a similar way, the motion of the equatorial plane with respect to the Andoyer plane is computed with the aid of

$$\begin{equation} \psi -\lambda =(\psi ^{\ast }-\lambda ^{\ast }) +\Delta (\psi -\lambda), \theta -I=(\theta ^{\ast }-I^{\ast }) +\Delta \left(\theta -I\right). \end{equation} \tag{ 55 }$$

Proceeding as in the case of the Andoyer plane, we find that

$$\begin{equation} \frac{d(\psi ^{\ast }-\lambda ^{\ast })}{dt}=0+O(\sigma ^{\ast })\,\ \frac{d(\theta ^{\ast }-I^{\ast })}{dt}=0+O(\sigma ^{\ast }).\ \end{equation} \tag{ 56 }$$

Therefore, the secular evolution of the equatorial plane and Andoyer plane is the same at the zero order in σ*. Following the nomenclature introduced in the Section 2, this part of the motion of the angles ψ and θ will be the precession, providing Equation (51) the rates of precession in longitude and in obliquity.

With respect to the Oppolzer terms of Δ(ψ − λ) and Δ(θ − I), we find that

$$\begin{eqnarray} \Delta (\psi -\lambda) &=&{\sum _{p=S,M} }\frac{k_{p}}{\sin I}\nonumber\\ &&\times \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\left[ \big(\mathcal {N}_{1}+ \mathcal {N}_{E}^{\rm in}\big)\tau \sin \Theta _{i}+\mathcal {N}_{E}^{\rm out}\cos \Theta _{i}\right] \nonumber \\ &&+\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin ^{2}I}\Bigg[ \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\mathcal {N}_{2}^{s}\ \tau \sin \Theta _{i} \nonumber \\ &&+\displaystyle \sum _{i}\displaystyle \sum _{j\ne 0}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\displaystyle \underbrace{\mathcal { N}_{2}^{p1}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+\mathcal {N}_{3}^{1L})\sin (\tau \Theta _{i}-\rho \Theta _{j}\big) \nonumber \\ && +\displaystyle \sum _{i}\displaystyle \sum _{j}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\mathcal {N}_{2}^{p2}+\mathcal {N}_{3}^{2}+\mathcal {N}_{3}^{3L}\big)\nonumber\\ &&\times\sin (\tau \Theta _{i}-\rho \Theta _{j})\Bigg]+ O(\sigma). \end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray} \Delta (\theta -I) &=&{\sum _{p=S,M} }k_{p}\displaystyle \sum _{i} \displaystyle \sum _{\tau =\pm 1}\big[ \big(\mathcal {N}_{1}+\mathcal {N} _{E}^{\rm in}\big)\cos \Theta _{i}\nonumber\\ &&-\mathcal {N}_{E}^{\rm out}\ \tau \sin \Theta _{i} \big] \nonumber \\ &&+\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin I}\Bigg[ \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\mathcal {N}_{2}^{s}\ \cos \Theta _{i} \nonumber \\ &&+\displaystyle \sum _{i}\displaystyle \sum _{j\ne 0}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\displaystyle \underbrace{\mathcal { N}_{2}^{p1}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+\mathcal {N}_{3}^{1O}\big)\cos (\tau \Theta _{i}-\rho \Theta _{j}) \nonumber \\ && +\displaystyle \sum _{i}\displaystyle \sum _{j}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\mathcal {N}_{2}^{p2}+\mathcal {N}_{3}^{2}+\mathcal {N}_{3}^{3O}\big)\nonumber\\ &&\times\cos (\tau \Theta _{i}-\rho \Theta _{j})\Bigg] +O(\sigma).\vspace*{-20pt} \end{eqnarray} \tag{ 58 }$$

The meaning of functions $\mathcal {N}_{\alpha }^{\beta }$ is similar to that of the Andoyer plane case. It is worth to point out that the amplitudes $\mathcal {N}_{3}^{3L}$ and $\mathcal {N}_{3}^{3O}$ are exclusively due to the second order in σ in the expansion of the Euler angles in terms of the Andoyer variables (Equation (18)). Expressions of these functions can be found in Appendix E. The sum of Equations (53) and (57), and (54) and (58), provides the short-period part of the motion of the angles ψ and θ, that is to say, the nutation of these angles.

4.2. Time Evolution of $\protect \phi$

The angle ϕ can also be written into the form (Equation (18))

$$\begin{eqnarray} \phi &=&\mu +\nu +\left(\phi -\mu -\nu \right), \nonumber \\ \phi -\mu -\nu &=&{-}\sigma \frac{\cos I}{\sin I}\sin \mu\nonumber\\ &&+\sigma ^{2}\frac{ 1+\cos ^{2}I}{2\sin ^{2}I}\sin \mu \cos \mu +O(\sigma ^{3}). \end{eqnarray} \tag{ 59 }$$

Following a similar procedure as in the equatorial plane situation, we firstly get the solution of the part independent of σ

$$\begin{eqnarray} &&\frac{d\left(\mu ^{\ast }+\nu ^{\ast }\right) }{dt} =\frac{M^{\ast }-N^{\ast }}{A}+\frac{N^{\ast }}{C}+S_{E}^{M}+{\sum _{p=S,M} }\frac{ k_{p}}{\sin I^{\ast }}S_{1}^{M}\nonumber \\ &&\quad\! +\!\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin ^{2}I^{\ast }}\displaystyle \sum _{\tau,\rho =\pm 1} \left[ \displaystyle \sum _{i,j\ne 0}S_{2a}^{M}+\displaystyle \sum _{ i,j} S_{2b}^{M}\!\right] _{\tau \Theta _{i}=\rho \Theta _{j}}\!\!\!\! +O(\sigma ^{\ast }),\nonumber\\ \end{eqnarray} \tag{ 60 }$$

and

$$\begin{eqnarray} \Delta (\mu +\nu) &= &{\sum _{p=S,M} }\frac{k_{p}}{\sin I}\displaystyle \sum _{ i\ne 0}\left[ \big(\mathcal {M}_{1}+\mathcal {M}_{E}^{\rm in}\big)\sin \Theta _{i}+ \mathcal {M}_{E}^{\rm out}\cos \Theta _{i}\right] \nonumber \\ &&+\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin ^{2}I}\Bigg[ \displaystyle \sum _{ i\ne 0}\mathcal {M}_{2}^{s}\sin \Theta _{i} \nonumber \\ &&+\displaystyle \sum _{\tau,\rho =\pm 1}\Bigg(\displaystyle \sum _{i,j\ne 0}\big(\displaystyle \underbrace{\mathcal { M}_{2}^{p1}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+\mathcal {M}_{3}^{1}\big)\nonumber\\ &&+\displaystyle \sum _{i,j}\big(\displaystyle \underbrace{\mathcal {M}_{2}^{p2}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+ \mathcal {M}_{3}^{2}\big)\Bigg) \sin (\tau \Theta _{i}-\rho \Theta _{j})\Bigg] +O(\sigma).\nonumber\\ \end{eqnarray} \tag{ 61 }$$

As for the remaining terms in Equation (59), we have

$$\begin{equation} \frac{d(\phi ^{\ast }-\mu ^{\ast }-\nu ^{\ast })}{dt}=0+O(\sigma ^{\ast }). \end{equation} \tag{ 62 }$$

and

$$\begin{eqnarray} &&\Delta (\phi -\mu -\nu) =-\!\cos I{\sum _{p=S,M} }\frac{k_{p}}{\sin I }\nonumber\\ &&\quad\times\displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\left[ \big(\mathcal {N}_{1}+ \mathcal {N}_{E}^{\rm in}\big)\tau \sin \Theta _{i}+\mathcal {N}_{E}^{\rm out}\cos \Theta _{i}\right] \nonumber \\ &&\quad-\cos I\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin ^{2}I}\Bigg[ \displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\mathcal {N}_{2}^{s}\ \tau \sin \Theta _{i} \nonumber \\ &&\quad+\displaystyle \sum _{i}\displaystyle \sum _{j\ne 0}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\displaystyle \underbrace{\mathcal { N}_{2}^{p1}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+\mathcal {N}_{3}^{1M}\big)\sin (\tau \Theta _{i}-\rho \Theta _{j}) \nonumber \\ &&\quad +\displaystyle \sum _{i}\displaystyle \sum _{j}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\displaystyle \underbrace{\mathcal { N}_{2}^{p2}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+\mathcal {N}_{3}^{2}-\mathcal {N}_{3}^{3L}\big)\nonumber\\ &&\quad\times\sin (\tau \Theta _{i}-\rho \Theta _{j})\Bigg] +O(\sigma). \end{eqnarray} \tag{ 63 }$$

The functions $S_{\alpha }^{M},\; \mathcal {M}_{\alpha }^{\beta }$, and $\mathcal {N}_{\alpha }^{\beta }$ have the same meaning as in the case of the motion of the equatorial plane; their expressions are included, respectively, in Appendices C, D, and E.

4.3. Polar Motion

The components ω_x and ω_y of the angular velocity vector in the terrestrial frame determine the polar motion. It is convenient to introduce the dimensionless complex magnitude

$$\begin{equation} U=\frac{1}{\omega _{E}}(\omega _{x}+i\omega _{y}), \end{equation} \tag{ 64 }$$

where ω_E is the mean value of the angular velocity of the Earth. By considering Equation (19), we can determine the expression of the function U in terms of the Andoyer variables

$$\begin{equation} U=i\frac{M}{A\omega _{E}}\sin \sigma e^{-i\nu }. \end{equation} \tag{ 65 }$$

As in the previous sections, the evolution of this function is obtained with the aid of the equation U = U* + ΔU. The solution of the first term is given by

$$\begin{equation} \frac{dU^{\ast }}{dt}=\lbrace U^{\ast },\mathcal {H}^{\ast }\rbrace =i\ \mathcal {F}\ U^{\ast }. \end{equation} \tag{ 66 }$$

Therefore, the part of the solution independent of Θ_i harmonically evolves with a frequency $\mathcal {F}$. This frequency characterizing the polar motion can be cast as the sum of three terms

$$\begin{equation} \mathcal {F=F}_{0}+\mathcal {F}_{1}+\mathcal {F}_{2},\vspace*{-12pt} \end{equation} \tag{ 67 }$$

with

$$\begin{eqnarray} &&\mathcal {F}_{0}=\omega _{E}\displaystyle \ \frac{C-A}{A} \rightarrow \rm {Euler\hbox{-}free\ frequency,} \nonumber \\ &&\mathcal {F}_{1}=-3{\displaystyle \sum _{p=S,M} }k_{p}B_{0} \rightarrow \rm { First\hbox{-}order\ perturbation,} \nonumber \\ &&\mathcal {F}_{2}=\displaystyle \sum _{p,q=S,M}\displaystyle \frac{k_{p}k_{q}}{\sin I^{\ast }}\Bigg\lbrace \displaystyle \sum _{\tau,\rho =\pm 1}\Bigg[ \displaystyle -6\sum _{i,j\ne 0} \mathcal {C}_{p0}^{a}\nonumber\\ &&+2\displaystyle \sum _{i,j}\mathcal {C}_{p2}^{c1}\Bigg] \Bigg\rbrace _{\tau \Theta _{i}=\rho \Theta _{j}}\!\!\!\!\! \rightarrow \rm {Second\hbox{-}order\ perturbation,}\quad \end{eqnarray} \tag{ 68 }$$

where $\mathcal {C}_{p0}^{a}$ and $\mathcal {C}_{p2}^{c1}$ are components of the Poisson bracket $\mathcal {C}_{P}$ given in Appendix B. Let us note that $\mathcal {F}$ is the real frequency of the polar motion in presence of the perturbations. As it can be seen, its value is shifted by an amount $\mathcal {F}_{1}+\mathcal {F}_{2}$ with respect to the Eulerian free frequency $\mathcal {F}_{0}$.

As for the short-period solution, from Equations (27) and (65), and separating real and imaginary parts, it turns out to be

$$\begin{eqnarray} &&\left\lbrace \begin{array}{@{}l@{}} \Delta \omega _{x} \nonumber \\ \Delta \omega _{y} \end{array} \right\rbrace &=&\frac{C}{A}\Bigg[ {\sum _{p=S,M} }k_{p}\displaystyle \sum _{i} \displaystyle \sum _{\tau =\pm 1}\big(\mathcal {N}_{1}+\mathcal {N}_{E}^{\rm in}\big)\nonumber\\ &&\times\left\lbrace \begin{array}{@{}l@{}} \sin (\mu +\nu -\tau \Theta _{i}) \nonumber \\ \cos (\mu +\nu -\tau \Theta _{i}) \end{array} \right\rbrace \nonumber \\ &&+{\sum _{p=S,M} }k_{p}\displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\mathcal {N}_{E}^{\rm out}\nonumber\\ &&\times\left\lbrace \begin{array}{@{}r@{}} -\cos (\mu +\nu -\tau \Theta _{i}) \nonumber \\ \sin (\mu +\nu -\tau \Theta _{i}) \end{array} \right\rbrace \nonumber \\ && +\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin I}\displaystyle \sum _{i}\displaystyle \sum _{\tau =\pm 1}\mathcal {N}_{2}^{s} \nonumber \\ &&\times\left\lbrace \begin{array}{@{}l@{}} \sin (\mu +\nu -\tau \Theta _{i}) \nonumber \\ \cos (\mu +\nu -\tau \Theta _{i}) \end{array} \right\rbrace \nonumber \\ && +\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin I}\displaystyle \sum _{i}\displaystyle \sum _{j\ne 0}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\displaystyle \underbrace{\mathcal {N}_{2}^{p1}}_{ \tau \Theta _{i}\ne \rho \Theta _{j}}+ \mathcal {N}_{3}^{1P}\big) \nonumber \\ &&\times\left\lbrace \begin{array}{@{}l@{}} \sin (\mu +\nu -\tau \Theta _{i}+\rho \Theta _{j}) \nonumber \\ \cos (\mu +\nu -\tau \Theta _{i}+\rho \Theta _{j}) \end{array} \right\rbrace \nonumber \\ && +\displaystyle \sum _{p,q=S,M}\frac{k_{p}k_{q}}{\sin I}\displaystyle \sum _{i}\displaystyle \sum _{j}\displaystyle \sum _{\tau,\rho =\pm 1}\big(\displaystyle \underbrace{\mathcal {N}_{2}^{p2}}_{\tau \Theta _{i}\ne \rho \Theta _{j}}+ \mathcal {N}_{3}^{2}\big) \nonumber \\ &&\times \left\lbrace \begin{array}{@{}l@{}} \sin (\mu +\nu -\tau \Theta _{i}+\rho \Theta _{j}) \\ \cos (\mu +\nu -\tau \Theta _{i}+\rho \Theta _{j}) \end{array} \right\rbrace \Bigg] +O(\sigma). \end{eqnarray} \tag{ 69 }$$

The explicit expressions of the functions $\mathcal {N}_{\alpha }^{\beta }$ are listed in Appendix E.

In this investigation, we have derived a complete analytical solution of the rotation of the rigid Earth up to the second order, which accounts for the spin–spin coupling. We have considered the perturbation on the rotational motion caused by the J₂ term of the gravitational potential energy expansion (Equation (41)), as well as an additional term E that emerges when referring the motion to the moving ecliptic. As a first result, first-order contributions in the small parameters k_p and e_i coincide with those determined in Kinoshita (1977). With respect to the second-order contributions, and as far as we know, we have performed for the first time a complete study of the crossed terms coming from the interaction of the J₂ term with the moving ecliptic term E and from the interaction of the J₂ term with itself. These crossed terms originate the different amplitudes appearing in the analytical formulae given in Section 4.

In Table 1, we present those amplitudes indicating their origin and the associated crossed terms causing them. The column "Origin" refers to the second-order part of the Poisson bracket from which the amplitude is derived according to Equation (27) for the short-period parts and to Equations (44) and (50) in the case of the long-period contributions. On the other hand, the column "Crossed terms" indicates the crossed interaction responsible of the amplitude. For further comparisons, we have expanded any function of the Andoyer variables into the form

$$\begin{equation} f=f^{(0)}+f^{(1)}+f^{(2)}+O(\sigma ^{3}), \end{equation} \tag{ 70 }$$

Table 1. Second-order Contributions to the Earth Rotational Motion

Function	Amplitude	Origin	Crossed terms
λ,  μ + ν	$\mathcal {S}_{2a}^{L},\,\mathcal {S}_{2a}^{M}$	$\left\lbrace f; \left\lbrace \mathcal {H}_{1\rm per};\mathcal {W}_{1}\right\rbrace _{\rm sec}\right\rbrace$a	V(0)per–$\mathcal {W}_{1}^{(0)}$
λ,  μ + ν	$\mathcal {S}_{2b}^{L},\,\mathcal {S}_{2b}^{M}$	$\left\lbrace f; \left\lbrace \mathcal {H}_{1\rm per};\mathcal {W}_{1}\right\rbrace _{\rm sec}\right\rbrace$a	V(1)per–$\mathcal {W}_{1}^{(1)}$
U(1)	$\mathcal {F}_{2}$	$\left\lbrace f; \left\lbrace \mathcal {H}_{1\rm per};\mathcal {W}_{1}\right\rbrace _{\rm sec}\right\rbrace$b	V(0,1)per–$\mathcal {W}_{1}^{(0,1)}$
λ,  I,  μ + ν	$\mathcal {L}_{E}^{\rm in},\,\mathcal {O}_{E}^{\rm in},\,\mathcal {M}_{E}^{\rm in}$	$\left\lbrace f;\mathcal {W}_{2s}\right\rbrace$	E(0)–$\mathcal {W}_{1}^{(0)}$
λ,  I,  μ + ν	$\mathcal {L}_{E}^{\rm out},\,\mathcal {O}_{E}^{\rm out},\,\mathcal {M}_{E}^{\rm out}$	$\left\lbrace f;\mathcal {W}_{2s}\right\rbrace$	E(0)–$\mathcal {W}_{1}^{(0)}$
λ,  I,  μ + ν	$\mathcal {L}_{2}^{s},\,\mathcal {O}_{2}^{s},\,\mathcal {M}_{2}^{s}$	$\left\lbrace f;\mathcal {W}_{2s}\right\rbrace$	V(0)sec–$\mathcal {W}_{1}^{(0)}$
λ,  I,  μ + ν	$\mathcal {L}_{2}^{p1},\,\mathcal {O}_{2}^{p1},\,\mathcal {M}_{2}^{p1}$	$\left\lbrace f;\mathcal {W}_{2p}\right\rbrace$	V(0)per–$\mathcal {W}_{1}^{(0)}$
λ,  I,  μ + ν	$\mathcal {L}_{2}^{p2},\,\mathcal {O}_{2}^{p2},\,\mathcal {M}_{2}^{p2}$	$\left\lbrace f;\mathcal {W}_{2p}\right\rbrace$	V(1)per–$\mathcal {W}_{1}^{(1)}$
λ,  I,  μ + ν	$\mathcal {L}_{3}^{1},\,\mathcal {O}_{3}^{1},\,\mathcal {M}_{3}^{1}$	$\left\lbrace \left\lbrace f;\mathcal {W}_{1}\right\rbrace ;\mathcal {W}_{1}\right\rbrace$	$\mathcal {W}_{1}^{(0)}$–$\mathcal {W}_{1}^{(0)}$
λ,  I,  μ + ν	$\mathcal {L}_{3}^{2},\,\mathcal {O}_{3}^{2},\,\mathcal {M}_{3}^{2}$	$\left\lbrace \left\lbrace f;\mathcal {W} _{1}\right\rbrace ;\mathcal {W}_{1}\right\rbrace$	$\mathcal {W}_{1}^{(1)}$–$\mathcal {W}_{1}^{(1)}$
(ψ − λ)(1), (θ − I)(1), (ϕ − μ − ν)(1), U(1)	$\mathcal {N}_{E}^{\rm in},\,\mathcal {N}_{E}^{\rm in},\,\mathcal {N}_{E}^{\rm in},\,\mathcal {N}_{E}^{\rm in}$	$\left\lbrace f;\mathcal {W}_{2s}\right\rbrace$	E(0)–$\mathcal {W}_{1}^{(1)}$
(ψ − λ)(1), (θ − I)(1), (ϕ − μ − ν)(1), U(1)	$\mathcal {N}_{E}^{\rm out},\,\mathcal {N}_{E}^{\rm out},\,\mathcal {N}_{E}^{\rm out},\,\mathcal {N}_{E}^{\rm out}$	$\left\lbrace f;\mathcal {W}_{2s}\right\rbrace$	E(0)–$\mathcal {W}_{1}^{(1)}$
(ψ − λ)(1), (θ − I)(1), (ϕ − μ − ν)(1), U(1)	$\mathcal {N}_{2}^{s},\,\mathcal {N}_{2}^{s},\,\mathcal {N}_{2}^{s},\,\mathcal {N}_{2}^{s}$	$\left\lbrace f;\mathcal {W}_{2s}\right\rbrace$	$V_{\rm sec}^{(0,2_B)}$–$\mathcal {W}_{1}^{(1)}$
(ψ − λ)(1), (θ − I)(1), (ϕ − μ − ν)(1), U(1)	$\mathcal {N}_{2}^{p1},\,\mathcal {N}_{2}^{p1},\,\mathcal {N}_{2}^{p1},\,\mathcal {N}_{2}^{p1}$	$\left\lbrace f;\mathcal {W}_{2p}\right\rbrace$	$V_{\rm per}^{(0,1,2_B)}$–$\mathcal {W}_{1}^{(0,1,2_B)}$
(ψ − λ)(1), (θ − I)(1), (ϕ − μ − ν)(1), U(1)	$\mathcal {N}_{2}^{p2},\,\mathcal {N}_{2}^{p2},\,\mathcal {N}_{2}^{p2},\,\mathcal {N}_{2}^{p2}$	$\left\lbrace f;\mathcal {W}_{2p}\right\rbrace$	$V_{\rm per}^{(1,2_D)}$–$\mathcal {W}_{1}^{(1,2_D)}$
(ψ − λ)(1), (θ − I)(1), (ϕ − μ − ν)(1), U(1)	$\mathcal {N}_{3}^{1L},\,\mathcal {N}_{3}^{1O},\,\mathcal {N}_{3}^{1M},\,\mathcal {N}_{3}^{1P}$	$\left\lbrace \left\lbrace f;\mathcal {W}_{1}\right\rbrace ;\mathcal {W}_{1}\right\rbrace$	$\mathcal {W}_{1}^{(0,1,2_B)}$–$\mathcal {W}_{1}^{(0,1)}$
(ψ − λ)(1), (θ − I)(1), (ϕ − μ − ν)(1), U(1)	$\mathcal {N}_{3}^{2},\,\mathcal {N}_{3}^{2},\,\mathcal {N}_{3}^{2},\,\mathcal {N}_{3}^{2}$	$\left\lbrace \left\lbrace f;\mathcal {W}_{1}\right\rbrace ;\mathcal {W}_{1}\right\rbrace$	$\mathcal {W}_{1}^{(2_D)}$–$\mathcal {W}_{1}^{(1)}$
(ψ − λ)(2), (θ − I)(2), (ϕ − μ − ν)(2)	$\mathcal {N}_{3}^{3L},\,\mathcal {N}_{3}^{3O},\,\mathcal {N}_{3}^{3L}$	$\left\lbrace \left\lbrace f;\mathcal {W}_{1}\right\rbrace ;\mathcal {W}_{1}\right\rbrace$	$\mathcal {W}_{1}^{(1)}$–$\mathcal {W}_{1}^{(1)}$

Notes. ^aThese amplitudes provide the second-order long-period contributions for λ, μ + ν, and also for ψ and ϕ according to Equation (56). The remaining amplitudes in this table, with the exception of $\mathcal {F}_{2}$, provide second-order short-period contributions. See the main text for a discussion. ^bThe amplitude $\mathcal {F}_{2}$ provides the second-order contribution to the Eulerian free frequency $\mathcal {F}_{0}$ of the polar motion (Equation (68)).

Download table as:  ASCIITypeset image

where f⁽ⁱ⁾ is proportional to σⁱ. For example, from Equation (18) the function (ψ − λ) is expressed as

$$\begin{eqnarray} \left(\psi -\lambda \right) &=&\sigma \left[ \frac{1}{\sin I}\sin \mu \right]\nonumber\\ && +\sigma ^{2}\left[ -\frac{\cos I}{\sin ^{2}I}\sin \mu \cos \mu \right] +O(\sigma ^{3}) \nonumber \\ &=&\left(\psi -\lambda \right) ^{\left(1\right) }+\left(\psi -\lambda \right) ^{\left(2\right) }+O(\sigma ^{3}). \end{eqnarray} \tag{ 71 }$$

In the same way, we write the functions (θ − I), (ϕ − μ − ν), etc. Let us recall that the σ zero-order part of the variables ψ, θ, and ϕ, denoted as ψ⁽⁰⁾, θ⁽⁰⁾, and ϕ⁽⁰⁾, coincides, respectively, with the Andoyer variables λ, I, and μ + ν. In addition, in the case of the functions V and $\mathcal {W}_{1}$, we have also separated, when necessary, f⁽ⁱ⁾ into different parts according to its dependence with the coefficients B_i, C_i,τ, and D_i,τ. With this nomenclature the gravitational potential energy (Equation (5)) is represented as

$$\begin{eqnarray} V &=&{\sum _{p=S,M} }k_{p}^{\prime }\sum _{i}\left\lbrace \left[ B_{i}\cos \Theta _{i}\right] -\sigma \left[ \sum _{\tau =\pm 1}C_{i,\tau }\cos (\mu -\tau \Theta _{i})\right] \right. \nonumber \\ &+ &\left.\!\!\sigma ^{2}\left[\! -\frac{3}{2}B_{i}\cos \Theta _{i}+\frac{1}{4} \sum _{\tau =\pm 1}\! D_{i,\tau }\cos (2\mu -\tau \Theta _{i})\!\right] \right\rbrace +O(\sigma ^{3}) \nonumber \\ &=&V^{(0)}+V^{(1)}+[ V^{(2_{B})}+V^{(2_{D})}] +O(\sigma ^{3}). \end{eqnarray} \tag{ 72 }$$

A similar representation has been used for the first-order generating function $\mathcal {W}_{1}$.

It is interesting to perform some comparison with the Hamiltonian treatment given in Souchay et al. (1999), who re-calculated second-order contributions to catch all the coupling effects down to 0.1 μas instead of 5 μas, as it was performed in the precedent work by Kinoshita & Souchay (1990). Anyway, this comparison is not a direct task since we are faced with the difficulty emerging from the lack of explicit analytical formulae for the different aspects of the Earth's rotational motion. In particular, although in Souchay et al. (1999) some analytical formulae are given, those expressions are not presented in their final, simplified form as is the case in Kinoshita (1977) or in this work, but instead only some intermediate formula that could allow the numerical computation of the final results with the aid of specific software procedures is provided. On the other hand, in these works the authors focused their studies solely on the motion of the equatorial plane, in contrast to the second-order analytical development of the complete rotational motion that we have reported here. However, we can discuss some general aspects of both approaches and numerically compare some second-order amplitudes concerning the motion of the equatorial plane. In this regard, let us point out that to numerically compute the amplitudes presented in this investigation, we have employed the values of the parameters listed in Tables 8, 9, and 10, which have been obtained from Kinoshita & Souchay (1990), Souchay et al. (1999), and references therein.

With respect to the considerations of general nature, let us underline that, as we have pointed out in the Introduction, in this work we have limited ourselves to the investigation of the spin–spin coupling, sometimes referred to as the crossed-nutation effect (Souchay et al. 1999); that is to say, we have only been concerned with the rotational problem without considering the coupling effect between the orbital motion of the Moon and the rotational motion of the Earth, i.e., the spin–orbit coupling effect. In this sense, this work can be considered a comprehensive extension of the theory of the rotation of the rigid Earth developed in Kinoshita (1977) at the second order. Hence, our dynamical system is described with the help of six canonical variables (Andoyer variables), instead of the twelve canonical variables that would arise from adding the six Delaunay canonical variables associated with the orbital motion of the Moon. This fact sometimes entails further complications for numerical comparisons, since some second-order contributions merge these two effects; for example, this is the case for the precession rate in longitude (Kinoshita & Souchay 1990).

Regarding the spin–spin coupling treatment, there is an important difference when implementing the perturbation procedure between this investigation and that by Souchay et al. (1999): we include the terms E and $V_{\rm {sec}}$ in the secular part of the first-order Hamiltonian $\mathcal {H}_{1}$, whereas Souchay et al. (1999), as well as Kinoshita & Souchay (1990) and Kinoshita (1977), seem to incorporate⁶ E in the unperturbed part of the Hamiltonian $\mathcal {H}_{0}$ and $V_{\rm {sec}}$ in the secular part of the first-order Hamiltonian $\mathcal {H}_{1}$. From our point of view, our approach makes the effects of the moving ecliptic on the rotation of the Earth clear, both for their in-phase and out-phase components, and allows a more direct computation of the generating function since the auxiliary angle I is strictly constant. In any case, the numerical results for the out-of-phase components of the leading nutation term agree with those reported by Souchay et al. (1999), although we find more terms beyond the threshold of 0.1 μ as established by those authors, as it can be seen in Table 2. Let us recall that the existence of those out-of-phase terms was established in Williams (1994), and as we can deduce immediately from our analytical formulae, their origin is due to the moving ecliptic and not to any kind of dissipative effect, absent on this model.

Table 2. Second-order Out-of-phase Nutations of the Andoyer Plane Due to the Ecliptic Term E and Comparison With Souchay et al. (1999), REN–2000

Argument					Period	Δλ(cos )		ΔI(sin )
lM	lS	F	D	Ω	Days	$\mathcal {L}_{E}^{\rm out}$	REN–2000	$\mathcal {O}_{E}^{\rm out}$	REN-2000
+0	+0	+0	+0	+1	−6798.36	−247.17	−248	−26.79	−27
+0	+0	+0	+0	+2	−3399.18	−1.30	...	−0.70	...
+0	+1	+0	+0	+0	365.26	0.09	...	0.00	...
+0	+0	+2	−2	+2	182.62	−0.43	...	−0.23	...

Notes. The unit of the amplitudes is 1 μas, having listed the terms whose values are equal or greater than 0.05 μas. In all the tables displayed in this section, the given values correspond to the nutations themselves, their sign being not adapted to the conventional usage (to match them, the sign of each term should be reversed). In a similar way, the amplitude denote all the expressions necessary to compute the final numerical values. For example, with $\mathcal {L}_{E}^{\rm out}$ we denote the amplitude arising from each argument in the formula $\sum _{p=S,M}\frac{k_p}{\sin I} \sum _{i\ne 0}\mathcal {L}_{E}^{\rm out}\cos \Theta _i$ of the angle Δλ.

Download table as:  ASCIITypeset image

Likewise, as a consequence of the inclusion of $V_{\rm {sec}}$ and E as first-order perturbations, the mean motion of the variable λ in the unperturbed problem is exactly 0. It means that the frequency n_i = dΘ_i/dt is independent of the canonical momenta (compare with Equation (6.17) in Kinoshita 1977) and must be taken as constant when computing the Poisson brackets. In our treatment, the mean motion of λ, and also of I, couples with the rotation at the second order, hence the amplitudes $\mathcal {L}_{E}^{\rm in}$, $\mathcal {O}_{E}^{\rm in}$, $\mathcal {L}_{2}^{s}$, and $\mathcal {O}_{2}^{s}$ are equivalent to the nutations derived from a first-order generating function obtained by considering the mean motions of λ and I, which would arise if the terms E and $V_{\rm {sec}}$ were incorporated in the unperturbed Hamiltonian $\mathcal {H}_{0}$. In the case of the first-order periodic variation of λ, this equivalence is established if we expand Equation (6.17) in Kinoshita (1977) in powers of the products k_pk_q and k_pe_i, recovering the amplitudes $\mathcal {L}_{2}^{s}$ and $\mathcal {L} _{E}^{\rm in}$. In a similar way, it would also be the origin of the amplitudes $\mathcal {N}_{E}^{\rm in}$ and $\mathcal {N}_{2}^{s}$ of the short-period variations of (ψ − λ), (θ − I), (ϕ − μ − ν), and U.

Yet another very substantial difference of this investigation with respect to the spin–spin coupling studied in Souchay et al. (1999) is the treatment of the auxiliary variable σ in the expression of the angles ψ and θ as functions of the Andoyer variables, in the expression of the gravitational potential energy V and of the generating function $\mathcal {W}_{1}$, and in the computation of the Poisson brackets. Specifically, in Souchay et al. (1999), second-order contributions of the motion of the equatorial plane are constructed considering

$$\begin{eqnarray} \theta &=&I+O(\sigma),\,\quad \psi =\lambda +O(\sigma), \nonumber \\ V &=&{\sum _{p=S,M} }k_{p}^{\prime }\sum _{i}B_{i}\cos \Theta _{i}+O(\sigma),\,\,\nonumber\\ \mathcal {W}_{1} &=&{\sum _{p=S,M} }k_{p}^{\prime }\sum _{i\ne 0}\frac{B_{i}}{n_{i}}\sin \Theta _{i}+O(\sigma), \nonumber \\ \left\lbrace f;\mathcal {H}\right\rbrace &\simeq &\frac{\partial f}{\partial \lambda } \frac{\partial \mathcal {H}}{\partial \Lambda }-\frac{\partial \mathcal {H}}{ \partial \lambda }\frac{\partial f}{\partial \Lambda }. \end{eqnarray} \tag{ 73 }$$

This formulation is inherited from Kinoshita & Souchay (1990) and entails the identification of the rotational motion, at the second order, of the equatorial plane with that of the Andoyer plane, and according to the notation introduced previously, the reduction of V to V⁽⁰⁾ and of $\mathcal {W}_{1}$ to $\mathcal {W}_{1}^{\left(0\right) }$. Therefore, looking at Table 1, it is followed that in our treatment these second-order contributions correspond with the amplitudes $\mathcal {L}_{2}^{p1}$, $\mathcal {O}_{2}^{p1}$, $\mathcal {L}_{3}^{1}$, and $\mathcal {O}_{3}^{1}$. Since, as we have explained above, Souchay et al. (1999) do not provide explicit, final formulae for these terms, we can only compare them numerically. In Table 3, we have displayed the total contributions of these amplitudes, as well as those reported in Souchay et al. (1999). As it can be seen there is a good agreement between both approaches, although some terms present discrepancies probably due to the fact that Souchay et al. (1999) considered a large set of arguments Θ_i in the expansion of the orbital motion of the Moon and the Sun.

Table 3. Second-order Nutations of the Andoyer Plane: V⁽⁰⁾_per–$\mathcal {W}_{1}^{(0)}$ and $\mathcal {W}_{1}^{(0)}$–$\mathcal {W}_{1}^{(0)}$ Contributions and Comparison with Souchay et al. (1999), REN–2000

Argument					Period	Δλ(sin )				ΔI(cos )
lM	lS	F	D	Ω	Days	$\mathcal {L}_{2}^{p1}$	$\mathcal {L}_{3}^{1}$	Total	REN–2000	$\mathcal {O}_{2}^{p1}$	$\mathcal {O}_{3}^{1}$	Total	REN–2000
+0	+0	+0	+0	+1	−6798.38	−18.01	−11.63	−29.65	−28.80	+20.68	+7.92	+28.60	+27.70
+0	+0	+0	+0	+2	−3399.19	*	−1212.65	−1212.65	−1220.60	*	+236.47	+236.47	+238.10
+0	+0	+0	+0	+3	−2266.13	−1.00	+22.31	+21.31	+21.60	+0.43	−4.18	−3.75	−4.20
+0	+0	+0	+0	+4	−1699.60	*	−0.11	−0.11	−0.10	*	+0.02	+0.02
+0	+1	+0	+0	+1	+386.00	−2.24	+3.23	+0.98	+1.00	+1.56	−1.40	+0.16	+0.40
+0	+1	−2	+2	−3	−385.96	+0.34	−2.29	−1.95	−1.90	+0.15	−0.43	−0.28	−0.50
+0	+1	+0	+0	+0	+365.26	+1.02	*	+1.02	+1.00	*	−0.11	−0.11	−0.20
+0	+1	−2	+2	−2	−365.22	+1.05	+0.42	+1.47	+1.50	+0.50	+0.17	+0.67	+1.70
+0	+1	+0	+0	−1	+346.64	−1.81	+3.23	+1.42	+1.40	−1.26	+1.40	+0.14	+0.20
+0	+1	−2	+2	−1	−346.60	+0.34	+1.20	+1.53	+1.50	+0.39	+0.81	+1.20	+0.80
+0	+0	+2	−2	+4	+192.99	*	+1.40	+1.40	+1.40	+0.00	−0.28	−0.28	−0.30
+0	+0	+2	−2	+3	+187.66	+19.43	−136.75	−117.32	−117.70	−8.42	+25.65	+17.23	+17.30
+0	+0	+2	−2	+2	+182.62	−0.03	−0.02	−0.05	+1.70	+0.01	+0.01	+0.02	−0.30
+0	+0	+2	−2	+1	+177.84	+21.17	+71.45	+92.62	+92.80	−24.31	−48.65	−72.96	−73.10
+0	+0	+2	−2	+0	+173.31	−1.04	*	−1.04	−0.40	*	+0.83	+0.83	+1.00
+0	+1	+2	−2	+4	+126.27	+0.00	+0.05	+0.05	–	*	−0.01	−0.01	–
+0	+1	+2	−2	+3	+123.97	+0.74	−5.33	−4.59	−4.60	−0.32	+1.00	+0.68	+1.00
+0	+1	+2	−2	+2	+121.75	+0.12	+0.42	+0.54	+0.60	−0.06	−0.17	−0.22	−0.20
+0	+1	−4	+4	−4	−121.74	+0.00	−0.14	−0.14	−0.10	−0.00	−0.03	−0.03
+0	+1	+2	−2	+1	+119.61	+0.84	+2.78	+3.62	+3.60	−0.96	−1.90	−2.86	−3.80
+0	+0	+4	−4	+4	+91.31	*	−4.27	−4.27	−4.30	+0.00	+0.85	+0.85	+0.80
+0	+1	+4	−4	+4	+73.05	−0.00	−0.33	−0.33	−0.30	+0.00	+0.07	+0.07
+1	+0	+0	+0	+1	+27.67	−1.09	+1.74	+0.65	+0.60	+0.75	−0.75	−0.00
+1	+0	+0	+0	+0	+27.56	+0.12	*	+0.12	+0.80	*	−0.02	−0.02
+1	+0	+0	+0	−1	+27.44	−1.07	+1.74	+0.66	+0.70	−0.75	+0.75	+0.01
+1	+0	−2	+0	−2	−27.09	+0.09	+0.04	+0.13	−0.30	+0.04	+0.01	+0.06
+1	+0	+2	−2	+2	+23.94	−0.14	+0.23	+0.09	+0.10	+0.07	−0.09	−0.02
+0	+1	+0	−2	+0	−15.39	−0.06	−0.00	−0.06	–	*	−0.03	−0.03	–
+0	+0	+0	+2	+0	+14.76	+1.32	*	+1.32	+1.30	*	−0.82	−0.82	−0.80
+0	+0	+0	+2	−1	+14.73	+0.05	−0.14	−0.09	–	+0.06	−0.10	−0.04	–
+0	+1	−2	+0	−2	−14.19	+0.06	+0.07	+0.13	–	+0.03	+0.03	+0.06	–
+0	+0	+2	+0	+4	+13.72	*	+0.22	+0.22	+0.20	*	−0.04	−0.04
+0	+0	+2	+0	+3	+13.69	+2.97	−21.88	−18.91	−19.00	−1.29	+4.10	+2.82	+2.80
+0	+0	+2	+0	+2	+13.66	−0.01	−4.82	−4.83	−4.80	+0.00	+0.94	+0.94	+0.90
+0	+0	+2	+0	+1	+13.63	+3.57	+11.45	+15.02	+15.20	−4.10	−7.80	−11.90	−12.00
+0	+0	+2	+0	+0	+13.61	−2.13	+0.00	−2.13	−2.10	*	−0.45	−0.45	−0.50
+0	+1	+2	+0	+2	+13.17	+0.05	+0.07	+0.12	–	−0.02	−0.03	−0.05	–
+0	+0	+4	−2	+4	+12.71	*	−1.37	−1.37	−1.40	−0.00	+0.27	+0.27	+0.30
+0	+0	+4	−2	+3	+12.69	−0.03	−0.27	−0.30	−0.30	+0.01	+0.05	+0.06
+0	+1	+4	−2	+4	+12.28	+0.00	−0.05	−0.05	–	+0.00	+0.01	+0.01	–
+1	+0	+0	+2	+0	+9.61	+0.16	*	+0.16	−0.20	*	−0.10	−0.10	−0.10
+1	+0	+2	+0	+3	+9.14	+0.38	−2.81	−2.43	−2.40	−0.16	+0.53	+0.36	+0.60
+1	+0	+2	+0	+1	+9.12	+0.46	+1.47	+1.93	+1.90	−0.53	−1.00	−1.53	−1.50
+1	+0	+4	−2	+4	+8.70	+0.00	−0.18	−0.18	−0.20	−0.00	+0.03	+0.03
+0	+0	+4	+0	+4	+6.83	*	−0.11	−0.11	−0.10	*	+0.02	+0.02

Notes. The unit of the amplitudes is 1 μas, having listed the terms whose values are equal or greater than 0.05 μas. The symbol "*" represents the fact that, with the internal accuracy used in the computations, the amplitude does not appear in the corresponding argument. In addition, we have performed a comparison with the values given in Souchay et al. (1999). Let us notice that in Souchay et al. (1999) the amplitudes are given with a precision of 0.1 μas (Table 1 in Souchay et al. 1999) and probably computed with a complete list of arguments Θ_i in the expansion of the gravitational potential energy of the J₂ term, instead of the eleven arguments used in this work. The characters "–" reflects that those values are not given in Souchay et al. (1999).

Download table as:  ASCIITypeset image

With respect to the long-period evolution of the angles ψ and θ, that is to say, the precession, it is computed with the same premises as before in Kinoshita & Souchay (1990), so it would correspond with the amplitude $\mathcal {S}_{2a}^{L}$ of Table 1. Regrettably, in this case we cannot perform a direct numerical comparison, since the numerical contribution reported in that work also contains a part coming from the spin–orbit coupling effect. However, this effect is partially worked out in Kinoshita (1977), who only considered the largest term with the argument Θ_i = Ω in the development of V and $\mathcal {W}_{1}$, providing an analytical expression for the precession rate in longitude. If we particularize the amplitude $\mathcal {S}_{2a}^{L}$ to this situation (consider Equation (51) and adapt the notations), we recover the same expression as given in that investigation.⁷

It must be emphasized, however, that the simplifications assumed in Equation (73), based on the smallness of both the angle σ and second-order contributions, depend drastically on the truncation level established for the amplitude of the nutation terms. In this way, if the threshold is fixed in 5 μas, as it was established in Kinoshita & Souchay (1990), the simplification is valid, whereas if we bring it down up to 0.1 μas (Souchay et al. 1999), we would be omitting some nutation amplitudes whose value is greater than 0.1 μas, as we will show below. In any case, it is difficult to give a priori some method that ensures that those kind of simplifications allow catching all the amplitudes with a definite truncation level.⁸ Therefore, a precise determination of the amplitudes requires an explicit computation of the different contributions, as we have performed in this investigation.

By so doing, to catch up all the final second-order contributions arising from the spin–spin coupling at order O(σ⁰), it is necessary to remove the simplifications assumed in Equation (73) by Souchay et al. (1999). In particular, we must keep the expansions up to the second order in σ in the Euler–Andoyer relationships (Equation (18)), in the generating function, and in the gravitational potential energy. It is a similar situation as in the case of a first-order theory, where the computation of the Oppolzer terms (Kinoshita 1977) needs to keep those expansions up to the first order in σ. The reason is clear if we consider the following property of the Poisson brackets

$$\begin{eqnarray} &&\lbrace \sigma ^{n}h_{1}\left(\mu \right),\sigma ^{m}h_{2}\left(\mu \right) \rbrace = \frac{1}{M}\Bigg(mh_{2}\left(\mu \right) \frac{\partial h_{1}\left(\mu \right) }{\partial \mu }\nonumber\\ &&\qquad -nh_{1}\left(\mu \right) \frac{\partial h_{2}\left(\mu \right) }{\partial \mu }\Bigg) \sigma ^{n+m-2}+O(\sigma ^{n+m}), \end{eqnarray} \tag{ 74 }$$

where n, m are two non-negative integers with n + m ⩾ 2. Therefore, taking into account that some second-order contributions (Equation (27)) arise from computing a Poisson bracket of a Poisson bracket of the form

$$\begin{equation} \left\lbrace \left\lbrace s_{1};s_{2}\right\rbrace ;s_{3}\right\rbrace, \end{equation} \tag{ 75 }$$

it is followed that the combination of σ² terms in s₁ with σ terms in s₂, and vice versa, gives rise to a term of the order O(σ¹), whose combination with the σ part of s₃ produces a final contribution at order O(σ⁰). Hence, it is necessary to go beyond the zero-order terms in σ in the functions s₁, s₂, and s₂. This fact explains the origin of the amplitudes $\mathcal {L}_{2}^{p2}$, $\mathcal {O}_{2}^{p2}$, $\mathcal {M} _{2}^{p2}$, $\mathcal {L}_{3}^{2}$, $\mathcal {O}_{3}^{2}$, and $\mathcal {M} _{3}^{2}$, and of all the amplitudes of the form $\mathcal {N}_{\alpha }^{\beta }$ (see Table 1). It is interesting to point out that even a precise computation of the second-order contributions of the motion of the Andoyer plane requires the consideration of the terms proportional to σ in V and $\mathcal {W}_{1}$. It is an important difference with respect to the first-order situation, where all the final contributions, at order O(σ⁰), to the motion of the Andoyer plane come from V⁽⁰⁾ and $\mathcal {W}_{1}^{\left(0\right) }$. In a similar way, some second-order contributions to the motion of the equatorial plane need the consideration of the terms proportional to σ² in the relationships between Euler and Andoyer variables (Equation (18)), whereas in the first-order case it was only necessary to consider the terms proportional to σ (Kinoshita 1977). Summarizing up, the additional second-order contributions to the motion of the equatorial plane have a double origin: one part comes from the additional second-order contributions to the motion of the Andoyer plane, the other one arises from the second-order contributions to the motion of the equatorial plane itself with respect to the Andoyer plane. In Table 4, we show the different amplitudes of the motion of the Andoyer and equatorial planes indicating whether or not they are included in Souchay et al. (1999).

Table 4. Different Amplitudes of the Rotational Motion of the Andoyer Plane and the Equatorial Plane at the Second-order and Comparison with Souchay et al. (1999), REN–2000

Plane	Amplitude	REN–2000
Andoyer (precession)	$\mathcal {S}_{2a}^{L}$	Yes
Andoyer (precession)	$\mathcal {S}_{2b}^{L}$	No
Andoyer	$\mathcal {L}_{E}^{\rm in},\,\mathcal {O}_{E}^{\rm in},\,\mathcal {L}_{E}^{\rm out},\,\mathcal {O}_{E}^{\rm out},\,\mathcal {L}_{2}^{s},\,\mathcal {O}_{2}^{s}$	Yes (as first order)
Andoyer	$\mathcal {L}_{2}^{p1},\,\mathcal {O}_{2}^{p1},\,\mathcal {L}_{3}^{1},\,\mathcal {O}_{3}^{1}$	Yes
Andoyer	$\mathcal {L}_{2}^{p2},\,\mathcal {O}_{2}^{p2},\,\mathcal {L}_{3}^{2},\,\mathcal {O}_{3}^{2}$	No
Equatorial	$\mathcal {N}_{E}^{\rm in},\,\mathcal {N}_{E}^{\rm out},\,\mathcal {N}_{2}^{s}$	Yes (as first order)
Equatorial	$\mathcal {N}_{2}^{p1},\,\mathcal {N}_{2}^{p2}$	No
Equatorial	$\mathcal {N}_{3}^{1L},\,\mathcal {N}_{3}^{1O},\,\mathcal {N}_{3}^{2},\,\mathcal {N}_{3}^{3L},\,\mathcal {N}_{3}^{3O}$	No

Notes. The rows referred to the equatorial plane gives the motion of this plane with respect to the Andoyer one. As it is displayed in the table, some second-order amplitudes of this work are included as first-order amplitudes in Souchay et al. (1999) and in Kinoshita (1977). It is a consequence of the inclusion of V_sec and E in the secular part of the first-order Hamiltonian $\mathcal {H}_1$. The remaining amplitudes determined in this work relative to the evolution of the angle ϕ and the polar motion are not tackled in Souchay et al. (1999).

Download table as:  ASCIITypeset image

We have computed the numerical contributions of those new additional second-order terms to the motion of the equatorial plane. As expected, their amplitudes are small, although some of them are greater than the truncation level of 0.1 μas established in Souchay et al. (1999). In particular, as it is displayed in Table 5, we have found four additional terms greater than 0.05 μas that contribute at the second order to the motion of the Andoyer plane (three terms in longitude and three terms in obliquity exceed 0.1 μas). Those contributions exclusively come from the amplitudes $\mathcal {L} _{2}^{p2}$ and $\mathcal {O}_{2}^{p2}$, the values of $\mathcal {L}_{2}^{3}$ and $\mathcal {O}_{2}^{3}$ being negligible. This fact can be partially understood from a qualitative point of view by considering the analytical expressions of those functions (Equation (D5)) since, unlike $\mathcal {L} _{2}^{p2}$and $\mathcal {O}_{2}^{p2}$, in the denominator both frequencies n_i and n_j are summed to the mean motion n_μ, so, roughly speaking, those terms are smaller than the amplitudes $\mathcal {L}_{2}^{p2}$ and $\mathcal {O}_{2}^{p2}$ by the factor

$$\begin{equation} \frac{\tau n_{i}-\rho n_{j}}{n_{\mu }-\rho n_{j}}\simeq \frac{n_{k}}{\omega _{E}}, \end{equation} \tag{ 76 }$$

where n_k is the frequency associated to the argument τΘ_i − ρΘ_j, and we have taken into account that n_μ ≃ ω_EC/A ≃ ω_E. This factor is small for most of the terms, since the frequencies in the expansion of the orbital motion of the Moon and the Sun are greater than the mean value of the angular velocity of the Earth. In a similar way, the numerical values of the amplitudes $\mathcal {N}_{\alpha }^{\beta }$, which determine the second-order motion of the equatorial plane with respect to the Andoyer plane, are shown in Table 6. We have found 10 additional terms greater than 0.05 μas (nine terms in longitude and eight terms in obliquity exceed 0.1 μas). The contributions arise from the amplitudes $\mathcal {N}_{2}^{p1}$, $\mathcal {N}_{3}^{1L}$, and $\mathcal {N}_{3}^{1O}$, whereas the remaining terms $\mathcal {N}_{2}^{p2}$, $\mathcal {N}_{3}^{2}$, $\mathcal {N}_{3}^{3L}$, and $\mathcal {N}_{3}^{3O}$ do not provide any significant contribution. The explanation is similar to the case of the Andoyer plane, as it can be derived from the analytical expressions of these amplitudes (Equation (E4)).

Table 5. Second-order Nutations of the Andoyer Plane: V⁽¹⁾_per–$\mathcal {W}_{1}^{(1)}$ and $\mathcal {W}_{1}^{(1)}$–$\mathcal {W}_{1}^{(1)}$ Contributions

Argument					Period	Δλ(sin )			ΔI(cos )
lM	lS	F	D	Ω	Days	$\mathcal {L}_{2}^{p2}$	$\mathcal {L}_{3}^{2}$	Total	$\mathcal {O}_{2}^{p2}$	$\mathcal {O}_{3}^{2}$	Total
+0	+0	+0	+0	+1	−6798.38	−0.45	−0.00	−0.45	+1.37	+0.00	+1.37
+0	+0	+0	+0	+2	−3399.19	+0.06	−0.00	+0.06	−0.03	+0.00	−0.03
+0	+0	+2	−2	+2	+182.62	−0.20	+0.00	−0.20	+0.11	−0.00	+0.11

Notes. The unit of the amplitudes is 1 μas, having listed the terms whose values are equal or greater than 0.05 μas. These terms are not present in Souchay et al. (1999).

Download table as:  ASCIITypeset image

Table 6. Second-order Nutations of the Equatorial Plane with Respect to the Andoyer Plane

Argument					Period	Δ(ψ − λ)(sin )						Δ(θ − I)(cos )
lM	lS	F	D	Ω	Days	$\mathcal {N}_{2}^{p1}$	$\mathcal {N}_{3}^{1L}$	$\mathcal {N}_{2}^{p2}$	$\mathcal {N}_{3}^{2}$	$\mathcal {N}_{3}^{3L}$	Total	$\mathcal {N}_{2}^{p1}$	$\mathcal {N}_{3}^{1O}$	$\mathcal {N}_{2}^{p2}$	$\mathcal {N}_{3}^{2}$	$\mathcal {N}_{3}^{3O}$	Total
+0	+0	+0	+0	+1	−6798.38	−0.25	+0.26	+0.00	+0.00	−0.00	+0.01	+0.36	+0.36	−0.00	−0.00	+0.00	+0.72
+0	+0	+0	+0	+2	−3399.19	+0.08	+0.27	−0.00	−0.00	−0.00	+0.35	−0.04	−0.06	+0.00	+0.00	−0.00	−0.10
+0	+0	+2	−2	+3	+187.66	−0.08	−0.53	+0.00	+0.00	+0.00	−0.61	+0.04	+0.07	−0.00	−0.00	+0.00	+0.11
+0	+0	+2	−2	+1	+177.84	+0.47	+0.16	−0.00	−0.00	−0.00	+0.63	−0.16	−0.19	+0.00	+0.00	−0.00	−0.35
+0	+0	+0	+2	+0	+14.76	+0.06	+0.04	+0.00	−0.00	+0.00	+0.10	−0.02	−0.04	−0.00	−0.00	+0.00	−0.06
+0	+0	+2	+0	+3	+13.69	−0.19	−1.26	+0.00	+0.00	+0.00	−1.45	+0.09	+0.17	−0.00	−0.00	+0.00	+0.26
+0	+0	+2	+0	+2	+13.66	−0.09	−0.29	−0.00	−0.00	+0.00	−0.37	+0.04	+0.06	+0.00	+0.00	+0.00	+0.10
+0	+0	+2	+0	+1	+13.63	+1.07	+0.37	−0.00	−0.00	+0.00	+1.44	−0.36	−0.44	+0.00	+0.00	−0.00	−0.80
+1	+0	+2	+0	+3	+9.14	−0.04	−0.25	+0.00	+0.00	+0.00	−0.29	+0.02	+0.03	−0.00	−0.00	+0.00	+0.05
+1	+0	+2	+0	+1	+9.12	+0.21	+0.07	−0.00	−0.00	+0.00	+0.28	−0.07	−0.09	+0.00	+0.00	−0.00	−0.16

Notes. The unit of the amplitudes is 1 μas, having listed the terms whose values are equal or greater than 0.05 μas. These terms are not present in Souchay et al. (1999).

Download table as:  ASCIITypeset image

Consequently, to reach the truncation level of 0.1 μas established in Souchay et al. (1999), in order to model the rotational motion of the equatorial plane at the second order, we should take into account the numerical contributions related with the Andoyer plane (Table 5) and those ones related with the equatorial plane (Table 6). Their sum is presented in Table 7 following the same criteria as in Souchay et al. (1999). Let us note that we have found three terms in longitude greater or equal than 0.1 μas and nine terms in obliquity. In addition, we should also consider the contribution to the precession rate in longitude due to the amplitude $\mathcal {S}_{2b}^{L}$, which is not present in Kinoshita & Souchay (1990), and whose value is 0.692 mas cy⁻¹.

However, from the perspective of extending this investigation to more realistic Earth models, i.e., non-rigid models, the importance of the new terms, derived from removing the simplifications taken in Equation (73) by Souchay et al. (1999), is much more relevant. It is due to the fact that the amplitudes $\mathcal {L}_{2}^{p2},$ $\mathcal {O} _{2}^{p2}$, $\mathcal {L}_{2}^{3},$ $\mathcal {O}_{2}^{3}$, $\mathcal {N} _{2}^{p1}$, $\mathcal {N}_{3}^{1L}$, $\mathcal {N}_{3}^{1O}$, $\mathcal {N} _{2}^{p2}$, $\mathcal {N}_{3}^{2}$, $\mathcal {N}_{3}^{3L}$, and $\mathcal {N} _{3}^{3O}$ (Equations (D5) and (E4)) depend on the mean motion n_μ, what reflects that they are influenced by the internal mechanical properties of the Earth, since n_μ is related with the Euler-free frequency. Hence, when considering a non-rigid Earth model—in particular, a two-layer Earth model composed of a mantle that encloses a fluid core—these terms might play an important role since the internal structure of the Earth could amplify their amplitudes due to the resonance caused by the fluid core. As a matter of fact, in the case of the precessional motion, it has been shown in Ferrándiz et al. (2004) that the presence of the fluid core does affect the precession rate in longitude, contrary to the extended statement that holds that the internal structure of the Earth does not affect this precessional motion. Specifically, the rigid contribution coming from $\mathcal {S}_{2b}^{L}$, whose value is 0.692 mas cy⁻¹, is amplified in the above-mentioned non-rigid Earth model providing a value of about 23 mas cy⁻¹ (Ferrándiz et al. 2004).

Therefore, the availability of a systematic second-order theory for the rotation of a non-rigid Earth becomes necessary, since that model might provide new second-order contributions that cannot be longer ignored within nowadays accuracy requirements. In this first stage, we are only involved with the second-order effects arising from the rotation–rotation interaction. Even in this situation, to extend the present investigation to the non-rigid Earth, it is necessary to develop new and complex mathematical methods, since the intrinsic difficulties of the non-rigid model are much greater than in the rigid case. In this sense, the analytical formulae that we have obtained in this investigation will be a benchmark to check and compare the correctness of the expressions derived in the non-rigid case. This work is in progress and will be presented in a forthcoming paper.

We acknowledge an anonymous referee for critical comments that improved the manuscript. This work has been partially supported by the Spanish projects I+D+I, AYA2004-07970, and AYA2007-67546; and Junta de Castilla y León project VA070A07. Its final version was carried out while A. Escapa was on sabbatical leave from the University of Alicante at the National Astronomical Observatory of Japan (NAOJ), supported by the Spanish Ministerio de Educación, project PR2009-0379, within the Programa Nacional de Movilidad de Recursos Humanos I−D+i 2008–2011. The generous hospitality of the NAOJ staff is gratefully acknowledged. D. Miguel's contributions form part of his PhD thesis under the supervision of J. Getino and A. Escapa at the University of Valladolid.

The relationship between the Euler and Andoyer variables can be derived by considering that the rotation matrix that brings the celestial frame to the terrestrial frame can be parameterized in both sets (Figure 1). In particular, we have

$$\begin{equation} R_{3}(\phi) R_{1}(\theta) R_{3}(\psi) =R_{3}(\nu) R_{1}(\sigma) R_{3}(\mu) R_{1}(I) R_{3}(\lambda), \end{equation} \tag{ A1 }$$

where the rotation (counterclockwise) about the first (x) axis or the third (z) axis through an angle ϑ is given, respectively, by

$$\begin{eqnarray} R_{1}&=&\left(\begin{array}{@{}lll@{}} 1 & 0 & 0 \\ 0 & \cos \vartheta & \sin \vartheta \\ 0 & -\sin \vartheta & \cos \vartheta \end{array} \right),\nonumber\\ R_{3}&=&\left(\begin{array}{@{}lll@{}} \cos \vartheta & \sin \vartheta & 0 \\ -\sin \vartheta & \cos \vartheta & 0 \\ 0 & 0 & 1 \end{array} \right). \end{eqnarray} \tag{ A2 }$$

This equality can be cast into the form R₃(ϕ − ν)R₁(θ)R₃(ψ − λ) = R₁(σ)R₃(μ)R₁(I), giving different equations that implicitly define the Euler angles as functions of the Andoyer set. For example, we have

$$\begin{eqnarray} \cos \theta & = & \cos I\cos \sigma -\sin I\sin \sigma \cos \mu, \nonumber\\ \sin \theta \sin (\psi -\lambda) & = & \sin \sigma \sin \mu, \nonumber\\ \sin \theta \cos (\psi -\lambda) & = & \sin I\cos \sigma +\cos I\sin \sigma \cos \mu, \nonumber\\ \sin \theta \sin (\phi -\nu) & = & \sin I\sin \mu, \nonumber\\ \sin \theta \cos (\phi -\nu) & = & \cos I\sin \sigma +\sin I\cos \sigma \cos \mu. \nonumber\\ \end{eqnarray} \tag{ A3 }$$

From them, we can obtain the explicit relations

$$\begin{eqnarray} \theta & = & \arccos \left[ \cos I\cos \sigma -\sin I\sin \sigma \cos \mu \right], \nonumber\\ \psi & = & \lambda +\arctan \left[ \displaystyle \ \frac{\sin \sigma \sin \mu }{\sin I\cos \sigma +\cos I\sin \sigma \cos \mu }\right], \nonumber\\ \phi & = & \nu +\arctan \left[ \displaystyle \ \frac{\sin I\sin \mu }{\cos I\sin \sigma +\sin I\cos \sigma \cos \mu }\right]. \end{eqnarray} \tag{ A4 }$$

Let us note that when σ goes to 0 these equalities provide

$$\begin{equation} \theta \rightarrow I,\quad \psi \rightarrow \lambda,\quad \phi \rightarrow \mu +\nu. \end{equation} \tag{ A5 }$$

Since the numerical value of the variable σ is small (about 10⁻⁶ rad), instead of employing directly the Equations (A4) we will consider their Taylor series expansion around the value σ = 0, what allows expressing the explicit dependence of the Euler angles on the Andoyer variables in a simpler way.

With the help of a symbolic processor of general purposes like, for example, Maple or Mathematica, we can derive those expansions up to the second order. They turn out to be

$$\begin{eqnarray} \theta &=&I+\sigma \cos \mu +\sigma ^{2}\frac{\cos I}{2\sin I}\sin ^{2}\mu +O(\sigma ^{3}), \nonumber \\ \psi &=&\lambda +\sigma \frac{1}{\sin I}\sin \mu -\sigma ^{2}\frac{\cos I}{ \sin ^{2}I}\sin \mu \cos \mu +O(\sigma ^{3}), \nonumber \\ \phi &=&\mu +\nu -\sigma \frac{\cos I}{\sin I}\sin \mu\nonumber\\ &&+\sigma ^{2}\frac{ (1+\cos ^{2}I) }{2\sin ^{2}I}\sin \mu \cos \mu +O(\sigma ^{3}). \end{eqnarray} \tag{ A6 }$$

The coefficients of the Poisson bracket $\mathcal {C}_{S}=\left\lbrace \mathcal {H} _{1\,{\rm sec} };W_{1}\right\rbrace$ (Section 3.3) are given by

$$\begin{eqnarray} \mathcal {C}_{s0}^{\rm in} &=&\left[ e_{1}-(e_{2}\cos \lambda +e_{3}\sin \lambda) \frac{\cos I}{\sin I}\right] m_{5i}\frac{B_{i}}{n_{i}}, \nonumber \\ \mathcal {C}_{s1}^{\rm in} &=&\left[ \tau m_{5i} e_{1}+(e_{2}\cos \lambda +e_{3}\sin \lambda)\frac{1-\tau m_{5i}\cos I}{\sin I}\right]\nonumber\\ &&\times\frac{ C_{i,\tau }}{n_{\mu }-\tau n_{i}}, \nonumber \\ \mathcal {C}_{s2}^{\rm in} &=&-\!\frac{1}{4}\left[ \tau m_{5i} e_{1}+(e_{2}\cos \lambda +e_{3}\sin \lambda)\frac{2-\tau m_{5i}\cos I}{\sin I}\right]\nonumber\\ &&\times \frac{ D_{i,\tau }}{2n_{\mu }-\tau n_{i}}, \nonumber \\ \mathcal {C}_{s0}^{\rm out} &=&(e_{2}\sin \lambda -e_{3}\cos \lambda)\frac{ B_{i}^{\prime }}{n_{i}}, \nonumber \\ \mathcal {C}_{s1}^{\rm out} &=&-\!(e_{2}\sin \lambda -e_{3}\cos \lambda)\frac{ C_{i,\tau }^{\prime }}{n_{\mu }-\tau n_{i}}, \nonumber \\ \mathcal {C}_{s2}^{\rm out} &=&\frac{1}{4}(e_{2}\sin \lambda -e_{3}\cos \lambda) \frac{D_{i,\tau }^{\prime }}{2n_{\mu }-\tau n_{i}}, \nonumber \\ \mathcal {C}_{s0}^{a} &=&-\!m_{5i}\frac{B_{i}}{n_{i}}B_{0}^{\prime }, \nonumber \\ \mathcal {C}_{s1}^{a} &=&\frac{C_{i,\tau }}{n_{\mu }-\tau n_{i}}[ (\cos I-\tau m_{5i})B_{0}^{\prime }-3\sin I B_{0}], \nonumber \\ \mathcal {C}_{s2}^{a} &=&-\!\frac{1}{4}\frac{D_{i,\tau }}{2n_{\mu }-\tau n_{i}} [ (2\cos I-\tau m_{5i})B_{0}^{\prime }-6\sin I B_{0}],\nonumber\\ \end{eqnarray} \tag{ B1 }$$

where we have used the notation

$$\begin{equation} B_{i}^{\prime }=\frac{\partial B_{i}}{\partial I},\quad C_{i,\tau }^{\prime \prime }=\frac{\partial ^{2}C_{i,\tau }}{\partial I^{2}}\,\ \ldots. \end{equation} \tag{ B2 }$$

With regard to the term $\mathcal {C}_{P}=\frac{1}{2}\lbrace \mathcal {H} _{1\rm per};W_{1}\rbrace$ (Section 3.3), we have that

$$\begin{eqnarray} \mathcal {C}_{p0}^{a} &=&-\!\frac{1}{8}\frac{1}{\tau n_i}(\tau m_{5i} B_{i} B_{j}^{\prime }+\rho m_{5j} B_{i}^{\prime } B_{j}), \nonumber \\ \mathcal {C}_{p1}^{b} &=&-\frac{1}{4}\left(\frac{1}{n_{\mu }-\tau n_{i}}- \frac{1}{\rho n_{j}}\right) \nonumber \\ &&\times [ C_{i,\tau}(3\sin I\ B_{j}+(\tau m_{5i}-\cos I)B_{j}^{\prime })\nonumber\\ &&+\rho m_{5j} C_{i,\tau }^{\prime } B_{j}], \nonumber \\ \mathcal {C}_{p2}^{b} &=&\frac{1}{16}\left(\frac{1}{2n_{\mu }-\tau n_{i}}- \frac{1}{\rho n_{j}}\right)\nonumber \\ &&\times [ D_{i,\tau }(6\sin I\ B_{j}+(\tau m_{5i}-2\cos I)B_{j}^{\prime })\nonumber\\ &&+\,\rho m_{5j} D_{i,\tau }^{\prime } B_{j}], \nonumber \\ \mathcal {C}_{p0}^{c} &=&-\!\frac{1}{2}\frac{1}{n_{\mu }-\tau n_{i}}\sin I C_{i,\tau } C_{j,\rho }, \nonumber \\ \mathcal {C}_{p1}^{c} &=&\frac{1}{4}\left(\frac{1}{2n_{\mu }-\tau n_{i}}+ \frac{1}{n_{\mu }-\rho n_{j}}\right) \sin I D_{i,\tau } C_{j,\rho }, \nonumber \\ \mathcal {C}_{p2}^{c1} &=&\frac{1}{4}\frac{1}{n_{\mu }-\tau n_{i}}[ C_{i,\tau }(6\sin I\ C_{j,\rho }+(\tau m_{5i}-\cos I)C_{j,\rho }^{\prime }) \nonumber \\ &&+(\rho m_{5j}-\cos I) C_{i,\tau }^{\prime } C_{j,\rho }]\nonumber\\ &&-\frac{1}{8}\frac{1}{2n_{\mu }-\tau n_{i}}\sin I D_{i,\tau } D_{j,\rho }, \nonumber \\ \mathcal {C}_{p2}^{c2} &=&\frac{1}{4}\frac{1}{n_{\mu }-\tau n_{i}}[ (\tau m_{5i}-\cos I)C_{i,\tau } C_{j,\rho }^{\prime }\nonumber\\ &&-\,(\rho m_{5j}-\cos I) C_{i,\tau }^{\prime } C_{j,\rho }]. \end{eqnarray} \tag{ B3 }$$

Equations (51) and (60) depend on the functions S^β_α that are associated to the ecliptic term, to the first-order integration, and to the second-order integration. In particular, the effect of the ecliptic (Equation (11)) is given by

$$\begin{eqnarray} S_{E}^{L} &=&e_{1}-(e_{2}\cos \lambda +e_{3}\sin \lambda)\frac{\cos I}{\sin I}, \nonumber \\ S_{E}^{O} &=&e_{3}\cos \lambda -e_{2}\sin \lambda, \nonumber \\ S_{E}^{M} &=&(e_{2}\cos \lambda +e_{3}\sin \lambda)\frac{1}{\sin I}. \end{eqnarray} \tag{ C1 }$$

Corresponding to first-order terms due to $\mathcal {W}_{1},$ we have

$$\begin{eqnarray} S_{1}^{L} &=&-\!B_{0,p}^{\prime }, \nonumber \\ S_{1}^{M} &=&\cos I B_{0,p}^{\prime }. \end{eqnarray} \tag{ C2 }$$

Finally, the second-order integration gives rise to

$$\begin{eqnarray} S_{2a}^{L} &=&\frac{1}{4}\frac{m_{5i}}{n_{i}}\left[ B_{i}\left(B_{j}^{\prime \prime }-\frac{\cos I}{\sin I}B_{j}^{\prime }\right) +B_{i}^{\prime } B_{j}^{\prime }\right], \nonumber \\ S_{2b}^{L} &=&\frac{1}{2}\sin I\frac{C_{i,\tau }^{\prime }\ C_{j,\rho }+C_{i,\tau }\ C_{j,\rho }^{\prime }}{n_{\mu }-\tau n_{i}}, \nonumber \\ S_{2a}^{M} &=&\frac{1}{4}\frac{m_{5i}}{n_{i}}\left[ B_{i}\left(\frac{ B_{j}^{\prime }}{\sin I}-\cos I B_{j}^{\prime \prime }\right) -\cos I B_{i}^{\prime } B_{j}^{\prime }\right], \nonumber \\ S_{2b}^{M} &=& -\!\frac{1}{2}\sin I\nonumber\\ &\times &\frac{\cos I (C_{i,\tau }^{\prime }\ C_{j,\rho }+C_{i,\tau }\ C_{j,\rho }^{\prime })-\sin I\ C_{i,\tau }\ C_{j,\rho }}{ n_{\mu }-\tau n_{i}}. \end{eqnarray} \tag{ C3 }$$

To numerically evaluate these formulae, it is necessary to recall that these functions depend on the starred canonical variables, whose values are given in Table 10. This observation is also extended to Appendices D and E.

Table 10. Numerical Parameters Used in this Work (Taken from Souchay et al. 1999 and References Therein)

Parameter	Value
I*	−0.4090928041 rad
$\omega _E\,(\simeq \dot{\Phi })$	230121.67526278 rad cy−1
Hd	.0032737548
kM	7546.717329 arcsec cy−1
kS	3475.413512 arcsec cy−1
e1	0 arcsec cy−1
e2	5.341 arcsec cy−1
e3	46.82 arcsec cy−1
$\lambda ^{\ast }_{\rm {J2000}}$	0 rad

Note. From these values, we can obtain the numerical value of n*_μ, since n*_μ = M*/A ≃ ω_EC/A = ω_E(1 − H_d)⁻¹, where the dynamical ellipticity of the Earth is given by H_d = (C − A)/C.

Download table as:  ASCIITypeset image

Equations (53), (54), and (61) contain different functions, their expressions being

• 1.  
  Ecliptic effect (in-phase)

  $$\begin{eqnarray} \mathcal {L}_{E}^{\rm in} &=&-\!\frac{m_{5i}}{n_{i}^{2}}\left[ e_{1}B_{i}^{\prime }+ \frac{(e_{2}\cos \lambda +e_{3}\sin \lambda)}{\sin I}\right.\nonumber\\ &&\left.\times\left(\frac{B_{i}}{ \sin I}-\cos I B_{i}^{\prime }\right) \right], \nonumber \\ \mathcal {O}_{E}^{\rm in} &=&-\!\frac{1}{n_{i}^{2}}\Bigg[ m_{5i}^{2}\ e_{1}B_{i}+(e_{2}\cos \lambda +e_{3}\sin \lambda)\nonumber\\ &&\times\left(B_{i}^{\prime }-m_{5i}^{2}\frac{\cos I }{\sin I}B_{i}\right) \Bigg], \nonumber \\ \mathcal {M}_{E}^{\rm in} &=&\cos I\frac{m_{5i}}{n_{i}^{2}}\left[ e_{1}B_{i}^{\prime }+\frac{(e_{2}\cos \lambda +e_{3}\sin \lambda)}{\sin I} \right.\nonumber\\ &&\left.\times\left(\frac{B_{i}}{\sin I}-\cos I B_{i}^{\prime }\right) \right]. \end{eqnarray} \tag{ D1 }$$
• 2.  
  Ecliptic effect (out-of-phase)

  $$\begin{eqnarray} \mathcal {L}_{E}^{\rm out} &=&\frac{1}{n_{i}^{2}}(e_{2}\sin \lambda -e_{3}\cos \lambda)B_{i}^{\prime \prime }, \nonumber \\ \mathcal {O}_{E}^{\rm out} &=&\frac{m_{5i}}{n_{i}^{2}}(e_{2}\sin \lambda -e_{3}\cos \lambda)\left(\frac{\cos I }{\sin I}B_{i}-B_{i}^{\prime }\right), \nonumber \\ \mathcal {M}_{E}^{\rm out} &=&-\cos I \frac{1}{n_{i}^{2}}(e_{2}\sin \lambda -e_{3}\cos \lambda)B_{i}^{\prime \prime }. \end{eqnarray} \tag{ D2 }$$
• 3.  
  First-order effect ($\mathcal {W}_{1}$)

  $$\begin{equation} \mathcal {L}_{1}=-\frac{B_{i}^{\prime }}{n_{i}},\,\quad \mathcal {O}_{1}=-m_{5i}\frac{B_{i}}{n_{i}},\, \quad \mathcal {M}_{1}=\cos I\frac{B_{i}^{\prime }}{n_{i}}. \end{equation} \tag{ D3 }$$
• 4.  
  Second-order effect ($\mathcal {W}_{2s}$)

  $$\begin{eqnarray} \mathcal {L}_{2}^{s} &=&\frac{m_{5i}}{n_{i}^{2}}\left[ B_{i}\left(B_{0}^{\prime \prime }-\frac{\cos I }{\sin I}B_{0}^{\prime }\right) +B_{i}^{\prime } B_{0}^{\prime }\right], \nonumber \\ \mathcal {O}_{2}^{s} &=&\frac{m_{5i}^{2}}{n_{i}^{2}}B_{i}\ B_{0}^{\prime }, \nonumber \\ \mathcal {M}_{2}^{s} &=&\frac{m_{5i}}{n_{i}^{2}}\left[ B_{i}\left(\frac{ B_{0}^{\prime }}{\sin I}-\cos I B_{0}^{\prime \prime }\right) -\cos I B_{i}^{\prime } B_{0}^{\prime }\right].\nonumber\\ \end{eqnarray} \tag{ D4 }$$
• 5.  
  Second-order effect ($\mathcal {W}_{2p}$)

  $$\begin{eqnarray} \mathcal {L}_{2}^{p1} &=&\frac{1}{8}\frac{1}{\tau n_{i}-\rho n_{j}}\left(\frac{1}{\tau n_{i}}+\frac{1}{\rho n_{j}}\right)\nonumber\\ &&\times \left[ \tau m_{5i}B_{i}\left(B_{j}^{\prime \prime }-\frac{\cos I }{\sin I}B_{j}^{\prime }\right) +\rho m_{5j}B_{i}^{\prime } B_{j}^{\prime }\right], \nonumber \\ \mathcal {O}_{2}^{p1} &=&\frac{1}{8}\frac{\tau m_{5i}}{\tau n_{i}-\rho n_{j}} \left(\frac{1}{\tau n_{i}}+\frac{1}{\rho n_{j}}\right)\nonumber\\ &&\times [ \tau m_{5i }B_{i}\ B_{j}^{\prime }+\rho m_{5j}B_{i}^{\prime } B_{j}^{\prime } ], \nonumber \\ \mathcal {M}_{2}^{p1} &=&\frac{1}{8}\frac{\tau m_{5i}}{\tau n_{i}-\rho n_{j}} \left(\frac{1}{\tau n_{i}}+\frac{1}{\rho n_{j}}\right)\nonumber\\ &&\times \left[ B_{i}\left(\frac{B_{j}^{\prime }}{\sin I}-\cos I B_{j}^{\prime \prime }\right) -\cos I B_{i}^{\prime } B_{j}^{\prime }\right], \nonumber \\ \mathcal {L}_{2}^{p2} &=&\frac{\sin I}{2}\frac{C_{i,\tau }^{\prime }\ C_{j,\rho }+C_{i,\tau }\ C_{j,\rho }^{\prime }}{(\tau n_{i}-\rho n_{j})(n_{\mu }-\tau n_{i})}, \nonumber \\ \mathcal {O}_{2}^{p2} &=&\frac{\sin I}{2}\frac{(\tau m_{5i}-\rho m_{5j})\ C_{i,\tau }\ C_{j,\rho }}{(\tau n_{i}-\rho n_{j})(n_{\mu }-\tau n_{i})}, \nonumber \\ \mathcal {M}_{2}^{p2} &=&\!\frac{\sin I}{2}\frac{\sin I\ C_{i,\tau }\ C_{j,\rho }-\cos I(C_{i,\tau }^{\prime }\ C_{j,\rho }+C_{i,\tau }\ C_{j,\rho }^{\prime })}{(\tau n_{i}-\rho n_{j})(n_{\mu }-\tau n_{i})}.\nonumber\\ \end{eqnarray} \tag{ D5 }$$
• 6.  
  Second-order effect ($\left\lbrace \left\lbrace f;\mathcal {W}_{1}\right\rbrace ; \mathcal {W}_{1}\right\rbrace$)

  $$\begin{eqnarray} \mathcal {L}_{3}^{1} &=&\frac{1}{8}\frac{1}{\tau n_{i}\ \rho n_{j}}\Bigg[ \tau m_{5i}\ B_{i}^{\prime } B_{j}^{\prime }\nonumber\\ && +\rho m_{5j}B_{j}\left(B_{i}^{\prime \prime }-\frac{\cos I }{\sin I}B_{i}^{\prime }\right) \Bigg], \nonumber \\ \mathcal {O}_{3}^{1} &=&\frac{1}{8}\frac{\tau m_{5i}}{\tau n_{i}\ \rho n_{j}} \left[ \tau m_{5i}\ B_{i} B_{j}^{\prime }+\rho m_{5j}B_{i}^{\prime }B_{j} \right], \nonumber \\ \mathcal {M}_{3}^{1} &=&-\!\frac{1}{8}\frac{1}{\tau n_{i}\ \rho n_{j}}\Bigg[ \tau m_{5i}\ \cos I\ B_{i}^{\prime } B_{j}^{\prime }\nonumber\\ && +\rho m_{5j}B_{j}\left(\cos I\ B_{i}^{\prime \prime }-\frac{B_{i}^{\prime } }{\sin I}\right) \Bigg], \nonumber \\ \mathcal {L}_{3}^{2} &=&\frac{1}{2}\sin I\frac{C_{i,\tau }\ C_{j,\rho }^{\prime }}{(n_{\mu }-\tau n_{i})(n_{\mu }-\rho n_{j})}, \nonumber \\ \mathcal {O}_{3}^{2} &=&\frac{\sin I}{2}\frac{(\cos I-\tau m_{5i})\ C_{i,\tau }\ C_{j,\rho }}{(n_{\mu }-\tau n_{i})(n_{\mu }-\rho n_{j})}, \nonumber \\ \mathcal {M}_{3}^{2} &=&\frac{1}{2}\sin I\frac{C_{i,\tau }\big(\frac{1}{2} \sin I\ \ C_{j,\rho }-\cos I\ C_{j,\rho }^{\prime }\big) }{(n_{\mu }-\tau n_{i})(n_{\mu }-\rho n_{j})}.\nonumber\\ \end{eqnarray} \tag{ D6 }$$

These terms are given by Equations (57), (58), and (63), where the different functions appearing there are

• 1.  
  Ecliptic effect (in-phase and out-of-phase)

  $$\begin{eqnarray} \mathcal {N}_{E}^{\rm in} &=&-\!\frac{C_{i,\tau }}{(n_{\mu }-\tau n_{i})^{2}}\Bigg[ \tau m_{5i}\ e_{1}+(1-\tau m_{5i}\cos I)\nonumber\\ &&\times\frac{(e_{2}\cos \lambda +e_{3}\sin \lambda)}{\sin I}\Bigg], \nonumber \\ \mathcal {N}_{E}^{\rm out} &=&\frac{C_{i,\tau }^{\prime }}{(n_{\mu }-\tau n_{i})^{2}}(e_{2}\sin \lambda -e_{3}\cos \lambda). \end{eqnarray} \tag{ E1 }$$
• 2.  
  First-order effect ($\mathcal {W}_{1}$)

  $$\begin{equation} \mathcal {N}_{1}=\frac{C_{i,\tau }}{(n_{\mu }-\tau n_{i})}. \end{equation} \tag{ E2 }$$
• 3.  
  Second-order effect ($\mathcal {W}_{2s}$ and $\mathcal {W}_{2p}$)

  $$\begin{eqnarray} \mathcal {N}_{2}^{s} &=&\frac{C_{i,\tau }}{(n_{\mu }-\tau n_{i})^{2}}[ (\tau m_{5i}-\cos I)B_{0}^{\prime }+3\sin I\ B_{0}], \nonumber \\ \mathcal {N}_{2}^{p1} &=&\frac{1}{4}\left(\frac{1}{n_{\mu }-\tau n_{i}}- \frac{1}{\rho n_{j}}\right) \frac{1}{n_{\mu }-\tau n_{i}+\rho n_{j}} \nonumber \\ &&\times [ C_{i,\tau }[ (\tau m_{5i}-\cos I)B_{j}^{\prime }+3\sin I\ B_{j}]\nonumber\\ && +\rho m_{5j}C_{i,\tau \ }^{\prime }B_{j}], \nonumber \\ \mathcal {N}_{2}^{p2} &=&-\!\frac{\sin I}{4}\left(\frac{1}{2n_{\mu }-\tau n_{i} }+\frac{1}{n_{\mu }-\rho n_{j}}\right)\nonumber\\ &&\times \frac{\ D_{i,\tau }\ C_{j,\rho }}{ n_{\mu }-\tau n_{i}+\rho n_{j}}. \end{eqnarray} \tag{ E3 }$$
• 4.  
  Second-order effect ($\left\lbrace \left\lbrace f;\mathcal {W}_{1}\right\rbrace ; \mathcal {W}_{1}\right\rbrace$)

  $$\begin{eqnarray} \mathcal {N}_{3}^{1L} &=&\frac{1}{4}\frac{1}{(n_{\mu }-\tau n_{i})\ \rho n_{j} }\Bigg[ \rho m_{5j}\ B_{j}\Bigg(2\frac{\cos I }{\sin I}C_{i,\tau }-C_{i,\tau }^{\prime }\Bigg)\nonumber \\ && +3\sin I\ C_{i,\tau } B_{j}-(\tau m_{5i}+\cos I)C_{i,\tau } B_{j}^{\prime }\Bigg], \nonumber \\ \mathcal {N}_{3}^{1O} &=&\frac{1}{4}\frac{1}{(n_{\mu }-\tau n_{i})\ \rho n_{j} }[ -\rho m_{5j}\ C_{i,\tau }^{\prime }\ B_{j} \nonumber \\ && +\,3\sin I\ C_{i,\tau } B_{j}-(\tau m_{5i}+\cos I)C_{i,\tau } B_{j}^{\prime }], \nonumber \\ \mathcal {N}_{3}^{1M} &=&\frac{1}{4}\frac{1}{(n_{\mu }-\tau n_{i})\ \rho n_{j} }\Bigg[ \rho m_{5j}\ B_{j}\nonumber\\ &&\times\Bigg(2\frac{C_{i,\tau } }{\sin I\cos I} -C_{i,\tau }^{\prime }\Bigg) \nonumber \\ && +3\sin I\ C_{i,\tau } B_{j}-(\tau m_{5i}+\cos I)C_{i,\tau } B_{j}^{\prime }\Bigg], \nonumber \\ \mathcal {N}_{3}^{1P} &=&-\!\frac{1}{4}\frac{1}{(n_{\mu }-\tau n_{i})\ \rho n_{j}}[ \rho m_{5j}C_{i,\tau }^{\prime }\ B_{j} \nonumber \\ &&-\,3\sin I\ C_{i,\tau } B_{j} +(\tau m_{5i}-\cos I)C_{i,\tau } B_{j}^{\prime }], \nonumber \\ \mathcal {N}_{3}^{2} &=&-\!\frac{\sin I}{4}\frac{D_{i,\tau }\ C_{j,\rho }}{ (2n_{\mu }-\tau n_{i})\ (n_{\mu }-\rho n_{j})}, \nonumber \\ \mathcal {N}_{3}^{3L} &=&-\!\frac{1}{2}\frac{C_{i,\tau }\ C_{j,-\rho }}{(n_{\mu }-\tau n_{i})\ (n_{\mu }+\rho n_{j})}, \nonumber \\ \mathcal {N}_{3}^{3O} &=&\frac{1}{4}\frac{C_{i,\tau }\ C_{j,-\rho }}{(n_{\mu }-\tau n_{i})\ (n_{\mu }+\rho n_{j})}\nonumber\\ && -\frac{1}{4}\frac{C_{i,\tau }\ C_{j,\rho }}{(n_{\mu }-\tau n_{i}) (n_{\mu }-\rho n_{j})}. \end{eqnarray} \tag{ E4 }$$