Issue 
A&A
Volume 604, August 2017



Article Number  A92  
Number of page(s)  13  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201730490  
Published online  14 August 2017 
Dynamical adjustments in IAU 2000A nutation series arising from IAU 2006 precession
^{1} Dept. of Aerospace Engineering, University of León, 24071 León, Spain
^{2} Dept. of Applied Mathematics, University of Alicante, PO Box 99, 03080 Alicante, Spain
email: Alberto.Escapa@ua.es
^{3} Dept. of Applied Mathematics, University of Valladolid, 47011 Valladolid, Spain
Received: 20 January 2017
Accepted: 23 May 2017
The adoption of International Astronomical Union (IAU) 2006 precession model, IAU 2006 precession, requires IAU 2000A nutation to be adjusted to ensure compatibility between both theories. This consists of adding small terms to some nutation amplitudes relevant at the microarcsecond level. Those contributions were derived in previously published articles and are incorporated into current astronomical standards. They are due to the estimation process of nutation amplitudes by Very Long Baseline Interferometry (VLBI) and to the changes induced by the J_{2} rate present in the precession theory. We focus on the second kind of those adjustments, and develop a simple model of the Earth nutation capable of determining all the changes arising in the theoretical construction of the nutation series in a dynamical consistent way. This entails the consideration of three main classes of effects: the J_{2} rate, the orbital coefficients rate, and the variations induced by the update of some IAU 2006 precession quantities. With this aim, we construct a first order model for the nutations of the angular momentum axis of the nonrigid Earth. Our treatment is based on a Hamiltonian formalism and leads to analytical formulae for the nutation amplitudes in the form of inphase, outofphase, and mixed secular terms. They allow numerical evaluation of the contributions of the former effects. We conclude that the accepted corrections associated with the J_{2} rate must be supplemented with new, hitherto unconsidered terms of the same order of magnitude, and that these should be incorporated into present standards.
Key words: astrometry / ephemerides / reference systems / methods: analytical / celestial mechanics
© ESO, 2017
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The present International Astronomical Union (IAU) model of the Earth nutationprecession is comprised of two parts. The first part describes the Earth nutation and is based on the work by Mathews et al. (2002) and the second considers the precessional motion and was developed in Capitaine et al. (2003, 2005). They are commonly referred to as IAU 2000A nutation, for its more precise model, and IAU 2006 precession. They were adopted in Resolution B1.6 and Resolution B1 of the XXIVth and XXVIth IAU general assemblies, respectively, held in Manchester (2000) and Prague (2006).
Although the forced rotation of the Celestial Intermediate Pole (CIP) is kinematically a single motion, the approach followed by IAU, which consists of separating precession from nutation, is very convenient. This is due to the high complexity involved in the right modeling of the Earth’s rotation. Indeed, former IAU theories, such as the pairs IAU 1976 precession (Lieske 1977) and IAU 1980 nutation (Seidelmann 1980) or IAU 1976 precession and IAU 2000A nutation^{1}, followed a similar scheme.
When pursuing utmost accuracies, however, this twofold approach needs to be complemented with the inclusion of some amendments in the nutation or precession component, an aspect not considered yet in IAU resolutions. By doing so, we can guarantee that both theories remain consistent and compatible, a central aim in the current development of the theories of the Earth rotation as recognized in the last IAU working groups devoted to this issue (e.g., Ferrándiz & Gross 2016).
In fact, accuracy requirements at the microarcsecond level (μas) demand the adoption of such adjustments. This is the case of the small changes considered in IAU 2000A nutation series as a consequence of adopting IAU 2006 precession (e.g., Capitaine et al. 2005, and Capitaine & Wallace 2006). They are incorporated into the most important standards, such as IERS Conventions (2010), the Explanatory Supplement to the Astronomical Almanac (2013), and Standards of Fundamental Astronomy (SOFA) routines (e.g., Hohenkerk 2012)^{2}.
The origin of these small adjustments is twofold (Capitaine et al. 2005). Some of them stem from the nutation amplitude estimation process related to Very Long Baseline Interferometry (VLBI) technique, with the same root as that discussed for the VLBI estimate for the precession rate on the obliquity of the ecliptic (Capitaine et al. 2004); this only affects nutation in longitude, being equal to its IAU 2000A counterpart multiplied by a common constant for all the amplitudes. That constant depends on the IAU 1976 and IAU 2006 precession obliquity values, which are different. The relevant contributions at the μas level provide inphase terms.
Others are due to the J_{2} rate introduced in IAU2006 precession, but not considered in the IAU 2000A nutation model. They entail additional corrections both in longitude and obliquity. Since nutation amplitudes are roughly proportional to J_{2}, and in view of the smallness of the ratio of the J_{2} rate to J_{2}, they were accounted for in a simple way by multiplying IAU 2000A nutation series by that ratio and time (Capitaine et al. 2005). In this way, they gave rise to additional mixed secular terms^{3}.
Here, we aim at revisiting the contributions to the IAU 2000A nutation stemming from IAU 2006 precession, completing and extending the work initiated by Escapa et al. (2014). We focus on the changes arising from the theoretical construction of the nutation series, not considering the aforementioned contributions related to VLBI.
Our approach, however, is different from that followed by Capitaine et al. (2005) or Capitaine & Wallace (2006). Specifically, we develop a consistent dynamical model from which we obtain all the changes entailed by the adoption of IAU 2006 precession by reconstructing the analytical solutions of the nutational motion.
Those changes comprise three kind of effects. The first one arises from the inclusion of the J_{2} rate in the Earth model, equivalent to an Earth dynamical ellipticity H_{d} rate. In addition to mixed secular terms, it will be derived that this rate also originates outofphase contributions.
For consistency, it is also necessary to include in our discussion other time rates present in Earth rotation models, similar to that of J_{2}. This leads to the full consideration of the orbital coefficients rate, which is currently only taken into account in the determination of ordinary mixed secular terms. We show that this gives raise to outofphase nutations, whose magnitude must be considered at the μas level.
Finally, the variation on the nutations caused by the changes of some IAU 2006 precession quantities with respect to the IAU 1976 ones are analyzed; these latter ones being employed when evaluating IAU 2000A nutation amplitudes. This accounts for new inphase and mixed secular term contributions.
The effects to be considered are small. Therefore, their impact in the nutations are also expected to be small, typically of some tens of μas and tens of μas per century (e.g., Capitaine et al. 2005; Escapa et al. 2014). This fact makes it possible to consider a simplified, but dynamically coherent, treatment.
Specifically, we restrict our research to the first order modeling of the nutations of the Earth angular momentum axis, that is, to the Poisson terms (Kinoshita 1977). At this perturbation order, these terms are independent of the internal structure of the Earth (Moritz & Mueller 1987).
Although there are several possible ways to tackle this problem, we follow the Hamiltonian approach. There are two main reasons for this. One is on the basis of the rigid model employed in IAU 2000A nutation (Souchay el al. 1999). Following this approach guarantees that the derived contributions can be properly viewed as a means of ensuring compatibility between IAU 2006 precession and IAU 2000A nutation models, as was already done in Capitaine et al. (2005).
On the other one, the Hamiltonian formalism can be extended to consider other nonrigid Earth models (e.g., Getino & Ferrándiz 1995; Getino & Ferrándiz 2001; Escapa et al. 2001, Escapa 2011). Hence, the developments made here can be directly incorporated to those theories, providing an integrated framework to study the Earth’s rotational motion.
The developed analytical model for the Poisson terms has allowed us to comprehensively discuss the influence of IAU 2006 precession in the IAU 2000A nutation series. Explicitly, we conclude that current corrections associated to the J_{2} rate (Capitaine et al. 2005) must be supplemented with new, hitherto unconsidered terms of the same order of magnitude.
The structure of the paper is as follows. In Sect. 2 we formulate the Hamiltonian that describes the nutations of the angular momentum axis of a nonrigid Earth model with a secular variation of its dynamical ellipticity, due to a J_{2} rate of nontidal origin (Williams 1994). Those nutations are induced by the gravitational action of the Moon and the Sun and in their modeling the time rate of their associated orbital coefficients has been taken into account.
The resulting equations of motion are solved, at the first order, in Sect. 3 with the aid of a canonical perturbation method (Hori 1966). These expressions provide analytical formulae for the nutation amplitudes. They are transformed to the form of inphase, outofphase, and mixed secular terms, following the common practice of numerical standards (e.g., IERS Conventions 2010).
The nutation terms are numerically evaluated in Sect. 4, where we display the amplitudes greater than 0.01μ as arising from the H_{d} rate, the orbital coefficients rate, and the change in the values of some precession quantities with respect to IAU 1976 precession ones. Besides, we make a comparison with the contributions found in Capitaine et al. (2005) and considered in IERS Conventions (2010).
Finally, in Sect. 5 we draw some conclusions about the robustness of our Earth model and provide the final value of the adjustments greater that 1 μas to IAU 2000A nutation series. Their inclusion in current standards would guarantee the compatibility of IAU 2000A nutation and IAU 2006 precession. The paper is completed with an appendix where we update, by means of a simple model, the numerical value of the orbital coefficients rate.
2. Hamiltonian of the nonrigid Earth model with secular varying dynamical ellipticity
2.1. Rotation of the angular momentum axis
Although the fundamentals of the Hamiltonian theory are well known, it is necessary to present some details of the developments. In doing that, we emphasize the dependence of the nutations on the precession quantities and incorporate some contributions that are usually discarded but are needed to ensure the consistency of the model. Further explanations can be found, for example, in Kinoshita (1977); Souchay et al. (1999); Efroimsky & Escapa (2007); and Getino et al. (2010) for rigid models and in Getino & Ferrándiz (2001); Escapa et al. (2001); and Escapa (2011) for nonrigid ones.
The rotational motion of the Earth around its center of mass O is characterized by relating a quasiinertial reference system OXYZ with a reference system Oxyz attached to the Earth in some prescribed way; for example by Tisserand mean axes, Oz =e_{z} being the symmetry axis, also named the Earth figure axis.
That rotational motion can be described by a timedependent rotation matrix R that is usually parameterized by the 3–1–3 Euler sequence $\mathit{R}\mathrm{=}\mathit{R}\left(\mathit{\phi ,\epsilon ,\psi}\right)\mathrm{=}{\mathit{R}}_{\mathrm{3}}\left(\mathit{\phi}\right){\mathit{R}}_{\mathrm{1}}\mathrm{(}{\mathit{\epsilon}}^{\mathrm{)}}{\mathit{R}}_{\mathrm{3}}\left(\mathit{\psi}\right)\mathit{,}$(1)where R_{1,3}are the rotation matrices with respect to the first and the third axes. From a dynamical point of view, the problem is solved when the dependence of the Euler angles on time is given explicitly.
Kinoshita’s theory (Kinoshita 1977) mathematically tackles the rotational motion not by means of Euler angles and their time derivatives, but through Andoyer canonical variables given by the canonical pairs (M,μ); ${\mathrm{(}}^{}\mathit{N}\mathrm{,}{\mathit{\nu}}^{\mathrm{)}}$; and ${\mathrm{(}}^{}\mathrm{\Lambda}\mathrm{,}{\mathit{\lambda}}^{\mathrm{)}}$. The canonical momenta M, N, Λ are related to the rotational angular momentum of the Earth L by $\mathit{M}\mathrm{=}\left{L}\right\mathrm{=}\mathit{L,}\mathit{N}\mathrm{=}{L}\mathrm{\xb7}{{e}}_{\mathit{z}}\mathit{,}\mathrm{\Lambda}\mathrm{=}{L}\mathrm{\xb7}{{e}}_{\mathit{Z}}\mathit{,}$(2)where e_{Z} denotes the third vector of the system OXYZ. By introducing the angles that L makes with the vectors e_{Z} and e_{z}, denoted as I and σ, respectively, in Eqs. (2), we can also write $\mathrm{\Lambda}\mathrm{=}\mathit{M}\mathrm{cos}\mathit{I}\mathrm{,}\mathit{N}\mathrm{=}\mathit{M}\mathrm{cos}\mathit{\sigma}\mathit{.}$(3)In terms of the Andoyer variables, the rotation matrix R is described through $\mathit{R}\mathrm{=}{\mathit{R}}_{\mathrm{3}}\mathrm{(}{\mathit{\nu}}^{\mathrm{)}}{\mathit{R}}_{\mathrm{1}}\mathrm{(}{\mathit{\sigma}}^{\mathrm{)}}{\mathit{R}}_{\mathrm{3}}\left(\mathit{\mu}\right){\mathit{R}}_{\mathrm{1}}\mathrm{(}{\mathit{I}}^{\mathrm{)}}{\mathit{R}}_{\mathrm{3}}\mathrm{(}{\mathit{\lambda}}^{\mathrm{)}}\mathit{.}$(4)Combining Eq. (1) with Eq. (4) we can relate Euler angles to Andoyer variables. Since in the Earth case the angle σ is about 10^{6} rad, it is possible to obtain simple polynomial expansions in σ for those expressions. In particular, to obtain an accurate representation of the motion of the figure axis it is necessary to reach up to the first order in σ (Kinoshita 1977) or even second order (Getino et al. 2010). These orders are mainly related to Oppolzer terms.
However, since the contributions discussed in this study are small, we can maintain a zero order expansion in σ. By doing so, the angles fixing the position of the figure axis, that is, the longitude and the obliquity ^{4} read as $\mathit{\psi}\mathrm{=}\mathit{\lambda}\mathrm{+}\mathit{O}\mathrm{\left(}\mathit{\sigma}\mathrm{\right)}\mathit{,\epsilon}\mathrm{=}\mathit{I}\mathrm{+}\mathit{O}\mathrm{\left(}\mathit{\sigma}\mathrm{\right)}\mathit{,}$(5)that is to say, its evolution is the same as that of the angular momentum axis e_{L} = LL^{1}. This temporal evolution arises from three of the differential equations of motion, or Hamiltonian equations, $\frac{\mathrm{d}\mathit{\lambda}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{\partial}\mathrm{\mathscr{H}}}{\mathit{\partial}\mathrm{\Lambda}}\mathit{,}\frac{\mathrm{d\Lambda}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{}\frac{\mathit{\partial}\mathrm{\mathscr{H}}}{\mathit{\partial \lambda}}\mathit{,}\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{}\frac{\mathit{\partial}\mathrm{\mathscr{H}}}{\mathit{\partial \mu}}\mathit{,}$(6)with the proper initial conditions at t = t_{0}. Alternatively, we can write $\frac{\mathrm{d}\mathit{\lambda}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{}\frac{\mathrm{1}}{\mathit{M}\mathrm{sin}\mathit{I}}\frac{\mathit{\partial}\mathrm{\mathscr{H}}}{\mathit{\partial I}}\mathit{,}\frac{\mathrm{d}\mathit{I}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{1}}{\mathit{M}\mathrm{sin}\mathit{I}}\frac{\mathit{\partial}\mathrm{\mathscr{H}}}{\mathit{\partial \lambda}}\mathrm{}\frac{\mathrm{cot}\mathit{I}}{\mathit{M}}\frac{\mathit{\partial}\mathrm{\mathscr{H}}}{\mathit{\partial \mu}}\mathrm{\xb7}$(7)In these expressions, ℋ is the Hamiltonian function of the system that depends on the canonical variables and on time.
The above reasonings can also be applied to nonrigid models of one, two, and three layers. For some of these models it is necessary to extend the number of canonical variables by defining Andoyerlike sets. However, the subset M, N, Λ, μ, ν, and λ is still present in all of them and keeps their relationship with the total angular momentum of the Earth L. Hence, Eqs. (5) and (7) also hold for nonrigid models.
2.2. Reference to the ecliptic of date
Before formulating the Hamiltonian of the Earth model, it is necessary to point out that Kinoshitalike theories do not refer Andoyerlike variables, and hence Euler angles, to the quasiinerial reference system OXYZ, but to a noninertial reference system OX′Y′Z′. That quasiinertial system is denominated as ecliptic of epoch, and the noninertial one as ecliptic of date also named as moving plane (Kinoshita 1977). This circumstance was recognized by Kinoshita himself (Kinoshita 1977) as an advantage of his theory, since by doing that, the development of the disturbing function, due to the gravitational interaction of Moon and Sun, is greatly simplified.
The motion of the ecliptic of date is specified by the angles π and Π. The angle π provides the inclination of the ecliptic of date on the ecliptic of epoch and the angle Π the longitude of the ascending node of the ecliptic of date reckoned from the origin of longitudes of the ecliptic of epoch, that is, the mean equinox of epoch.
The secular variations of these angles, denoted with a subscript A, are assumed to be known functions of time, that is to say, their secular motion with respect to the ecliptic of epoch is determined. This time dependence is provided by some precession theory and given in the form $\begin{array}{ccc}{\mathit{P}}_{\mathit{A}}& \mathrm{=}& \mathrm{sin}{\mathit{\pi}}_{\mathit{A}}\mathrm{cos}{\mathrm{\Pi}}_{\mathit{A}}\mathrm{=}{\mathit{s}}_{\mathrm{1}}\mathit{t}\mathrm{+}{\mathit{s}}_{\mathrm{2}}{\mathit{t}}^{\mathrm{2}}\mathrm{+}{\mathit{s}}_{\mathrm{3}}{\mathit{t}}^{\mathrm{3}}\mathrm{+}\mathit{...}\mathit{,}\\ {\mathit{Q}}_{\mathit{A}}& \mathrm{=}& \end{array}$(8)Throughout our discussions, t denotes the time expressed in Julian centuries of 36 525 days (in practice TT) from some epoch; J2000.0 in our case. After IAU 2006 General Assembly Resolution B1 (e.g., Hilton et al. 2006), this motion was denominated as precession of the ecliptic and the recommended values of the coefficients c_{i} and s_{i} to be used were those arising from the model developed in Capitaine et al. (2003). They are provided for the J2000.0 ecliptic of epoch and its equinox.
When referring the rotational motion to the ecliptic of date, the geometrical meaning of Euler angles, Andoyer variables, their associated angles, and their relationships are analogous to the former ones but the roles of the systems OXYZ and OX′Y′Z′ change. The angles ψ and λ, however, are still reckoned from the OX positive axis^{5}. In this way, the secular part of λ can be identified with the general precession in longitude.
2.3. Hamiltonian of the nonrigid Earth model
2.3.1. Kinetic and potential energies
The Hamiltonian function of our nonrigid Earth model is given by $\mathrm{\mathscr{H}}\mathrm{=}\mathrm{\mathcal{T}}\mathrm{+}\mathrm{\mathcal{V}}\mathrm{+}\mathrm{\mathcal{E}}\mathit{,}$(9)where is the rotational kinetic energy of the model, its potential energy, and ℰ is a complementary term arising from the rotation of the reference system OX′Y′Z′ with respect to the OXYZ one.
The expression of the kinetic energy does depend on the Earth model and provides the torquefree motion, which is the leading term in Eq. (9) because the Earth is a fast rotator. However, in that free motion, the conservation of the angular momentum of the Earth implies that does not depend on the canonical variables μ, Λ, and λ. Therefore, will not provide any contribution to Eqs. (7). Hence, as far we are just concerned with obtaining first order Poisson terms of the rotation of the Earth, it is not necessary to work out the explicit form of .
The function in Eq. (9) is the gravitational potential energy of the system. Here, this potential energy is due to the gravitational interaction caused by the Moon and the Sun, and also the planets, on the nonspherical Earth. It disturbs the rotation of the Earth, inducing small forced deviations with respect to the torquefree motion situation.
Because the effects investigated here are small, we can simply keep the second degree term in the multipolar expansion of the potential due to the Moon and the Sun. This term provides the main contribution to the perturbations of the free rotational motion and has the form $\mathrm{\mathcal{V}}\mathrm{=}\frac{{\mathit{\kappa}}^{\mathrm{2}}{\mathit{m}}^{\mathrm{\prime}}}{{\mathit{r}}^{\mathrm{3}}}\frac{\mathrm{2}\mathit{C}\mathrm{}\mathit{A}\mathrm{}\mathit{B}}{\mathrm{2}}{\mathit{P}}_{\mathrm{2}}\left(\mathrm{sin}\mathit{\delta}\right)\mathit{,}$(10)where κ^{2} is the universal constant of gravitation; m′ is the mass of the perturbing body; r its distance to O, and δ its latitude with respect to the system Oxyz.
The principal moments of inertia A, B, C of the Earth along the axes Ox, Oy, Oz are assumed to have longtime changes induced by the J_{2} rate of nontidal origin, mainly from postglacial rebound due to the nonrigidity of mantle, with ${\mathit{J}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}\mathit{C}\mathrm{}\mathrm{(}\mathit{A}\mathrm{+}{\mathit{B}}^{\mathrm{)}}}{\mathrm{2}\mathit{m}{\mathit{R}}^{\mathrm{2}}}\mathit{,}$(11)where m is the mass of the Earth and R its equatorial radius. The influence of this J_{2} rate on Earth’s precession was first considered in Williams ( 1994) and then incorporated to IAU 2006 precession.
Longtime variations in the moments of inertia stemming from J_{2} rate are given by (Burša et al. 2008) $\frac{\mathrm{d}\mathit{C}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathrm{d}{\mathit{J}}_{\mathrm{2}}}{\mathrm{d}\mathit{t}}\mathit{m}{\mathit{R}}^{\mathrm{2}}\mathit{,}\frac{\mathrm{d}\mathit{A}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{d}\mathit{B}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{}\frac{\mathrm{1}}{\mathrm{3}}\frac{\mathrm{d}{\mathit{J}}_{\mathrm{2}}}{\mathrm{d}\mathit{t}}\mathit{m}{\mathit{R}}^{\mathrm{2}}\mathit{.}$(12)Thus, the rate of the dynamical ellipticity of the Earth, ${\mathit{H}}_{\mathrm{d}}\mathrm{=}\frac{\mathrm{2}\mathit{C}\mathrm{}\mathrm{(}\mathit{A}\mathrm{+}{\mathit{B}}^{\mathrm{)}}}{\mathrm{2}\mathit{C}}\mathit{,}$(13)can be also computed (Burša et al. 2008).
The J_{2} secular evolution as a function of time is slow, hence we express it simply as a first order polynomial on time. In particular, with reference to some epoch t_{0} we consider expansions of the form ${\mathit{J}}_{\mathrm{2}}\mathrm{=}{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}\mathrm{+}\mathit{t}{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{1}}\mathrm{,}$(14)where ${\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}\mathrm{=}{\mathit{J}}_{\mathrm{2}}{}^{\mathrm{(}}\mathit{t}_{\mathrm{0}}{}^{\mathrm{)}}$ and ${\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{1}}\mathrm{=}\mathrm{d}{\mathit{J}}_{\mathrm{2}}{}^{\mathrm{(}}\mathit{t}_{\mathrm{0}}{}^{\mathrm{)}}\mathit{/}\mathrm{d}\mathit{t}$, also denoted as . According to the values given in Williams (1994) the ratio J_{2,1}/J_{2,0} is −2.771 × 10^{6}cy^{1}, so we simply keep firstorder terms of this parameter.
In a similar way, we make use of the expansions $\begin{array}{ccc}{\mathit{k}}^{\mathrm{\prime}}& \mathrm{=}& \mathrm{3}\frac{{\mathit{\kappa}}^{\mathrm{2}}{\mathit{m}}^{\mathrm{\prime}}}{{\mathit{a}}^{\mathrm{3}}}\frac{\mathrm{2}\mathit{C}\mathrm{}\mathit{A}\mathrm{}\mathit{B}}{\mathrm{2}}\mathrm{=}{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}\mathrm{+}\mathit{t}{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}\mathit{,}\\ {\mathit{H}}_{\mathrm{d}}& \mathrm{=}& \end{array}$(15)with a the semimajor axis derived through the Keplerian equation from the corresponding mean motion of the perturber. Equations (11), (12), and (13) combined with expansions (14) and (15) give rise to $\frac{{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{1}}}{{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}}\mathrm{=}\frac{{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{1}}}{{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}}\mathrm{}\frac{\mathrm{1}}{{\mathit{C}}_{\mathrm{0}}}\frac{\mathit{dC}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{1}}}{{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}}\left(\mathrm{1}\mathrm{}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}\right)\mathit{.}$(16)Since 2H_{d,0}/ 3 is about 10^{3} (Williams 1994), within our order of approximation we can neglect this term when compared with 1 and take $\frac{{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{1}}}{{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}}\mathrm{=}\frac{{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{1}}}{{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}\mathit{.}$(17)As a consequence of the above considerations, Eq. (10) can be cast into the form $\mathrm{\mathcal{V}}\mathrm{=}{\mathit{k}}^{\mathrm{\prime}}{\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{3}}{\mathit{P}}_{\mathrm{2}}\left(\mathrm{sin}\mathit{\delta}\right)\mathit{.}$(18)To formulate the differential equations (Eq. (7)), it is necessary to determine the Andoyer variables’ dependence on Eq. (18). This procedure can be accomplished following Kinoshita (1977). In this way, zero order in σ gives $\mathrm{\mathcal{V}}\mathrm{=}{\mathit{k}}^{\mathrm{\prime}}\sum _{\mathit{i}}{\mathit{B}}_{\mathit{i}}\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}\mathit{.}$(19)This contribution to the gravitational potential energy must be considered for each perturber, that is, the Moon and the Sun. However, to lighten the notation we do not make any explicit reference to this fact unless there is some risk of confusion.
In Eq. (19) the index i denotes a quintuplet of integers $\mathit{i}\mathrm{=}{}^{\mathrm{(}}\mathit{m}_{\mathit{i}\mathrm{1}}\mathit{,}{\mathit{m}}_{\mathit{i}\mathrm{2}}\mathit{,}{\mathit{m}}_{\mathit{i}\mathrm{3}}\mathit{,}{\mathit{m}}_{\mathit{i}\mathrm{4}}\mathit{,}{\mathit{m}}_{\mathit{i}\mathrm{5}}{}^{\mathrm{)}}$and the argument Θ_{i} is given by $\begin{array}{ccc}{\mathrm{\Theta}}_{\mathit{i}}& \mathrm{=}& {\mathit{m}}_{\mathit{i}\mathrm{1}}{\mathit{l}}_{\mathrm{M}}\mathrm{+}{\mathit{m}}_{\mathit{i}\mathrm{2}}{\mathit{l}}_{\mathrm{S}}\mathrm{+}{\mathit{m}}_{\mathit{i}\mathrm{3}}\mathit{F}\mathrm{+}{\mathit{m}}_{\mathit{i}\mathrm{4}}\mathit{D}\mathrm{+}{\mathit{m}}_{\mathit{i}\mathrm{5}}\left(\mathrm{\Omega \u0305}\mathrm{}\mathit{\lambda}\right)\\ & \mathrm{=}& \end{array}$(20)where l_{M}, l_{S}, F, D are related to the Delaunay variables of the Moon and the Sun and is the Moon mean longitude of the node referring to the mean equinox of epoch^{6} (e.g., Simon et al. 1994). For our purposes, can be considered as an affine function of time of the form $\mathrm{\Theta \u0305}\mathit{i}\mathrm{=}\mathit{n\u0305}\mathit{i}\mathit{t}\mathrm{+}\mathrm{\Theta \u0305}\mathit{i}\mathrm{0}\mathit{.}$(21)The functions B_{i}, or in short, the orbital functions, are defined (Kinoshita 1977) as ${\mathit{B}}_{\mathit{i}}\mathrm{=}\mathrm{}\frac{\mathrm{1}}{\mathrm{6}}\mathrm{(}\mathrm{3}{\mathrm{cos}}^{\mathrm{2}}\mathit{I}\mathrm{}{\mathrm{1}}^{\mathrm{)}}\mathit{A}\begin{array}{c}\left(\mathrm{0}\right)\\ \mathit{i}\end{array}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\mathit{A}\begin{array}{c}\left(\mathrm{1}\right)\\ \mathit{i}\end{array}\mathrm{sin}\mathrm{2}\mathit{I}\mathrm{}\frac{\mathrm{1}}{\mathrm{4}}\mathit{A}\begin{array}{c}\left(\mathrm{2}\right)\\ \mathit{i}\end{array}{\mathrm{sin}}^{\mathrm{2}}\mathit{I,}$(22)with $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i}\end{array}\mathrm{=}\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{0}\end{array}\mathrm{+}\mathit{t}\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}\mathit{.}$(23)The rates $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$, which render $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i}\end{array}$ timedependent, are due to the secular change of the value of the Sun’s eccentricity (Kinoshita 1977). Although they are also present in the Moon orbital expansions, due to an indirect effect, their main contribution is related to the Sun terms. They induce a time rate on B_{i}, so we must write ${\mathit{B}}_{\mathit{i}}\mathrm{=}{\mathit{B}}_{\mathit{i,}\mathrm{0}}\mathrm{+}\mathit{t}{\mathit{B}}_{\mathit{i,}\mathrm{1}}\mathit{.}$(24)Typically, the largest value of the ratio B_{i,1}/B_{i,0} is of the order 10^{3}cy^{1}, thus their influence on the nutations was only partially considered (e.g., Kinoshita 1977). To obtain a consistent treatment of the Earth model considered, however, it is necessary to take into account the above longterm evolution of the orbital functions due to that of the orbital coefficients. This is clear if we examine the time dependence of k′ (Eq. (15)).
A list of the arguments i and the numerical values of coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{0}\end{array}$ and $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ can be found in Appendix A. There, we also provide the time expansions of l_{M}, l_{S}, F, D, and .
2.3.2. Complementary term
The complementary term ℰ obeys the expression (Efroimsky & Escapa 2007) $\mathrm{\mathcal{E}}\mathrm{=}\mathrm{}{L}\mathrm{\xb7}\varkappa \mathit{,}$(25)where ϰ is the angular velocity vector of the reference system OX′Y′Z′ with respect to the OXYZ one. Writing both vectors in the reference system OX′Y′Z′ gives $\mathrm{\mathcal{E}}\mathrm{=}\mathit{M}\mathrm{sin}\mathit{I}\mathrm{(}{\varkappa}_{\mathrm{2}}\mathrm{cos}\mathit{\lambda}\mathrm{}{\varkappa}_{\mathrm{1}}\mathrm{sin}\mathit{\lambda}\mathrm{)}\mathrm{+}\mathrm{\Lambda}{\varkappa}_{\mathrm{3}}\mathit{.}$(26)The angular velocity ϰ depends on the functions of time π_{A} and Π_{A} through $\begin{array}{c}\\ {\varkappa}_{\mathrm{1}}& \mathrm{=}\mathrm{cos}{\mathrm{\Pi}}_{\mathit{A}}\frac{\mathrm{d}{\mathit{\pi}}_{\mathit{A}}}{\mathrm{d}\mathit{t}}\mathrm{}\mathrm{sin}{\mathit{\pi}}_{\mathit{A}}\mathrm{sin}{\mathrm{\Pi}}_{\mathit{A}}\frac{\mathrm{d}{\mathrm{\Pi}}_{\mathit{A}}}{\mathrm{d}\mathit{t}}\mathit{,}\\ \\ {\varkappa}_{\mathrm{2}}& \mathrm{=}\mathrm{sin}{\mathrm{\Pi}}_{\mathit{A}}\frac{\mathrm{d}{\mathit{\pi}}_{\mathit{A}}}{\mathrm{d}\mathit{t}}\mathrm{+}\mathrm{sin}{\mathit{\pi}}_{\mathit{A}}\mathrm{cos}{\mathrm{\Pi}}_{\mathit{A}}\frac{\mathrm{d}{\mathrm{\Pi}}_{\mathit{A}}}{\mathrm{d}\mathit{t}}\mathit{,}\\ \\ \mathrm{and}\u2001{\varkappa}_{\mathrm{3}}& \mathrm{=}\mathrm{(}\mathrm{1}\mathrm{}\mathrm{cos}{\mathit{\pi}}_{\mathit{A}}\mathrm{)}\frac{\mathrm{d}{\mathrm{\Pi}}_{\mathit{A}}}{\mathrm{d}\mathit{t}}\mathrm{\xb7}\end{array}$(27)The consideration of Eqs. (8) in the above relationships leads to a polynomial expansion in t of ϰ_{i}. To the first order in t we have ${\varkappa}_{\mathrm{1}}\mathrm{=}{\mathit{s}}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{s}}_{\mathrm{2}}\mathit{t,}{\varkappa}_{\mathrm{2}}\mathrm{=}{\mathit{c}}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{c}}_{\mathrm{2}}\mathit{t,}{\varkappa}_{\mathrm{3}}\mathrm{=}\mathrm{0.}$(28)
3. Analytical formulae for the nutations in longitude and obliquity
3.1. First order solution
Once constructed, the Hamiltonian of the rotation of our nonrigid Earth model, we have to solve Eq. (7). Nevertheless, the integration of those differential equations in closed form is not possible. Since the action of the external bodies on the rotation of the free Earth can be viewed properly as a disturbance, perturbation theories are the best candidates to accomplish that integration, yielding analytical approximated solutions.
In the case of Kinoshitalike theories, the perturbation method is implemented via Hori’s algorithm (Hori 1966), based on Lie canonical transformations and an averaging method. We sketch the main features of the method in the context of our research. A detailed review of its fundamentals can be found in FerrazMello (2007) and of its application in Kinoshita (1977).
The relative magnitude of the different parts of the Hamiltonian function allow it to be separated into an unperturbed part ℋ_{0}, related to the torquefree motion, plus a perturbation ℋ_{1} characterized by the coefficients k′ and ϰ_{i} appearing in Eqs. (19) and (26). That is to say, we can write^{7}${\mathrm{\mathscr{H}}\mathrm{=}\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{+}{\mathrm{\mathscr{H}}}_{\mathrm{1}}\mathrm{=}\mathrm{\mathcal{T}}\mathrm{+}\left(\mathrm{\mathcal{V}}\mathrm{+}\mathrm{\mathcal{E}}\right)\mathit{.}$(29)First order nutations are computed from a function that transforms the original canonical set into a new one, denoted with an asterisk, that is easier to integrate. This generating or determining function is given by ${\mathrm{\mathcal{W}}}^{\mathrm{\ast}}\mathrm{=}{\mathrm{\int}}_{\mathrm{UP}}{\mathrm{\mathscr{H}}}_{\mathrm{1}\mathrm{per}}\mathrm{d}\mathit{t,}$(30)where ℋ_{1per} is the periodic part of the perturbation. In our case, it arises from the periodic part of the gravitational potential energy (Kinoshita 1977), since the complementary term ℰ does not have shortperiod terms. Therefore, we have ${\mathrm{\mathscr{H}}}_{\mathrm{1}\mathrm{per}}\mathrm{=}{\mathrm{\mathcal{V}}}_{\mathrm{per}}\mathrm{=}{\mathit{k}}^{\mathrm{\prime}}\sum _{\mathit{i}\ne \mathrm{0}}{\mathit{B}}_{\mathit{i}}\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}\mathrm{=}{\mathrm{\mathscr{H}}}_{\mathrm{1}\mathrm{per}}\mathrm{(}\mathit{I,}{\mathit{\lambda}}^{\mathrm{)}}\mathit{.}$(31)The integral in Eq. (30) is evaluated along the trajectories of the unperturbed problem (UP). Here, this problem is given by the torquefree motion. Since from Eq. (31), ℋ_{1per} is a function of I and λ, only the evolution of these variables is needed. As we point out in Sect. 2.3.1, in the torquefree motion they are constant. Once performed the integral, I and λ have to be substituted by their counterparts I^{∗} and λ^{∗} of the new canonical set.
Following this procedure we can explicitly calculate the generating function by means of Eq. (30). That integration is direct and allows us to obtain without performing approximations or neglecting any contribution. However, in view of the magnitude of the ratios ${\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}\mathit{/}{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}\mathrm{=}{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{1}}\mathit{/}{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}$ and B_{i,1}/B_{i,0}, we simply keep terms up to the first order in these parameters ${\mathrm{\mathcal{V}}}_{\mathrm{per}}\mathrm{=}\sum _{\mathit{i}\ne \mathrm{0}}{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}{\mathit{B}}_{\mathit{i,}\mathrm{0}}\left[\mathrm{1}\mathrm{+}\left(\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}\mathrm{+}\frac{{\mathit{B}}_{\mathit{i,}\mathrm{1}}}{{\mathit{B}}_{\mathit{i,}\mathrm{0}}}\right)\mathit{t}\right]\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}\mathrm{+}\mathit{O}\mathrm{\left(}{\mathit{t}}^{\mathrm{2}}\mathrm{\right)}\mathit{.}$(32)Therefore, considering the constancy of I and λ in the unperturbed problem, and the time dependence of Θ_{i} and ${\mathit{B}}_{\mathit{i}}{}^{\mathrm{\right(}}\mathit{I}^{\mathrm{\left)}}$ (Eqs. (20), (21), and (24)) the integrands reduce to integrals of the form $\begin{array}{ccc}\mathrm{\int}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{a}{\mathit{t}}^{\mathrm{)}}\mathrm{cos}\left(\mathit{\alpha t}\mathrm{+}\mathit{\beta}\right)& \mathrm{=}& \frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{a}{\mathit{t}}^{\mathrm{)}}}{\mathit{\alpha}}\mathrm{sin}\left(\mathit{\alpha t}\mathrm{+}\mathit{\beta}\right)\\ & & \end{array}$(33)where a, α ≠ 0, and β are real constants. By doing so, the generating function can be written as ${\mathrm{\mathcal{W}}}^{\mathrm{\ast}}{\mathrm{=}\mathrm{\mathcal{W}}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{+}\mathit{t}{\mathrm{\mathcal{W}}}_{\mathrm{1}}^{\mathrm{\ast}}\mathit{,}$(34)with $\begin{array}{ccc}{\mathrm{\mathcal{W}}}_{\mathrm{0}}^{\mathrm{\ast}}& \mathrm{=}& \sum _{\mathit{i}\ne \mathrm{0}}\frac{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{n\u0305}\mathit{i}}\mathrm{sin}{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}\mathrm{+}\sum _{\mathit{i}\ne \mathrm{0}}\frac{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{n\u0305}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}}\left(\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}\mathrm{+}\frac{{\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}}{{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}\right)\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}\mathit{,}\\ {\mathrm{\mathcal{W}}}_{\mathrm{1}}^{\mathrm{\ast}}& \mathrm{=}& \end{array}$(35)The first term in ${\mathrm{\mathcal{W}}}_{\mathrm{0}}^{\mathrm{\ast}}$ was computed in Kinoshita (1977) and coincides with his expression, except for , due to the different choice of the unperturbed problem^{8}. It is also the case for the term in ${\mathrm{\mathcal{W}}}_{\mathrm{1}}^{\mathrm{\ast}}$ proportional to ${\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}\mathit{/}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}$. The parts linear in ${\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}\mathit{/}{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}$ appear in the generating function as a consequence of incorporating a J_{2} rate into the Earth model. They were not considered in Kinoshita (1977) nor Souchay et al. (1999), since they arise from the nonrigidity.
The parameters ${\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}\mathit{/}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}$ are also present in ${\mathrm{\mathcal{W}}}_{\mathrm{0}}^{\mathrm{\ast}}$ as outofphase terms. Since their contribution to the periodic part of the nutations is very small, those outofphase nutations coming from the orbital coefficients rate were not considered in Kinoshita (1977), although it was explicitly stated that they must be taken into account in an exact integration of the generating function.
3.2. Nutations in longitude and obliquity
At the first order, the periodic variation of any canonical function f of the original variables can be found through the relationship (Hori 1966) $\mathrm{\Delta}\mathit{f}\mathrm{=}\left\{{\mathit{f}}^{\mathrm{\ast}}\mathit{,}{\mathrm{\mathcal{W}}}^{\mathrm{\ast}}\right\}\mathit{,}$(36)where f^{∗} is equal to f but expressed in the transformed set by literal substitution and {− , −} stands for the Poisson bracket of two functions of the canonical variables. To obtain Δf as an explicit function of time we have to substitute, after computing the Poisson bracket, the transformed variables in terms of time.
In principle, that would require solving the canonical equations generated by the Hamiltonian , free of shortperiod terms, resulting from applying Hori’s perturbation method with generating function .
For our purposes those needed variables reduce to λ^{∗}, I^{∗}, and M^{∗}. Since their evolution is secular or involves longperiod terms, they are usually expressed as ${\mathit{\lambda}}^{\mathrm{\ast}}\mathrm{=}{\mathit{\lambda}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{+}\mathit{t}\hspace{0.17em}{\mathit{n}}_{{\mathit{\lambda}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{+}\mathit{...}\mathit{,}{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{+}\mathit{t}\hspace{0.17em}{\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{+}\mathit{...}\mathit{,}{\mathit{M}}^{\mathrm{\ast}}\mathrm{=}{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{+}\mathit{t}{\mathit{M}}_{\mathrm{1}}^{\mathrm{\ast}}\mathrm{+}\mathit{...}$(37)where ${\mathit{\lambda}}_{\mathrm{0}}^{\mathrm{\ast}}$, ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$, and ${\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}$ denote the initial conditions at J2000.0. The function λ^{∗} is related to the general precession in longitude and I^{∗} with the precession in obliquity. The numerical values of ${\mathit{\lambda}}_{\mathrm{0}}^{\mathrm{\ast}}$, ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$, ${\mathit{n}}_{{\mathit{\lambda}}_{\mathrm{0}}^{\mathrm{\ast}}}$, and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$ are borrowed from some precessional model (e.g., Lieske et al. 1977; Fukushima 2003; Capitaine et al. 2003).
The new canonical variable M^{∗} can be approximated by ${\mathit{M}}^{\mathrm{\ast}}\mathrm{\simeq}\mathit{C}{\mathit{\omega}}_{\mathrm{E}}^{\mathrm{\ast}}$, ω_{E} being the angular velocity of the Earth (Kinoshita 1977). As pointed out by Williams (1994), processes leading to a secular variation of J_{2} leave M^{∗} unaffected, that is, ${\mathit{M}}^{\mathrm{\ast}}\mathrm{=}{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{=}{\mathit{C}}_{\mathrm{0}}{\mathit{\omega}}_{\mathrm{E}\mathit{,}\mathrm{0}}^{\mathrm{\ast}}$, although both C and ${\mathit{\omega}}_{\mathrm{E}}^{\mathrm{\ast}}$ depend on time^{9}.
This approach could appear as unsatisfactory, since from a dynamical perspective the complete evolution of the canonical variables should be obtained from the same Hamiltonian function. From a practical point of view, however, this procedure is adequate considering current levels of accuracy.
Therefore, the periodic evolution of the figure axis, which at our level of approximation is the same as that of the angular momentum axis, is computed from $\mathrm{\Delta}\mathit{\psi}\mathrm{=}\mathrm{\Delta}\mathit{\lambda}\mathrm{=}\left\{{\mathit{\lambda}}^{\mathrm{\ast}}\mathit{,}{\mathrm{\mathcal{W}}}^{\mathrm{\ast}}\right\}\mathit{,}\mathrm{\Delta}\mathit{\epsilon}\mathrm{=}\mathrm{\Delta}\mathit{I}\mathrm{=}\left\{{\mathit{I}}^{\mathrm{\ast}}\mathit{,}{\mathrm{\mathcal{W}}}^{\mathrm{\ast}}\right\}\mathit{.}$(38)Considering Eqs. (3), (20), and (35), we have $\mathrm{\Delta}\mathit{\psi}\mathrm{=}\mathrm{}\frac{\mathrm{1}}{{\mathit{M}}^{\mathrm{\ast}}\mathrm{sin}{\mathit{I}}^{\mathrm{\ast}}}\frac{\mathit{\partial}{\mathrm{\mathcal{W}}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\mathit{,}\mathrm{\Delta}\mathit{\epsilon}\mathrm{=}\frac{\mathrm{1}}{{\mathit{M}}^{\mathrm{\ast}}\mathrm{sin}{\mathit{I}}^{\mathrm{\ast}}}\frac{\mathit{\partial}{\mathrm{\mathcal{W}}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{\lambda}}^{\mathrm{\ast}}}\mathrm{\xb7}$(39)Specifically, we get $\begin{array}{ccc}\mathrm{\Delta}\mathit{\psi}& \mathrm{=}& \frac{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}}\sum _{\mathit{i}\ne \mathrm{0}}\mathrm{[}{\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{in}}\mathrm{sin}{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}\mathrm{+}{\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{out}}\mathrm{cos}{{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}}^{\mathrm{]}}\mathit{,}\\ \mathrm{\Delta}\mathit{\epsilon}& \mathrm{=}& \end{array}$(40)with $\begin{array}{ccc}{\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{in}}& \mathrm{=}& \mathrm{}\frac{\mathrm{1}}{\mathit{n\u0305}\mathit{i}\mathrm{sin}{\mathit{I}}^{\mathrm{\ast}}}\left[\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\mathrm{+}\left(\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\mathrm{+}\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)\mathit{t}\right]\mathit{,}\\ {\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{out}}& \mathrm{=}& \mathrm{}\frac{\mathrm{1}}{\mathit{n\u0305}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}\mathrm{sin}{\mathit{I}}^{\mathrm{\ast}}}\left(\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\mathrm{+}\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)\mathit{,}\\ {\mathrm{\mathcal{O}}}_{\mathit{i}}^{\mathrm{in}}& \mathrm{=}& \mathrm{}\frac{{\mathit{m}}_{\mathit{i}\mathrm{5}}}{\mathit{n\u0305}\mathit{i}\mathrm{sin}{\mathit{I}}^{\mathrm{\ast}}}\left[{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}\mathrm{+}\left(\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}\mathrm{+}{\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}\right)\mathit{t}\right]\mathit{,}\\ {\mathrm{\mathcal{O}}}_{\mathit{i}}^{\mathrm{out}}& \mathrm{=}& \end{array}$(41)These analytical formulae clearly reflect the influence of the J_{2} rate and the orbital coefficients rate on the nutation model. These formulae also depend on the precession quantities I^{∗} and λ^{∗} (through the argument ${\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}$, Eq. (20)).
When taking ${\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}\mathrm{=}\mathrm{0}$ in Eqs. (41), the amplitudes ${\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{in}}$ and ${\mathrm{\mathcal{O}}}_{\mathit{i}}^{\mathrm{in}}$ coincide ^{10} with those of Kinoshita (1977) except for , which is substituted by n_{i}. This is due to the different choice of the unperturbed problem followed in Kinoshita (1977). When computing the nutations in Kinoshitas’s way, one must be aware of the dependence of n_{i} on the canonical variables I^{∗} and λ^{∗}, since it is the case of n_{λ∗}. With our approach, these additional terms arise as a secondorder effect coming from the coupling of the complementary term and the secular part of the gravitational potential energy with the generating function (Getino et al. 2010). The outofphase amplitudes ${\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{out}}$ and ${\mathrm{\mathcal{O}}}_{\mathit{i}}^{\mathrm{out}}$ arise as a consequence of the J_{2} and orbital coefficients rate and they are not present in other studies.
Finally, the computation of the nutations as quasipolynomials in t is achieved by substituting I^{∗} as an explicit function of time, in addition to the time dependence of the arguments ${\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}$(Eq. (20)). Within the scope of this study, we consider only the first order expansion ${\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{+}\mathit{t}\hspace{0.17em}{\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$ to evaluate the amplitudes ${\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{in}}$, ${\mathrm{\mathcal{O}}}_{\mathit{i}}^{\mathrm{in}}$, ${\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{out}}$, and ${\mathrm{\mathcal{O}}}_{\mathit{i}}^{\mathrm{out}}$, with the constants ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$ and${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$ provided by IAU 2006 precession.
3.3. Final expressions
The common use of numerical standards (e.g., Kaplan 2005; IERS Conventions 2010; or Explanatory Supplement to the Astronomicals Almanac 2013) employs a firstorder time expansion of the amplitudes provided in Eqs. (40).
If we denote generically any of those amplitudes as ℱ_{i} and take into account Eq. (37), we can write at the first order in t$\begin{array}{ccc}{\mathrm{\mathcal{F}}}_{\mathit{i}}& \mathrm{=}& {\left({\mathrm{\mathcal{F}}}_{\mathit{i}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{+}\mathit{t}\left[{\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}{\left(\frac{\mathit{\partial}{\mathrm{\mathcal{F}}}_{\mathit{i}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{+}{\left(\frac{\mathit{\partial}{\mathrm{\mathcal{F}}}_{\mathit{i}}}{\mathit{\partial t}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\right]\\ & \mathrm{=}& \end{array}$(42)The coefficients ℱ_{i,1} have a double origin: one part is due to the rate of I^{∗}, ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$, and the other comes from the rates of k′ and B_{i}.
With this relationship in mind, and keeping only firstorder terms in the related parameters associated to the rates ${\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}$, B_{i,1} (and its derivatives with respect to I^{∗}), ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$, we can finally write the nutations in longitude and obliquity as $\mathrm{\Delta}\mathit{\psi}\mathrm{=}{\mathrm{\Delta}}_{\mathrm{0}}\mathit{\psi}\mathrm{+}{\mathrm{\Delta}}_{\mathrm{1}}\mathit{\psi ,}\mathrm{\Delta}\mathit{\epsilon}\mathrm{=}{\mathrm{\Delta}}_{\mathrm{0}}\mathit{\epsilon}\mathrm{+}{\mathrm{\Delta}}_{\mathrm{1}}\mathit{\epsilon}\mathit{.}$(43)The terms with subscript 0 are trigonometrical polynomials with constant coefficients of the form $\begin{array}{ccc}{\mathrm{\Delta}}_{\mathrm{0}}\mathit{\psi}& \mathrm{=}& \frac{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}}\sum _{\mathit{i}\ne \mathrm{0}}\mathrm{[}{\mathrm{\mathcal{L}}}_{\mathit{i,}\mathrm{0}}^{\mathrm{in}}\mathrm{sin}{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}\mathrm{+}{\mathrm{\mathcal{L}}}_{\mathit{i,}\mathrm{0}}^{\mathrm{out}}\mathrm{cos}{{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}}^{\mathrm{]}}\mathit{,}\\ {\mathrm{\Delta}}_{\mathrm{0}}\mathit{\epsilon}& \mathrm{=}& \end{array}$(44)where $\begin{array}{ccc}{\mathrm{\mathcal{L}}}_{\mathit{i}\mathit{.0}}^{\mathrm{in}}& \mathrm{=}& \mathrm{}\frac{\mathrm{1}}{\mathit{n\u0305}\mathit{i}\mathrm{sin}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}{\left(\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathit{,}{\mathrm{\mathcal{L}}}_{\mathit{i,}\mathrm{0}}^{\mathrm{out}}\mathrm{=}{\mathrm{(}{{\mathrm{\mathcal{L}}}_{\mathit{i}}^{\mathrm{out}}}^{\mathrm{)}}}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathit{,}\\ {\mathrm{\mathcal{O}}}_{\mathit{i,}\mathrm{0}}^{\mathrm{in}}& \mathrm{=}& \end{array}$(45)The terms with subscript 1 contain t in their amplitudes, that is, they are mixed secular terms^{11}, whose expressions are given by ${\mathrm{\Delta}}_{\mathrm{1}}\mathit{\psi}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}}\sum _{\mathit{i}\ne \mathrm{0}}{\mathrm{\mathcal{L}}}_{\mathit{i,}\mathrm{1}}^{\mathrm{in}}\mathit{t}\mathrm{sin}{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}\mathit{,}{\mathrm{\Delta}}_{\mathrm{1}}\mathit{\epsilon}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}}\sum _{\mathit{i}\ne \mathrm{0}}{\mathrm{\mathcal{O}}}_{\mathit{i,}\mathrm{1}}^{\mathrm{in}}\mathit{t}\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}^{\mathrm{\ast}}\mathit{,}$(46)where $\begin{array}{c}\\ {\mathrm{\mathcal{L}}}_{\mathit{i,}\mathrm{1}}^{\mathrm{in}}\mathrm{=}& \mathrm{}\frac{\mathrm{1}}{\mathit{n\u0305}\mathit{i}\mathrm{sin}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\{\left[\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{\left(\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{+}{\left(\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\right]\\ & \\ & \mathrm{}{\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\left[\mathrm{cot}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}{\left(\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{}{\left(\frac{{\mathit{\partial}}^{\mathrm{2}}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{{\mathrm{\ast}}^{\mathrm{2}}}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\right]\begin{array}{c}\mathrm{\u23ab}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ac}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ad}\end{array}\mathit{,}\\ & \\ {\mathrm{\mathcal{O}}}_{\mathit{i,}\mathrm{1}}^{\mathrm{in}}\mathrm{=}& \mathrm{}\frac{{\mathit{m}}_{\mathit{i}\mathrm{5}}}{\mathit{n\u0305}\mathit{i}\mathrm{sin}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\{\left[\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{\mathrm{(}{{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}^{\mathrm{)}}}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{+}{\mathrm{(}{{\mathit{B}}_{\mathit{i,}\mathrm{1}}^{\mathrm{\ast}}}^{\mathrm{)}}}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\right]\\ & \\ & \mathrm{}{\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\left[\mathrm{cot}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}{\mathrm{(}{{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}^{\mathrm{)}}}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{}{\left(\frac{\mathit{\partial}{\mathit{B}}_{\mathit{i,}\mathrm{0}}^{\mathrm{\ast}}}{\mathit{\partial}{\mathit{I}}^{\mathrm{\ast}}}\right)}_{{\mathit{I}}^{\mathrm{\ast}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\right]\begin{array}{c}\mathrm{\u23ab}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ac}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ad}\end{array}\mathit{.}\end{array}$(47)
4. Discussion
4.1. Adopted parameters
The evaluation of Δ_{0}ψ, Δ_{0}ε, Δ_{1}ψ, and Δ_{1}ε requires the adoption of different parameters. The main ones, affecting all the contributions tackled in this work, are related to the orbital functions B_{i}, the orbital frequency , ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$, and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$.
The orbital functions and their derivatives can be computed from their definition (Eq. (22)), with the i arguments and orbital coefficients listed in Table A.4 of Appendix A. For each argument its unperturbed frequency is obtained considering the time evolution, at the first order in t, of l_{M}, l_{S}, F, D, and given in Table A.1. Finally, the values ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$ and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$ are taken as those adopted in IAU 2006 precession for ϵ_{A} expansion, with the proper change of sign, that is, ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{=}\mathrm{}{\mathrm{84381.4060}}^{\mathrm{\prime \prime}}$ and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{=}{\mathrm{46.836769}}^{\mathrm{\prime \prime}}\mathit{/}\mathrm{cy}$ (Capitaine et al. 2005).
Besides, it is necessary to consider the value of the scaling factor ${\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}\mathrm{}\mathrm{1}}$ for each perturber. This can be written as ${\mathit{k}}_{\mathrm{0}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}}{{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{=}\mathrm{3}\frac{{\mathit{\kappa}}^{\mathrm{2}}{\mathit{m}}^{\mathrm{\prime}}}{{\mathit{a}}^{\mathrm{3}}{\mathit{\omega}}_{\mathrm{E}\mathit{,}\mathrm{0}}^{\mathrm{\ast}}}\frac{\mathrm{2}{\mathit{C}}_{\mathrm{0}}\mathrm{}{\mathit{A}}_{\mathrm{0}}\mathrm{}{\mathit{B}}_{\mathrm{0}}}{\mathrm{2}{\mathit{C}}_{\mathrm{0}}}\mathrm{=}\left(\mathrm{3}\frac{{\mathit{\kappa}}^{\mathrm{2}}{\mathit{m}}^{\mathrm{\prime}}}{{\mathit{a}}^{\mathrm{3}}{\mathit{\omega}}_{\mathrm{E}\mathit{,}\mathrm{0}}^{\mathrm{\ast}}}\right){\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}\mathit{.}$(48)The values of k_{0} for the Moon and the Sun were provided in Souchay et al. (1999) for a rigid Earth model. The nonrigidity of the Earth entails a change of the value of H_{d,0}, so we have to consider ${\mathit{k}}_{\mathrm{0}}\mathrm{=}{\mathit{k}}_{\mathrm{0}}^{\mathit{R}}\left(\frac{{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}}{{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}^{\mathit{R}}}\right)\mathrm{\xb7}$(49)In particular, for the rigid Earth we take the values ${\mathit{k}}_{\mathrm{0}\mathit{,}\mathrm{M}}^{\mathit{R}}\mathrm{=}{\mathrm{7546.717329}}^{\mathrm{\prime \prime}}$/cy, and ${\mathit{k}}_{\mathrm{0}\mathit{,S}}^{\mathit{R}}\mathrm{=}{\mathrm{3475.413512}}^{\mathrm{\prime \prime}}\mathit{/}$cy, for the Moon and the Sun, respectively, and ${\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}^{\mathit{R}}\mathrm{=}\mathrm{0.0032737548}$ (Souchay et al. 1999). As regards the dynamical ellipticity of the nonrigid Earth, we use the value associated with IAU 2000A nutation, that is to say, H_{d,0} = 0.0032737949 with an uncertainty under four parts in 10^{7} (Mathews et al. 2002).
4.2. Contributions from the H_{d} rate and the orbital coefficients rate
The influence of the H_{d} rate appears in Δ_{0}ψ, Δ_{0}ε, Δ_{1}ψ, and Δ_{1}ε through the parameter ${\mathit{k}}_{\mathrm{1}}^{\mathrm{\prime}}\mathit{/}{\mathit{k}}_{\mathrm{0}}^{\mathrm{\prime}}\mathrm{=}{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{1}}\mathit{/}{\mathit{H}}_{\mathrm{d}\mathit{,}\mathrm{0}}\mathrm{=}{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{1}}\mathit{/}{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}$ (Eq. (17)). It induces outofphase and mixed secular terms both in longitude and obliquity.
The numerical value of this parameter was fixed in IAU 2006 precession from J_{2,0} = 1.0826358 × 10^{3} and J_{2,1} = −3.001 × 10^{9}/cy (Capitaine et al. 2005), providing^{12}H_{d,1}/H_{d,0} = −2.7719 × 10^{6}/cy. This value is very close to that computed from Williams’ data (Williams 1994), −2.7710 × 10^{6}/cy. More recently, IAU adopted an updated System of Astronomical Constants in its 2009 General Assembly (e.g., Luzum et al. 2011 ). From this system we also get H_{d,1}/H_{d,0} = −2.7710 × 10^{6}/cy, which is the value that we have employed in our computations.
The numerical influence of the H_{d} rate on the nutations is displayed in Table 1 where, as for the other tables of this section, we have only shown those arguments i having any amplitude equal to or greater than 0.01μas or 0.01μas/cy.
Orbital coefficients rate contributions (in brackets values already included in IAU 2000A nutation series).
The mixed secular terms were previously computed in Capitaine et al. (2005) or Capitaine & Wallace (2006) and are considered in IERS Conventions (2010), the Explanatory Supplement to the Astronomical Almanac (2013), or Standards of Fundamental Astronomy (SOFA) routines (e.g., Hohenkerk 2012). Their values are displayed in parentheses in Cols. 9 and 12 of Table 1, according to Table 5.2f^{13} provided in IERS Conventions (2010, Sect. 5.6.3) that has a cutoff of 0.1μas.
Our calculations show good agreement with those results, the larger differences being below 0.2μas/cy. They are mainly related to the different way in obtaining the changes associated with the H_{d} rate. In particular, in Capitaine et al. (2005) it was assumed that all the nutation amplitudes were proportional to H_{d}. Strictly speaking, this is only true for firstorder contributions of gravitational origin arising from the J_{2} term in the geopotential expansion. However, it is a valid approximation for obtaining the small adjustments associated with H_{d} rate, as we have checked with the quite different and comprehensive approach developed in this research.
In contrast, the outofphase terms have not previously been considered with the exception of some preliminary computations reported in Escapa et al. (2014). Indeed, their dynamical origin is the same as that giving raise to mixed secular terms. Hence, a congruent treatment to determine the influence of H_{d} rate on the nutations should also take them into account.
To be consistent, the inclusion of the H_{d} rate entails the consideration of the orbital coefficients rate, since its influence on the nutations is of the same nature as that discussed previously for H_{d} rate. It is shown in Table 2. The mixed secular terms due to that rate, shown in brackets in Cols. 8 and 10, are included in current standards (e.g., IERS Conventions 2010), since they were implicitly considered in rigid Earth theories (e.g., Kinoshita 1977; or Kinoshita & Souchay 1990).
However, the derivations leading to outofphase terms are given here for the first time (Cols. 7 and 9 of Table 2). These contributions would also be present for the rigid Earth, their origin stemming from the secular evolution of Sun eccentricity and not from any dissipative process. They must be considered to reach the truncation level of 0.1μas established in some rigid Earth models (e.g., Souchay et al. 1999).
4.3. Contributions from the change of precession quantities
As a consequence of the adoption of IAU 2006 precession, the value of some precession quantities changed with respect to those adopted in IAU 1976 precession. The nutation amplitudes are affected by these changes, since the underlying rigid Earth theory of IAU 2000A nutation (Souchay et al. 1999) was based on IAU 1976 precession^{14}.
The relevant precession quantities entering into Δ_{0}ψ, Δ_{0}ε, Δ_{1}ψ, and Δ_{1}ε are ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$ and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$. In principle, we should also consider the dependence of the nutations with the Earth dynamical ellipticity H_{d,0}. However, the value given in IAU 2006 precession, 0.00327379448 (Capitaine et al. 2003), belongs to the uncertainty interval of the dynamical ellipticity given in IAU 2000A nutation that runs from 0.00327379372 to 0.00327379612 (Mathews et al. 2002). So, we do not take into account the variation, or indirect effect (Escapa et al. 2016), on the nutations due to this small change of the dynamical ellipticity^{15}.
To compute the variations of the nutation amplitudes associated with the changes in ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$ and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$, we have calculated the difference among the amplitudes, evaluated with IAU 2006 parameters, and their counterparts obtained from IAU 1976 values ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{=}\mathrm{}{\mathrm{84381.448}}^{\mathrm{\prime \prime}}$and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}\mathrm{=}{\mathrm{46.8150}}^{\mathrm{\prime \prime}}$/cy provided in Lieske et al. (1977). The results are displayed in Table 3.
Precession quantities variation contributions.
The changes of precession quantities provide some inphase and mixed secular terms both in longitude and obliquity. There are no outofphase contributions (Eqs. (45)) greater that the threshold fixed in the construction of our tables (0.01μas). Of the discussed effects in this study, it is the only contribution that originates from a change of inphase components, especially relevant for the 18.6yr term in longitude, as was first pointed out in Escapa et al. (2014).
Those inphase changes can be viewed as a reevaluation of rigid Earth nutation series for a new set of parameters, compatible with IAU 2006 precession, affecting in this way all the amplitudes in longitude and obliquity, according to their different dependence on ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$. Since these amplitudes are common with the nutations of the angular momentum axis, they can be directly incorporated to the nonrigid Earth nutation series.
The adjustments considered in Capitaine et al. (2005) affecting the inphase components in longitude, included in current standards, are not related to the variation of the precession quatities discussed here, as was implicitly considered in Escapa et al. (2014).
They are related to the VLBI technique, which is not sensitive to an ecliptic (Capitaine et al. 2003, 2004). Hence, the quantity estimated from the observation is not directly the longitude but the longitude multiplied by the sinus of the obliquity at a certain epoch. This value depends on the precession theory, therefore it is necessary to make a global rescaling when moving from one precession theory to another as indicated in Capitaine et al. (2005).
5. Conclusion
The adoption of the new IAU 2006 precession theory (Capitaine et al. 2003, 2005) entails the introduction of small changes in IAU 2000A nutation theory (Mathews et al. 2002) to ensure compatibility between both models. In particular, we have focused our attention on the changes induced in the theoretical construction of the nutational model, independently of those arising from the estimation process associated to it (Capitaine et al. 2005).
They are of a double nature. On the one hand, it is necessary to take into account in the formulation the effects of the J_{2} rate or, equivalently, the H_{d} rate, present in IAU 2006 precession. In turn, for consistency, this inclusion motivates the full consideration of the orbital coefficients rate. On the other hand, some precession quantities of IAU 2006 precession differ from those adopted in IAU 1976 precession. These variations modify the Earth nutation model, since IAU 2000A nutation was constructed using IAU 1976 precession quantities.
We found several corrections to be included in the series of IAU 2000A nutation. They affect inphase, outofphase, and mixed secular terms both in longitude and obliquity. Specifically, the adjustments larger than 1μas are given by $\begin{array}{ccc}\left(\mathrm{}\mathrm{d\Delta}\mathit{\psi}\right)& \mathrm{=}& \mathrm{}\mathrm{15.6}\mathrm{sin}\mathrm{\Omega}\mathrm{}\mathrm{1.4}\mathrm{cos}\mathrm{\Omega}\mathrm{}\mathrm{0.5}\mathrm{cos}{\mathit{l}}_{\mathrm{S}}\\ & & \u2001\mathrm{+}\mathrm{39.8}\mathit{t}\mathrm{sin}\mathrm{\Omega}\mathrm{}\mathrm{0.6}\mathit{t}\mathrm{sin}\mathrm{2}\mathrm{\Omega}\\ & & \u2001\mathrm{+}\mathrm{3.5}\mathit{t}\mathrm{sin}\left(\mathrm{2}\mathit{F}\mathrm{}\mathrm{2}\mathit{D}\mathrm{+}\mathrm{2}\mathrm{\Omega}\right)\mathrm{+}\mathrm{0.6}\mathit{t}\mathrm{sin}\left(\mathrm{2}\mathit{F}\mathrm{+}\mathrm{2}\mathrm{\Omega}\right)\mathit{,}\\ \left(\mathrm{}\mathrm{d\Delta}\mathit{\epsilon}\right)& \mathrm{=}& \mathrm{+}\mathrm{0.8}\mathrm{cos}\mathrm{\Omega}\mathrm{}\mathrm{0.8}\mathrm{sin}\mathrm{\Omega}\mathrm{}\mathrm{25.1}\mathit{t}\mathrm{cos}\mathrm{\Omega}\\ & & \end{array}$(50)where we considered the sign corresponding to the conventional use of longitude and obliquity. In these expressions the amplitudes are given in μas and μas/cy, t being the number of Julian centuries from J2000.0.
The former corrections were derived from an analytical simplified model of the rotation of the nonrigid Earth developed within the Hamiltonian framework. It was constructed by considering the firstorder modeling of the nutations of the angular momentum axis, that is, the Poisson terms. In view of the numerical magnitude of the obtained contributions (Tables 1–3), the incorporation of other features into our model does not seem to be relevant at the 1μas level.
For example, we could consider the influence of Earth internal structure or of secondorder terms. Earth internal structure affects mainly the nutations through the Oppolzer terms. For a twolayer Earth model, the largest contribution is about 80 mas (milliarcseconds) for the 18.6yr term in longitude (Getino & Ferrándiz 2001). If we take as proxy the value of H_{d,1}/H_{d,0} = −2.7710 × 10^{6}/cy, they produce adjustments smaller than 3 × 10^{1}μas in absolute value.
Reasoning in a similar way, with our scheme of perturbation, secondorder terms are proportional to k^{′2} and k′ϰ_{i}. The largest contribution associated with k^{′2} is about 1250μas for the 9.3yr term in longitude (e.g., Souchay et al. 1999; or Getino et al. 2010). So, we get a change of the order of 6 × 10^{3}μas. The nutations due to k′ϰ_{i} are outofphase terms of the order of 250μas (Getino et al. 2010). In addition to ${\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}$and ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$, they also depend on the precessional quantities s_{1}, s_{2}, c_{1}, and c_{2} through ϰ_{i} (Eq. (28)), especially on c_{1} (Getino et al. 2010). Their IAU 2006 precession values are different from the IAU 1976 precession ones. The relative change in c_{1} is of the order of 1 × 10^{5}, hence the induced variations would be below 3 × 10^{2}μas.
At present, the compatibility between IAU 2006 precession and IAU 2000A nutation is considered in the main sources of Astronomical standards, such as IERS Conventions (2010), the Explanatory Supplement to the Astronomical Almanac (2013), and Standards of Fundamental Astronomy (SOFA) software (e.g., Hohenkerk 2012). The adjustments taken into account (Capitaine et al. 2005) originate from the different value of the obliquity when estimating IAU 2000A nutation amplitudes and from the introduction of the J_{2} rate into the model.
The first effect is not related to this research, since it has its roots in the particular characteristics of VLBI observations. The second one was also considered here, where we showed that the mixed secular terms derived in Capitaine et al. (2005) must be supplemented with some outofphase terms. Besides, our study has unveiled other contributions relevant at the 1μas level, all of them gathered in Eqs. (50) and related to the analytical dynamical development of the nutational model. They should also be incorporated into current standards.
IAU Resolution B1.6 in 2000 referred to as IAU 2000 PrecessionNutation model. In practise, the precession part was limited to some correction of the precession rates in longitude and obliquity of IAU 1976 precession. The limitations from the dynamical point of view of this approach gave raise to the development of new precession theories (e.g., Fukushima 2003; Capitaine et al. 2003, 2005), leading to the adoption of IAU 2006 precession.
To avoid confusions in the context of Earth rotation studies, we prefer the denomination mixed secular terms instead of Poisson terms. In orbital theories, Poisson terms (e.g., Simon et al. 1994) refer to expressions of the form ${\mathit{t}}^{\mathit{n}}\mathrm{sin}\mathit{\alpha}{}^{\mathrm{\right(}}\mathit{t}^{\mathrm{\left)}}$ or ${\mathit{t}}^{\mathit{n}}\mathrm{cos}\mathit{\alpha}{}^{\mathrm{\right(}}\mathit{t}^{\mathrm{\left)}}$, n being a nonnegative integer. They are related to Poisson series (e.g. Danby et al. 1965). However, in Earth rotation theories the name Poisson terms is commonly employed to describe the nutations of angular momentum axis (e.g., Kinoshita 1977). These nutations are the solution of the Poisson equations (Woolard 1953).
Strictly speaking, Eq. (20) refers to the secular part of λ not to λ itself.
Specifically, in Kinoshita (1977) instead of it is considered n_{i}, both related through
It is not the case for other phenomena like those involving tides with energy dissipation. They entail a secular variation of M^{∗}. However, since ${\mathit{M}}_{\mathrm{1}}^{\mathrm{\ast}}\mathit{/}{\mathit{M}}_{\mathrm{0}}^{\mathrm{\ast}}\mathrm{=}\mathrm{}\mathrm{2.20}\mathrm{\times}{\mathrm{10}}^{8}\mathrm{c}{\mathrm{y}}^{1}$ (Williams 1994) their effects are two order of magnitude smaller than those due to J_{2,1}/J_{2,0}, which themselves are small.
We note that the righthand side of Eqs. (6.14) and (6.15) in Kinoshita (1977) must have the opposite sign.
The mixed secular terms arising from ${\mathit{n}}_{{\mathit{I}}_{\mathrm{0}}^{\mathrm{\ast}}}$ are incorporated into current standards (e.g., IERS Conventions 2010). However, to our knowledge, there is no reference providing analytical formulae similar to those derived in this work.
There is a missprint in Capitaine et al. (2005) concerning J_{2} rate which appears as −0.3001 × 10^{9}/cy instead of −3.001 × 10^{9}/cy, a typo also present in IAU 2006 Resolution B1 (e.g., IERS Conventions 2010, Appendix B). In that work the ratio J_{2,1}/J_{2,0} is assigned with the value −2.7774 × 10^{6}/cy from the rate to the acceleration of precession in longitude (Bourda & Capitaine 2004). These differences lead to variations of about 0.1μas in our computations.
Mathews et al. (2002, par. [57]) transformed the rigid nutation amplitudes to the prograde and retrograde form. In this process, they employed a slightly different value of I_{0} from that of IAU 1976 precession and IAU 2006 precession, that is, −arcsin(0.3977769687). As long as the same value of I_{0} is used to reverse the process, it does not influence our discussion.
Acknowledgments
We are grateful to N. Capitaine for helpful comments which enhanced the manuscript and to the valuable suggestions and advices of an anonymous referee which helped to improve this paper. This work has been partially supported by the Spanish projects I+D+I AYA201022039C0202 and AYA201679775P (AEI/FEDER, UE).
References
 Bourda, G., & Capitaine, N. 2004, A&A, 428, 691 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bretagnon, P., & Francou, G. 1988, A&A, 202, 309 [Google Scholar]
 Brower, D., & Clemence, G. M. 1961, Methods of Celestial Mechanics (Academic Press) [Google Scholar]
 Burša, M., Groten, E., & Šíma, Z. 2008, AJ, 135, 1021 [NASA ADS] [CrossRef] [Google Scholar]
 Capitaine, N., & Wallace, P. T. 2006, A&A, 450, 855 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Capitaine, N., Wallace, P. T., & Chapront, J. 2003, A&A, 412, 567 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Capitaine, N., Wallace, P. T., & Chapront, J. 2004, A&A, 421, 365 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Capitaine, N., Wallace, P. T., & Chapront, J. 2005, A&A, 432, 355 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Danby, J. M. A., Deprit, A., & Rom, A. R. M. 1965, Boeing Scientific Research Laboratories, Mathematical Note No. 432 [Google Scholar]
 Efroimsky, M., & Escapa, A. 2007, Celest. Mech. Dyn. Astron, 98, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Escapa, A. 2011, Celest. Mech. Dyn. Astron, 110, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Escapa, A., Getino, J., & Ferrándiz, J.M. 2001, J. Geophys. Res. 106, 11387 [NASA ADS] [CrossRef] [Google Scholar]
 Escapa, A., Getino, J., Ferrándiz, J. M., & Baenas, T. 2014, in Proc. of the Journées 2013, ed. N. Capitaine, Observatoire de Paris, 148 [Google Scholar]
 Escapa, A., Ferrándiz, J. M., Baenas, T., et al. 2016, Pure Appl. Geophys. 173, 861 [NASA ADS] [CrossRef] [Google Scholar]
 Explanatory Supplement to the Astronomical Almanac 2012, third edition, eds. S. E. Urban, & P. K. Seidelmann (University Science Books), 734 [Google Scholar]
 Ferrándiz, J. M., & Gross, R. 2016, International Association of Geodesy Symposia, 143, 533 [Google Scholar]
 FerrazMello, S. 2007, Canonical Perturbation Theories: Degenerate Systems and Resonance. (New York: Springer) [Google Scholar]
 Fukushima, T. 2003, AJ, 126, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Getino, J., & Ferrándiz, J. M. 1995, Celest. Mech. Dyn. Astron., 61, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Getino, J., & Ferrándiz, J. M. 2001, MNRAS, 322, 785 [NASA ADS] [CrossRef] [Google Scholar]
 Getino, J., Escapa, A., & Miguel, D. 2010, AJ, 139, 1916 [NASA ADS] [CrossRef] [Google Scholar]
 Hilton, J., Capitaine, N., Chapront, J., et al. 2006, Celest. Mech. Dyn. Astron., 94, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Hohenkerk, C. Y. 2012, in Proc. of the Journées 2011, eds. H. Schuh, S. Böhm, T. Nilsson, & N. Capitaine, Vienna University of Technology, 21 [Google Scholar]
 Hori, G. 1966, PASJ, 18, 287 [NASA ADS] [Google Scholar]
 IERS Conventions 2010, IERS Technical Note 36, eds. G. Petit, & B. Luzum, 179 [Google Scholar]
 Kaplan, G. H. 2005, Circular No. 179, US Naval Observatory, 103 [Google Scholar]
 Kinoshita, H. 1977, Celest. Mech. Dyn. Astron., 15, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Kinoshita, H., & Souchay, J. 1990, Celest. Mech. Dyn. Astron., 48, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Lieske, J. H., Lederle, T., Fricke, W., & Morando, B. 1977, A&A, 58, 1 [NASA ADS] [Google Scholar]
 Luzum, B., Capitaine, N., Fienga, A., et al. 2011, Celest. Mech. Dyn. Astron., 110, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Mathews, P. M., Herring, T. A., & Buffett, B. A. 2002, J. Geophys. Res., 107, B4 [Google Scholar]
 Moritz, H., & Mueller, I. 1987, Earth Rotation (New York: Frederic Ungar) [Google Scholar]
 Seidelmann, P. K. 1982, Celest. Mech., 27, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. L., Bretagnon, P., Chapront, J., et al. 1994, A&A, 282, 663 [NASA ADS] [Google Scholar]
 Souchay, J., Losley, B., Kinoshita, H., & Folgueira, M. 1999, A&AS, 135, 111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Williams, J. G. 1994, AJ, 108, 711 [NASA ADS] [CrossRef] [Google Scholar]
 Woolard, E. W. 1953, Astr. Pap. Amer. Ephem. Naut. Almanach XV, I, 165 [Google Scholar]
Appendix A: Orbital coefficients rate
The orbital coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i}\end{array}$ arise from the development of the seconddegree disturbing potential of the Moon and the Sun (e.g., Kinoshita 1977; Kinoshita & Souchay 1990). They are defined from the expressions
$\begin{array}{ccc}& & {\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{3}}{\mathit{P}}_{\mathrm{2}}\left(\mathrm{sin}\mathit{\beta}\right)\mathrm{=}\sum _{\mathit{i}}\mathit{A}\begin{array}{c}\left(\mathrm{0}\right)\\ \mathit{i}\end{array}\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}\mathit{,}\\ & & {\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{3}}{\mathit{P}}_{\mathrm{2}}^{\mathrm{1}}\left(\mathrm{sin}\mathit{\beta}\right)\mathrm{sin}\mathrm{(}\mathrm{2}{\mathit{\alpha}}^{\mathrm{)}}\mathrm{=}\mathrm{3}\sum _{\mathit{i}}\mathit{A}\begin{array}{c}\left(\mathrm{1}\right)\\ \mathit{i}\end{array}\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}\mathit{,}\\ & & \end{array}$(A.1)where α and β denote the longitude and the latitude of the Moon or the Sun referring to the ecliptic of date and ${\mathit{P}}_{\mathit{n}}^{\mathit{m}}\left(\mathrm{cos}\mathit{\beta}\right)$ are the nthdegree, mthorder Legendre associated functions. It is assumed that the orbital motion of both perturbers, that is to say, r, α, and β, are known functions of time, provided by some analytical or numerical orbital ephemeris theory.
When focusing on the Sun coefficients and recalling that its orbital motion is given in the ecliptic of date, that is, β ≃ 0, the following simplifications can be adopted (Kinoshita & Souchay 1990) ${\mathit{P}}_{\mathrm{2}}\left(\mathrm{sin}\mathit{\beta}\right)\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}{\mathit{P}}_{\mathrm{2}}^{\mathrm{1}}\left(\mathrm{sin}\mathit{\beta}\right)\mathrm{=}\mathrm{0}\mathit{,}{\mathit{P}}_{\mathrm{2}}^{\mathrm{2}}\left(\mathrm{sin}\mathit{\beta}\right)\mathrm{=}\mathrm{3.}$(A.2)Hence, from Eqs. (A.1) and (A.2), $\mathit{A}\begin{array}{c}\left(\mathrm{0}\right)\\ \mathit{i}\end{array}$ and $\mathit{A}\begin{array}{c}\left(\mathrm{2}\right)\\ \mathit{i}\end{array}$ are determined through $\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{3}}\mathrm{=}\sum _{\mathit{i}}\mathit{A}\begin{array}{c}\left(\mathrm{0}\right)\\ \mathit{i}\end{array}\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}\mathit{,}\hspace{0.17em}{\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{3}}\mathrm{cos}\mathrm{(}\mathrm{2}{\mathit{\alpha}}^{\mathrm{)}}\mathrm{=}\sum _{\mathit{i}}\mathit{A}\begin{array}{c}\left(\mathrm{2}\right)\\ \mathit{i}\end{array}\mathrm{cos}{\mathrm{\Theta}}_{\mathit{i}}\mathit{,}$(A.3)the orbital coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{1}\right)\\ \mathit{i}\end{array}$ being nil.
The Sun orbital coefficients rate, $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ (Eq. (23)), was computed by Kinoshita (1977). Within the context of our work, it is possible to obtain a simple update of those values through the mean orbital elements of the Earth given in Simon et al. (1994), by employing expansion procedures valid for elliptic motion (e.g., Brower & Clemence 1961, Chaps. 1 and 2). Those mean orbital elements were derived from the orbital theory VSOP87 (Bretagnon & Francou 1988), the same used in IAU2006 precession theory (Capitaine et al. 2003). The x, y, z coordinates of a celestial body with respect to a given reference system can be obtained from its , , coordinates referring to the orbital reference system, with the axis in the direction of the periapsis and parallel to the orbital angular momentum vector. This process is accomplished with the introduction of the orbital elements, subscript o, and the transformation $\left(\begin{array}{c}\\ \mathit{x\u0305}\\ \mathit{y\u0305}\\ \mathit{z\u0305}\end{array}\right)\mathrm{=}{\mathit{R}}_{\mathrm{3}}\mathrm{\left(}{\mathit{\omega}}_{\mathrm{o}}\mathrm{\right)}{\mathit{R}}_{\mathrm{1}}\mathrm{\left(}{\mathit{I}}_{\mathrm{o}}\mathrm{\right)}{\mathit{R}}_{\mathrm{3}}\left({\mathrm{\Omega}}_{\mathrm{o}}\right)\left(\begin{array}{c}\\ \mathit{x}\\ \mathit{y}\\ \mathit{z}\end{array}\right)\mathit{,}$(A.4)R_{i} being the elemental rotation matrices and Ω_{o}, I_{o}, and ω_{o} the longitude of the ascending node, the inclination, and the argument of the periapsis.
In our case the xy plane of the given reference system has to be identified with the ecliptic of date. It follows that I_{o} = β = 0, so and z coordinates play no role with z = 0. Besides, since the line of nodes is not defined, the longitude of the periapsis $\begin{array}{c}{}_{\mathrm{\u02dc}}\\ {\mathit{\omega}}_{\mathrm{o}}\end{array}\mathrm{=}{\mathit{\omega}}_{\mathrm{o}}\mathrm{+}{\mathrm{\Omega}}_{\mathrm{o}}$ must be introduced. Hence, Eq. (A.4) leads to $\left(\begin{array}{c}\\ \mathit{x}\\ \mathit{y}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \mathrm{cos}\begin{array}{c}\mathrm{\u02dc}\\ {\mathit{\omega}}_{\mathrm{o}}\end{array}& \mathrm{}\mathrm{sin}\begin{array}{c}\mathrm{\u02dc}\\ {\mathit{\omega}}_{\mathrm{o}}\end{array}\\ \mathrm{sin}\begin{array}{c}\mathrm{\u02dc}\\ {\mathit{\omega}}_{\mathrm{o}}\end{array}& \mathrm{cos}\begin{array}{c}\mathrm{\u02dc}\\ {\mathit{\omega}}_{\mathrm{o}}\end{array}\end{array}\right)\left(\begin{array}{c}\\ \mathit{x\u0305}\\ \mathit{y\u0305}\end{array}\right)\mathit{.}$(A.5)These relationships allow the left hand sides of Eqs. (A.3 ) to be expressed in terms of the elliptic motion elements, since $\mathit{x}\mathrm{=}\mathit{r}\mathrm{cos}\mathit{\alpha ,}\mathit{y}\mathrm{=}\mathit{r}\mathrm{sin}\mathit{\alpha}\mathrm{=}\mathrm{\Rightarrow}\{\begin{array}{c}\\ \mathit{r}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}\\ \\ \mathrm{cos}\mathrm{(}\mathrm{2}{\mathit{\alpha}}^{\mathrm{)}}\mathrm{=}\frac{{\mathit{x}}^{\mathrm{2}}\mathrm{}{\mathit{y}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{\xb7}\end{array}$(A.6)To obtain those expansions, the evolution of the celestial body must be provided in the orbital plane. It is given by its , coordinates that can be written in terms of the eccentric anomaly, u, the orbit eccentricity, e, and the semimajor axis, a. Namely, we have $\mathit{x\u0305}\mathrm{=}\mathit{a}\mathrm{(}\mathrm{cos}\mathit{u}\mathrm{}{\mathit{e}}^{\mathrm{)}}\mathit{,}\mathit{y\u0305}\mathrm{=}\mathit{a}\sqrt{\mathrm{1}\mathrm{}{\mathit{e}}^{\mathrm{2}}}\mathrm{sin}\mathit{u}\mathit{.}$(A.7)Given that the mean anomaly l_{S} is the quantity entering into Θ_{i} (Eq. (20)), it is necessary to express the eccentric anomaly in terms of l_{S}. This is achieved through Kepler’s equation $\mathit{u}\mathrm{}\mathit{e}\mathrm{sin}\mathit{u}\mathrm{=}{\mathit{l}}_{\mathit{S}}\mathit{,}$(A.8)and by introducing the Bessel functions of the first kind and order s and their derivative, J_{s} and ${\mathit{J}}_{\mathit{s}}^{\mathrm{\prime}}$. Specifically, it can be deduced that $\begin{array}{ccc}\mathrm{cos}\mathit{u}& \mathrm{=}& \mathrm{}\frac{\mathit{e}}{\mathrm{2}}\mathrm{+}\mathrm{2}\sum _{\mathit{s}\mathrm{=}\mathrm{1}}^{\mathrm{+}\mathrm{\infty}}\frac{\mathrm{1}}{\mathit{s}}{\mathit{J}}_{\mathit{s}}^{\mathrm{\prime}}\mathrm{(}\mathit{s}{\mathit{e}}^{\mathrm{)}}\mathrm{cos}\mathrm{(}\mathit{s}{\mathit{l}}_{\mathrm{S}}\mathrm{)}\mathit{,}\\ \mathrm{sin}\mathit{u}& \mathrm{=}& \end{array}$(A.9)For small values of the eccentricity, as in our case, the former equations can be expanded with respect to e and truncated to some definite power. For example, to the third order in e we have $\begin{array}{ccc}\mathrm{cos}\mathit{u}& \mathrm{=}& \mathrm{}\frac{\mathit{e}}{\mathrm{2}}\mathrm{+}\left(\mathrm{1}\mathrm{}\frac{\mathrm{3}{\mathit{e}}^{\mathrm{2}}}{\mathrm{8}}\right)\mathrm{cos}{\mathit{l}}_{\mathrm{S}}\mathrm{+}\left(\frac{\mathit{e}}{\mathrm{2}}\mathrm{}\frac{{\mathit{e}}^{\mathrm{3}}}{\mathrm{3}}\right)\mathrm{cos}\mathrm{2}{\mathit{l}}_{\mathrm{S}}\\ & & \mathrm{+}\frac{\mathrm{3}{\mathit{e}}^{\mathrm{2}}}{\mathrm{8}}\mathrm{cos}\mathrm{3}{\mathit{l}}_{\mathrm{S}}\mathrm{+}\frac{{\mathit{e}}^{\mathrm{3}}}{\mathrm{3}}\mathrm{cos}\mathrm{4}{\mathit{l}}_{\mathrm{S}}\mathit{,}\\ \mathrm{sin}\mathit{u}& \mathrm{=}& \left(\mathrm{1}\mathrm{}\frac{{\mathit{e}}^{\mathrm{2}}}{\mathrm{8}}\right)\mathrm{sin}{\mathit{l}}_{\mathrm{S}}\mathrm{+}\left(\frac{\mathit{e}}{\mathrm{2}}\mathrm{}\frac{{\mathit{e}}^{\mathrm{3}}}{\mathrm{6}}\right)\mathrm{sin}\mathrm{2}{\mathit{l}}_{\mathrm{S}}\\ & & \end{array}$(A.10)By doing so, the combination of Eqs. (A.10), (A.7), (A.5), and (A.6) leads to the construction of the equations defining the orbital coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\right)\\ \mathit{i}\end{array}$ and $\mathit{A}\begin{array}{c}\left(\mathrm{2}\right)\\ \mathit{i}\end{array}\hspace{0.17em}$(Eqs. (A.3)). These equations must be expanded again with respect to the eccentricity in order to transform them into the form given by the righthand sides of Eqs. (A.3).
The resulting expressions of this expansion are trigonometric polynomials of l_{S} and $\begin{array}{c}{\mathrm{\u02dc}}_{}\\ {\mathit{\omega}}_{\mathrm{o}}\end{array}$, whose coefficients are polynomials of the eccentricity. In order to express them in terms of Θ_{i} it is necessary to take into account that, from the meaning of the variables defining Θ_{i} (e.g., Kinoshita 1977; Kinoshita & Souchay 1990), the mean longitude of the Sun, L_{S}, equals F + Ω−D (Eq. (20)). Therefore, the longitude of the periapsis^{16} can be written as $\begin{array}{c}\mathrm{\u02dc}\\ {\mathit{\omega}}_{\mathrm{o}}\end{array}\mathrm{=}{\mathit{L}}_{\mathrm{S}}\mathrm{}{\mathit{l}}_{\mathrm{S}}\mathrm{=}\mathit{F}\mathrm{+}\mathrm{\Omega}\mathrm{}\mathit{D}\mathrm{}{\mathit{l}}_{\mathrm{S}}\mathit{.}$(A.11)In this way, to the third order in the eccentricity, we have recomputed the orbital coefficients rate $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ for the same arguments as those provided in Kinoshita (1977) with the aid of Maple software. This task requires setting the numerical value of the mean eccentricity, which for our purposes can be taken as (Simon et al. 1994) $\mathit{e}\mathrm{=}{\mathit{e}}_{\mathrm{0}}\mathrm{+}\mathit{t}{\mathit{e}}_{\mathrm{1}}\mathrm{+}\mathit{...}\mathrm{=}\mathrm{0.0167086342}\mathrm{}\mathrm{0.00004203654}\mathit{t}\mathrm{+}\mathit{...}$(A.12)We have also taken this work (Simon et al. 1994) as the source providing the time evolution of l_{S} and the remaining lunisolar arguments l_{M}, F, D, Ω, and , necessary for other computations performed in this research, such as the value of . They are listed in Table A.1 to the first order in t.
Time evolution of lunisolar arguments (epoch J2000.0).
Sun orbital coefficients rate $\mathit{A}\begin{array}{c}\left(\mathrm{0}\right)\\ \mathit{i,}\mathrm{1}\end{array}$.
The analytical expressions and numerical values of orbital coefficients rate $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ are provided in Tables A.2 and A.3. Their values show little difference with respect to Kinoshita’s computations (Kinoshita 1977), taking into account that these coefficients are multiplied by 10^{8} rd/cy. Considering that for the Moon the effect of the secular change of Sun eccentricity is even smaller, we can keep for this perturber the original values provided in Kinoshita (1977).
For the dynamical adjustments computed in this research due to the orbital coefficients rate, the only argument providing contributions at the μas level, in particular an outofphase term in longitude, is the one with period 365.26 days in Table A.2. However, at the 0.1 μas level, the threshold established in the rigid theory by Souchay et al. (1999), the arguments of periods 365.22 and 121.75 days in Table A.3 must also be considered.
Sun orbital coefficients rate $\mathit{A}\begin{array}{c}\left(\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$.
To make the recalculation of the different numerical evaluations carried in this work more straightforward, we have gathered the used values of the coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i}\end{array}$ in Table A.4. Since we are tackling with small contributions, we have limited to the thirteen arguments providing the larger nutation amplitudes, in addition to the argument with infinite period. The orbital coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{0}\end{array}$ are taken from Kinoshita (1977) with the update provided in Kinoshita & Souchay (1990). The coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ are those from Tables A.2 and A.3 with the exception of the Moon ones that were borrowed from Kinoshita (1977). The epoch used in Kinoshita ( 1977) was J1900.0 and not J2000.0. This affects only $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{0}\end{array}$ and has no effect on $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ coefficients, as was implicitly assumed in the Kinoshita & Souchay (1990) update.
Main orbital coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i}\end{array}$ for Moon and Sun with $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{0}\end{array}$ in 10^{7} rd and $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ in 10^{8} rd/cy (epoch J2000.0).
All Tables
Orbital coefficients rate contributions (in brackets values already included in IAU 2000A nutation series).
Sun orbital coefficients rate $\mathit{A}\begin{array}{c}\left(\mathrm{0}\right)\\ \mathit{i,}\mathrm{1}\end{array}$.
Sun orbital coefficients rate $\mathit{A}\begin{array}{c}\left(\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$.
Main orbital coefficients $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i}\end{array}$ for Moon and Sun with $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{0}\end{array}$ in 10^{7} rd and $\mathit{A}\begin{array}{c}\left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\right)\\ \mathit{i,}\mathrm{1}\end{array}$ in 10^{8} rd/cy (epoch J2000.0).
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.