The Middle Pleistocene Transition by frequency locking and slow ramping of internal period

The increase in glacial cycle length from approximately $41$ to on average $100$ thousand years around $1$ million years ago, called the Middle Pleistocene Transition (MPT), lacks a conclusive explanation. We describe a dynamical mechanism which we call Ramping with Frequency Locking (RFL), that explains the transition by an interaction between the internal period of a self-sustained oscillator and forcing that contains periodic components. This mechanism naturally explains the abrupt increase in cycle length from approximately $40$ to $80$ thousand years observed in proxy data, unlike some previously proposed mechanisms for the MPT. A rapid increase in durations can be produced by a rapid change in an external parameter, but this assumes rather than explains the abruptness. In contrast, models relying on frequency locking can produce a rapid change in durations assuming only a slow change in an external parameter. We propose a scheme for detecting RFL in complex, computationally expensive models, and motivate the search for climate variables that can gradually increase the internal period of the glacial cycles.

epochs of icy and cold conditions on the one hand, and warm and ice-free conditions on the other (Fig. 1). While these glacial cycles are attested from geological records (EPICA community members, 2004;Huybers, 2007;Lisiecki and Raymo, 2005), there is no single conclusive theory of their origins. Historically, the focus has been to explain the approximately 100 thousand years (kyr) long glacial cycles that dominate the past 800 kyr (see (Imbrie and Imbrie, 1979) for a review). But as sediment records later revealed that these cycles were ca 40 kyr long prior to 1.2 Myr ago the question arose what caused the shift to approximately 100 kyr long cycles, called the Middle Pleistocene Transition (MPT) (Clark et al., 2006).
The main strategy to address these questions has been to replicate the palaeoclimatic records using simple models of glacial cycles with few variables, referred to as conceptual models (see (Crucifix, 2012) for a review). One reason for this is that the rather regular and cyclic variations in data suggest that the the main dynamics can be captured by a system of few degrees of freedom, even as the full climate system obviously has a large number of degrees of freedom. None of these models describe the climate system in detail, but they are useful for understanding underlying dynamics. Virtually all models involve insolation variations due to changes to Earth's orbital configuration relative to the sun, an idea heralded by Adhémar, Croll and Milankovitch (Imbrie and Imbrie, 1979). But the specific role played by insolation variations is still unknown and debated.
Several solutions to the cause of the MPT have been presented within the context of conceptual models. Some mechanisms rely on a bifurcation occurring in the unforced climate system which fundamentally changes how the system operates (Ashwin and Ditlevsen, 2015;Ditlevsen, 2009;Huybers and Langmuir, 2017;Maasch and Salzman, 1990;Tziperman and Gildor, 2003). Other mechanisms invoke a "spontaneous" change, such as a shift between attractors due to subtle changes in insolation (Omta et al., 2015 (Lisiecki and Raymo, 2005). Higher values means more ice. b) Durations between successive major terminations (black), showing an increase at the Middle Pleistocene Transition (MPT). Error bars indicate one standard deviation dating uncertainty. Contours (darker to lighter) show the amplitude of a wavelet spectral estimate. Each lighter contour corresponds to an increase of 7% of the maximum amplitude, starting at 30%. Outside the cone of influence (thick black line), edge effects are important. See Appendix B and C for details et al., 2018) or random fluctuations (Imbrie et al., 2011;Salzman and Verbitsky, 1993), or as a coincidence (Huybers, 2009). A third possible mechanism for the MPT assumes one essential mode of oscillation throughout the Pleistocene and relies crucially on the interaction between insolation variations and an increasing internal period. This mechanism, previously imprecisely referred to as phase/frequency locking and non-linear resonance -but here Ramping with Frequency Locking (RFL) -is the focus of this paper.
The main appeal of this mechanism is that nothing special had to occur in the climate system over the MPT (Huybers, 2007); it is only required that the internal period was ramped slowly -interactions with forcing are enough to cause an abrupt increase in durations between glacial terminations (Fig. 1, bottom panel).
The first publications where RFL was used (Paillard, 1998;Paillard and Parrenin, 2004) did not explain why the durations between glacial transitions increased abruptly over the MPT. Ashkenazy (2006); Ashkenazy and Tziperman (2004) hinted how frequency locking (therein called phase locking) could produce an abrupt increase in duration, by showing diagrams of average duration as a function of a system parameter (Devil's staircases). Huybers (2007) was first to both show a model trajectory of ice volume using the mechanism, and to attribute the effect to "skipping of obliquity cycles", a frequency locking effect. Recently, Daruka and Ditlevsen (2015); Feng and Bailer-Jones (2015); Mitsui et al. (2015); Tzedakis et al. (2017) alluded to the mech-anism, but neither emphasised that frequency locking can explain the MPT assuming only a slow linear change in a climate parameter. Instead, by ramping some parameter in a way that mimics the rapid change in durations over the MPT, they prescribe an abrupt increase in period over the MPT rather than explaining it. Here, we for the first time properly define RFL and emphasise its generality.
Following (Huybers, 2007), we question the common assumption that climate entered a stationary state in the late Pleistocene, and instead argue that the sequence of durations between glacial terminations is consistent with a slow increase of the internal period of the climate system until present (Fig. 1 b)). According to this view, the typical durations between transitions changed from ∼40 kyr to ∼80 kyr around 1200 kyr ago, after which they increased gradually in the mean to present time, with the last duration being ∼120 kyr long.
Here, we first aim to explain RFL in a clear way, using a harmonically forced simple model. We use harmonic (pure sine) forcing because it makes frequency locking concepts clearer, while still producing qualitatively similar behaviour to astronomical forcing curves. We should not expect model runs with such simplified forcing to agree well with data, however. We use forcing with period 41 kyr, corresponding to the main period of obliquity variations (Berger, 1978), which determine the total insolation integrated over the summer at Northern latitudes (Huybers, 2006). We then define RFL, specify a class of models able to reproduce the MPT using the mechanism, and propose a decomposition of model components to understand the abruptness of the MPT. We consider evidence in data for a 40 to 80 kyr shift in durations between terminations and a subsequent gradual increase, and why this supports RFL in favour of some other mechanisms for the MPT. We then discuss how insights from harmonic forcing relate to nonharmonic forcing, how RFL can be detected in complex and computationally expensive models, and some climate variables that can cause an increase in the internal period of the glacial cycles.

