© ESO 2019

1. Introduction

Most of the theoretical contributions to the precessional rate of the Earth’s equator, included in the current IAU2006 theory (Capitaine et al. 2003, 2005), were taken from previous works by Williams (1994) and Mathews et al. (2002). The effect of the tidal deformations on the precession rates in longitude and obliquity – the main scope of this work – is included in the IAU2006 set of contributions with the unspecific name of non-linear effect, and it is established as −0.960 mas cy⁻¹ for the precession rate in longitude, and 0.340 mas cy⁻¹ in obliquity.

A recent publication by Liu & Capitaine (2017) has evaluated the possibility of upgrading IAU2006 precession theory. In that process, the updated values of the non-linear effect to be considered would be −0.124 and 1.844 mas cy⁻¹, as reported in Lambert & Mathews (2008). Liu & Capitaine (2017) designate those contributions as an effect of the second-order torque on precession rates.

In a previous study by the authors (Baenas et al. 2017), the revision of such components of the secular motion of the equator was initiated, and it was shown that the second-order lunisolar contribution changes significantly when the assumed rigid Earth model is replaced by a two-layer one that incorporates some of the effects of the deformation one (achieving a difference of −8.51 mas cy⁻¹ in the precession in longitude with respect to the IAU2006 reference value).

In this study, the IAU2006 contributions arising from the variations of the Earth’s gravitational potential induced by the deformations resulting from the action of the lunisolar gravitational potential (also known as tidal potential or tide-raising potential) are revisited. The physical phenomenon is referred to as tidal deformation (Love 1906, Sect. 185) or bodily tide (Jeffreys 1976, Chap. 7), and the additional term in the gravitational energy of the system is known as redistribution tidal potential (as it is caused by a mass redistribution). We show that the changes of the precession rates are not negligible when realistic hypotheses about the Earth’s rheology and dynamics are assumed. The resulting contributions are numerically comparable to some of the components of the current precession theory, displayed, for example, in Table 2 of Baenas et al. (2017).

The effects of the redistribution potential on Earth nutation and precession have been partially discussed by Souchay & Folgueira (2000), Escapa et al. (2003, 2004), Ferrándiz et al. (2012), and Baenas (2014), within a Hamiltonian framework. Another approach to the problem, based on an application of modified SOS (Sasao et al. 1980) equations, can be found in Lambert & Capitaine (2004) and Lambert & Mathews (2006, 2008). Krasinsky (1999) and, recently, Williams & Boggs (2016), computed the obliquity rate from the vectorial Euler equations. In the works by Lambert & Capitaine (2004) and Souchay & Folgueira (2000), a simplification of the problem, consisting in considering only the zonal contribution of the redistribution tidal potential, was performed. This hypothesis must be improved since the other harmonic contributions of the potential, the tesseral and sectorial terms, are of the same magnitude. The formalism used in Lambert & Mathews (2006, 2008) does not provide simple analytical formulae that allow for the recomputing of the precession and nutation due to the redistribution of mass when the Love numbers model is changed. That is the case, for example, when considering a Love numbers set that includes, in addition to the ocean load, the ocean tide direct effect on the redistribution (Williams & Boggs 2016 and this work). That inclusion produces significant differences with respect to the values provided in Lambert & Mathews (2006, 2008) as is shown in Sect. 4.

The works by Escapa et al. (2003, 2004), Ferrándiz et al. (2012) and Baenas (2014) constitute a progressive line that leads to the present article, in which the problem of the redistribution of mass in the evolution of the Earth’s rotation will be thoroughly studied by means of the Hamiltonian formalism and the use of canonical perturbation methods, allowing an analytical approach. With these tools, closed-analytical formulae for the precession rates will be obtained, which enables a simple numerical evaluation of the redistribution effects under different elastic responses of the Earth (or rheologies), characterized by means of the Love numbers formalism.

The effects of the mass redistribution are particulary relevant in the Earth’s precession modeling, due to its connection with the Earth’s dynamical ellipticity¹H_d = 1 − A/C (expression for an axial symmetric Earth, A and C being the principal moments of inertia, equatorial and polar, respectively). It should be remarked that, in spite of the simplified notation, A is taken in most of the theories (including IAU2000) as the mean (A + B)/2 of the equatorial moments of inertia. This simplification suffices for this paper, although it is pertinent to cite that the triaxial problem has been treated in many different papers; some examples are Chao et al. (1991), Getino et al. (2000), Escapa et al. (2002), Shen et al. (2007), Chen & Shen (2010), Bizouard & Zotov (2013), and Chen et al. (2015).

Independently of the symmetry assumed in the background model, the H_d parameter is inferred from the precession rate in longitude adjusted from astronomical observations at an epoch, p_A (general precession in longitude). Hence, it must be readjusted when some of the smaller components of p_A, which are treated in the process like known constants, are modified. Besides, corrections to the dynamical ellipticity indirectly affect the nutation amplitudes at the tens of μas level, through a consistency factor (H_d + δH_d)/H_d (Escapa et al. 2016).

The structure of the paper is as follows. In Sect. 2, the fundamentals of the Earth redistribution potential theory are described, with special emphasis in the process leading to its analytical formulation and its generalization through Love numbers formalism. In Sect. 3, the analytical formulae for the precession rates are obtained by means of the Hamiltonian approach. Section 4 is composed of the computation of the precession rates for different Earth models, given by the corresponding Love numbers sets. In Sect. 5, the main implications of the numerical results derived formerly are discussed. They refer to the permanent tide, the change in the dynamical ellipticity of the Earth, and its effects on the nutations. We also draw some of the fundamental conclusions of this research. Finally, two appendices have been added for clarifications and additional material.

2. Earth redistribution tidal potential

2.1. Basic hypothesis

The external gravitational potential energy of an extended body can be expanded in a series of spherical harmonics (e.g., MacMillan 1958, Sect. 172). That procedure is especially suitable when both the geometry and the mass density functions of the body are not far from the spherical symmetry, or when the potential energy is evaluated at great distances compared with the dimensions of the body.

From the point of view of continuum mechanics, the spherical harmonics expansion depends only on the volume of the body and its volumic mass distribution². This fact is neatly appreciated when the expansion is written in terms of the inertial integrals (Ibid, Sect. 204) given by

$$\begin{array}{c}\hfill {\displaystyle {\int }_B{x}_1^l{x}_2^m{x}_3^n\hspace{0.25em}\rho \left(\mathbf{r}\right)d^3\mathbf{r},}\hspace{0.25em}l,m,n=0,1,2,\dots \end{array}$$(1)

In this expression, r = (x₁,x₂,x₃)^t are the Cartesian coordinates of a material particle of the extended body with volume B and volumic density function ρ(r), and d³r stands for the volume element. Since the effects to be investigated are small, it is enough to consider just the second degree terms in the expansion of the external gravitational potential energy.

The relationships of the I_ij matrix elements of the inertia tensor of the body³,

$$\begin{array}{c}\hfill I_{\mathit{ij}}={\displaystyle {\int }_B\rho \left(\mathbf{r}\right)(r^2{\delta }_{\mathit{ij}}-x_i{x}_j)}{d}^3\mathbf{r},\end{array}$$(2)

with the second degree inertial integrals (l + m + n = 2) lead to the well-known expression of the second degree terms in the multipolar expansion of the external gravitational potential (MacCullagh’s formula, MacCullagh 1844; Peale 1973)

$$\begin{array}{c}\hfill \mathcal{V}\left({\mathbf{r}}_p\right)=-\frac{G{m}_p}{2r_p^5}{\displaystyle \sum _{i,j=1(j\ge i)}^3{I}_{\mathit{ij}}[r_p^2{\delta }_{\mathit{ij}}-(3-2{\delta }_{\mathit{ij}})x_{i,p}{x}_{j,p}]}.\end{array}$$(3)

Here r_p = (x_1, p,x_2, p,x_3, p)^t and r_p = |r_p| represent the position vector and the modulus of a material particle of mass m, labeled as p, where the gravitational potential energy is evaluated, and G is the gravitational constant.

Possible temporal dependence has not been made explicit in the previous formulae. The matrix of inertia, I=(I_ij), or the inertial integrals, Eq. (1), might depend on time, inducing a similar dependence on the expansion of the external gravitational potential. In turn, they can also depend on time due to an explicit temporal variation of r_p. Whereas the second situation does not present special difficulties, it is not the case for the first one, so it will be analyzed with some detail.

In such a case, time dependence in I can be due to two facts. First, the reference system employed in the computation of the matrix of inertia can be animated with a motion independent of that of the extended body. Second, the extended body can experience changes of shape and volume (redistribution of mass). Both situations would make the volume and volumic density of the extended body time dependent, B = B(t) and ρ = ρ(r,t), however, they have a quite different nature. The first one has a kinematical origin, while the second one is due to the deformable character of the extended body. A convenient way to separate, to some extent, those two facts is to attach in some prescribed way a reference system to the extended body.

The procedure is best appreciated if the case in which a rigid extended body is first considered. Since by definition the distances between its material particles must keep constant (invariable form and shape), all the time dependence in B(t) and ρ(r,t) emerges as a consequence of the motion of the reference system used in their description. Alternatively, the extended body temporal evolution is a rigid motion given by a combination of a translation of the barycenter of the body O, and a rotation around it.

It entails that there exist reference systems, or body figure axes, where I is independent of time. Among them, it is especially convenient to employ the principal axes system, Oxyz, because I is a diagonal matrix in it (for specific problems, another choice of the reference system might be more suitable in order to simplify the formulation of the rotational dynamics problem.).

In this way, the motion of the extended body can be completely identified with the motion of the reference system Oxyz and vice versa, that motion referred to another reference system, OXYZ, not necessarily inertial but with known time evolution. In the study of the rotation of the Earth, those other reference systems are also assumed to be geocentric like, for example, the one defined by the ecliptic of date.

The question is more complex for deformable bodies as was early recognized in Tisserand (1891, Tome II, Sect. 213). Following a similar procedure to that used in the case of the rigid body, first a reference system Oxyz can be linked to the extended body. It will induce a rigid motion common for all its material particles.

However, since the body is deformable, in general there will exist other part of the motion not described by that of Oxyz, and depending on each material particle. This will be referred to as deformation. In this way, the motion of the extended body is split into two components: a rotation around O (rigid motion) plus a deformation.

It allows the unambiguous definition of the rotation of a deformable body through the rotation of Oxyz relative to OXYZ (see, e.g., Escapa 2011). In contrast to the rigid case situation, however, the tensor of inertia I will depend on time due to the deformation component of the motion, even when it is described from Oxyz (which entails a variation of the relative distance among the material particles of the body).

The main problem with this approach is that the decomposition of the motion is not unique. Different criteria can be employed in order to define the rigid and deformable components (Munk & MacDonald 1960, Chap. 3, Sect. 2).

