North Atlantic Drift Sediments Constrain Eocene Tidal Dissipation and the Evolution of the Earth‐Moon System

Cyclostratigraphy and astrochronology are now at the forefront of geologic timekeeping. While this technique heavily relies on the accuracy of astronomical calculations, solar system chaos limits how far back astronomical calculations can be performed with confidence. High‐resolution paleoclimate records with Milankovitch imprints now allow reversing the traditional cyclostratigraphic approach: Middle Eocene drift sediments from Newfoundland Ridge are well‐suited for this purpose, due to high sedimentation rates and distinct lithological cycles. Per contra, the stratigraphies of Integrated Ocean Drilling Program Sites U1408–U1410 are highly complex with several hiatuses. Here, we built a two‐site composite and constructed a conservative age‐depth model to provide a reliable chronology for this rhythmic, highly resolved (<1 kyr) sedimentary archive. Astronomical components (g‐terms and precession constant) are extracted from proxy time‐series using two different techniques, producing consistent results. We find astronomical frequencies up to 4% lower than reported in astronomical solution La04. This solution, however, was smoothed over 20‐Myr intervals, and our results therefore provide constraints on g‐term variability on shorter, million‐year timescales. We also report first evidence that the g4–g3 “grand eccentricity cycle” may have had a 1.2‐Myr period around 41 Ma, contrary to its 2.4‐Myr periodicity today. Our median precession constant estimate (51.28 ± 0.56″/year) confirms earlier indicators of a relatively low rate of tidal dissipation in the Paleogene. Newfoundland Ridge drift sediments thus enable a reliable reconstruction of astronomical components at the limit of validity of current astronomical calculations, extracted from geologic data, providing a new target for the next generation of astronomical calculations.


