On two formulations of polar motion and identification of its sources

Differences in interpretation may arise from differences in formulation of the equations of celestial mechanics. This paper focuses on the Liouville system of differential equations. In a"modern"presentation of the equations, variations in polar motion (PM) and variations in length of day (lod) are decoupled. Their source terms result from redistribution of masses and torques. Elasticity of the Earth, large earthquakes, or external forcing by the fluid envelopes have been successively invoked in the"modern"approach. In the"classical"presentation, PM is governed by the inclination of Earth's rotation pole and the derivative of its declination (close to the derivative of lod). These are coupled by the Liouville system and all"source"terms are astronomical. The"classical"approach also implies that there should be a link between the rotations and the torques exerted by the planets of the solar system. A sufficiently strong source of energy can be found. There is remarkable agreement between the sum of forces exerted by the four Jovian planets and components of Earth's polar motion. Since 1970, lod tends to decrease, but a recent acceleration of rotation velocity, that contradicts previous models, finds a simple explanation with the"classical"formalism. Analysis of lod finds nine components, all with physical sense: first comes a"trend", then pseudo-oscillations with periods 80 yr , 18.6 yr, 11 yr , 1 year and 0.5 yr, 27.54 days, 13.66 days, 13.63 days and 9.13 days. The longer periods, 1 yr, 11 yr, 18.6 yr and 80 yr are common to lod and PM. The Liouville system of equations in the"classical"analysis implies a strong link between the two: the straightening of the inclination of the axis of rotation should accompany a decrease in lod.Also, if this link is valid, all the components with extraterrestrial periods should be present in the series of sunspots and indeed they are.


