Why could ice ages be unpredictable?

It is commonly accepted that the variations of Earth's orbit and obliquity control the timing of Pleistocene glacial-interglacial cycles. Evidence comes from power spectrum analysis of palaeoclimate records and from inspection of the timing of glacial and deglacial transitions. However, we do not know how tight this control is. Is it, for example, conceivable that random climatic fluctuations could cause a delay in deglaciation, bad enough to skip a full precession or obliquity cycle and subsequently modify the sequence of ice ages? To address this question, seven previously published conceptual models of ice ages are analysed by reference to the notion of generalised synchronisation. Insight is being gained by comparing the effects of the astronomical forcing with idealised forcings composed of only one or two periodic components. In general, the richness of the astronomical forcing allows for synchronisation over a wider range of parameters, compared to periodic forcing. Hence, glacial cycles may conceivably have remained paced by the astronomical forcing throughout the Pleistocene. However, all the models examined here also show a range of parameters for which the structural stability of the ice age dynamics is weak. This means that small variations in parameters or random fluctuations may cause significant shifts in the succession of ice ages if the system were effectively in that parameter range. Whether or not the system has strong structural stability depends on the amplitude of the effects associated with the astronomical forcing, which significantly differ across the different models studied here. The possibility of synchronisation on eccentricity is also discussed and it is shown that a high Rayleigh number on eccentricity, as recently found in observations, is no guarantee of reliable synchronisation.


Introduction
showed that southern ocean climate benthic records exhibit spectral peaks around 19, 23-24, 42 and 100 thousand years (thousand years are henceforth denoted 'ka'). More or less concomitantly Berger (1977) showed, based on celestial mechanics, that the power spectrum of climatic precession was dominated by periods of 19, 22 and 24 ka, and that of obliquity was dominated by a period of 41 ka. These authors concluded that the succession of ice ages is somehow controlled by the astronomical forcing.
The much less cited paper by Birchfield and Weertman (1978) is, however, at least as important. These authors considered a dynamical ice sheet model, which they forced by astronomically-induced variations in incoming solar radiation (insolation). They managed to reproduce grossly the spectral signature found by Hays et al. (1976). However, subtle changes in the model parameters, well within the range allowed by physics, disturbed significantly the precise sequence of ice ages, without altering the power spectrum of ice volume variations.
It is this author's experience that some patient tuning is generally needed to reproduce the exact sequence of glacial-interglacial cycles with a conceptual model. On Figure 1 are shown two examples of ice volume history reproduced with models previously published in the palaeoclimate modelling literature (Saltzman and Maasch, 1991;Tziperman et al., 2006). In both cases, small changes in model parameters do, at some stage in the climate history, induce a shift in the sequence of ice ages. Sometimes this sensitivity is explicitly acknowledged by the authors (Paillard, 1998;Imbrie et al., 2011), but not always, and this may have given the false impression that these models unambiguously confirm the tight control of astronomical forcing on ice ages.  Saltzman and Maasch (1990) and Tziperman et al. (2006), forced by normalised insolation at 65 • N. The blue lines are obtained with the published parameters; the latter were slightly changed to obtain the red ones: p 0 = 0.262 Sv instead of 0.260 Sv in Tziperman et al. (2006), and w = 0.6 instead of w = 0.5 in Saltzman and Maasch (1990). While the qualitative aspect of the curves are preserved, the timing of ice ages is affected by the parameter changes.
Our purpose here is to understand the dynamical factors which may induce instability in the succession of ice ages. The approach is dynamics-oriented: we use tools from mathematics, and focus more on the understanding of the dynamics, than on the identification of physical mechanisms. Though, it will not be concluded whether or not glacial-interglacial cycles are indeed predictable or not. This requires an additional step of statistical inference, which is left for another article.
Which model to use? There are many models of ice ages, spanning different orders of complexity and based on different physical interpretations. We will therefore work with different models, but only of the class of the simplest ones. This choice offers us a greater flexibility in analysing model dynamics with computing intensive techniques, and it also allows us to keep our hypotheses to a minimum.
Indeed, most of the simplest models of ice ages Maasch, 1990, 1991;Paillard and Parrenin, 2004;Tziperman et al., 2006;Imbrie et al., 2011) share a number of characteristics: 1. These are dynamical systems: climate has a memory (in contrast to Milankovitch (1998)); 2. the astronomical forcing is introduced as an additive or quasi-additive forcing term, which involves a combination of precession and obliquity; 3. there are non-linear terms involved in the internal system dynamics, which induce episodically conditions of instability. The general hypothesis is that high glaciation levels are unstable. The instability conditions may lie implicitly in the system dynamics (Saltzman and Maasch, 1990), or postulated explicitly by means of a threshold criteria (as in Paillard, 1998;Imbrie et al., 2011). The threshold may be a function of precession and obliquity (Parrenin and Paillard, 2012), and a dependency on eccentricity was also proposed (Rial, 2004).
System instability is an important aspect of Pleistocene theory. It is a convenient starting point to explain the existence of large climatic fluctuations such as deglaciations, even when the astronomical forcing is weak. Termination V, which occurred 400 ka ago, is an often-cite example (Paillard, 2001). Instability may also explain the emergence of 100-ka climatic cycles independently of the effect of eccentricity (see, e.g. Crucifix (2012) for a review). The present article is structured as follows. The discussion starts with the van der Pol oscillator forced by astronomical forcing. As the other models cited so far, this is a dynamical system that combines the accumulative action of astronomical forcing with an instability mechanism causing regime changes. The van der Pol (1926) Pol oscillator was first introduced as a model of an electronic circuit and it has been studied for over 80 years. This gives us the possibility to anchor the present work in a long tradition of dynamical system theory. Next, the analysis techniques used with the van der Pol oscillator are applied to 6 other models previously published in the literature. We will then be able to determine which conclusions seem the most robust.