Introduction
Accurate and highly resolved age models are necessary for unraveling causal relationships in paleoclimatology. For this reason, cyclostratigraphy and astrochronology became indispensable tools over the last few decades. However, astrochronologies crucially depend on reliable orbital calculations of Solar System dynamics far back in time (Sinnesael et al., 2019). Variations in Earth's astronomical parameters of eccentricity, obliquity and precession significantly alter the distribution and amount of incoming solar radiation in space and time (Laskar et al., 2004;Milanković, 1941). Following Quinn et al. (1991), recent orbital models of the Solar System are constructed by direct numerical integration (Laskar et al., 2004, Laskar et al., 2011aLaskar, Gastineau, et al., 2011;Varadi et al., 2003;Zeebe, 2017;Zeebe & Lourens, 2019). These numerical integrations describe the orientation and shape of the orbital bodies considered, and can be characterized by a set of fundamental frequencies that describe the motion of orbits within, and perpendicular to their orbital planes (the so-called g i terms or eccentricity modulation terms, and s i terms or orbital inclination terms, respectively) (Pälike, 2005). Each of those frequencies is associated with the particular orbital movement of one of the nine (dwarf) planets in the Solar system (the "i" index), whereby i = 1 refers to Mercury and i = 9 to Pluto. These fundamental frequencies change slowly on million-year time scales (Laskar, 1990) and are therefore denoted "secular terms". To date, our knowledge on the variability and trends of these terms comes almost exclusively from astronomical models, and first attempts to reverse the cyclostratigraphic approach indicate that the past variability of these terms could have been greater than previously thought (Meyers & Malinverno, 2018). Moreover, the orbital calculations also exhibit chaotic behavior (Laskar, 1990) and a fuller characterization of the orbital system must involve resonant and non-linear terms. More recent orbital calculations that also include major asteroids like Ceres and Vesta (Laskar, Gastineau, et al., 2011c;Zeebe, 2017) have shown that the chaotic evolution of the solar system limits the time span for which valid calculations can be made of Earth's eccentricity and orbital inclination to around 60 Ma in the most optimistic case. However, different solutions diverge significantly beyond ∼48 Ma. Several other uncertainties, like, for example, the J2 quadrupole moment of the Sun (Laskar, 1999;Laskar et al., 2004;Zeebe & Lourens, 2019), further indicate that orbital solutions must be constrained by geological observations. Certain orbital amplitude modulation terms have been shown to be stable over long periods of time. These so-called astronomical metronomes can be used as a framework to construct detailed cyclostratigraphies in deep time. The 405-kyr long eccentricity term (g 2 -g 5 ), for example, has been shown to be stable beyond 60 Ma (Kent et al., 2018;Olsen et al., 2019). Likewise, the s 3 -s 6 orbital inclination term (∼173-kyr obliquity amplitude modulator) is still phase coherent between different orbital integrations back to 48 Ma (Boulila et al., 2018). While the above discussion relates to the orbital components, Earth's climatic precession and obliquity evolution also involves the precession frequency of the Earth p, which is influenced by the distribution of angular momentum within the Earth-Moon system. Tidal friction steadily slows the Earth's rotational speed, among other effects. To maintain angular momentum within the Earth-Moon system, the Earth-Moon distance increases as an immediate consequence. This effect results in a decrease of the precession frequency p over time and, therefore, in a progressive decline in the main frequency components of Earth's obliquity (p + s 1-4 and p + s 6 terms) and precession (p + g 1-5 ) periodicities. The current rate at which rotational energy dissipates is, however, a poor guide for the past. Present-day tidal dissipation would imply zero Earth-Moon distance at ∼1.5 Ga, whereas the Moon is ∼4.5 Gyr old (Green et al., 2017;Hansen, 1982;Waltham, 2015). While it is generally accepted that tidal dissipation must have been lower throughout much of Earth history, uncertainties on its temporal evolution remain (Daher et al., 2021). This ambiguity also holds for geologic Epochs as recent as the Eocene: Color reflectance (a*b*) data from Walvis Ridge ODP Site 1262 has been used both by Meyers and Malinverno (2018) and Zeebe and Lourens (2022) to reconstruct tidal dissipation around 55 Ma. Their approaches are very different, but both studies adopt iterative data-model fitting approaches that simultaneously derive floating age-depth models and the precession frequency p. Both teams likewise firmly conclude that early Eocene tidal dissipation must have been lower than the present-day value, but caution that more precise tidal dissipation reconstructions require independent age-depth model constraints. Thus, the lingering uncertainty on early Cenozoic tidal dissipation limits accurate insolation reconstructions to the Neogene (Lourens et al., 2001;Pälike & Shackleton, 2000; Zeeden motion and on the Earth-Moon interactions. These astronomical reconstructions based on geological data can now be used by astronomers to describe the evolution of the solar system further back in time than was previously possible. et al., 2014). To move forward, orbital terms extracted from accurately dated geological archives are now required to benchmark astronomical solutions regarding the evolution of the solar system.
The middle Eocene is a time of gradual global cooling, interrupted by a remarkable warm event called the Middle Eocene Climatic Optimum (MECO) (Bohaty & Zachos, 2003;Bohaty et al., 2009). The MECO is marked by large temperature, sea level and salinity variations in the North Atlantic, the Arctic Ocean and the North Sea (e.g., Cramwinckel et al., 2020;Marchegiano & John, 2022;van der Ploeg et al., 2023). The contourite drift deposit sequences of International Ocean Discovery Program (IODP) Site U1408 and Site U1410 (Newfoundland Ridge, North Atlantic, IODP Expedition 342, Figure 1) provide a high-resolution paleoceanographic archive throughout this time interval and have the potential to supply a test for Earth-Moon dynamics and the evolution of the solar system. The onset of contourite deposition occurred at both sites around 47 Ma, when an order-of-magnitude increase in terrigenous mass accumulation rate resulted in high overall sedimentation rates (>2 cm/kyr) (Boyle et al., 2017;Cappelli et al., 2019) and a sedimentary system sensitive to astronomical insolation forcing . Post-expedition, the science team collaborated to obtain cm-resolution benthic isotope data, as well as X-Ray Fluorescence (XRF) core scans for this unique Eocene sedimentary archive. However, the utility of this multi-proxy data set has been limited by problems related to obtaining a reliable composite depth-and time-scales, which in turn are due to complicated drift sediment lithologies at IODP Sites U1408 and U1410 (Boulila et al., 2018;Boulila & Hinnov, 2022;Cappelli et al., 2019;Vahlenkamp et al., 2018;Zeebe & Lourens, 2022). Hence, the first task for this study is to resolve the U1408-U1410 middle Eocene time-scale controversies. This is achieved by integrating high-resolution benthic foraminiferal stable isotope records (N = 3,424) and X-ray Fluorescence (XRF) derived elemental ratios (N = 9,662) into a new two-site composite section and a robust stratigraphic framework that acknowledges the fragmentary nature of the individual IODP Sites. With a revised and improved chronostratigraphic framework, these geochemical records provide the means to attain the second objective of extracting information on the long-term chaotic evolution of the solar system.

X-Ray Fluoresence and Benthic Stable Isotope Proxies
XRF-derived Ca/Fe ratios were measured at 2-cm resolution at MARUM (University of Bremen) and the Scripps Institution of Oceanography (University of California San Diego) on Avaatech core scanners. The studied interval consists of 9,662 Ca/Fe data points in the final two-site composite. Stable carbon and oxygen isotopes of benthic foraminifera Nuttallides truempyi (2-6 shells per analysis) were measured at 3-6 cm resolution by a consortium consisting of MARUM (Bremen University), National Oceanography Center (Southampton), University of California (San Diego), the University of California (Santa Cruz), and Yale University (New Haven). In all labs, 20cc sediment samples were dried and then washed through a 63-μm sieve. Two to six individuals of Nuttalides truempyi were picked and run for stable carbon and oxygen isotopes using a Kiel IV carbonate preparation device coupled to a ThermoFisher MAT 253 Plus isotope ratio mass spectrometer, using established dual-inlet techniques. The studied interval consists of 3,424 data points in the final two-site composite. Benthic Nuttallides truempyi δ 13 C and δ 18 O values are converted to Cibicidoides for comparison with the CENOGRID reference curve (Westerhold et al., 2020), using the conversion factors of Katz et al. (2003):