The idea behind Ramping with Frequency Locking
We illustrate Ramping with Frequency Locking (RFL) using a deterministic and continuous time version of the H07 model (Huybers, 2007) (see Fig. 2). The model is arguably the simplest to represent alternating stages of intrinsic growth and decay of ice sheets, with the growth state ending abruptly as a critical ice volume is reached. It is is an integrate-and-fire threshold model conceptually very similar to the models in (Ashkenazy and Tziperman, 2004;Glass and Mackey, 1979;Huybers, 2007;Imbrie and Imbrie, 1980;Imbrie et al., 2011;Paillard, 1998;Paillard, 2003, 2012;van der Pol and van der Mark, 1927;de Saedeleer et al., 2013).
We assume that ice volume x(t) grows at a constant rate µ in a glacial state until it reaches a threshold θ(t). Then deglaciation starts, whereby ice volume decays to 0 over a fixed time T decay = 10 kyr: linearly decrease x(t) to 0 over time T decay , repeat. (1) Small perturbations to the model, such as having a constant rate of decay instead of a fixed time, does not qualitatively affect its behaviour. We split θ(t) into a forcing term A · F(t) -a zero-mean sum of periodic components -and a ramping term R(t): In the limit of constant ramping R(t) = R 0 and zero forcing A = 0, the system has a constant internal period of oscillation T o = R 0 µ + T decay (subscript o for oscillator), see Fig. 2 a). But if the threshold increases slowly over time, for instance linearly θ(t) = R(t) = R 0 + R 1 t as in Fig. 3 a)i), then the internal period T o (t) = R 0 µ + T decay + R 1 µ t also increases slowly (Fig. 3 b)i)). As there is no forcing, the durations between glacial terminations follow T o (t) closely. (1)), a) unforced and b) periodically forced. Ice volume (blue) grows linearly at a rate µ until a threshold (red) is hit, after which ice volume is reset to 0 over a time T d e c a y . In a) the threshold of glacial termination θ(t) is constant θ(t) = R 0 , whereas in b) it oscillates periodically as θ(t) = R 0 + A sin (2πt/T f ). T o is the internal (unforced) period of the model However, with periodic forcing θ(t) = R 0 + R 1 t + A sin (2π/T f ), durations D i are near multiples of the forcing period D i ≈ NT f , N ∈ N (Fig. 3 b)i)). Roughly speaking, the multiple that is realised is the one closest to the internal period T o . This phenomenon, called frequency locking (Pikovsky et al., 2001), has been studied extensively over the past century (e.g. Cartwright and Littlewood (1945); Glass and Mackey (1979); Le Treut and Ghil (1983); van der Pol and van der Mark (1927); Tziperman et al. (2006)).
In Fig 3 a)i) the durations D i change abruptly from 1 ×T f to 2×T f and finally 3×T f . These abrupt changes in durations resulting from a gradual change in an underlying parameter is one possible dynamical mechanism behind the MPT.
We call the mechanism Ramping with Frequency Locking (RFL), rather than non-linear resonance, phase locking or frequency locking as it has previously been called. This we do to emphasise both that an internal period must increase gradually over time (ramping), and that the internal oscillations must be locked to external forcing. This is opposed to e.g. the mechanism in (Omta et al., 2015), which realises the MPT through jumps between coexisting frequency locked solutions.
We note that RFL is a special case of "slow passage through bifurcation" (e.g. (Baer et al., 1989;Do and Lopez, 2012)), for which the bifurcations typically are saddle-node bifurcations of limit cycles marking transitions in and out of frequency locking regions (Pikovsky et al., 2001). (However, see e.g. Guckenheimer et al. (2003); Levi (1990) for other relevant bifurcations).
Finally, we note that H07 is an illustrative example of RFL, and not representative of all glacial cycle models. However, the rapid jumps between frequency locking regions occur generically in a broad class of models, defined next. Fig. 2) is just one particular model capable of realising the MPT through RFL. We could simply call these self-

