A punctuated equilibrium analysis of the climate evolution of cenozoic exhibits a hierarchy of abrupt transitions

The Earth’s climate has experienced numerous critical transitions during its history, which have often been accompanied by massive and rapid changes in the biosphere. Such transitions are evidenced in various proxy records covering different timescales. The goal is then to identify, date, characterize, and rank past critical transitions in terms of importance, thus possibly yielding a more thorough perspective on climatic history. To illustrate such an approach, which is inspired by the punctuated equilibrium perspective on the theory of evolution, we have analyzed 2 key high-resolution datasets: the CENOGRID marine compilation (past 66 Myr), and North Atlantic U1308 record (past 3.3 Myr). By combining recurrence analysis of the individual time series with a multivariate representation of the system based on the theory of the quasi-potential, we identify the key abrupt transitions associated with major regime changes that separate various clusters of climate variability. This allows interpreting the time-evolution of the system as a trajectory taking place in a dynamical landscape, whose multiscale features describe a hierarchy of metastable states and associated tipping points.

Early evidence of abrupt transitions in Camp Century and Dye 3 Greenland ice cores 1,2 attracted a lot of attention from the paleoclimatic community before being well acknowledged and understood.These findings gave evidence of a sequence of previously unknown abrupt climatic variations.However, such transitions did not seem to agree with other marine and terrestrial records, which led to considerable debate in the field [3][4][5] .Over the course of several decades spent retrieving and studying more detailed paleorecords, the existence of these rapid climatic variations, known as Dansgaard-Oeschger events (DO), has become well accepted.Further support for these findings has recently been provided by the identification of additional abrupt transitions from the NGRIP ice core, which has been made possible by the increased temporal resolution of the record 6 .These additional events correspond to changes of either short duration or amplitude in the stable O 18 over O 16 isotopes ratio δ 18 O that visual or standard statistical inspections do not necessarily flag.We remind that the fractionation of oxygen isotopes can be used to reconstruct local temperature conditions.The Earth climate has experienced numerous abrupt and critical transitions during its long history, well beyond the specific examples above 7,8 .As discussed in Rothman 9 following Newell 10 , rapid variations change are more likely to lead to catastrophic consequences in the biosphere -as in the case extreme case of mass extinction eventsbecause it is hard for the evolutionary process to keep pace with shifting environmental conditions.Such transitions are often referred to as climatic tipping points (TPs), associated with possibly irreversible changes in the state of the system.The term TP was originally introduced in social sciences 11 and made popular more recently by Gladwell 12 .The study of TPs has recently gained broad interest and perspective in Earth and Environmental sciences, especially with regard to the future of our societies under the present climate warming scenarios [13][14][15][16] .The term Tipping Elements (TE) was introduced in 15,17 , and subsequently adopted by other authors 14,18 to characterize the specific components of the Earth System that are likely to experience a TP in the near future as a result of the ongoing climate crisis 16,[19][20][21] .Recently, the concept of TP has been used to define, in turn, rapid societal changes that might lead to potentially positive impacts in terms of addressing the climate crisis 22,23 .
Here, we want to investigate critical transitions in the Earth climate history by looking at two key high-resolution datasets that show evidence of abrupt transitions.The first dataset is the CENOGRID benthic δ 18 O and δ 13 C record corresponding to the compilation of 14 marine records over the past 66 Myr, from the Cretaceous-Paleogene (K-Pg) extinction event till present 24 .The second dataset comprises the North Atlantic U1308 benthic δ 18 O, δ 13 C and δ 18 O bulk carbonate time series covering the past 3.3 Myr 25 .While visual evidences of abrupt transitions have already been discussed for these datasets, we wish to identify key abrupt thresholds by applying the recurrence quantification analysis (RQA) to each individual univariate time series and supplementing it with the Kolmogorov-Smirnov (KS) test 26 , see 27 and discussion below.Then, the selected transitions are discussed in the context of the Earth climate history allowing the definition of dynamical succession of abrupt transitions.Such transitions are then interpreted taking into account the evolution of key climate factors such as CO 2 concentration, average global sea level, and depth of the carbonate compensation.
9][30][31] and discussion in 32,33 .The multistability of the climate system comes the presence of more than one possible climate states for a given set of boundary conditions 34 .While earlier analyses have mostly evidenced the possibility of bistable behaviour, multistability can indeed include multiple competing states [35][36][37][38] .Recently, it has been proposed that the metastability properties of the climate system can be understood by interpreting the climate evolution as a diffusion process taking place in an effective dynamical landscape 37,39,40 defined by the the Graham's quasi-potential 41 .The local minima of the quasi-potential indicate the competing metastable states, with the transitions between such states occurring preferentially through the saddles of the quasi-potential.Such a viewpoint mirrors earlier proposals for interpreting biological evolution, namely the Waddington's epigenetic landscape [42][43][44][45][46] , see also relevant literature associated with synthetic evolution models like Tangled Nature 49,50 , and foresee the climatic transitions associated with the TPs as a manifestation of a dynamics characterized by punctuated equilibria 47-50, The interplay between periods of stasis and rare transitions between competing metastable states seem to reflect the fact that climate fluctuations behave according to the dominance of stabilizing vs destabilizing feedbacks when considering short vs long time scales, respectively 51 .We will see that in most cases the RQA applied to the CENOGRID δ 18 O dataset flags TPs that disagree in terms of dating from those derived from the δ 13 C dataset.This is unsurprising because the two proxy data are sensitive to vastly different climatic processes.Nonetheless, the candidate TPs from both datasets come hand in hand with saddles of the bidimensional quasipotential estimated from the bivariate time series.This indicates consistency between separate ways of detecting TPs.Additionally, TPs featuring faster characteristic time scales are associated with smaller-scale decoration of the quasi-potential, in agreement with what conjectured in 37 .The analysis of the benthic δ 18 O and δ 13 C suggests that the evolution of the climate in Cenozoic is characterized by a hierarchy of TPs due to an underlying multiscale quasi-potential.

