Forced precession and nutation of Earth

This appendix gives a more precise treatment of the precession and nutation of the Earth's axis of rotation, forced by the joint gravitational influence of the Sun and the Moon, than that presented in Section 8.10.

Consider the approximately inertial Cartesian coordinate system, $x$, $y$, $z$, described in Section 8.10. Let

$\displaystyle \mbox{\boldmath$\omega$}$$\displaystyle = \omega\,(-\sin\theta\,\sin\phi,\,\sin\theta\,\cos\phi,\,\cos\theta)$ (F.1)

be the Earth's instantaneous angular velocity vector due to its daily rotation. Here, $\theta $ is the inclination angle of the Earth's rotation axis to the normal to the ecliptic plane, whereas $\phi $ is the angle subtended between the projection of the rotation axis onto the ecliptic plane and the $y$-axis. Note that $\theta $ and $\phi $ are Euler angles. (See Section 8.7.) Suppose that a body of mass $m$ is in a Keplerian orbit about the Earth. The coordinates of the body are thus (see Section 4.12)

$\displaystyle x$ $\displaystyle = r\left[\cos{\mit\Omega}\,\cos(\tilde{\omega}+f)-\sin{\mit\Omega}\,\sin(\tilde{\omega}+f)\,\cos I\right],$ (F.2)
$\displaystyle y$ $\displaystyle = r\left[\sin{\mit\Omega}\,\cos(\tilde{\omega}+f)+\cos{\mit\Omega}\,\sin(\tilde{\omega}+f)\,\cos I\right],$ (F.3)
$\displaystyle z$ $\displaystyle = r\,\sin(\tilde{\omega}+f)\,\sin I,$ (F.4)


$\displaystyle r = \frac{a\,(1-e^{\,2})}{1+e\,\cos f} = a\,(1-e\,\cos E).$ (F.5)

Here, $r=(x^{\,2}+y^{\,2}+z^{\,2})^{1/2}$, and $a$, $e$, $I$, $\tilde{\omega}$, ${\mit \Omega }$, $f$, and $E$ are the body's orbital major radius, eccentricity, inclination (to the ecliptic), argument of the perigee, longitude of the ascending node, true anomaly, and eccentric anomaly, respectively. (See Section 4.12.)

The potential energy of the Earth in the gravitational field of the orbiting body is

$\displaystyle U = - \frac{G\,m\,M}{r} + \frac{G\,m\,({\cal I}_\parallel-{\cal I}_\perp)}{r^{\,3}}\,P_2(\cos\gamma).$ (F.6)

Here, $M$, ${\cal I}_\parallel$, and ${\cal I}_\perp$ are the Earth's mass, moment of inertial about its rotation axis, and moment of inertia about an axis lying in its equatorial plane, respectively. (See Section 8.10.) Moreover,

$\displaystyle \cos\gamma = \frac{\mbox{\boldmath$\omega$}\cdot{\bf r}}{\vert\mbox{\boldmath$\omega$}\vert\,\vert{\bf r}\vert}$ $\displaystyle = \cos(\tilde{\omega}+f)\,\sin\theta\,\sin({\mit\Omega}-\phi)$    
  $\displaystyle \phantom{=}+ \sin(\tilde{\omega}+f)\left[\cos I\,\sin\theta\,\cos({\mit\Omega}-\phi)+\sin I\,\cos\theta\right],$ (F.7)

where ${\bf r}$ is the position vector of the orbiting body.

We are interested in the effect of the orbiting body on the Earth's axis of rotation on timescales that are much longer than the orbital period. We can concentrate on this effect, and filter out any relatively short-term oscillations, by averaging the potential energy function, $U$, over an orbital period. In other words, we replace $U$ by

$\displaystyle \overline{U} = (1-e^{\,2})^{-1/2}\oint\left(\frac{r}{a}\right)^2 U\,\frac{df}{2\pi}=\oint \frac{r}{a}\,U\,\frac{dE}{2\pi}.$ (F.8)

(See Section 10.5.) It follows that