H07 (
The Ramping with Frequency Locking (RFL) mechanism for the periodically forced H07 model. a) Ice volume (blue sawtooth) over time for a i) linear and ii) sigmoidal ramp of the upper threshold θ(t). Periodic threshold in orange and unforced threshold in grey. b) Average duration between glacial terminations D τ over frozen time τ (black solid) (known as Devil's staircases, Sect. 4.1) for the i) linearly and ii) sigmoidally ramped thresholds, and sample durations (red dotted lines) for the forced solutions in a). Magenta lines show the internal period T o (τ). The inset shows that durations with time-varying ramping R(t) do not agree perfectly with the average duration D τ , computed for R(t) = const (see Sect. 5). Parameters are µ = 1, T d e c a y = 10, T f = 41, A = 20 and ramp functions are i) R(τ) = 26 + 0.05 × (τ + 2000) and ii) R(t) = 26 + 50(tanh ( t +1000 300 ) + 1) respectively sustained oscillators, but we aim to be more precise and to establish notation.
First, we naturally require the model to be a dynamical system, such that there is an evolution rule f (t, x) taking a state x(t) forward in time t. We identify the model with the evolution rule and denote it f (without arguments) for brevity. f (t, x), x(t) and t can be very general, for instance; t can be continuous or discrete, x can be of any dimension, and f (t, x) can e.g. be a piecewise smooth ODE paired with a switching rule, as for H07.
We also require f to be forced by a continuous zeromean sum of periodic components A(t)F(t) with an amplitude A(t), called the forcing. We further require that f is parametrised by a set of parameters p(t), whose timevarying subset R(t) is called the ramping. Thus we can write f = f (t, x, R(t), A(t)F(t)).
We define the frozen system f τ := f (t, x, R(τ), A(τ)F(t)) as f with parameters frozen at time t = τ. Importantly, we require that f is a self-sustained oscillator with internal period T o (R(τ)), meaning that every solution to f τ with A(t) ≡ 0 tends asymptotically (as t → ∞) to a periodic solution with period T o (R(τ)). For RFL to be relevant we require that T o (R(τ)) increases as a function of τ. This is the ramping part of RFL.
The frequency locking part of RFL comes from the response of f to non-zero but constant forcing A(τ). For small and medium size A(τ), asymptotic solutions to the frozen system f (t, x, R(τ), A(τ)F(t)), are generically periodic with periods related rationally to the forcing periods (Pikovsky et al., 2001). The oscillator period can for instance be twice that of the forcing period. If so, the oscillator period (and therefore frequency) remains constant on open sets of parameters and we say that solutions are frequency locked to the forcing (we return to this in Section 4).
The essence of RFL is that the period of the frozen system can change rapidly as function of T o (R) when a ramped parameter causes the system to switch between frequency locking regions.
However, some remarks are in place. Firstly, the system with time varying parameters f (t, x(t), R(t), A(t)F(t)) is not the same as the frozen system f (t, x(t), R(t), A(τ)F(τ)) since solutions to the former cannot equilibrate to solutions of the latter in finite time. We return to differences between the two systems in Section 5 but until then we focus on the frozen system.
Secondly, the period of an oscillator is not the same as the length of individual "cycles". For instance, around −1350 kyr in Fig. 3, short and long "cycles" alternate. This makes the average time between terminations 61.5 kyr, whereas the period (time until repetition, two large peaks) is 123 kyr. Therefore, we instead characterise local behaviour with the average duration where D i,τ denotes the i:th duration between successive crossings of a fixed threshold for the frozen system f τ . For some models like H07 glacial terminations are natural such thresholds. For other models Poincaré sections can be considered instead (Pikovsky et al., 2001).

Breaking down the dependency of D τ on τ
Comparing Fig. 3 b)i) and b)ii) shows that D τ can rise steeply both from frequency locking effects under a gradual change of parameter (Fig. 3 b)i)) and from ramping of a climate parameter rapidly (Fig. 3 b)ii)).
We wish to break down the contribution to the local change in average duration from these effects and do so by considering the change ∆D τ under a small perturbation ∆τ: where e.g.
∆T o , and where we have neglected higher order terms. This approximation is generally better the smaller ∆τ is. As ∆τ → 0, (3) tends to the chain rule, but since ∆D τ (T o , A) ∆T o = 0 wherever differentiable (see Section 4.1) it is more appropriate to consider ∆D τ over short intervals of time ∆τ. In what follows we restrict ourselves to ∆A(τ) ∆τ = 0. Eq. (3) says that the rate of change (abruptness) in time of the average duration D τ is approximately the product of the rates at which R(t) changes with time, T o changes with R, and D τ changes with T o . Our point is that each of D τ , T o and R can contribute to an abrupt change of D τ at the MPT, but they  (Crucifix, 2013;Pikovsky et al., 2001;de Saedeleer et al., 2013). Inside major 1:N tongues, solutions are periodic with period N times the forcing period T f = 41 kyr, as evidenced in (Fig. 4 a)). Minor tongues emanate at A = 0 from other rationals of T f , and in between them are quasiperiodic solutions. (We show only M : N, M = {1, 2} Arnold tongues, defined numerically as sets for which |D τ − T f N M | < 0.5. D τ is estimated over 6 million years.) A change in τ that in turn leads to a change in T o (R(τ)) traces out a path in (A, T o ) space (black and magenta lines in Fig. 4 a). Such a path represents the change in system state as one or more parameters change in time over the MPT. The paths in Fig. 4 a) pass through the major 1:1, 1:2 and 1:3 locking tongues, in which there are respectively 1, 2 and 3 forcing periods per oscillator period. We learn that for larger A, a larger portion of the path stays inside the major 1:N tongues, an observation also made in (Ashkenazy, 2006).
Another way of visualising the change in average duration D τ as a function of T o are Devil's staircases (Pikovsky et al., 2001) (Fig. 4 b)), in which the forcing amplitude A is fixed. We see that the average duration D τ is constant within Arnold tongues and that the staircase for larger A contains longer steps of constant duration, as predicted from Fig. 4 a).
. For further details, see Fig. 4 and the text Hence, stronger forcing influence tends to cause more abrupt changes to the average period.