Detecting Critical Transitions of the past 66 Myr -3 Myr history of the Earth Climate
The augmented KS test of the benthic δ 18 O record of the past 66 Ma identifies six major abrupt transitions corresponding to two major warming events at about 58 Ma and 56 Ma, followed by four major coolings at 47 Ma, 34 Ma, 14 Ma and 2.8 Ma respectively (Fig. 1A).The competing metastable states associated with these transitions feature rather long temporal persistence (a long time-window of 1-4 Ma is used, see Suppl.Mat.).These events are classical ones described from the literature 52 , where the first two transitions led to warmer conditions, while the latter four led to colder conditions.The same transitions are identified by employing the recurrence plot (RP) and recurrence rate (RR) analyses 27,53 , which also identify three more events, occurring at around 63 Ma, 40 Ma and 9.7 Ma, see Fig. 1C, Suppl.Tab. 1.We have chronologically labeled these TP O 1 to TP O 9, where the lower index refers to the used proxy data.e older macrocluster corresponds to very warm climate conditions and, using Marwan et al.'s nomenclature 53 , a disrupted variability.The average global temperature was estimated to be between 8°C and 16°C above the present day one, with no apparent presence of any major continental ice bodies 24 (Fig. 1A).Five major transitions are found within this period.They include the 2 major abrupt warmings at 58 Ma (TP O 2)    climates in the earlier, warmer period, followed by the "Coolhouse" (TP O 5-TP O 6) and "Icehouse" climates (TP O 6 to present); see Fig. 1.The first two states alternated in a warm-hot-warm sequence under extremely high CO 2 concentrations 65 as compared to those measured over the past 800 Kyr in the Antarctic ice cores, which represent the reference states for the IPCC potential scenarios of climate change warming 67 (Fig. 2C).Before 34 Ma, one finds substantially larger values for GMSL, CCD depth, and CO 2 concentration whose average values are: +38 ± 15 m, 4600 ± 150 m, and 630 ± 300 ppmv respectively (Suppl.Tab. 2).Note that the CO 2 concentration is much higher than the levels observed over the past 800 kyr in the Antarctic ice cores, which represent the reference states for the potential scenarios of climate warming outlined by the IPCC 67 .
Conversely, from 34 Ma until present, the Earth experienced much lower CO 2 concentrations and GMSL, thus generating the classical climate trend towards the recent ice-age conditions 8,24,52 (Fig. 2A,C).Indeed, the last 34 Myr show average values in GMSL, CO 2 concentration, and CCD depth of -3.5 ± 13 m, 330 ± 160 ppmv, and 3500 ± 400 m, respectively (Suppl.Tab. 2), which are much lower than in the older interval.This second set of means has underestimated values since it is based on data from the Miller et al. 58 dataset, which ends at 0.9 Ma, These GMSL, CCD and CO 2 reconstructions show key transitions that fit with the CENOGRID thresholds deduced from the KS and RR analysis of the benthic δ 18 O as they signal an increase or a decrease in the global mean sea level corresponding roughly to warming or cooling episodes of the Earth history or strong variations in the concentration of atmospheric CO 2 .The variations observed in CO 2 , CCD and GMSL at the 9 identified TPs from the benthic δ 18 O record indicate heterogeneous characteristics prior to the TP6 major threshold (Suppl.Tab. 3).On the contrary, more homogenous features are noticed in the three climate proxies after TP O 6, translating the occurrence of major reorganizations in the climate system, which become interesting to test at shorter timescales; see discussion below.

Quasi-potential Landscape and Critical Transitions
As discussed in detail in the Methods section, taking inspiration from the use of the Waddington epigenetic landscape to describe evolution [42][43][44][45][46] and from the theory of punctuated equilibrium 49,50 , it has been proposed in 37,39,40 to study the global stability properties of the climate system by introducing an effective quasi-potential 41 , which generalizes the classical energy landscape framework often used to study metastable stochastic system.The quasi-potential formalism allows one to study general nonequilibrium system and to associate local maxima of the probability distribution function (pdf) of the system stable states.Additionally, saddles of the pdf are associated with Melancholia (M) states 34 , which are unstable states living in the boundary between different basins of attraction.In the weak-noise limit, noise-induced transitions between competing stable states are expected to go through such M states.What we have done here is to construct the empirical bivariate pdf of climate system in the projected (δ 13 C, δ 18 O) and check whether TP O 1-9 indicated in Fig. 1 correspond, at least approximately, to saddles of the pdf.We find -see Fig 3A -that, indeed, this is the case for TP O s 1, 2, 4, 5, 6, 8, 9, whereas no agreement is found for TP O s 3 and 7, which seem to take place in regions where the density is very small and very large, respectively.A very similar version of the bivariate pdf shown here had already been used in Westerhold et al. 24 to identify the Icehouse, the Coolhouse, the Warmhouse, and the Hothouse states, as groupings of nearby peaks corresponding to qualitatively similar climates.
A major improvement we propose here is to identify the transitions between such states, as well as less obvious transitions, by looking at the saddle points.TPo6 clearly emerges as the most important transitions, as it basically breaks the pdf into two separate parts.Note that the transitions shown in Fig. 1 are associated, because of the choice of a long time-window, with events that occur over long periods and that lead to persistent changes in the state of the system.Two key climatic features that appear as TPs in both records: the PETM is analogous to TP3 in both records, while TP O 6 and TP C 8 represent the EOT.