From the point of view of establishing the dynamics of the extended body, the prescription of the well-defined direct orthonormal Oxyz reference system called the Tisserand system is especially suitable (Ibid., Moritz & Mueller 1986; Kinoshita & Sasao 1977; see Escapa 2011 for clarifications on the equivalence of the different definitions).

In the Tisserand system⁴, the position r′ of every material particle in the deformed body (or deformed state of the body), with integration domain B′, will be obtained by means of a transformation defined by a Lagrangian particle displacement vector field, u(r), given by

$$\begin{array}{c}\hfill {\mathbf{r}}^{\prime }=\mathbf{r}+\mathbf{u}\left(\mathbf{r}\right).\end{array}$$(4)

The modulus of u(r) vectors will be assumed to be small enough to admit a first-order deformation theory, where any related magnitude is linear in the components of the displacement vector or its spatial derivatives. The Jacobian of the change of variable r′(r) related to the previous transformation, in the linear elasticity limit, is given by (Landau & Lifshitz 1959, Chap. 1)

$$\begin{array}{c}\hfill det\left(\frac{\partial {\mathbf{r}}^{\prime }}{\partial \mathbf{r}}\right)=1+\mathrm{\nabla }u\left(\mathbf{r}\right),\end{array}$$(5)

being the divergence of u known as volume dilatation, as it stands for the relative volume change in the deformation between volume elements, d³r′ = [1 + ∇u(r)]d³r. Besides, by means of the displacement field, the volumic density function ρ′(r′) in the deformed state can be also linked to its counterpart distribution in the undeformed domain (Moritz & Mueller 1986, Chap. 4),

$$\begin{array}{c}\hfill {\rho }^{\prime }\left({\mathbf{r}}^{\prime }\right)=\rho \left(\mathbf{r}\right)-\rho \left(\mathbf{r}\right)\mathrm{\nabla }u\left(\mathbf{r}\right).\end{array}$$(6)

The last result is not only a consequence of the definition of u(r), but of the continuity equation in the case of an Eulerian change of density due to an instantaneous deformation, δρ(r′) ≡ ρ′ (r′) − ρ(r′) = −∇[ρ(r)u(r)] (Backus 1967).

Equations (4)–(6) can be used to perform the calculation of the inertial integrals over the deformed state B′, ${\int }_{B^{\prime }}{x}_1^{\prime l}{x}_2^{\prime m}{x}_3^{\prime n}\hspace{0.25em}{\rho }^{\prime }\left({\mathbf{r}}^{\prime }\right)d^3{\mathbf{r}}^{\prime }$, and consequently the matrix elements of the inertia tensor of the deformed body, through an integration over the undeformed volume B. This is achieved by means of a first-order Taylor expansion of the product $x_1^{\prime l}{x}_2^{\prime m}{x}_3^{\prime n}$ induced by the transformations $x_i^{\prime }=x_i+u_i$ and the smallness of u_i in the linear elasticity limit, so that the following expression is obtained:

$$\begin{array}{cc}\hfill {\displaystyle {\int }_{B^{\prime }}{x}_1^{\prime l}{x}_2^{\prime m}{x}_3^{\prime n}\hspace{0.25em}{\rho }^{\prime }\left({\mathbf{r}}^{\prime }\right)}{d}^3{\mathbf{r}}^{\prime }=& {\displaystyle {\int }_B{x}_1^l{x}_2^m{x}_3^n\hspace{0.25em}\rho \left(\mathbf{r}\right)}{d}^3\mathbf{r}\hfill \\ \hfill & +{\displaystyle {\int }_B{x}_1^l{x}_2^m{x}_3^n\hspace{0.25em}(l\frac{u_1}{x_1}+m\frac{u_2}{x_2}+n\frac{u_3}{x_3})\rho \left(\mathbf{r}\right)}{d}^3\mathbf{r}.\hfill \end{array}$$(7)

The application of this result to the integrals of $I_{\mathit{ij}}^{\prime }$ elements (given by Eq. (2)), with primed variables in the case of a deformed extended body, allows us to decompose them into two parts: the matrix elements of the undeformed state, I_ij, and the perturbation due to the deformed one, ΔI_ij, performed over the undeformed B volume, namely,

$$\begin{array}{c}\hfill I_{\mathit{ij}}^{\prime }=I_{\mathit{ij}}+\mathrm{\Delta }{I}_{\mathit{ij}},i,j=1,2,3.\end{array}$$(8)

where (see, e.g., Chandrasekhar 1969, Sect. 15, or Wahr 1981)

$$\begin{array}{cc}& I_{\mathit{ij}}^{\prime }={\displaystyle {\int }_{B^{\prime }}{\rho }^{\prime }\left({\mathbf{r}}^{\prime }\right)(r^{\prime 2}{\delta }_{\mathit{ij}}-x_i^{\prime }{x}_j^{\prime })d^3{\mathbf{r}}^{\prime }},\hfill \\ \hfill & I_{\mathit{ij}}={\displaystyle {\int }_B\rho \left(\mathbf{r}\right)(r^2{\delta }_{\mathit{ij}}-x_i{x}_j)d^3\mathbf{r}},\hfill \\ \hfill & \mathrm{\Delta }{I}_{\mathit{ij}}={\displaystyle {\int }_B\rho \left(\mathbf{r}\right)[2(\mathbf{u}·r){\delta }_{\mathit{ij}}-(u_i{x}_j+u_j{x}_i)]d^3\mathbf{r}}.\hfill \end{array}$$(9)

In matrix form, Eq. (8) will be written as I′ = I + ΔI.

The substitution of Eq. (9) into Eq. (3) leads to a similar decomposition for the external gravitational potential of the extended body

$$\begin{array}{c}\hfill {\mathcal{V}}^{\prime }\left({\mathbf{r}}_p\right)=\mathcal{V}\left({\mathbf{r}}_p\right)+{\mathcal{V}}_t\left({\mathbf{r}}_p\right),\end{array}$$(10)

with

$$\begin{array}{cc}\hfill \mathcal{V}\left({\mathbf{r}}_p\right)& =-\frac{G{m}_p}{2r_p^5}{\displaystyle \sum _{i,j=1(j\ge i)}^3{I}_{\mathit{ij}}[r_p^2{\delta }_{\mathit{ij}}-(3-2{\delta }_{\mathit{ij}})x_{i,p}{x}_{j,p}]},\hfill \\ \hfill {\mathcal{V}}_t\left({\mathbf{r}}_p\right)& =-\frac{G{m}_p}{2r_p^5}{\displaystyle \sum _{i,j=1(j\ge i)}^3\mathrm{\Delta }{I}_{\mathit{ij}}[r_p^2{\delta }_{\mathit{ij}}-(3-2{\delta }_{\mathit{ij}})x_{i,p}{x}_{j,p}]}.\hfill \end{array}$$(11)

The additional term 𝒱_t (r_p) in the potential energy is the redistribution potential, since it arises from the redistribution of mass associated with the deformation field u(r). As it is derived from Eqs. (9) and (12), the form of 𝒱_t (r_p) is determined once u(r) is known. For example, in the case of the bodily tide of the Earth, it would be necessary to find the displacement vector field produced on its body by the attraction of the Moon and the Sun. The way of constructing the redistribution potential shown here is a generalization of the method proposed in Getino & Ferrándiz (1990, see also Escapa 2011).

Another widespread procedure to obtain 𝒱_t (r_p) is by solving directly the elastic problem, since the dynamical equations of motion of the displacement vector field must be supplemented with the Poisson equation for the redistribution potential. This is the approach followed by classical works on the bodily tide like, for example, those of Love (1906), Takeuchi (1950), and Alterman et al. (1959), leading to the introduction of Love numbers. This method allows the use of Eq. (3) in reverse: from the known form of the redistribution potential it is possible to obtain the variation of the inertia matrix due to deformation (e.g., Jeffreys 1976, Chap. 7, Peale 1973; Sasao et al. 1980).

2.2. Earth redistribution tidal potential for SNREI models

The SNREI (seismology acronym, Dahlen 1968) model refers to an (oceanless) Earth assuming a Spherical, Non-Rotating, Elastic and Isotropic body, in the undeformed state. SNREI Earth models have also been considered when determining the bodily tide (e.g., Love 1906; Takeuchi 1950). Although the real Earth departs from the SNREI model, those theories provide reasonable results for the tidal deformations (Wahr 1982), which makes them suitable to study the effects that the bodily tides have on other motions of the Earth.

It is the case when determining the effects of the tidal deformations on the Earth rotation (e.g., Jeffreys & Vicente 1957; Sasao et al. 1980; Getino & Ferrándiz 1995, or Escapa 2011). Since the variations of the Earth rotation, the displacement vector field associated with the tide, and the deviations of the SNREI Earth model with respect to the real Earth are small, ignoring rotational and elliptical effects on the bodily tide leads to a very good first approximation when studying Earth rotation⁵. This fact has been explicitly confirmed by theories of the Earth rotation that include first-order rotational and elliptical effects in the deformations, as it is stated in the very elucidating notes by Wahr (1982).

The features of a SNREI model imply the following assumptions. Firstly, the deformation is purely elastic, which implies instantaneity and reversibility⁶, both properties consistent with the discussion made in the previous section. Moreover, rotational and elliptical effects are not considered on the Earth’s deformation. In order to perform the ΔI_ij integrals given by Eq. (9), the density function in a SNREI model is assumed to be isotropic or spherical, ρ(r) = ρ(r).

In this situation, the case of a deformation induced by an external (gravitational) potential, proportional to solid spherical harmonics (of second degree), 𝒰₂, was studied by Love (1911) and Takeuchi (1950). In fact, let (r_q,θ_q,ϕ_q) be the spherical coordinates of a perturber q, the tidal (or tide-raising) degree-2 potential due to the gravitational action of the perturber q on a point (r,θ,ϕ) of the Earth is given by the well-known expression

$$\begin{array}{c}\hfill {\mathcal{U}}_2(\mathbf{r},{\mathbf{r}}_q)=-G\frac{m_q}{r_q}{\left(\frac{r}{r_q}\right)}^2{P}_2(cos\gamma ),\end{array}$$(12)

where γ is the angular separation between vectors r_q and r, and P₂ is the second degree Legendre polynomial. Then, the acceleration derives from the potential in the form −∇𝒰₂. The solution for the displacement field is given by (Love 1911, Sect. 109),

$$\begin{array}{c}\hfill \mathbf{u}\left(\mathbf{r}\right)=-F_2\left(r\right)\mathrm{\nabla }{\mathcal{U}}_2-G_2\left(r\right){\mathcal{U}}_2\mathbf{r},\end{array}$$(13)