Numerical Analysis of Fundamental Frequencies
Astronomical frequencies were extracted from the geological data using numerical analysis of fundamental frequencies (NAFF). NAFF performs high-resolution frequency analysis of time series as designed by Laskar (1990), and previously used by Olsen et al. (2019) to extract astronomical frequencies from a Triassic-Jurassic lake depth rank series. Here, we apply the NAFF technique to the log(Ca/Fe) time-series (39.47-42.81 Ma), as well as to the amplitude envelope of the precession cycles within the log(Ca/Fe) data , and to the benthic isotope time-series (39.47-40.82 Ma). The NAFF method has been implemented in a macOS command line tool by Pälike (2021b). The astrochron R package was used to calculate multi-taper method spectra, Taner bandpass filters and Hilbert transforms for amplitude envelope extraction (Meyers, 2014).

TimeOptMCMC
To verify the NAFF result, and to quantify the uncertainties in the extracted astronomical frequencies, we additionally applied the TimeOptMCMC method (Meyers & Malinverno, 2018) to the log(Ca/Fe) time-series (39.47-42.81 Ma). The TimeOptMCMC function is implemented in the astrochron R package (Meyers, 2014). Through Bayesian Markov Chain Monte Carlo simulations, different plausible combinations of the solar system secular frequencies g 1 to g 5 and the precession constant p are compared to the stratigraphic data, and their match is reformulated in terms of a likelihood function. For this study, we apply TimeOptMCMC to the log(Ca/Fe) time-series rather than the depth-series, yet we allowed for minor perturbations of our time-axis by ±10%. Time-OptMCMC was run twice: First, TimeOptMCMC is run with the same prior beliefs on g-term and precession constant averages as in the Eocene analysis in Meyers and Malinverno (2018). Second, TimeOptMCMC is run with updated prior beliefs, based on the NAFF results. Both analyses each consist of 106 Markov Chain Monte Carlo runs, each of which consists of 200,000 samples, sufficient to explore the entire parameter space. Every chain is characterized by a burn-in phase, during which the combinations of astronomical frequencies converge toward a high-probability region of the posterior probability density function. We determined the burn-in phase in the same way as Meyers and Malinverno (2018), first computing the median likelihood value in the second half of the chain and then defining the burn-in end as the first sample in the chain that reaches a likelihood value greater than this median value. Finally, we combine the results of the 106 independent chains to display posterior distributions.

Dynamical Ellipticity, Tidal Dissipation and the Precession Constant
The periods for precession and obliquity are known to have been shorter in the geologic past, compared to the present-day (Berger et al., 1992). This is the direct result of changes in the precession constant p, which appears in all precession (p + g i ) and obliquity (p + s i ) arguments. The precession constant p arises from the solution of the Poisson equation describing the Earth-Moon system, which is Equation 7 in Berger et al. (1992) and Equation 8 in Laskar et al. (2004). This equation entails both the rotational angular velocity of the Earth and the dynamical ellipticity of the Earth. Both parameters change slowly through geologic time. The rotational angular velocity of the Earth decreases due to tidal friction and other energy dissipative effects (e.g., core-mantle friction). The dynamical ellipticity of the Earth changes due to plate tectonics or the build-up of continental ice sheets. Several cyclostratigraphers have reconstructed tidal dissipation and dynamical ellipticity from Neogene geologic data without translating them to their respective precession constant p (Lourens et al., 2001;Pälike & Shackleton, 2000;Zeeden et al., 2014). We developed a Graphical User Interface (De Vleeschouwer, 2023) that approximates the effect of dynamical ellipticity and tidal dissipation on the precession constant p, using the formulation in Laskar et al. (1993) and assuming S 0 and ε 0 constant. The interface was used to translate the results of Pälike and Shackleton (2000)

A Two-Site U1408-U1410 Composite
The middle Eocene drift deposit sequences of IODP Site U1408 and Site U1410 (Newfoundland Ridge, North Atlantic) consist of nannofossil clay marked by rhythmical alternations between more clay-rich and more nannofossil-rich endmembers. The clay-content variability occurs on the meter-scale (0.6-1.2 m) and was initially interpreted as obliquity-paced cyclicity (Boulila et al., 2018;Boulila & Hinnov, 2022;Vahlenkamp et al., 2018). However, a site-to-site correlation by Cappelli et al. (2019) revealed numerous hiatuses within the individual records, as well as splicing mistakes that doubled strata. The recovered middle Eocene sequences therefore represent less time than originally thought, and doubts emerged on the original astronomical interpretation. The suspicion arose that the lithological rhythm is reflecting precession rather than obliquity. The  (Figure 4), maximizing stratigraphic completeness. We then convert the composite from depth-to-time and assess the hypothesis that the distinct lithological cycle in clay-content represents precession rather than obliquity.
The construction of the U1408-U1410 composite (called CCSF-X) is underpinned by an integrated stratigraphic approach, combining all available magneto-, bio-and chemostratigraphic information that constrains the two-site correlation in the depth-domain. While we were able to make a site-to-site correlation throughout the studied interval, we had to adopt a different composite-building strategy in the upper and lower part of the composite. In the upper part of the composite (Bohaty composite in Figure 4), the site-to-site correlation reveals up to 10-m-thick intervals that are present at one site but not at the other ( Figure 3). Hence, the upper part of the composite is constructed by filling in gaps at one site with sections from the other. This approach results in a composite section that is significantly more expanded compared to the individual sites, and that is therefore closer to a continuous sedimentary representation of geologic time. In the lower part of the composite (Cappelli composite in Figure 4), the Cappelli et al. (2019) site-to-site correlation revealed that Site U1408 is the more complete record, with gaps and condensed intervals occurring more frequently at Site U1410. For that reason, Site U1410 was mapped onto the Site U1408 composite depth scale (Cappelli et al., 2019). In contrast to the upper part of the composite, Site U1410 does not contain any sections that can fill in for stratigraphic gaps at Site U1408. Obviously, this does not imply that Site U1408 is continuous. Instead, we suspect that there are missing cycles at Site U1408, which are also missing from Site U1410. We therefore infer that the lower composite is likely not time-continuous, in contrast to the upper part of the composite. Nevertheless, the site-to-site correlation in the lower composite provides long stretches of high-resolution stable isotope data (∼1 kyr resolution) which may come in useful for a wide range of paleoceanographic investigations (Figure 2).