T o (R) and R(τ)
The function T o (R), if continuous and monotonic, stretches and squeezes Arnold tongues by scaling the independent variable T o of D τ (T o ). In the Ashkenazy model (Ashkenazy, 2006), for instance, a faster-than-linearly increasing T o (R) makes the 1:2 and 1:3 Arnold tongues, as a function of ice volume threshold, narrower and more closely spaced than the 1:1 tongue.
Ramping R(τ) continuously and monotonically, similarly stretches and squeezes Arnold tongues. For instance, in Fig. 5 a sigmoidal ramping R(t) makes the 1:2 Arnold tongue narrower compared to a linear change of R(t). Fig. 3 further illustrates this, showing model runs for either a sigmoidally (R(t) = 26 + 50 × (tanh( t+1000 300 ) + 1)) or a linearly R(t) = 26 + 0.05 × (t + 2000) ramped threshold. The sigmoidal ramping accelerates the increase in average duration around −1000 kyr, making the transition more abrupt. Note that the parameters in the functions R(τ) in Fig. 5  A model having an abrupt change due to D τ (T o , A) relies on frequency locking properties, and assumes only slowly varying functions T o (R) and R(τ). Hence, the internal period is assumed to change slowly with model parameters, and parameters are assumed to change slowly in time. Such a model, relying on few assumptions about the climate system, makes full use of the RFL mechanism. The models in (Huybers, 2007;Paillard, 1998;Paillard and Parrenin, 2004) and H07 are of this kind.
A model which relies predominantly on T o (R) for an abrupt change in average duration D τ is also consistent with a slowly changing external parameter R(τ), but a particular function T o (R) requires a physical explanation.
A model relying on a rapidly changed external parameter R(τ) does not need frequency locking properties of D τ (T o , A) or a non-linear response of internal dynamics to the parameter T o (R). However, such a model prescribes the abrupt change in average duration at the MPT rather than explaining the dynamics behind it. Therefore, such an explanation requires justification for the rapidly changed external parameter. The models in (Ashkenazy and Tziperman, 2004;Daruka and Ditlevsen, 2015;Mitsui et al., 2015;Tzedakis et al., 2017) can be said to fall under this category, although they also achieve some abruptness through D τ (T o , A).