with F₂ and G₂ being some functions of r only⁷. We note that the toroidal component (proportional to r × ∇𝒰₂) is null in this model. Therefore, using the addition theorem for the un-normalized degree-2 spherical harmonic (Whittaker & Watson 1950, Sect. 15.7) in Eq. (12),

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathcal{U}}_2(\mathbf{r},{\mathbf{r}}_q)=-G\frac{m_q}{r_q}{\left(\frac{r}{r_q}\right)}^2\times [{\mathcal{C}}_{20}{\mathcal{C}}_{20,q}}\hfill \\ {\displaystyle +\frac{1}{3}({\mathcal{C}}_{21}{\mathcal{C}}_{21,q}+{\mathcal{S}}_{21}{\mathcal{S}}_{21,q})+\frac{1}{12}({\mathcal{C}}_{22}{\mathcal{C}}_{22,q}+{\mathcal{S}}_{22}{\mathcal{S}}_{22,q})}],\hfill \end{array}\end{array}$$(14)

is obtained, 𝒞_2m and 𝒮_2m functions being real surface spherical harmonics of second degree, 𝒞_2m (θ,ϕ) = P_2m (cos θ) cos (mθ), 𝒮_2m (θ,ϕ) = P_2m (cos θ)sin(mθ), and P_2m the associate Legendre function. In Eq. (14) the arguments of the spherical harmonics have been omitted, for the sake of simplifying the notation. We note that those with q subscript are related to the spherical coordinates of the perturber.

Taking into account the previous considerations, the integral of the perturbation ΔI_ij to the inertia matrix elements (Eq. (9)) can be performed. This requires a rewriting of the Cartesian coordinates u_i of the displacement vector in spherical ones, given by

$$\begin{array}{cc}\hfill u_r& =-[2\frac{F_2\left(r\right)}{r}+G_2\left(r\right)r]{\mathcal{U}}_2,\hfill \\ \hfill u_{\theta }& =-\frac{F_2\left(r\right)}{r}\frac{\partial {\mathcal{U}}_2}{\partial \theta },\hfill \\ \hfill u_{\varphi }& =-\frac{F_2\left(r\right)}{rsin\theta }\frac{\partial {\mathcal{U}}_2}{\partial \varphi },\hfill \end{array}$$(15)

where the relation ∂𝒰₂/∂r = 2𝒰₂/r has been used, in addition to the spherical expression of ∇𝒰₂. Then, performing the integral in the angular variables (θ,ϕ), the following expression of the perturbation of the matrix of inertia due to the perturber q is obtained:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathbf{I}}_q=2{\kappa }_{t,q}{\left(\frac{a_q}{r_q}\right)}^3\left(\begin{array}{ccc}{\mathcal{C}}_{20,q}-\frac{1}{2}{\mathcal{C}}_{22,q}& -\frac{1}{2}{\mathcal{S}}_{22,q}& -{\mathcal{C}}_{21,q}\\ -\frac{1}{2}{\mathcal{S}}_{22,q}& {\mathcal{C}}_{20,q}+\frac{1}{2}{\mathcal{C}}_{22,q}& -{\mathcal{S}}_{21,q}\\ -{\mathcal{C}}_{21,q}& -{\mathcal{S}}_{21,q}& -2{\mathcal{C}}_{20,q}\end{array}\right),\end{array}$$(16)

with κ_t, q being a parameter describing the elastic response of the SNREI Earth to the perturbation caused by the perturber, whose integral form is

$$\begin{array}{c}\hfill {\kappa }_{t,q}=\frac{4\pi }{15}\frac{G{m}_q}{a_q^3}{\displaystyle {\int }_0^{a_{\mathrm{E}}}{r}^4\rho \left(r\right)[5F_2\left(r\right)+r^2{G}_2\left(r\right)]}dr.\end{array}$$(17)

In Eqs. (16) and (17), the mean distance between the geocenter and the perturber, a_q, has been introduced by normalization convenience; a_E is the Earth’s mean equatorial radius, used to define the integral in the radial coordinate. Equation (16) is equivalent to those of Peale (1973), Kubo (1991), Getino & Ferrándiz (1990, 1995), or Lambert & Mathews (2006), but the parameters describing the elastic behavior are different up to proportionality factors.

However, the process that we have followed in their computation is different. In those works Eq. (16) is obtained from MacCullagh’s formula (3), since the redistribution tidal potential is known from the solution of the elastic problem. Here, we construct the increment of the matrix of inertia through the displacement vector. This approach provides an alternative derivation to that presented in Getino & Ferrándiz (1990). It should be recalled that the trace(ΔI) is null, due to a general result of linear elasticity (Rochester & Smylie 1974, generalizing Darwin’s theorem, Darwin 1910, for a solenoidal deformation).

2.3. Expression of the redistribution tidal potential

MacCullagh’s formula, Eq. (3), allows the linking of the perturbation of the inertia tensor to the additional gravitational potential energy related to the mass redistribution. Therefore, from Eq. (11), the redistribution tidal potential energy, 𝒱_t; p, q, due to a couple of perturbed (p) and perturbing (q) bodies is written as

$$\begin{array}{c}\hfill {\mathcal{V}}_{t;p,q}=-G\frac{m_p}{2r_p^5}{\displaystyle \sum _{i,j=1(j\ge i)}^3\mathrm{\Delta }{I}_{ij,q}[r_p^2{\delta }_{\mathit{ij}}-3(2-{\delta }_{\mathit{ij}})x_{i,p}{x}_{j,p}]}.\end{array}$$(18)

Here r_p = (x_1, p,x_2, p,x_3, p) are the Cartesian coordinates of the point where the gravitational potential is evaluated; in our application, it will be the position of the perturbed body p. After some algebra, introducing the ΔI_ij, q elements given by Eq. (16), and taking into account that the solid spherical harmonics of the second degree can be expressed as linear combinations of products x_i, px_j, p, namely,

$$\begin{array}{cc}\hfill r_p^2{\mathcal{C}}_{20,p}& =\frac{1}{2}(3x_{3,p}^2-r_p^2),\hfill \\ \hfill r_p^2{C}_{22,p}& =3(x_{1,p}^2-x_{2,p}^2),\hfill \\ \hfill r_p^2{\mathcal{C}}_{21,p}& =3x_{1,p}{x}_{3,p},\hfill \\ \hfill r_p^2{\mathcal{S}}_{21,p}& =3x_{2,p}{x}_{3,p},\hfill \\ \hfill r_p^2{\mathcal{S}}_{22,p}& =6x_{1,p}{x}_{2,p},\hfill \end{array}$$(19)

the 𝒱_t; p, q potential becomes

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathcal{V}}_{t;p,q}=-G\frac{m_p}{a_p^3}{\left(\frac{a_p}{r_p}\right)}^3{\left(\frac{a_q}{r_q}\right)}^3{\kappa }_{t,q}\times [6{\mathcal{C}}_{20,p}{\mathcal{C}}_{20,q}}\hfill \\ {\displaystyle +2({\mathcal{C}}_{21,p}{\mathcal{C}}_{21,q}+{\mathcal{S}}_{21,p}{\mathcal{S}}_{21,q})+\frac{1}{2}({\mathcal{C}}_{22,p}{\mathcal{C}}_{22,q}+{\mathcal{S}}_{22,p}{\mathcal{S}}_{22,q})}].\hfill \end{array}\end{array}$$(20)

It will be useful for the purpose of this study to link the κ_t, q parameter with the Love number k₂ (Love 1911, Sect. 59), due to the widespread use of such formalism in the study of the Earth’s deformability, as will be shown in the next sections. By definition, k₂ Love number is given by the proportionality between the additional redistribution potential at the (displaced) surface, arising solely from the redistribution of mass (Munk & MacDonald 1960, Chap. 5).

Therefore, in order to use Eqs. (14) and (20), the following relation will be fulfilled, 𝒱_t; p, q(r_p = a_E)/m_p = k₂𝒰₂(r = a_E), which simplifies to

$$\begin{array}{c}\hfill k_2=6\frac{a_q^3}{m_q{a}_{\mathrm{E}}^5}{\kappa }_{t,q}.\end{array}$$(21)

Obviously, k₂ is independent from the perturber q, and characterizes the elastic response of Earth in a SNREI model. The well-known form of the redistribution potential is then recovered,

$$\begin{array}{c}\hfill \frac{{\mathcal{V}}_{t;p,q}}{m_p}=k_2{\left(\frac{a_{\mathrm{E}}}{r_p}\right)}^5{\mathcal{U}}_2({\mathbf{r}}_p,{\mathbf{r}}_q),\end{array}$$(22)

with 𝒰₂ given by Eq. (14). This general expression is used, for instance, in Munk & MacDonald (1960), Peale (1973) or, in recent works, by Efroimsky (2012a,b) or Williams & Boggs (2016).

The tide-raising potential is commonly referred to the Earth surface (e.g., Kaula 1964; Efroimsky 2012a,b, or Williams & Boggs 2016). By writing r_p = r_pe_p and with the aid of Eq. (14), we derive that

$$\begin{array}{c}\hfill {\mathcal{U}}_2({\mathbf{r}}_p,{\mathbf{r}}_q)={\left(\frac{r_p}{a_{\mathrm{E}}}\right)}^2{\mathcal{U}}_2(a_{\mathrm{E}}{\mathbf{e}}_p,{\mathbf{r}}_q)={\left(\frac{r_p}{a_{\mathrm{E}}}\right)}^2{\mathcal{U}}_2\left({\mathbf{r}}_q\right),\end{array}$$(23)

where we have introduced the abridged notation 𝒰₂(a_Ee_p,r_q) = 𝒰₂(r_q). In this way, we can write Eq. (22) in the form

$$\begin{array}{c}\hfill \frac{{\mathcal{V}}_{t;p,q}}{m_p}=k_2{\left(\frac{a_{\mathrm{E}}}{r_p}\right)}^3{\mathcal{U}}_2\left({\mathbf{r}}_q\right).\end{array}$$(24)

At the Earth surface r_p = a_E; the former equation provides the Love number k₂ as the constant of proportionality between the redistribution potential and the tide-raising potential on the spherical Earth surface. In what follows, the criterion given by Eq. (24) will be adopted, due to its wider use.

From Eq. (17), an integral form for the k₂ Love number in the SNREI Earth model can be provided,

$$\begin{array}{c}\hfill k_2=\frac{8\pi }{5}\frac{G}{a_{\mathrm{E}}^5}{\displaystyle {\int }_0^{a_{\mathrm{E}}}{r}^4\rho \left(r\right)[5F_2\left(r\right)+r^2{G}_2\left(r\right)]}\mathrm{d}r.\end{array}$$(25)