The recent past
The past 3.3 Myr record from North Atlantic core U1308 can be considered as a blowup of the CENOGRID dataset (Fig. 4).As previously mentioned, the last 3.3 Myr have been defined as an Icehouse climate state, with the appearance, development, and variations of the NHIS 24 , whilst the Antarctic ice sheets had already mostly reached their maximal expansion.The variations in the deep-water temperature, as expressed by the benthic δ 18 O, are interpreted as an indicator of the continental ice volume with clear interglacial-glacial successions [68][69][70] .The Icehouse state is characterized by a change of the interplay between benthic δ 13 C and δ 18 O, which corresponds to a new relationship between the carbon cycle and climate 71 .Indeed, one finds a very strong correlation between the two records (Pearson's coefficient being approximately -0.6).The correlation mainly results from the fact that the time series have approximately a common quasi-periodic behavior due to amplified response to the astronomical forcing as dictated by the Milankovich theory., and therefore illustrates the dynamics of the Northern Hemisphere ice sheets (NHIS), yields similar dates to those obtained for the benthic δ 18 O record (see Suppl.Fig. 3, Suppl.Tab. 1) 72 .Indeed, one finds abrupt transitions at 2.75 Ma, at 1.5 Ma, at 1.25 Ma, at 0.9 Ma and 0.65 Ma.Finally, as opposed to the CENOGRID δ 13 C RP, U1308 δ 13 C RP shows a drifting pattern similar to that of benthic δ 18 O, with only 2 key transitions at 2.52 Ma and 0.48 Ma (see Suppl.Fig. 3).Note that the 0.48 Ma transition does not have any equivalent in the benthic δ 18 O records.

Discussion
Studying the same CENOGRID dataset, Boettner et al. 78 78 , which is instead absent in the case for the EOT key transition.Based on the results of both the RQA and the KS test of the δ 18 O and δ 13 C time series considered in this study, and the bivariate analysis performed using the framework of the quasi-potential theory, we propose a succession of critical transitions as described in Fig. 5.The critical transitions TP O 1 to TP O 9 shaped the Earth climate towards the onset and development of the Southern ice sheets and the later build-up of the NHIS.TP O 9 is followed by four more recent RTPs during the Quaternary, which steered the evolution of the ice sheets and of the climate as a whole until the present day.The climatic evolution during the Cenozoic until about 3 Ma seems to conform to a punctuated equilibrium framework, where the TP O s are associated with rapid transitions between rather different metastable modes of operation of the climate system.In particular, a key step that separates two rather diverse sets of climatic states occurred around 34 Ma at the EOT (TP O 6).Without the major drop in GMSL, in CO 2 concentrations, and in CCD, the Earth climate could have been different.However, after TP O 6, the Earth climate entered new dynamical regimes marked by much lower CO 2 concentrations, a lower GMSL, and a lower CCD.The remodeling of oceanic basins and mountain uplifts changed the marine and atmospheric circulations patterns, which played a crucial role in initiating and shaping the development of the NHIS.
Interestingly, the analysis of the δ 13 C time series identifies a different set of critical transitions, with the exception of the PET (TP O 2 and TP C C) and the EOT (TP O 6 and TP C 8), which are identified for both analyzed proxies.Looking into the bivariate pdf in the projected (δ 13 C, δ 18 O) space allows a better understanding of the nature and the origin of the TPs separately detected by studying the recurrence properties of the univariate time series.Indeed, we are in most cases able to associate both the TP O s and the TP C s to transitions across saddles of the effective quasi-potential.Clearly, some critical transitions might be more easily detectable when examining one time series rather than the other, because different proxies might be more sensitive to the active climatic processes, yet the approach taken here allows placing all TPs within a common ground., Northern Hemisphere ice sheet maps are from Batchelor et al. 80 .The paleogeographic maps have been generated using the Ocean Drilling Stratigraphic Network (ODSN) plate tectonic reconstruction service: <https://www.odsn.de/odsn/services/paleomap/paleomap.html>.The red arrows on the tectonic maps indicate the key events that correspond to the identified abrupt transition.
More recently, variations in the extent and volume of the NHIS have contributed to the occurrence of the millennial variability marked by the Bond cycles, which have been best described during the last climate cycle.However, the onset of these cycles has been proposed to date back to 0.9 Ma 72 .Human activity is now rapidly pushing the Earth system towards the limits of its safe operating space associated with the occurrence of TPs; see also Rothman 9 for a geochemical perspective focusing on the ongoing perturbation imposed on the marine carbon cycle..This concern is supported by actual observations impacting numerous tipping elements (see Lenton 14 ) through very drastic tipping cascades 16 .This paves the way for a possible upcoming major transition, which might lead us to climate conditions that are fundamentally different from what has been observed in the recent or more distant past, and, at the very least, could bring us into a climate state with a much reduced or absent NHIS 38 .This potential major transition, leading to de facto irreversible changes for the climate and the biosphere, could become a boundary between the Cenozoic icehouse and a new, warmer and radically different climate state compared to Pleistocene conditions.

Cool -Ice house landscape
New Anthropocene landscape?

Recurrence Plots and Kolmogov-Smirnov Test
The Kolmogorov-Smirnov test is a robust method for accurately detecting discontinuities in a particular time series and is therefore a very precise way for determining the timing of abrupt transitions.Our method, described by Bagniewski et al. 26 is modified from the two-sample KS test and has been successfully applied to various geological time series [24][25]69 . Themethod uses the KS statistic, D KS , to compare two sample distributions taken before and after a potential transition point within a sliding window.If the samples do not belong to the same continuous distribution, i.e., D KS is greater than a predefined threshold, a transition point is identified.This classic two-sample KS test is augmented by additional criteria to refine the results and pinpoint significant abrupt transitions indicative of a true climatic shift.This involves discarding transitions below a rate-of-change threshold since smaller changes in the time series might be due to data errors or short-term variability within intervals shorter than the proxy record's sampling resolution.Furthermore, as the frequency with which the KS test detects transitions largely depends on the chosen window length, D KS is calculated for different window length, within a range that corresponds to the desired time scale at which a given paleorecord is to be investigated.The analysis starts by identifying transitions using the longest window, which has the largest sample size and thus carries greater statistical significance.Subsequently, the method incorporates transitions detected using shorter windows to capture transitions occurring on shorter time scales.For more details see Bagniewski et al. 26 , Rousseau et al. 72 .
To gain further insight into the climate story revealed by proxy records, we performed an analysis based on recurrence plots (RPs), first introduced by Eckmann 81 in the study of dynamical systems and later popularized in the climate sciences by Marwan et al. 53,82 .
The RP for a time series { x !: i = 1, … , N} is constructed as a square matrix in a Cartesian plane with one copy { x !} of the series on the abscissa and another copy { x !} on the ordinate, with both axes representing time.A dot is entered into a position (i, j) of the matrix when x ! is sufficiently close to x ! .For the details -such as how "sufficiently close" is determined -we refer to Marwan et al. 50.The square matrix of dots that is the visual result of RP exhibits a characteristic pattern of vertical and horizontal lines, indicating recurrences.These lines sometimes form clusters that represent specific periodic patterns.Eckmann 81 distinguished between large-scale typology and small-scale texture in the interpretation of RPs.The most interesting typologies in RP applications are associated with recurrent patterns that are not strictly periodic and are thus challenging to detect by purely spectral approaches to time series analysis.RPs offer an important advantage by enabling a visual identification of nonlinear relationships and dynamic patterns that characterize high-dimensional systems subject to time-dependent forcing, such as the climate system.However, in periodically forced systems, the recurrence structure is dominated by repetitive, periodic patterns rather than the underlying dynamics.This may pose challenges when attempting to identify meaningful transitions or regime shifts in climate at time scales strongly affected by orbital forcing.Marwan et al. 82 extensively discussed the various measures used to objectively quantify the typologies observed in RPs.Collectively known as recurrence quantification analysis (RQA), these measures include the Recurrence Rate (RR), which represents the probability of a specific state recurring within a given time interval.The RR is calculated by determining the density of recurrence points along the diagonal of the RP within a sliding window.Since low RR values correspond to an unstable behavior of the system, the minima of the RR may be used to identify abrupt transitions.Here we follow the approach by Bagniewski et al. 26 and select local minima based on their prominence.

Quasi-potential, Melancholia States, and Critical Transitions
Traditionally, tipping points are schematically represented as being associated with the bifurcation occurring for a system described by a one-dimensional effective potential when a change in the value of a certain parameter leads to a change in the number of stable equilibriums.Hence, conditions describing the nearing of a tipping point can be related to the presence of slower decay of correlations (critical slowing down).This viewpoint, while attractive, suffers from many mathematical issues due to the fact that the true dynamics of the system occurs in a possibly very high dimensional space.Tantet et al. 83,84 have introduced a mathematically rigorous framework for the occurrence of tipping points that clarifies the link between rate of decay of correlations, sensitivity of the system to perturbations, and robustness of the unperturbed dynamics.Here, we wish to take a different angle on the problem.Instead of focusing on the individual tipping points, we attempt to capture the global stability properties of the system.We take inspiration from the application of the Waddington epigenetic landscape to describe morphological evolution [42][43][44][45][46] and from the theory of punctuated equilibrium 49,50 , which associates periods of stasis (characterized by relatively stable morphology) with the emergence of new species through abrupt changes (named cladogenesis) from a previous species.Lucarini & Bodai 39,40 proposed describing the global properties of the climate system using the formalism of quasi-potential 41 .Roughly speaking, assuming that the system lives in ℝ !, and its dynamics is described by a stochastic differential equation the probability that its state is within the volume  around the poin  ∈ ℝ ! is given by  ,  =    where   ≈  !!! ! ! ! is the probability distribution function (pdf) and Φ  is the quasipotential.The function Φ  depends in a nontrivial way on the drift term and noise law defining the stochastic differential equation.This setting generalizes the classical energy landscape and applies to a fairly large class of stochastic dynamical systems.One can see the dynamics of the system as being driven towards lower values Φ  (plus an extra rotational effect that is typical of non-equilibrium systems), while the stochastic forcing noise makes sure that the system is erratically pushed around.Hence, the minima of Φ  correspond to local maxima of the pdf, and the saddles (which coincide for both Φ  and   ) coincide with the Melancholia (M) states 34 .Such M states are unstable states of the system that live at the boundary between basins of attraction and are the gateways for the noise-induced transitions between competing stable states.Margazoglou et al. 37 applied this method to the investigation of the metastability properties of an intermediate complexity climate model and suggested that the presence of decorations of the quasi-potential at different scales could be interpreted as being associated with a hierarchy of tipping points.Indeed, passing near M states is intimately associated with the occurrence of critical transitions.Hence, the construction of the quasi-potential Φ  can be seen as the structural counterpart of the investigation of the time-evolution of the system and of its critical transitions.Zhou et al. 85 give a complete overview of different methods applied to perform such an analysis.Global Mean Sea Level in meters from Miller et al. 58 .Identification of particular warm and of glaciation events.The Laurentide, GIW-WAIS and Ice free lines are from Miller et al. 58 ; b) Carbonate Compensation Depth (CCD) in meters from Pälike et al. 64 .The purple circles identify the TPs on this record; c) Estimate of the CO 2 concentration in parts per million volume (ppmv) from Beerling & Royer 65 ., Northern Hemisphere ice sheet maps are from Batchelor et al. 80 .
The paleogeographic maps have been generated using the Ocean Drilling Stratigraphic Network (ODSN) plate tectonic reconstruction service: <https://www.odsn.de/odsn/services/paleomap/paleomap.html>.The red arrows on the tectonic maps indicate the key events that correspond to abrupt transition.