Validity of the quasi-static approximation f ∼ f τ
The quasistatic approximation is the approximation that parameters R(τ) change so slowly that the local average duration of f at time t = τ is equal to D τ . I(τ) is the set of indices of durations D i,τ within a time interval [τ−τ 0 , τ+τ 1 ] around τ, with τ 0 , τ 1 > 0. If I(τ) = ∅, then we define D τ,loc = 0. If the quasistatic approximation holds, then the average duration, Arnold tongue diagrams and Devil's staircases calculated for the frozen system f τ provide accurate information about local dynamics of f . However, if R(t) and/or T o (R(t)) change rapidly around t = τ, then there are two sources of discrepancy between D τ,loc and D τ .
The first comes from that the length of the interval of time needed for a good average may be long relative to the local change of D τ for the system f τ . Fig. 6 illustrates that for a 9 × 41 = 369 kyr-periodic solution, a long interval is needed to get a local average duration D τ,loc in agreement with D τ . At the same time, a long averaging interval fails to capture abrupt changes to D τ .
The second is that solutions to f may fail to track solutions to f τ . This occurs if the "frozen" attractor of f τ changes (in some sense) at a fast rate, and if solutions attract to the frozen attractor at a slow rate. Quantifying these rates in a coordinate-and model-independent way seems difficult, however. A candidate measure of rate of attraction is the maximal Lyapunov exponent of the return map mapping one transition time to another (Pikovsky et al., 2001). This can be normalised to a common time scale between models and is coordinate independent. However, since it is only a local measure it neglects the time it takes to enter a small neighbourhood of the attractor. This time can in practice dominate, as is the case in the standard circle model (not shown, model described in (Pikovsky et al., 2001)).
The local change in average duration | ∆D τ ∆τ | is a candidate measure of rate of change of an attractor of f τ , since it exists in all models f and is coordinate-independent. It is ambiguous how large ∆τ should be, however. Furthermore, the average duration D τ is only a proxy for the position of an attractor in phase space; an attractor can move even if ∆D τ ∆τ = 0. This explains the consistent deviation of single durations from the predicted and locally constant D τ = 82kyr in the inset of Fig. 3 b) i).

Is there a 100 kyr world?
The late Pleistocene (∼800 -0 kyr) is sometimes referred to as the "100 kyr world", carrying the implicit notion that the Earth system has settled in a stationary mode with a dominant time scale of 100 kyr (Fig. 1 a)). This view, originating from the closeness to the 100 kyr component of eccentri-city (an astronomical parameter), is supported by the rate of increase of mean ice volume seemingly levelling off (Clark et al., 2006;Mudelsee and Schulz, 1997), and that the Fourier spectrum over the last ∼800 kyr is centred around 100 kyr.
We propose on the contrary, following (Huybers, 2007), that the glacial period increased gradually from ∼80 kyr around −1200 kyr to ∼120 kyr at present day. The change from ∼40 to ∼80 kyr long cycles at −1200 kyr can be a shift from 1 × 41 to 2 × 41 kyr obliquity frequency locking, and/or 2 ×21 to 4 ×21 kyr precession locking. We base this claim on durations between major glacial terminations and a wavelet spectrum of the LR04 stack (see Fig. 1 b)); both quantities increase rather rapidly around −1200 kyr and show a steady but irregular increase towards present time.

Identifying the shift to longer periods
While Huybers (2007) observed that the mean period of global ice volume variations increases over time, we make the stronger claim that an abrupt shift from 40 to 80 kyr long durations occurred around −1200 kyr. We base this claim on our identification of major glacial terminations, which unlike spectral decomposition ignores glacial cycle shape and is unaffected by time-frequency resolution.
A disadvantage of using glacial termination events is that it is unclear what constitutes a major termination, and whether it is meaningful to characterise glacial cycles by termination events. Nevertheless, we believe that our identification of major terminations is sufficiently robust to support the claim that the duration shifted abruptly from ∼40 to ∼80 kyr around −1200 kyr.
6.2 Testing for trend after the MPT It appears that the durations between successive glacial terminations are increasing over time starting at the onset of the MPT around −1200 kyr.
We evaluate whether this trend is statistically significant, using a variation on the Mann-Kendall test (Kendall, 1955;Mann, 1945). Our null hypothesis H0 is that the sequence of thirteen durations from −1126 kyr until present is generated by a process with stationary mean, and that any observed monotonicity is by chance. SinceD i = D i − D mean , successive deviations from the mean duration D mean = 91 kyr are correlated, we immediately reject a white noise process as assumed in the standard Mann-Kendall test. Instead, we model them as an AR(1) process, such thatD i+1 = αD i + σ d ξ i , where ξ i are independent Gaussian zero mean and unit variance elements. The parameters α = 0.6 and σ d = 14.5kyr are the standard estimates of lag 1 and 0 autocorrelation coefficients respectively.
We test the hypothesis using the Kendall τ test statistic for monotonicity τ K , based on the number of ordered and disordered pairs in a sequence. τ K = 1 for a perfectly ordered sequence and τ K = 0 for sequence with equally many ordered and disordered pairs (see Appendix A for a definition of τ K ). We evaluate τ K for 2 · 10 4 samples of the AR(1) process. As indicated in Fig. 7, it is unlikely (p < 0.05) to observe the test statistic in durations from data, assuming that the durations follow an AR(1) process. Therefore, we reject the null hypothesis of no trend.
Adding age model uncertainty to the Monte Carlo sequences of durations only makes it more difficult to reject H0. Furthermore, slightly different choices of major glacial terminations, or the use of an untuned record, does not influence the conclusion of the test.