This expression allows a numerical integration in order to obtain the k₂ value for different Earth rheologies (within the SNREI modeling). The Takeuchi (1950) computational procedure to obtain F₂ and G₂ functions was followed, as an example of a similar procedure, by Getino & Ferrándiz (1991) ⁸ – with Takeuchi’s Earth model – or Getino (1993) – with 1066A and 1066B Gilbert & Dziewonski (1975) models. However, this is not the purpose of this research, and the numerical values of k₂ (and those of its generalizations for more realistic Earth models) will be considered instead, as inputs of the theory, as occurs with the ephemeris describing the motion of the perturbers.

The time dependence of the redistribution tidal potential (Eq. (20)) is due to that of the Love number and the relative motion of the perturber q. It is better described in the frequency domain. The dependence of k₂ with the tidal forcing frequencies, due to the relative motion of the perturber q in a SNREI Earth model, is very weak at least when they are smaller than the free elastic oscillations of the Earth model as is the case (e.g., Wahr 1982 or Moritz & Mueller 1986, Chap. 4). It means that when deriving the expression of the displacement vector, it is possible to assume that the tide-raising potential is a constant function, that is, its forcing frequency equals zero. In other words, the elastic response of the SNREI Earth model is an equilibrium response giving rise to what is usually referred to as the static, or quasi-static, approximation (Moritz & Mueller 1986, Chap. 4). In entails that with those provisos the Love number k₂ is constant.

With respect to the part related to the perturber q, its resulting frequency set arises from a combination of the orbital evolution features of q with the Earth rotational motion. This kind of Fourier decomposition is common in tidal studies, leading to tide generating potential (TGP) expansions, and can be accomplished from multiple approaches (e.g., Doodson 1922; Kaula 1964, or Kudryavtsev 2004).

For the purposes of this study, however, it is more convenient to obtain the Fourier decomposition with the same method as that introduced by Kinoshita (1977) when developing the disturbing function in his theory of the Hamiltonian rotation of the rigid Earth. In this way, all the considered perturbations on the Earth rotation coming from the Moon and the Sun are developed from an identical baseline expansion, which enhances the consistency of the canonical approaches.

2.4. Generalization to the Earth redistribution potential

For the consideration of a more general elastic response of the Earth, mainly two kinds of effects can be introduced within the scope of the elastic-viscoelastic analogy⁹: the anelastic behavior intrinsic to the rheological model, and the contributions of the ellipticity and rotation to the deformation, both involving the transition from the stationary Love number k₂ to its dynamic version (time dependent).

2.4.1. Anelastic response (rheological model)

The anelastic behavior implies that 𝒱_t; p, q and 𝒰₂ in Eq. (24) are not determined at the same time, but with delay (without loss of reversibility). This fact can be written, with explicit arguments, as (Efroimsky 2012b)

$$\begin{array}{c}\hfill \frac{{\mathcal{V}}_{t;p,q}\left(t\right)}{m_p}={\left(\frac{a_{\mathrm{E}}}{r_p}\right)}^3{\displaystyle {\int }_{-\infty }^t{\dot{k}}_2(t-t^{\prime }){\mathcal{U}}_2({\mathbf{r}}_q,t^{\prime })\mathrm{d}{t}^{\prime }},\end{array}$$(26)

with k₂(t − t′) being the kernel of an integral representation of the dynamic counterpart of the Love number (the dot represents the time derivative), which acts as an operator mapping 𝒰₂ onto 𝒱_t; p, q. This formalism generalizes the classic idea of the constant time lag for the description of the delayed response relative to the perturbing potential, which is inconsistent in certain aspects as discussed in Efroimsky & Makarov (2013). We note that if k₂ (t − t′) = k₂Θ(t − t′) (Θ being the Heaviside step-function), Eq. (24) is recovered, describing an instantaneous reaction or, in other words, a perfectly elastic response.

In such an approach, the frequency-dependent complex Love functions, k̄₂(ω), are defined as Fourier transforms of ${\dot{k}}_2(t-t^{\prime })$ by means of

$$\begin{array}{c}\hfill {\dot{k}}_2(t-t^{\prime })={\displaystyle {\int }_{-\infty }^{\infty }{\overline{k}}_2\left(\omega \right){\mathrm{e}}^{i\omega (t-t^{\prime })}}\mathrm{d}\omega ={\displaystyle {\int }_{-\infty }^{\infty }\left|{\overline{k}}_2\left(\omega \right)\right|{\mathrm{e}}^{i[\omega (t-t^{\prime })+\epsilon \left(\omega \right)]}}\mathrm{d}\omega ,\end{array}$$(27)

where, in the second equality, the module-argument expression of the Love function, ${\overline{k}}_2\left(\omega \right)=\left|{\overline{k}}_2\left(\omega \right)\right|{\mathrm{e}}^{i\epsilon \left(\omega \right)}$, has been used.

2.4.2. Rotation and ellipticity (dynamical model)

The influence of rotation and ellipticity on the bodily tides has been long studied, for instance, within the rheological theories where the undeformed state is assumed to be ellipsoidal (or oblate) and rotating. This implies forgetting the simplifying hypotheses of the SNREI model. This is the case for Smith (1974); Shen & Mansinha (1976), and Wahr (1981) among others. Other approaches to the problem are those of Mathews et al. (1995), Dehant et al. (1999), or its generalization in the IERS Conventions (2010), both of them on the basis of the Sasao et al. (1980) dynamical theory (known as SOS formalism).

As these works show, such type of behavior can be built with a frequency (or mode) dependence of the Love numbers or, more precisely, Love functions (Efroimsky 2012b), which are introduced in Eq. (26). Strictly speaking, this equation involves any mass redistribution effect that can be formally described with a 2-degree 𝒰₂-like potential, which can be incorporated into the external tide-generating potential. This is the case, for instance, for oceanic loads and tides (see, e.g., the Love number model for an Earth with oceans in Williams & Boggs 2016, which will be considered in Sect. 4.3).

In the presence of rotation and ellipticity there appear normal modes related to the Earth rotation (e.g., Escapa et al. 2001). Those modes depend on the unperturbed rotation state and the internal structure of the Earth. For example, for a two-layer Earth model with a solid mantle and a fluid core, the most important ones are the Chandler wobble (CW) and the retrograde free core nutation (RFCN, or simply FCN). If we consider that the core contains a solid part, the solid inner core, we have two additional modes: the inner core wobble (ICW) and the prograde free core nutation (PFCN).

Among all those modes the FCN has a particularly notable effect on the Earth rotation. The reason is that its associated frequency relative to the Earth system, also referred to as nearly-diurnal free wobble (NDFW), is in the diurnal band¹⁰, as a consequence of the value of the ellipticity in the core-mantle boundary. Since there are also some frequencies of the tide-generating potential in the same band, a resonance phenomenon appears. Its presence has been modeled in different Earth rotation models (e.g., Poincaré 1910; Wahr 1981, or Getino & Ferrándiz 2001) and its effects on the nutations observed (e.g., IERS Conventions 2010, Chap. 5).

The excitation of the FCN mode by the forcing frequencies, related to the perturber q, couples with the Earth tidal deformation. It makes the Love number k₂ resonant at the FCN frequency (e.g., Jeffreys 1976, Chap. 7, Wahr 1982; Moritz & Mueller 1986, Chap. 4). That resonance of k₂ in the diurnal band has emerged in many studies using quite different approaches (e.g., Jeffreys & Vicente 1957; Sasao et al. 1980, or Wahr 1981).

The mechanisms entering into this process are complex and lead to a resonance expansion of the Love number k₂, as it is explained in the comprehensive work by Mathews et al. (1995). Recently, this problem has been revisited within the Hamiltonian formalism by Baenas et al. (2017), achieving closed-analytical formulae compatible with the previous works. That expansion can be written in the general form (e.g., Mathews et al. 1995) in the frequency domain as

$$\begin{array}{c}\hfill {\overline{k}}_2\left(\omega \right)=L_0+{\displaystyle \sum _i\frac{L_i}{\omega -{\omega }_i}},\end{array}$$(28)

where the sum runs over the rotation normal modes of the Earth model. The numbers L₀ and L_i are constant and ω_i is the frequency of the i normal mode relative to the Earth system. In the context of this research, the most relevant contribution to that resonance formula comes from the FCN mode.

The effects of rotation and ellipticity, by themselves, are conservative, so the numbers in Eq. (28) are real. However, if some dissipative torque is incorporated into the Earth rotation model, those constant numbers become complex. That induces a small imaginary part in ${\overline{k}}_2\left(\omega \right)$ or, equivalently, a phase ε(ω) as in the anelastic situation. That is the case, for example, when considering a viscous-electromagnetic torque at the core-mantle boundary (e.g., Getino & Ferrándiz 2001).

There are different resonance expansions generalizing Eq. (28) according to the features of the Earth rotation model under consideration. Typical examples are those provided in Mathews et al. (1995) or in IERS Conventions (2010, Chap. 6).

2.4.3. Earth redistribution potential

Both effects described in Sects. 2.4.1 and 2.4.2 can be incorporated into a sole formalism for the Love functions by using complex Love functions (similar to the complex compliances in SOS formalism), within the language of complex Fourier expansion of the involved potentials used by the IERS Conventions (2010). Therefore, Eq. (24) can be written in terms of the Fourier complex components of 𝒱_t; p, q and 𝒰₂ in the frequency domain as

$$\begin{array}{c}\hfill \frac{{\overline{\mathcal{V}}}_{t;p,q}\left(\omega \right)}{m_p}={\overline{k}}_2\left(\omega \right){\left(\frac{a_{\mathrm{E}}}{r_p}\right)}^3{\overline{\mathrm{U}}}_2({\mathbf{r}}_q,\omega ).\end{array}$$(29)

Equation (27) and (29) come from a standard Fourier expansion and details can be consulted in Efroimsky (2012b). In our formalism, the concise notation ω for the excitation frequency will be substituted by mω_E ± n_j, m being the frequency band (spherical harmonics order), ω_E the mean value of the Earth’s angular velocity, and n_j the j-th orbital frequency. These frequencies are obtained from a Fourier-like expansion of the orbital motion of the perturbing bodies, referred to the terrestial reference system (Kinoshita 1977, see Appendix A). Complex Love functions, usually named as complex Love numbers, will be denoted as

$$\begin{array}{c}\hfill {\overline{k}}_{2m,j}=\left|{\overline{k}}_{2m,j}\right|{\mathrm{e}}^{i{\epsilon }_{2m,j}}.\end{array}$$(30)

It should be remembered that the frequency dependencies ${\overline{k}}_{2m,j}$ must be given by the addition of a rheological and a dynamical model. Within a first-order theory, the solution offered by both models is a known function of time. In other words, ${\overline{k}}_{2m,j}$ values and the related frequencies are not affected by their own corrections on the Earth’s rotational motion, as happens with the orbital motion of the perturbers.

