© ESO, 2017

1. Introduction

The present International Astronomical Union (IAU) model of the Earth nutation-precession 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¹, 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)².

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 in-phase terms.

Others are due to the J₂ 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₂, and in view of the smallness of the ratio of the J₂ rate to J₂, 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³.

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₂ 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 out-of-phase 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₂. 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 out-of-phase 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 in-phase 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 non-rigid 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₂ 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 non-rigid Earth model with a secular variation of its dynamical ellipticity, due to a J₂ 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 in-phase, out-of-phase, 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 non-rigid 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 non-rigid ones.

The rotational motion of the Earth around its center of mass O is characterized by relating a quasi-inertial 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 time-dependent rotation matrix R that is usually parameterized by the 3–1–3 Euler sequence (1)where R_1,3are 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,μ); ; and . The canonical momenta M, N, Λ are related to the rotational angular momentum of the Earth L by (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 (3)In terms of the Andoyer variables, the rotation matrix R is described through (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 ⁴ read as (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, (6)with the proper initial conditions at t = t₀. Alternatively, we can write (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 non-rigid models of one, two, and three layers. For some of these models it is necessary to extend the number of canonical variables by defining Andoyer-like 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 non-rigid models.

2.2. Reference to the ecliptic of date

Before formulating the Hamiltonian of the Earth model, it is necessary to point out that Kinoshita-like theories do not refer Andoyer-like variables, and hence Euler angles, to the quasi-inerial reference system OXYZ, but to a non-inertial reference system OX′Y′Z′. That quasi-inertial system is denominated as ecliptic of epoch, and the non-inertial 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 (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⁵. In this way, the secular part of λ can be identified with the general precession in longitude.

2.3. Hamiltonian of the non-rigid Earth model

2.3.1. Kinetic and potential energies

The Hamiltonian function of our non-rigid Earth model is given by (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 torque-free 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 non-spherical Earth. It disturbs the rotation of the Earth, inducing small forced deviations with respect to the torque-free 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 (10)where κ² 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 long-time changes induced by the J₂ rate of nontidal origin, mainly from post-glacial rebound due to the non-rigidity of mantle, with (11)where m is the mass of the Earth and R its equatorial radius. The influence of this J₂ rate on Earth’s precession was first considered in Williams ( 1994) and then incorporated to IAU 2006 precession.

Long-time variations in the moments of inertia stemming from J₂ rate are given by (Burša et al. 2008) (12)Thus, the rate of the dynamical ellipticity of the Earth, (13)can be also computed (Burša et al. 2008).

The J₂ 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₀ we consider expansions of the form (14)where and , also denoted as . According to the values given in Williams (1994) the ratio J_2,1/J_2,0 is −2.771 × 10^-6cy^-1, so we simply keep first-order terms of this parameter.

In a similar way, we make use of the expansions (15)with a the semi-major 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 (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 (17)As a consequence of the above considerations, Eq. (10) can be cast into the form (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 (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 and the argument Θ_i is given by (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⁶ (e.g., Simon et al. 1994). For our purposes, can be considered as an affine function of time of the form (21)The functions B_i, or in short, the orbital functions, are defined (Kinoshita 1977) as (22)with (23)The rates , which render time-dependent, 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 (24)Typically, the largest value of the ratio B_i,1/B_i,0 is of the order 10^-3cy^-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 long-term 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 and 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) (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 (26)The angular velocity ϰ depends on the functions of time π_A and Π_A through (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 (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 non-rigid 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 Kinoshita-like 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 Ferraz-Mello (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 ℋ₀, related to the torque-free motion, plus a perturbation ℋ₁ characterized by the coefficients k′ and ϰ_i appearing in Eqs. (19) and (26). That is to say, we can write⁷(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 (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 short-period terms. Therefore, we have (31)The integral in Eq. (30) is evaluated along the trajectories of the unperturbed problem (UP). Here, this problem is given by the torque-free 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 torque-free 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 and B_i,1/B_i,0, we simply keep terms up to the first order in these parameters (32)Therefore, considering the constancy of I and λ in the unperturbed problem, and the time dependence of Θ_i and (Eqs. (20), (21), and (24)) the integrands reduce to integrals of the form (33)where a, α ≠ 0, and β are real constants. By doing so, the generating function can be written as (34)with (35)The first term in was computed in Kinoshita (1977) and coincides with his expression, except for , due to the different choice of the unperturbed problem⁸. It is also the case for the term in proportional to . The parts linear in appear in the generating function as a consequence of incorporating a J₂ rate into the Earth model. They were not considered in Kinoshita (1977) nor Souchay et al. (1999), since they arise from the non-rigidity.

The parameters are also present in as out-of-phase terms. Since their contribution to the periodic part of the nutations is very small, those out-of-phase 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) (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 short-period 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 long-period terms, they are usually expressed as (37)where , , and 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 , , , and 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 , ω_E being the angular velocity of the Earth (Kinoshita 1977). As pointed out by Williams (1994), processes leading to a secular variation of J₂ leave M^∗ unaffected, that is, , although both C and depend on time⁹.

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 (38)Considering Eqs. (3), (20), and (35), we have (39)Specifically, we get (40)with (41)These analytical formulae clearly reflect the influence of the J₂ rate and the orbital coefficients rate on the nutation model. These formulae also depend on the precession quantities I^∗ and λ^∗ (through the argument , Eq. (20)).

When taking in Eqs. (41), the amplitudes and coincide ¹⁰ 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 second-order 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 out-of-phase amplitudes and arise as a consequence of the J₂ and orbital coefficients rate and they are not present in other studies.

Finally, the computation of the nutations as quasi-polynomials in t is achieved by substituting I^∗ as an explicit function of time, in addition to the time dependence of the arguments (Eq. (20)). Within the scope of this study, we consider only the first order expansion to evaluate the amplitudes , , , and , with the constants and 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 first-order 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(42)The coefficients ℱ_i,1 have a double origin: one part is due to the rate of I^∗, , and the other comes from the rates of k′ and B_i.

With this relationship in mind, and keeping only first-order terms in the related parameters associated to the rates , B_i,1 (and its derivatives with respect to I^∗), , we can finally write the nutations in longitude and obliquity as (43)The terms with subscript 0 are trigonometrical polynomials with constant coefficients of the form (44)where (45)The terms with subscript 1 contain t in their amplitudes, that is, they are mixed secular terms¹¹, whose expressions are given by (46)where (47)

4. Discussion

4.1. Adopted parameters

The evaluation of Δ₀ψ, Δ₀ε, Δ₁ψ, and Δ₁ε 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 , , and .

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 and are taken as those adopted in IAU 2006 precession for ϵ_A expansion, with the proper change of sign, that is, and (Capitaine et al. 2005).

Besides, it is necessary to consider the value of the scaling factor for each perturber. This can be written as (48)The values of k₀ for the Moon and the Sun were provided in Souchay et al. (1999) for a rigid Earth model. The non-rigidity of the Earth entails a change of the value of H_d,0, so we have to consider (49)In particular, for the rigid Earth we take the values /cy, and cy, for the Moon and the Sun, respectively, and (Souchay et al. 1999). As regards the dynamical ellipticity of the non-rigid 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⁷ (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 Δ₀ψ, Δ₀ε, Δ₁ψ, and Δ₁ε through the parameter (Eq. (17)). It induces out-of-phase 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¹²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.

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¹³ 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 first-order contributions of gravitational origin arising from the J₂ 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 out-of-phase 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 out-of-phase 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¹⁴.

The relevant precession quantities entering into Δ₀ψ, Δ₀ε, Δ₁ψ, and Δ₁ε are and . 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¹⁵.

To compute the variations of the nutation amplitudes associated with the changes in and , we have calculated the difference among the amplitudes, evaluated with IAU 2006 parameters, and their counterparts obtained from IAU 1976 values and /cy provided in Lieske et al. (1977). The results are displayed in Table 3.

The changes of precession quantities provide some in-phase and mixed secular terms both in longitude and obliquity. There are no out-of-phase 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 in-phase components, especially relevant for the 18.6-yr term in longitude, as was first pointed out in Escapa et al. (2014).

Those in-phase changes can be viewed as a re-evaluation 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 and . Since these amplitudes are common with the nutations of the angular momentum axis, they can be directly incorporated to the non-rigid Earth nutation series.

The adjustments considered in Capitaine et al. (2005) affecting the in-phase 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₂ 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 in-phase, out-of-phase, and mixed secular terms both in longitude and obliquity. Specifically, the adjustments larger than 1μas are given by (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 non-rigid Earth developed within the Hamiltonian framework. It was constructed by considering the first-order 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 second-order terms. Earth internal structure affects mainly the nutations through the Oppolzer terms. For a two-layer Earth model, the largest contribution is about 80 mas (milliarcseconds) for the 18.6-yr 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, second-order terms are proportional to k^′2 and k′ϰ_i. The largest contribution associated with k^′2 is about 1250μas for the 9.3-yr 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 out-of-phase terms of the order of 250μas (Getino et al. 2010). In addition to and , they also depend on the precessional quantities s₁, s₂, c₁, and c₂ through ϰ_i (Eq. (28)), especially on c₁ (Getino et al. 2010). Their IAU 2006 precession values are different from the IAU 1976 precession ones. The relative change in c₁ 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₂ 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 out-of-phase 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.

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 AYA2010-22039-C02-02 and AYA2016-79775-P (AEI/FEDER, UE).

References

Appendix A: Orbital coefficients rate

The orbital coefficients arise from the development of the second-degree disturbing potential of the Moon and the Sun (e.g., Kinoshita 1977; Kinoshita & Souchay 1990). They are defined from the expressions

(A.1)where α and β denote the longitude and the latitude of the Moon or the Sun referring to the ecliptic of date and are the nth-degree, mth-order 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) (A.2)Hence, from Eqs. (A.1) and (A.2), and are determined through (A.3)the orbital coefficients being nil.

The Sun orbital coefficients rate, (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 (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 must be introduced. Hence, Eq. (A.4) leads to (A.5)These relationships allow the left hand sides of Eqs. (A.3 ) to be expressed in terms of the elliptic motion elements, since (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 semi-major axis, a. Namely, we have (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 (A.8)and by introducing the Bessel functions of the first kind and order s and their derivative, J_s and . Specifically, it can be deduced that (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 (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 and (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 right-hand sides of Eqs. (A.3).

The resulting expressions of this expansion are trigonometric polynomials of l_S and , 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¹⁶ can be written as (A.11)In this way, to the third order in the eccentricity, we have recomputed the orbital coefficients rate 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) (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.

The analytical expressions and numerical values of orbital coefficients rate 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 out-of-phase 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.

To make the re-calculation of the different numerical evaluations carried in this work more straightforward, we have gathered the used values of the coefficients 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 are taken from Kinoshita (1977) with the update provided in Kinoshita & Souchay (1990). The coefficients 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 and has no effect on coefficients, as was implicitly assumed in the Kinoshita & Souchay (1990) update.

All Tables