A Secular Solar System Resonance that Disrupts the Dominant Cycle in Earth’s Orbital Eccentricity (g 2 − g 5): Implications for Astrochronology

The planets’ gravitational interaction causes rhythmic changes in Earth’s orbital parameters (also called Milanković cycles), which have powerful applications in geology and astrochronology. For instance, the primary astronomical eccentricity cycle due to the secular frequency term (g 2−g 5) (∼405 kyr in the recent past) utilized in deep-time analyses is dominated by the orbits of Venus and Jupiter, i.e., long eccentricity cycle. The widely accepted and long-held view is that (g 2−g 5) was practically stable in the past and may hence be used as a “metronome” to reconstruct accurate geologic ages and chronologies. However, using state-of-the-art integrations of the solar system, we show here that (g 2−g 5) can become unstable over long timescales, without major changes in, or destabilization of, planetary orbits. The (g 2−g 5) disruption is due to the secular resonance σ 12 = (g 1 − g 2) + (s 1 − s 2), a major contributor to solar system chaos. We demonstrate that entering/exiting the σ 12 resonance is a common phenomenon on long timescales, occurring in ∼40% of our solutions. During σ 12-resonance episodes, (g 2−g 5) is very weak or absent and Earth’s orbital eccentricity and climate-forcing spectrum are unrecognizable compared to the recent past. Our results have fundamental implications for geology and astrochronology, as well as climate forcing, because the paradigm that the long eccentricity cycle is stable, dominates Earth's orbital eccentricity spectrum, and has a period of ∼405 kyr requires revision.


