Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations

. The relative role of external forcing and of intrinsic variability is a key question of climate variability in general and of our planet’s paleoclimatic past in particular. Over the last 100 years since Milankovic’s contributions, the importance of orbital forcing has been established for the period covering the last 2.6 Myr and the Quaternary glaciation cycles that took place during that time. A convincing case has also been made for the role of several internal mechanisms that are active on timescales both shorter and longer than the orbital ones. Such mechanisms clearly have a causal role in Dansgaard–Oeschger and Heinrich events, as well as in the mid-Pleistocene transition. We introduce herein a uniﬁed framework for the understanding of the orbital forcing’s effects on the climate system’s internal variability on timescales from thousands to millions of years. This framework relies on the fairly recent theory of non-autonomous and random dynamical systems, and it has so far been suc-cessfully applied in the climate sciences for problems like the El Niño–Southern Oscillation, the oceans’ wind-driven circulation, and other problems on interannual to interdecadal timescales. Finally, we provide further examples of climate applications and present preliminary results of interest for the Quaternary glaciation cycles in general and the mid-Pleistocene transition in particular.


Introduction and motivation
In the early 20th century, Milutin Milankovic presented his theory of ice ages (Milankovitch, 1920). Based on his own calculations and on insightful suggestions from Wladimir Köppen and Alfred Wegener (Imbrie and Imbrie, 1986), he proposed that the transitions between glacial and interglacial climate conditions were primarily caused by variations of incoming solar radiation, which by that time was known to vary in a quasi-periodic manner on slow timescales of tens of thousands to hundreds of thousands of years (Poincaré, 1892(Poincaré, -1899. These variations of insolation, which arise as a consequence of the gravitational interaction of the Earth with the other planets and with its own Moon, are typically referred to as orbital forcing. The orbital forcing comprises variations in (i) the eccentricity of the Earth's orbit around the sun with dominant spectral peaks around 400 and 100 kyr; (ii) the obliquity, or axial tilt, i.e., the angle between the Earth's rotational and its orbital axis, with dominant periodicity around 41 kyr; and (iii) the climatic precession, which determines the phase of the summer solstice along the Earth's orbit and has its most pronounced spectral power around 23 and 19 kyr (Berger, 1978).
For 2 centuries or so of modern geology, records of our planet's physical and biological past were merely discrete sequences of strata with specific properties, like coloration 864 K. Riechers et al.: Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations and composition (Imbrie and Imbrie, 1986). This state of affairs led, after the initial success of the Milankovitch (1920) theory of the ice ages, to severe criticism of the temporal mismatch between insolation minima and glaciation maxima (e.g., Flint, 1971).
The advent of marine-sediment cores after World War II led, for the first time, to the availability of records that were more or less continuous in time. Like all climate records, these cores covered limited time intervals and did so with limited resolution and with inaccuracies in absolute dating, as well as in the quantities being measured. Moreover, they posed the problem of inverting proxy records of isotopic and microbiotic counts to physical quantities like temperature and precipitation.
In spite of these limitations, the spectral analysis of deepsea records allowed Hays et al. (1976) to overcome the difficulties previously encountered by the orbital theory of Quaternary glaciations, in particular the absence of the imprint of precessional and obliquity peaks in glaciation proxy records. Specifically, Hays et al. (1976) were able to create a composite record -back to over 400 kyr b2k, i.e., over 400 000 years before the year 2000 CE -from two relatively long marinesediment records of the best quality available in the early 1970s. The authors demonstrated therewith that significant precessional and obliquity peaks near 20 and 40 kyr were present in this record's spectral analysis; see Fig. 1. The power spectrum in the figure also made it quite clear that these peaks were superimposed on a continuous background -the stippled area in the figure -whose total variance much exceeded the sum of the variances present in the peaks.
The work of Hays et al. (1976) and of the subsequent CLIMAP and SPECMAP projects resulted in a much more detailed spatiotemporal mapping of the Quaternary and extended the belief in the pacemaking role of orbital variations into the more remote past. The spectral peaks near 20 and 40 kyr have been widely interpreted within the geological community as evidence for a linear response of the climate system to the orbital forcing (Imbrie and Imbrie, 1986). A third spectral peak at 100 kyr was, however, the most pronounced but much more difficult to reconcile with the orbital theory of Quaternary glaciations. Since no sufficiently pronounced counterpart can be found in the spectra of the seasonal insolation forcing, Hays et al. (1976) hypothesized a nonlinear response of the climate system in order to explain this dominant periodicity of the Late Pleistocene glacialinterglacial cycles. At the same time, the advent of higherresolution marine cores and especially ice cores from both Greenland and Antarctica, led to the discovery of Heinrich events (Heinrich, 1988), Dansgaard-Oeschger (D-O) events (e.g., Dansgaard et al., 1993), and Bond cycles (Bond et al., 1997), which were hard to explain by orbital forcing given their shorter timescales.
In fact, interest in past climates was heightened not only by these striking observational discoveries but also by the growing concerns about humanity's impact on the climate (SMIC, Figure 1. Power spectrum of a composite δ 18 O record using deepsea cores RC11-120 and E49-18. This figure is based on the work of Hays et al. (1976), as presented by Imbrie and Imbrie (1986). Reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer Nature. Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics by Ghil and Childress (1987), © 1987 by Springer Sci-ence+Business Media New York. All rights reserved. 1987All rights reserved. . 1971National Research Council, 1975). Given the declining temperatures between the 1940s and 1970s, on the one hand, as shown in Fig. 2 (see also Ghil and Vautard, 1991, Fig. 3), and the substantial advances in the description of the Quaternary glaciations, on the other hand, it is clear that interest was mainly on the planet's falling into another ice age (e.g., National Research Council, 1975).
As a result of the twofold stimulation provided by data about past glaciations and concern about future ones, a number of researchers in the early-to-mid 1970s worked on energy balance models (EBMs) of climate with multiple stable steady states (Held and Suarez, 1974;North, 1975;Ghil, 1976). Two such stable "equilibria" corresponded to the present climate and to a "deep-freeze", as it was called at the time, i.e., to a totally ice-covered Earth. At the time there was some disbelief about this second climate, as its calculated temperatures were much lower than those associated with the Quaternary glaciations and incompatible with paleoclimatic evidence available in the 1970s.
New geochemical evidence, however, led in the early 1990s to the discovery of a "snowball" or, at least, "slushball" Earth prior to the emergence of multicellular life, sometime before 650 Myr b2k (Hoffman et al., 1998). It thus turned out that this climate state -predicted by several EBMs and confirmed by a general circulation model (GCM) with much higher spatial resolution (Wetherald and Manabe, 1975) -had actually occurred and is now being modeled in much greater detail (Pierrehumbert, 2004;Ghil and Lucarini, 2020).
On the other hand, it also became clear that these early models, whose only stable solutions were stationary, could not reproduce the wealth of variability that the proxy records were describing very well, not even in the presence of stochastic forcing (e.g., Ghil, 1994). Certain theoretical paleoclimatologists therefore turned to coupling a "climate" equation, with temperature as its only dependent variable, with an ice sheet equation (Källén et al., 1979;Ghil and Le Treut, 1981) or a carbon dioxide equation (Saltzman et al., 1981;Saltzman and Maasch, 1988). These coupled climate models, albeit highly idealized, did produce oscillatory solutions that captured some of the features of the Quaternary glaciation cycles as known at that time. For instance, the models of Ghil and associates (Källén et al., 1979;Ghil and Le Treut, 1981) captured the phase differences between peak ice sheet extent and minimum temperatures in the North Atlantic suggested by Ruddiman and McIntyre (1981), while the work of Saltzman and associates (e.g., Saltzman and Maasch, 1988) captured the asymmetry of the glaciation cycles with their more rapid "terminations" (Broecker and Van Donk, 1970).
The stable self-sustained oscillations of these coupled models, however, were totally independent of any orbital or other time-dependent forcing, i.e., the solar input to their radiative budget was constant in time. Hence, they could not capture the wealth of spectral features, with their orbital and other peaks, of the paleo-records available by the 1980s. The basic quandary of the Quaternary glaciation cycles -at least from the point of view of theoretical climate dynamics (Ghil and Childress, 1987, Part IV) -is formulated in Fig. 3 below; see also Ghil (1994): How does the quasi-periodic orbital forcing, with its relatively narrow spectral peaks, act on the climate's internal variability to produce a response that is characterized by significant spectral peaks superimposed on a broad continuous background? In addition, how does the climate's spectral peak at 100 kyr arise given its absence in the power spectrum of the forcing?
In this paper, we try to show a path toward resolving the four fundamental questions listed in Box 1. In the next section, we summarize existing results on how the climate system's intrinsic variability arises on Quaternary timescales and on how this variability is modified by the time-dependent orbital forcing, which was added to the previously autonomous climate models as the next step in paleoclimate modeling evolution; see, for instance, Le Treut and Ghil (1983) and Le Treut et al. (1988) vs. Ghil and Le Treut (1981). In Sect. 3, we outline a more general framework for the study of such mechanisms, as given by the theory of non-autonomous and random dynamical systems (NDSs and RDSs), and sketch an application of this theory to other climate problems. An application to the problem at hand is proposed in Sect. 4 and conclusions follow in Sect. 5.

866
K. Riechers et al.: Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations Figure 3. The theoretical quandary of modeling the Quaternary glaciation cycles, illustrated here by schematic diagrams of the composite power spectral densities (PSD) of (a) the paleo-records and (b) the orbital forcing. In panel (a), the dominant peak for the Late Pleistocene is near 100 kyr, while in panel (b) eccentricity forcing is distributed over several spectral lines. The peaks at 6-8 and 1-2 kyr in panel (a) correspond to Bond cycles (Bond et al., 1992 and to the mean recurrence of D-O events, and they lack a match in the forcing lines of panel (b). Box 1. Fundamental questions regarding the Quaternary glacialinterglacial cycles.

A simple mechanism for climate oscillations
We follow Ghil (1994) in sketching the simplest physical mechanism for a self-sustained climate oscillation at fixed insolation forcing. Consider the Källén et al. (1979, hereafter KCG) oscillator, to the best of our knowledge the first such self-sustained climate oscillator. The model itself was built from the ground up, coupling a scalar version of the Ghil (1976) energy balance model (EBM) with a simplified, scalar version of the Weertman (1964Weertman ( , 1976 ice sheet model (ISM). The model's details and further analyses of its ingredients and variants can be found in several references (e.g., Crafoord and Källén, 1978;Ghil and Tavantzis, 1983;Ghil, 1984;Bódai et al., 2015).
The basic workings of this climate oscillator can be represented by two coupled ordinary differential equations (ODEs), written symbolically as follows: Here T stands for global temperature and V for global ice volume, while Eq. (1a) is an EBM and Eq. (1b) is an ice sheet model (ISM). The " " symbol stands for a binary relation of rough proportionality and is intended to neglect the details of the equation's right-hand side (RHS), including its nonlinearities. The EBM represents the well-known ice-albedo feedback used by both Budyko (1969) and Sellers (1969), while the ISM relies on the precipitation-temperature feedback postulated by KCG and used also by Ghil and Le Treut (1981), who coined the term. The latter feedback can be better understood by writing the following equations: Here p is net precipitation on the single ice sheet of the globally integrated model, given by the difference in Eq. (2b) between the accumulation p ac and the ablation p ab (KCG). As first observed by George C. Simpson -the meteorologist of Robert F. Scott's Terra Nova expedition to Antarctica in 1910Antarctica in -1912 and later the longest-serving Director of the U.K. Meteorological Office -warmer winters have more snow, and hence, at least in central Antarctica, the increase of p ac with T exceeds the more obvious increase of p ab with T . Hence, p T , and we have derived therewith Eq. (1b), i.e., V p T . For more recent studies of the precipitationtemperature feedback, see Tziperman and Gildor (2002).
More generally, the presence of feedbacks of opposite sign in a system of two linear coupled ODEs, leads to an oscillation, with the solution given by two trigonometric functions in quadrature with each other, x(t) = sin(t), y(t) = cos(t) = sin(t + π/2), and the trajectory describing a circle in the (x, y) phase plane, x 2 (t) + y 2 (t) = 1. In a nonlinear system, however -like the full KCG model or any other climate oscillator mentioned so far -the possibility of an oscillation, as indicated by the system (1), is actually realized in the explicit, full set of equations only for certain parameter values and not for others. This can be understood by considering the so-called normal form of a Hopf bifurcation, which leads from a stable steady state, called a fixed point in dynamical systems theory, to a stable oscillatory solution, called a limit cycle. The easiest way to see this transition is by writing the normal form in polar coordinates, as in Arnold (1988) and in Ghil and Childress (1987, Sect. 12.2), namely, Here z = x + ıy is complex, where ı = √ −1 is the imaginary unit, while µ is a real bifurcation parameter, and c, ω are real and nonzero. Note that the KCG model per se is not in the normal form above, and we will discuss its bifurcation parameter µ * in the next subsection.
A very natural transformation of variables, leads from the complex ODE (3) to the system of two real and decoupled ODEs, Equation (5b) simply provides an angular rotation around the origin ρ = 0 = x = y, since the complex exponential in Eq. (4) is periodic with period 2π . Equation (5a) is quadratic in ρ, and thus it can have two real roots, ρ = ρ 0 = 0 and ρ = ρ * = −µ/c. But ρ has to be positive, and thus in the case in which c < 0, the only possible solution for µ < 0 is the origin, and it is stable since ρ(µ + cρ) is negative for ρ > 0 in this case; hence, ρ has to be monotonically decreasing, i.e., all the solutions of Eq. (5) spiral into the origin. The Hopf bifurcation from this stable steady state to a periodic solution, i.e., a limit cycle with radius ρ 1/2 = ρ 1/2 * , occurs as µ crosses 0. Since now ρ(µ+cρ) > 0 for 0 < ρ < ρ * and ρ(µ+cρ) < 0 for ρ > ρ * , the limit cycle is stable and trajectories spiral out from inside this cycle and into it from outside; see Fig. 4 for the so-called supercritical (or soft) Hopf bifurcation case with c < 0.

Intrinsic climate oscillations and the mid-Pleistocene transition (MPT)
In this subsection, we present an argument for the role of intrinsic oscillations in the mid-Pleistocene transition (MPT). The first point to be made is that while orbital forcing clearly plays a major role in the power spectrum of the Quaternary's climatic variability, it cannot be, in and of itself, the cause of the MPT. Indeed, changes in the solar system's orbital periodicities only occur on much longer timescales than the entire Quaternary's duration (Varadi et al., 2003;Laskar et al., 2004a). Our argument continues with a further analysis of the Hopf bifurcation presented in the previous subsection. Such an analysis was carried out for the KCG model by Ghil and Tavantzis (1983). Physically speaking, the presence or absence of the regular, purely periodic oscillations obtained by KCG and illustrated in Ghil and Childress (1987, Fig. 12.6) depends on whether c≷0 in Eq. (5a). The KCG model's bifurcation parameter is µ * = c T /c L , where c T is the heat capacity in its EBM, while c L is the "mass capacity" in its ISM (Ghil and Tavantzis, 1983). Large µ * corresponds physically to a very small, possibly pre-Pleistocene ice cap (Ghil, 1984;Saltzman and Sutera, 1987). At these values of µ * , the KCG model's nullclines and fixed points -the latter being given by the intersection of the former -are very different from those that are obtained for Quaternary-sized ice sheets, for which c L is comparable in value to c T ; see Ghil and Tavantzis (1983, Figs. 3-5). As µ * decreases to O(1), i.e., as we proceed from very small to more substantial ice sheets, the fixed point transfers its stability to a branch of periodic solutions via a subcritical Hopf bifurcation (Källén et al., 1979;Ghil and Tavantzis, 1983); see also Ghil and Childress (1987, Figs. 12.8 and 12.9).
To clarify the simple physical concepts that underlie subcritical and supercritical Hop bifurcations, let us consider a purely mechanical oscillator with mass m, a spring kx, and a dashpot αẋ (Landau and Lifshitz, 1960;Jordan and Smith, 1987), If k = const. and α = const., we have the simplest, linear setup, but we will be interested in the nonlinear cases here. Normalizing by the mass m and not changing notation otherwise, by rearranging terms and adding a periodic forcing we get the following equation: x + α(x)ẋ + k(x)x = F cos(ωt).
Two classical nonlinear cases are those of the Duffing (1918) equation, in which k(x) = x 2 and α = const., and of the Van der Pol (1926) equation, in which k = const. and α(x) = ν(x 2 − 1). The fully nonlinear case in which both the spring and the damping are nonlinear, with k(x) = x 2 and α(x) = ν(x 2 − 1), is known as the Van der Pol-Duffing oscillator In this case, the origin constitutes the only stable fixed point, and all trajectories will spiral into this point, as illustrated by the single brown trajectory.
(b) Vector fieldż(z) and three trajectories for µ = 1. In this case, the origin is an unstable fixed point, while the limit cycle with radius ρ 1/2 * = 1 constitutes the only stable solution. Trajectories that start inside this limit cycle, with ρ(t 0 ) < ρ * , tend to it by spiraling out -as illustrated by the magenta trajectory -while trajectories that start outside the limit cycle, with ρ(t 0 ) > ρ * , approach it by spiraling inward, as illustrated by the gray trajectory. (c) Dependence of the stable solution of Eq. (3) on the parameter µ, for c = −1. For µ ≤ 0, the single stable solution is the fixed point located at the origin, ρ ≡ 0. For µ > 0, the stable solution is the limit cycle given by ρ = µ.
(e.g., Jackson, 1991;Pierini et al., 2018). Note that all three types of nonlinear oscillators can exhibit chaotic behavior even in the presence of simple periodic forcing (e.g., Guckenheimer and Holmes, 1983;Pierini et al., 2018, and references therein). The idea of using such simple, classical oscillators in modeling Quaternary glaciation cycles goes back to Saltzman et al. (1981). Jordan and Smith (1987, Sect. 5.6) specifically discuss the case of a soft and a hard spring for a generalized Duffing equation, with k(x) = k 0 + h(x), where 0 < | | 1, h(−x) = h(x), and h > 0. A spring is soft if it is sublinear, < 0, and hard if it is superlinear, > 0; see their Eq. (5.37) and Fig. 5.4a, b, with h(x) = x 2 .
The supercritical Hopf bifurcation in the absence of forcing is analogous to the nonlinear response of a soft, sublinear spring to periodic forcing in which the oscillations in the position x of the mass m increase gradually in amplitude as the spring constant k 0 increases past a critical value k * , while the subcritical case is analogous to the response of a hard, superlinear spring, for which the oscillations in x jump suddenly from zero amplitude to a finite amplitude as the spring constant k 0 crosses the value k * . Please compare the behavior of supercritical and subcritical Hopf bifurcations in Ghil and Childress (1987, Sect. 12.2 and Figs. 12.7-12.9) and see Jordan and Smith (1987, Fig. 5.7) for the change in the nonlinear response of a Duffing oscillator as its spring changes from soft, with < 0, to hard, with > 0.
There is a clear-cut analogy with the mid-Pleistocene transition, occurring at roughly 0.8 Ma b2k, at which smallamplitude climate variability with a dominant periodicity near 40 kyr becomes larger, dominated by a periodicity that is close to 100 kyr, as well as being more irregular (e.g., Huybers, 2009;Quinn et al., 2018;Rousseau et al., 2020). A fair number of distinct dynamical theories for this transition have been formulated (e.g., Maasch and Ghil, 1994;Ashwin and Ditlevsen, 2015;Daruka and Ditlevsen, 2016;Omta et al., 2016;Ditlevsen and Ashwin, 2018). A rather obvious one is that a Hopf bifurcation occurs at that point, which leads to a more vigorous response to the multi-periodic orbital forcing; thus, the latter does not need to change in the least in order to explain the observed phenomena. In Saltzman and Sutera (1987), there is only a comment on the likely role of a Hopf bifurcation in the transition, but their Fig. 3 suggests that in their model such a bifurcation would have to be of the supercritical type and lead to a fairly gradual transition.
In contrast, the subcritical Hopf bifurcation of the KCG and Ghil and Le Treut (1981) oscillators would have to lead to a more abrupt transition, as suggested by Ghil (1984). Later, Crucifix (2012) showed that the models by Maasch (1990, 1991) exhibit MPT-like behavior via supercritical or subcritical Hopf bifurcations, depending on the parameter values. The existing δ 18 O and δ 13 C records might or might not have sufficient resolution back in time up to 1.2 Ma to settle this question about the abruptness of the transition. In case some of them do, an objective test of suddenness -as proposed by Bagniewski et al. (2021) for the highresolution NGRIP record (North Greenland Ice Core Project members, 2004) -will have to be applied to such records.

Orbital forcing of a climate oscillator
We start this section by describing some fairly simple ways in which the orbital forcing might have modified intrinsic climate variability, thus helping to solve the mismatch between Fig. 3a and b in Sect. 1. To explore this possibility, Le Treut and Ghil (1983) used somewhat simplified insolation forcing based on the calculations of Berger (1978) and applied it to a slightly modified version of the Ghil and Le Treut (1981) oscillator. These authors found that, as expected for a nonlinear oscillator, its internal frequency f 0 is affected strongly by the forcing ones, {f 1 , . . ., f 5 }, resulting in both nonlinear resonance and combination tones (Landau and Lifshitz, 1960).
Depending on parameter values, the periodicity P 0 = 1/f 0 of the Ghil and Le Treut (1981) oscillator is P 0 6-7 kyr. The lines in the simplified insolation spectrum used by Le Treut and Ghil (1983) had the following periodicities: P 1 = 19 kyr, P 2 = 23 kyr, P 3 = 41 kyr, P 4 = 100 kyr, and P 5 = 400 kyr. These periodicities correspond to the two precessional ones, the obliquity one, and the two eccentricity ones. The actual celestial-mechanics calculations that Berger (1978) based his insolation calculations on have been substantially updated since (e.g., Varadi et al., 2003;Fienga et al., 2015). However, these advances have not modified the spectrum of the planetary-orbit solutions over the 2.6 Myr of the Quaternary very much, which is a rather short interval in celestial-mechanics terms.
The results of the Le Treut and Ghil (1983) model on the evolution of the primary climate variables T and V were converted to δ 18 O values in simulated isotopic records of marine sediment and ice cores by Le Treut et al. (1988); the spectra of the latter are plotted in Fig. 5. The values on the abscissa of Fig. 5a are values of the logarithm of frequency, while those in Fig. 5b are values of the frequency itself; the values on the ordinate of both panels are powers of 10. One refers to such figures as being in (a) log-log coordinates vs. (b) log-linear coordinates for short.
Aside from the spectral features noted in the figure caption and discussed in greater detail by Ghil (1994), it is important to realize (i) that the large continuous background in Fig. 5a is purely of deterministically chaotic origin since there is no stochastic element whatsoever in the Le Treut and Ghil (1983) model or in its forcing and (ii) that the dominant peak at 109 kyr is not directly forced by the f 4 = 1/100 kyr −1 eccentricity line but rather it is due to the difference tone between the two precessional frequencies, f 1 and f 2 . Finally, it is the nonlinear, broad resonance of the model's f 0 frequency with the quasi-periodic forcing that produces the bump in the spectrum of Fig. 5a to the right of the orbital frequencies.
In returning to the "fundamental question no. 2" in Box 1, one must recall that on the paleoclimatic timescales of interest -apart from deterministic chaos à la Lorenz (1963), as obtained by Le Treut and Ghil (1983) and Le Treut et al. (1988) and shown here in Fig. 5a -stochastic contributions à la Hasselmann (1976) to the continuous part of the spectrum must also play an important role. In fact, RDS theory, as outlined in the next subsection, provides an excellent framework for a "grand unification" of these two complementary points of view (Ghil, 2014(Ghil, , 2019. In the paleoclimatic context, Ditlevsen et al. (2020) have suggested that, aside from red noise processes, dating uncertainties in the proxy records from which the spectra are derived may contribute, in all likelihood, to this background; see also Boers et al. (2017a, b). In this context, Verbitsky and Crucifix (2020) also provide a simple theory that addresses scaling properties in the glacial cycles and their spectra, based on the so-called Buckingham π -theorem (e.g., Barenblatt, 1996).

Basic facts of NDS and RDS life
The highly preliminary results summarized in Sect. 3.1 encourage us to pursue the effects of the orbital forcing on intrinsic climatic variability in a more systematic way, effects that may have contributed to generate the rich paleoclimate spectrum on Quaternary and longer timescales (e.g., Westerhold et al., 2020). In fact, several research groups in the climate sciences have carried out during the past 2 decades an important extension of the dynamical systems and model hierarchy framework presented by Ghil and Childress (1987) and by Ghil (2001), from deterministically autonomous to non-autonomous and random dynamical systems (NDSs and RDSs, e.g., Ghil et al., 2008;Bódai and Tél, 2012).
On the road to including deterministically time-dependent and random effects, one needs to realize first that the climate system -as well as any of its subsystems on any timescaleis not closed: it exchanges energy, mass, and momentum with its surroundings, whether this pertains to other subsystems or the interplanetary space and the solid earth. The typical applications of dynamical systems theory to climate variability until not so long ago have only taken into account exchanges that are constant in time, thus keeping the model -whether governed by ordinary, partial, or other differential equations -autonomous; i.e., the models had coefficients and forcings that were constant in time.
Alternatively, the external forcing or the parameters were assumed to change either much more slowly than a model's internal variability, meaning that the changes could be assumed to be quasi-adiabatic, or much faster, meaning that they could be approximated by stochastic processes. Some of these issues are covered in much greater detail by Ghil and Lucarini (2020, Sect. III.G). The key concepts and tools of NDSs and RDSs go beyond such approaches that rely in an essential way on a scale separation between the characteristic times of the forcing and the internal variability of a given system; such a separation is rarely if ever actually present in the climate sciences.
The presentation of the key NDS and RDS concepts and tools in this subsection is aimed at as large a readership as possible and follows Ghil (2014). Slightly more in-depth but still fairly expository presentations can be found in Crauel and Kloeden (2015) and in Caraballo and Han (2017). Readers who are less interested in this mathematical frameworkwhich allows a truly thorough understanding of the way that orbital forcing acts on intrinsic climate variability on Quaternary timescales -may skip at a first reading the remainder of this section and continue with Sect. 4.  Treut et al. (1988) with idealized orbital forcing. Panel (a) is in log-log coordinates and clearly shows a dominant peak at 109 kyr and a very large amount of variance in the continuous spectrum, which has a roughly −2 power law slope. The orbital peaks at P 1 = 19 kyr, P 2 = 23 kyr and P 3 = 41 kyr are also present, along with peaks at 14.7 and 10.4 kyr, which correspond to the sum tones f 2 +f 3 and f 1 +f 2 , respectively. Panel (b) is in log-linear coordinates and shows additional harmonics and sum tones, as well as the difference tone f 1 − f 2 , which corresponds to the dominant peak at 109 kyr. Inside the figure, the notation ka instead of kyr has been kept unchanged from the original publication.

Autonomous and non-autonomous systems
Succinctly, one can write an autonomous system aṡ where X stands for any state vector or climate field. While F is a smooth function of X and of the parameter µ, it does not depend explicitly on time. This autonomous character of Eq. (8) greatly facilitates the analysis of its solutions' properties. For instance, two distinct trajectories, X 1 (t) and X 2 (t), of a well-behaved, smooth autonomous system cannot intersect -i.e., they cannot pass through the same point in phase space -because of the uniqueness of solutions. This property helps one draw the phase portrait of an autonomous system, as does the fact that we only need to consider the behavior of solutions X(t) as time t tends to +∞. The sets of points so obtained are (possibly multiple) equilibria, periodic and quasiperiodic solutions, and chaotic sets. In the language of dynamical systems theory, these are called, respectively: fixed points, limit cycles, tori, and strange attractors.
We know only too well, however, that the seasonal cycle plays a key role in climate variability on interannual timescales, while orbital forcing is crucial on the Quaternary timescales of many millennia. In addition, more recently it has become obvious that anthropogenic forcing is of utmost importance on the interdecadal timescales between these two extremes.
How can one take into account these types of timedependent forcings and analyze the non-autonomous systems that they lead us to formulate? One succinctly writes such a system as follows: In Eq. (9), the dependence of F on t may be periodic, F (X, t + P ) = F (X, t), as in various El Niño-Southern Oscillation (ENSO) models, where the period P = 12 months, or monotonic, |F |(X, t + τ ) ≥ |F |(X, t), as it is when studying scenarios of anthropogenic climate forcing. An even more general situation includes time dependence in one or more parameters {µ 1 , . . ., µ p }.
To illustrate the fundamental distinction between an autonomous system like Eq. (8) and a non-autonomous one like Eq. (9), consider the simple scalar version of these two equations: anḋ respectively. We assume that both systems are dissipative, i.e., β > 0, and that the forcing is monotonically increasing, γ ≥ 0, as would be the case for anthropogenic forcing in the industrial era. Lorenz (1963) pointed out the key role of dissipativity in giving rise to strange but attracting solution behavior, while Ghil and Childress (1987) emphasized its importance and pervasive character in climate dynamics. Clearly the only attractor for the solutions of Eq. (10a), given any initial point X(0) = X 0 , is the stable fixed point X * = 0, attained as t → +∞.
In the case of Eq. (10b), however, this forward-in-time approach yields blow-up from t → +∞ for any initial point. To make sense of what happens in the case of time-dependent forcing, one instead introduces the pullback approach, in which solutions are allowed to still depend on the time t at which we observe them but also on a time s from which the solution is started, with X(s) = X 0 and s t. With this little change of approach, one can easily verify that for all t and X 0 , where A t = (γ /β)(t−1/β). We obtain therewith, in this pullback sense, the intuitively obvious result that the solutions, if we start them far enough in the past, all approach the family of attracting sets A t ; this family follows the forcing γ t, and it thus has a linear growth in time t. Hence, the fixed point X * of Eq. (10b) is in fact a moving target and it is given by X * = γ t/β. Due to the system's inertia, the set A t that is approached by the trajectories lags this time-dependent fixed point by a constant offset of γ /β 2 .

Pullback attractor (PBA)
Formally, the indexed family of all pullback attracting sets A t , is termed the pullback attractor (PBA) of an NDS if the following two conditions are fulfilled: i. each snapshot A t is compact and the family of snapshots {A(t)} t∈R is invariant with respect to the dynamics, and ii. the pullback attraction occurs for all times, To further improve the reader's intuition for PBAs, we provide a second illustrative example here. A system defined in polar coordinates bẏ with ρ, µ ∈ R + , and φ ∈ R/2π can easily be seen to exhibit a limit cycle in the (x, y) plane with (x = ρ cos φ, y = ρ sin φ). An initial deviation of ρ from µ will decay exponentially, and the system converges to an oscillation of radius µ with the angular velocity ω. Here, we transform this autonomous dynamical system into a non-autonomous one by modulating the target radius µ with a sinusoidal forcing, where the modulation is moderate to guarantee that µ 0 + β sin(νt) > 0 for all t.
Since the dynamics of the phase φ and of the radius ρ are decoupled, the corresponding equations can be solved and analyzed separately. While the temporal development of the phase is trivial, the pullback-invariant attracting set of the radius for the initial condition ρ(s) = ρ 0 is given by as shown in Appendix B. Note that in the limit s → −∞, the dependence on the initial value ρ 0 vanishes and the attracting set A (ρ) t performs an oscillation of the same frequency as the forcing. It lags the phase of the time-dependent fixed point µ(t) by the constant ϑ, while its amplitude is amplified by the factor α. Since ρ is restricted to positive values, this solution requires αβ < µ 0 .
The PBA with respect to the coordinate ρ is comprised of the family of all the sets A (ρ) t as defined in Eq. (18) and thus reads Since the pullback limit for the phase φ does not exist, no constraints on it other than φ ∈ [0, 2π) are imposed by the dynamics. Hence, for the system (15) comprised of radius and phase, we find that the distance of any trajectory at time t -i.e. (ρ(t; t 0 , ρ 0 ), φ(t; t 0 , φ 0 )) -to the set A t = {(αβ sin(νt + ϑ) + µ 0 , ϕ) : ϕ ∈ [0, 2π )} tends to zero as we pullback the initial time t 0 to −∞. The pullback attracting sets A t at time t are circles in the (x, y) plane with a radius that oscillates in time, and the system's PBA is given by the family of these circles Figure 6 shows trajectories of the system starting from different points in the past. In Fig. 6a, the trajectories are depicted in the three-dimensional (3-D) space spanned by the two Cartesian coordinates (x, y) and the time t, where the usual transformation from polar to Cartesian coordinates was applied. The shaded surface in this panel represents the PBA of the system. Corresponding trajectories of ρ(t) and their convergence to the PBA A   (16); in all the panels the parameter values are µ 0 = 15, α = 10, ν = 0.2, and ω = 0.3, unless stated otherwise. (a) Trajectories (ρ(t), φ(t)) of the system starting from different times in the past in the 3-D space spanned by (x, y) and time t; the system's PBA lies on the red-shaded surface. (b) Temporal evolution of the radius (solid colors) of the three trajectories shown in (a) together with its PBA (dotted red). (c-e) Trajectories integrated from t 0 = −100 to t f = 200 for ν = 0.3, 0.2, 0.5 1/2 in panels (c), (d), and (e), respectively. The values are chosen such that the ratio between ν and ω is ω/ν = 1, ω/ν ∈ Z\{1} and ω/ν ∈ R\Z in panels (c) to (e). Correspondingly, in (c) the trajectory quickly converges to a circle, whose center is slightly shifted from the origin. In panel (d) a quasi-closed three-fold loop can be observed, since ω/ν = 3/2. In panel (e) the trajectories will densely fill the annular disk defined in the video supplement statement. (f) Heat map of numerous trajectories projected onto the (x, y) plane. The trajectories start at random points in state space and are integrated from t i = −200 to t f = 200. A video of the heat map filling up as more and more trajectories with different initial conditions are added is provided in the video supplement. The heat map shown here is a snapshot of the video taken at time t = 0.20; see details in the video supplement. a measure in NDSs and RDSs, please see Caraballo and Han (2017), , Crauel and Kloeden (2015), or Ghil et al. (2008). Essentially, the heat map here counts the number of times that 100 trajectories integrated from t = −200 to t = 200, with randomized initial conditions like the ones shown in Fig. 6a, cross small pixels in the (x, y) plane. Figure 6c-e demonstrate a particularity of this system, which is characteristic of dynamics confined to a torus. Namely, the structure of the system's trajectories depends on the frequency ratio ω/ν, and three different cases must be distinguished. If the radius is modulated with the same frequency as the oscillation itself, i.e. ω = ν as in panel (c), the forcing and the system have a fixed phase relation. That is, for a given phase of the system, its radius is always attracted by the same fixed point. Hence, the system practically repeats its orbit after a short time. More precisely, the radius of the oscillation does differ from one "round trip" around the torus to the next, but this difference tends to zero as ρ(t) approaches the PBA A (ρ) t .
If ω and ν are rationally related, i.e., m ω = n ν with n, m ∈ N, as in Fig. 6d, then -after n periods of the radial modulation and m periods of the system's oscillation -the phase relation between the system and its forcing will repeat itself, and hence we observe the same quasi-repetition of the orbit after the time n 2π/ν = m2π/ω. That is, such a trajectory will appear as an n-fold quasi-closed loop.
Finally, if ω/ν ∈ Z, as in Fig. 6e, then a given phase of the system will never coincide with the same phase of the radius modulation more than once. Hence, the trajectory does not repeat itself but instead densely covers the annular disk

Random attractor
Let us return now to the more general, nonlinear case of Eq. (9) and add not only deterministic time dependence F (X, t) but also random forcing, where η = η(t, ω) represents a Wiener process -with dη commonly referred to as "white noise" -and ω now labels the particular realization of this random process. When G = const., the noise is additive, while for ∂G/∂X = 0 we speak of multiplicative noise. The distinction between dt and dη in the stochastic differential Eq. (20) is necessary since, roughly speaking and following the Einstein (1905) paper on Brownian motion, it is the variance of a Wiener process that is proportional to time and thus dη ∝ (dt) 1/2 . In Eq. (20), we dropped the dependence on a parameter µ for the sake of simplicity. The noise processes may include "weather" and volcanic eruptions when X(t) is "climate", thus generalizing the linear model of Hasselmann (1976), or cloud processes when dealing with the weather itself: one person's signal is another person's noise, as the saying goes. In the case of random forcing of Eq. (20), the concepts introduced by the simple deterministic examples of Eq. (10b) and Eqs. (15) and (16) above can be illustrated by the random attractor (ω) in Fig. 7.
A key feature of the pullback point of view on noiseperturbed dynamical systems that characterizes RDS theory is the use of a single noise realization, as opposed to the traditional, forward viewpoint of the Fokker-Planck equation and associated concepts, in which multiple noise realizations play a role. For a precise definition of a random attractor -as well as the commonalities and differences between the deterministic and random cases of time-dependent forcingplease see Caraballo and Han (2017).  studied a specific case of such a random attractor for the paradigmatic, climate-related Lorenz (1963) convection model. The authors introduced multiplicative noise into each of the ODEs of the original, deterministically chaotic system, as shown below: where r = 28, P r = 10, and b = 8/3 are the standard parameter values for chaotic behavior in the absence of noise and σ is a constant variance of the Wiener process that is not Figure 7. Schematic diagram of a random attractor (ω) and of the pullback attraction to it; here ω labels the particular realization of the random process θ (t)ω that drives the system. We illustrate the evolution in time t of the random process θ (t)ω (light solid black line at the bottom), the random attractor (ω) itself (yellow band in the middle), with the snapshots A 0 (ω) = A(ω; t = 0)) and A(ω; t) (the two vertical sections, heavy solid lines), and the flow of an arbitrary compact set B from "pullback times" t = −τ 2 and t = −τ 1 onto the attractor (heavy blue arrows). See Ghil et al. (2008, Appendix A) for the requisite properties of the random process θ (t)ω that drives the RDS formulated by Eq. (20). Adapted from Physica D: Nonlinear Phenomena, 237, Ghil et al. (2008), Climate dynamics and fluid mechanics: natural variability and related uncertainties, 2111-2126, © 2008 Elsevier B.V. All rights reserved. 2008, with permission from Elsevier. necessarily small; both the noise realization η(t) and σ are the same in all three equations. The well-known strange attractor of the deterministic case is replaced by the Lorenz model's random attractor, dubbed LORA by the authors. Four snapshots A t (ω) of LORA are plotted in Fig. 8 and a video of its evolution in time (ω) = {A t (ω)} t∈R is available as supplementary material to  at https://doi.org/10.1016/j.physd.2011.06.005. Charó et al. (2021) have further analyzed the striking effects of the noise on the nonlinear dynamics that are visible in Fig. 8 and in the video of  and gathered further insights into the abrupt changes of the snapshots' topology at critical points in time. These remarkable changes suggest the possibility of random processes giving rise to qualitative jumps, such as the MPT, in paleoclimatic variability.

Application to Dansgaard-Oeschger (D-O) events
Before discussing conceptual glacial cycle models, we take a little detour and introduce a simpler -yet interesting and at the same time highly instructive -application of NDS theory to another important climate phenomenon. During past glacial periods, Greenland experienced a series of sudden decadal-scale warming events that left a clear trace in ice core records (Dansgaard et al., 1993). These so-called D-O events were followed by intervals of steady moderate cooling, be- fore a short phase of enhanced cooling brought the temperatures back to their pre-event levels (e.g., Rasmussen et al., 2014). This pattern is very clearly apparent in NGRIP δ 18 O records (North Greenland Ice Core Project members, 2004), and it can be qualitatively simulated by the fast component of a FitzHugh-Nagumo (FHN) model (FitzHugh, 1961;Nagumo et al., 1962), as pointed out by Mitsui and Crucifix (2017), among others; see also Kwasniok (2013), Rial and Yang (2007), Roberts and Saha (2017), and Vettoretti et al. (2022).
We discuss the example of the FHN model at some length in order to illustrate how external forcing can act on a system's internal variability and thereby give rise to more complex dynamics. This model's concise mathematical formulation and its widespread application in paleoclimate modeling and other fields make it ideally suited for this goal. We start with a description of the autonomous model with no timedependent forcing. Subsequently, we introduce a simple sinusoidal forcing and numerically compute the corresponding PBA. We then extend these consideration into the realm of random dynamical systems by adding stochastic forcing and discuss the resulting random attractor. Finally, we replace the synthetic forcings by one that corresponds to a paleoclimate proxy record of past CO 2 concentrations retrieved from Antarctic ice cores (Bereiter et al., 2015a) and show that this setup brings the model's trajectories into good qualitative agreement with the D-O patterns observed in δ 18 O records from Greenland ice cores. In doing so, we pay less attention to the physical interpretation of the model's variables, while focusing on the detailed explanation of model behavior and on the role of the forcing in the resulting dynamics.

The FitzHugh-Nagumo (FHN) model of fast-slow oscillations
The FHN model consists of two coupled ODEs that govern behavior alternating between slow evolutions and fast transitions. Typically, the timescales of the two variables are separated by introducing the parameters τ x and τ y , with x(t) being the slow component and y(t) the fast one: In order to develop an understanding for the way that such a model can simulate the rapid D-O warmings, followed by slow coolings, we start by discussing its autonomous behavior for time-independent γ . First, consider the case of large timescale separation τ x τ y .
This choice guarantees that the fast y component adjusts adiabatically to quasi-static changes of the slow x component. The time derivative of y(t), as shown in Fig. 10a, exhibits either three real roots or a single one, depending on the value of x(t), which shifts the graph of the cubic polynomial P 3 (x, y) = α(y −y 3 )−x globally upwards or downwards. Of the three potential roots, the outer two are stable fixed points for y(t) at a fixed value of x, while the inner one is unstable. Note that the two stable fixed points are always located either left or right of the local minimum or maximum of P 3 (x, y), respectively. Accordingly, we label them y l and y r . The positions of the local extrema, namely y min = − √ 1/3 and y max = √ 1/3, provide an upper and a lower bound for the left and right stable fixed point, respectively. Now, let us investigate the coupled dynamics of the slow and fast variables x(t) and y(t). Assume we are in a state where x > 2α/ √ 27 such that y l is the only root of P 3 (x, y). Provided that |γ | < √ 1/3, the time derivativeẋ and x itself will have opposite signs, with the consequence that a slow adjustment process of x(t) sets in. This will shift the polynomial P 3 (x, y) upwards in Fig. 10a, which in turn revives the other two roots of P 3 (x, y). The fast y(t) closely follows y l , which is shifted accordingly. Once x(t) = −2α/ √ 27, the root y l stops existing and a fast transition to the neighborhood of the opposite root y r ensues. This transition switches the sign ofẋ(t) and the slow adjustment of x(t) now happens in the opposite direction.
In the (x, y) plane, this behavior manifests itself as a stable limit cycle. However, any value |γ | > √ 1/3 preventṡ x(t) from switching sign and therefore interrupts the cyclic destruction and revival of the respective opposite root of P 3 (x, y). Instead, |γ | > √ 1/3 gives rise to a single stable fixed point for the entire system in the (x, y) plane and in this case both variables relax towards it.
In fact, in the autonomous setting, the system's qualitative behavior is controlled by the value of the parameter γ , which decides between internal oscillations along a limit cycle and relaxation towards a fixed point in the (x, y) plane. Note that previously we referred to the roots y l and y r also as stable fixed points of y for a given value of x. Here, the term stable fixed point refers to the entire system defined by the coupled ODEs (22a)  This behavior can be better understood by considering the nullclines of Eqs. (22b) and (22a) in the (x, y) plane, as shown in Fig. 9. If the branches of the y nullcline that correspond to y l and y r , and thus to stable fixed points of Eq. (22b) for a given value of x, intersect with the x nullcline given by y = γ , then this intersection constitutes a stable fixed point for the entire system.
If they do not, the system first relaxes along the fast direction toward the y nullcline. Only then does the adjustment of the slow component start to drag the system along the y nullcline in the direction where the distance to the x nullcline decreases. At the point where the y nullcline reverses, the fast component is immediately attracted by the other branch of the fast nullcline and the same process starts all over again. So far we have described the formation of the limit cycle in the FHN model under the assumption of clear timescale separation and the independence of the x nullcline {y = γ } from x. See Rocsoreanu et al. (2000) for the emergence of the limit cycle in the more general FHN model.
The highly nonlinear, two-time behavior of the FHN model somewhat modifies the way that stable limit cycles arise in it. While we saw the oscillation's radius grow with the square root of the bifurcation parameter in the case of the normal form given by Eqs. (3)-(5), in the case of the FHN model, the radius of the oscillations actually grows exponentially over a small range of γ values right after the bifurcation point and then stabilizes. These exponentially growing stable limit cycles have been termed "canard cycles" (Benoît, 1983). Roberts and Saha (2017) have pointed out a possible link between canard cycles and D-O cycles, and they play a role in other excitable climate models; see Pierini and Ghil (2022) and references therein.

A pullback attractor of a periodically forced FHN model
Introducing a sinusoidal time dependence into the slow Eq. (22a) makes the system non-autonomous, as discussed in general terms in Sect. 3.2, and it periodically switches the self-oscillatory behavior on and off. We consider here the case in which the variations of the external forcing γ (t; τ f ) occur on a slower timescale than the entire internal dynamics, i.e., Note that we give up here on the restriction of strict timescale separation between x(t) and y(t), as expressed by Eq. (23). The trajectories plotted in Fig. 10c-e and f-h represent the solutions of such a periodically forced FHN model starting at different points in the past. The time units are arbitrary and the two sets of solutions are for two different forcing timescales τ f . These trajectories illustrate the applicability of the pullback perspective suggested by Sect. 3.2 and Fig. 6 to the periodically forced FHN model. In contrast to the illustrative examples of Sect. 3.2, an analytical solution is not available in this case. However, this small sample of numerically computed trajectories, together with the phenomenological discussion herein, is quite sufficient for an intuitive understanding of the system's pullback behavior.
For Fig. 10c-e, the forcing's timescale is much slower then the internal timescales. The amplitude of the forcing γ (t) was chosen such that it repeatedly crosses the two thresholds ± √ 1/3 and thus induces a sequence of Hopf bifurcations by switching between intervals of self-sustained oscillation and attraction to a stable fixed point. Crossing such a bifurcation point due to slow changes in the forcing is referred to as a bifurcation-induced tipping or "B-tipping" (Ashwin et al., 2012;Ghil, 2019;Ghil and Lucarini, 2020).
Strikingly, all trajectories converge to one another during non-oscillatory time intervals, when they are simultaneously attracted by the single existing fixed point. During oscillatory intervals, phase differences between individual trajectories may, in principle, persist. Still, convergence during a single non-oscillatory interval is so strong that after it the trajectories can no longer be discriminated visually. Numerically, however, the distance between trajectories only tends to zero but never reaches it. At the end of non-oscillatory intervals, the trajectories always re-enter the oscillatory regime from the same location in the (x, y) plane -to within negligible numerical differences -and hence very nearly repeat themselves. Qualitatively speaking, the PBA -i.e., the family of invariant snapshots {A(t)} t∈R of Eq. (13) -is an infinite repetition of the common trajectory structure that can be observed in Fig. 10d, e between −5000 and 15 000 time units. In the case at hand, each snapshot A(t) consists of a single point. (a-c) One illustrative trajectory (red dotted) for γ = 0, 0.5, and 1.0, respectively. The upper branch and the lower branch of the y nullcline (solid blue) correspond to the roots y r and y l , respectively, as defined in the text. For fixed values of x, they constitute stable fixed points for Eq. (22b). The middle branch (dashed blue) corresponds to the unstable fixed point for Eq. (22b) for a given x. The trajectories for panels (a) and (b) follow a limit cycle, as the x nullcline does not intersect with either of the stable branches of the y nullcline, given that |γ | < √ 1/3, while the trajectory in panel (c) approaches a single stable fixed point for the coupled system formed by the intersection of the x nullcline with the y l branch of the y nullcline.
For Fig. 10f-h, the timescale separation between the forcing and the internal dynamics is reduced, resulting in a qualitatively different behavior of the non-autonomous system. The frequency of occurrence of B-tipping points is much higher, and hence the trajectories do not even execute a full oscillation during a single time interval that permits oscillations. As a result, two stable patterns of trajectories are formed. These two patterns can be brought into agreement by switching the sign of one pattern and shifting it in time by τ f /2. This symmetry reflects the symmetry of the stable nullcline of the fast system component, as shown in Fig. 9.
Again, the PBA of this non-autonomous system can be thought of as an infinite repetition of the common trajectory structure that can be observed in Fig. 10g, h between −5000 and 15 000 time units. In contrast to the slow-forcing case, each snapshot A(t) now is comprised of two points in the (x, y) plane. This example illustrates how the action of an external force on an autonomous system can give rise to considerably richer dynamics, which crucially depends on both the system's internal variability and the nature of the forcing.

A random attractor of a periodically and stochastically forced FHN model
Based on the brief introduction to RDSs in Sect. 3.2, we take our investigation of the periodically forced FHN model one step further and include a random component into the external forcing, acting on the model's fast y component: Here, η denotes a Wiener process, as in Eq. (20) -i.e., a continuous stochastic process whose increments dη are inde-pendently and normally distributed, with mean zero and unit variance -and γ (t) remains a periodic forcing of the x component proportional to sin(τ f t).
In order to study the random attractor of this system, we compute trajectories with random initial conditions over a time span long enough to reveal the asymptotic behavior, as shown in Fig. 11. For both the slower and faster deterministic forcing -i.e., τ f = 1000 and τ f = 350, respectively, as studied in the previous paragraphs -several random attractors are approximated for increasing noise variance values σ in Eq. (26b). For each attractor approximation, we use 20 trajectories with the same noise realization. Each random attractor (red) is shown together with the corresponding PBA of the FHN system subject to purely periodic forcing (blue and green).
For the long periodic forcing with τ f = 1000, the trajectories in Fig. 11a converge fairly rapidly to a single one, as in the case with no noise in Fig. 10d, e. Furthermore, there is a clear similarity of pattern and proximity in phase between the PBA of the deterministic system and the random attractor. However, the deviations of the random attractor from the PBA increase in both pattern and phase as the noise variance increases. These deviations are most striking during the oscillatory intervals because the noise can induce phase shifts in the oscillations, which then persist for the duration of the oscillatory interval. During non-oscillatory intervals, the random attractor is less susceptible to the noise because the resulting perturbations decay in the presence of a stable fixed point of the underlying deterministic system.
For the shorter period forcing with τ f = 350, the PBA of the deterministic system features two stable branches, which could in principle persist in the random attractor. As for the case of τ f = 1000, a rapid convergence of the trajectories can be observed for all noise variances. However, two separate 27. The purple lines labeled y min and y max mark the right and left boundaries for the roots y l and y r , respectively: y l and y r can never be located in between the two purple lines. (b) Trajectories of the nonautonomous model with γ (t) = 0.8 sin(t/τ f ) and τ f = 1000, plotted in the (x, y) phase plane; the trajectories are colored by their starting times {t 0 = −20 000, t 1 = −16 000, t 2 = −13 000, t 3 = −7000}, and the initial positions were drawn from a standard Gaussian bivariate distribution. (c) The slow time-dependent forcing γ (t) = 0.8 sin(t/1000). (d, e) The same trajectories as in (b) but plotted in time as y = y(t) and x = x(t), respectively. Panels (f-h) are same as panels (c)-(e) but for the fast time-dependent forcing γ (t) = 0.8 sin(t/350). The gray shading in panels (c)-(h) indicates intervals during which |γ | > √ 1/3, and the internal oscillation is hence suppressed.
bundles of trajectories persist, each associated with one of the two separate branches of the deterministic PBA. Only after some time do the two bundles merge and subsequently follow either one of the two branches of the PBA or the other. This phenomenon is called noise-induced synchronization (e.g., Arnold, 1998;. Which of the two PBA branches the random attractor follows depends on the exact noise realization. The random attractor may also switch irregularly from one branch of the PBA to the other (not shown). Already in the simple setup explored here, one notices that the higher the noise variance, the faster the synchronization of the trajectories. This statement must be understood probabilistically: it may certainly happen for two given noise realizations that the random attractor with the higher noise level takes longer to synchronize, as is the case for σ = 0.7 and σ = 0.8 in Fig. 11b. The investigation carried out herein assesses a very special case of an FHN model's random attractor. Random attractors of FHN-type models have been studied intensively (e.g., Wang, 2009;Yamakou et al., 2019). Most studies so far, however, have concentrated on the excitable regime of the stochastically forced FHN model, i.e., the deterministic model's parameters were chosen so as to exhibit a single stable fixed point in the absence of the noise. The noise variance was then chosen sufficiently high to ensure that random fluctuations can push the system out of equilibrium and allow it to take a round trip on what would be the deterministic model's limit cycle, given different parameter values. forcing with period τ f = 1000. The top graph shows the periodic forcing, with non-oscillatory intervals marked by gray shading. The next three graphs show the fast y component of the approximate random attractors, as per Eq. (26b); the noise variance σ increases from top to bottom. Each random attractor is approximated by integrating 20 trajectories (solid red) with different initial conditions over time and using the same Wiener process as their common stochastic forcing; the corresponding deterministic PBA is shown in blue. The random attractors and the PBA are very similar for small noise variances, but they differ more and more as the noise variance increases. (b) Results for faster periodic forcing with period τ f = 350. The forcing is only shown for the first 8000 time units, with non-oscillatory intervals again shaded gray. The left part of panel (b) shows five approximate random attractors, computed as in panel (a), on a common time axis. The panel's righthand side shows their continuations on individual time axes in order to display the moment when full noise-induced synchronization of the trajectories takes place. Prior to this point, the random attractor is split into two branches, which closely follow the PBA of the deterministic system (blue and green). For σ = 0.7, the synchronization already takes place during the first 8000 time units.
The situation studied here is rather different, as shown across Figs. 9-11.

An FHN model of the NGRIP record
Readers who are familiar with the NGRIP δ 18 O record (North Greenland Ice Core Project members, 2004) might have already realized the qualitative resemblance between the proxy data and the fast component's trajectory of the periodically forced FHN model in Fig. 10d. In particular, the prominent sawtooth pattern of the data is satisfactorily captured by the fast and slow dynamics of the model. Figure 12 shows a trajectory of the FHN model for which the sinusoidal forcing used in Fig. 10 was replaced by a rescaled time series of atmospheric CO 2 concentrations re-trieved from Antarctic ice cores (Bereiter et al., 2015a): Is it remarkable how well this simple forcing brings the oscillatory intervals of the FHN model into agreement with the time intervals of the record that are dominated by D-O cycles without any systematic tuning of the model parameters. Clearly, the CO 2 -forced FHN model fails to reproduce the exact waiting times between D-O events. However, with these waiting times being at least in part stochastically determined (Ditlevsen et al., 2007), the purely deterministic FHN model is not meant to reproduce the exact pattern of D-O events. Vettoretti et al. (2022) have carried out a detailed In fact, Rousseau et al. (2022, Fig . 6) describe in detail a somewhat more complex, proxy-record-based picture of the interaction between episodes full of D-O events, (Heinrich, 1988) events, and longer-term cooling trends. It is quite possible that a simple model like the one in this section but including explicitly continental ice sheets could capture such a detailed picture.
In the present framework, the FHN model's fast variable y(t) may be interpreted as the intensity of the Atlantic meridional overturning circulation (AMOC), which switches between on and off states during self-oscillatory behavior; see, for instance, Henry et al. (2016), Ghil (1994, Table 5), and Ghil and Lucarini (2020, Table I). The slower x variable that drives the transition between the on and off states of the AMOC may then be taken, for instance, as the waxing and waning of Northern Hemisphere ice sheets (e.g., , linked in turn to varying ice shelf extent (e.g., Boers et al., 2018), or as the weakening and strengthening of Antarctic bottom water production (Vettoretti et al., 2022).
The interaction between the fast variable and the slow one happens here in the presence of a climate forcing represented by CO 2 concentration. On the much slower timescales of Quaternary glaciations, an interplay between the CO 2 concentration and mean global temperature might also occur, as we shall see in the next section.

An NDS for the Quaternary glaciations
Apparently, it was Crucifix (2013) who first applied pullback ideas to the problems of Quaternary glaciations, independently of earlier work on the topic in the climate literature (Ghil et al., 2008;Bódai and Tél, 2012). His work concentrated mainly on the connection between the pacemaking role of the orbital forcing and the observed irregularity of the glacial terminations during the Late Pleistocene; cf. Broecker and Van Donk (1970) and Ghil and Childress (1987, Fig. 11.2). Based on the considerable success of NDS and RDS applications to other climate problems -such as ENSO (Ghil et al., 2008;Ghil and Zaliapin, 2015;Chekroun et al., 2018;Marangio et al., 2019), the wind-driven ocean circulation (Pierini et al., 2016(Pierini et al., , 2018, or the evaluation of the ensemble simulations routinely performed in support of the Assessment Reports of the Intergovernmental Panel on Climate Change (Drótos et al., 2015;Vissio et al., 2020) -it would appear worthwhile to proceed further along these lines.
Box 2 summarizes open questions with respect to Quaternary glaciations whose investigation should be aided by NDS and RDS theory. Of course, each of these problems requires one or more distinct climate models, as well as very careful modeling of the kinds of time-dependent changes in forcing and parameters that are most enlightening, as well as most relevant and plausible. A good way would be to start testing ideas with relatively simple models and pursue the investigation systematically across a hierarchy of models -through intermediate ones and on to the most detailed ones -in order to further increase understanding of the climate system and of its predictability on the various paleoclimatic time scales mentioned in the list above.
Such an approach can usefully complement the more common one of merely pushing onwards to higher and higher model resolution in order to achieve ever more detailed simulations of the system's behavior for a limited set of semiempirical parameter values. Ghil (2001) and Held (2005), among others, have emphasized the need to pursue such a model hierarchy, as originally proposed by Schneider and Dickinson (1974), in order to balance the need to broaden the number of plausible hypotheses vs. the need for confronting them with spatiotemporal details derived from observations.
In this section, we illustrate how the PBA concept can help shed more light upon the dynamics of ice age models. As pointed out in Sect. 3.1 and elsewhere in this paper, there is a long history of modeling the climate of the Quaternary by means of conceptual models, and many non-autonomous models have been proposed to simulate glacial-interglacial cycles of the last 400 kyr to 2.6 Myr based on the orbital forcing. In Appendix A, we provide a long but still not exhaustive list of glacial-cycle models and specify some of their key characteristics, including the degree of their success at simulating the MPT; see also the discussion in Sect. 2.2.
Among these glacial-cycle models, the model of Daruka and Ditlevsen (2016, DD16 hereafter) belongs to the more abstract ones, as it is not derived from detailed physical considerations. Still, its concise form, interesting nonlinear dynamics, and ability to simulate glacial cycles, as well as the MPT, make the DD16 model well suited for our illustrative purposes. We first modify this model slightly from its original formulation. We do so mainly in order for the model to better approximate the benthic δ 18 O proxy reconstruction of glacial-interglacial cycles due to Lisiecki and Raymo (2005a), especially the timing of glacial terminations; compare our Fig. 13 with Fig. 1 in DD16. Thereafter, we compute the PBAs of the modified DD16 model (hereafter M-DD16) to investigate the dynamical stability of its glacial cycles over the past 2.6 Myr.
Our model's variables, following DD16, are a global temperature anomaly y that is proportional to minus the global ice volume and an effective climatic memory term x that represents the internal degrees of freedom. In the deterministic case, the governing equations of the M-DD16 model are given by In the original DD16 model, MPT-like behavior was produced by a slow sigmoid variation of the parameter κ in Eq. (28b), In our M-DD16 model, we introduce instead a slow change in the parameters α(t) and β(t) of Eq. (28b) as follows: The functions α(t) and β(t) so defined are plotted in Fig. 13b, and they induce, as we shall see forthwith, a change in model behavior that not only resembles the MPT but also shows correct timings for most of the terminations. Moreover, to simulate δ 18 O model , we add a linear trend to the slow variable y to mimic the overall cooling at time scales of millions of years, and thus δ 18 O model = 4.3 − 1.4y + 0.0003t. Equation (28) generalizes a first-order ODE system that is equivalent to the Duffing form of the nonlinear spring Eq (7) discussed in Sect. 2.2. The use of such classical first-order systems -like the Duffing and Van der Pol ones -in paleoclimate models was initiated by Saltzman et al. (1981); see also Sect. 2.2 herein. ; the latter dataset, provided by Lisiecki and Raymo (2005a), is also known as the LR04 stack. Figure 13c shows a time series of simulated glacialinterglacial changes δ 18 O model (red) in comparison with the benthic δ 18 O (blue) of Lisiecki and Raymo (2005a). The model's initial condition is taken to be x = −1 and y = 0 at t = 10 000 kyr b2k and, since the insolation forcing in Fig. 13a is prescribed as a time series with 100-year sampling, we solved Eq. (28) using Heun's predictor-corrector method (Isaacson and Keller, 2012, chap. 8) with a step size of 100 years. A large spin-up time is chosen to guarantee that transient effects caused by the initial conditions have abated by the year 2.6 Ma b2k, which is the starting point for the time interval under study. The correlation between the model simulation and the proxy record is 0.75 for the time interval from 2600 to 0 kyr b2k and 0.72 over the interval from 1000 to 0 kyr b2k. Varying the parameters slowly across the time interval of interest, as shown in Fig. 13b, leads to a change in the frequency -from a dominant 41 kyr periodicity prior to the MPT, at roughly 1.2 Ma kyr b2k, to a dominant 100 kyr periodicity after the MPT, at roughly 800 kyr b2k -and a substantial increase in the amplitude.
We next approximate the PBA by taking 40 random initial conditions at 10 Ma b2k and integrating the model of Eqs. (28) and (30) up to the present time. The PBA in this case is simply a moving fixed point, as plotted in Fig. 14a, since the model dynamics is predominantly stable in the long time interval prior to the MPT. It would thus appear that the orbital forcing simply moves this fixed point around and fully determines Earth's climate. This agrees with the clear statement in DD16 that "first and foremost, our model does not have any internal periods of oscillation".
However, when keeping the parameters α and β fixed at their post-MPT values α = 0.7 and β = 3.9 throughout the simulation interval and repeating the computation of the PBA, a more complex picture arises. In the latter case, Fig. 14b shows a bunching of trajectories into separate fuzzy clusters, subject to the quasi-periodic orbital forcing of Fig. 13a.
There are two interesting inferences to be drawn. First, post-MPT dynamics is much more irregular and unstable than the more stable dynamics prior to the MPT. The robustness of the 40 kyr glacial cycles and instability of 100 kyr glacial cycles against perturbations is in line with the conclusions of previous studies (Mitsui et al., 2015;Quinn et al., 2018). It also appears to be consistent with the Willeit et al. (2019) simulation of the last 3 Myr of Earth's history that used an Earth system model of intermediate complexity (CLIMBER-2) that included ice sheets and a carbon cycle, along with atmosphere-ocean dynamics. In their work, trajectories starting from different initial states tended to converge to a single attracting trajectory in the Early Pleistocene,  Fig. 13, with the slowly changing parameters α(t) and β(t) given by Eq. (30). (b) The PBA for α(t) and β(t) kept constant at the post-MPT values of α = 0.7 and β = 3.9. In panel (b), the PBA is shown over a shorter time interval of the last 1000 kyr so that the detailed structure is more clearly visible. Without changes in α and β in time, the overall structure of the PBA is similar before and after 1000 kyr b2k.
while several distinct trajectories survived in the Late Pleistocene after the MPT.
Second, the separate bundles or "ropes" of trajectories in Fig. 14b seem to point to the type of generalized synchronization discussed in the paleoclimate context by De Saedeleer et al. (2013) and in the context of interannual and interdecadal climate variability by Pierini and Ghil (2022) and Vannitsem et al. (2021). Generalized synchronization in the strict sense of the existence of a map between a time-dependent control and the system's asymptotic behavior has only been shown to hold for non-chaotic systems. Work is under way, however, to further generalize this concept to chaotic systems as well (e.g., Rulkov et al., 1995;Zhang et al., 2007).

Conclusions
In this review and research paper, we have covered the contributions of the 1970s to the rebirth of the Milankovitch (1920) theory of the ice ages in Sect. 1 and the 1980s advances in modeling the Quaternary climate's intrinsic variability in Sect. 2. In Sect. 3, we presented first results on the effects of the orbital insolation forcing of the data discussed in Sect. 1 on the intrinsic variability of the data discussed in Sect. 2, and proceeded to introduce the novel concepts and tools of the theory of non-autonomous and random dynamical systems (NDSs and RDSs) that can help to better model and understand these effects. Section 3 concluded with the formulation and study of a FitzHugh-Nagumo (FHN)-type model of recurrent Dansgaard-Oeschger (D-O) events, in which historical CO 2 concentrations induced episodes of D-O events alternating with episodes of their absence in excellent qualitative agreement with NGRIP δ 18 O data; see Fig. 12.
Finally, in Sect. 4, we listed a number of open issues on Quaternary and longer paleoclimate timescales and proposed to address them by using the tools of Sect. 3.2. This approach was illustrated by a Duffing-type model of Daruka and Ditlevsen (2016), modified to include slow changes in the parameters that mimic such changes in the Earth system over the duration of the Quaternary period.
When the parameters are gradually changed in time so as to exhibit the mid-Pleistocene transition (MPT), the PBA is simply a moving fixed point. However, when the parameters are fixed at their post-MPT values, the PBA so obtained is chaotic and exhibits clusters of trajectories that we termed ropes. This suggests (a) that the stability of the system is gradually lost while crossing the MPT and (b) that the Late Pleistocene climate, albeit chaotic, may well be subject to a kind of generalized synchronization (cf. De Saedeleer et al., 2013;Pierini and Ghil, 2022;Vannitsem et al., 2021) with the orbital forcing that is illustrated in Fig. 5 of Sect. 3.1 herein. In the specific situation at hand, separate ropes may be associated with various combination tones of the forcing frequencies.
In a broader perspective -and leaving aside various finer points of the MPT conversation outlined in Sect. 2.2 -one can see the work that was reviewed and extended in this paper as a confirmation of the fine intuition of Emiliani and Geiss (1959), 6 decades ago, as summarized in and further expanded by Ghil and Childress (1987, Sect. 12, pp. 446-447).
Hence the following scenario (compare Emiliani and Geiss, 1959) suggests itself for the successive climatic transitions from Pliocene to Pleistocene and from Early to Late Pleistocene: As land masses moved towards more northerly positions, small ice caps formed on mountain chains and at high latitudes. These ice caps, due to their feedback on albedo, made climate more sensitive to insolation variations than it was in the total absence of ice. The response of the climatic system to such variations during the Early Pleistocene (2000 [kyr]-1000 [kyr] ago) was still relatively weak, of a fraction of a degree centigrade in global temperature perhaps, in agreement with the quasi-equilibrium results of Sect. 10.2.
As ice caps passed, about 1000 [kyr] ago, a certain critical size, the unforced system jumped from its stable equilibrium to its stable limit-cycle state (Figures 12.5 and 12.9), increasing dramatically the climate's total variability, to a few degrees centigrade in global temperature. Furthermore, resonant response became possible (see also Oerlemans (1984] and Sergin (1979)), enhancing abruptly the amplitude of the peak at 100 kyr, among others.
The take-home message is that slow and fast processes, both intrinsic and extrinsic, interact on all paleoclimatic timescales and that we are mastering the art of modeling such interactions.

Appendix A: Low-order dynamical-system models of glacial cycles
The dynamical modeling of glacial cycles dates back to the 1970s. Calder (1974) proposed a model of global ice volume changes that had different sensitivities to the insolation when ice sheets were waxing and waning, respectively. His model can be written as an NDS, according to Paillard (2001). Subsequently, conceptual dynamical models were further developed (e.g., Imbrie and Imbrie, 1980;Berger, 1999). Some of the more recent models simulate the proxy records of glacial cycles remarkably well (Paillard, 1998;Imbrie et al., 2011;Parrenin and Paillard, 2012).
Shortly after Calder (1974) presented his work, Weertman (1976) proposed a simple ice sheet model based on the flow law of a perfectly plastic solid. Next, researchers extended this simple ice sheet model by coupling it with an energy balance model (Källén et al., 1979) and further with the isostatic response of the underlying bedrock (Oerlemans, 1980;Ghil and Le Treut, 1981;Le Treut and Ghil, 1983;Pollard, 1983). Källén et al. (1979) found self-sustained oscillations in their simple coupled ice sheet-energy balance model. Le Treut and Ghil (1983) showed that the dominant 100 kyr periodicity of glacial cycles is generated -in their simple oscillator coupling ice sheet volume with the bedrock's isostatic rebound, on the one hand, and with the atmosphere and ocean's energy balance, on the other hand -via nonlinear resonance with the multi-periodic orbital forcing. More recently, Verbitsky et al. (2018) developed a simple physical model through a scaling argument that respects the underlying physics. Another branch of simple models explicitly includes the carbon cycle as an essential ingredient Paillard and Parrenin, 2004;Hogg, 2008;Toggweiler, 2008;Omta et al., 2016;Talento and Ganopolski, 2021).
A deeper understanding of glaciation cycles cannot be obtained without process-based models that focus on the detailed physics and biogeochemical phenomena involved Ganopolski and Calov, 2011;Abe-Ouchi et al., 2013;Ganopolski and Brovkin, 2017;Willeit et al., 2019). Still, simple dynamical models like the ones mentioned above, as well as models based on more mathematical considerations, are also useful for understanding the climate system's behavior and changes therein, since complex systems can sometimes exhibit familiar dynamics, regardless of the details (Nicolis and Nicolis, 2012;Crucifix, 2011). As Henri Poincaré pointed out, "mathematics is the art of giving the same name to different things".
For example, coupled nonlinear oscillators frequently exhibit synchronization with simple frequency ratios, either with each other (Pikovsky et al., 2001) or with the forcing (Ghil and Childress, 1987). Thus, mathematical models that ignore many physical details may also help us elucidate emergent properties in paleoclimatic dynamics (Tziperman et al., 2006;Crucifix, 2011De Saedeleer et al., 2013;Mitsui et al., 2015;Ashwin and Ditlevsen, 2015;Daruka and Ditlevsen, 2016) or test different orbital hypotheses (Huybers, 2011). Many low-order models of glacial dynamics are listed in Table A1, although it is by no means exhaustive.

Appendix B: PBA for a limit cycle with sinusoidally modulated radius
Here we study the system of two formally decoupled ODEs that were introduced in Sect. 3.2 and analytically derive their invariant sets as well as the corresponding pullback attractor (PBA). Following Crauel and Kloeden (2015), the PBA is given by the family First, we define ρ(t) = ρ(t) − µ, which gives rise to This is an inhomogeneous ODE and can thus be solved by the variation of parameters method (e.g., Boyce and DiPrima, 2005). The ansatz  His model can be written as a non-autonomous dynamical system (NDS) (Paillard, 2001). Weertman (1976) Simplified ice sheet model based on the flow law of a perfectly plastic solid.

Non-chaotic.
No selfsustained oscillation in the absence of forcing.
Chaotic in the presence of orbital forcing (Ghil, 1994). Exhibits self-sustained oscillations with a period of ≈ 10 kyr; the basic mechanism for getting ∼ 100 kyr cycle is nonlinear resonance at a difference tone between the 19 and 23 kyr precession cycles (Le Treut and Ghil, 1983); see Sect. 2.1 herein.
Not explicitly. It can exhibit, however, the MPT in principle, since 41 kyr oscillations and dominant 100 kyr oscillations can be generated depending on the model's parameter values. Källén et al. (1979) is an ice sheet-temperature coupled oscillator. Ghil and Le Treut (1981) introduced the bedrock and Le Treut and Ghil (1983) added the orbital forcing. Imbrie and Imbrie (1980) Piecewise linear model with two different timescales for ice sheet waxing and waning.

Non-chaotic.
No selfsustained oscillation in the absence of forcing.
Non-chaotic. After the MPT, it exhibits self-sustained oscillations with a period of ∼ 100 kyr in the absence of forcing.
Yes, the MPT arises via a Hopf bifurcation in the underlying system; see also Crucifix (2012).
It has a strange non-chaotic attractor, as well as a chaotic one, depending on the parameter setting (Mitsui and Aihara, 2014). Paillard (1998), Paillard (2003, 2012) Hybrid dynamical system with discrete states that switch when conditions are satisfied.
Non-chaotic but the dynamics is sensitive to parameter changes near the switching boundaries (Paillard, 2001).
No self-sustained oscillations in the absence of forcing.
Yes, in Paillard (1998). The authors call these models relaxation oscillators, in spite of their discrete-state formulation.
Berger (1999) Delayed differential equation for ice volume, Robustness of trajectories against random perturbations is mentioned. It exhibits damped oscillations in the absence of forcing.
Not specifically addressed. Paillard and Parrenin (2004) Coupled oscillator based on Northern Hemisphere ice volume-Antarctic ice extent-CO 2 coupling.
Non-chaotic. It exhibits selfsustained oscillations in the absence of forcing.
Yes, the MPT is induced by a slow drift in the bottom water formation efficiency around the Antarctic. Ashwin et al. (2018) showed that this model can exhibit chaotic dynamics when it is slightly modified. Ashkenazy and Tziperman (2004), Tziperman et al. (2006) An ice mass model that simplifies the sea ice switch model of Gildor and Tziperman (2000).
It exhibits self-sustained oscillations in the absence of forcing. Non-chaotic under the orbital forcing.
Yes, in Ashkenazy and Tziperman (2004), where the maximum ice volume threshold is increased in time according to the regolith hypothesis.
Uses the temperatureprecipitation feedback introduced by Källén et al. (1979). Huybers (2007Huybers ( , 2011 The ice mass grows monotonically and collapses to zero when it exceeds a threshold modulated by the obliquity (Huybers, 2007) or by a hybrid measure consisting of the obliquity and climatic precession (Huybers, 2011).
Non-chaotic. It exhibits selfsustained oscillations in the absence of forcing.
Yes, in Huybers (2007). Due to the presence of the threshold, the model is sensitive to changes in parameters or in its position in phase space.
Yes, MPT-like behavior is produced solely by changes in orbital parameters.
Crucifix (2012) Non-chaotic (Crucifix, 2012) but sensitive to the noise. It exhibits self-sustained oscillations in the absence of forcing.
Yes, MPT-like behavior is generated via a Hopf bifurcation with an explosive character (Crucifix, 2012). Ashwin et al. (2018) showed that this model has only a small parameter region corresponding to chaotic dynamics but may have a wider chaotic region when it is generalized to the VdP-Duffing system. Mitsui et al. (2015) One-dimensional phase oscillator model.

Non-chaotic.
Yes, MPT-like frequency change accompanies a smooth or non-smooth saddle-node bifurcation of tori.
It has a strange non-chaotic attractor or a quasi-periodic attractor in a classical sense, i.e., not in the pullback sense.
Non-chaotic; overall stability of simulated glacial cycles against dynamical noise is reported.
It exhibits MPT-like behavior via a trans-critical bifurcation of the slow manifold in the fast dynamics. Omta et al. (2016) Forced two-dimensional oscillator consisting of ocean alkalinity and calcifier population.
Non-chaotic or chaotic given periodic forcing, depending on the parameters. It exhibits self-sustained oscillations in the absence of forcing.
No self-sustained oscillations in the absence of forcing. Can be chaotic depending on parameter values; the authors call this the "climatic butterfly effect".
Yes, MPT-like behavior is induced by a slow change in the damping coefficient parameter κ.
See also Sect. 4 herein. (2017) Coupled system of ice volume, temperature, and atmospheric CO 2 incorporating a delayed CO 2 contribution from ocean ridge volcanism. It is reduced to a forced 1-D delay differential equation.

Huybers and Langmuir
A phase-locking property is reported. The underlying system has two stable and one unstable stationary states. However, it is close to a Hopf bifurcation point and is excitable by the forcing.
Yes, the MPT is modeled as a switch from smallamplitude oscillations to large-amplitude ones, which is triggered by the amplitude modulation of obliquity cycles. Quinn et al. (2018) Scalar delay differential equation for ice volume, derived from the Saltzman and Maasch (1988) model via the linear chain approximation.
Non-chaotic before the MPT and temporarily chaotic after the MPT around 800 kyr b2k.
In the absence of astronomical forcing, the delayed feedback leads to bistable behavior, in which stable large-amplitude oscillations and an equilibrium coexist.
Yes, MPT-like behavior is induced by summer insolation forcing as a transition from small-amplitude ∼ 41 kyr cycles to large-amplitude ∼ 100 kyr cycles. Ashwin et al. (2018) Study of several low-order dynamical systems models, including Paillard and Parrenin (2004), Crucifix (2012), and a generalized VdP-Duffing oscillator.
Chaotic or non-chaotic depending on the model and its parameter values.
Not explicitly discussed. The models can, however, exhibit MPT-like behavior in principle because 41 kyr oscillations and dominant 100 kyr oscillations are generated depending on parameter values.
The authors stress the possibility of chaotic dynamics occurring. Yes, MPT-like behavior is induced by an enhancement of positive feedbacks against negative feedbacks in the model. Talento and Ganopolski (2021) Ice volume-CO 2 coupled conceptual model.

Not specifically addressed.
This model was developed to simulate the last 800 kyr of glacial cycles and Earth's future climate on the millionyear timescale.
Low-order dynamical models in which stochastic processes are of the essence Benzi et al. (1981Benzi et al. ( , 1982, Nicolis (1981) Two independent formulations of a stochasticresonance model. Additive noise favors the bimodal response of a periodically modified double-well potential to global annual insolation variations.
Additive stochastic perturbations play a key role.
Not specifically addressed. An MPT-like transition could be obtained by subjecting both the depth of the deterministic double well and the pure periodicity of its modulation to a more or less gradual modification. Benzi et al. (1981Benzi et al. ( , 1982 and Nicolis (1981) proposed the stochastic-resonance idea for ice age cyclicity independently of each other. Matteucci (1989) Generalized stochastic resonance model with obliquity and precessional forcing.
Additive stochastic perturbations play a key role.
Not specifically addressed.
Pelletier (2003) Model based on the coherence resonance mechanism (Pikovsky and Kurths, 1997) uses a single temperature ODE with a discrete delay.
Additive stochastic perturbations play a key role.
Yes. The delay feedback in the temperature is based on ice sheet extent reconstructions and affects the model's albedo.
Additive stochastic perturbations play a key role.
Video supplement. The video supplement to this article (https://doi.org/10.5281/zenodo.6346211, Riechers, 2022) illustrates the pullback attractor (PBA) associated with the simple system governed by Eqs. (15) and (16). It shows a heat map of the phase plane, derived from an increasing number of trajectories with common initial time t i = −200 and final time t f = 200. The initial radius and phase of each trajectory are randomly sampled from Gaussian distributions centered at µ ρ = 20 and µ φ = 0 with standard deviations of σ ρ = 5 and σ φ = 10, respectively. Over the course of the video, 100 trajectories are continuously added to the heat map and the annular disk D = {(ρ, φ) : ρ ∈ [µ − αβ, µ + αβ] and φ ∈ [0, 2π)} fills up. The heat map in Fig. 6b is a snapshot from this video at time t = 0.2.
Author contributions. MG conceived and designed the study. KR and TM carried out the major part of the article's new research in close interaction with MG and NB. All authors interpreted and discussed the results and wrote the manuscript.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement.
This article is part of the special issue "A century of Milankovic's theory of climate changes: achievements and challenges (NPG/CP inter-journal SI)". It is a result of the conference "One Hundred Years of Milankovic's Theory of Climate Changes: synergy of the achievements and challenges of the next century", 17-18 November 2020.
to give a broader perspective of relevant work on the Quaternary's glacial cycles and the mid-Pleistocene transition (MPT The publication of this article was funded by the Open Access Fund of the Leibniz Association.
Review statement. This paper was edited by Marie-France Loutre and reviewed by Tamas Bodai and two anonymous referees.