A Punctuated Equilibrium Analysis of the Climate Evolution of Cenozoic exhibits a Hierarchy of Abrupt Transitions
Denis-Didier Rousseau, Witold Bagniewski, Valerio Lucarini

A Historical Account of the Critical Transitions
In addition to the Chicxulub meteor impact, which injected a considerable amount of CO 2 into the atmosphere 1,2 , Deccan traps were already spreading at the beginning of the Cenozoic, contributing to the release of massive amount of CO 2 3 .CO 2 concentrations continued to rise until reaching about 500 ppmv, coinciding with a very active period of the North Atlantic Igneous Province around 58 Ma -56 Ma (TP O 2 and TP O 3), which was associated with the opening of the North Atlantic Ocean 4 .During this time, the Northern Hemisphere plates were connected and did not experience the present Arctic conditions, allowing faunal and vegetal dispersion.In contrast, other plates were undergoing reorganization, such as India moving northeastward towards the Asian continent.Equatorial Pacific carbonate compensation depth (CCD) reached its minimum value a bit later at about 54 Ma when CO 2 concentrations were at a maximum for the whole Cenozoic, above 1,100 ppmv 5 .The decrease in GMSL and the CCD at TP O 4 (about 40 Ma) has been interpreted as the start of the icing of Antarctica through dating mountain glacier deposits using K-Ar dating of lava flows 6 , as well as other glacial evidence, i.e. from the Gamburtsev subglacial mountains 7 .During this period, the Northern Hemisphere plates remained connected.By TP O 4 India was approaching the Asian plate while Northern Hemisphere plates remained connected, allowing continental migrations of mammals at high latitudes 8,9 .Around 34 Ma, the EOT (TP O 6) is associated with the opening (in several steps) of the Drake and Tasmanian passages 10 , which led to a drastic change in the global ocean circulation.This resulted in a reduced strength of deep water formation in the Southern Hemisphere, and a decline in deep-sea temperatures, accompanied by a significant drop in relative sea level and a decrease in CCD (see Fig. 2).Based on paleoaltimetry estimates using oxygen isotopes, Rowley & Currie 11 indicate that the Tibetan Plateau had an elevation of about 4000m, favoring therefore the physical weathering of rocks, resulting in the consumption of CO 2 , and the enhanced burial of carbon through high sedimentation rates in the nearby seas.This may have contributed to the major threshold in the variation of the CO 2 concentration 12 , which also drops very strongly 5 (Fig. 2).Recent analysis of the evolution of sea surface temperatures in the North Atlantic Ocean during the Eocene-Oligocene Transition (EOT) suggests that the cooling during this period was primarily triggered by a decrease in CO2 concentration, with paleogeographic changes playing a secondary role 13 .Comparisons between marine and terrestrial data and models reveal a more abrupt decrease in pCO2 in marine proxies, as evidenced by the decrease in the CCD, contrasting with a more gradual decrease in terrestrial indicators.The gradual decrease observed in terrestrial proxies remains a topic of debate.Additionally, these studies suggest that paleogeographic changes played a significant role in initiating the Atlantic Meridional Overturning Circulation (AMOC), which became an important ocean circulation pattern in the growth of polar ice sheets (Hutchinson et al. 14 , and references herein).
As mentioned, The EOT transition is the major boundary between two different climate landscapes dominated by intensive plate tectonics and strong volcanism for the older one, and by major ice sheet coverage in both Hemispheres, with still very active plate tectonics (closure of seaways, orogenies), for the younger one 9 .Following the EOT, the climate experienced the build-up of the East Antarctic ice sheet (Oi-1 glaciation), which can be regarded as the onset of the cold world in which we are presently living.During this time, India has almost ended its transfer to the Asian plate.Between TP O 6 and TP O 7, i.e., 34 Ma and 14 Ma respectively, the East Antarctic ice sheet underwent waxing and waning with several major glaciations occurring before the 17 Ma to 14.5 Ma interval, characterized by rising sea levels, a severe shoaling episode of CCD, and higher CO 2 concentrations (Fig. 2).Such high CO 2 concentrations may have been fueled by the Columbia River major volcanism, which ended by TP O 7 at about 14 Ma 3 .Plate tectonics still played a very active role, with the closure of both the Indonesian gateway and the Tethyan seaway, and contributing to the start of the development of the Mediterranean 75 .Eurasia is now separated from North America and India colliding with the Asian continent, and the Andes are uplifting, thus modifying the geometry of the marine basins and the global oceanic circulation.West Antarctica is beginning to form while the East Antarctic ice sheet continues to strengthen and expand.TP O 8 at about 9 Ma sees a strong lowering of the GMSL and of the CO 2 concentrations (Fig. 2).TP O 9 is associated with a final major tectonic event corresponding to the closure of the Panama Isthmus, which, as a result of shutting down the exchange of water between the two oceans, led to re-routing large-scale flows in both the Pacific and the Atlantic Ocean, and configured a oceanic circulation that is very similar to present day 15 .This is associated with a significant decrease in the global sea level, a deepening of the CCD, and lower CO 2 concentrations (Fig. 2).These conditions will shape the climate history during the Quaternary, leading to the build-up of the Northern Hemisphere glaciers and ice sheets.It is worth noting that all the major transitions that have been identified during the interval between 66 Ma and 2.9 Ma -2.5 Ma are associated with the build-up, waxing, and waning of the Northern and Southern Hemisphere ice sheets.
above the 750 ppmv considered as the Antarctic glaciation threshold 3 , and marked by a relative minimum at about TP o 5 at 40 Ma (Fig. 3B,C).This first interval of the GMSL, which ends with a strong deepening of about 1000m, is associated with a strong decrease in CO2 concentration and an almost completely ice free Earth, with no major ice sheet in either the southern or northern hemisphere.This is deduced primarily from the high GMSL, mostly remaining above 12m (Fig. 2A), which would correspond to the mutual contribution of the Greenland and the West Antarctic ice sheets.Additionally, high CO2 concentrations estimates are in agreement with a CCD generally lower than 4000 m (Fig 2B,C).
The second main interval shows a completely different scenario with the GMSL varying between +30m and -80m without considering the late Quaternary interval, and much lower CCD and CO2 concentration (Fig. 2B,C).First the strong decrease in the GMSL at EOT reaches negative values of about -25m at about 33.5Ma.It is interpreted as evidence of the first continental scale Antarctic ice sheet (first Oligocene isotope maximum -Oi1 -Fig.2A) 4 .After a return to values similar to present mean sea level, about +2m, between 32 Ma and 30 Ma, a new decrease of about 25m in GMSL occurs between 29.5 Ma and 27 Ma corresponding to another continental scale Antarctic ice sheet extent labeled Oi2 (Fig. 2A).After a two-step increase in GMSL of about 40m between 27 Ma and 24 Ma, a new sharp decrease of about 40m is noticed at about 23 Ma.It corresponds to the Middle Oligocene Maximum (Mi1 -Fig.2A), another Antarctic ice sheet wide expansion 4,5 2A).The CCD indicates about 600 m shoaling which lasted around 2.5 Myr linked to high estimates of CO2 concentration from paleosols and stomata 2 (Fig. 2B,C).This strong increase corresponds to the Miocene Climatic optimum (MCO -Fig.2A) between 17 Ma and TP o 7 at 13.9 Ma, which is the last interval during which GMSL reaches such high values higher than +20m above the present day ones (Fig. 2).The interval between TP o 7 at 13.9 Ma and 13 Ma corresponds to the Middle Miocene transition during which GMSL once more decreases significantly by about 35m.
Such lowering is associated with the growth of the East Antarctic ice sheet to near its present state, remaining a perennial ice body that is thereafter impacting the Earth climate 6,7 .Although GMSL remains relatively stable between 13 Ma and 12 Ma, another strong decrease, again of about 30m, occurs between TP o 8 at about 9 Ma and 8.5 Ma, associated with the strongest deepening recorded by the CCD, around 4800m.GMSL increases again by about 20m until 7.5 Ma to remain relatively stable until 5.5 Ma when GMSL increases by about 20m between 5.5 Ma and 3.5 Ma, corresponding to the Pliocene Climatic Optimum (PCO -Fig.2A).Between 3.5 Ma until TP o 9 at about 2.7 Ma GMSL shows a sharp decreasing trend of about 35m (Mi2a, 3, 4 -Fig.2A).This is associated with the development of the large northern ice sheets between 2.9 Ma (TP o 9a) and 2.5 Ma (TP o 9b -Fig.2), especially the Laurentide ice sheet corresponding to about 50m decrease of GMSL with regards to the present day value (Fig. 2A).From TP o 8 at about 9 onward, CCD shows significant fluctuations, although it remains at around 4500m with two strong deepening events at about TP o 9a 2.9 Ma and TP o 9b (Fig. 2) at 2.5 Ma.
In another global sea level reconstruction 8 , several major thresholds were identified agreeing with the Miller et al 10 10 .These two key dates correspond to respectively to the first major iceberg calving in the Nordic Seas 11 and from the Laurentide ice sheet 12 although this interpretation was rejected by Naafs et al. 13 who instead attributed the IRD  ) detected by the RQA.