$\displaystyle \overline{U}$ $\displaystyle = -\frac{G\,m\,M}{a} - \frac{G\,m\,({\cal I}_\parallel-{\cal I}_\...
... I}_\parallel\left\{-\left(1-\frac{3}{2}\,\sin^2 I\right)\cos(2\,\theta)\right.$    
  $\displaystyle \phantom{=}\left. + \sin(2\,I)\,\sin(2\,\theta)\,\cos ({\mit\Omega}-\phi)
-\sin^2 I\, \sin^2\theta\,\cos[2\,({\mit\Omega}-\phi)]\right\},$ (F.9)


$\displaystyle \alpha$ $\displaystyle = \frac{3}{8}\,\frac{\skew{3}\tilde{\epsilon}}{(1-e^{\,2})^{\,3/2}}\,\frac{G\,m}{a^{\,3}},$ (F.10)
$\displaystyle \skew{3}\tilde{\epsilon}$ $\displaystyle = \frac{{\cal I}_\parallel-{\cal I}_\perp}{{\cal I}_\parallel}.$ (F.11)

The Earth's orbit-averaged rotational Lagrangian is

$\displaystyle {\cal L} = \frac{1}{2}\left({\cal I}_\perp\,\skew{5}\dot{\theta}^...
{\cal I}_\parallel\,\omega^{\,3}\right)-\overline{U},$ (F.12)


$\displaystyle \omega =\cos\theta\,\skew{5}\dot{\phi}+\skew{5}\dot{\psi}.$ (F.13)

(See Section 8.10.) Here, $\psi$ is the Earth's third Euler angle. The two non-trivial equations of motion derived from the preceding Lagrangian are

$\displaystyle \frac{d}{dt}\left(\frac{\partial{\cal L}}{\partial\skew{5}\dot{\theta}}\right)-\frac{\partial{\cal L}}{\partial \theta}$ $\displaystyle = 0,$ (F.14)
$\displaystyle \frac{d}{dt}\left(\frac{\partial{\cal L}}{\partial\skew{5}\dot{\phi}}\right)- \frac{\partial{\cal L}}{\partial \phi}$ $\displaystyle = 0.$ (F.15)

(The trivial equation merely confirms that the Earth's axial rotation rate, $\omega$, is a constant.) Assuming that $\partial/\partial t\ll \omega$, we obtain

$\displaystyle \omega\,\sin\theta\,\,\skew{5}\dot{\phi}$ $\displaystyle \simeq -2\,\alpha\,[1-(3/2)\,\sin^2 I]\,\sin(2\,\theta)-2\,\alpha\,\sin(2\,I)\,\cos(2\,\theta)\,\cos({\mit\Omega}-\phi)$    
  $\displaystyle \phantom{=}+ \alpha\,\sin^2 I\,\sin(2\,\theta)\,\cos[2\,({\mit\Omega}-\phi)],$ (F.16)
$\displaystyle \omega\,\sin\theta\,\,\skew{5}\dot{\theta}$ $\displaystyle \simeq \alpha\,\sin(2\,I)\,\sin(2\,\theta)\,\sin({\mit\Omega}-\phi)-2\,\alpha\,\sin^2 I\,\sin^2\theta\,\sin[2\,({\mit\Omega}-\phi)].$ (F.17)

Finally, assuming that $\phi =0$ at $t=0$, and also that

$\displaystyle {\mit\Omega}(t) = {\mit\Omega}_0+\skew{5}\dot{\mit\Omega}\,\,t,$ (F.18)

the lowest-order (in $\skew{3}\tilde{\epsilon}$) solutions to Equations (F.16) and (F.17) are written

$\displaystyle \phi$ $\displaystyle = -{\mit\Omega}_\phi\,t +\delta\phi_1\left(\sin{\mit\Omega}-\sin{...
...ght) - \delta\phi_2
\left[\sin(2\,{\mit\Omega})-\sin(2\,{\mit\Omega}_0)\right],$ (F.19)
$\displaystyle \theta$ $\displaystyle = \theta_0 + \delta\theta_1\,\cos{\mit\Omega}-\delta\theta_2\,\cos(2\,{\mit\Omega}),$ (F.20)


$\displaystyle {\mit\Omega}_\phi$ $\displaystyle =\frac{3}{2}\,\frac{\skew{3}\tilde{\epsilon}}{\omega\,(1-e^{\,2})^{\,3/2}}\,\frac{G\,m}{a^{\,3}}\left(1-\frac{3}{2}\,\sin^2 I\right)\cos\theta_0,$ (F.21)
$\displaystyle \delta\theta_1$ $\displaystyle =-\frac{3}{4}\,\frac{\skew{3}\tilde{\epsilon}}{\omega\,\skew{5}\d...
...t\Omega}\,(1-e^{\,2})^{\,3/2}}\,\frac{G\,m}{a^{\,3}}\,\sin(2\,I)\,\cos\theta_0,$ (F.22)
$\displaystyle \delta\theta_2$ $\displaystyle = -\frac{3}{8}\,\frac{\skew{3}\tilde{\epsilon}}{\omega\,\skew{5}\dot{\mit\Omega}\,(1-e^{\,2})^{\,3/2}}\,\frac{G\,m}{a^{\,3}}\sin^2 I\,\sin\theta_0,$ (F.23)
$\displaystyle \delta\phi_1$ $\displaystyle = -\frac{3}{4}\,\frac{\skew{3}\tilde{\epsilon}}{\omega\,\skew{5}\...
.../2}}\,\frac{G\,m}{a^{\,3}}\,\sin(2\,I)\,\frac{\cos(2\,\theta_0)}{\sin\theta_0},$ (F.24)
$\displaystyle \delta\theta_2$ $\displaystyle =- \frac{3}{8}\,\frac{\skew{3}\tilde{\epsilon}}{\omega\,\skew{5}\dot{\mit\Omega}\,(1-e^{\,2})^{\,3/2}}\,\frac{G\,m}{a^{\,3}}\sin^2 I\,\cos\theta_0.$ (F.25)

Here, ${\mit\Omega}_\phi$ is the precession rate of the Earth's axis of rotation (a positive rate corresponds to retrograde precession) induced by the orbiting body. Furthermore, $\delta\theta_1$ and $\delta \theta_2$ are the amplitudes of the nutation in the polar angle, $\theta $, induced by the body. The first amplitude corresponds to nutation with a period equal to the body's nodal precession period, whereas the second amplitude corresponds to nutation with half this period. Finally, $\delta\phi_1$ and $\delta\phi_2$ are the corresponding amplitudes of the nutation in the azimuthal angle, $\phi $. Note that, in the preceding analysis, there is an implicit assumption that $\vert\delta\theta_1\vert$, $\vert\delta\theta_2\vert$, $\vert\delta\phi_1\vert$, $\vert\delta\phi_2\vert\ll 1$, and, hence, that $\vert\skew{5}\dot{\mit\Omega}\vert\gg {\mit\Omega}_\phi$.

The effects of the Sun and the Moon on the Earth's axis of rotation are additive (because the Sun and Moon's gravitational fields are additive). For the Sun, we can write