Consequences for modelling the MPT
Some explanations for the MPT do not reproduce the sequence of successively longer durations between glacial terminations in data as naturally as RFL. Instead, they produce long period cycles at the onset of the MPT which shorten towards the present as a parameter is ramped.
The Maasch and Salzman 1990 model in Fig. 8 is one such model (Maasch and Salzman, 1990). The inconsistency with data is evident when comparing the model durations with those in the LR04 stack (Fig. 1). Another such model Fig. 7 Histogram shows a Monte Carlo distribution of the Kendall tau (τ K ) test statistic, under the null hypothesis H0 that the sequence of durations between glacial terminations from −1126 kyr follow an AR(1) process. Larger τ K indicates a more monotonic sequence. Black line shows the test statistic τ K for durations in an ice volume proxy (Fig. 1). τ K under H0 exceeds the observed τ K only in 5% of the cases. For details, see Section 6.2 is the Tziperman and Gildor 2003 model (Tziperman and Gildor, 2003).
Although different dynamical mechanisms are at play in these models, they have in common that a long period limit unforced cycle emerges near a region of slow motion in phase space. As a parameter is varied, the limit cycle moves farther from this region, shortening the internal period.
RFL on the other hand naturally explains both a sudden shift from 40 kyr to 80 kyr cycles and a gradual increase towards longer cycles, since the system can respond both smoothly and abruptly to an increasing internal period, due to the Devil's staircase structure (e.g. Fig. 4). We interpret the progression of durations as evidence against models like Maasch andSalzman 1990 andGildor 2003, and for mechanisms that naturally produce increasing glacial cycle length, such as RFL.

Non-harmonic forcing
RFL is not restricted to harmonic forcing, but occurs also for astronomical, non-harmonic forcing. This is for instance the case for the Paillard and Parennin 2004 model (Fig. 10), forced by summer solstice insolation at 65 degrees North (65Nss, Fig. 9 a)). As a parameter is increased linearly, durations first cluster around 41 kyr, then shift abruptly to cluster around 82 kyr at −1000 kyr, after which they increase gradually until present. The shift to 80 kyr durations is later than in proxy data (Fig. 1) and there are some short and long durations not clear in the proxy record, but overall the glacial terminations coincide well.
Multi-frequency forcing generally produces Devil's staircases with shorter steps of constant duration, making them look "smooth" (e.g. Fig. 11). This is apparently a problem for RFL since it relies on rapid jumps in durations. However, RFL can still be relevant as demonstrated by H07 forced by an equal amount of obliquity and precession (Fig. 9 b), Fig. 12 and Fig. 11). Such forcing corresponds well to e.g. (Maasch and Salzman, 1990) forced by Summer solstice insolation at 65 degrees North. a) Global ice volume over time (black), with glacial terminations (red dots) at peaks chosen for simplicity to be above 1.2 normalised ice volume units and spaced at least 60 kyr apart. Self-sustained cycles emerge around -800 kyr and shorten towards the present. b) Durations and wavelets as in Fig. 1, except that contours start at 10% of the maximum wavelet amplitude 65Nss (au) a) Time (kyr) OblPrec (au) b) Fig. 9 Astronomical insolation curves. a) Summer solstice insolation at 65 degrees North (65Nss), normalised to zero mean and unit variance (Laskar et al., 2004). The signal is approximately a linear combination of 33% normalised obliquity and 77% normalised precession, two modulated sinusoidal signals with central frequencies 41 and 22 kyr (Crucifix, 2013). b) A normalised insolation curve consisting of 50% obliquity and 50% precession caloric summer insolation or integrated insolation above a threshold (Huybers, 2011;Tzedakis et al., 2017). The median and mode of the distribution of durations change more abruptly than the mean, which reflects that the gradual increase in average duration is caused by a gradual redistribution of durations between clusters, rather than a gradual increase of the most typical durations. In a simulation with timedependent ramping parameter, the local-in-time distribution of durations cannot be sampled well. Therefore the majority of the realised durations come from the dominant clusters of durations, which can give the impression that durations shift rapidly, in spite of the average duration changing gradually ( Fig. 12 and Fig. 11)).

Fig. 8 Simulation of the Maasch and Salzman model in
Multi-frequency forcing gives rise to many interesting phenomena regarding predictability of solutions, see for instance (Ashwin et al., 2018;Crucifix, 2013;Grebogi et al., 1984;Imbrie and Imbrie, 1980;Le Treut and Ghil, 1983;Mitsui et al., 2015;de Saedeleer et al., 2013;Tziperman et al., 2006). Importantly, however, these phenomena are not essential to RFL. Whether solutions are truly frequency locked or depend on initial conditions is irrelevant, as long as durations undergo abrupt change and tend to cluster.
We conclude that RFL, clearly understood under periodic forcing, also is relevant for astronomical forcing. Indeed, recent studies provide evidence for the long-standing hypothesis that a combination of precession and obliquity paces the glacial cycles (Feng and Bailer-Jones, 2015;Huybers, 2011;Tzedakis et al., 2017). Differences between periodic and multi-frequency forcing exist, but are not crucial for modelling the MPT with RFL.