Figure. 1 .
Figure. 1. KS test and Recurrence Quantification Analysis (RQA) of CENOGRID benthic δ 18 O.a) Time series in Ma BP with difference of the reconstructed and present Mean Global Temperature in pink).KS test identifying abrupt transitions towards warmer conditions in red and cooler or colder conditions in blue; b) Recurrence plot (RP) with identification of the main two clusters prior and after 34 Ma.The main abrupt transitions identified are highlighted by red circles, and c) Recurrence rate (RR).The pink crosses and vertical green lines indicate the abrupt transitions (TP O ) detected by the RQA.CENOGRID benthic δ 18 O data are from Westerhold et al.24

Page 5 of 21 and 56
Ma (TP O 3).These two transitions are thresholds towards much warmer oceanic deep water characterizing the first late Paleocene-Eocene hyperthermal

Figure 2 .
Figure 2. Variation through time of three main climate factors and comparison with the identified abrupt transitions (TP) in the CENOGRID benthic δ 18 O.a) Global Mean Sea Level in meters from Miller et al.58 .Identification of specific warm and of glaciation events.The Laurentide, GIW-WAIS and Ice free lines are from Miller et al.58 ; b) Carbonate Compensation Depth (CCD) in meters from Pälike et al.64 .The purple circles identify the TPs on this record; c) Estimate of the CO 2 concentration in parts per million volume (ppmv) from Beerling & Royer65 .The Antarctica glaciation threshold at 750 ppmv and the NH glaciation threshold at 280 ppmv lines are from DeConto et al.66