INTRODUCTION
Laying the foundations of chaos theory, Henri Poincaré wrote: "It may happen that small differences in the initial conditions produce very great ones in the final phenomena.A small error in the former will produce an enormous error in the latter.Prediction becomes impossible . .." (Poincaré 1914).In reference to the solar system, the sensitivity to initial conditions is indeed a key feature of the large-scale dynamical chaos in the system, which has been confirmed numerically (Sussman & Wisdom 1988;Laskar 1989;Ito & Tanikawa 2002;Morbidelli 2002;Varadi et al. 2003;Batygin & zeebe@soest.hawaii.eduLaughlin 2008;Zeebe 2015;Brown & Rein 2020;Abbot et al. 2023).Dynamical chaos affects the secular frequencies g i and s i (see Appendix A), where the terms (g 4 −g 3 ) and (s 4 −s 3 ), for instance, show chaotic behavior already on a 50-Myr time scale.As a result, astronomical solutions diverge around t = ±50 Myr, which fundamentally prevents identifying a unique solution on time scales > ∼ 10 8 y (Laskar et al. 2004;Zeebe 2017;Zeebe & Lourens 2019).The chaos therefore not only severely limits our understanding and ability to reconstruct and predict the solar system's history and longterm future, it also imposes strict limits on geological and astrochronological applications such as developing a fully calibrated astronomical time scale beyond ∼50 Ma (for recent efforts, see Zeebe & Lourens 2019, 2022).In contrast to these limitations (largely due to unstable terms such as (g 4 −g 3 ) and (s 4 −s 3 )), another frequency term appears more promising as it shows more stable behavior.For example, it has hitherto been assumed that (g 2 −g 5 ) was practically stable in the past and has been suggested for use as a "metronome" in deep-time geological applications, i.e., far exceeding 50 Ma (Laskar et al. 2004;Kent et al. 2018;Spalding et al. 2018;Meyers & Malinverno 2018;Montenari 2018;Lantink et al. 2019;De Vleeschouwer et al. 2024).The (g 2 −g 5 ) cycle, which is the dominant term in Earth's orbital eccentricity in the recent past (∼405 kyr, see Fig. 1a) may thus have been regarded as an island of stability in a sea of chaos.However, we show in this contribution that also (g 2 −g 5 ) can become unstable over long time scales.
Solar system dynamics affect Earth's orbital evolution, as well as Earth's climate, which is paced by astronomical cycles on time scales > ∼ 10 kyr.The cycles include precession and obliquity cycles of Earth's spin axis (∼20 and 40 kyr in the recent past), and the short and long eccentricity cycles (∼100 and 405 kyr, see Figs. 1a and 6) (Milanković 1941;Montenari 2018).Orbital eccentricity is controlled by the solar system's orbital dynamics and is the focus of this study.The primary tuning target used in astrochronology and cyclostratigraphy for deep-time stratigraphic age models is the long eccentricity cycle (LEC) because it is widely assumed to be stable in the past (see above).Dynamically, Earth's orbital eccentricity and inclination cycles originate from combinations of the solar system's fundamental frequencies, called g-and s-frequencies (or modes), loosely related to the apsidal and nodal precession of the planetary orbits (Fig. A.1).The LEC is dominated by Venus' and Jupiter's orbits, viz., (g 2 − g 5 ), or g 25 for short, and represents the strongest cycle in Earth's eccentricity spectrum in the recent past (see Figs. 1 and 6).Assuming a stable LEC may appear plausible because Jupiter (dominating g 5 ) is the most massive planet and less susceptible to perturbations.Astronomical computations have confirmed g 5 's stability and hitherto did not indicate instabilities in g 25 (Berger 1984;Quinn et al. 1991;Varadi et al. 2003;Laskar et al. 2004;Zeebe 2017;Spalding et al. 2018;Zeebe & Lourens 2022).However, compared to g 5 , g 2 's long-term stability is less certain but has been overlooked so far.The long-term stability of g 25 and the LEC is critical for, e.g., climate forcing/insolation, constructing accurate age models and chronologies, expanding the evidence for the astronomical theory of climate into the more distant past, extending the astronomically calibrated geological time scale into deep time, and more (see discussion).Below we show that the LEC can become unstable over long-time scales owing to g 2 , without major changes in, or desta-bilization of, planetary orbits.Orbital destabilization is a known, separate dynamical phenomenon relevant to the future, see below.
For the present study, we performed state-of-the-art solar system integrations, including the eight planets and Pluto, a lunar contribution, general relativity, the solar quadrupole moment, and solar mass loss (see Section 2) (Zeebe 2017;Zeebe & Lourens 2019;Zeebe 2022Zeebe , 2023)).Initial conditions at time t 0 were taken from the latest JPL ephemeris DE441 (Park et al. 2021) and the equations of motion were numerically integrated to t = −3.5 Gyr (beyond −3.5 Gyr parameters such as the lunar distance have large uncertainties, see Section 2).Owing to solar system chaos, the solutions diverge around t = −50 Myr, which prevents identifying a unique solution on time scales > ∼ 10 8 y (see above).Hence we present results from long-term ensemble integrations to explore the possible solution/phase space of the system (see Section 2).Importantly, because of the chaos, each ∼10 8 y interval of the integrations represents a snapshot of the system's general/possible behavior that is largely independent of the actual numerical time of a particular solution (provided here that t < −τ 12 , where τ 12 is of order 10 8 to 10 9 y, see below).In other words, a numerical solution's behavior around, say, t = −1.5 Gyr may represent the actual solar system around t = −600 Myr and so on.

Solar System Integrations
Solar system integrations were performed following our earlier work (Zeebe 2017;Zeebe & Lourens 2019;Zeebe 2022) with our integrator package orbitN (v1.0) (Zeebe 2023), using a symplectic integrator and Jacobi coordinates (Wisdom & Holman 1991;Zeebe 2015a).The open source code is available at zenodo.org/records/8021040 and github.com/rezeebe/orbitN.The methods used here and our integrator package have been extensively tested and compared against other studies (Zeebe 2017;Zeebe & Lourens 2019;Zeebe 2022Zeebe , 2023)).For the present study, we also included simulations with an independent integrator package (HNBody) (Rauch & Hamilton 2002) and found the same dynamical behavior.All simulations include contributions from general relativity, available in orbitN as Post-Newtonian effects due to the dominant mass.The Earth-Moon system was modeled as a gravitational quadrupole (Quinn et al. 1991;Varadi et al. 2003;Zeebe 2017Zeebe , 2023)).Initial conditions for the positions and velocities of the planets and Pluto were generated from the latest JPL ephemeris DE441 (Park et al. 2021) using the SPICE toolkit for Matlab.Coordinates were obtained at JD2451545.0 in the ECLIPJ2000 and (c) Spectra of solutions R06 (Run 06) and R45 during σ12-resonance episodes (see Section 3) centered at t = −1180 Myr and t = −3260 Myr, respectively.Note the unrecognizable spectrum pattern in (c) compared to (a) and the reduced/absence of power around the LEC frequency of ∼3.2" y −1 (∼405 kyr) in (b) and (c).The peaks around 10 to 13" y −1 are due to the short eccentricity cycle.