A Precessional Rhythm for the Middle Eocene Newfoundland Ridge Lithological Cycles
The conversion from depth to time consists of twenty tie-points ( Figure 5). Linear interpolation is adopted in-between tie-points. The age-depth model was constructed by aligning the composite's isotope records with the CENOGRID stratigraphic backbone (Westerhold et al., 2020), while simultaneously aligning pronounced lithological cycles (i.e., high-amplitude Ca/Fe cycles) with eccentricity maxima in the La11 eccentricity solution (Laskar et al., 2011b), and by considering the available magneto-and biostratigraphic datums (Table 1). The age-depth model came about through a trial-and-error approach, until a satisfactory fit was achieved with as few age-depth tie-points as possible ( Figure 6). This low number of tie-points is essential to avoid the importation of astronomical terms from CENOGRID or the La11 solution into the U1408/U1410 composite, and thus to avoid circular reasoning during the extraction of astronomical components from that composite. Tie-points are typically spaced several hundred thousand years apart ( Figure 5), meaning there could only be an influence on the 405-kyr g 2 -g 5 eccentricity term. This potential influence is considered acceptable since the g 2 -g 5 term is extremely stable throughout geological history ("astronomical metronome," Hinnov, 2013), with both g 2 and g 5 exhibiting at least one order of magnitude less variability compared to other g-terms (Laskar et al., 2004).
The upper part of the composite (50-160 m CCSF-X) is essentially time-continuous on astronomical time-scales. We make this judgment based on the good fit in terms of bio-, magneto-, and chemostratigraphy ( Figures 5  and 6). In the lower part of the composite (160-250 m CCSF-X), magneto-and biostratigraphic datums indicate that more geological time is represented by less stratigraphy compared to the upper part. This could indicate a lower sedimentation rate, but we consider this option unlikely because the lithologic cycles in the lower part are slightly thicker, certainly not thinner, and of the same shape (i.e., the waveform of individual lithological cycles) as in the upper section ( Figure 7). Instead, we infer sedimentation rates to be slightly higher in the lower part of the composite, with hiatuses causing more time to be represented by less stratigraphy. We propose four hiatuses throughout the lower part of the composite: at 160. 61, 190.21, 214.55, and 221.12 m CCSF-X ( Table 2). The upper two hiatuses were placed at abrupt shifts in δ 13 C (>0.25‰ over a few centimeters) that coincided with an disturbed amplitude modulation pattern in Ca/Fe. The lower two hiatuses were placed to improve the bioand chemostratigraphic match between 190.21-214.55 m and 221.12-250 m CCSF-X. For the small section in-between, only limited confidence can be put in the proposed age-depth model, as it is solely based on the chemostratigraphic correlation between the U1408-U1410 composite and the CENOGRID reference curve (i.e., ODP Site 1,263 in this interval).
Together, the two-site composite and its accompanying age-depth model are unequivocal about the astronomical origin of the meter-scale lithological cycles: They are chiefly related to precession (∼20 kyr periodicity = ∼50 cycles/Myr frequency). Indeed, the log(Ca/Fe) multi-taper method spectra display dense clusters of significant frequencies between 40 and 60 cycles/Myr, corresponding to precession (Figures 7e-7h). A smaller cluster also occurs at the obliquity frequency (∼25 cycles/Myr), but it is clearly subordinate to the precession cluster. The δ 18 O and δ 13 C NAFF analyses, on the other hand, show spectral peaks that can be related to obliquity (∼25 cycles/Myr), short eccentricity (∼10 cycles/Myr) and long eccentricity (∼2.5 cycles/Myr) (Figures 8c  and 8d). The clear-cut amplitude modulation patterns in Ca/Fe have also been subjected to NAFF analysis and are marked by frequencies that are related to short and long eccentricity, in agreement with astronomical theory (Figure 8b). Based on these elements, we rebut the original obliquity interpretation (Boulila et al., 2018;Boulila & Hinnov, 2022;Vahlenkamp et al., 2018), which was erroneously arrived at for two main underlying reasons. First, the common occurrence of hiatuses in the studied Newfoundland drifts system makes for a complex stratigraphic system. When hiatuses remain undetected, sedimentation rates are underestimated and cycle durations overestimated. Second, Fourier-Transform-based power spectra of Ca/Fe depth-series are dominated by a single peak (or cluster of peaks) that is related to the basic lithological rhythm (Figure 7). Additional peaks at lower frequencies are absent or modestly developed. The latter hampers the application of the classic cyclostratigraphic frequency-ratio method (Sinnesael et al., 2019). There is, of course, the clear bundling of the lithological cycles into groups of 3-5, but this bundling is not sufficiently distinct to discriminate between eccentricity-modulated precession (1:5 ratio) or the amplitude modulation of obliquity by the ∼173-kyr s 3 -s 6 term (1:4 ratio).
The precession interpretation sheds new light on the mechanistic pathway between astronomical insolation forcing and drift deposition, as well as on sedimentation rates at Sites U1408 and U1410.  is associated with strong cooling of surface waters in the Greenland-Norwegian Sea during precession maxima (when Earth is in the aphelion during northern hemisphere summer). This in turn allows for a vigorous Deep Western Boundary Current, transporting more clay from sources along the northeast Canadian margin to Newfoundland Ridge, leaving behind a low Ca/Fe ratio in the sedimentary record. The precession interpretation also yields sedimentation rate estimates twice as high as previously thought, and thus an exceptionally high time-resolution for the multi-proxy datasets.