The van der Pol oscillator 2.1 Model definition
The van der Pol model can be introduced as a dynamical system of two coupled ordinary differential equations: with: (x, y) the climate state vector, τ a time constant, α a time-decoupling factor, β a bifurcation parameter and F (t) the forcing. The parameter β does not appear in the 5 original van der Pol equations. The present variant is sometimes referred to as the biased van der Pol model. The autonomous (i.e. F (t) = 0) model displays self-sustained oscillations as long as |β| < 1. For later reference the period of the unforced oscillator is denoted T n (τ ). The variable x follows then a saw-tooth periodic cycle, the asymmetry of which is controlled by β.
Variable y varies the more abruptly that α is high. Here we chose α = 30, so that y may be termed the 'fast' variable. In climate terms, x may be interpreted as a glaciation index, which accumulates slowly the effects of the astronomical forcing F (t), while y, which shifts between approximately −1 and 1, might be interpreted as some representation of the ocean dynamics. This is admittedly arguable and other interpretations would be possible. Keep in mind that the van der Pol model is only used here to identify phenomena emerging from the combination of limit cycle dynamics and astronomical forcing, and models with better physical justifications are analysed section 3.
In ice age models the forcing function is generally one or several insolation curves, computed for specific seasons and latitudes. The rationale behind this choice is that whichever insolation is used it is, to a very good approximation, a linear combination of climatic precession and obliquity (Loutre (1993), see also Appendix A). The choice of one specific insolation curve may be viewed as a modelling decision about the effective forcing phase of climatic precession, and the relative amplitudes of the forcings due to precession and obliquity. In turn, climatic precession and obliquity can be expressed as a sum of sines and cosines of various amplitudes and frequencies (Berger, 1978), so that F (t) can be modelled as a linear combination of about a dozen of dominant periodic signals, plus a series of smaller amplitude components. They are shown on Figure 2. More details are given in section 2.4. With these hypotheses the van der Pol model may be tuned so that the benthic curve over the last 800,000 years is reasonably well reproduced ( Figure 3).
An abundant literature analyses the response of the van der Pol oscillator to a periodic forcing (e.g. Mettin et al., 1993;Guckenheimer et al., 2003, and ref. therein). The response of oscillators to the sum of two periodic forcings has been the focus of attention because it leads to the emergence of 'strange non-chaotic attractors', on which we will come back (Wiggins, 1987;Romeiras and Ott, 1987;Wojewoda, 1990, 1993;Belogortsev, 1992;Feudel et al., 1997;Glendinning et al., 2000). To our knowledge, however, there is no systematic study of the response of an oscillator to a signal of the form of the astronomical forcing, except for preliminary work of our group (De Saedeleer et al., 2012). Le Treut and Ghil (1983), for example, represented the astronomical forcing as a sum of only two or three periodic components and it will be shown here that it matters to consider the astronomical forcing with all its complexity.