It should be recalled that, by relying on Eq. (29), the redistribution potential energy is therefore generalized by extracting the real part of the complex ${\overline{\mathcal{V}}}_{t;p,q}$ function in Eq. (29). Then, the following expression including arguments of the trigonometric functions within the spherical harmonics holds:

$$\begin{array}{cc}\hfill {\mathcal{V}}_{t;p,q}& =-G{m}_p{m}_q{\displaystyle \sum _j\frac{a_{\mathrm{E}}^5}{r_q^3{r}_p^3}}\hfill \\ \hfill & \times \{\left|{\overline{k}}_{20,j}\right|cos{\epsilon }_{20}{\mathcal{C}}_{20,p}({\theta }_p,{\varphi }_p){\mathcal{C}}_{20,q}({\theta }_q,{\varphi }_q)\hfill \\ \hfill & +\frac{1}{3}\left|{\overline{k}}_{21,j}\right|\left[\begin{array}{c}{\mathcal{C}}_{21,p}({\theta }_p,{\varphi }_p){\mathcal{C}}_{21,q}({\theta }_q,{\varphi }_q-{\epsilon }_{21,j})\\ +{\mathcal{S}}_{21,p}({\theta }_p,{\varphi }_p){\mathcal{S}}_{21,q}({\theta }_q,{\varphi }_q-{\epsilon }_{21,j})\end{array}\right]\hfill \\ \hfill & +\frac{1}{12}\left|{\overline{k}}_{22,j}\right|\left[\begin{array}{c}{\mathcal{C}}_{22,p}({\theta }_p,{\varphi }_p){\mathcal{C}}_{22,q}({\theta }_q,{\varphi }_q-\frac{{\epsilon }_{22,j}}{2})\\ +{\mathcal{S}}_{22,p}({\theta }_p,{\varphi }_p){\mathcal{S}}_{22,q}({\theta }_q,{\varphi }_q-\frac{{\epsilon }_{22,j}}{2})\end{array}\right]\}.\hfill \end{array}$$(31)

We note that the complete expressions of the real spherical harmonics also depend on the orbital frequencies, once they are provided by some calculation of ephemeris (for perturbed and perturbing bodies). This fact is implicit in Eq. (31). The total redistribution potential energy, 𝒱_t, is the sum over p and q – both representing either Moon or Sun – of terms 𝒱_t; p, q. Each portion can be also partitioned in the sum of zonal (long-period, m = 0), tesseral (diurnal, m = 1), and sectorial (semidiurnal, m = 2) harmonic components of the tide-raising potential.

3. Precession formulae

Once the expression of the redistribution potential energy has been derived in a quite general form (Eq. (31)), it is necessary to quantify its contribution to the Earth precession rates in longitude and obliquity. Although other approaches are possible, the formulation of this problem in the Hamiltonian formalism of the rotation of the non-rigid Earth (e.g., Getino & Ferrándiz 1995, 2001; Escapa et al. 2001, or Escapa 2011) presents some advantages (e.g., Baenas et al. 2017).

Two of the most relevant are the achievement of a consistent theory because all the effects are derived from a sole function, the Hamiltonian of the system, and the suitability of this framework to obtain analytical expressions even for the second-order terms, that is, terms that are quadratic in the small parameters. This is specially useful from our perspective, since the contributions of the redistribution of mass to the precession can be easily recomputed for different rheologies and dynamical Earth models. Next, the main elements necessary for this construction will be sketched. Further details can be found in two recent works by the authors, Baenas et al. (2017) and Escapa et al. (2017), and references therein.

3.1. Canonical formulation of the rotational problem of a non-rigid two-layer Earth

An Andoyer-like set of canonical variables (Getino 1995) {λ,μ,ν,λ_c,μ_c,ν_c;Λ,M,N,Λ_c,M_c,N_c} is used to describe the rotation that transforms the non-rotating (OXYZ) system into the (Oxyz) system linked – in Tisserand’s sense – to the Earth, which is referred to as a terrestrial system. The coordinates {λ,μ,ν} and their conjugated moments {Λ,M,N} stand for the Andoyer variables of the whole Earth, and the rest, with subscript c, for the fluid outer core (FOC). The auxiliary angles σ (between the angular momentum and the figure axis) and I (between the angular momentum vector and the Z axis of the non-rotating or quasi-inertial system) are defined through Λ = Mcos I and N = Mcos σ. The order of magnitude of σ is about 10⁻⁶ rad (Kinoshita 1977). More details on these variables within the precession context can be found in Ferrándiz et al. (2004) or Baenas et al. (2017).

The Hamiltonian of the problem is given by

$$\begin{array}{c}\hfill \mathcal{H}={\mathcal{T}}_0+{\mathcal{V}}_0+\mathcal{E}+{\mathcal{T}}_t+{\mathcal{V}}_t,\end{array}$$(32)

where 𝒯₀ is the kinetic energy corresponding to the rigid mantle case (Poincaré model, Getino & Ferrándiz 2001); 𝒱₀ is the degree-2 gravitational perturbation exerted by Moon (M) and Sun (S) (𝒱₀ = 𝒰_2; M + 𝒰_2; S);ℰ is the non-inertial complementary term (due to ecliptic motion, Kinoshita 1977); 𝒯_t is the redistribution kinetic energy (Getino & Ferrándiz 2001), and 𝒱_t is the redistribution potential energy, given in Sect. 2.4.3.

The precessional motion is caused by the secular part of the Hamiltonian. This is obtained from Expression (32) by means of the elimination of the periodic terms through a Lie-Hori canonical transformation (Hori 1966; Ferraz-Mello 2007), which leads to the averaged Hamiltonian, responsible for the secular motion of the equatorial plane.

The contribution of each term can be studied separately, obtaining the final values of the precession rates from summing up all the effects. The precession rate in longitude corresponding to 𝒯₀ + 𝒱₀ + ℰ (Poincaré model) was computed up to the second-order of perturbation by Ferrándiz et al. (2004). The effect of the mantle elasticity, introduced through the 𝒯_t term, was performed by Baenas et al. (2017), which is a second-order effect¹¹ in the scope of the perturbation methods, and updates the value of the general precession in longitude given by the previous work.

Therefore, this work’s aims are restricted to the secular Hamiltonian

$$\begin{array}{c}\hfill {\mathcal{H}}_{\sec }={\mathcal{V}}_{t,\sec },\end{array}$$(33)

whose order of magnitude, O(𝒱_t/𝒯₀) ∼ 1.7 × 10⁻¹², allows the use of a first-order perturbation method¹².

The contribution of 𝒱_t, sec to the precessional motion is determined through the variation of the derivatives n_λ = dλ/dt, for the precession in longitude, and n_I = dI/dt, in obliquity. We will denote by δn_λ and δn_I the new additive contributions, stemming from the application of the dynamical equations to the secular Hamiltonian.

These equations are given by δ(df/dt) = {f,ℋ_sec}, f being a smooth function of the canonical set, while {−,−} stands for the Poisson bracket written down in the Andoyer-like canonical variables. In this case,

$$\begin{array}{cc}\hfill \delta n_{\lambda }& =\{\lambda ,{\mathcal{H}}_{\sec }\}=-\frac{1}{MsinI}\frac{\partial {\mathcal{V}}_{t,\sec }}{\partial I},\hfill \\ \hfill \delta n_I& =\{I,{\mathcal{H}}_{\sec }\}=-\frac{1}{MsinI}(cosI\frac{\partial {\mathcal{V}}_{t,\sec }}{\partial \mu }-\frac{\partial {\mathcal{V}}_{t,\sec }}{\partial \lambda }),\hfill \end{array}$$(34)

is obtained.

It should be noted that the contributions δn_λ and δn_I are strictly related to the angular momentum axis motion. Within our order of approximation, those rates are the same as those of the longitude and obliquity of the figure axis (Kinoshita 1977). Besides, the commonly used astronomical longitude and obliquity angles, ψ = −λ_f and ε = −I_f, will be henceforth employed in this work to express the precessional rates, δn_ψ = −δn_λ and δn_ε = −δn_I.

3.2. Canonical expression of the secular redistribution potential

Considering Eq. (31), the real spherical harmonics must be expressed in the Andoyer-like canonical set of variables of the Earth, as a combination of spherical harmonics defined in the OXYZ system. In this work, the Moon and the Sun are considered as the sole perturbing bodies. Although some planets have an appreciable direct effect on Earth precession, as first demonstrated by Vondrák (1982), their direct influence on the Earth’s mass redistribution precession can be neglected within the current accuracy threshold¹³.

In this way, when the orbital motions of the Moon and the Sun are provided by some ephemeris (Kinoshita 1977), the redistribution tidal potential energy can be re-expressed in terms of Andoyer-like canonical variables. This process, at first order in σ, is given in Eq. (A.1). The related expressions are applicable to both perturbing and perturbed bodies, then a different notation will be used in order to distinguish them: dummy indexes i, τ for perturbed bodies, and j, ϵ for perturbing ones.

As the orbital motion of the perturbing bodies is considered to be decoupled from the Earth’s rotational motion, the related Andoyer coordinates, $\stackrel{\sim }{\lambda },\stackrel{\sim }{\mu },\stackrel{\sim }{\nu }$, and the fundamental argument, ${\stackrel{\sim }{\mathrm{\Theta }}}_j={\mathrm{\Theta }}_{0j}-m_{5i}\stackrel{\sim }{\lambda }$ (defined in Appendix A), are not canonical variables, but explicit time-dependent functions. The tilde symbol (∼) will be used from now on to denote time-dependent functions or, in other words, magnitudes depending on the variables or coordinates of the perturbers q.

The forcing orbital frequencies are given by ${\stackrel{\sim }{n}}_j=\mathrm{d}{\stackrel{\sim }{\mathrm{\Theta }}}_j/\mathrm{d}t$. The spherical harmonics expansion in such case (with the tilde symbol) can be truncated at the zero order in $\stackrel{\sim }{\sigma }$ because of the small magnitude of this angle as a function of time, and the fact that it does not participate in the construction of the equations of motion (e.g., Kubo 1991, appendix).

The secular part of the redistribution potential energy, 𝒱_t, sec, is computed in order to be free from short-period terms, following a similar procedure to that of Kinoshita (1977). This will be achieved by making explicit the canonical variables in Eq. (31) through Eq. (A.1), and then retaining only the terms whose frequency vanishes. In this procedure, some canonical variables can be substituted by their values in the steady rotation state. These are (see, e.g., Escapa 2011)

$$\begin{array}{c}\hfill \lambda ={\lambda }_0,N=N_0,M=M_0,\mu +\nu ={\omega }_{\mathrm{E}}t+{(\mu +\nu )}_0,\end{array}$$(35)

