The Mid-Pleistocene Transition induced by delayed feedback and bistability

The Mid-Pleistocene Transition, the shift from 41 kyr to 100 kyr glacial-interglacial cycles that occurred roughly 1 Myr ago, is often considered as a change in internal climate dynamics. Here we revisit the model of Quaternary climate dynamics that was proposed by Saltzman and Maasch (1988). We show that it is quantitatively similar to a scalar equation for the ice dynamics only when combining the remaining components into a single delayed feedback term. The delay is the sum of the internal times scales of ocean transport and ice sheet dynamics, which is on the order of 10 kyr. We find that, in the absence of astronomical forcing, the delayed feedback leads to bistable behaviour, where stable large-amplitude oscillations of ice volume and an equilibrium coexist over a large range of values for the delay. We then apply astronomical forcing. We perform a systematic study to show how the system response depends on the forcing amplitude. We find that over a wide range of forcing amplitudes the forcing leads to a switch from small-scale oscillations of 41 kyr to large-amplitude oscillations of roughly 100 kyr without any change of other parameters. The transition in the forced model consistently occurs near the time of the Mid-Pleistocene Transition as observed in data records. This provides evidence that the MPT could have been primarily a forcing-induced switch between attractors of the internal dynamics. Small additional random disturbances make the forcing-induced transition near 800 kyr BP even more robust. We also find that the forced system forgets its initial history during the small-scale oscillations, in particular, nearby initial conditions converge prior to transitioning. In contrast to this, in the regime of large-amplitude oscillations, the oscillation phase is very sensitive to random perturbations, which has a strong effect on the timing of the deglaciation events.