Relevance for complex models and physical mechanisms
We see two practical uses of our description of RFL: To guide modelling of the MPT in complex models, and to drive the search for slowly changing climate variables.

Relevance for complex models
While climate physics are highly simplified in conceptual models like H07, their dynamics are well understood. The opposite holds true for Earth System Models (ESMs), which resolve multiple processes of climate in detail. To learn if the dynamical mechanism of RFL applies to such a model, we could in theory produce an Arnold tongue diagram as for H07 (Fig. 4). However, since running ESMs is computationally expensive this is presently not possible. Nevertheless, it might be possible to detect signatures of RFL from only few model runs.
First, one should investigate whether the glacial cycles are self-sustained by fixing model parameters at plausible values and fixing the insolation field at its mean value. This is the case if, after a transient time, variations in ice volume on the order of 10-100 kyr persist. The next step is to sparsely sample an Arnold tongue diagram. First, a ramping parameter must be chosen. This does not have to be a scalar, but can be a function like a parametrisation, as long as its change over time is well defined. The parameter should be one that feasibly could influence the internal period of glacial cycles.
If changing the parameter changes the internal period, then one can compare the average duration of (insolation  (Paillard and Parrenin, 2004) forced by Summer solstice insolation at 65 degrees North (Fig. 9 a)). a) Model ice volume over time (black) contrasted with the LR04 stack (blue) (Fig. 1), with glacial terminations (red dots) at times when a switch in Southern ocean circulation occurs. b) Durations and wavelets as in Fig. 1   Fig. 11 Devil's staircases of H07 forced by equal amounts of obliquity and precession (Fig. 9 b)). The average duration D (black line) as function of internal period T o is gradually increasing, whereas the the median D me di a n and the mode D mo d e are more step-like. Blue dots are the population of durations for fixed internal period; darker colours indicate higher density of durations. D mo d e is defined from binning the durations; D mo d e is the mean of the edges of the 4-kyr bin with the highest frequency. All quantities are evaluated from −2000 kyr to the present. Model parameters as in Fig. 12 variation) forced and unforced solutions. If the average duration of the forced solutions is close to either 40 or 80 kyr and remains close even under parameter perturbations that change the internal period, then this is an indication that the system is frequency locked to insolation in a way relevant for the glacial cycles. In that case, there is good reason to research RFL more closely in the model. Ashkenazy (2006) suggested that synchronisation can be detected by running the system from multiple initial conditions and see if solutions converge. This procedure is not enough for us; we need to know if the internal period can be shifted appropriately with a change in parameter, and we need to know if the durations can robustly cluster on 40 and 80 kyr.

Ramped climate variables
To evaluate whether RFL caused the MPT one must identify slowly changing climate variables. Two such candidate variables are atmospheric CO 2 and atmospheric or oceanic temperatures. Since less CO 2 leads to a generally cooler atmosphere, it can be viewed as a proxy for global average atmospheric temperature. Local cooling can occur for other reasons, however.
There are currently no direct measurements of atmospheric CO 2 across the MPT, but a recent reconstruction back to −2000 kyr suggests that the mean CO 2 did not change in the mean until at least −1300 kyr (Hönisch et al., 2009). Since the reconstruction implies that CO 2 fell 31 ppm by −700, the decrease in CO 2 must either have been rapid and driving the MPT, or a consequence of it. A rapid change in CO 2 is still consistent with RFL, but in that case RFL does not explain the abrupt increase in cycle length at the MPT, and instead one must find an explanation for the rapid increase in CO 2 . However, the planned European BEOIC deep ice core drilling in Antarctica can hopefully improve estimates of CO 2 across the MPT.
There is evidence of a gradual deep ocean cooling since the onset of northern Hemisphere glaciation 2.7 Million years ago (Lisiecki and Raymo, 2005). How much of this a) b) MPT Fig. 12 Simulation of the H07 model forced by a sum of equal amounts obliquity and precession ( Fig. 9 b)). a) Model ice volume (black) shown with the LR04 stack (blue) (Fig. 1), with glacial terminations (red dots) at times when a threshold of deglaciation is reached. b) Durations and wavelets as in Fig. 1. The threshold of deglaciation increases linearly as R(t) = 40 + 0.04(t + 2000), and forcing amplitude is A = 26. All other parameters are as in Section 2 cooling occurred across the MPT is not known, however. The reconstruction of deep water temperatures by (Elderfield et al., 2012) indicates a gradual cooling in the mean from −1300 kyr until present, but also a puzzling warming from −1500 kyr to −1300 kyr. Therefore, glacial cycle length does not appear to have a direct relation with mean deep ocean temperature. However, it may be that sea surface temperatures in the vicinity of major ice sheets are more relevant for glacial dynamics. If so, detailed and reliable reconstruction of such temperatures is necessary to evaluate whether they act as ramped climate variables in RFL.
Another slowly varying parameter could be the erosion of regolith. According to this hypothesis soft material under ice sheets eroded throughout the Pleistocene, enabling them to grow larger before collapsing Clark and Pollard (1998). The hypothesis is difficult to test empirically, however.
In addition to the candidate ramping climate variables mentioned, there may be others that are relevant for the MPT. RFL motivates the search for other such climate variables. These might not only be relevant for RFL, but for any mechanism of the MPT invoking deterministic bifurcation.