Periodic forcing
Consider a sine-wave forcing (F (t) = γ sin(2π/P 1 t + φ P 1 )), with period P 1 = 23, 716 years and φ P 1 = 32.01 • . This is the first component of the harmonic development of climatic precession (Berger, 1978). If certain conditions are met-they will soon be given-the van der Pol oscillator may become synchronised on the forcing. Synchronised means, in this particular context, that the response of the system displays p cycles within q forcing periods, where p and q are integers. It is said that the system is in a p : q synchronisation regime (Pikovski et al., 2001, p. 66-67). The output is periodic, and its period is equal to q × P . There are several ways to identify the synchronisation in the output of a dynamical system. One method is to plot the state of the system at a given time t, and then superimpose on that plot the state of the system at every time t + nP , where n is integer. The system is synchronised if only q distinct points appear on the graph, discarding transient effects associated with initial conditions. These correspond to the stable fixed points of the iteration bringing the system from t to t + qP . In the following we refer to this kind of plot as a "stroboscopic section" of period P (Figure 4a). If the system is not synchronised, then there are two options: the stroboscopic section is a closed curve (the response is quasi-periodic), or a figure with a strange geometry (the response is a-periodic).
There is another, equivalent way to identify synchronisation. Suppose that the system is started from arbitrary initial conditions. Then, plot the system state at a given time t, long enough after the initial conditions. Repeat the experiment with another set of initial conditions, superimpose the result to the plot, and so on with a very large number of different initial conditions. In doing so one constructs the section of the global pullback attractor at time t (henceforth referred to as the pullback section); the pullback attractor itself being the continuation of this figure over all times ( Figure 4b). Each component of the global pullback attractor (two aro illustrated on Figure 4b) is a local pullback attractor. Rasmussen (2000, chap. 2) reviews all the relevant mathematical formalism.
In the particular case of a periodic forcing, the stroboscopic section and the pullback section are often identical ( Figure 5). 1 The number of points on the pullback section may then be estimated for different combinations of parameters and we can use this as a criteria to detect synchronisation. This is done on Figure 6 for a range of γ and τ . It turns out that synchronisation regimes are organised in the form of triangles, known in the dynamical system literature as Arnol'd tongues (Pikovski et al., 2001, p. 52). A p : q synchronisation regime appears when the ratio between the natural period (T n ) and the forcing period is close to q/p. The tolerance, i.e., how distant this ratio can afford to be with respect to q/p, increases with the forcing amplitude and decreases with p and q. Synchronisation is weakest (least reliable) near the edge of the tongues. Unreliable synchronisation characterises a system that is synchronised, but in which small fluctuations may cause episodes of desynchronisation. In the particular case of periodic forcing the episode of desynchronisation is called a phase slip, as is well explained in Pikovski et al. (2001, p. 54) Using the pullback attractor to identify synchronisation is not a very efficient method in the periodic forcing case. Arc-length continuation methods are faster and more accurate (e.g., Schilder and Peckham, 2007, and ref. therein). It is shown however in De Saedeleer et al. (2012) that the pullback section method gives results that are acceptable enough for our purpose, and it is adopted here because it provides a more intuitive starting point to characterise synchronisation with multi-periodic forcings. Periodic forcing (P 1 ) Figure 5: Pullback (at t 0 = 0) and stroboscopic sections (t = 0 + nP 1 ) obtained with the van der Pol oscillator, with parameters α = 30, β = 0.7, τ = 36 forced by F (t) = γ sin(2π/P 1 t + φ P 1 ), P 1 = 23.7 ka, φ P 1 = 32.01 • and γ = 0.6. The two plots indicate a case of 4 : 1 synchronisation, and they are identical because the forcing is periodic.