$\displaystyle \frac{G\,m}{a^{\,3}}$ $\displaystyle \simeq n_s^{\,2},$ (F.26)
$\displaystyle e$ $\displaystyle =e_s,$ (F.27)
$\displaystyle I$ $\displaystyle = I_s=0,$ (F.28)

where $n_s$ and $e_s$ are the mean orbital angular velocity, and eccentricity, of the Sun's apparent orbit about the Earth, respectively. For the Moon, we can write

$\displaystyle \frac{G\,m}{a^{\,3}}$ $\displaystyle = \mu_m\,n_m^{\,2},$ (F.29)
$\displaystyle e$ $\displaystyle =e_m,$ (F.30)
$\displaystyle I$ $\displaystyle = I_m,$ (F.31)
$\displaystyle \skew{5}\dot{\mit\Omega}$ $\displaystyle = -\vert\skew{5}\dot{\mit\Omega}_m\vert,$ (F.32)

where $n_m$, $e_m$, $I_m$, and $\skew{5}\dot{\mit\Omega}_m$ are the mean orbital angular velocity, eccentricity, inclination (to the ecliptic), and (negative) nodal precession rate, of the Moon's orbit about the Earth, respectively. Furthermore,

$\displaystyle \mu_m = \frac{M_m}{M_E+M_m},$ (F.33)

where $M_m$ and $M_E$ are the masses of the Moon and Earth, respectively.

The net luni-solar precession rate of the Earth's rotation axis (which is the sum of the precession rates due to the Sun and the Moon) is written

$\displaystyle {\mit\Omega}_\phi = \frac{3}{2}\,\frac{\skew{3}\tilde{\epsilon}\,...
\,\left(1-\frac{3}{2}\,\sin^2 I_m\right)\right]$ (F.34)

Now, $\omega/n_s = T_s({\rm day})$, $n_s/n_m = T_m({\rm yr})$, and $n_s/{\mit\Omega}_\phi = T_\phi({\rm yr})$, where $T_s({\rm day})$ is the number of stellar (sidereal) days in a sidereal year, $T_m({\rm yr})$ is the length of a sidereal month in units of sidereal years, and $T_\phi({\rm yr})$ is the precession period of the Earth's axis of rotation in sidereal years. Thus,