Figure 3 .
Figure 3. Probability density of the climate system in the projected CENOGRID benthic δ 18 O and δ 13 C. space.a) Chronologically ordered TP O s (diamonds) selected according to the KS methodology for δ 18 O with time window 1-4 My are shown.The two extra TP O s found via RP are indicated with a *.b) Chronologically ordered fast TP O s (FTP O s, diamonds) selected according to the KS methodology for δ 18 O

Figure 4 .
Figure 4. RQA of U1308 benthic δ 18 O.a) Time series in ka BP; b) RP; and c) RR.Pink crosses and green lines as in Fig. S1.TP O 9 and RTP O 1-4 abrupt transitions are identified from the RR.U1308 benthic δ 18 O data are from Hodell & Channell 25

Figure 5 .
Figure 5. Evolution of the Earth Climate history among 2 different dynamical landscapes and a potential third one.The first dynamical landscape, in light red, corresponds to the Hot-Warm House time interval.The second one, in light blue, represents the Cold-Ice house time interval.The third one, in light green, highlights the potential new dynamical landscape represented by the Anthropocene time interval.The different abrupt transitions identified in the present study are reported as TPos or RTPos to distinguish the major tipping points up to the early Quaternary period from the more recent TPs characterizing less drastic climate changes during the Quaternary.Various plate tectonic and ice sheet events are indicated and supported by maps of plate movements and North and South Hemisphere ice sheets.The Antarctica maps are from Pollard & DeConto79 , Northern Hemisphere ice sheet maps are from Batchelor et al.80 .The paleogeographic maps have been generated using the Ocean Drilling Stratigraphic Network (ODSN) plate tectonic reconstruction service: <https://www.odsn.de/odsn/services/paleomap/paleomap.html>.The red arrows on the tectonic maps indicate the key events that correspond to the identified abrupt transition.