Numerical Analysis of Fundamental Frequencies
To distil astronomical components from the Newfoundland Ridge geologic data, we adopt a strategy that is largely similar to the strategy applied by Olsen et al. (2019). The main difference between Olsen et al. (2019) and this work is that the U1408/U1410 age model was obtained by aligning the benthic isotope records with CENOG-RID (Figure 6), whereas Olsen et al. (2019) tuned their time series using a pure 405-kyr g 2 -g 5 cosine model. Simple age-depth model for the Middle Eocene composite section from Newfoundland Ridge. Twenty tie-points relate the two-site composite depth scale (m CCSF-X) to geologic age. In-between tie-points linear interpolation is adopted. The age-depth model was constructed by aligning the composite's isotope records with the CENOGRID stratigraphic backbone (Westerhold et al., 2020). Thereby, we considered all available magneto-and biostratigraphic constraints (from both Sites, Table 1), while we explicitly did not assume the sedimentary composite to be time-continuous throughout. The drift Sites U1408-U1410 are characterized by exceptionally high sedimentation rates in comparison with the pelagic reference sites (Walvis Ridge Site 1263 for this part of the Eocene) incorporated into CENOGRID (Westerhold et al., 2020).
The Olsen et al. (2019) tuning completely excludes the inheritance of g-term variability from the astronomical solution, which is not true for our approach. However, the conservative tuning with only 20 tie-points is designed to restrict possible inheritance effects to a minimum (through CENOGRID which is in turn tuned to eccentricity solution La10b). For this study, we deliberately adopt the 20-tie-points age model for NAFF analysis because a tuning of the U1408/U1410 records to a pure 405-kyr cosine model (as in Olsen et al., 2019) led to unsatisfactory age models. The NAFF analysis presented here thus represents a balancing between geochronological accuracy and possible inheritance.
The NAFF method extracts the dominant periodic components from a data series and sorts them by decreasing amplitude. Hence, NAFF provides a tool to objectively detect consistent periodic components in a signal, even when they are buried in noise. An additional convenience of NAFF is that multiple proxies can be used in tandem. This is useful when different astronomical parameters are best registered in different proxies. Here, we apply NAFF to the time-continuous upper part of the middle Eocene composite: between 39.47 and 42.81 Ma for the Ca/Fe time-series, and between 39.47 and 40.82 Ma for the isotopic time-series. The first 70 NAFF components (23 components for the Ca/Fe precession envelope) are shown in Figure 8, but it becomes readily clear that periodic components reminiscent of the periodic components of Earth's orbital eccentricity and spin-axis obliquity and precession occur within the first few high-amplitude components. We employ the NAFF result of the Ca/Fe time-series to select periodic components that correspond to the 5 main precession components (p + g 1-5 ). Similarly, the NAFF result of the Ca/Fe precession envelope is used to select periodic components that correspond to various eccentricity components. The long 405-kyr eccentricity component g 2 -g 5 constitutes an exception to this  (Katz et al., 2003). The curve in the background of our benthic isotope time-series is the Cenogrid stratigraphic reference (Westerhold et al., 2020). The eccentricity curve in the background of the Ca/Fe time-series is the La11 solution (Laskar et al., 2011b). rule though: While a periodic component associated with g 2 -g 5 could be discerned in the NAFF result of the Ca/ Fe precession envelope, this component is even more prevalent in the NAFF result of the δ 18 O time-series. Moreover, the frequency of the g 2 -g 5 component in δ 18 O (2.49 cycles/Myr) is closer to the g 2 -g 5 frequency predicted by astronomical models (2.46 cycles/Myr), compared to the g 2 -g 5 component in the precession amplitude envelope of Ca/Fe (2.85 cycles/Myr). It should be noted that an imprint of g 2 -g 5 eccentricity can be discerned in all four NAFF results in Figure 8. Finally, obliquity-related components (p + s 2-4 and p + s 6 ) are extracted from the δ 13 C NAFF result.
We thus follow a strategy of extracting individual astronomical parameters (precession, obliquity, short eccentricity and 405-kyr eccentricity) from the proxy that most clearly carries its signal. Furthermore, there is also a paleoclimate rationale behind this strategy: precession is the driving force behind local changes in Deep Western Boundary Current velocity, with Ca/Fe as a proxy. Hence, we reconstruct precession by means of Ca/Fe. The impact of precession on insolation and climate, in turn, is amplitude-modulated by eccentricity. Therefore, it is logical to determine the 100-kyr eccentricity component by using the Ca/Fe precession envelope. An obliquity imprint in the U1408/U1410 composite is discernible in the Ca/Fe record, as well as in both isotope records.
Here, we work with δ 13 C benthic to reconstruct obliquity because it displays the full suite of obliquity frequencies with clear high-amplitude peaks, and because it is not affected by precession-obliquity interference peaks.  astronomical components (red bars) allow the reconstruction of fundamental astronomical frequencies at ∼41 Ma (i.e., the mid-point of the analyzed interval), just by assuming that the outer Solar System is stable over the age of the Earth. Concretely, we assume the g 5 and s 6 frequencies to be invariant through geologic time, and adopt their present-day values as reported in Laskar et al. (2004). These fundamental frequencies are related to Jupiter and Saturn, respectively, and their assumed stability is due to the large mass of these planets (317.8 and 95.2 times the mass of the Earth). With these assumptions, we first calculate g 2 from the δ 18 O-extracted g 2 -g 5 periodic component. Subsequently, we calculate g 1 from g 2 -g 1 . After these two steps, we obtain an overdetermined system and can directly calculate g 3 from g 3 -g 5 as well as from g 3 -g 2 . The same goes for g 4 , which is calculated from g 4 -g 5 , g 4 -g 2 , and (g 4 -g 2 )−(g 2 -g 5 ). The resulting g i frequencies are compared to the La04 astronomical model, as well as to an independent reconstruction from Walvis Ridge early Eocene data by Meyers and Malinverno (2018). Our results for g 4 are consistent with La04, but our estimates for the g 1 and g 3 frequencies are up to 4% lower than the minimum value predicted by the La04 model, and our g 2 frequency estimate is ∼0.5% higher. Interestingly, Meyers and Malinverno (2018) observed an analogous data-model mismatch. These authors explain the mismatch by pointing out that Laskar et al. (2004) adopts 20-Myr averaging intervals before plotting the variation in secular frequencies g 1-4 , whereas the geology-based reconstructions for the Eocene span much shorter time intervals.
The NAFF result for g 3 is remarkable because a 4% lower g 3 (∼16.705"/year) in combination with a roughly constant g 4 (∼17.769"/year) would imply a g 4 -g 3 eccentricity term of ∼1.064"/year. This frequency is equivalent to a g 4 -g 3 eccentricity term with a periodicity of 1.22 Myr. In the present day, the g 4 -g 3 term is an eccentricity cycle with a periodicity around 2.4 Myr (∼0.544″/year). Together with the 1.2-Myr-long s 4 -s 3 obliquity modulation cycle, the g 4 -g 3 eccentricity term is one of the "Grand Cycles" in cyclostratigraphy (Hinnov, 2013). However, this 2:1 ratio between g 4 -g 3 and s 4 -s 3 may not have been steadfast through geologic time. Astronomical models predict that the g 4 -g 3 term can quickly transition from its "normal" frequency ∼0.544″/year to a value >0.9″/ year, creating a new 1:1 resonance between the two "Grand Cycles" (Hoang et al., 2021). In other words, the g 4 -g 3 periodicity has the potential to temporarily change from a ∼2.4 Myr period (libration) to a ∼1.2 Myr period (circulation). While theoretic astronomical models do predict that such g 4 -g 3 transitions occurred in the geologic past, they consider it unlikely that this would have occurred in the last 50 million years (Hoang et al., 2021). That being said, the exact timing of such chaotic transitions remains an open question. This is because, on the one hand, astronomical models are strongly dependent on initial conditions, and on the other hand, it is difficult to extract these chaotic transitions from the geologic record. To date, there are only two publications that report on tentative geologic indications of chaotic transitions around 52 Ma (Westerhold et al., 2017) and 87 Ma (Ma et al., 2017). In case our reconstructed values for g 3 and g 4 are accurate, we provide a third possible timing for a chaotic resonance transition in Earth's history, around the Middle Eocene Climatic Optimum. In this study, the robustness of this interpretation is further scrutinized with TimeOptMCMC (see Section 5.2).
The middle Eocene precession constant is estimated in nine different ways, using five different precession arguments (p + g 1-5 ), one obliquity argument (p + s 6 ), and multiple estimates for g 3 and g 4 ( Table 3). The nine estimates range between 50.088″/year and 51.62″/year, with a median value of 51.20″/year and an interquartile range between 50.88 Figure 8. Extraction of fundamental astronomical frequencies with NAFF. Precession (p + g i ), obliquity (p + s i ) and eccentricity (g i -g j ) arguments are differently expressed in the different proxies. We use log(Ca/Fe) to extract precession, δ 13 C for obliquity, the precession envelope of Ca/Fe for short eccentricity and δ 18 O for long eccentricity. Selected NAFF frequencies (bar chart) are indicated in red and labeled with their associated astronomical argument (listed in Table 3). Multi-taper method spectra are shown in the background. and 51.42″/year (boxplot in Figure 9e). The spread between the nine precession constant p estimates gives a sense of the uncertainty related to its reconstruction, yet it does not constitute a full-blown error propagation analysis. Uncertainties arise from several assumptions that went into age-depth modeling, while associating NAFF frequencies to astronomical parameters also requires user decisions. These uncertainties are hardly quantifiable and the propagation of uncertainties into NAFF results is not attempted in the framework of this work. Instead, the TimeOptMCMC technique is adopted in §5.2. This algorithm extracts the g-terms and precession constant p from the U1408/U1410 composite, while simultaneously estimating their uncertainties with a Bayesian Monte Carlo approach.