Ensemble Integrations
We performed ensemble integrations of the solar system with a total of N = 64 members.Note that a larger N is unnecessary for the current problem, which samples a common phenomenon (∼40% of solutions), not a rare event, which requires large N (Abbot et al. 2023).Different solutions were obtained by offsetting Earth's initial position by a small distance (largest offset ∆x 0 ≃ 1 × 10 −12 au), which is within observational uncertainties (Zeebe 2015(Zeebe , 2017)).The different ∆x 0 lead to complete randomization of solutions on a time scale of ∼50 Myr due to solar system chaos.We also tested different histories of the Earth-Moon distance (R), which has little effect on our results (see Section 2.3).Because of the large uncertainties in R prior to ∼3 Ga, we restrict our integrations to t = −3.5 Gyr.Our solutions are available at www.ncdc.noaa.gov/paleo/study/xxxxxand www2.hawaii.edu/∼ zeebe/Astro.html.

Past Earth-Moon distance
Our integrations included a lunar contribution, i.e., a gravitational quadrupole model of the Earth-Moon system (Quinn et al. 1991;Varadi et al. 2003;Zeebe 2017Zeebe , 2023)).In the present context, the lunar contribution has a relatively small effect on the overall dynamics, yet the integration requires the Earth-Moon distance (R) as parameter at a given time in the past.We tested two approaches, both avoiding the known problem of unrealistically small R at −3.5 Gyr (see Fig. 3).(i) A linear extrapolation of R into the past starting with dR 0 /dt close to the present rate and (ii) a 3rd-order polynomial fit to observations.The two approaches made essentially no difference in our computations and both yielded solutions including σ 12 -resonance intervals at a similar frequency (see below).For the observational constraints on R, we selected robust data sets based on the reconstruction of Earth's axial precession frequency obtained by cyclostratigraphic studies (Meyers & Malinverno 2018;Lantink et al. 2022;Sørensen et al. 2020;De Vleeschouwer, D. et al 2023) (see Fig. 3).The classical integration of precession equations starting at the present rate dR 0 /dt (see Fig. 3, green dashed line) follows MacDonald (1964); Goldreich (1966); Touma & Wisdom (1994).

Time series analysis of astronomical solutions
The solar system's fundamental g-and s-frequencies were determined from the output of our numerical integrations using fast Fourier transform (FFT) over consecutive 20-Myr intervals.For the spectral analysis (see Figs. 4 and 5), we used Earth's orbital elements and the classic variables: where e, I, ϖ, and Ω are eccentricity, inclination, longitude of perihelion, and longitude of ascending node, respectively.The spectra for Earth's k and q, for example, show strong peaks at nearly all g-and s-frequencies, respectively (see Fig. 5).The g-and s-modes are loosely related to the apsidal and nodal precession of the planetary orbits (see Fig. A.1).However, there is generally no simple one-to-one relationship between a single mode and a single planet, particularly for the inner planets.The system's motion is a superposition of all modes, although for the outer planets, some modes are dominated by a single planet.