Synchronisation on two periods
Consider now a forcing function that is the sum of two periodic signals. Two cases are considered here: the two forcing periods differ by a factor of about 2, and the two forcing periods are close.
2.3.1 P 1 = 23.716ka and O 1 = 41.000ka We adopt F (t) = γ[sin(2π/P 1 t + φ P 1 )) + cos(2π/O 1 t + φ O1 )], with P 1 = 23.716ka, O 1 = 41.000ka and φ P 1 = 32.01 • and φ O1 = 251.09 • . P 1 is the first period in the development of precession, O 1 is the first period in the development of obliquity and φ P 1 ,φ O1 the corresponding phases given by Berger (1978), so that F (t) may already be viewed as a very rough representation of the astronomical forcing. Let us begin with τ = 36 ka, which corresponds to a limit cycle in the van der Pol oscillator of period T n = 98.79 ka, and consider the stroboscopic section on P 1 (Figure 7, line 1). Forcing amplitude γ is set to 0.6. Due to the presence of the O 1 forcing, the four points of the periodic case seen on Figure 5 have mutated into four local attractors, which appear as closed curves (some are very flat). Every time P 1 elapses, the system visits a different local attractor. They are attractors in the sense that they attract solutions of the iteration bringing the system from t to t + 4 · P 1 . In this particular example, the system is said to be phase-or frequency-locked on P 1 with ratio 1:4 (Pikovski et al., 2001, p. 68), because on average, one ice age cycle takes four precession cycles, even though the response is no longer periodic. In this example, the curves on the P 1 stroboscopic section nearly touch each other. This implies that synchronisation is not reliable since a solution captured by one of these attractors could easily escape and fall into the basin of attraction of another local attractor. One can also see that it is not synchronised on O 1 since the stroboscopic section of period O 1 shows one closed curve englobing all possible phases. It may also be said that the system is synchronised in the generalised sense (Pikovski et al., 2001, p. 150), because the pullback section is made of only four points : starting from arbitrary conditions, the system converges to only a small number of solutions at any time t. It is also stable in the Lyapunov sense, a point that will not be further discussed here, but see De Saedeleer et al. (2012).
Consider now τ = 41 ka. The four closed curves on the P1-stroboscopic section have collided and merged into one attractor with strange geometry. A similar figure appears on the E1-stroboscopic section. The phenomenon of strange non-chaotic attractor has been described since Romeiras and Ott (1987), its occurrence in the van der Pol oscillator is discussed in Kapitaniak and Wojewoda (1990), and its relevance to climate dynamics was suggested by Sonechkin and Ivachtchenko (2001). In our specific example, the system is neither frequency-locked on P 1 nor on O 1 , but it is synchronised in the generalised sense: the pullback section has two points. Finally, with τ = 44 ka there is frequency-locking on O 1 (regime 3:1) but not on P 1 . Clearly, the system underwent changes in synchronisation regimes as τ increased from 36 to 44 ka. Further insight may be had by considering the τ − γ plot (Figure 8). The frequency locking regime on P 1 lies in the relic of the 1:4 tongue visible in the periodic forcing case (Figure 6). Frequency locking on O 1 belongs to the 1:3 tongue associated with O 1 . The strange non-chaotic regime occurs where the tongues associated with these different forcing components merge.
The word bifurcation has been defined for non-autonomous dynamical systems (Rasmussen, 2000, chap. 2). This is a complex subject and we will admit here the rather informal notion that there is a bifurcation when a local pullback attractor appears or ceases to exist (adapted from Def. 2.42, in Rasmussen, 2000). With this definition, there is a bifurcation at least every time color changes on Figure 8 (assuming t back is far enough in the past).
Another view on the bifurcation structure may be obtained by plotting the x and y solutions of the system at t 0 = 0, initiated from a grid of initial conditions at t back = −5 Ma, as a function of τ , still with γ = 0.6 ( Figure 9). This plot outlines a region of weak structural stability, where the system shows pretty strong dependence on the value of the parameters. It occurs here between τ = 37 and 42 ka. Quasi-periodic forcing (P 1 +E 1 ) Quasi-periodic forcing (P 1 and E 1 ) These observations have two important consequences for our understanding of the phenomena illustrated on Figure 1. To see this it is useful to refer to general considerations about autonomous dynamical systems. A bifurcation generally separates two distinct (technically: non-homeomorphic) attractors, which control the asymptotic dynamics of the system. As the bifurcation is being approached, the convergence to the attractor is slower, while the attractor that exists on the other side of the bifurcation may already take some temporary control on the transient dynamics of the system. This is, namely, one possible mechanism of excitable systems. One sometimes refers to 'remnant' or 'ghost attractors' to refer to these attractors that exist on the other side of the bifurcation and may take control on the dynamics of the system over significant time intervals (e.g. Nayfeh and Balachandran, 2004, p. 206) The idea may be generalised to non-autonomous systems. Consider Figure 10. The upper plot shows the two local pullback attractors of the system obtained with τ = 41 ka. The middle panel displays one local attractor obtained with τ = 40. The two τ = 41 attractors are reproduced with thin lines for comparison. Observe that this τ = 40 attractor is qualitatively similar to the τ = 41 attractors, and most of its time is spent on a path that is nearly undistinguishable from those obtained with τ = 41. However, on a portion of the time interval displayed it follows a sequence of ice ages that is distinct from those obtained with τ = 41. In fact, there are four pullback attractors at τ = 40, which clarifies the fact that there is at least one bifurcation between τ = 40 and τ = 41 ka.
Let us now consider a third scenario. Parameter τ = 41, but an additive stochastic term (σ dω dt , σ 2 = 0.25ka −1 , and ω symbolises a Wiener process) is added to the second system equation. This is thus a slightly noisy version of the original system. Shown here is one realisation of this stochastic equation, among the infinity of solutions that could be obtained with this system. Expectedly, the system spends large fractions of time near one or the other of the two pullback attractors. However, it also spends a significant time on a distinct path. Speculatively, this distinct path is under the influence of a 'ghost' pullback attractor. As the bifurcation structure is complex and dense, as shown on Figure 9, we expect a host of ghosts to lie around, ready to take control of the system over significant fractions of time, and this is what happens in this particular case.
To further support this hypothesis, consider a second experiment. Figure 11 displays the number of distinct solutions counted at time t 0 = 0, when the system is started from 121 distinct initial conditions at a time back t back , as a function of t back . Surprisingly, one needs to go back to −30 Ma (million years) to identify the true pullback attractor. Obviously 30 Ma is a very long time compared to the Pleistocene and so this solution is in practice no more or relevant than the 4 or 8 solutions that can be identified by only starting the system back in time 1 or 2 Ma ago. They may be interpreted as ghost ('almost alive') pullback attractors, and following the preceding discussion they are likely to be visited by a system forced by large enough random external fluctuations.  The two periods now being combined are the second and third components of precession, still according to Berger (1978). These two periods were selected for two reasons. The first one is that the addition of the two periodic signal produces an interference beating with period 123 319 yr, not too far away from the usual 100-ka cycle that characterises Late Pleistocene climatic cycles. Second, the period of the beating is not close to an integer number of the two original periods (this occurs, accidentally, when using P 1 and P 3 ). This was important to be able to clearly distinguish a synchronisation to the beating from a higher resonance harmonic to either forcing components.
It is known from astronomical theory that the periodicity of eccentricity is mechanically related to the beatings of the precession signal (Berger, 1978).The scientific question considered here is whether the correspondence between the period of ice age cycles and eccentricity is coincidental, or whether a phenomenon of synchronisation of climate on eccentricity developed.
To address this question we need a marker of synchronisation on the precession beating. The Rayleigh number has already been used to this end in palaeoclimate applications (Huybers and Wunsch, 2004;Lisiecki, 2010). Let P b be the beating period, and X i the system state snapshotted every t = t 0 + iP b , the Rayleigh number R is defined as | X i − X|/ |X i −X|, where the overbar denotes an average. R is strictly equal to 1 when the solution is synchronised with a periodic forcing of period P b , assuming no other source of fluctuations. As a reference, Lisiecki (2010) estimated 0.94 the Rayleigh number of a stacked benthic δ 18 O signal with respect to eccentricity over the last million years.
The bifurcation diagram showing the number of pullback solutions is displayed on Figure 12. The frequency-locking tongues on P 2 and P 3 are easily identified at low forcing amplitude; as forcing amplitude increases the tongues merge and produce generalised synchronisation regimes. Regions of synchronisation on P b , identified as R > 0.95, are hashed. They occur when the system natural period is close to P b but narrow bands also appear near T N = P b /2. Observe also that these synchronisation regimes are generally not unique (several pullback attractors co-exist), and additional sensitivity experiments show that convergence is quite slow. More specifically, the synchronisation diagram was computed here using t back = −10 Ma. With shorter backward time horizons, the number of remaining solutions identified in the beating-synchronisation regime often exceeds 6 and could not been seen on the graph, while the Rayleigh number was still beyond 0.95. Hence, a high Rayleigh number is not necessarily a good indicator of reliable synchronisation. 2-period forcing (P 2 and P 3 ) Figure 12: As Figure 8 but with : F (t) = γ(sin(2π/P 2 t + φ P 2 ) + sin(2π/P 3 t + φ P 3 )).
Hashes indicate regions of Rayleigh number > 0.95 on the beating associated with P 2 and P 3 , the period of which is called E 3 in the Berger (1978) nomenclature (third component of eccentricity).