$\displaystyle T_\phi({\rm yr}) = \frac{2}{3}\,\frac{T_s({\rm day})}{\skew{3}\ti...
\frac{[1-(3/2)\,\sin^2 I_m]}{(1-e_m^{\,2})^{\,3/2}}\right\}^{-1}.$ (F.35)

Using the following observed values, $T_s({\rm day})=366.26$, $T_m({\rm yr})=0.07480$, $e_s=0.016711$, $e_m=0.05488$, $I_m=5.16^\circ$, $\mu_m=0.01215$, and $\theta_0=23.44^\circ$, (Yoder 1995) we obtain

$\displaystyle T_\phi({\rm yr}) = \frac{84.34}{\skew{3}\tilde{\epsilon}}.$ (F.36)

In fact, the observed precession period of the Earth's axis of rotation is $T_\phi({\rm yr})=25,772\,{\rm yr}$ (Yoder 1995). This suggests that the Earth's dynamical ellipticity, $\skew{3}\tilde{\epsilon}$, which cannot be measured directly, takes the value

$\displaystyle \skew{3}\tilde{\epsilon} = 0.003273.$ (F.37)

(In fact, this is how the Earth's dynamical ellipticity is most accurately determined.) The previous value is less than the value, $\skew{3}\tilde{\epsilon}=0.00335$, used in Section 8.10, which was derived from the Earth's observed flattening on the simplistic assumption that it is a homogeneous body. [See Equation (8.95).] We can write

$\displaystyle \frac{{\cal I}_\parallel}{M\,R^{\,2}}=\frac{J_2}{\skew{3}\tilde{\epsilon}},$ (F.38)

where $R$ is the Earth's equatorial radius, and $J_2= 1.083\times 10^{-3}$ is its measured dimensionless gravitational quadrupole moment. (See Section 10.5.) Thus, we deduce that

$\displaystyle {\cal I}_\parallel = 0.331\,M\,R^{\,2},$ (F.39)

which confirms that the Earth's internal mass distribution is not homogeneous (because, if it were then ${\cal I}_\parallel/M\,R^{\,2}$ would take the value $0.4$), but is instead centrally concentrated. In fact, the previous value for ${\cal I}_\parallel/(M\,R^{\,2})$ is very similar to that predicted by the so-called Darwin-Radau equation. (See Appendix D.) This suggests that the Earth's response to its centrifugal potential is essentially fluid-like.

Now, $n_s/\vert\skew{5}\dot{\mit\Omega}_m\vert=T_n({\rm yr})$, where $T_n({\rm yr})=18.61$ is the Moon's observed nodal precession period in sidereal years (Yoder 1965). Thus, making use of the deduced value, given in Equation (F.37), for the Earth's dynamical ellipticity, we find that the nutation amplitudes induced by the Moon are

$\displaystyle \delta\theta_1$ $\displaystyle = \frac{3}{4}\,\frac{T_n({\rm yr})}{T_s({\rm day})\,[T_m({\rm yr}...
  $\displaystyle \phantom{=} = 9.2'',$ (F.40)
$\displaystyle \delta\theta_2$ $\displaystyle = \frac{3}{8}\,\frac{T_n({\rm yr})}{T_s({\rm day})\,[T_m({\rm yr}...
...kew{3}\tilde{\epsilon}\,\mu_m}{(1-e_m^{\,2})^{\,3/2}}\,\sin^2 I_m\,\sin\theta_0$    
  $\displaystyle \phantom{=}=0.090'',$ (F.41)
$\displaystyle \delta\phi_1$ $\displaystyle = \frac{3}{4}\,\frac{T_n({\rm yr})}{T_s({\rm day})\,[T_m({\rm yr}...
  $\displaystyle \phantom{=}= 17.2'',$ (F.42)
$\displaystyle \delta\phi_2$ $\displaystyle = \frac{3}{8}\,\frac{T_n({\rm yr})}{T_s({\rm day})\,[T_m({\rm yr}...
...kew{3}\tilde{\epsilon}\,\mu_m}{(1-e_m^{\,2})^{\,3/2}}\,\sin^2 I_m\,\cos\theta_0$    
  $\displaystyle \phantom{=}= 0.21''.$ (F.43)

(The corresponding nutation amplitudes included by the Sun are all zero, because $I_s=0^\circ$.) The observed nutation amplitudes are $\delta\theta_1=9.2''$, $\delta\theta_2 = 0.090''$, $\delta\phi_1=17.2''$, and $\delta\phi_2 = 0.21''$ (Meeus 2005). The fact that the calculated nutation amplitudes agree so well with the observed amplitudes validates our general approach.