Figure 3 .
Figure 3. Probability density of the climate system in the projected CENOGRID benthic δ 18 O and δ 13 C. space.a) Chronologically ordered TP O s (diamonds) selected according to the KS methodology for δ 18 O with time window 1-4 My are shown.The two extra TP O s found via RP are indicated with a *.b) Chronologically ordered fast TP O s (FTP O s, diamonds) selected according to the KS methodology for δ 18 O with time window 0.25-1 Myr are shown.c) Same as a), but for the δ 13 C record.The approximate timing of the TPs is indicated (rounded to .01My).The 5 Ky-long portions of trajectories before and after each TP are also plotted.Figure 4. RQA of U1308 benthic δ 18 O.a) Time series in ka BP; b) RP; and c) RR.Pink crosses and green lines as in Fig. 1.TPo9-10 and RTPo1-4 abrupt transitions are identified from the RR.U1308 benthic δ 18 O data are from Hodell & Channell 25 Figure 5. Evolution of the Earth Climate history among 2 different tipping dynamical landscapes and proposal for a potential third one.The first dynamical landscape, in light red, corresponds to the Hot-Warm House time interval.The second one, in light blue, represents the Cold-Ice house time interval.The third one, in light green, highlights the potential new dynamical landscape represented by the Anthropocene time interval.The different abrupt transitions identified in the present study are reported as TPos or RTPos to differentiate the major tipping points from the critical transitions characterizing transitions of lighter significance in the climate history.Various plate tectonic and ice sheet events are indicated and supported by maps of plate movements and North and South Hemisphere ice sheets.The Antarctica maps are from Pollard & DeConto 79 , Northern Hemisphere ice sheet maps are from Batchelor et al. 80 .

Figure 4 .Channell 25 Figure 5 .
Figure 3. Probability density of the climate system in the projected CENOGRID benthic δ 18 O and δ 13 C. space.a) Chronologically ordered TP O s (diamonds) selected according to the KS methodology for δ 18 O with time window 1-4 My are shown.The two extra TP O s found via RP are indicated with a *.b) Chronologically ordered fast TP O s (FTP O s, diamonds) selected according to the KS methodology for δ 18 O with time window 0.25-1 Myr are shown.c) Same as a), but for the δ 13 C record.The approximate timing of the TPs is indicated (rounded to .01My).The 5 Ky-long portions of trajectories before and after each TP are also plotted.Figure 4. RQA of U1308 benthic δ 18 O.a) Time series in ka BP; b) RP; and c) RR.Pink crosses and green lines as in Fig. 1.TPo9-10 and RTPo1-4 abrupt transitions are identified from the RR.U1308 benthic δ 18 O data are from Hodell & Channell 25 Figure 5. Evolution of the Earth Climate history among 2 different tipping dynamical landscapes and proposal for a potential third one.The first dynamical landscape, in light red, corresponds to the Hot-Warm House time interval.The second one, in light blue, represents the Cold-Ice house time interval.The third one, in light green, highlights the potential new dynamical landscape represented by the Anthropocene time interval.The different abrupt transitions identified in the present study are reported as TPos or RTPos to differentiate the major tipping points from the critical transitions characterizing transitions of lighter significance in the climate history.Various plate tectonic and ice sheet events are indicated and supported by maps of plate movements and North and South Hemisphere ice sheets.The Antarctica maps are from Pollard & DeConto 79 , Northern Hemisphere ice sheet maps are from Batchelor et al. 80 .
reconstruction and our present analysis.This reconstruction estimates the EOT global sea level drop at 34 Ma (TP o 6) as about 30m, while previous reconstructions by Houben et al. 9 , and by Miller et al. 10 found a decrease of 70m-80m.Miller et al 10 also estimate this drop in global sea level being associated with a 2.5°C cooling interpreted previously as above the onset of the Antarctic ice sheet glaciation.Rohling et al. 8 also identify 14 Ma (TP o 7) threshold as the end of the last intermittently ice free period in the Earth history of the last 40 Ma with only the southern hemisphere ice sheets impacting the Earth climate.Miller et al. 10 indicates a 35m lowering at that particular transition.Indeed Rohling et al. 8 indicate a slight negative sea level shift at around 10 Ma (TP o 8), of about 10m, , which they interpret as the onset of partial or ephemeral northern Hemisphere ice bodies with two other sea level thresholds at about 3 Ma (close to TP o 9a) and 2.75 Ma (TP o 9b), also observed in Miller et al.