Full astronomical forcing
The next step is to consider the full astronomical forcing, as the sum of standardized climatic precession (Π) and the deviation of obliquity with respect to its standard value (O): where The various coefficients are taken from Berger (1978). We take N p = N o = 34, so that the signal includes in total 68 harmonic components. With this choice the BER78 solution (Berger, 1978) is almost perfectly reproduced. BER78 is still used in many palaeoclimate applications. Compared to a state-of-the-art solution such as La04 (Laskar et al., 2004), the error on amplitude is between 0 and 25 %, and the error on phase is generally much less than 20 • .
The bifurcation diagram representing the number of pullback solutions as a function of forcing amplitude and τ is shown on Figure 13. We have taken γ = γ p = γ o . One recognises the tongues originating from the individual components merging gradually into a complex pattern. The number of attractors settles to 1 as the amplitude of the forcing is further increased. Let us call this the 1-pullback attractor regime. We already know that synchronisation is generally not reliable in the region characterised by the complex and dense bifurcation region, where more than one attractor exist. The remaining problem is to characterise the reliability of synchronisation in the 1-pullback attractor regime.
The literature says little about systematic approaches to quantify the reliability of generalised synchronisation with quasi-periodic forcings. To develop further the ideas developed in section 2.3.1, one can plot pullback solutions at a certain time t as a function of one or several parameters. This is done Figure 14. Here γ is kept constant (= 1.0) and τ is varied between 25 and 40. There is a brief episode of 2-solution regime between 31 and 32 ka. Within the 1-solution regime there is a number of abrupt transitions (at 26, 27, 34 and 37 ka).
As we have seen above, the density of bifurcations in the parameter space is an indicator of the structural stability of the system. Changes in the number of pullback attractors are clearly bifurcations. Abrupt variations such as near τ = 26, 27, 34 and 37 ka could well be bifurcations because one observes slower convergence near the transitions, with co-existence of two or more solutions when t back is only 1 Ma (not shown). It may also be verified that near these transition zones scenarios similar to those depicted on Figure 10 may be observed. Most of the transition zones disappear as γ is further increased. Hence, being in the 1-pullback attractor zone is not quite enough to guarantee a reliable synchronisation: one need to be deep into that zone.