The Mid-Pleistocene Transition, the shift from 41 kyr to 100 kyr glacialinterglacial cycles that occurred roughly 1 Myr ago, is often considered as a change in internal climate dynamics. Here we revisit the model of Quaternary climate dynamics that was proposed by Saltzman and Maasch ('Carbon cycle instability as a cause of the late Pleistocene ice age oscillations: modelling the asymmetric response'. Glob Biogeochem Cycle 1988; 2: 177-185-from this point, referred to as SM88). We show that it is quantitatively similar to a scalar equation for the ice dynamics only when combining the remaining components into a single delayed feedback term. The delay is the sum of the internal time scales of ocean transport and ice sheet dynamics, which is on the order of 10 kyr.
We find that, in the absence of astronomical forcing, the delayed feedback leads to bistable behaviour, where stable large-amplitude oscillations of ice volume and an equilibrium coexist over a large range of values for the delay. We then apply astronomical forcing using the forcing data for integrated summer insolation at 65 degrees North from Eisenman (Integrated Summer Insolation Calculations. NOAA/NCDC Paleoclimatology Program Data Contribution # 2006-079, 2006). Since the precise scaling of the forcing amplitude is not known, we perform a systematic study to show how the system response depends on the forcing amplitude. We find that over a wide range of forcing amplitudes the forcing leads to a switch from smallscale oscillations of 41 kyr to large-amplitude oscillations of roughly 100 kyr without any change of other parameters. The transition in the forced model consistently occurs at about the same time as the Mid-Pleistocene Transition (between 1200 and 800 kyr BP) as observed in the data records from Lisiecki and Raymo ('A Pliocene-Pleistocene stack of 57 globally distributed benthic δ 18 O records '. Paleoceanography 2005;20:1-17).. This provides evidence that the MPT could have been primarily forcing-invoked. Small additional random disturbances make the forcing-induced transition near 800 kyr BP even more robust.
We also find that the forced system forgets its initial history during the small-scale oscillations, in particular, nearby initial conditions converge prior to transitioning. In contrast to this, in the regime of large-amplitude oscillations, the oscillation phase is very sensitive to random perturbations, which has a strong effect on the timing of the deglaciation events.  (Lisiecki and Raymo, 2005).

I. INTRODUCTION
The Pleistocene is characterised by its successive glacial-interglacial cycles whose documented periodicities has been a subject of interest for decades. Hays, Imbrie, and Shackleton (1976); Paillard (1998), and Ganopolski and Calov (2011) have all presented evidence that orbital forcing plays a large role in the onset and termination of glacial periods, but the full extent of the variations cannot be explained solely by forcing (Ridgwell, Watson, and Raymo, 1999;Shackleton, 2000;Paillard, 2001). One actively researched aspect of the Pleistocene is the abrupt change in frequency of major glaciations known as the Mid-Pleistocene Transition (MPT) (Pisias and Moore, 1981;Maasch, 1988;Maasch and Saltzman, 1990;Paillard, 1998;Clark et al., 2006;Crucifix, 2012;Ashwin and Ditlevsen, 2015). Fig. 1 shows a proxy record compiled by Lisiecki and Raymo (2005) for global temperature taken from ocean core sediments which can be directly related to global ice cover. Spectral analysis has been performed on this time series, and it has been observed that the signal of the 100-kyr cycle began to rise 1250 kyr BP and was established as the dominant cycle by 700 kyr BP (Clark et al., 2006;Lisiecki and Raymo, 2007;Dijkstra, 2013). Many researchers have tried to determine not only the cause of this switch, but also the driving force behind the 100 kyr cycles themselves (Maasch and Saltzman, 1990;Saltzman and Maasch, 1991;Paillard, 1998;Gildor and Tziperman, 2001;Paillard and Parrenin, 2004;Ganopolski and Calov, 2011;Ashwin and Ditlevsen, 2015). While the 41-kyr oscillations have been attributed to external forcing (Paillard, 2001), the 100-kyr cycles have been proposed to be a result of nonlinear responses of the climate system itself (Imbrie et al., 1993;Yiou et al., 1994;Gildor and Tziperman, 2001).
The major glacial-interglacial cycles were originally attributted to changes in insolation due to orbital variations (i.e. the Milankovitch cycles). Milankovitch argued that summer insolation at high northern latitudes determined the main pacing of glaciations, as these changes control the length of seasons and amount of solar energy being recieved in high latitudes where the ice resides. The most prominent frequencies of these changes are 19 and 23 kyr due to precession and 41 kyr due to obliquity (Paillard, 2001). Precession is the orientation of the rotational axis of the earth (axial) and the rotation of the orbital axis (apsidal), and obliquity is the angle between the rotational axis and the orbital axis. Milankovitch (1941) argued that the main driving forces of climate fluctuations are obliquity and precession, which agrees with the records prior to the mid-Pleistocene transition. There is a smaller signal of orbital forcing that correlates with the 100-kyr oscillations (Imbrie et al., 1993). Eccentricity, the amount the earth's orbital ellipse deviates from a circle, shifts with a main period of 413 kyr, but there are components that vary with periods between 95 kyr and 125 kyr. The amplitude of these forcing signal bands, however, are an order of magnitude smaller than the signals of the 23-and 41-kyr bands (Hays, Imbrie, and Shackleton, 1976).
Since Milankovitch's theory, there have been many other attempts to identify the cause of the MPT (see review paper by Crucifix (2012)). The starting point of our analysis is the collection of Saltzman models from the late twentieth century. They were a collection of attempts to model the transition as a bifurcation in low-order dynamical systems due to climate feedbacks. The models typically only had three dynamic variables usually representing some measure of global ice volume, global mean temperature, and ocean circulation strength. Other modelling attempts, after the Saltzman models, include (in order of in-creasing complexity) conceptual models (Paillard, 2017), threshold models (Tzedakis et al., 2017), nonsmooth models (Paillard, 1998), relaxation oscillators (Ashwin and Ditlevsen, 2015), box models of intermediate complexity (Gildor and Tziperman, 2001;Tziperman and Gildor, 2003), and Earth system models of intermediate complexity (EMICs) (Ganopolski and Calov, 2011). While all of these studies were able to reproduce a transition similar to the MPT, they are all based on the hypothesis that the MPT occured because of some significant change in the dynamics of the system. Huybers (2009) proposed an alternative hypothesis that the MPT was a spontaneous response to the astronomical forcing. In this study we focus on this alternative hypothesis, and provide more evidence using a previously established model by Saltzman and Maasch (1988).
The paper is arranged as follows. Section 2 revisits and analyses the three-variable Saltzman and Maasch 1988 model. We observe that two of the variables mostly act as delays in feedback loops such that we may reduce the model to a scalar equation with delayed feedback. We also identify a region of bistability at the parameter values of interest which was not explored in the original model analysis. Section 3 shows the model response in the bistable region when astronomical forcing is introduced, which includes MPT-like realisations with no change in parameters. In section 4 we investigate the sensitivity of the model. We observe that the system forgets its initial condition and has low sensitivity whenever it exhibits a small-amplitude oscillatory response. It is much more sensitive to disturbances when exhibiting a large-amplitude oscillatory response. This sensitivity affects mostly the timing of major deglaciation events. Section 5 discusses our results.

II. ADAPTATION OF THE SALTZMAN AND MAASCH MODEL
We begin by considering the original model of Saltzman and Maasch (1988), which we will refer to as SM88,Ẋ Here X, Y, Z are scaled versions of global ice mass, atmospheric CO 2 , and North Atlantic Deep Water (NADW) respectively. More precisely, the variables represent anomalies from a background state. The first term in (1a) represents the feedbacks of global ice mass. The combined effect of damping and negative feedback from NADW production is taken to outweigh the positive feedback from ablation (melting and loss of ice through icebergs). The second term in (1a) is the direct effect of higher temperatures (caused by higher CO 2 levels) leading to loss of ice mass. In equation (1c) the negative sign of ice anomaly X is motivated by the reduction of NADW production with loss of ice mass. This equation has a time scaling q which is the ratio of ice sheet time constant (10 kyr, which equals t = 1 in the scaling of (1a)-(1c)) to the response time of the deep ocean. The equation for atmospheric CO 2 (1b) has a negative term −pZ for the negative effect due to NADW-more NADW leads to more ventilation of the deep ocean (stronger mixing between surface and deep waters) and an increased downdraw of atmospheric CO 2 into the deep ocean. The term rY accounts for the positive feedbacks of atmospheric CO 2 , related to changes in sea surface temperatures, sea ice extent, and sea level, which outweigh negative feedbacks. The nonlinearity (s−Y )Z 2 is motivated by the appearance of locally enhanced instabilities in the Southern Ocean due to increased NADW production. The NADW meets colder, denser Antarctic Bottom Water and induces vertical mixing. This brings CO 2 -rich water to the surface which releases CO 2 into the atmosphere depending on the current level of atmospheric CO 2 (Y ). Each of these effects (except damping of the nonlinear terms) have an associated parameter which controls their relative strengths. These were the feedbacks between global ice mass, atmospheric CO 2 , and NADW deemed most important by Saltzman and Maasch (1988). Present day studies continue to stress the interaction between these three variables as being essential for the MPT, including a study by Chalk et al. (2017) which uses a geochemical box model to confirm that carbon cycle feedbacks related to ocean circulation are necessary to sustain the long glacial cycles in the late Pleistoncene.
where V is the negative of the global ice mass perturbations, replaces the two linear equations (1a) for X and (1c) for Z by a chain of linear first-order filterṡ (2) A linear chain of filters is well known to approximate a delay (Smith, 2010): Appendix A gives further comments on the linear chain approximation for delays. Using this connection between linear chains and delays, we can rewrite the model as a scalar delay-differential equation (DDE) for Y . To facilitate the comparison to data at a later stage, we shift time and exploit that in the DDE, the ice mass anomaly V = −X is just the delayed (by time 1) CO 2 anomaly: . Thus, we may express the DDE model in terms of ice mass anomaly X: In DDE (4) one of the parameters r, p, s, τ is redundant and could be removed by a rescaling of X and time. Thus, (4) is objectively simpler than the original SM88 model (1) as it has fewer free parameters. For the purpose of model comparison, the DDE system will be left as above.

A. Dynamics of DDE, compared to SM88
We initially perform a systematic analysis of DDE (4) and compare it to the original SM88. We find that models (1) and (4) have the same long-time behaviour (using the relation q = 1 τ −1 implied by (3)). Since delays, whether in the form of chains as in SM88 or explicit, don't affect location of equilibria, i.e. X(t) = X(t − τ ) (or X(t) = −Y (t) = −Z(t)), equilibria for both models, (1) and (4), (called X eq,j ) satisfy 0 = −pX eq + rX eq − sX 2 eq − X 3 eq .
Thus, we have a trivial equilibrium X eq,1 = 0 for all parameter values and non-zero equilibria X eq,2,3 for s 2 > 4(p − r): The stability of the equilibria may differ, however, since eigenvalues of the linearisation may cross the imaginary axis (leading to oscillations in a Hopf bifurcation (Kuznetsov, 2013)) at different parameter values for (1) and (4). All other trajectories have to be studied numerically (for both, SM88 (1) and DDE(4)). Since the scaling between the two models are the exactly same, we will use values for p, r, and s from Saltzman and Maasch (1988) for all numerical studies: while varying the timescale (or delay) of the feedback processes (τ for DDE (4), q for SM88 model (1)). A value of q = 1.2 was used in SM88, which would correspond to τ ≈ 1.83. The only constraint on q was q ≥ 1 (q was defined as the ratio of ice sheet time constant, 10 kyr, to the slow response time of deep ocean FIG. 2: Bifurcation diagrams of both models for delay parameter τ . Other parameters: ocean response time of around 8.3 kyr. This is long for ocean timescales as is shown through simulations of general circulation models. For example, Yang and Zhu (2011) studied the response of the deep ocean under varying climate scenarios and found it to be 1.5 to 2 kyr. While arguments for a slightly longer ice sheet response or ocean response could be made, τ < 1.6 is more realistic that the suggested τ = 1.83 in SM88. The comments below (4) also show that the change τ → τ /α, p → αp, r → αr, s → √ αs leaves the system unchanged such that all phenomena reported in this paper are observable for any delay and appropriately rescaled parameters. Figure 2 shows the bifurcation diagrams for varying τ in DDE (4) and SM88 model (1), respectively. Although the locations of bifurcations are slightly shifted, qualitatively the two figures agree nicely. In particular, even though the space of possible initial conditions for DDE (4) is infinite-dimensional (every possible history on [−τ, 0] gives a different trajectory), the long-time behaviour of trajectories in the (X(t), X(t − τ ))-plane follows a two-dimensional ordinary differential equation (ODE) in the range we explored. The diagrams in Fig. 2  ] can be split further into two sections: one containing unstable periodic orbits and the other not. The effect of the unstable periodic orbits is a drastically reduced basin of attraction for the stable equilibrium. However, since there is no change in the attractors, we will consider these as one region. Region [r el ] was not discussed in a later parameter study of the original model (Maasch and Saltzman, 1990). For the remainder of this paper we will discuss only the DDE model, but similar results can be seen for the ODE model.

B. The bistable region [r el ]
The focus of this study is the bistability seen for τ ∈ [1.295, 1.625] in the DDE system. Not only is this a novel dynamical region, but it also fits within more realistic timescales of the ice sheet and deep ocean response. In this region there are two possible stable solutions (shown in Fig. 3): a stable equilibrium and stable large amplitude periodic orbits. The time profile of the periodic orbits has the assymetrical shape observed in the ice age cycles of the late Pleistocene. The assymetry is attributed to a slow accumulation of ice mass followed by rapid melting. In addition, the period remains between 109 and 120 kyr throughout the bistable region (with an exception of τ very close to the [r e -r el ] boundary where the period approaches infinity). This cycle length agrees with what is seen in the data (even more so when one adds the external forcing; see Sec. III). In previous studies of SM88 model (1), a transition between two stable states was enforced by a parameter shift through a Hopf bifurcation (Saltzman and Maasch, 1988;Maasch and Saltzman, 1990). We will demonstrate in the next section that the bistability in region [r el ] makes transitions between the two states possible without any change of parameter when the model is subjected to external forcing.

III. THE MID-PLEISTOCENE TRANSITION UNDER THE FORCED MODEL
In this section we show that when adding astronomical forcing to model (4), a transition typically occurs at the same time as the MPT is seen in recorded data without further parameter tuning. In section III A we describe the astronomical forcing considered and how we include it in the model. We then examine the responses for different delays and different forcing strengths in sections III B and III C.
A. Astronomical forcing Hays, Imbrie, and Shackleton (1976) have provided evidence that the glacial cycles during the Pleistocene are driven primarily by variations in the earth's orbital cycle. This includes FIG. 4: Normalised integrated July insolation at 65 • N, adapted from Huybers and Eisenman (2006). changes in precession (orientation of the rotational axis), obliquity (angle between the rotational axis and orbital axis), and eccentricity (orbital ellipse's deviation from a circle). These modes vary approximately periodically, with cycle lengths of 19/23 kyr, 41 kyr, and 100 kyr respectively (Milankovitch, 1941;Berger, 1978;Huybers and Eisenman, 2006). Fig. 4 shows a time series of average daily summer insolation at 65 • N computed by Huybers and Eisenman (2006) based on the model introduced by Huybers (2006). To investigate the effect of this forcing on DDE (4), we include the astronomical forcing in the same way as in Maasch and Saltzman (1990). We add forcing signal M (t) shown in Fig. 4 with negative amplitude u,Ẋ The precise procedure for extracting M (t) from the publicly available data source can be found in Appendix B. We are interested in how the bistable region responds to this external forcing in dependence of two parameters: the delay τ and the forcing amplitude u.
For small u the system is expected to exhibit two types of responses to forcing in the bistable region. Each of them is a perturbation of an attractor of the unforced system, namely the equilibrium and the large-amplitude periodic orbit, which persist for small u.
We will refer to these responses as the small-amplitude and the large-amplitude response (compare red and blue time profiles in Fig. 5c). For increasing u we expect to observe increasingly frequent transitions between these responses. For large forcing amplitudes u the internal dynamics of the model will be dominated by the forcing. The first comparison is made close to the boundary [r e ]-[r el ] at τ 1 = 1.3. We see in figure 5a that both trajectories change in synchrony, exhibiting only the small-amplitude response. In the middle of the bistable region (τ 2 = 1.45, shown in Fig. 5b), the solution shows a small-amplitude response in most of the first half of the time window and a largeamplitude response in the second half. Close to the Hopf bifurcation of the autonomous stable equilibrium (τ 3 = 1.6, shown in Fig. 5c) the solution exhibits primarily a largeamplitude response. This is expected due to the weakening attraction of the autonomous stable equilibrium and its shrinking basin of attraction. We will demonstrate in the following section that these transitions generate dynamics with time profiles that are qualitatively similar to the records of the MPT.