Resonant angle
The resonant angle θ 12 associated with the σ 12 resonance (see Eq. ( 6)) was determined following Lithwick & Wu (2011).The method is numerically efficient and easy to implement.Consider Eqs. ( 1) and ( 2), and use sin(I/2) ≃ I/2 (applicable to small I, as in our solutions).The variable pairs (h, k) and 2(p, q) can then be combined into two complex variables for each planet (î = √ −1): where index k = M, V, E, . . ., N refers to the planets.The (z k , ζ k ) for k = M, V , for instance, were determined from Mercury's and Venus' computed orbital elements.Next, we applied a simple bandpass filter (rectangular window) centered on the fundamental frequencies g i and s i of interest (index i).For example, for i = 1, 2, the passed frequency range was set to the present rate.Cyan: 3rd-order polynomial fit to observations.Using R/R E based on the blue and cyan lines made essentially no difference in our computations; both approaches yielded solutions including σ12-resonance intervals at a similar frequency (see Section 3). [ , where ℑ and ℜ denote the imaginary and real part of a complex number.Finally, the resonant angle θ 12 associated with σ 12 (see Eq. ( 6)) is calculated as: (5) 3. RESULTS

Secular frequencies and eccentricity
Contrary to expectations, we found in ∼40% of the solutions that g 25 was not stable at a period P 25 ≃ 405 kyr but shifted abruptly due to shifts in g 2 (Fig. 4).Importantly, in those cases the g 2 spectral peak usually split into two peaks at significantly reduced power (Fig. 5), resulting in a very weak or absent LEC (Fig. 6).Note that for, e.g., geological applications, the weak/absent LEC is crucial, not the actual value of the P 25 shift (Fig. B.1), which is immaterial because it would be unidentifiable in a stratigraphic record owing to the low g 25 power.Time series analysis of the solutions (see Section 2) revealed that the weak LEC intervals are associated with a secular resonance between the g-and s-modes dominated by Mercury and Venus (|g Several observations lend confidence to the robustness of our astronomical computations.(i) The methods used here and our integrator package have been extensively tested and compared against other studies (Zeebe 2017;Zeebe & Lourens 2019;Zeebe 2022Zeebe , 2023)).(ii) The σ 12 resonance was recognized previously, although to our knowledge only by two studies (Lithwick & Wu 2011;Mogavero & Laskar 2022) and not its effect on g 25 /LEC (see below).(iii) We tested an independent integrator package (HNBody Rauch & Hamilton (2002), see Section 2) and found the same dynamical behavior.(iv) Re-examination of previous 5-Gyr future integrations from this group (Zeebe 2015) also revealed various solutions with σ 12 -resonance intervals.(v) Total energy and angular momentum errors were small throughout our present 3.5-Gyr integrations (relative errors in test runs: < ∼ 6 × 10 −10 and < ∼ 7 × 10 −12 , Fig. 2) and our numerical timestep sufficiently resolves Mercury's perihe-Figure 4. Evolution of fundamental solar system frequencies.The g-and s-frequencies (in arcsec y −1 = " y −1 ) were determined from our solar system integrations using fast Fourier transform (FFT) over consecutive 20-Myr intervals and Earth's k and q variables (see Section 2).The g-and s-modes are loosely related to the apsidal and nodal precession of the planetary orbits (see Fig. A.1). Solutions including σ12-resonance intervals (∼40%) are highlighted in color, the remaining solutions are displayed in gray.The frequencies g1, s1, and s2 drift most strongly over time owing to chaotic diffusion.In addition, g2 shows large and rapid shifts (spikes) at specific times when the spectral g2 peak splits into two peaks at significantly reduced power during σ12-resonance episodes (see Fig. 5).Alternating maximum power between the two peaks then causes the spikes in g2.As a result, g25 = (g2 − g5) is unstable and weak/absent during σ12-resonance intervals (see Fig. B.1). g5, g6, and s6 (dominated by Jupiter and Saturn) are practically stable over 3.5 Gyr (s5 is zero due to conservation of total angular momentum/existence of an invariable plane, see  lion (Wisdom 2015;Hernandez et al. 2022;Abbot et al. 2023) (see Section 2).
During σ 12 -resonance episodes, Earth's orbital eccentricity pattern and hence Earth's climate forcing spectrum due to eccentricity becomes unrecognizable compared to the recent past (Figs. 1 and 6).For example, a geologist examining a climate record exhibiting Milanković cycles within the σ 12 resonance (e.g., Fig. 6bd) would fail to identify the rhythm as eccentricity cycles, given the currently known pattern (Fig. 6a).The resonance motifs are so different that (coincidentally) some frequency and amplitude modulation (AM) features (Fig. 6c) show more similarities with Mars' orbital inclination in the recent past (Fig. C.1) (Zeebe 2022) than Earth's eccentricity (Fig. 6a).The estimated time scale τ 12 for a possible σ 12 -resonance occurrence (τ 12 = temporal distance to the present) is of order 10 8 to 10 9 y.In several solutions, we found reduced g 2 and g 25 power (lower than short eccentricity power), as well as unusual eccentricity patterns, at t < ∼ −500 Myr (see e.g., Fig. 6d) and in one solution at t ≃ −420 Myr.However, we have so far tested only 64 solutions (see Section 2), hence the youngest possible age of a σ 12 -resonance interval that could be detectable in the geologic record is yet unknown.The duration of a σ 12 -resonance episode may range from a few Myr to tens of millions of years (multiple entries/exits often occurring over several 100 Myr, see Figs. 6 and 9).

Insolation and climatic precession
The total mean annual insolation (or energy W ) Earth's receives is a function of its orbital eccentricity (e ♁ ): In the recent past, 0 < ∼ e ♁ < ∼ 0.06, whereas during σ 12resonance episodes, max{e ♁ } may be as low as ∼0.04 (Fig. 6c).Thus, the relative variation/difference in W between eccentricity maxima and minima ((1 − 0 2 ) − 1 2 = 1) is substantially reduced by the factor: Thus, in addition to a weak/absent LEC, both the total variation in eccentricity climate forcing and the extreme values are diminished on a 10 6 -year time scale during σ 12 intervals.Moreover, the σ 12 resonance causes major changes in climatic precession (p), the primary climate driver on the shortest Milanković time scale (∼20 kyr in the recent past).The main p frequencies are given by Ψ + g i , where Ψ is the lunisolar precession frequency.
The disruption of g 2 (and hence Ψ + g 2 ) causes major changes in p's total amplitude and AM (see Figs. 7 and 8).For example, during σ 12 -resonance intervals, the forcing power at the Ψ + g 2 precession frequency may drop by orders of magnitudes compared to the recent past (Fig. 8).The altered forcing in both, eccentricity and climatic precession, would scale down the climate response to orbital forcing, and hence its expression in geological sequences, as well as affect threshold behavior for triggering orbitally forced climate events (for recent examples such as the Paleocene-Eocene Thermal Maximum and the Eocene hyperthermals, see Zeebe & Lourens ( 2019)).

σ 12 resonance
In a quasi-periodic (non-chaotic) system, the fundamental frequencies are constant over time.In contrast, the solar system's frequencies change over time owing to dynamical chaos (Fig. 4), priming the system to enter/exit the σ 12 resonance over long time scales.Several solar system resonances have been studied previously, including (g 1 − g 5 ) − (s 1 − s 2 ) and 2(g 4 − g 3 ) − (s 4 − s 3 ) (Laskar 1990;Sussman & Wisdom 1992;Batygin et al. 2015;Ma et al. 2017;Zeebe 2022;Abbot et al. 2023).However, to our knowledge only two studies recognized σ 12 , yet did not investigate its consequences for g 25 and Earth's orbital eccentricity (Lithwick & Wu 2011;Mogavero & Laskar 2022).σ 12 may be characterized by the resonant angle θ 12 (see Eq. ( 5)), where ϖ * and Ω * are associated with the g-and s-modes (analog to longitude of perihelion and ascending node, but not of individual planets, see Appendix A).Chaos is often associated with resonant angles that alternate between circulation and libration (Fig. 9).Generally, θ 12 circulated in our solutions without σ 12 -resonance episodes, whereas θ 12 -circulation and libration occurred in solutions that showed σ 12 -resonance episodes and a weak/absent LEC.The latter case was usually associated with intervals of slightly elevated eccentricity in Mercury's orbit (0.25 < ∼ e < ∼ 0.35, see Fig. 9).Importantly, none of our solutions showed high eccentricities (e > ∼ 0.4), which could indicate progressing chaotic behavior or a potential destabilization of the inner solar system -known, separate dynamical phenomena, most relevant to future chaos, studied previously (Laskar 1990;Sussman & Wisdom 1992;Lithwick & Wu 2011;Batygin et al. 2015;Zeebe 2015;Brown & Rein 2020;Abbot et al. 2023).Solutions displaying any significant destabilization in the past can of course be excluded (incompatible with the solar system's known history).Furthermore, given that all our solutions showed at most slightly elevated e demonstrates that the sys-  tem can enter/exit the σ 12 resonance without major changes in, or destabilization of, planetary orbits.Thus, past σ 12 -resonance episodes and a weak LEC are a possible and likely dynamical phenomenon, present in ∼40% of our solutions.

IMPLICATIONS
We anticipate far-reaching consequences of our findings for (i) exploring the effects of secular resonances (particularly σ 12 ) on the long-term dynamical evolution, chaos, and planetary climates in the solar system, (ii) understanding and unraveling Earth's past climate forcing and climate change via parameters including eccentricity (total insolation) and climatic precession (see Eq. ( 7) and Figs. 1,6,7,8), (iii) reconstructing the solar system's chaotic dynamics constrained by geologic evidence (Ma et al. 2017;Meyers & Malinverno 2018;Zeebe & Lourens 2019;Olsen et al. 2019), (iv) expanding the evidence for the astronomical theory of climate in yet understudied parts of Earth's history (e.g., the Precambrian), (v) studying effects of deep-time Milanković forcing on Earth's environmental and climatic long-term evolution, and (vi) extending the astronomically calibrated geological time scale (Montenari 2018; Zeebe & Lourens 2022) into deep time.
It appears that the σ 12 secular resonance and its effect on solar system dynamics and planetary climates has been understudied thus far.To our knowledge, only two studies recognized σ 12 , yet did not investigate its consequences on, for instance, (g 2 −g 5 ) and Earth's or- Note different lunisolar precession frequency Ψ in the past owing to the Earth-Moon system's evolution (for details and code, see (Zeebe 2022), github.com/rezeebe/snvec,www2.hawaii.edu/∼ zeebe/Astro.html,and Section 2.3).Top: Standard spectrum in the recent past centered at t = −3 Myr (nearly identical in all solutions).Bottom: Spectrum based on solution R06 centered at t = −1180 Myr.Note the split and reduced power of Ψ + g2 in bottom vs. top panel (red circles).
bital eccentricity (Lithwick & Wu 2011;Mogavero & Laskar 2022).Several of the secular modes (or terms) related to the g i and s i (or differences between pairs) show multiple, strong interactions for i = 1, . . . 4. In other words, secular resonances usually affect multiple frequency pairs.For example, as shown here, the σ 12 resonance has a major impact on (g 1 −g 2 ) and (s 1 −s 2 ), but also on (g 2 −g 5 ).Moreover, because there is no simple one-to-one relationship between a single mode and a single inner planet (the motion is a superposition of all modes), resonances (say σ ij with i, j = 1, . . .4, i ̸ = j) affect the entire inner solar system.Here, we have only investigated σ 12 's effect on (g 2 −g 5 ) and Earth's orbital eccentricity.Future work should explore whether there are other important, yet unknown, effects of σ 12 (and other resonances) on the dynamics and planetary climates in the inner solar system.
Our results have fundamental implications for, e.g., current astrochronologic and cyclostratigraphic prac-tices, which are based on the paradigm that the LEC is stable, dominates the eccentricity spectrum, and has a period of ∼405 kyr.Given our findings that the σ 12 resonance is a common phenomenon (occurring in ∼40% of our solutions), the assumption of a stable 405 kyr cycle in deep time can no longer be made.Specifically, the possibility of an unstable period and weakened LEC amplitude requires rethinking of currently employed strategies for building accurate and high-resolution ('floating' or radio-isotopically anchored) astrochronologic age models.The presumed 405-kyr "metronome" was particularly important for constructing pre-Cenozoic age models, where reliable astronomical solutions are absent owing to solar system chaos.Deep-time astrochronologies have thus far critically relied on identification of the LEC because it is the only Milanković cycle whose period has been widely regarded as stable (Laskar et al. 2004;Kent et al. 2018;Spalding et al. 2018) and of sufficiently large amplitude to be typically expressed in sed- imentary sequences.Our results indicate that, prior to several hundred Myrs in the past, the LEC may have become unstable over multi-million year intervals.Notably, on these time scales the periods of other critical Milanković parameters (climatic precession and obliquity) are also more uncertain due to changes in the Earth-Moon system's tidal evolution.
Does the presently explored geologic record provide examples consistent with our astronomical calculations?We note here a recently discovered section in the ∼2.46 Ga Joffre Member of the Brockman Iron Formation (Joffre Falls, Western Australia), in which a dominant short eccentricity cycle was reported (∼100 kyr), compared to a relatively weak ex-pression at the scale of the interpreted LEC (Lantink et al. 2022) (see Fig. D.1).The interpreted eccentricity modulation pattern in the Joffre Falls section differs from typical Cenozoic precession-eccentricity dominated records, which often display a strong 1:4 hierarchy of long vs. short eccentricity cycles.One first-order interpretation for the unusual bundling pattern is a complex nonlinear response of the paleoclimate and/or sedimentary system to orbital forcing, which is still poorly understood for the ancient deposits.Yet, given our findings, a fundamentally different pattern of Earth's orbital eccentricity variations (i.e., a weakened LEC) at the time of deposition provides an alternative explanation (see Fig. D.1).Further investigations are required to confirm past σ 12 -resonance episodes in geologic sequences.We propose that future exploration of highquality and rhythmic sediment successions (especially of Precambrian age) will be critical in constraining the LEC's past stability and hence the history of the solar system's chaotic evolution.
Generally, the possibility of an unstable/weak LEC argues strongly for internal consistency checks and tests of eccentricity-related cycles interpreted in stratigraphic sequences at multiple levels.For example, future studies need to include consistency checks between the period of short eccentricity and associated g-frequencies and the period of the interpreted (g 4 − g 3 ) cycle and/or other very long period eccentricity modulations.At a more advanced level, the internal consistency of all g-and s-frequencies that can be extracted from the sequence need to be examined, for which algorithms are already available (Meyers & Malinverno 2018;Olsen et al. 2019).Furthermore, the uncertainty in LEC stability substantially increases the ambiguity in interpreting cycle ratios of the eccentricity-precession forcing.Hence, independent sedimentation rate checks (preferably from accurate radiometric ages) will become inevitable to verify deep-time Milanković interpretations based on observed stratigraphic cycle hierarchy.Moreover, when significant obliquity signals are present, eccentricity-related cycle pattern will be more difficult to distinguish from those expected for obliquity.
hence for circulating Ω, the time-averaged nodal precession is retrograde (i.e., in the opposite direction as the orbital motion, see large colored arrows).Given conservation of total angular momentum (L), there exists an invariable plane perpendicular to L, which is fixed in space.It follows that one of the s frequencies is zero (s5, see main text).
B. PERIOD OF (g 2 −g 5 ) The secular frequency g2 shows large and rapid shifts (spikes, see Fig. 4) at specific times when the spectral g2 peak splits into two peaks at significantly reduced power during σ12-resonance episodes (see Fig. 5).Alternating maximum power between the two peaks then causes the spikes in g2 and hence in g25 (Fig. B.1).As a result, g25 is unstable and weak/absent during σ12-resonance intervals.Note that for, e.g., geological applications, the weak/absent LEC is crucial, not the actual value of the P25 shift, which is immaterial because it would be unidentifiable in a stratigraphic record owing to the low g25 power.

C. MARS' INCLINATION
The σ12-resonance motifs are fundamentally different from the recent past (Fig. 6), some of which (coincidentally) show more similarities with Mars' orbital inclination in the recent past (  (Lantink et al. 2022) shown on the right.Note the larger-scale modulations in the thickness and relief of the shaley beds (degree of weathering) forming two distinctive darker 'bundles' defined by cycles 17-19 and 23-28.This pattern deviates from an expected strong ∼1:4 bundling pattern or ratio between the medium-and large-scale cycles in case of a strong and stable LEC.(However, note the ∼1:4 ratio visible in the filter amplitude of the bandpass filtered cyclicity on the right.)Higher up in the stratigraphy, the shaley layers become weaker overall and thus any larger-scale modulations are more difficult to recognize.Nevertheless, we count at least 6 medium-scale cycles until the next more distinctive shaley interval, i.e., again no clear 1:4 ratio as would be expected in case of a strong and stable LEC.

Figure 3 .
Figure 3. Past Earth-Moon distance (R) in units of Earth radii (R E ). Green dashed line: Integration of precession equations starting at present rate dR0/dt (see text), yielding (well-known) unrealistic past R/R E .Red diamonds: Observational estimates based on robust data sets from cyclostratigraphic studies.Blue and cyan lines: Used in the present study.Blue: linear extrapolation into the past starting with dR0/dt close to

Figure 6 .
Figure 6.Earth's orbital eccentricity from ensemble integrations.(a) Eccentricity pattern in the recent past (last 20 Myr) with short and long (∼100 and ∼405 kyr) cycles (nearly identical in all solutions).Note the strong bundling of four short cycles into one long cycle (highlighted by 405 kyr filter).Panels (b, c, d) display examples of solutions during σ12-resonance intervals.(b) Solution R28 (Run 28) over a 20-Myr interval centered at t = −1450 Myr.The LEC is virtually absent and the eccentricity pattern is unrecognizable compared to (a).(c) Solution R06 over a 20-Myr interval centered at t = −1180 Myr.In addition to a weak LEC, the maximum eccentricity is reduced to ∼0.04, which affects the total insolation Earth receives over one year (see text).(d) Solution R44 over a 20-Myr interval centered at t = −528 Myr; the LEC is present but weak relative to the short eccentricity cycle.

Figure 7 .
Figure 7. Climatic precession p = e sin ω, where ω is the longitude of perihelion measured from the moving equinox (for details and code, see (Zeebe 2022), github.com/rezeebe/snvec,www2.hawaii.edu/∼ zeebe/Astro.html,and Section 2.3.(a) p in the recent past (last 1 Myr), nearly identical in all solutions.Panels (b, c) display examples of p based on orbital solutions during σ12-resonance intervals.Note the different lunisolar precession frequencies (carrier frequency) in the past owing to the Earth-Moon system's evolution.(b) p based on solution R28 over a 1-Myr interval centered at t = −1450 Myr.The amplitude modulation (AM) pattern fundementally differs from (a).(c) p based on R06 over a 1-Myr interval centered at t = −1180 Myr.In addition to an altered AM pattern, the total amplitude is reduced compared to (a).
Fig. C.1) than Earth's eccentricity.D. SECTION AT JOFFRE FALLS, WESTERN AUSTRALIA In the recently discovered Joffre Falls section (∼2.46 Ga, Western Australia), a dominant short eccentricity cycle was reported (∼100 kyr) and a relatively weak expression of the LEC (Lantink et al. 2022) (see Fig. D.1).The regular medium-scale (∼85 cm) alternations of thicker units of banded iron formation and thinner, softer interval of a more shaley lithology (Fig. D.1, left) have been interpreted as the expression of short eccentricity (Lantink et al. 2022).Horizons highlighted in blue (Fig. D.1, left) correspond to cycle numbers in the original log Figure C.1.Mars' inclination over the past 20 Myr (nearly identical in all solutions, reference frame: ECLIPJ2000).