Other models
We now consider 6 previously published models. Mathematical details are given in the Appendix and the codes are available on-line at https://github.com/mcrucifix.

SM90
: This is a model with three ordinary differential equations representing the dynamics of ice volume, carbon dioxide concentration and deep-ocean temperature. The astronomical forcing is linearly introduced in the ice volume equation, under the form of insolation at 65 • North on the day of summer solstice. Only the carbon dioxide equation is non-linear, and this non-linearity induces the existence of a limit-cycle solution-spontaneous glaciation and deglaciation-in the corresponding autonomous system. The SM90 model is thus a mathematical transcription of the hypothesis according to which the origin 100,000 year cycle is to be found in the biological components of Earth's climate.

SM91
: This model is identical to SM90 except for a difference in the carbon cycle equation.

PP04:
The Paillard-Parrenin model (Paillard and Parrenin, 2004) is also a 3-differentialequation system, featuring Northern Hemisphere ice volume, Antarctic ice area and carbon dioxide concentration. The carbon dioxide equation includes one non-linear term associated to a switch on/off of the southern ocean ventilation. Astronomical forcing is injected linearly at three places in the model: in the ice-volume equation, in the carbon dioxide equation, and in the ocean ventilation parameterisation. The autonomous version of the model also features a limit cycle. As in SM90 and SM91 the non-linearity introduced in the carbon cycle equation plays a key role but the bifurcation structure of this model differs from SM90 and SM91 (Crucifix, 2012).

