Crossover and peaks in the Pleistocene climate spectrum; understanding from simple ice age models

The power spectrum provides a compact representation of the scale dependence of the variability in time series. At multi-millennial time scales the spectrum of the Pleistocene climate is composed of a set of narrow band spectral modes attributed to the regularly varying changes in insolation from the astronomical change in Earth’s orbit and rotation superimposed on a continuous background generally attributed to stochastic variations. Quantitative analyses of paleoclimatic records indicate that the continuous part comprises a dominant part of the variance. It exhibits a power-law dependency typical of stochastic, self-similar processes, but with a scale break at the frequency of glacial-interglacial cycles. Here we discuss possible origins of this scale break, the connection between the continuous background and the narrow bands, and the apparently modest spectral power above the continuum at these scales. We demonstrate that the observed scale break around 100 ka can have a variety of different origins and does not imply an internal time scale of correlation as implied by the simplest linear stochastic model.


Introduction
The schematic picture of the climate spectrum proposed by Mitchell (1976), reproduced in Fig. 1, provides a classical, standard view of climate variability. Mitchell made a distinction between the continuous stochastic background spectrum and the discrete (line) spectrum superimposed, comprised by the periodic diurnal -, tidal -, and annual cycles.
Obviously, Mitchell's figure is schematic and not intended to be interpreted too literally. For example, if we could compute the Fourier transform of temperature records available over the entire geological history, the spectral line of the diurnal cycle would be broadened due to the tidally-forced slowing down of the Earth's rotation. Likewise, the 100 ka cycle, which is characteristic of the Pleistocene (2.6 Ma-12 ka BP), would not appear clearly on the spectrum of a 1 billion-year-long record. In other words, the low frequency part of Mitchell's spectrum is not an average over geological times, but rather a depiction of the dynamics averaged over a characteristic time in the recent Earth history. Yet, it offers an important and key insight by combining two aspects of climate variability. One is associated with stochastic processes, where the accumulation of random fluctuations generates the background spectrum. The other aspect is associated with the periodic or quasi-periodic climate variations. These are the peaks on the spectrum.
The underlying theory and interpretation of Mitchell's background spectrum were given by Hasselmann (1976). The basic assumption is that fast (chaotic) climatic fluctuations (atmosphere temperature) decorrelate quickly enough to be effectively seen as noise by the slower components (ocean temperature) of the climate system. The influence of the noise on the slow climate variables is then integrated over time, counteracted by stabilising feedback processes. The latter brings the slow variables back to an equilibrium state with some characteristic relaxation time. In the Hasselmann picture, the fluctuations are white noise, and the relaxation process is linear. It is a standard example of an Ornstein-Uhlenbeck (OU) process. Its power spectrum is easily calculated, as we shall do shortly. It behaves as (P(f ) ∼ f −2 ) for frequencies larger than the inverse relaxation time and has a flat spectrum (P(f ) ∼ f 0 ) for frequencies smaller than the inverse relaxation time.
The origin of peaks at 20, 40 and 100 ka are, in Mitchell's spectrum, related to the astronomical forcing of climate. From celestial mechanics it is indeed established (Milankovitch 1941) that over time scales of tens of thousands of years, the precession of Earth's axis of rotation, variation in the inclination of the axis (obliquity) and in eccentricity affect the seasonal and latitudinal distributions of insolation. These orbital-or astronomical-variations constitute the dominant external forcing to the climate system. They are due to many-body perturbations of Earth's Keplerian orbit and steady rotation, and can be represented as a series of harmonic functions of time, with leading components having periods of 19, 21, 23 (precession), 41 (obliquity), 100 and 400 ka (eccentricity) (Berger 1978). Evidence that the astronomical forcing somehow controls the succession of ice ages became largely acknowledged when Hays et al. (1976) recognised the astronomical frequencies in the power spectra of marine isotopic records. Investigators have therefore focused on the mechanisms by which the astronomical forcing may generate, or at least set the timing, of ice ages. An approach has been to model the ice age phenomenon as governed by a determinstic dynamical system, featuring a small number of essential interacting climatic variables. A variety of models were proposed, involving interactions between ice volume/ ocean temperature/atmospheric CO 2 (Saltzman and Sutera 1987;Saltzman and Maasch 1990), ice sheet extent/planetary albedo/atmospheric CO 2 (Ghil et al. 1985), calcifying plankton/ocean alkalinity/CO 2 (Omta et al. 2013), and ocean temperature/ice sheet area/basal temperature (Verbitsky et al. 2018), to cite but a few. Other models are more conceptual in the sense that they do not immediately identify variables with components of the climate system, but provide more generic results about the response to the astronomical forcing of a non-linear multi-stable system (Paillard 1998;Ditlevsen 2009;Daruka and Ditlevsen 2015), bi-stable potential with stochastic forcing (Benzi et al. 1982), internal oscillator (De Saedeleer et al. 2013;Mitsui et al. 2015;Ashwin and Ditlevsen 2015;Nyman and Ditlevsen 2019), or a threshold crossing process (Imbrie and Imbrie 1979;Huybers and Wunsch 2005). Some models require the astronomical forcing for generating glacial-interglacial cycles, while others exhibit a self-sustained cycle which can be synchronized on the astronomical forcing. Yet, they all generate power spectra in which the power is narrowly concentrated in the region 10-100 ka, and exhibit pronounced spectral peaks.
Although the Milankovitch theory relating the glacial cycles to the narrow band astronomical forcing was confirmed, the spectral variance in the bands of astronomical forcing seems to capture only a limited fraction of the total variance of the record. This is seen in Fig. 2, with the spectra from the past 2 Ma of the Lisiecki and Raymo (2005) (LR04, green) and the Huybers and Wunsch (2004) (HW05, blue) records. These are composite records, obtained by averaging ("stacking") several individual records. The two stacks differ mainly in the dating assumptions: The LR04 is tuned to a simple climate model using the insolation curve, while the HW05 is dated using a constant sedimentation rate model independent from the astronomical forcing. As LR04 is tuned, it displays slightly more power in the obliquity and precession bands. Even with that tuning, though, the background spectrum contains most of the total variance. Before us, Wunsch (2003) estimated the fraction of the spectral power of the precession and obliquity bands (excluding 100 ka eccentricity band) in the paleoclimatic records to be less than 15%. When assuming a somewhat broader band spectrum around the astronomical periods, Tziperman et al. (2006) proposed that 65% of the spectral power exists over the full astronomical bands of precession (8%), obliquity (18%) and eccentricity (39%). Meyers et al. (2008) gave yet different estimates: 28% for precession and obliquity bands, 41% for eccentricity band, and 31% for the background continuum. The differences originate from the fact that there is no unique way of separating the spectrum into the two parts. In the latter estimate the spectrum is assumed to be generated by a linear sum of a red noise process and the narrow band astronomical forcing. The relative weights of the two components are obtained by estimating the parameters for the red noise process, an AR(1) process, from the continuous part of the spectrum away from the narrow peak spectrum rather than estimating from the  Mitchell (1976). The dashed lines within the gray area below the continuous spectrum indicate the contributions from processes on regional scales below 1000 km and 100 km respectively. Why variability on small spatial scales should be suppressed on long temporal scales is not well explained by Michell and not clear to us. This is, however, not the focus here, where we shall only discuss the full spectrum lag-one autocorrelation, which is heavily influenced by the narrow peak spectrum (Mann and Lees 1996).
The Hasselmann conceptual framework for interpreting the background also needs to be questioned. Figures 2 and 3 show spectra calculated from the oxygen isotopic ratio of benthic foraminifera, sampled in deep sea cores. The records reflect changes in ice volume and deep-sea temperature, and are considered as reasonable indicators of climate variations during the Cenozoic (66 Ma-present). The spectra of the shortest records, covering the last million years, display a scale break around 100 ka. This scale break at 100 ka has been recognized by several previous studies (Wunsch 2003;Pelletier 2003;Lovejoy and Schertzer 2013;Shao and Ditlevsen 2016) and will hereafter be referred to as the "crossover". As seen by the red curve in Fig. 2, the power spectra of the benthic records are reasonably well accommodated by a Lorentzian (defined in Sect. 3), suggesting that the OU-process with a relaxation time scale of 100 ka is a priori a reasonable model for this part of the spectrum. However, interpreting this in terms of climatic feedback raises the problem: There are no obvious components of the climate system described by linear feedback processes with a relaxation time scale of the order of 100 ka.
The longest records (Fig. 3) have spectra scaling in the low frequency end of the spectrum roughly as P(f ) ∼ f − with ≈ 1 , which is between the two scaling regimes for the OU-process. One approach is to explain a ≈ 1 scaling spectrum from a multitude of OU-processes by extending the stochastic relaxation model to account for multiple relaxation time scales. The spectrum obtained from observations of the system will then be composed of several Lorentzians, which can be arranged to accommodate almost any background spectrum with a slope (in the log-log plot) varying between 0 and − 2, including the scaling spectrum with = 1 (Rypdal and Rypdal 2014). However, this solution exacerbates the problem already mentioned: at time scales of several millions of years the dynamics are better seen as an evolutionary process than as a form of linear relaxation.
Some alternative interpretations of the spectral background were suggested. A constant slope (in a log-log plot) may reveal some form of scaling invariance related to the symmetry of the equations governing the dynamics of the system being observed (Lovejoy 2015). The view is inspired by turbulence theory, in which scale analysis of the governing Navier-Stokes equations provides a constraint on the slope of the expected spectrum (Lovejoy and Schertzer 2013). Similar scale invariance has recently been detected in a low-order ice-sheet climate model, in which it was shown that forcing at millennial time scale can excite low-frequency response following an upward cascade which is constrained by an approximate scale invariance of the dynamical equations (Verbitsky et al. 2018(Verbitsky et al. , 2019. In this view, a crossover, or a change in slope, is therefore related to a change in the nature of dominant processes which control the system evolution over a certain time scale. Lovejoy and Schertzer (2013) used the fluctuation spectrum as an empirical basis to define different such regimes: "weather", "macroweather", "climate", and "macro-climate". The 100-ka crossover in The dating of the former is obtained by fitting to an ice-sheet model forced by the 65N summer insolation, while the latter is dated based on sedimentation rates, independent from astronomical forcing. The spectra are shown in the right panel.
As expected the tuned record (LR04) has slightly higher weights in the astronomical bands than the untuned (HW05) (indicated by grey bars). The discrepancy between the spectra at high frequencies is due to different procedures for smoothing out non-climatic noise in the stacked records. The spectra of the three astronomical parameters are shown at the bottom: eccentricity (magenta), obliquity (purple) and precession (red). The smooth red curve is the Lorentzian indicating a spectral crossover around the astronomical time scales. Here, as well as in the rest of the paper, the spectra are simply estimated by the Fourier periodogram this picture then corresponds to the transition between the "climate" and "macroclimate" regimes. Finally, in the schematic spectrum in Fig. 1 the spectral density varies less than an order of magnitude over the range 10-10 5 years. By contrast, the spectra in Fig. 3 vary by four to five orders of magnitude in the same frequency range.
Our objective here is, on one hand, to reconcile the lack of strong astronomical peaks in the climate spectrum obtained from paleoclimatic observations at glacial time scales with the theories of the dominant role of astronomical forcing in shaping the glacial cycles. On the other hand, we suggest that the spectral crossover around 100 ka should not be interpreted as a relaxation time scale as suggested in the linear stochastic climate model approach.
We shall argue that the interpretation of the spectrum as a superposition of independent discrete (deterministic) and continuous (stochastic) parts, also suggested by the linear stochastic climate model, is insufficient in accounting for the climate spectrum. The observed spectrum strongly suggests that the dynamics implies a broad band response to the narrow band forcing. To approach this problem, we will illustrate the interplay between scale breaks in the continuous spectrum and spectral peaks by focusing on the range of scales imprinting the Pleistocene glacial cycles using the simplest available models of the glacial climate.
After briefly reviewing some basic relations for power spectra and stochastic models, we perform simple numerical experiments with the Saltzman and Maasch (1990) limit cycle glacial model (SM90) to show that even moderate dating uncertainties, realistic for paleoclimatic records, result in substantial spectral broadening of a discrete spectrum. This broadening leads to a spectrum with a scale break similar to what is seen in the estimated spectra for the paleoclimatic observations. For specific parameter choices the SM90 model does show chaos (Mitsui and Aihara 2014), but for illustrating the effect of chaos, we proceed with the classical Lorenz '63 system to show that the interplay between spectral peaks and scale breaks is subtle. The Lorenz '63 model is not relevant as a model of glacial cycles, but it is the most well-known and well studied low-dimensional chaotic system, and the spectral features are similar to any other simple chaotic model. The characteristic time scale of running through one lobe of the Lorenz attractor shows up as a scale break in the x-and y-components and as spectral peaks in the z-component. The different records give a consistent estimate for the spectra. From top to bottom it is seen that for the records covering the past approximately 5 Myrs the spectrum is flat (P(f ) ∼ 1) for time scales longer than the astronomical periods, indicated by the vertical dashed lines. When including the warmer and older parts of the Paleocene, the spectrum is not flat but rather scaling with a slope close to − 1 (P(f ) ∼ f −1 ) , indicated by the black lines.
Here the spectra are calculated with the multitaper method (Thomson 1982), where the bandwidth parameter is 2 and 3 tapers as suggested by Ghil et al. (2002) We then consider some non-intuitive effects of the forcing. Resonance theory tells us that a periodically forced, damped linear oscillator responds with a spectral peak at the forcing frequency. If, however, the forcing also has a stochastic component, a broad peak around the internal frequency appears if the damping is small, and a crossover at a much lower frequency that the internal frequency appears when the damping coefficient is large. To complete the picture, we consider several cases of purely stochastic models with drift constrained by a threshold, all previously proposed as models of the glacial cycles. In these stochastic models the spectrum shows a crossover at a time scale determined by a characteristic transit time.

The autocorrelation and the power spectrum
For later reference, we briefly review the basics of the spectrum and the standard OU-process in the next two sections. We, furthermore derive a few central mathematical results in the Appendix. This can be skipped without loss of continuity, it is included for the presentation to be self-contained. In our discussion of the spectrum so far we have simply obtained the spectrum from the Fourier transform calculated from a time series, without considering the nature of the time series. This empirical approach is natural when considering observations, but in the following we will need to be a little more specific because we compare estimated spectra with analytic results and realizations of different models. Firstly, the finite time series can be considered as "windowed data" from an infinite time series, thus we consider x(t) over −T∕2 < t < T∕2 and assume x(t) = 0 elsewhere. This is a convenient way of handling convergence problems. When relevant we will implicitly assume that the limit T → ∞ is taken. Then, the Fourier transform of x(t) may be written as with the inverse Fourier transformation The power spectral density, here just denoted the power spectrum, is defined in two fundamentally different ways, which are sometimes confused. For a deterministic time series x(t), sometimes called a signal, the power spectrum is (3) while for a stochastic time series, x(t), also called a stochastic process, the power spectrum is defined as where ⟨.⟩ denotes an ensemble mean. Thus the power spectrum is considered as the mean (which is also the expectation value) over many realizations of the stochastic process. An ergodic process is defined by the property that the temporal mean (3) and the ensemble mean (4) (in the limit T → ∞ ) are identical. The autocorrelation function c( ) for the process x(t) is defined as the ensemble mean, Here it is implicitly assumed that the process is stationary, which implies that the statistics is time independent, thus the autocorrelation in (5) does not depend on t. By substituting t → t + it follows that the autocorrelation is a symmetric function, c(− ) = c( ) . Again for an ergodic process we can as well define the autocorrelation as the temporal mean: which is obviously independent of the integration variable t. A very important result, the Wiener-Khinchin theorem, states that the power spectrum and the autocorrelation are a Fourier transformed pair (see Appendix B): Since both c( ) and P( ) are symmetric functions, the integrals are just twice the integrals between 0 and ∞ . The autocorrelation, like any correlation function, is inherently related to stochastic variables. However, for any (square integrable) deterministic function, we can of course formally apply (6). Thus for the harmonic function cos t we get The harmonic function is not a stationary process, but for a rigorous treatment, we can introduce a stochastic phase, , uniformly distributed over [ 0, 2 ∕ ] so that cos (t + ) becomes a stationary stochastic process. By taking the ensemble average over we get a result identical to (8). This trick can easily be applied to any deterministic function. By using the spectral representation of the Dirac delta distribution (see Appendix A), the power spectrum for cos 0 t is obtained as the Fourier transform of (8): which is a single discrete peak in the positive spectrum ( > 0).

Stochastic climate
As argued by Hasselmann (1976), the climate is described by a slow variable x(t) submitted to the effects of relaxation ( − x(t) ) and forced by a fast process, (t) , which is idealized as a standard (unit intensity) white noise process with the correlation function ⟨ (t) (t + )⟩ = ( ) . These assumptions yield the linear stochastic differential equation (SDE): where is the intensity of the noise and the dot represents the derivative with respect to time. This defines the OU-process. An equation for the autocorrelation can be derived and solved directly from (11): where we have used the symmetry, c(− ) = c( ) to obtain the result for < 0 . Physically, the parameter represents stabilising feedbacks pulling the climate state toward the equilibrium state (normalized to x = 0 ). As determines the exponential decay of perturbations, it determines a time scale = 1∕ related to the physical processes governing the stabilising feedback.
The variance is c(0) = ⟨x 2 ⟩ = 2 ∕(2 ) (see Appendix D). The power spectrum is then immediately obtained: This functional form of the power spectrum is called the Lorentzian. For later use, we shall calculate the spectrum for the OU-process in a slightly different but illuminating way by expressing (11) in Fourier space: In a rigorous mathematical sense the Fourier transformed of the white noise, ̂( ) , is not defined, since the white noise is delta-correlated with infinite variance. By multiplying with the complex conjugate, multiplying by 2 ∕T and taking the average on both sides in (13), the correct result is obtained by applying Ito isometry for the white noise (see Appendix E). Thus we obtain the power spectra directly from (4): where the power spectrum of the unit variance white noise is P ( ) = 1∕(2 ) (See Appendix C).
We shall use angular frequency or oscillation frequency whenever convenient ( = 2 f ). A realization of the OU-process is shown in the left panel of Fig. 4. The power spectrum, a Lorentzian (in orange), is shown together with the power spectrum calculated from the realization (in blue) in the right panel of Fig. 4. The power spectrum of this process is flat (white noise) on long time scales; ≪ ⇒ P( ) ∼ 2 ∕(2 2 ) , while on short time scales the spectrum scales as −2 (red noise); ≫ ⇒ P( ) ∼ 2 ∕(2 2 ) . It follows that the crossover frequency, ≈ is the inverse of the characteristic autocorrelation time of the process, which is shown by the arrow in the second panel of Fig. 4.
Since the Fourier transform is linear, it directly follows that the dynamics resulting in the spectrum (14) does not involve any transfer of power between different modes. The dominance of the low frequencies ("reddening of the spectrum"), is merely due to the slower damping at lower frequencies. Had there, say, been an additional periodic component in the forcing, 2A cos̃t on the right hand side in (11), there would be a peak in the spectrum, P ( ) = ( −̃)A 2 ∕( 2 + 2 ) , corresponding to the square of the amplitude of the solution to the simple ordinary differential equation, Equations (11) to (14) do not rely on a hypothesis of scale separation between fast and slow processes in the stochastic climate model. In the case (t) is itself a correlated noise process, with correlation time 1 = 1∕ 1 and variance ⟨ 2 ⟩ = 2 ∕(2 1 ) , the spectrum for x(t) in (11) is readily obtained from (14): This holds even in the case 1 > . There are thus two crossover times and an even steeper −4 spectrum for ≫ max( , 1 ).
A red-noise stochostic forcing (with its own time scale) is one of several mechanisms which may explain several time scales visible as several crossovers on the spectrum. Several crossovers may appear in a linear model with a 2-dimensional (or more) state vector, associated with relaxation equations featuring distinct time scales. This is the framework implicit in the presentation of Mitchell (1976), and which is easily generalised to an arbitrary (up to countable infinite) number of time scales (Fredriksen and Rypdal 2017). Multiple crossover may also occur as a natural nonlinear extension of the linear-relaxation stochastic climate model. In a bistable climate, the climate jumps between two stable states, either because of the external forcing or because of an internal stochastic noise. The latter case can be described as a process similar to Eq. (11) but with a nonlinear drift, where is a double-well potential.
The power spectrum can not be calculated analytically but be estimated as a spectrum composed by the spectrum of the jump process, which is a slow Poisson process of jumping between the two states called a telegraph process (Goldstein 1951), and the spectrum of a fast OU-process close to either of the stable states. The telegraph process also has an exponential autocorrelation, c 1 ( ) = ⟨x 2 ⟩e −� �∕T , where T is the escape time from one well across the potential barrier. This escape time is obtained with Kramer's formula, where U a is the potential at the minimum, U b is the potential at the barrier and H = U b − U a is the height of the barrier. Since the two equilibria are at x = ±1 , the slow telegraph-like process has variance ⟨x 2 ⟩ ≈ 1 . The fast OU-process between jumps can be obtained in the small noise limit, by linear expansion of the potential around either of the equilibria, which with the notation in (16) is exactly the process (11). There is thus a fast correlation time scale −1 . The two Lorentzian power spectra are shown in the second panel in Fig. 5 as red and orange curves, while the sum of the two and the cross spectrum is shown in purple. The blue curve is the spectrum of the realization shown in the first panel. The average slope in the spectrum, which in this case contains two crossovers separated by a plateau, is then somewhere in between − 2 and 0. This is slightly different from the composite spectrum of a sum of independent OU-processes suggested by Mitchell (1976).

The spectrum from low order dynamical models
From these purely stochastic models, we now turn our attention to models with specific identifiable and dominating time scales of variation. Discrete peaks in the spectrum can be related to either periodic or quasiperiodic external forcing or periodic oscillations in the internal dynamics. However, the power in the narrow bands covering the astronomical forcing barely exceeds the continuous part of the paleoclimatic spectrum. A continuous spectrum can arise from either chaotic dynamics or spill out from peaks in imperfectly dated observations. Chaotic dynamics with a characteristic time scale can produce a spectrum with spectral peaks on top of a continuous background. To illustrate this we examine two low-dimensional dynamical models, one glacial cycle model with an internal oscillation and a limit cycle, and one with chaotic dynamics and a strange attractor. The first model is the three variable dynamical model of glacial cycles proposed by Saltzmann and Maasch (SM90): Here x, y and z represent global ice volume, atmospheric CO 2 concentration and ocean temperature in suitable dimensionless form, and parameters are: (p 1 , … , p 5 ) = (0.0075, 0.006, 0.0075, 0.006, 0.009) . Note that in this simple form of the model the glacial cycles are internally generated, completely independent from the astronomical periods.
The x component (ice volume) is shown in the first panel of Fig. 6 (blue). The corresponding spectrum is shown in second panel (blue).
Since the solution is strictly periodic-a limit cyclethe spectrum is discrete, and the first discrete peak corresponds to the period. The rest are over-tones of that. The dating uncertainty in the paleoclimatic records will perturb a strictly periodic signal. We model this uncertainty by assigning the time series x(t i ), i = 1, … , n to n uniformly distributed stochastic time points in interval [ (t 1 , t n ) ], and (18) Ice volume:̇x = − p 1 x − p 1 y Atmospheric CO 2 concentration:̇y = p 2 y − p 3 z + p 4 z 2 − p 1 yz 2 Deep ocean temperature:̇z = − p 5 x − p 5 z then interpolate the time series back to evenly spaced times. With this stochastic uncertainty in dating, realistic for the dating uncertainty in the paleoclimatic records, the signal (black curve in left panel) is no longer periodic and only the main period and the first overtone show up in the spectrum. The spectrum at centennial to millennial time scales is generated by dating uncertainty with a continuous tail not very different from an f −2 spectrum.
The second model is the chaotic Lorenz '63 system; with parameters, ( , , ) = (10, 28, 8∕3) . The attractor is the well-known two-lobe strange attractor (Fig. 7, left panel). The spectrum of the chaotic system is continuous, but the dominant time scale of variation turns up as peaks in the spectrum. The Lorenz system has three unstable fixed points, one at the origin and two at the centers of the lobes (red crosses). The x and y components are essentially identical and monitor in which lobe the system is while the z component contains no information on which lobe the system is in. Realizations of the x(t) (blue) and z(t) (green) components are shown in Fig. 7, middle panel. The corresponding power spectra are shown in the right panel. The spectrum for x(t) (blue) only shows a crossover from a fast decaying high frequency part to a flat spectrum above a crossover time scale. In the spectrum for z(t) the crossover time scale shows up as a strong discrete peak. The gray sine curve on top of the z(t) curve (green) in the middle panel corresponds exactly to the peak, thus it is seen that this is simply the period of rotating around the center of one lobe. The gray  (16). The power spectrum (right panel) can be understood as the sum of two Lorentzian spectra, and the cross spectrum, from the jump process (a random telegraph process) and the OU-process in either of the steady states. The Lorentzians of a telegraph process and of a fast OU-process are shown in orange and red respectively. The sum of the two and the cross spectrum between the two is shown in purple sine curve on top of the x(t) curve (blue) is obtained from the imaginary part ( 0 ) of the eigenvalue obtained from a linear perturbation around the central unstable fixed point. The two gray sine curves have slightly different frequencies.
Though x(t) clearly oscillates at the frequency obtained from the stability analysis, this frequency does not show up in the corresponding power spectrum. The reason is that the projection of x(t) onto sin 0 t from one lobe cancels the projection from the other lobe, leaving only the continuous part related to the chaotic shifts between the two lobes. In fact, a spectrum estimate from a time series too short for allowing inter-lobe transitions will actually show a peak (Fig. 8).
Indeed, when there is no inter-lobe transition, no phase cancellation occurs. Another way to describe the phenomenon is to observe that the z-component, which shows a peak in the spectrum, is degenerate in the sense that a Takens' delay embedding-a plot of [z(t), z(t + ), z(t + 2 )]-will not for any value of reproduce the two-lobe structure of the original attractor.
These simple examples illustrate firstly, that the discrete spectrum of a limit cycle can be smoothed into a continuous spectrum by dating uncertainty and secondly, that for autonomous chaotic systems the power spectrum contains both a continuous part and spectral peaks. Both cases evidence difficulties for separating the discrete and continuous part in the estimated power spectrum. Dating uncertainties will result in an underestimation of the discrete part, while the Lorenz model case points to an entirely different effect: The non-degenerate x and y coordinates in the Lorenz model

Crossover and resonance
We now turn the attention to the crossover in the spectral slope in the paleoclimatic records around the time scales of the astronomical forcing. In the stochastic climate model perspective, a crossover in the background spectrum is interpreted as the typical relaxation time scale in the climate system. This is, however, an unlikely explanation for the crossover observed at 100 ka, since there is no relaxation process known to operate at that time scale. Over the last million years, glacial cycles have been characterized by an oscillation around 100 ka. The question is whether this oscillation is causing the crossover in the spectrum. To this end, we represent the climate signal by a linear damped oscillator, forced by the sum of white noise and a single harmonic function: This stochastic differential equation is easily solved in Fourier space, in analogy to the solution of the OU-process above: where P F ( ) is the power spectrum of the forcing F(t). The spectrum for the forcing is readily obtained from this autocorrelation: c F ( ) = 2A 2 cos f + 2 ( ) , which is obtained (21) using (8) and the independence between the two components. Thus, , and we have Due to the combined deterministic and stochastic forcings this is a composite spectrum, with a spectral peak at the forcing frequency f , depending on the strength of the damping and how far the forcing frequency f is from the resonance 0 . Furthermore, there is a continuous background proportional to the variance of the noise 2 and a broad peak around the resonance.
Two realizations of this process are shown in Fig. 9. In both cases the forcing frequency is f = 2 ∕41 kyr −1 corresponding to the obliquity cycle and the internal frequency is taken to be 0 = 2 ∕10 kyr −1 . The spectra are shown in the right panels with the orange curves obtained from (22). The continuous spectrum depends substantially on the strength of the damping: For small damping ( = 0.5 kyr −1 ) the spectrum has a resonance peak around 0 . This is indicated by the arrow in the top right panel. It has a steep P( ) ∼ −4 spectrum for shorter time scales. In the other case of a strongly damped oscillator ( = 10 kyr −1 ), the crossover moves to a much longer time scale ∼ ∕ 2 0 , indicated by the arrow in the lower right panel.
Other forms of resonance have been discussed in the context of Quaternary dynamics. The one which is arguably the most physically relevant in this context is coherence resonance, proposed by Pelletier (2003). Stochastic resonance was historically suggested as a possible mechanism for glacial-interglacial cycles (Benzi et al. 1982). We will discuss it briefly as well, even though its physical motivation for ice ages is probably no longer tenable. Coherence resonance, in an excitable system, is a phenomenon where stochastic noise drives the system away from a stable fixed point into a large excursion in phase space returning back into the vicinity of the fixed point. Two time scales are at play here: the Kramer waiting time for escaping from the fixed point and the excursion time. For small noise, excursions are rarely excited and the waiting time becomes irregular, while for large noise, the trajectory is disturbed and the excursion time becomes irregular. For a certain amount of noise, the stochastic system becomes more coherent (more periodic). Pikovsky and Kurths (1997) described the coherence resonance phenomenon in the noisy FitzHugh-Nagumo model. They suggested the following form of the autocorrelation function: where is a decay rate and 0 is the resonance frequency. The corresponding power spectrum is given by the Fourier transform of C(t) by the Wiener-Khinchin theorem: The power spectrum P( ) is shown in the left panel in Fig. 10 with different parameter values. As before there are two scaling regimes, for ≫ 0 + , the power spectrum scales as P( ) ∼ −2 , while for ≪ 0 + , it is constant such that P( ) = ∕[ ( 2 + 2 0 )] . The crossover frequency = 0 + is related to the peak of the spectrum for strong resonance: ≪ 0 ⇒ ≈ 0 . Following Pikovsky and Kurths (1997), we integrate a stochastic FitzHugh-Nagumo model: where a broad resonance is seen at the internal frequency. Bottom row is for large damping, = 10 kyr −1 , where the spectral crossover is approximately at 2 0 ∕ , indicated by the yellow arrow. The orange curves are given by (22) with the parameters, ( , a, ) = (0.01, 1.05, 0.25) . We scale time, such that the resonance frequency 0 is 2 ∕100 kyr −1 . The y component is shown in the middle panel in Fig. 10, with the corresponding power spectrum in the right panel. This spectrum is well fitted by (24), which is multiplied by the variance of y in order to satisfy Parseval's theorem. Benzi et al. (1982) introduced the stochastic resonance as a possible explanation of the 100 ka spectral peak. The original premise is that changes in globally averaged, annual mean insolation, which is controlled by excentricity only (here assumed to be represented by a pure 100 ka oscillation) are too weak to cause glacial-interglacial transitions. An amplifying mechanism is thus needed, this is achieved by modelling the global temperature x as a mass-less (no acceleration term) viscous particle moving in a double-well potential, with stochastic forcing, and with a weak periodic forcing with angular velocity f : The double-well potential, V(x), is such that with no forcing the system has two steady states x = ±1 . Furthermore, the amplitude of the periodic forcing is smaller than needed to destabilize either of the steady states (A < 2 √ 3∕9) . The forcing term can be included in the potential by redefining V(x, t) = x 4 ∕4 − x 2 ∕2 − Ax cos( f t) . This potential has one deep and one shallow well when the forcing is either maximum or minimum, determining two distinct Kramers rates r s→d and r d→s for noise assisted escapes from the shallow and from the deep well, respectively. Resonance happens when the intensity of the noise is such that r d→s ≪ f ≪ r s→d .
An approximate expression for the power spectrum of the stochastic resonance phenomenon in the periodically-forced double-well potential system is (Gammaitoni et al. 1998): where f is the frequency of the external forcing and r K is the Kramers rate of the noise-induced hopping in the symmetric double-well potential. The spectrum consists of a discrete peak and a Lorentzian part analogous to the case of the linear damped and forced oscillator. Note that the resonance frequency = f is unrelated to the crossover frequency 2r K of the Lorentzian part. Here, two times the escape rate 2r K plays the role of the internal frequency 0 (Fig. 11).
In summary, we have shown that for the resonance oscillator models the crossover frequency is not directly linked to the external forcing frequency. For the linear damped oscillator the crossover time scale depends on the linear damping and the internal frequency: ∼ ∕ 2 0 , while in the coherence resonance case it is inversely related to the sum of the resonance frequency and a damping rate: ∼ 1∕( 0 + ) . Finally, in the case of the stochastic resonance, the crossover time is related to the Kramers rate of noise-induced transition, ∼ 1∕r K .

Threshold models
A different perspective, focusing on the 100-ka late Pleistocene time scale of glaciation, was presented by Wunsch (2003). He suggested that the local maximum in 100-ka variability in the power spectrum should not be interpreted as a resonance peak. It is "a reflection of the absence of energy at the adjacent lower frequencies". In this interpretation, the existence of a crossover is not related to a feedback time scale in the system. The hypothesis is that the ∼ 100 ka glaciation is rather governed by constraints in ice volume independent from the astronomical forcing. Our first model encoding this assumption is that the ice growth is a purely random walk process, ̇x = bounded by a maximum volume x max . The constraint can be modeled as reflecting boundaries for the random walk, making x a stationary process. A characteristic time scale in this case would be the diffusion time for diffusing the distance x max : ∼ (x max ∕ ) 2 . Such a process is shown in Fig. 12 top left panel. The power spectrum is shown in blue in the top right panel together with a power spectrum (orange) obtained from a 10 5 ka realization. The diffusion time is indicated by the arrow in the top right panel.
A slightly more realistic threshold model was proposed by Huybers and Wunsch (2004) (HW04). This is a random walk with a constant drift ̇x(t) = + (t) . The ice volume x(t) collapses to 0 when reaching a threshold x max , corresponding to glacial terminations. This process is shown in Fig. 12 middle left panel. Though a distinct time scale for glacial cycles can be derived from the linear drift, = x max ∕ , the spectrum has a weak spectral peak coinciding with the crossover, indicated by the arrow, as is seen in middle right panel.
With no noise, the purely deterministic version of the HW04 model is simply a periodic sawtooth curve. It has a simple Fourier series, with f = ∕x max . Thus its discrete power spectrum S(f n ), (f n = nf ) is given as which, by coincidence, is the same spectral slope as we have for the OU-process. This discrete spectrum is indicated by the grey bars in the middle right panel of Fig. 12.
In a slightly more advanced version of the HW04 model, the astronomical forcing is introduced via a time varying threshold, x max (t) = x 0 max + F(t) where F(t) is a suitably rescaled version of the forcing. In this case (lower panel in Fig. 12), using the obliquity curve (Berger 1978), the response is an approximate 2:1 frequency locking to the 41 ka almost harmonic obliquity curve. The 2 × 41 ka = 82 ka period shows up as broad peak in the spectrum, while the 41 ka period shows up as a more narrow peak (Fig. 12,lower right panel).
To summarize, for the threshold models the crossover time scale is either set by a diffusion time scale, or given by a linear drift, which in the low noise limit corresponds to the time needed for establishing a glacial maximum. When assuming that the astronomical forcing governs the threshold, the spectrum shows peaks at multiples of the forcing frequencies in accordance with a frequency locking scenario (Nyman and Ditlevsen 2019).

Summary and conclusions
Spectrum estimates of paleoclimate records tend to show weak astronomical modes over a dominant continuous stochastic background. While this has been noted since the  Fig. 12 Three models, where left panels are the realizations, right panels the power spectra of the realization (blue) and the mean spectra are obtained from 1000 realizations of the model (orange): Top row is a random walk with reflecting boundaries at x max = ± 0.9 . The yellow arrow in the right panel corresponds to the diffusion time for moving between x = 0 and x = x max . The middle row is a linear drift with threshold x max = 0.9 , where x is reset to 0. This is the Huybers-Wunsch model of glacial terminations. The broad peak indicated by the yellow arrow corresponds to the mean drift time to reach x = x max after a resetting. The bottom row shows is the Huybers-Wunsch model with the time varying threshold determined by the obliquity parameter, which is close to a 41 ka harmonic curve. With the linear drift being the same as above, the system shows an approximate 2:1 frequency locking to the period of the forcing. This is reflected in the power spectrum, where the broad peak is shifted towards 1/82 ka −1 , indicated by the yellow arrow. The peak at 1/41 ka −1 is visible in the mean spectrum seventies, it still challenges the widespread conception that glacial cycles are generated by changes in Earth's orbit and inclination. Any satisfactory theory of Quaternary variability must explain both the weakness of astronomical modes in the power spectrum, and the existence of a crossover around 100 ka, yet be compatible with the fact that deglaciations and glacial inceptions are paced by the orbital forcing. On the basis of the different simple dynamical and stochastic models analysed in this study, we can list possible keys to the interpretation of the crossover time scale. Most deterministic models supporting ice age theories present one strong mode around the 100 ka period with little or no background power. We used the SM90 as an archetypical representative of this class of models. The model has a limit cycle solution for the ice ages with a dominant mode. However, realistic dating uncertainty in proxy records weakens such modes, generating a scaling background spectrum with appealing properties; a ramp up to 100 ka followed by a scale break. Hence, part of the dispersion from the mode to a background continuum may conceivably originate from variability in the sedimentary deposition rate and other dating uncertainties.
In contrast to oscillator models, chaotic models have a continuous spectral background. This is illustrated by the Lorenz '63 model, which also for the z-component has spectral modes over the continuous background. The crossover to a flat spectrum occurs as it should be expected at the correlation (Lyapunov) time scale set by the chaotic character of the dynamics. In this case the spectral peaks and the continuous background are parts of the same dynamics, thus estimating weight of a peak above a background in the spectrum of a specific observable is not relevant. In this simple case the high frequency end of the spectrum falls off exponentially, reflecting the smooth dynamics, which of course does not provide a realistic representation of the multi-scale climate dynamics. In the case of the Lorenz '63 model, a spectral peak corresponding to the oscillatory mode around the lobes does not show up in the spectrum of variable x, thus for any given dynamical model an analysis and cautious interpretation of the power spectrum is required.
At the other extreme, in purely linear stochastic climate models, fluctuations of the stochastic components will generate a Lorentzian spectrum. The basic problem is the physical interpretation that a crossover near 100 ka would require a relaxation time scale of that order of magnitude, which is difficult to identify in known components of the climate system. What we show here is that this is a too simplified view. We examined alternative simple linear and non-linear processes, generating a crossover which need not to be directly related to a time scale of a process. This means that on the one hand, one should not rule out certain models solely based on the observed power spectrum. On the other hand estimates of crossovers and spectral peaks could perhaps be interpreted in such a way as to identify physical mechanisms. Unexpectedly, a system with linear resonance can generate a crossover at a much lower frequency than the resonance frequency if it has a very strong damping (crossover frequency at 2 0 ∕ , cf. Fig. 9). In coherent resonance, the position of the crossover time scale relative to the resonance period is related to a decorrelation time scale, which itself depends on the noise level. In stochastic resonance, the crossover time scale is related to the typical exit time from one potential well to the other.
Finally, models associated with threshold-triggered transitions also generate crossovers, the time scale of which is associated to the typical transit time-physically interpreted as the time needed to build an ice sheet until it reaches levels of instability triggering collapse. Incorporating the astronomical forcing through a varying threshold can generate a spectrum with spectral peaks similar to the spectra generated from the paleoclimatic records.
To conclude, in the spirit of Mitchell's original paper, we used the simplest possible models of the ice age climate, concentrating on the millennial to astronomical time scales in the Pleistocene climate spectrum. We have shown that the interplay between the astronomical peaks and the crossover observed in the climate records can have several origins. As the time series of the past climate do not permit discrimination between the various proposed models, these obviously do not constitute a final theory of the Pleistocene climate. In our view, the different stochastic, deterministic, self-oscillating, and forced models describe different aspects of the climate system, which need not be mutually exclusive.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.