C. Variable forcing strength
For the exploration of the effect of different forcing strengths, we take the same values for τ as in Fig. 5, covering the range of the bistable region. For each τ we compute trajectories for different u, ranging from u = 0 to u = 0.6, with the initial history X(t) = −0.5 for t ∈ [2 + τ, 2]Myr BP. Figure 6 shows the difference of these trajectories to a respective reference trajectory for τ ref = 1.25 and the same initial condition and forcing strength u, averaged over a window of length τ . We note that differences to a reference trajectory at u = 0 and identical delay τ (using same initial condition) give qualitatively similar results.
Remarkably, figures 6a and 6b show that there is a distinct period around 700-800 kyr BP where the solutions diverge from the reference trajectory in a large range of forcing strengths u. This suggests that some aspect of the forcing around this time kicks the trajectories into the basin of attraction of the large-amplitude response. An additional area like this is seen in figure 6b near 1600 kyr BP.
The second interesting feature in all three figures is a threshold behaviour in the forcing strength u: below a certain value of u specific to each τ , transitions do not occur such that This forcing has a period of 41 kyr. If u > 0, equation (8) has two stable solutions: a sinusoidal signal with a 41 kyr period, or a large-amplitude asymmetrical quasiperiodic response. When recreating figure 6b with this periodic forcing, we see the same threshold transition at u = 0.08. We can show that this transition is due to moving basins of attraction for the two stable solutions. The other sharp transition we saw in Fig. 6 would not be possible with this forcing. Either the solution would start to transition immediately, or it would not transition at all. This is also the case if an envelope of modulated amplitude describing changes of eccentricity at 1.1Myr periods is added for the periodic forcing. This leads us to believe that the quasiperiodicity is necessary for a transition of this type.

