On the Slow Drift of Solstices: Milankovic Cycles and Mean Global Temperature

The Earth's revolution is modified by changes in inclination of its rotation axis. Despite the fact that the gravity field is central, the Earth's trajectory is not closed and the equinoxes drift. Milankovic (1920) argued that the shortest precession period of solstices is 20,7kyr: the Summer solstice in one hemisphere takes place alternately every 11kyr at perihelion and at aphelion. We have submitted the time series for the Earth's pole of rotation, global mean surface temperature and ephemeris to iterative Singular Spectrum Analysis. iSSA extracts from each a trend, a 1yr and a 60yr component. Both the apparent drift of solstices of Earth around the Sun and the global mean temperature exhibit a strong 60yr oscillation. The"fixed dates"of solstices actually drift. Comparing the time evolution of the Winter and Summer solstices positions of the rotation pole and the first iSSA component (trend) of the temperature allows one to recognize some common features. A basic equation from Milankovic links the derivative of heat received at a given location on Earth to solar insolation, known functions of the location coordinates, solar declination and hour angle, with an inverse square dependence on the Sun-Earth distance. We have translated the drift of solstices as a function of distance to the Sun into the geometrical insolation theory of Milankovic. Shifting the inverse square of the 60yr iSSA drift of solstices by 15 years with respect to the first derivative of the 60yr iSSA trend of temperature, that is exactly a quadrature in time, puts the two curves in quasi-exact superimposition. The probability of a chance coincidence appears very low. Correlation does not imply causality when there is no accompanying model. Here Milankovic's equation can be considered as a model that is widely accepted. This paper identifies a case of agreement between observations and a mathematical formulation.