Introduction
Over more than two centuries, scientists have attempted to measure and to explain the variations in the length of the day (lod), or the equivalent rotation velocity of Earth, and changes in the geographical location of the pole of rotation, that is the place where the rotation axis intersects the Earth surface. A thorough treatment is in the Treatise of celestial mechanics of Pierre-Simon de Laplace (1749Laplace ( -1827Laplace (1799)) where the great scientist derives the system of differential equations that fully describes the motions of the rotation axis of any celestial body, among others Earth. This system has come to be known as Liouville-Euler after mathematicians Leonhard Euler (1707-1783; the first to establish this system for pole motion) and Joseph Liouville (1809-1882; he complemented the the-ory of partial differential equations, showing in a theorem that the volume of phase space of a system is constant along trajectories of the system -in modern terms). The theory has been confirmed and elaborated on by a number of authors, Poincaré (1899) among them. Recent formulations are found in many papers and textbooks. In the present note, we focus on that of Lambeck (2005). In a first section, we recall the theoretical derivations of Laplace and Lambeck and show aspects in which they are formally identical, but also differences that can be significant and need to be explained. The main point is the identification of the sources (or excitation functions) of polar motion and length of day (lod). In a second section, we illustrate these differences by applying the theory to modern data of polar motion and lod. We discuss them and draw some conclusions in the final section. Lopes et al. (2021) and Courtillot et al. (2021) have recently recalled in some detail how Laplace (1799) derived the system of differential equations that was later to be named Liouville-Euler. We call "full" polar motion the vector consisting of the two spherical coordinates of the pole m 1 and m 2 (according to the X 1 OX 2 plane) and the third coordinate m 3 , linked to the length of the day (cf. Fig.1). The motion of the Earth's rotation axis (ω) can be seen as the combination of three Euler angles ω 1 , ω 2 and ω 3 . The rotation axis moves only very small distances from its mean position (at least over the past century of continuous measurements) and one can write:

Two formulations of the Liouville-Euler equations
where Ω (=7.292115*10 −5 rad/s) is the Earth's mean rotation velocity today computed on the last 3 decades. Applying the theorem of kinetic momentum to the rotation of a non rigid body and following Lambeck (2005), chapter 3, equations (1) lead to the set of Liouville -Euler equations (system 3.2.9 in Lambeck (2005)): where i = √ −1, m = m 1 + im 2 , σ r is the Euler frequency and f 3 are the so-called excitation functions (e.g. Lambeck (2005), chapter 4). C and A are respectively, the axial and equatorial moments of inertia of Earth (8.0365*10 3 7 kg.m 2 and 8.010*10 3 7 kg.m 2 , see Chen and Shen (2010)). In this derivation, the behavior of the pole position (m 1 ,m 2 ) and m 3 have been fully separated. Through σ r , (m 1 ,m 2 ) involve the (internal) terrestrial data C and A. Lambeck (2005), page 34, writes: "m 1 and m 2 are the components of the polar motion or wobble and Ω dω 3 dt is nearly the acceleration in diurnal rotation". The generally accepted reading (physical interpretation) of this formulation is that polar motion (m 1 , m 2 ) is linked to geophysical excitation such as atmospheric or oceanic circulation, lithospheric and mantle convection or electromagnetic coupling, and that the m 3 component is linked to astronomical phenomena such as tides. Lambeck (2005), page 36) concludes: " Equations (3.2.6) clearly separate the astronomical and geophysical problems." This is one reason for which on long time scales mechanical properties of the mantle are called upon. Some of these excitations can vary along with climate variations. Such is the case in the theory of isostasy (e.g. Peltier and Andrews (1976); Nakiboglu and Lambeck (1980)). It is also the reason why one calls upon the ephemerids of the Moon and Sun in order to compute Earth tides (e.g. Melchior (1966); Ray and Erofeeva (2014)) or to evaluate their influence on lod variations (e.g. Wahr et al. (1981); Le Mouël et al. (2019)).
In the classical problem of the motion of a Lagrange top (e.g. Pérez (1997)), the weight of the top plays the role of a perturbation of the full motion of the rotation axis. This is replaced by the astronomical torques of the Moon and Sun as excitations of the Earth's rotation axis. This formalism has allowed one to compute the period of equinoctial precession (26,000 yr). This is how one is led to the Milankovitch theory of climate variations (Milanković (1920)). The role of other bodies in the solar system, mainly the large and remote Jovian (ice) planets, must be taken into account (e.g. Laskar et al. (2004); Laskar et al. (2011)), for they excite responses in eccentricity at very long periods (e.g. 400 kyr). As is the case for the top, changes in Earth rotation perturb its revolution about the Sun.
Laplace's "classical" formulation of the theory is formally identical to the "modern" formulation recalled above (Lambeck (2005)). But it leads to a different interpretation (reading). Laplace (1799), page 74, system D, deduces from the fundamental law of dynamics a system of equations that can be shown to be identical to (2). It reads (Laplace (1799), page 74, system D): On the left side of (3a), (p,q,r) stand for the Euler angles (ω 1 , ω 2 and ω 3 ) of equations (1) and (p',q',r') = (Cp, Aq, Br), still with (A,B,C) being the Earth's moments of inertia. Let (x, y, z) be the coordinates of the Earth's center of gravity, (x',y',z') the coordinates of a mass element dm of Earth. Then: are classical expressions for the second order inertia tensor. θ and ψ are as defined is Fig.1. How can N, N' and N" be evaluated without knowing the motions of all particles dm in and out of the Earth? In a bold move, Laplace (1799), volume 1, book 5, page 305, paragraph 3, assumes that dN, dN' and dN" can be computed from the positions and motions of the celestial bodies that act on it.
Laplace starts with the Sun, its mass L and distance r 1 . (3a) becomes: This system is close to that proposed by Guinot (1976), page 530, system 1: Changing the directions of the axes and since B ∼ A, (3c) and (3d) are equivalent. But for Laplace, all terms on the right side are celestial (astronomical), whereas for Guinot the torques L t , M t and N t (not to be mistaken for N in 3a) can be external or internal to the Earth. For the Laplace formulation, one must take into account all planets that can produce effects one wants to account for. For instance, Laplace gives the full equations with the Sun and Moon included: (1 + λ).m.
The inclination θ of the rotation axis has the current value h in (4a). dψ dt is linked to the Earth's rotation, therefore to the lod. On the right side of (4a and 4b) are the ephemerids and masses of the Moon and Sun that enter the classical theory of gravitation (see Appendix A in Lopes et al. (2021) for more details). Length of day and polar inclination are clearly connected by equations (4a and 4b). Thus, Laplace reduces the problem to a system of two equations for the inclination and time derivative of the declination of the Earth's rotation axis. θ and dψ dt (and the norm that can be considered as a known constant) give the direction of the polar rotation axis and its variations. The time difference (in ms) between the theoretical and measured Earth rotation is proportional to ψ v , v being the rotation velocity (and the Earth's radius is a constant). Either ψ alone, or v alone, or both can vary. We assume the former, since the mean rotation rate apparently remains constant, as was already noted above, and equation (4b) implies studying the time derivative of declination of the rotation axis, thus studying a quantity that is linearly related to the derivative of lod. Stephenson and Morrison (1984) have compiled the reference data for lod from 700AD onwards. This has allowed Gross (2001) to build a monthly data set of lod from 1832 to 1997 (LUNAR97). We have combined it with the daily data provided by IERS (International Earth Rotation Service) since 1962 1 . Fig.(2) shows the two data sets, and their (smoother) trends (the red curve). In order to determine these trends, we have applied Singular Spectrum Analysis (SSA; see Golyandina and Zhigljavsky (2013); Lemmerling and Van Huffel (2001); Golub and Reinsch (1971)). The trends are the first, leading (in terms of pseudo-period and amplitude) components of the data series. The trends of the two series are smoothly continuous where they meet (1962).

Confronting the theory with the observations
The data for the m 1 and m 2 components of polar motion since 1846 are also available from IERS and are shown in Fig.3a. Their respective trends, extracted by SSA, are shown in Fig.3b. m 1 and m 2 are used to compute a global trend of m (= m 1 + im 2 ), called the Markowitz drift (Markowitz (1968)). This is displayed in Fig.4a as a thinner gray curve, and compared to the Morrison/IERS trend of lod (thicker black curve). These curves are in excellent agreement with previous determinations (by e.g. Stoyko (1968); Hulot et al. (1996)), though they are smoother due to SSA extraction. We recall that the Markowitz drift is one of the three main components of polar motion along with the Chandler free oscillation and the forced annual oscillation (e.g. Gibert et al. (1998); Zotov and Bizouard (2012); Lopes et al. (2021)). The magnitude of the Markowitz drift is on the same order as plate tectonic velocities, that originally made it quite difficult to detect.

Equations (4a) and (4b) imply an integrative link between
θ and dψ dt , that is polar motion and length of day. On Fig.4b, we show again the lod curve pictured in Fig.4a (thicker black curve) and the derivative of the Markowitz drift (thinner gray curve). Despite extraction by the SSA method, the lod curve is still affected by higher frequency variations. This may be due to the fact that the sampling rate jumps from monthly to daily in 1962 and to the presence of a derivative (that is a high-pass filter). We smooth further (moving average windows of 10 years) the lod curve (thicker black curve; Fig.4c) that is displayed with the derivative of the Markowitz drift already shown in Fig.4b (thinner gray curve).
The two curves of Fig.4c are quite similar, with a lag of about a decade (lod and polar motion being in quadrature). This is consistent with the integration link between the two. We have computed the relative phase and amplitude variations that would put the two curves in best agreement (normalized to take care of the different units). We have done this by applying the simulated annealing technique (Kirkpatrick et al. (1983)) in order to bring the pole motion curve of Fig.4c (gray curve) into superposition with lod (black curve).     .5 shows the results, the amplitude and phase shift respectively at the top and bottom. We do not display the results of this inversion on waveforms since they are practically superposable (the dynamic time warping algorithm calculates the modulations of phase and amplitude of the operator that transforms the black curve of Fig.4c into the gray curve; the end result is two curves that are identical to 99%). The scaling factor (or amplitude of the operator) required to transform the normalized polar motion to the normalized lod ranges between 0.6 and 1.6 and averages 1.03 ± 0.30: to first order, the amplitudes of the two geophysical quantities m and lod evolve in parallel. The larger differences occur between 1870 and 1890, and a century later between 1970 and 1990 (Fig.5 top). The largest "jump" in both lod (4 ms) and pole motion (8.10-3 "/yr) takes place between 1870 and 1900 (Fig.4c). The phase shift of that operator ranges between 6 and 16 years and averages 10.7 ± 3.0 yr. Note that during the period between 1920 and 1940, when the amplitude factor is close to a minimum (∼ 0.7), the Chandler free oscillation of polar motion suffers a well-known phase jump of π.

Discussion and concluding remarks
Passing from the "classical" system of equations (4a and 4b, Laplace (1799)) to the "modern" system of equations (2, e.g. Lambeck (2005)) is not complicated but is rather lengthy. There is no contradiction between the two formulations. However, each can be read with its own emphasis. The main difference between what Laplace and Poincaré write, and what Lambeck and Guinot write arises at the step of system (3a), leading either to (2) or to (4a and 4b). How does one account for system (3b), that is for the distribution and motions of masses inside or outside Earth, at any place and instant? Laplace solves the point with a hypothesis: if one does not have access to masses and accelerations, one can still determine the force budget. For Laplace, these forces are exclusively external, leading to system (3c), and then (4a and 4b) with the Sun and Moon included. As a consequence, internal and external mass motions only serve to dissipate the energy received by polar motion from all celestial bodies (see Lopes et al. (2021), Appendix 1). With Lambeck and Guinot's formulation (2), one has to resort to purely terrestrial forces to explain mantle motions and changes in climate. Emphasis is on the separation between the polar coordinates (polar motion) and the third coordinate (linked to lod), and on determining the excitation functions, that can be external or internal to the Earth. These are mathematical functions that "include all factors that perturb the rotation motion" (Lambeck (2005), page 36). These excitation functions are listed by Lambeck (2005), page 47: "the excitation functions consist of contributions from (i) redistribution of mass, (ii) relative motion of matter, and (iii) torques", the difference between (i) and (ii) being tenuous. For Guinot (system 3d) the torques can be one or more external and internal actions.
In short, in the "modern" presentation, polar motion and length of day are decoupled and correspond to different physical mechanisms, whereas in the "classical" presentation polar motion involves coupling of the inclination of the rotation pole and the derivative of its declination.
In a previous SSA analysis of lod using more than 50 years of IERS observations, Le Mouël et al. (2019) find nine components, all likely astronomical -thus external to the Earth ("trend", ∼ 80 yr, 18.6 yr, 11 yr, 1 year, 0.5 yr, 27.54 days, 13.66 days, 13.63 days, 9.13 days). The QBO at 2.36 yr is interpreted as a Sun-related oscillation. The lunar components at 13.63 and 13.66 days could contain a solar contribution. The longer periods, 1 yr, 11 yr, 18.6 yr and 80 yr (Markowitz,Fig.4c) are common to lod and polar motion. If the link between the two is as established by Laplace (1799), then the straightening of the inclination of the axis of rotation would accompany a decrease in lod. The Earth indeed straightens up (cf. Stoyko (1968) and Fig.3b). Also, if this link is valid, all the components with extraterrestrial periods should be present in the series of sunspots; and indeed they are (e.g. Courtillot et al. (2021)). As far as quasi annual components of sunspots are concerned, there is no exact 1 yr line, but two nearby lines with periods 0.93 and 1.05 yr (Le Mouël et al. (2020),  (2022)). Can a sufficiently strong source of energy be found? Indeed it has been known for some time (e.g. Dickman (1989); Chao et al. (1995); Varga et al. (2005)) that the components with luni-solar periods found above account for 95% of the total variance of the lod signal (Le Mouël et al. (2019)).
Another benefit of using the original Laplace approach is the determination of the period of the Euler free oscillation. In one of the first analyses of actual observations of polar motion, Chandler (1891a, 1891b) discovered the wobble that now bears his name. Chandler wrote that "the general result of a preliminary discussion is to show a revolution of the earth's pole in a period of 427 days." The observed period of this free oscillation is much larger than the theoretical value of 306 days, even more so now that the period has reached about 433 days (e.g. Zotov and Bizouard (2012); Lopes et al. (2017)).
In the more classical reading of Laplace, still considering an elastic (or plastic) Earth, we have seen that the equations (4a) and (4b) allow one to calculate the period of the Euler free oscillation. This oscillation actually varies with the pole inclination θ from 306 to 578 days. A transition to a double period (∼ 430 and (∼ 433 days:) has taken place at the time of the 1930 phase jump of the Chandler oscillation; this can be accounted for by variations in polar inclination (cf. Fig.4c). As is the case for a top, as recalled in section 2, the various rotations (precession, etc. . . ) vary with the inclination of the rotation axis.
Laplace's (1799) calculations and conclusions have been confirmed, prominently by Poincaré (1899). Based on these equations, a link between the rotations and the torques exerted by the planets of our solar system is expected. Indeed, we have shown elsewhere the influence of Jovian planets on the Sun (sunspots, cf. Courtillot et al. (2021)), and Earth (Lopes et al. (2021), Figure 11). In the latter paper, we show the similarity between the envelope of the Chandler oscillation and the ephemerids of Neptune. Figure 03 of Lopes et al. (2021), reproduced here as Fig.6, shows the remarkable agreement between the sum of forces exerted by the four Jovian planets and the m 1 component of polar motion.
As far as the length of day is concerned, we have seen in Fig.2 that since 1970 it tends to decrease, i.e. the Earth's velocity of rotation tends to increase. The same phenomena envisioned above (earthquakes, variations in the fluid envelopes,. . . ) are seen as potential causes of this trend (e.g Rochester (1984) 2021)). Also, the recent acceleration of rotation velocity, that contradicts previous models, may find a simple explanation with Laplace's formalism (e.g. Trofimov et al. (2021)). In closing, we emphasize that all these results are based on actual data and the validity of Laplace's astronomical hypothesis, not on a model.