T06:
The Tziperman et al. (2006) model is a mathematical idealisation of more complex versions previously published by Gildor and Tziperman (2000). T06 features the concept of sea-ice switch, according to which sea-ice growth in the Northern Hemisphere inhibits accumulation of snow over the ice sheets, and vice-versa. Mathematically, T06 is presented as a hybrid model, which is the combination of a differential equation in which the astronomical forcing is introduced linearly as a summer insolation forcing term, and a discrete variable, which may be 0 or 1 to represent the absence or presence of sea-ice in the northern hemisphere.

I11:
The Imbrie et al. (2011) was introduced by its authors as a "phase-space" model. It is a 2-D model, of which the equations were designed to distinguish an 'ice accumulation phase' and an 'abrupt deglaciation' phase, which is triggered when a threshold defined in the phase space is crossed. I11 was specifically tuned to reproduce the phase-space characteristics of the benthic oxygen isotopic dynamics. A particularity of this model is that the phasing and amplitude of the forcing depend on the level of glaciation.
PP12: Similar to Imbrie et al. (2011), the Parrenin and Paillard (2012) model distinguishes accumulation and deglaciation phases. Accumulation is a linear accumulation of insolation, without restoring force (hence similar to equation 1 of the van der Pol oscillator); deglaciation accumulates insolation forcing but a negative relaxation towards deglaciation is added. Contrarily to Imbrie et al. (2011), the trigger function, which determines the regime change, is mainly a function of astronomical parameters. An ice volume term only appears in the function controlling the shift from 'accumulation' to 'deglaciation' regime.
The ice volumes (or, equivalently, glaciation index or sea-level) simulated by each of these models are shown on Figure 15. Shown here are estimates of the pullback attractors. More specifically, the trajectories obtained with an ensemble of initial conditions at t back = −20 Ma; in some cases the curves actually published (in particular in SM90) are not pullback attractors, but ghost trajectories in the sense illustrated on Figure 11. In some cases (PP04 and PP12) the parameters had to be slightly adjusted to reproduce the published version satsifactorily. Details are given in appendix.
Some of these models include as much as 14 adjustable parameters (e.g.: PP04) and a full dynamical investigation of each of them is beyond the scope of the study. Rather, we proceeded as follows. Every model responds to a state equation, which may be written, in general (assuming a numerical implementation), as: where t i is the discretized time, x i the climate state (a 2-D or 3-D vector) at t i and F (t) is the astronomical forcing, which is specific to each model because the different authors made different choices about the respective weights and phases of precession and obliquity.
In all generality, the equation (or its numerical approximation) may be rewritten as follows, posing τ = 1 and γ = 1: The parameters τ and γ introduced this way have a similar meaning as in the van der Pol oscillator, since τ controls the characteristic response time of the model, while γ controls the forcing amplitude.
Bifurcation diagrams, similar to Figure 13 are then shown on Figure 16. Remember that γ = τ = 1 corresponds to the model as originally published.
A first group of four models appears (SM90, SM91, T06 and PP04), on which one recognises a similar tongue-synchronisation structure as in the van der Pol oscillator. This was expected since these models are also oscillators with additive astronomical forcing. Depending on parameter choices synchronisation may be reliable or not. Synchronisation is clearly not reliable in the standard parameters used for SM90 and SM91. In T06 the standard parameters are not far away from the complex multi-pullback-attractor regime and this explains why transitions such as displayed on Figure 1 could be obtained under small parameter changes (or, equivalently, with some noise, as discussed in the original Tziperman et al. (2006) study). The standard parameter choice of PP04 is further into the stability zone and, indeed, experimenting with this model shows that instabilities such as displayed on Figure 1 are harder to obtain. I11, with standard parameters, is also fairly deep in the stable zone. One has to consider much smaller forcing values than published to recognize the synchronisation tongue structure that characterises oscillators ( Figure 17).
PP12 finally turns out to be the only case not showing a tongue-like structure. This may be surprising because this model also has limit cycle dynamics (self-sustained oscillations in absence of astronomical forcing). Some aspects of its design resemble the van der Pol oscillator. The role of the variable y in the van der Pol oscillator is here played by the mode, which may either be g (glaciation) or d (deglaciation). Also similar to the van der Pol, the direct effect of the astronomical forcing on the ice volume (x in the van der Pol; v in PP04) is additive. T06 has also similar characteristics. The distinctive feature of PP12 is that the transition between deglaciation and glaciation modes is controlled by the astronomical forcing and not by the system state, as in the other models. To assess the importance of this element of design we considered a hacked version of PP12, where the d → g occurs when for v < v 1 (cf. appendix B.6 for model details). In this case the tongue synchronisation structure is recovered, with standard parameters marginally in the reliable synchronisation regime (Figure 17).