where the subscript 0 stands for integration constants, N₀ ≃ M₀ = Cω_E, and consequently, I = I₀, σ = σ₀ ≃ 0 (note that the remaining unperturbed solutions are not needed because they do not participate at the first order in σ, inasmuch as no derivatives with respect to σ are present in the dynamical equations for the precession).

Finally, the application of this technique results in

$$\begin{array}{c}\hfill {\mathcal{V}}_{t,\sec }={\mathcal{V}}_{t,\sec }^{(0)}+{\mathcal{V}}_{t,\sec }^{(1)}+{\mathcal{V}}_{t,\sec }^{(2)},\end{array}$$(36)

where the (0), (1), and (2) superscripts stand for the zonal, tesseral, and sectorial contributions of the redistribution potential, respectively, and are given by

$$\begin{array}{cc}\hfill {\mathcal{V}}_{t,\sec }^{(0)}& =-\frac{9}{4}C{\omega }_{\mathrm{E}}{\displaystyle \sum _{p,q=M,S}}{f}_q{k}_p{\displaystyle \sum _{i,j;\tau ,ϵ\in \mathcal{I}}}\left|{\overline{k}}_{20,j}\right|B_{i;p}{\stackrel{\sim }{B}}_{j;q}\hfill \\ \hfill & \times cos(\tau {\mathrm{\Theta }}_i-ϵ{\stackrel{\sim }{\mathrm{\Theta }}}_j+{\epsilon }_{20,j}),\hfill \\ \hfill {\mathcal{V}}_{t,\sec }^{(1)}& =-3C{\omega }_{\mathrm{E}}{\displaystyle \sum _{p,q=M,S}}{f}_q{k}_p{\displaystyle \sum _{i,j;\tau ,ϵ\in \mathcal{I}}}\left|{\overline{k}}_{21,j}\right|C_{i;p}{\stackrel{\sim }{C}}_{j;q}\hfill \\ \hfill & \times cos(\mu +\nu -\tau {\mathrm{\Theta }}_i-\stackrel{\sim }{\mu }-\stackrel{\sim }{\nu }+ϵ{\stackrel{\sim }{\mathrm{\Theta }}}_j+{\epsilon }_{21,j}),\hfill \\ \hfill {\mathcal{V}}_{t,\sec }^{(2)}& =-\frac{3}{4}C{\omega }_{\mathrm{E}}{\displaystyle \sum _{p,q=M,S}}{f}_q{k}_p{\displaystyle \sum _{i,j;\tau ,ϵ\in \mathcal{I}}}\left|{\overline{k}}_{22,j}\right|k_p{D}_{i;p}{\stackrel{\sim }{D}}_{j;q}\hfill \\ \hfill & \times cos(2\mu +2\nu -\tau {\mathrm{\Theta }}_i-2\stackrel{\sim }{\mu }-2\stackrel{\sim }{\nu }+ϵ{\stackrel{\sim }{\mathrm{\Theta }}}_j+{\epsilon }_{22,j}).\hfill \end{array}$$(37)

Here, ℐ set has been defined for the summation conditions

$$\begin{array}{c}\hfill \mathcal{I}=\{\tau ,ϵ\in \{-1,+1\}\mid \tau {\mathrm{\Theta }}_i-ϵ{\stackrel{\sim }{\mathrm{\Theta }}}_j=0\},\end{array}$$(38)

equivalent to the aforesaid cancelation of the frequency of the trigonometric arguments of 𝒱_t. In the previous expression and hereinafter, j subscript of Love and Kinoshita’s functions refers to an abridged notation for ϵ = ±1 and jth orbital frequency, which in fact should be denoted with a different subindex. However, this abuse of notation is kept due to the absence of likelihood of confusion.

The k_p parameter is the one defined by Kinoshita (1977), and the dimensionless f_q parameter has been introduced for notational convenience, namely,

$$\begin{array}{cc}\hfill k_p& =\frac{3G{m}_p}{{\omega }_{\mathrm{E}}{a}_p^3}{H}_{\mathrm{d}},\hfill \\ \hfill f_q& =\frac{m_q{a}_{\mathrm{E}}^2}{3C{H}_{\mathrm{d}}}{\left(\frac{a_{\mathrm{E}}}{a_q}\right)}^3.\hfill \end{array}$$(39)

3.3. Resulting precession formulae

Once Lie derivatives of Eqs. (34) have been computed, it is possible to identify the (μ,ν,λ) canonical variables with the $(\stackrel{\sim }{\mu },\stackrel{\sim }{\nu },\stackrel{\sim }{\lambda })$ time functions, given the coincidence between perturbed and perturbing bodies. This fact involves identifying the arguments Θ_i and ${\stackrel{\sim }{\mathrm{\Theta }}}_i$ likewise. The Kinoshita’s orbital functions, B_i, C_i, and D_i, can be evaluated at some epoch wherein I = Ĩ ≡ I₀, keeping Eq. (35) as a good approximation for this purpose (see, e.g., Escapa 2011). Finally, the precession rates read as

$$\begin{array}{cc}\hfill \delta n_{\psi }& =-\frac{1}{sinI}{\displaystyle \sum _{p,q=M,S}}{f}_q{k}_p{\displaystyle \sum _{\begin{array}{c}i,j;\tau ,ϵ\in \mathcal{I}\\ m=0,1,2\end{array}}}\left|{\overline{k}}_{2m,j}\right|T_{ijpq,m}^{\left(n_{\psi }\right)}cos{\epsilon }_{2m,j},\hfill \\ \hfill \delta n_{\epsilon }& =-\frac{1}{sinI}{\displaystyle \sum _{p,q=M,S}}{f}_q{k}_p{\displaystyle \sum _{\begin{array}{c}i,j;\tau ,ϵ\in \mathcal{I}\\ m=0,1,2\end{array}}}\left|{\overline{k}}_{2m,j}\right|T_{ijpq,m}^{\left(n_{\epsilon }\right)}sin{\epsilon }_{2m,j},\hfill \end{array}$$(40)

where the following functions have been defined

$$\begin{array}{cc}\hfill T_{ijpq,m}^{\left(n_{\psi }\right)}& =\frac{9}{4}\frac{\partial B_{i;p}}{\partial I}{\stackrel{\sim }{B}}_{j;q}{\delta }_{m0}+3\frac{\partial C_{i;p}}{\partial I}{\stackrel{\sim }{C}}_{j;q}{\delta }_{m1}+\frac{3}{4}\frac{\partial D_{i;p}}{\partial I}{\stackrel{\sim }{D}}_{j;q}{\delta }_{m2},\hfill \\ \hfill T_{ijpq,m}^{\left(n_{\epsilon }\right)}& =-\frac{9}{4}{B}_{i;p}{\stackrel{\sim }{B}}_{j;q}\tau m_{5i}{\delta }_{m0}+3C_{i;p}{\stackrel{\sim }{C}}_{j;q}(\tau m_{5i}-cosI){\delta }_{m1}\hfill \\ \hfill & +\frac{3}{4}{D}_{i;p}{\stackrel{\sim }{D}}_{j;q}(\tau m_{5i}-2cosI){\delta }_{m2}.\hfill \end{array}$$(41)

In these expressions, the ∼ symbol has been kept, although now it is strictly redundant.

Equations (40) and (41) are one of the main achievements of this work. They allow the numerical computation of the precessional rates due to the redistribution of mass once given the corresponding Love numbers of the model. In this way there is no need to recompute the redistribution potential energy itself.

4. Redistribution precession for different Earth Love numbers sets

In order to obtain numerical values for the precessional rates, the dependence of the Love numbers on the forcing frequencies is needed. In this research, some different models of Earth’s rheology have been considered.

The main parameters used to perform the calculations are included in Table 1. In turn, Table 2 separately displays the contributions coming from the different harmonic terms of the tidal potential – zonal, tesseral, and sectorial – for each Earth model.

In the different models under consideration, the contribution of the time-independent part of the redistribution potential will be discounted (further details will be discussed in Sect. 5.1). From Eq. (16) –the form of ΔI_q– and Eq. (A.1), the time-independent part in the Earth’s inertia tensor comes from the zonal contribution and is given by

$$\begin{array}{cc}\hfill \mathrm{\Delta }{C}_q& =-2C{H}_{\mathrm{d}}{f}_q\left|{\overline{k}}_{20,j}\right|{\left(\frac{a_q}{r_q}\right)}^3{\stackrel{\sim }{\mathcal{C}}}_{20}\hfill \\ \hfill & =-6C{H}_{\mathrm{d}}{f}_q{\displaystyle \sum _j\left|{\overline{k}}_{20,j}\right|{\stackrel{\sim }{B}}_{j;q}cos{\stackrel{\sim }{\mathrm{\Theta }}}_j},\hfill \end{array}$$(42)

when ${\stackrel{\sim }{\mathrm{\Theta }}}_j=0$ (zero frequency), or equivalently, with index j = (0, 0, 0, 0, 0). The portion of the deformation coming from this time-independent zonal contribution of the Earth’s inertia tensor is named permanent tide.

The permanent tide contribution to precession in Table 2 has been computed using the fluid (limit) Love number, k_2f = 0.94 (Burša & KarelPěč 1993, Sect. 4.9; Lambeck 1980, Sect. 2.4) and k_2f = 0.30 = k₂ for comparative purposes.

The fluid Love number arises from considering the Earth response to forces that stay constant over very long periods of the order of billions of years, not to tidal periodic ones. It entails considering a completely fluid Earth model instead of a typical SNREI Earth model, hence the difference in the values of tidal-effective k₂ and fluid¹⁴k_f Love numbers (e.g., Munk & MacDonald 1960, Chap. 5; Lambeck 1980, Chap. 2; or Moritz & Mueller 1986, Chap. 3).

When making comparisons between different studies involving Earth parameters like its inertia moments, ellipticity, J₂ Stokes coefficient, and others, it is very important to take into account in what way the effects of the permanent tide have been considered. It leads to the denominations zero-tide (including the permanent tide) or tide-free (excluding the permanent tide) systems, recommended in IERS Conventions (2010, Chap. 1). The zero-tide system is also known as the zero-frequency one (e.g., Burša & KarelPěč 1993, Chap. 4). Its use is recommended by the Resolutions of the International Association of Geodesy, as noted in the IERS Conventions (2010, Chap. 6).

In this regard, since tide-free values are unobservable and the value of k_f is uncertain, some authors prefer to consider the parameters in the zero-frequency system (Burša 1995). In our case, it would imply that the permanent tide contributions would be contained in the values of the moments of inertia of the undeformed state. In Eq. (42) and Kinoshita’s constant, Eq. (39), the dynamical ellipticity H_d is a zero-tide magnitude. The corresponding parameter of IAU2000 also includes the permanent tide (Groten 2004).

4.1. SNREI Earth