A. Finite-time Lyapunov Exponents
To analyse how trajectories depend on their history at specific instances of the forcing, we compute finite-time Lyapunov exponents (FTLEs) using a QR decomposition method (details in Appendix C). This method was previous used in De Saedeleer, Crucifix, and Wieczorek (2013) to illustrate the desynchronisation of nearby trajectories in a van der Pol-type oscillator model. We will apply the ideas presentated in that study to our forced model.
The difference between FTLEs and classical Lyapunov exponents is that FTLEs are recorded for a family of time windows of finite length w rather than over the entire long time run. Thus, FTLEs are time-dependent functions instead of real numbers. A positive FTLE along a given trajectory X(t) at time t 0 indicates that some nearby trajectories diverge exponentially from X(t) in the time window [t 0 −w, t 0 ]. Therefore, X(t) is sensitive to small perturbations in the time window [t 0 − w, t 0 ]. While a trajectory could be asymptotically stable, a positive FTLE at a time t 0 indicates temporary amplification of perturbations from the attracting trajectory (observed as temporary desynchronisation, see (De Saedeleer, Crucifix, and Wieczorek, 2013)).

B. FTLE implications on MPT and timing of major deglaciations
We analyse one trajectory, showing a forcing-induced MPT (τ = 1.45, u = 0.15,shown in Fig. 7) . We compute the FTLE over a sliding window of our example MPT tragectory with a window length of w = 250 kyr. This window length is chosen in order to filter out the dominant frequencies of the forcing. Similar results are seen with any window lengths w from 150-500 kyr. Fig. 8 shows the time profile of the largest FTLE along the trajectory. Before the transition around 800 kyr BP, the FTLE generally remains negative apart from a few short excusions above zero. At 1000 kyr BP the FTLE approaches zero and remains there for some time. Just before 800 kyr BP it goes positive and then on average stays positive for the remainder of the trajectory.
The negative FTLEs leading up to the transition confirm that the trajectory forgets its initial history and the effect of past disturbances. In particular, this implies that the infinite-dimensional nature of the DDE's possible initial conditions does not play a role for the MPT. Whenever the system is exhibiting the small-amplitude response we observed negative FTLEs.
The positive FTLEs indicate a sensitivity of the trajectory during the large-amplitude response. This sensitivity affects the precise timing of the deglaciations, as we now demonstrate by noise-induced desynchronisation of nearby trajectories.
To explore further the desynchronisation phenomenon outlined in De Saedeleer, Crucifix, and Wieczorek (2013), we add noise to our system. The stochastic delay differential equation (SDDE) is then given as Here, W (t) is standard white noise and σ is the noise amplitude. We will always consider σ as a fraction of the forcing amplitude u, i.e. σ = u 30 . We set the deterministic forcing strength u = 0.15 and the noise amplitude σ = 0.005. We ran 500 realisations of the model, all with the same initial history. The results of 10 randomly selected realisations can be seen in Fig. 9. Up until 1000 kyr BP all trajectories generally track the same solution with only minor short desynchronisations. Shortly after 1000 kyr BP the trajectories begin to diverge, making the transition to the large amplitude oscillation state at different times. The trajectories then stay desynchronised, which corresponds to being at a different phase along the large amplitude oscillation. We illustrate this phenomenon in Fig. 10, where we show the distribution of trajectories in the (X(t), X(t − τ ))-plane along the unforced periodic orbit, and in Fig. 11. Figure 11 shows the transition from synchronization to desynchronization in a density plot. At approximately 800 kyr BP (indicated by the gray line) the probability density makes a sharp transition from a small-variance to a large-variance-low-maximum density. The timing of  Figures 10 and 11 also show that the distribution by desynchronisation is far from uniform. There are still distinct deglaciation times favoured for many realisations. For this reason the desynchronisation observed in figures 10 and 11 does not contradict studies that have argued for phase-locking to different components of the orbital forcing. One of the original hypotheses of phase-locking in the late Pleistocene was suggested to be related to eccentricity, specifically associated to events of low eccentricity (Hays, Imbrie, and Shackleton, 1976;Paillard, 2015). More recent studies have looked at the possibility of locking to precession or obliquity. In Ridgwell, Watson, and Raymo (1999) the authors argue locking to every 4th or 5th precession cycle, while Huybers and Wunsch (2005) argue for locking to every 2nd or 3rd obliquity cycle. Later, Huybers (2011) attributes locking to a combination of precession and obliquity, and states strongly that the pacing cannot be attributed solely to one of these components. Robust evidence for locking requires testing on a range of initial conditions or noise realizations (or long time series), but the short data record corresponds to only one realization of noise disturbance and one initial condition. Thus, data available may be insufficient to distinguish the higher-order locking proposed in the literature from the level of desynchronization we report in figures 10 and 11. C. Effect of noise with stronger forcing Figure 12 shows that the addition of noise enhances the MPT for strong forcing. When adding noise, we observe transitions for forcing strengths u for which there was no transition in the deterministic case (compare Figures 6b and 13). In Fig. 6b the last persistent transition occurs at u = 0.33. Fixing the forcing strength at u = 0.45, we compute ten realisations. Here, two realisations exhibit transitions that persist until the end of the run. Although the transitions don't appear as commonly as in the weaker forcing case, with noise it is possible to observe an MPT-like transition across a wider range of forcing amplitudes. Figure 13 shows a systematic overview of the transition enhancing effect of noise. We added noise of amplitude σ = u 30 and compute the distance diagram as in Fig 6b. We observe transitions occuring above the maximal value of u for which transitions occured in the deterministic case.