TimeOptMCMC
The NAFF technique has the advantage of being intuitive. It also allows working with multiple proxies when different astronomical parameters are best recorded by a range of different proxies. Major drawbacks for the NAFF technique reside in the fact that it is difficult to quantify uncertainties and the selection of astronomical frequencies depends on expert judgment. To circumvent both disadvantages and to scrutinize the NAFF results, we apply TimeOptMCMC on the time-continuous log(Ca/Fe) time-series between 39.47 and 42.81 Ma (allowing for time-axis perturbations of ±10%). The precession-dominated log(Ca/Fe) time-series and its strong amplitude modulation are well-suited for the TimeOptMCMC approach because this method is simultaneously evaluating the concentration of spectral power at the expected astronomical frequencies, and the amplitude modulation of precession by eccentricity. First, we utilized prior distributions (gray distributions in Figure 10) that have identical mean values compared to the Eocene Walvis Ridge analysis in Meyers and Malinverno (2018). The standard deviations of the prior distributions of the g-terms are however twice as large as in Meyers and Malinverno (2018). This choice was made to give TimeOptMCMC the freedom to test astronomical configurations that deviate further from the nominal astronomical solution, like for example, the 4% lower g 3 frequencies suggested by NAFF.
The posterior distribution for g 1 is bimodal, with the lower mode being in agreement with the NAFF result ( Figure 10a). The g 2 posterior distribution is similar to the prior distribution. The NAFF g 2 estimate occurs within the high probability range of the Bayesian approach (Figure 10b). The similarity between g 2 prior and posterior distributions is an important observation, as it underlines the stability and invariability of the g 5 -g 2 405-kyr eccentricity component. Indeed, our results once again affirm the status of the 405-kyr eccentricity term as the prime astronomical metronome for geologic time-keeping. The g 3 posterior distribution is multimodal with the first mode occurring close to the predicted g 3 frequency in the La04 astronomical solution (Figure 10c). The second mode, however, is in good agreement with the 4%-lower g 3 result obtained through NAFF analysis. We note that this second mode is at the lowermost end of the prior distribution, which implies that such low g 3 frequencies are relatively underexplored by the TimeOptMCMC algorithm. Nevertheless, the second g 3 mode illustrates that, when the algorithm examines relatively low g 3 frequencies, high likelihoods are obtained. The g 4 posterior distribution is significantly narrower than the prior distribution with two modes close to the predicted value in the La04 solution ( Figure 10d). We came to a similar observation based on the g 4 NAFF results and we thus consider the g 4 results of both techniques to be well-aligned. The same goes for the NAFF and TimeOptMCMC results for the precession parameter p (Figure 10e).
By plotting the posterior distribution of the g 4 -g 3 term in the TimeOptMCMC analyses (Figure 10f), we examined a possible chaotic resonance transition around 41 Ma as it was suggested by NAFF analysis. The g 4 -g 3 term in the TimeOptMCMC simulations is bimodally distributed, with a first mode around 2.75 Myr periodicities and a narrow but dense mode around 1.25 Myr. The TimeOptMCMC simulations thus provide support for a possible chaotic resonance transition between libration (∼2.4 Myr g 4 -g 3 period) and circulation (∼1.2 Myr g 4 -g 3 period) around 41 Ma.
The NAFF and TimeOptMCMC results are broadly compatible in the sense that they both hint at the possibility for a relatively low g 3 frequency around 41 Ma. To explore this possibility further, we updated our prior beliefs on the g-term frequencies with new information observed through NAFF analysis (Table 3). We reran TimeOptMCMC, now shifting the average values for the g 1 to g 4 fundamental frequencies to 5.2649, 7.4905, 16.705, and 17.769″/ year, respectively. These values represent the mean values of the updated prior distributions, depicted in gray on Figure 11. The resulting posterior g-term and precession constant p posterior distributions are shown as red distributions in Figure 11, and amount to 5.2956 ± 0.2485"/year for g 1 , 7.4828 ± 0.0289"/year for g 2 , 16.6421 ± 0.2786"/ year for g 3 , 17.6000 ± 0.2608"/year for g 4 and 51.2805 ± 0.5564"/year for the precession constant (uncertainties reflect ±1σ). We consider these values to be the best estimate of these frequencies from the between 39.5 and  Table 3). (f) The distribution of the very-long (g 4 -g 3 ) eccentricity period in the TimeOptMCMC simulations shows a bimodal distribution with one population indicating libration (∼2.75 Myr) and another population indicating circulation (∼1.2 Myr).
42.8 Ma dataset studied, and thus serve as a possible target for subsequent generations of astronomical solution at that time (Table 4, Figure 12). Moreover, the posterior g-term and precession distributions ( Figure 11) have single or multiple modes that are in excellent agreement with the NAFF reconstructions (Figures 9 and 12). We note that such agreement is to be expected as both methods were applied to the same proxy time-series. Yet, this is the first time the two techniques are compared, and confirmed to be compatible despite their major methodological differences. Thus, our conclusion that the evolution of the fundamental frequencies through geologic time exhibited much greater Myr-scale variability than previously assumed is corroborated by both NAFF and TimeOptMCMC results.