The tidal response of a SNREI model constitutes a topic of the literature that addresses the definition of Love numbers (e.g., Jeffreys 1976; Moritz & Mueller 1986, among many others). For the purposes of this study, its features can be characterized by taking a single constant value for the Love number set, ${\overline{k}}_{2m;j}=k_2$, and it can be analytically proved (Krasinsky 1999; Escapa et al. 2003, 2004; Lambert & Mathews 2006; Baenas 2014) that the contributions given by Eq. (40) to the precessional motion are cancelled out, δn_ψ = δn_ε = 0. A succinct demonstration of this fact can also be found in Appendix B. The first column in Table 2 shows numerically the cancelation of the harmonic contributions of the secular rates, on the last row of the table. It should be remarked that the cancelation is only achieved by adding the contribution of the permanent tide, assuming the same Love number as for the remaining tides, k_2f = k₂ = 0.30.

As displayed in Table 2, the obliquity rate for this model is zero. This result comes from Eq. (40), which show that this effect is only present when the Love numbers are complex, with any ε_2m, j ≠ 0. This is not the case for a perfect elastic, oceanless Earth model.

4.2. IERS Conventions (2010) models for solid Earth tides

4.2.1. Constant Love numbers per frequency band

In order to progressively show the influence on the precessional rates of the frequency dependence in the Earth’s anelastic response, a model with a nominal, reference value of the Love numbers for each harmonic order (or frequency band, long-period, diurnal, and semidiurnal) m, will be considered, that is, ${\overline{k}}_{2m;j}={\overline{k}}_{2m}$ (Col. 2 in Table 2). Basically, the band structure reflects the influence of rotation and ellipticity, whereas the imaginary part is related to the anelastic Earth response. This model is taken from Table 6.3 of IERS Conventions (2010) with ${\overline{k}}_{20}=0.30190-0.00000i$, ${\overline{k}}_{21}=0.29830-0.00144i$, and ${\overline{k}}_{22}=0.30102-0.00130i$.

4.2.2. IERS Conventions frequency-dependent Love numbers

In the IERS Conventions (2010), the frequency dependence in the Love numbers set (related to the rotation and ellipticity and the anelastic response of the mantle, Sect. 2.4) is introduced via the resonance formulae (Ibid, Eq. (6.9)), which update and summarize previous works by Wahr (1981), Mathews et al. (1995), and Dehant et al. (1999), among others. The modeling takes into account the retrograde free core nutation (FCN), the Chandler wobble (CW), and the free inner core nutation resonance processes in the diurnal band of frequencies for the solid Earth tides. In a second step, an effective ocean tide contribution (OT), $\delta k_{21;j}^{\mathrm{OT}}$ related to ${\overline{k}}_{21;j}$, following Wahr & Sasao (1981), is considered. This contribution incorporates the ocean load on the solid Earth tides. The ${\overline{k}}_{21;j}$ value is directly calculated from the nominal ${\overline{k}}_{21}=0.29830-i0.00144$ (Ibid., Table 6.3) and the tabulated δk_f; j values (Ibid., Table 6.5a), which are defined as

$$\begin{array}{cc}\hfill \delta k_{f;j}& =(L_0+{\displaystyle \sum _{k=1}^3\frac{L_k}{{\sigma }_j-{\sigma }_k}}-k_{21})+\delta k_{21;j}^{\mathrm{OT}},\hfill \\ \hfill {\overline{k}}_{21;j}& =k_{21}+\delta k_{f;j},\hfill \end{array}$$(43)

with L_k being the complex constants of the resonance formula for solid tides (Ibid., Eq. (69) and Table 6.4). In the zonal band, the IERS Conventions (2010) provide the following model, due to the mantle anelasticity (based in Wahr & Bergen 1986),

$$\begin{array}{c}\hfill {\overline{k}}_{20;j}=0.29525-5.795\times {10}^{-4}\{cot\frac{\alpha \pi }{2}[1-{\left(\frac{{\sigma }_r}{{\sigma }_j}\right)}^{\alpha }]+i{\left(\frac{{\sigma }_r}{{\sigma }_j}\right)}^{\alpha }\},\end{array}$$(44)

where σ_j is the excitation frequency (mω_E − ϵñ_j), σ_k and σ_r are reference frequencies, and α = 0.15 (Ibid., Table 6.4). The ${\overline{k}}_{22;j}$ complex Love number can be accurately approximated by its nominal value 0.30102 − i0.00130 (Ibid., Table 6.3). The related values of the precession rates components are shown in the third column (resonance formula) and fourth column (including the ocean load correction) of Table 2.

4.3. Williams & Boggs (2016) frequency-dependent Love numbers for the Earth with oceans

Finally, the last column in Table 2 shows the secular rates values obtained by using the Love numbers model of the Earth with oceans from Williams & Boggs (2016). This model is built over the IERS Conventions (2010) one for solid Earth tides (previous section), and the ocean model FES2004 (Lyard et al. 2006), giving a very complete description of the anelastic behavior of the tide effect. The Love numbers set can be consulted in Williams & Boggs (2016), Table 6, and the algorithm described on page 105.

This model allows us to compare the results obtained for the precessional rate in obliquity. Specifically, Williams & Boggs (2016) got the value for the precessional rate in obliquity (Ibid., Eq. (29)). They used a quite different approach from that followed here. In particular, they employed the Euler equations of vectorial mechanics, and constructed the Fourier expansion of the tide-generating potential with the same orbital ephemeris as us, but without introducing canonical variables in their derivations. By means of that procedure, Williams & Boggs (2016) derived a total tidal obliquity rate of 0.92 mas cy⁻¹ (which is split into −0.07 mas cy⁻¹, 0.1 mas cy⁻¹, and 0.89 mas cy⁻¹ for zonal, diurnal, and semidiurnal terms, respectively), while in this work an obliquity rate of 0.9341 mas cy⁻¹ (which is split into −0.0636 mas cy⁻¹, 0.0948 mas cy⁻¹, and 0.9030 mas cy⁻¹ for zonal, tesseral, and sectorial contributions, respectively) is obtained, showing an excellent agreement between both studies.

These equivalent results for the obliquity rate suppose a mutual validation of the procedures, due to them having been achieved through quite different approaches. In particular, it supports the correctness of Eqs. (40) and (41). Williams & Boggs (2016) do not provide a value for the longitude rate, as it is out of the scope of their work.

Lambert & Mathews (2006, 2008) also computed precession rates for a model close to that of Williams & Boggs (2016), but using the ocean tide model CSR4. The discrepancies with their numerical values are quite significant (see next section).

5. Discussion and conclusions

5.1. Contribution of the permanent tide

For the sake of convenience, the precession rate in longitude given by Eq. (40) will be decomposed in two parts, δn_ψ = p_P + p_P̄, p_P being the contribution of the permanent tide (zonal zero-frequency terms, i.e., those associated with B_0; p function), and p_P̄, the remaining non-permanent part (the “p” standardized notation is used for precession contributions). Therefore, from Eqs. (40) and (41), p_P can be written as

$$\begin{array}{cc}\hfill p_{\mathrm{P}}& =-\frac{1}{sinI}{\displaystyle \sum _{p,q=M,S}}{k}_{2f}{f}_q{k}_p{\displaystyle \sum _{\tau ,ϵ=\pm 1}}\frac{9}{4}\frac{\partial B_{0;p}}{\partial I}{\stackrel{\sim }{B}}_{0;q}\hfill \\ \hfill & =\left[\frac{1}{sinI}{\displaystyle \sum _{p=M,S}}{k}_p\frac{\partial B_{0;p}}{\partial I}\right][-9k_{2f}{\displaystyle \sum _{q=M,S}}{f}_q{\stackrel{\sim }{B}}_{0;q}].\hfill \end{array}$$(45)

Here, the $\left|{\overline{k}}_{20,0}\right|$ parameter has been simplified by k_2f fluid Love number, corresponding to permanent deformation. The first bracket of the decomposition coincides with the first-order lunisolar precession term, denoted as $p_A^{\prime }$ in Souchay & Kinoshita (1996) (an expanded form can be found in Kinoshita & Souchay 1990, Eq. (8.7)). The second one is related to the variation of the Earth dynamic ellipticity due to the permanent tide.

In fact, considering the tide-free dynamic ellipticity of a symmetric Earth (previous to tidal deformation), H_d, 0 = (C−A)/C, and its zero-tide counterpart, due to the permanent variation of the inertia moments,

$$\begin{array}{c}\hfill H_{\mathrm{d},\mathrm{P}}=\frac{C+\mathrm{\Delta }C-(A+\mathrm{\Delta }A)}{C+\mathrm{\Delta }C},\end{array}$$(46)

a first-order expansion in the small magnitude ΔC/C leads to the relation

$$\begin{array}{c}\hfill H_{\mathrm{d},\mathrm{P}}=H_{\mathrm{d},0}+\delta H_{\mathrm{d},\mathrm{P}}=H_{\mathrm{d},0}{\displaystyle (1-\frac{\mathrm{\Delta }C}{C})}+\frac{\mathrm{\Delta }C-\mathrm{\Delta }A}{C}·\end{array}$$(47)

Therefore, the δH_d, p following variation of the dynamic ellipticity holds

$$\begin{array}{c}\hfill \delta H_{\mathrm{d},\mathrm{P}}=\frac{\mathrm{\Delta }C(1-H_{\mathrm{d},0})-\mathrm{\Delta }A}{C}=\frac{\mathrm{\Delta }C}{C}{\displaystyle (\frac{3}{2}-H_{\mathrm{d},0})},\end{array}$$(48)

where, in the last equality, the relation ΔC + 2ΔA = 0 has been applied (Sect. 2.2). Recovering Eq. (42), with ΔC = ∑_qΔC_q,

$$\begin{array}{c}\hfill \frac{\delta H_{\mathrm{d},\mathrm{P}}}{H_{\mathrm{d},0}}=(-9+6H_{\mathrm{d},0})k_{2f}{\displaystyle \sum _q{f}_q{\stackrel{\sim }{B}}_{0;q}\simeq -9k_{2f}}{\displaystyle \sum _q{f}_q{\stackrel{\sim }{B}}_{0;q}},\end{array}$$(49)

which matches with the second bracket in Eq. (45). Hence, within a first-order precession theory (i.e., proportional to H_d, p), p_P can be written as

$$\begin{array}{c}\hfill p_{\mathrm{P}}=p_A^{\prime }\frac{\delta H_{\mathrm{d},\mathrm{P}}}{H_{\mathrm{d},0}}·\end{array}$$(50)

It should be recalled that H_d, P is a primary parameter (H_d, P = H_d, 0 + δH_d, P decomposition is not observationally distinguishable), while the tide-free H_d, 0 is a derived one. Accordingly, the p_P contribution should not be considered in the zero-tide system, since its effects are included within the H_d, P Earth dynamic ellipticity.