V. DISCUSSION
We have presented an analysis of the well known Saltzman and Maasch (1988) model of Quaternary climate dynamics. We have shown that it can be reduced to a scalar equation for the ice mass anomaly X(t), by collecting all other variables into a delayed feedback term. Through this formulation we are able to explore the dependence on the time scale τ of the delay in this feedback. We discovered a region in which a stable equilibrium and stable large-amplitude periodic orbit coexist for the unforced system. This parameter region is more consistent with physically relevant timescales for the ice sheet and deep ocean response than those explored in SM88.
We observe two different responses to astronomical forcing of the type propsed by Huybers and Eisenman (2006): small-amplitude response (tracking the equilibrium) and largeamplitude response (tracking the periodic orbit) In the deterministic bistable region the model exhibits a switching behaviour between these two responses. The timings of the switches are robust for different values of the delay and focing strength. The most promi-nent switch occurs around 800 kyr BP which is within the range of when the Mid-Pleistocene Transition is believed to have occurred.
One trajectory showing the MPT is analysed using finite-time Lyapunov exponents. Using a window of 250 kyr we observe that the largest FTLE is negative up to 1000 kyr BP, remains around zero during the transition, and is mostly positive after the transition. We demonstrate that this positivity induces phase desynchronisation which is then confirmed by adding noise to the model. The added noise also allows for the transition to be seen with stronger forcing.
The novelty of this study is in the observation that transitions similar to the MPT occur robustly without any additional climate forcing, such as a nonstationary background state. The transition arises from the bistability rather than a shift of a system parameter through a bifurcation. We also highlight that with noise the transition is possible but doesn't always occur, and that the phase (or timings of the deglaciations) are variable (however still influenced by the forcing frequencies). While our demonstration was done for the specific low-order model SM88, we conjecture many of the effects (instability caused by delayed feedback, transition induced but not caused by the forcing) to be present in other models with similar bistability.
Many models assume a slowly varying background state, typically related to slowly declining CO 2 levels, and attribute this as the primary mechanism for the MPT (Saltzman and Maasch, 1991;Tziperman and Gildor, 2003;Ashwin and Ditlevsen, 2015;Tzedakis et al., 2017;Paillard, 2017). Our proposed mechanism does not exclude the possibility of this changing background state. However, it does not require the slow change of background state to cause a bifurcation. Instead, we consider a new null hypothesis for the MPT which involves only bistability and external (orbital and/or stochastic) forcing. Future work could include an additional varying background state within the parameter regions mapped in Fig. 2a and analyse the effects on the transition.
We reported results for one specific astronomical forcing time series (65 • N summer integrated insolation), which is a weighted average version of the typical insolation threshold forcing (see Appendix B). Qualitatively identical results can also be found when choosing a particular insolation threshold. However the particular threshold at which the behaviour is observed changes for different delays. The weighted average over thresholds from Appendix B does not depend on an additional parameter (forcing threshold) and the dynamical behaviour is preserved.
Some questions that arise from this study will need to be addessed by deeper mathematical analysis. Why does the model prefer certain times to switch, and what causes the sharp threshold in forcing strength where switching becomes possible, visible in Fig. 6? The latter can be answered by considering simple harmonic forcing with the dominant period of 41 kyr seen in the astronomical forcing. This threshold should be analysed further in a follow-up study. Another question to consider is the effect of human activities on the system. Is it possible anthropogenic contributions could knock the system out of the oscillations again, and perhaps to a completely different dynamical region? Addressing this question would be useful for understanding current and future climate scenarios.
Appendix A: The linear chain approximation A linear chain of first-order ODEs approximates a delay in the following sense (see Smith (2010) for a precise derivation). Here we adapt the derivation to our particular system. The 3-dimensional SM88 model (1), using (2), has a 2-dimensional linear chain and is of the following form: Extracting the linear chain system we have where y = [V, Z] T . We have A and b(t) as follows: Since A is invertible, the fundamental solution matrix of (A2) is then given by where P is the matrix of eigenvectors associated with eigenmatrix Λ, We have the following basis of A and it's inverse, From Perko (2013) we have the fundamental solution Assuming the solution has existed arbitrarily far in the past and using a change of variables we can extract an equation for Z(t), where the kernel in (A8) is Approximating the true distributed delay kernel K in (A8) with a delta distribution at its expected value we approximate (A1) by a DDE with a discrete delay Appendix B: Extraction of integrated summer insolation forcing M (t) from data We use the Integrated Summer Insolation Calculations data set provided by Huybers and Eisenman (2006) (found at https://www1.ncdc.noaa.gov/pub/data/paleo/climate_ forcing/orbital_variations/huybers2006insolation/j_65north.txt), last accessed: 17 October 2018. The data is provided as daily average summer energy in GJ/m 2 calculated from number of days the insolation was above a given threshold (W/m 2 ). We used the data provided for the latitude of 65 • North. Rather than selecting a particular threshold, c. QR method We use a continuous QR algorithm for computing the Lyapunov exponents (see Dieci, Russell, and Van Vleck (1997)). The QR algorithm is based on the numerical linear algebra factorization of a matrix M into an orthogonal matrix Q and an upper triangular matrix R. Although DDEs have an infinite number of Lyapunov exponents, this method allows us to compute up to m + 1 through the discretisation. We can choose to compute only the largest l Lyapunov exponents by using a non-square Q. We initialise an arbitrary othogonal Q 0 as an m × l matrix of form Q 0 = I l 0 (m−l)×l .
Note that I l is the l×l identity matrix and 0 (m−l)×l is an (m−l)×l zero matrix. We define interatively Q i and R i by the QR decomposition of M i Q i−1 (matlab command qr(M i Q i−1 , 0)), which produces a square l × l upper triangular matrix R i with eigenvalues R i,jj > 0 (j = 1, ..., l). We store R i for each timestep. After N timesteps we have the relation The infinite Lyapunov exponents can then be approximated by For FTLEs we truncate the above sum after W timesteps, and view how this truncation changes in time, i.e. λ j,n = 1 hW n i=n−W ln R i,jj .
Here W = w/h, where w is the desired time window length. Note that n initialises at time W .