Figure S1 :
Figure S1: KS test and Recurrence Quantification Analysis (RQA) of CENOGRID benthic δ 13 C. A) KS test identifying abrupt transitions towards warmer conditions in red and cooler or colder conditions in blue; B) Recurrence plot (RP), and C) Recurrence rate (RR).The pink crosses and vertical green lines indicate the abrupt transitions (Table)detected by the RQA.

1.1 Discussion of the Detected Critical Transitions
61,60respectively.The third one at about 47 Ma (TP O 3) represents a transition towards cooler deep waters and named the Early-Middle Eocene cooling58. The l one at about 40 Ma (TP O 4) inaugurates the period of continuous cooling due to the decrease of CO 2 concentrations that eventually leads to the TP O 6 event61, which is associated with the build-up of the Antarctica ice sheets.

5 2 with
37me window 0.25-1 Myr are shown.c)Sameasa),butfortheδ13Crecord.The approximate timing of the TPs is indicated (rounded to .01My).The 5 Ky-long portions of trajectories before and after each TP are also plotted.By combining the recurrence analysis with the quasi-potential formalism we can extract further relevant information.It is natural to ask ourselves what happens if, instead, we consider a catalogue of transitions for the δ 18 O record that are detected by considering shorter time windows (0.25 Ma-1Ma) in the KS procedure.One finds -see Fig.3B-11 of such transitions (see Suppl.Mat.).Once we report such fast TP O s (FTP O s) into the empirical bivariate pdf of the climate system in the projected (δ 13 C, δ 18 O) space, we find that they correspond to finer and smaller structures of the pdf as opposed to the case of the TPs.Hence, events that are associated with faster time scales are associated with smaller jumps between secondary maxima in the pdf belong to a hierarchically lower rung than those occurring over longer time scales.This seems to support the proposal made in37.As additional step, we repeat the same analysis leading to Fig 3A by using, instead, the δ 13 C record, which features, as mentioned above, a total of 15 TP C s. Figure3Cshows clearly that the TP C s are, for the most part, saddles that had not been flagged by the TP O s.We have evidence that the same saddle is crossed more than once in a back-andforth fashion few millions of years apart (e.g.TP C s 4 & 5; TP C s 6 & 7; TP C s 12 & 13), which supports the dynamical interpretation discussed here.Comparing Figs 3A and 3C, one discovers that the same saddle is responsible for TP O 8 and for TP C 15 even if the dating is different.Similarly, the same saddle is responsible for TP O 5 and TP C s 6 & 7.
[73][74][75][76][77]tages (MIS) as already observed in numerous records covering the same interval (see Rousseau et al.72and references.therein).The first two transitions broadly match the previously discussed CENOGRID's TP O 9 associated with the onset of the Pleistocene (note that the interval between the two is much smaller than the resolution needed to separate to TP O s) and are followed by four more recent TPs associated with the benthic δ18O record (RTP O s).There is clear evidence of the Mid-Pleistocene (MPT) critical transition, between 1.25 Ma and 0.8 Ma, during which a shift occurred from climate cycles dominated by a 40-kyr periodicity (due to obliquity) to 100-Kyr periodicity (due to eccentricity) dominated ones[73][74][75][76][77].The 1.25 Ma date is particularly significant, since it is followed by an Increase in the amplitude of glacial-interglacial fluctuations.A complementary RQA of the δ 18 O bulk carbonate record from U1308, which characterizes episodes of iceberg calving into the North Atlantic Ocean IRD released into the North Atlantic Ocean The KS augmented test and the RQ of the benthic δ 18 O agree in identifying six abrupt transitions dated at approximately 2.93 Ma, at 2.52 Ma, at 1.51 Ma, at 1.25 Ma, at 0.61 Ma and at 0.35 Ma (Fig.4, Suppl.Tab. 1).They characterize the dynamics of North Hemisphere ice sheets (elevation and spatial expansion) and agree with the classical transitions characterizing . From TP o 6, at 34 Ma, and 23 Ma, the CO2 concentration decreases associated with a deepening trend in the CCD down to about -4600m.From 23 Ma until about 19 Ma, GMSL shows oscillations but with lower values than present day at about -20m, whereas from 19 Ma until 17 Ma, GMSL increases by about 50m to indicate high values around +30 m above present day value (Fig.

Table S1 :
Recurrence Quantification Analyses (RQA) of CENOGRID benthic δ 18 O and δ 13 C, and U1308 benthic δ 18 O, δ 13 C, and bulk carbonate δ 18 O.For each, listed are the dates of the abrupt transitions, the corresponding RR prominence values used to select the abrupt transitions, and the significance of those transitions.X denotes prominence values above the RR standard deviation, while $ denotes values below the standard deviation.Ages are in Ma for the CENOGRID and in ka BP for the U1308 records respectively.

Table S2 :
Statistics of the 66-34 Ma and 34 Ma-present intervals, from top to bottom: the Global Mean sea level (GMSL) in meters from Miller et al. 10 , the CO2 concentration in ppmv from Beerling and Royer 2 , and for the CCD depth in meters from Palike et al. 1 .

Table S3 :
Summary of the GMSL, CCD and CO2 concentration trends at the identified abrupt transitions TP1 to TP10.