Implications for Earth-Moon Dynamics
Our precession constant estimates are markedly lower compared to the corresponding value in the La04(1,1) astronomical solution (Figure 12e). Yet, it is in excellent agreement with the Waltham (2015) calculations that include time-fluctuating tidal friction. The early Eocene precession constant reconstruction based on Walvis Ridge data (Meyers & Malinverno, 2018) provides additional support for the Waltham tidal friction model: Both their and our Eocene reconstructions suggest that tidal friction is overestimated Figure 11. Summary of TimeOptMCMC results, ran with the updated prior belief that the g 3 frequency was ∼4% lower around 41 Ma. (a-e) Prior (gray distribution) and posterior (red distribution) for the different g-terms and precession constant are compared to the results of the NAFF approach (red diamonds, see also Table 3). (f) The distribution of the very-long (g 4 -g 3 ) eccentricity period in the TimeOptMCMC simulations shows a bimodal distribution with the main mode indicating circulation (∼1.2 Myr). in the nominal La93(1,1) and La04(1,1) solutions (Laskar, Joutel, et al., 1993;Lasker, Robutel, et al., 2004), at least when going back into the Paleogene. When only these two Eocene data points are considered, we find a much better model-data fit when considering an astronomical solutions with a tidal dissipation that is considerably lower (∼20% lower) compared to the present-day value (as used in La93(1,1); see Materials and Methods; Figures 9e and 12e). However, geology-based reconstructions of tidal dissipation and dynamic ellipticity for the Neogene (Lourens et al., 2001;Pälike & Shackleton, 2000;Zeeden et al., 2014) demonstrate good agreement with the evolution of the precession constant as implemented in the nominal Laskar solutions and the Farhat et al. (2022) model ( Figure 12e). Therefore, we infer that the tidal dissipation of rotational energy must have occurred at a relatively low pace throughout the Paleogene, in agreement with numerical tidal models (e.g., Green et al., 2017), after which a marked increase in tidal drag caused a much more rapid decrease of the precession constant during the Neogene and Quaternary.
The methodology we present here is unique in that it starts from a conservative age-depth model that was constructed prior to, and independent of, the extraction of astronomical components. While we find larger-than-expected variability in astronomical g-terms, our tidal dissipation reconstructions support similar results reported by previous studies (Meyers & Malinverno, 2018;Zeebe & Lourens, 2022). Therewith, our results urge Paleogene workers to move on from the default (1,1) setting for dynamical ellipticity and tidal dissipation in astronomical solutions when assessing detailed obliquity-precession interference patterns.

Conclusion
The extraction of astronomical components presented here exploits the particularly cyclic middle Eocene drift deposits on Newfoundland Ridge, cored during IODP Expedition 342. The lithologic clay-content cycles are now firmly identified as the imprint of climatic precession. This interpretation was made based on a carefully constructed two-site composite, in combination with an age-depth model that solely consists of only 20 age-depth tie-points. This feature of our analysis allows the classic cyclostratigraphic approach to be reversed while minimizing the risk for circular reasoning: Different rhythmic components in high-resolution proxy-series served the reconstruction of four g-terms and the precession constant p. These reconstructions provide novel constraints on the Cenozoic evolution of our solar system. First, the variability in g-term frequencies on million-year timescales has previously been underestimated. Second, the internally consistent evidence for a relatively slow tidal energy dissipation throughout the Paleogene and a strong increase during the Neogene caps a long-standing debate. Both pieces of information (g-terms and p, Table 4 and Figure 12) constitute targets that astronomers can use to extend the reliability of astronomical insolation models in geologic time.