Criteria for RFL
Having demonstrated RFL in H07 and Paillard and Parennin 2004, we ask in which models RFL is most likely to be relevant.
We expect RFL in all models similar to H07, that is, models with a critical threshold of deglaciation (explicit or not), two intrinsic growth and decay states, additive forcing, and a climate variable that naturally controls the internal period.
Furthermore, RFL is facilitated by dynamics focussed on a single strongly attracting limit cycle. This is because solutions to the system f with ramped parameters then track frozen solutions of f τ well, and because it is difficult for perturbations to bring solutions away from the neighbourhood of the attractor. Crucially, a model using RFL needs a parameter that can increase the internal period by 100 kyr. The models in e.g. (Le Treut and Ghil, 1983) and (Maasch and Salzman, 1990) are therefore difficult to reconcile with RFL since the internal periods are on the order of 10 and 100 kyr respectively, and do not change much within the physical range of model parameters.
Lastly, we note that e.g. excitable systems and dissipative resonant oscillators (Crucifix, 2012) also can undergo a rapid change in durations due to frequency locking related phenomena, although they are not self-sustained oscillators. Self-sustained oscillators are distinguished by having an internal period T o (R) through which we can define Arnold tongue diagrams and Devil's staircases; for non-selfsustained oscillators we have to define these through parameters R directly. Furthermore, it has been argued that the term frequency locking should be restricted to self-sustained oscillators (Marchionne et al., 2018;Pikovsky et al., 2001), why it makes sense to define RFL for self-sustained oscillators only.

Conclusions
The glacial cycles did not enter a stationary 100 kyr world at the MPT; instead, durations between glacial terminations shifted abruptly from approximately ∼40 to ∼80 kyr around −1200 kyr, followed by a gradual increase (Fig. 1). The dynamical mechanism Ramping with Frequency Locking (RFL) naturally explains this progression of durations. As the internal period of a model glacial cycle model increases gradually, frequency locking to insolation variations causes the durations between glacial terminations to increase sometimes abruptly and sometimes gradually.
Here we described how RFL can be understood in terms of a dynamical system f and a frozen system f τ with parameters R(t) fixed at times t = τ. The average duration D τ defined for f τ provides some information about single durations in solutions x(t) to f around t = τ, but since D τ is defined asymptotically, one must interpret solutions to f in terms of f τ with care.
Model behaviour can be understood from considering parameter paths through Arnold tongue diagrams and corresponding Devil's staircases (Fig. 3, 4 and 5). These diagrams as functions of frozen time τ depend on the change in average duration D τ (T o ) as function of T o , the change in internal period T o (R) as function of R, the change in parameters R(τ) as function of time τ, and the amplitude A of the forcing.
This decomposition clarifies different ways in which the average duration can change abruptly in models of the class f . For instance, the abruptness of the change in D τ can be adjusted either by changing the forcing amplitude A or the ramping of R(t). While the effects of changing A or the ramping of R(t) are typically easy to guess, we are not aware of any general rules dictating the widths of particular Arnold tongues. Such understanding may be researched further.
RFL is relevant also for multi-frequency astronomical forcing. Multi-frequency forcing tends to make Devil's staircases less abrupt, but durations can still increase rapidly when a model parameter is slowly ramped.
The RFL mechanism provides an explanation for the MPT without the climate system entering a new mode of operation. A shift from ∼40 kyr long to ∼80 kyr long cycles due to frequency locking to obliquity and precession, is consistent with data ( Fig. 1), and is used in models (Huybers, 2007;Paillard and Parrenin, 2004). This warrants further study of frequency locking characteristics of models throughout the model hierarchy, as well as a search for gradually increasing climate parameters. Some models use a rapidly ramped parameter to accelerate the increase in durations between glacial terminations at the MPT, but such a ramping begs for justification that a model relying solely on frequency locking does not require.

B Wavelets
Wavelet spectra are estimated with the MATLAB®function cwt, using Morlet basis functions with bandwidth parameter ω 0 = 6 (Torrence and Compo, 1998). Contours show wavelet amplitude (square root of variance) relative to the maximum, incremented in evenly spaced percentage units. The cone of influence marks the e-folding time of the amplitude of a discontinuity at the edge of the time interval. Inside the cone of influence edge effects are negligible (Torrence and Compo, 1998).