Conclusion
The present article is built around the paradigm of the 'pacemaker', that is, the timing of ice ages arises as a combination of climate's internal dynamics with the variations of incoming solar radiation induced by the variations of our planet's orbit and obliquity. This is not the only explanation of ice ages, but this is certainly one of the most plausible.
In this study we paid attention to the dynamical aspects that may affect the stability of the ice age sequence and its predictability. First, the astronomical forcing has a rich harmonic structure. We showed that a system like the van der Pol oscillator is more likely to be synchronised on the astronomical forcing as Nature provides it than on a periodic forcing, because the fraction of the parameter space corresponding to synchronisation is larger in the former case. A synchronised system is Lyapunov stable, so that at face value this would imply that the sequence of ice ages is stable. However, -this is the second point-even if the dynamical structure of the Pleistocene climate was correctly identified, there would be at least two sources of uncertainties : random fluctuations associated with the chaotic atmosphere and ocean and other statistically random forcings like volcanoes; and uncertainty on system parameters. In theory both types of uncertainties point to different mathematical concepts: path-wise stability to random fluctuations in the first case, and the structural stability in the second. In practice, however, lack of either form of stability will result in similar consequences: quantum skips of insolation cycles in the succession of ice ages. This was the lesson of Figure 10.
It was shown here that, compared to periodic forcing, the richness of the harmonic structure of astronomical forcing favours situations of weak structural stability. To preserve stability, the richness of the astronomical forcing has to be compensated for by large enough forcing amplitude.
Out of the seven models tested here, we ignore which one best captures ice ages dynamics. The overwhelming complexity of the climate system does not allow us to securely select the most plausible model on the sole basis of our knowledge of physics, biology and chemistry. Consequently, while we have understood here how and why the sequence of ice ages could be unstable in spite of available evidence (astronomical spectral signature; Rayleigh number), estimating the stability of the sequence ice ages and quantifying our ability to predict ice ages is also a problem of statistical inference : calibrating and selecting stochastic dynamical systems based on both theory and observations, which are sparse and characterised by chronological uncertainties. A conclusive demonstration of our ability to reach this objective is still awaited.
whereẏ := dy dt andÿ := dẏ dt . Note that the polynomial on the right-hand-side of the equation forẏ is a continuous fit to the piece-wise function used in the original Imbrie et al. (2011) publication. Time units are here ka.

B.6 PP12 model
This is a hybrid dynamical system, with ice volume v (expressed in m of equivalent sea-level) and state, which may be g (glaciation) or d (deglaciation).
Define first with a = 0.68034. Then define the following quantities, standardized as follows: Π = (f (Π) − 0.148)/0.808 Π = (f (Π) − 0.148)/0.808 the threshold θ = k ΠΠ + kΠΠ + k OŌ , and finally the following rule controlling the transition between state g and d: Ice volume v, expressed in sea-level equivalent, responds to the following equation: with the following parameter values: a Π = 1.456m/ka, aΠ = 0.387m/ka, a O = 1.137m/ka, a g = 0.978 m/ka, a d = −0.747m/ka, τ = 0.834ka, k Π = 14.635m, kΠ = 2.281m, k O = 23.5162m, v 0 = 122.918m and v 1 = 3.1031m, assuming that one time unit= 10ka. This parameter set is the one originally presented by Parrenin and Paillard in Climate of the Past Discussion (which differs from the final version in Climate of the Past), except that k O is 18.5162m in the original paper. This modification was needed to reproduce the exact sequence of terminations shown by the authors. Subtle details, such as the numerical scheme or the choice of the astronomical solution might explain the difference. All codes and scripts are available from GitHub at https://github.com/mcrucifix.