However, the precession rate due to the non-permanent tides, p_P̄, must be taken into account in correcting the reference value of H_d, P, as will be shown in the next section. In order to have a numerical estimate, p_P = 136.5964 mas cy⁻¹ (Table 2), and $p_A^{\prime }=5040684.593$ mas cy⁻¹ (obtained by subtraction in Capitaine et al. 2005) can be considered, thus

$$\begin{array}{c}\hfill \delta H_{\mathrm{d},\mathrm{P}}\simeq (8.8716\times {10}^{-8})k_{2f}.\end{array}$$(51)

We note that the approximately equal sign has been used as a reminder that the H_d, p value has been taken (Table 1) instead of the unobservable tide-free H_d, 0. Due to factor $p_{\mathrm{P}}/p_A^{\prime }\sim {10}^{-5}$ and the smallness of δH_d, P, this is an admissible approximation.

The second zonal Stokes parameter for a symmetric Earth, $J_{2,0}=(C-A)/\left(m_{\mathrm{E}}{a}_{\mathrm{E}}^2\right)$, is also consequently sensitive to the permanent tide. This effect has been previously studied, for instance, in Burša (1995), solving a boundary-value problem and leading to a numerical estimate of δJ_2, P = (3.08 × 10⁻⁸)k_2f (notation of this work).

In our case, by means of Eq. (50) it is only necessary to link the δJ_2, P variation with δH_d, P, which can be performed in a similar way to Escapa et al. (2017) (Eq. (16)) or combining Eqs. (48) and (50) with the definition of J_2, 0, H_d, 0, and $\delta J_{2,\mathrm{P}}=(\mathrm{\Delta }C-\mathrm{\Delta }A)/\left(m_{\mathrm{E}}{a}_{\mathrm{E}}^2\right)$. This relationship reads as

$$\begin{array}{c}\hfill \delta J_{2,\mathrm{P}}=J_{2,0}{{\displaystyle (1-\frac{2}{3}{H}_{\mathrm{d},0})}}^{-1}\frac{p_{\mathrm{P}}}{p_A^{\prime }}\simeq (3.1279\times {10}^{-8})k_{2f}.\end{array}$$(52)

Here, the zero-tide values for H_d, p and J_2, P = 1082635.9 × 10⁻⁹ (Groten 2004) have been used. This value can be considered in good agreement (a relative error of 1.5%) with that of Burša (1995), taking into account the modeling differences and the updating of physical parameters. In this sense, this fact also confirms the suitability of Eqs. (40) and (41) to treat the effects of the permanent tide.

5.2. Dynamical ellipticity variations due to non-permanent tides

In the general precession in longitude, p_A, the main component is the first-order lunisolar precession $p_A^{\prime }$ (see, for instance, Table 2 in Baenas et al. 2017). It is obtained by substraction of a set of higher order contributions, denoted as p_S, to the observed value of the precession, that is, $p_A=p_A^{\prime }+p_S$. The $p_A^{\prime }$ derived value allows the determination of the H_d, P dynamic ellipticity, through $p_A^{\prime }=\sum _p{k}_p{B}_{0;p}^{\prime }\left(I\right)/sinI$, recalling Eq. (39) for the linear relationship between k_p and H_d, P (see Kinoshita & Souchay 1990; Fukushima 2003 or Escapa et al. 2016, for more details on this procedure). The last expression will be formally and briefly written as $p_A^{\prime }=F\left(I\right)H_{\mathrm{d},\mathrm{P}}$.

If any change in the higher order part p_S is performed with respect to the reference value (in practice relative to IAU2006 precession one), the $p_A^{\prime }$ must be compensated by an equal and opposite change, and thus the H_d, P has to be recalculated. This is the case, for instance, for the recomputation of previously considered effects, like the second-order lunisolar precession (Baenas et al. 2017), the geodesic precession (Fukushima 2003), or the non-permanent p_P̄ component due to the redistribution of mass performed in this work, among others.

Let δp_S be a correction of p_S, the new value for the first-order lunisolar precession will be $p_A^{\prime }-\delta p_S$, and the variation the dynamical ellipticity, δH_d, δpS will be obtained from the equation $p_A^{\prime }-\delta p_S=F\left(I\right)(H_{\mathrm{d},\mathrm{P}}+\delta H_{\mathrm{d},\delta {\mathrm{p}}_{\mathrm{S}}}).$ Using that $F\left(I\right)=p_A^{\prime }/H_{\mathrm{d},\mathrm{P}}$, this equation yields the general expression

$$\begin{array}{c}\hfill \delta H_{\mathrm{d},\delta {\mathrm{p}}_{\mathrm{S}}}=-\frac{\delta p_S}{p_A^{\prime }}{H}_{\mathrm{d},\mathrm{p}},\end{array}$$(53)

pointed out, for instance, in Ferrándiz et al. (2012) or Baenas et al. (2017).

Considering δp_S = p_P̄ = −37.8472 mas cy⁻¹ (Table 2), the deviation from the IAU2006 reference value (−0.960 mas cy⁻¹, Ibid.) has an associate δH_d, δpS = 2.396 × 10⁻⁸, which implies a relative H_d, P change of +7.3 ppm. Therefore, $H_{\mathrm{d},\mathrm{P}}^{\mathrm{new}}=H_{\mathrm{d},\mathrm{P}}+\delta H_{\mathrm{d},\delta {\mathrm{p}}_{\mathrm{S}}}=0.00327381844$ is obtained. This value is somewhat higher if the Love numbers set for an Earth with oceans according to Williams & Boggs (2016) is considered. In particular, we have $H_{\mathrm{d},\mathrm{P}}^{\mathrm{new}}=0.00327382240$, with δH_d, δpS = 2.792 × 10⁻⁸ and a relative change of +8.5 ppm.

5.3. Indirect effects on nutation

Given a δp_S correction of p_S, the indirect effects on nutational motion emerge from the changes in the value of the dynamical ellipticity, δH_d, δpS (Escapa et al. 2016). Computed in the linear approximation, these effects are obtained by multiplying the nutation amplitudes by the factor (H_d, P+δH_d, δpS)/H_d, P.

The main terms arising from factor δH_d, δpS/H_d, P have been displayed in Table 3, considering the corrections to the dynamical ellipticity due to the second-order lunisolar precession (Baenas et al. 2017), and the redistribution potential (this work). The original nutation terms have been taken from Getino & Ferrándiz (2001).

The relative H_d, P change used in IERS2010 and WB2016 columns are those of the previous section. Some of the corrections are in the tens of μas level, and so they are relevant within the pursued accuracy of the current Earth rotation theories. Strictly speaking, these effects must be considered in the reconstruction of the Earth rotation theory, refitting the basic Earth parameters. However, the magnitudes provided in Table 3 are a clear indication that the effects of the redistribution tidal potential must be taken into account.

5.4. Conclusions

In this research the contributions to the IAU2006 precession of the mass redistribution processes induced on Earth by the gravitational action of the Moon and the Sun have been revisited. Closed analytical formulae for the precession rates in longitude and obliquity have been obtained, by means of the Hamiltonian formalism of the Earth rotation, through the consideration of the redistribution of tidal potential energy, whose derivation has been addressed in a comprehensive way from the changes in the Earth’s matrix of inertia. The analytical formulation enables the numerical computation of the precession rates in longitude and obliquity through Eqs. (40) and (41), under diverse Earth models or rheologies, characterized by different Love numbers that run from a SNREI model to those of an anelastic, elliptic and rotating Earth with oceans models.

The validation of Eqs. (40) and (41) has been shown in different contexts. For example, the precession rates formulae show the analytical cancelation in the SNREI model when the Love number is assumed to be constant regardless of the forcing frequency (including the long period through the permanent tide contribution). This result is consistent with other research (Krasinsky 1999; Lambert & Mathews 2006) that employed quite different methods. In addition, an excellent numerical agreement with Williams & Boggs (2016) in terms of the precession rate in obliquity for an Earth with oceans (with agreement in every frequency band) and Burša (1995) in relation to the permanent tide contribution to the zonal Stokes parameter, has been found.

The values displayed in Col. 4 of Table 2 are quite different from those computed Lambert & Mathews (2006, erratum Lambert & Mathews 2008). Although the authors do not specify which sets of Love numbers they used, it may be inferred that the basic elements used in their approach are similar to those employed in Williams & Boggs (2016). They got the value 1.84 mas cy⁻¹ (Lambert & Mathews 2008) replacing the former estimation of 0.13 mas cy⁻¹ (Lambert & Mathews 2006) for the obliquity rate. The value of 1.84 mas cy⁻¹ is quite close to 2.44 mas cy⁻¹ determined by Williams (1994). However, Williams & Boggs (2016) stated that the value 2.44 mas cy⁻¹ (Williams 1994) was derived under incorrect assumptions, because tidal forces have some components out of the ecliptic. Their new value of 0.92 mas cy⁻¹ is also consistent with similar results by Krasinsky (1999) and that obtained in this work.

In addition to the obliquity rate, our approach has allowed us to obtain the precession rate for the same model. Its numeric value also differs from that of Lambert & Mathews (2006, 2008). Since in those works explicit analytical formulae are not provided to compute the precession rates, it has not been possible to determine the source of those discrepancies. As we have pointed out, our Eqs. (40) and (41) are consistent with the results provided by other authors.

In conclusion, significant differences (−43.945 mas cy⁻¹ in longitude, 0.934 mas cy⁻¹ in obliquity) with respect to the values considered in IAU2006 precession (−0.960 mas cy⁻¹ in longitude, 0.340 mas cy⁻¹ in obliquity) have been found. This circumstance should be carefully considered when assessing a possible update of the IAU precession model (Liu & Capitaine 2017).

On the other hand, the indirect effects on nutations, due to the change of the dynamical ellipticity of some parts per million, must be taken into account for full consistency between the precession and nutation theories, required by the IAU Resolutions and also a main goal of the Global Geodetic Observing System (GGOS) of the International Association of Geodesy (IAG) and the IAU/IAG Joint Working Group on theory of Earth rotation and validation (Ferrándiz & Gross 2015; Transactions IAU 2015; IAG Travaux Reports 2017).

In turn, the obtained numerical results mean that much attention must be paid to the ocean tide models because of their relevant influence in precession and nutation. Our computation relies on the developments by to Williams & Boggs (2016) based on the ocean tide model FES2004 (IERS Conventions 2010, Sect. 6.3.2). Other models might be explored and their results given in a clear form adapted to the Love number formalism, in a similar way to Williams & Boggs (2016). The usefulness of that approach has been proven in this work.

Acknowledgments

The authors thank the anonymous referee and the Editor for their valuable comments and suggestions that helped to improve the manuscript. This work has been partially supported by the Spanish project AYA2016-79775-P (AEI/FEDER, UE).

References