Introduction
For over a century and a half, geological evidence has been put forward to argue that Earth had undergone cyclical changes in climate (Agassiz 1837). Adhémar (1860) described the long periodicities associated with glacial and interglacial cycles that were explained sixty years later by the mathematical theory of climate due to Milanković (1920). A geometrical derivation shows that the insolation W received at a location with coordinates ϕ (latitude) and ψ (longitude) is given by the (fundamental, yet simple) equation (eq. 20, page 15, Milanković (1920)): dW dt = I 0 ρ 2 [sin ϕ sin δ + cos ϕ cos δ cos(ω + ψ)] where ρ is the Sun-Earth distance, δ the Sun's declination and ω its hour angle. It is generally considered that the shortest Milanković cycle is the 19kyr precession cycle (eg. Laskar et al. (2004); Lopes et al. (2021a)). The Earth's eccentricity being quite small (e = 0.016), it is legitimate at first to consider that, over such long durations ρ and δ are constant. Thus, Milanković (1920) assumes that all climate variations arise because of movements of the rotation axis (ie. of the pole).
To first order, the rotation axis is perturbed by the conjugate effects of the Moon and Sun (eg. Laplace (1799); Poincaré (1899)), resulting for instance in the luni-solar tides (eg. Ray et Erofeeva (2014)). Over longer periods, the Jovian planets are the main source of perturbations of the rotation axis (cf. Laskar et al. (1993Laskar et al. ( , 2004), giving rise to variations in precession, obliquity and eccentricity.
Another source of perturbation of the Earth's rotation axis is the precession of equinoxes, associated with Kepler's second law (conservation of areas). There are four sources of precession of Earth's equinoxes, that is rotation of the whole orbit (revolution).
1) The first is associated with Kepler's laws. In the case of a central field and an elliptical orbit, for the orbit to be closed it is necessary and sufficient that the orbit's angular change after n revolutions be of the form ∆ϕ = 2πm/n, where m is the number of full revolutions necessary for the planet to recover its initial position. There are only two central fields in which ∆ϕ is a rational fraction of 2π, ensuring closed orbits, that is fields in r 2 and 1/r the latter being the case of our solar system (cf. Landau et Lifchitz (1964)).
2) The second involves the joint effects of the Moon and Sun. Let us quote d'Alembert (1749, page 14): "Enfin, l'inclinaison de l'axe terrestre au plan de l'ecliptique doit modifier aussi l'action du Soleil; car selon que cet axe sera différemment incliné, il feraà chaque point de l'ecliptique un angle différent avec la ligne qui joint les centres de la Terre et du Soleil; par conséquent la quantité et la loi de l'action du Soleil, dépend de l'inclinaison de l'axe, et c'est aussi ce que l'analyse apprend." As clearly stated by Laplace (1799) and later Poincaré (1899), the planet's revolution is modified by changes in inclination of the rotation axis, principally due to the joint actions of the Moon (for 2/3rds) and Sun (for 1/3rd). One is no longer in case 1: despite the fact that the field is a central one, the trajectory is not closed anymore. d'Alembert (1749, page 52) estimates the drift of the equinoxes to be about 50 " per year, that is a precession period of 25920 yr.
The two other processes are relativistic.
3) The Sun containing 99% of the total mass of the solar system, Schwarzschild (1916) shows that the planet's revolution about the Sun produces an additional precession of about 3.8" per century, or a period of some 33 million years. 4) Because the Sun is actually a huge rotating mass, there is an additional relativistic component of precession, with a period on the order of 5.8 million years Lense et Thirring (1918).
For Newton, planetary bodies attract (or repel) the oceans (and atmospheres) and this re-organization of masses modifies the rotation axis. For d'Alembert, Lagrange, Laplace and Poincaré, Changes in polar motion and revolution are coupled and involve the luni-solar torques, as in a top. Re-organization of Earth's fluid envelopes (eg. tides) follows.
As Milanković (1920) writes, precession of the equinoxes is actually due to the joint attraction of the Moon and Sun on the Earth's equatorial bulge and its period is in theory 26,000 yr. Because the Earth itself rotates, areolar velocities vary between perihelion and aphelion (Kepler's second law); because of centrifugal forces, a precession with a period of 19,000 yr appears. So the precession of equinoxes undergoes a double periodicity, with a mean of 22,500 yr (half period 11,250 yr). Indeed, Milankovic (1920, p.221) writes that the first precession of perihelion (for us solstices) is 20,700 Julian years and that the consequence of this precession is that the Summer solstice in one hemisphere (when that hemisphere receives maximum insolation) takes place alternately every 11,000 yr at perihelion (thus a warmer Summer) and aphelion (thus a cooler Summer). The difference in insolation (energy received by Earth) between maximum and minimum is a function of eccentricity.
The 26 kyr period of precession has first been determined in the frame of Newtonian physics by d'Alembert (1749). It is rather close to the first precession cycle of 19 kyr in Milanković's (1920) theory. When Milanković makes the assumption that the planetary distances to the Sun and the solar ephemerids are constant, he can estimate climate maxima but not their smooth transitions between equinoxes and solstices. Today, we have access to observations that allow one to drop the hypotheses that ρ, δ and ω in equation 1 are constant. Thus, we can evaluate the consequences of changes in the position of the rotation axis on, for instance, atmospheric temperature, that is the main parameter in Milanković's (1920) theory of climate.
2. The data: temperature, pole motion, and solar ephemerids 2-1: Mean global temperatures: We have used the data series maintained by the Hadley Center for Climate Prediction and Research under the name HadCrut. In order to have an idea of the reliability of the data, we have selected five successive sets of HadCrut data: HadCrutv * , 1870-2000Rayner et al. (1996HadCrut2 † , 1856-2006Rayner et al. (2003HadCrut3 ‡ , 1850Brohan et al. (2006; HadCrut4 § , 1850-2021 Osborn et Jones (2014) and HadCrut5 ¶ 1850-2022 Osborn et al. (2021). Figure 1a shows all the data, and Figure 1b their Fourier transforms. There are rather significant differences between the data series, for instance between 1940 and 2020 in HadCrut3 (yellow curve) vs HadCrut5 (blue curve). Differences become larger after 1950, to the point that HadCrut3 has a plateau after 2000 when HadCrut5 grows linearly since 1960. We have already worked on these data sets Le Mouël et al. (2020) and pointed out these differences Courtillot et al. (2013), Figure 4. These differences of course result in differences in the Fourier spectra of Figure 1. As a result, the dominant spectral peak shifts from 60 to 80 yr, a topic discussed in several papers (Mazzarella et Scafetta 2012;Courtillot et al. 2013;Gervais 2016;Veretenenko et Ogurtsov 2019;Scafetta et al. 2020).  We do not present a figure with the raw data: the Earth orbit's eccentricity is so small that an annual oscillation since 1846 would transform into an unreadable quasi-sinus with 176 oscillations on the 15cm (or so) width of the figure, that is 14 oscillations per cm.

2-3: Rotation pole and length of day:
The motions of the rotation pole and variations in rotation velocity are available at the International Earth Rotation Reference Systems Service ** (IERS). They consist in the couple of coordinates (m 1 , m 2 ) of motion of the rotation pole PM (see Lambeck (2005); Lopes et al. (2021b)) and the length of day lod (eg. Le Mouël et al. (2019a)). We have selected data set EOP C01 IAU19801. Figures 2a and 2b respectively show the evolution of the couple (m 1 , m 2 ) since 1946 and of lod since 1962. We have used the semi-annual lod data provided by Stephenson et Morrison (1984) for the period 1832-1997, combined with the IERS data, resulting in the mean curve between 1832 and 2022 shown in Figure 2c (Lopes et al. 2022  (c) Black curve: lod semi-annual data since 1832 (from Stephenson et Morrison (1984)); Gray: daily lod data since 1962 from IERS; Red curve: The median of the trend for these two combined data sets.

Extraction and Analysis of the trends and annual oscillations
As in a number of previous papers eg. Le Mouël et al. (2019b, 2020; Lopes et al. (2021b), we have submitted time series to iterative Singular Spectrum Analysis (iSSA) and we will now do the same for the rotation, temperature and ephemeris time series presented in the previous section. We refer the reader to these papers and to Golyandina et Zhigljavsky (2013) for the SSA method, to Lemmerling et Van Huffel (2001) for properties of the Hankel and Toeplitz matrices that it uses, and to Golub et Reinsch (1998) for the singular value decomposition algorithm SVD). Figure 3a shows the first oscillatory iSSA component of lod that has a period of 1 yr. In order to check the quality of iSSA, we have compared some of the iSSA results with those obtained with the method of continuous wavelets, that is widely used in the literature. Figure 3b shows the scalogram (ie. the continuous wavelet transform) of lod since 1962 (Figure 2b). We have selected a Morse wavelet (cf. Olhede et Walden (2002); Lilly et Olhede (2012)), an analytic wavelet that is well adapted to time series with variable spectra. Between the dashed red lines are the ordinates of the wavelet coefficients associated with the 1 yr oscillation. Given the property of information redundancy of the wavelet kernels, one can extract the ridge in the scales covered by the 1 yr periodicity and reconstruct the signal Gibert et al. (1998). The (wavelet) reconstructed annual component of lod is shown as the black curve in Figure 03a. The comparison is good but not perfect: the modulation of the iSSA curve (in red) is smoother than that of the wavelet reconstruction. The similarity of the two curves in Figure 3a together give us confidence that iSSA can correctly extract the annual components of the three time series introduced in Section 2.
A difficulty seen clearly in Figure 3b is that the signal's energy can spread and diffuse over several scales. The reconstruction could be optimized by applying correction methods, such as the reassignment method (cf. Auger et Flandrin (1995)), but this would require an additional step that is not necessary in the present study.

The Lissajous diagrams
The main points in space that can be used to monitor the precession of an elliptical orbit are the solstices and equinoxes. The fixed dates of their occurrences are December 21 for the Winter solstice, March 21 for the Spring equinox, June 21 for the Summer solstice, and September 23 for the Fall equinox. Since the legal and astronomical calendars are not exactly the same, this entails an error on their positions that can be estimated. The variations of these positions are actually very small: the perimeter of the Earth's orbit being close to 6.28 astronomical units (a.u.), the December 21 positions are at a distance of 0.98412 ± 3.46 10 −5 a.u., that is an error of 5.5 10 −4 between 1844 and 2022. The error is 5.8 10-4 for the other solstice and close to 7.9 10 −4 and 8.7 10 −4 for the two equinoxes. Figure 5a shows the evolution of the trajectories of Figure  4a as a function of Sun-Earth distance (ephemerids), that is in other words as a function of time. We call these by analogy to Lissajous orbits in astronomy "Lissajous diagrams". The Lissajous diagram of (m 1 , m 2 ) is shown in two perspectives in order to gain some insight on its topology. The locations of the four remarkable points (equinoxes and solstices) are shown in four different colors (the same are used throughout of the paper). We can see that the closer the Earth is to the Sun, the more the rotation axis straightens; the farther it is, the larger the amplitudes of motions and the flatter the rotation axis, i.e. the larger its declination (see also Figure 5b). Note that the "fixed  dates" of equinoxes and solstices appear to drift as a function of time.
The conservation of momentum of the orbiting planet implies that its areolar velocity is constant (the areolar velocity is the area spanned by the vector radius -the Sun to Earth vector -per unit time). As a consequence of this "law of areas", the orbital velocity varies from a maximum of 30.29 km/s at perihelion to a minimum of 29.29 km/s at aphelion. As explicitly stated by Lagrange (1781Lagrange ( , 1782; Laplace (1799) and Poincaré (1899)  In order to emphasize the relative amplitudes of the drift and the butterfly-like shape of the diagram, we have actually multiplied (m 1 , m 2 ) by the centered value d * S E = d S E -mean(d S E ). This also makes it clear that the drift of solstices is larger than that of equinoxes. Figure 6a shows the Lissajous diagram equivalent to that in Figure 05a for the couple (m * 1 ,m * 2 ). The Lissajous diagrams of Figure 6a represent the geodesic evolution of Kepler areas.  One sees clearly in Figure 6a that polar motion reaches a minimum at the equinoxes (red and orange dots) when solar attraction is the weakest, and a maximum at the solstices (green and blue dots) when the Sun, Earth and focus of the ellipse are aligned. We shall see below that the same applies to lod.
These results actually require another small correction, in relation with Kepler's second law. The conservation of momentum of the orbiting planet implies that its areolar velocity is constant (the areolar velocity is the area spanned by the vector radius -the Sun to Earth vector -per unit time). As a consequence of this "law of areas", the orbital velocity varies from a maximum of 30.29 km/s at perihelion to a minimum of 29.29 km/s at aphelion. We introduce new more physical variables by multiplying the polar motion coordinates by the Sun-Earth distance d S E : We display the drift of the four reference points (solstices and equinoxes) in the (m * 1 , m * 2 .) plane in Figure 6b. The time evolution of the separate coordinates m * 1 and m * 2 . for the solstices are shown in Figure 7a and 7b.      Milanković (1920) knew that eccentricity, precession and obliquity evolve slowly in time leading to a transition from a warmer to a cooler climate every 11 kyr (the period of precession of equinoxes being the shortest). Figures 5 to 7 show that one cannot consider that either the Sun-Earth distance or the hour angle, or the declination, or the daily variation are constant, contrary to what was done by Milanković (1920). Both polar motions and length of day are affected by their position on the elliptical orbit, on rather short periods on the order of less than 10 yr up to 1 century or more. The data shown in the present paper demonstrate that parameters in equation (1) evolve in time over these shorter time scales. They may therefore imply some forcing of climate on these same time scales. Figure 1a illustrates the differences between five data sets that were supposed when they were compiled to represent the same physical (or quasi-physical) datum, that is mean global surface temperature anomaly. Figure 08a displays the iSSA trends (component 1) of the five HadCrut temperature data sets introduced in section 2-1. The median of the five curves is shown as an inset in Figure 8a. Figures 8b and 8c represent the two other major iSSA components of global mean temperature anomaly, the annual and 60 yr oscillations. The frequencies of the five series are consistent except for the modulated amplitudes of the annual components (Figure 8b), which is puzzling.

Discussion
In Figure 9, we superimpose the time evolution of the Winter (blue) and Summer (green) solstices for component m * 1 from Figure 7a and the first iSSA component (trend) of the median of the five HadCrut Center temperature series (red) from Figure 8a (note: We have simplified the discussion by using only the m1 polar motion component, the one most clearly linked to the pole's inclination). With some familiarity with global temperature curves, one can recognize some common features with the evolution of the solstices Le Mouël et al. (2020). Based on the annual oscillations of the full polar motion (m 1 , m 2 and lod), we have seen that the celestial positions of solstices moved significantly in the past 180 years. To our knowledge, this is the first time one evidences in the observational data what d'Alembert (1749) called the apparent precession of solstices (that is seen from Earth) and that Milankovic (1920, part 2, chapter 2) worked on but on much longer time periods of millions of yearsWe recalled in the introduction the basic equation from Milankovic's thesis (1920) that links time variations of heat received at a given location on Earth to solar insolation, known functions of the location coordinates, solar declination and hour angle, with an inverse square dependence on the Sun-Earth distance. We can let W play the role of the heat/energy term in equation (1). The goal is to translate the drift of solstices as a function of distance to the Sun into the geometrical insolation theory of Milanković (1920).

Conclusion
The Earth's revolution is modified by changes in inclination of its rotation axis, principally due to the actions of the Moon and Sun. Despite the fact that the gravity field is central, the Earth's trajectory is not closed and the equinoxes drift (by a little less than one minute of arc per year, that is a precession period of some 26 kyr -d'Alembert (1749)). For d'Alembert, Lagrange, Laplace and Poincaré who based their reasoning on the action of torques, as in a top, changes in polar motion and revolution are coupled (through the Liouville-Euler system of equations). Re-organization of Earth's fluid envelopes follows. For Newton, who gives a central role to the inertia tensor, planetary bodies attract (or repel) the oceans (and atmospheres) and this re-organization of masses modifies the rotation axis. Milanković (1920, p.221) argued that the shortest precession period of perihelion (for us solstices) is 20,700 Julian years and that the consequence of this precession is that the Summer solstice in one hemisphere takes place alternately every 11 kyr at perihelion (a warmer Summer) and at aphelion (a cooler Summer). The difference in insolation between maximum and minimum is a function of eccentricity. Milanković assumed that the planetary distances to the Sun and the solar ephemerids are constant. There are now observations that allow one to drop the hypothesis that Sun-Earth distance, the Sun's declination and hourly angle (equation 1) are constant. Thus, we can evaluate whether changes in the position of the rotation axis affect, for instance, atmospheric temperature, that is the main parameter in Milanković's (1920) theory of climate.
(a) The first iSSA components (trends) of the five HadCrut Center temperature series introduced in section 2-1, and their median shown in black among the 5 trends and alone above as an inset.
(b) Annual iSSA components of the same series as in Figure 8a.  Both the apparent drift of solstices of Earth around the Sun and the global mean temperature exhibit a strong 60yr oscillation. In the present paper, we have confirmed the finding of a strong iSSA component with 60yr period in global temperatures and in the drift of solstices (Figure 8c and 10), hence the rotation axis, that has already been encountered in Lau et Weng (1995)  In pursuit of this goal, we obtained the Sun's ephemerids from 1846 to the present from Institut de Mécanique Céleste et du Calcul des Ephémérides (data set EOP C01 IAU1980). The motions of the rotation pole and variations in rotation velocity were taken from the International Earth Rotation Reference Systems Service. They consist in the couple of coordinates (m 1 , m 2 ) of motion of the rotation pole PM and the length of day lod. We used the semi-annual lod data provided by Stephenson et Morrison (1984) for the period 1832-1997, combined with the IERS data, resulting in a mean curve between 1832 and 2022 ( Figure 2c). For the mean global temperatures, we used five HadCrut data series from the Hadley Center for Climate Prediction and Research in order to estimate the consistence of the data, There are rather significant differences between the five data series, resulting in differences in the dominant spectral peak that shifts from 60 to 80 yr.
We have submitted the time series for the rotation, temperature and ephemeris to iterative Singular Spectrum Analysis (iSSA), a method which we have used extensively in a number of recent studies (Le Mouël et al. (2019aMouël et al. ( ,b, 2020). We have checked that the results obtained with iSSA match those obtained with the more common method of continuous wavelets that is widely used in the literature.
The main points in space that can be used to monitor the precession of an elliptical orbit are the solstices and equinoxes. The legal and astronomical calendars not being exactly the same, this entails a very small error on their positions (relative variations of their positions between 5 and 9 10-4 ) between 1844 and 2022. Figure 5a shows the evolution of the locations of the equinoxes and solstices trajectories as a function of Sun-Earth distance (we call these "Lissajous diagrams"). The closer the Earth is to the Sun, the more the rotation axis straightens; the farther it is, the larger the amplitudes of motions and the flatter the rotation axis. The "fixed dates" of equinoxes and solstices actually drift as a function of time. In order to emphasize the relative amplitudes of the drift and the butterfly-like shape of the Lissajous diagram, we have actually multiplied (m1, m2) by the centered value d * S E = d S E -mean(d S E ), yielding (m * 1 , m * 2 ). This also makes it clear that the drift of solstices is larger than that of equinoxes. Polar motion reaches a minimum at the equinoxes when solar attraction is the weakest, and a maximum at the solstices when the Sun, Earth and focus of the ellipse are aligned.
Both polar motions and length of day are affected by their position on the elliptical orbit, on rather short periods on the order of less than 10 yr up to 1 to a few centuries. Parameters in equation (1) evolve over these time scales. Some forcing of climate on these same time scales may be expected. Despite differences between the five HadCrut temperature data sets, the median of the first iSSA component (trend) appears to be representative (Figure 8a). The two other major iSSA components of the five global mean temperature anomaly series, the annual and 60 yr oscillations, are consistent; the modulated amplitudes of the annual components do not agree as well (Figure 8b).
We have superimposed the time evolution of the Winter and Summer solstices for component m * 1 and the first iSSA component (trend) of the median of the five HadCrut Center temperature series (Figure 9). Based on the annual oscillations of the full polar motion (m1, m2 and lod), we have shown that the celestial positions of solstices moved significantly in the past 180 years. To our knowledge, this is the first time one evidences in the observational data what d'Alembert (1749) called the apparent precession of solstices (that is seen from Earth) and that Milankovic (1920, part 2, chapter 2) worked on, but on much longer time periods of millions of years.
One can recognize in Figure 9 some common features between the evolution of the solstices and the trend of temperatures. We recalled in the introduction the basic equation from Milankovic's thesis (1920) that links time variations of heat W received at a given location on Earth to solar insolation, known functions of the location coordinates, solar declination and hour angle, with an inverse square dependence on the Sun-Earth distance. We translate the drift of solstices as a function of distance to the Sun into the geometrical insolation theory of Milankovic (1920). Both the apparent drift of solstices of Earth around the Sun and the global mean temperature exhibit a strong 60yr oscillation (Figure 8c and 10).
It may seem that we navigate between two pitfalls: remembering that correlation does not imply causality at the risk of discounting a potentially interesting relationship or accepting the causality when it does not exist, at the risk of pursuing a non existent theory. We do acknowledge that one should not jump too fast to conclusions, yet the probability of a chance coincidence in Figure 10 appears very low. Correlation certainly does not imply causality when there is no accompanying model. But in the case studied in this paper, equation (1) can be considered as a model that is widely accepted. Equation (1) links the time derivative of insolation with the inverse square of the Sun-Earth distance. On Figure 10 shifting the inverse square of the 60yr iSSA drift of solstices by 15 years with respect to the first derivative of the 60yr iSSA trend of temperature, that is exactly a quadrature in time, puts the two curves in quasi-exact agreement. This is a case of agreement between observations and a mathematical formulation. This new finding joins a host of recent results that argue in the same direction Bank et Scafetta (2022); the hypothesis proposes that a forcing by the giant Jovian planets is exerted on a vast number of solar and terrestrial phenomena, including as shown in this paper global surface temperature. This forcing induces a number of responses from the Earth's rotation axis, hence on climate at many time scales. This is a sort of extension of the Milankovic theory of climate to a range of periods that are much shorter than the ∼20 kyr minimum often associated with this theory.