Is the astronomical forcing a reliable and unique pacemaker for climate? A conceptual model study

There is evidence that ice age cycles are paced by astronomical forcing, suggesting some kind of synchronisation phenomenon. Here, we identify the type of such synchronisation and explore systematically its uniqueness and robustness using a simple paleoclimate model akin to the van der Pol relaxation oscillator and dynamical system theory. As the insolation is quite a complex quasiperiodic signal involving different frequencies, the traditional concepts used to define synchronisation to periodic forcing are no longer applicable. Instead, we explore a different concept of generalised synchronisation in terms of (coexisting) synchronised solutions for the forced system, their basins of attraction and instabilities. We propose a clustering technique to compute the number of synchronised solutions, each of which corresponds to a different paleoclimate history. In this way, we uncover multistable synchronisation (reminiscent of phase- or frequency-locking to individual periodic components of astronomical forcing) at low forcing strength, and monostable or unique synchronisation at stronger forcing. In the multistable regime, different initial conditions may lead to different paleoclimate histories. To study their robustness, we analyse Lyapunov exponents that quantify the rate of convergence towards each synchronised solution (local stability), and basins of attraction that indicate critical levels of external perturbations (global stability). We find that even though synchronised solutions are stable on a long term, there exist short episodes of desynchronisation where nearby climate trajectories diverge temporarily (for about 50 kyr). (...)


Introduction
This article is a contribution to the field of paleoclimate dynamics theory, which has experienced many developments in terms of ice age models since many years notably by [Le Treut andGhil, 1983, Saltzman andMaasch, 1990], and many others, and remains an active research field. Paleoclimate modeling is a complex problem, hence an uncomfortable situation for a scientist. On one hand, the data are scarce and marred by uncertainties. On the other hand, there is not a single well established model, the problem is non autonomous, the forcing is aperiodic, and stochastic ef-  [Lisiecki and Raymo, 2005] [Luethi et al., 2008]. High values of CO 2 correspond to a warmer climate (interglacial state). Figure 1: The long-term climatic signals reveal a slow-fast dynamics (only the most recent 500 kyr of the data are displayed here). The areas in grey (width τ slow ) are wider than the areas in white (width τ f ast ): while glaciation is a slow process of ice build-up, deglaciation occurs much more rapidly (τ slow > τ f ast ). One also recognizes the last deglaciation which started some 20 kyr ago, up to the present time (t = 0). fects are present.
Here, we focus on the slow variations of climate over the few last million years, which include the phenomenon of ice ages [Hays et al., 1976], that is, the repeated growth and decay of ice sheets in the Northern Hemisphere of a total mass as big as modern Antarctica's. When examining long-term climatic signals like the 5.3 Myr-long stack produced in [Lisiecki and Raymo, 2005], or the 800 kyr-long EPICA Dome C Ice Core from [Luethi et al., 2008], plotted respectively in Fig. 1(a) and Fig. 1(b) for the last 500 kyr, one immediately identifies three clearly visible features of the climatic time series: (i) oscillations: the signal oscillates between higher and lower values of ice volume corresponding to the glacial and interglacial states, (ii) asymmetry: in Fig. 1(a) typical transitions from a minimum to a maximum take much longer than transitions from a maximum to a minimum: deglaciations occur much more rapidly (τ f ast ≈ 10 kyr) than glaciations (τ slow ≈ 80 kyr), giving a distinctive saw-tooth structure in the glacial/interglacial (G/I) cycles, especially pronounced over the last 500 kyrs, (iii) 100-kyr dominant period: this has been identified by many authors since [Broecker and van Donk, 1970].
Note that the G/I cycles are not periodic.
The asymmetry in the oscillations has been studied by many authors. In order to reproduce it, some authors use underlying physical principles to build phenomenological models that exhibit slow-fast dynamics reasonably mimicking the climatic proxies [Saltzman, 2002]. Others assume this asymmetry by explicitly defining 2 different parameters such as the time intervals τ up = τ slow and τ down = τ f ast [Ashkenazy, 2006] or time constants (τ R and τ F in [Paillard, 1998] and T w and T c in [Imbrie and Imbrie, 1980]). Whatever the model, it has to ultimately exhibit asymmetric oscillations under the effect of the forcing, as it is aimed to mimic the oscillations between G/I states. Relaxation oscillators are therefore very straightforward natural candidates of ice age models. In this article, we will consider a slightly modified van der Pol oscillator model to illustrate the new contributions of our synchronization concepts.
In this paper, we concentrate on the influence of the astronomical forcing on Earth's climate. This forcing is induced by the slow variations in the spatial and seasonal distributions of incoming solar radiation (insola-tion) at the top of the atmosphere, associated with the slow variations of the Earth's astronomical elements: eccentricity (e), true solar longitude of the perihelion measured with respect to the moving vernal equinox (ϖ), and Earth obliquity (ε E ). These quantities are now accurately known over several tens of millions of years [Laskar et al., 2004], but analytical approximations of e, e sin ϖ, and ε E valid back to one million years have been known since [Berger, 1978]. They take the form of d'Alembert series (∑ A i sin[ω i t + φ i ]).
The external forcing F(t) used throughout this article is the insolation at 65 • N latitude on the day of the summer solstice. That specific insolation quantity is commonly related to the Milankovitch theory and can be thought of as a measure of how much ice may melt over summer. It can be written under the following compact form: (1) where the value of the 3 × 35 parameters (including ω i , s i , and c i ) are given in the Table 1 of Appendix A. The coefficients were extracted from [Berger, 1978] by performing a linear regression of the insolation on the ω i . The validity range of this approximation is [-1 Myr, 0 Myr], and its mean error (mean of the absolute value of the difference, compared to [Laskar et al., 2004]) is 6.7 W/m 2 with peaks at 27.5 W/m 2 . Note that the mean value (494.2447 W/m 2 ) has been removed; the theoretical framework that allows to work with anomalies was justified by [Saltzman and Maasch, 1991]. In short, this can be done as we are interested in oscillations, and not in the mean values themselves. The quasiperiodic 1 nature of the insolation forcing F(t) is illustrated in its spectrum decomposition in Fig. 2. Precession is dominated by two harmonics around 19 and 23 kyrs (1 kyr = 1,000 years) and obliquity is dominated by an harmonic with a period of 41 kyrs but it bears periods as long as 1,200 kyrs.

Synchronization
There is ample evidence that the astronomical forcing influences the climate system. The phrase 'pacemaker of ice ages' was coined in a seminal paper [Hays et al., 1976] to express the idea that the timing of ice ages is controlled by the astronomical forcing, 1 A quasiperiodic signal is the superposition of several periodic signals with uncommensurate periods.  1)). This is a graphical representation of Table 1 in Appendix A, with a 2 i = s 2 i + c 2 i and T i = 2π/ω i . The eight largest parameters a i , which represent already 80% of the signal, are the major components of the insolation; they clearly come from the precession (19 and 23 kyr), and from the obliquity (41 kyr) associated series. This insolation is the forcing used to construct all Figures subject to the quasiperiodic astronomical forcing. The main harmonic ε 1 associated with obliquity has an angular velocity ω ε 1 = 0.1532 rad/kyr (T ε 1 = 41.0 kyr) and an amplitude of a ε 1 = 11.77 W/m 2 . The three main harmonics associated with precession are denoted p 1 , p 2 and p 3 . while the ice age cycle itself is shaped by internal system dynamics. The paradigm has prevailed since then and is it still supported by the most recent analyses of palaeoclimate records Raymo, 2007, Huybers, 2007]. The notion of 'pacemaker' naturally evokes some sort of synchronization. However, despite some attempts, the actual type of synchronization has not been clearly identified or demonstrated to date.
For example, [Ashkenazy, 2006, Tziperman et al., 2006 speak of "nonlinear phaselocking" although they do not define suitable "phase variables" that can be used to demonstrate a fixed-intime relationship between phases of the forcing and the oscillator response.
Synchronization, as a universal nonlinear phenomenon, is a pervasive process in Nature, as it is associated with rhythmic processes. It is therefore not surprising to have synchronization also in Paleoclimatic Sciences. Depending on the forcing type (periodic, chaotic, stochastic), one can distinguish many types of synchronization including complete, lag, phase, frequency, identical, generalized [Rulkov et al., 1995], achronal and isochronous [Wu et al., 2006], and even noise synchronization. For a review, the reader is referred to [Balanov et al., 2009, Pikovsky et al., 2001, among others. While some terminology is still debated, [Brown and Kocarev, 2000] proposed an unified definition of synchronization for dynamical systems-there is synchronization if there exists a relationship h between the measured properties of the forcing, g( u), and those of the oscillator, g( v): that is fixed-in-time, meaning that h is time independent. Because we are interested in synchronization that is stable, for arbitrary initial conditions u(0) and v(0) that do not satisfy (2), we require that [Brown and Kocarev, 2000]: For example, if g( u) = u, g( v) = v, u and v have the same dimension, and (2) can be written as u = v, we speak of identical synchronization. More generally, if vectors u and v have different dimensions and (2) cannot be reduced to more than a functional relationship u = H( v), we speak of generalized synchronization; see also [Abarbanel et al., 1996, Rulkov et al., 1995, Pikovsky et al., 2001. Note that the relationship (2) need not be unique. If there are two or more relationships (2) for the same parameter settings, we speak of multistable synchronization [Pikovsky et al., 2001, Ch.15.3.2]. Then, which of the relationships the system settles to will depend on initial conditions.
In this paper, we use a simple van der Pol oscillator model to identify and illustrate for the first time the phenomenon of generalized synchronization between ice age cycles and astronomical forcing. The dynamical systems approach outlined in the next section (i) allows for stability analysis of such synchronization, (ii) uncovers interesting effects relating to the robustness of the synchronization with respect to external perturbations, and (iii) uncovers the phenomenon of multistable synchronization that has been overlooked by previous studies. We show that, in contrast to claims in [Tziperman et al., 2006], synchronization needs not be unique.
The article is structured as follows. Section 2 introduces a slightly modified version of the van der Pol oscillator as a suitable model for studying synchronization of ice ages to astronomical forcing. In Section 3, we analyse synchronization to periodic forcing and quasiperiodic astronomical forcing in terms of largest Lyapunov exponents. Section 4 is dedicated to the study of multistable synchronization in terms of attracting trajectories in the phase space of the forced system, and the associated basins of attraction. In Section 5, we investigate effects of the symmetry-breaking parameter β for the van der Pol oscillator model. Section 6 is concerned with the robustness of the synchronization and focuses on two aspects relating to predictability. Firstly, it shows that the local stability can be lost temporarily causing divergence of nearby climatic trajectories. Secondly, it demonstrates that in the multistable regime external perturbations (such as noise) may cause jumps between coexisting synchronized solutions when these solutions come close to their basin boundary. To be clear, all the treatment below is deterministic, except for Figs. 14 and 16.
This article requires some basics of Dynamical Systems theory (dynamical systems, nonlinear oscillations, limit cycles, bifurcations of vector fields, etc.), for which we refer the reader to [Guckenheimer and Holmes, 1983, Arnold, 1983, Strogatz, 1994, to [Saltzman, 2002] for dynamical paleoclimatology, and also to [Savi, 2005]  The hypothesis at the basis of the work by Milankovitch [Milankovitch, 1941] is that changes in total amount of continental ice (say: x) are driven by summer insolation F(t) already described in Eq. (1). One straightforward interpretation of this hypothesis is a simple differential where dΨ(x)/dx is the derivative of a climatic potential and γ is the forcing efficiency. However, models of this form fail in practice to correctly capture the rapid deglaciation phenomenon.
We therefore propose to model the paleoclimatic dynamical system with a dissipative self-sustained oscillator resembling the classical van der Pol oscillator 2 : where Φ (y) = y 3 /3 − y. Note that this system is nonautonomous because the right-hand side depends explicitly on time.
The physical interpretation of the model is as follows. Ice volume x integrates the external forcing F(t) over time but with a drift y + β . Assuming α 1, y is the faster variable whose dynamics is controlled by a two-well potential Φ(y). For example, there are arguments that the dynamics of the Atlantic ocean circulation may be approximated by an equation similar to Eq. 4b [Rahmstorf et al., 2005, Dijkstra et al., 2003.
Further interpretation and discussion of the fast variable can be found in [Saltzman et al., 1984, Tziperman and Gildor, 2003, Paillard and Parrenin, 2004, Tziperman et al., 2006, Crucifix, 2011. The parameter τ sets the slow time scale. The coupled system Eq. (4) has one stable equilibrium solution for |β | > 1 and a stable periodic orbit for |β | < 1. The ratio of time spent near the two stable branches of the slow manifold given by Φ (y) = x depends on β (see the paragraph "Time Spent" in Appendix B). We use T ULC to denote the period of the stable periodic orbit and ω ULC = 2π/T ULC to denote the corresponding angular velocity.
Relaxation oscillators have been proposed previously to study ice ages [Saltzman et al., 1984, Tziperman and Gildor, 2003, Paillard and Parrenin, 2004 although, to our knowledge, in a less general form than here. We adopted this form 3 because it is very close to the well-studied van der Pol oscillator, and a good agreement (timing of glaciations and deglaciations, and their amplitude) with ice volume proxies was easily found for well chosen values of α, β , γ and τ (Fig. 3). We note, though, that small changes in parameters or additive fluctuations may easily shift the timings of ice-age terminations for reasons that will be clarified later in the paper. The definition of synchronization can be applied to our model Eq. (4) as follows. The astronomical forcing F(t) corresponds to u(t), and the state vector whose two components are the slowly-varying ice volume x and the faster variable y corresponds to v(t). For nonperiodic forcing, relationship (2) can be very complicated (non-functional or even fractal-like) and hence difficult to detect. Therefore, other methods of detecting (2) had to be developed. As suggested by the auxiliary system approach [Abarbanel et al., 1996], relationships (2) and (3) are implied by an (invariant) attracting trajectory in the (x, y,t) phase space of the nonautonomous forced system (4) [Wieczorek, 2011]. In the remainder of the paper, such an attracting trajectory is denoted with AT and referred to as an attracting climatic trajectory or synchronized solution. All other solutions to Eq. (4) will be referred to as climatic trajectories.
They uncovered interesting dynamics including Arnol'd or mode-locked tongues consisting of 'interlocking' bubbles and open regions of multistability, nonsmooth bifurcations, and strange nonchaotic attractors. Here, we consider quasiperiodic forcing with 35 frequency components and focus on the regions of mode locking. Our approach is based on instabilities of attracting trajectories in the (x, y,t) phase space of the continoustime forced system because they relate directly to the concept of generalized synchronization. We can provide a systematic study of generalized synchronization to astronomical forcing by demonstrating existence of such trajectories and exploring their local and global stability properties. More specifically, we perform three types of calculations. Firstly, a clustering detection q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −2 −1 0 1 2 Y X −4 0 2 4 Forcing Figure 3: Top: the insolation forcing F(t), scaled by a ε 1 in order to work dimensionless. Bottom: the x and y climatic trajectories obtained using system Eq. (4) with α = 11.11, β = 0.25, γ = 0.75 and τ = 35.09. With these parameters, ω ε1 = 2.5ω ULC , where ω ε1 is the angular velocity associated with the dominant harmonic of obliquity and ω ULC is the angular velocity associated with the unforced system's periodic orbit. Blue dots correspond to the Lisiecki and Raymo stack (LR04) described in Fig. 1(a). Time t = 0 corresponds conventionally to the year 1950. The model is in good agreement with the ice volume proxy.
technique uncovers parameter regions with monostable (unique) and multistable (non-unique) synchronization. Secondly, the largest Lyapunov exponent along AT quantifies its long-and short-term local (linear) stability. Thirdly, a basin of attraction of AT quantifies its global (nonlinear) stability. Finally, we remark that in the theory of nonautonomous dynamical systems, attracting trajectories in the (x, y,t) phase space are linked to a modern and more general concept of a pullback attractor [Kloeden, 2000, Langa et al., 2002, Kosmidis and Pakdaman, 2003, Wiggins, 2003].
3 Synchronization of the paleoclimatic system to the insolation forcing

Illustration of the synchronization phenomenon
A typical climatic trajectory for γ = 0 is shown in Fig. 4, from two different points of view: the time series and the phase space portrait. In the time series, we recognize the slow variable x (the ice volume), while y exhibits slow-fast dynamics. This climatic trajectory is also shown in the two-dimensional (x, y) phase space of the autonomous system where arrows indicate direction of the flow. The trajectory converges to the limit cycle with slow-fast dynamics (the speed along the trajectory can be visually assessed by the circles of the residence plot). Let us now consider a set of 70 random initial conditions in the (x, y)-plane at time t 0 = 0, and study the resulting climatic trajectories in the threedimensional phase space (x, y,t) of the nonautonomous system ( Fig. 5(a)) for time t > t 0 . One clearly sees that all trajectories converge to a cylinder-the attracting set in the (x, y,t) space.
However, if we consider now an external forcing (γ > 0) then synchronization onto this forcing may occur under certain conditions [Ashkenazy, 2006, Tziperman et al., 2006. According to our definition (2-3), synchronization is represented by an attracting climatic trajectory in the (x, y,t) phase space.
Consider first the case of a purely periodic forcing with a period of T F = 41 kyr and strength γ = 3.33. The 70 initial conditions give rise to climatic trajectories that, after a sufficiently long integration time, converge to two attracting trajectories (see Fig. 5   other. This phenomenon is described in the literature as 2 : 1 phase-locking or frequency-locking. Generally speaking a n : m f requency-locking is defined as a fixed-in-time relation between the frequencies of the forcing (ω F ) and the oscillator response (ω R ) of the form n ω R = m ω F , where m and n are integers [Pikovsky et al., 2001, p.52].
Then consider the case of the quasiperiodic insolation forcing described in Eq. (1) with γ = 0.75 and τ = 43.86. Figure 5(e) shows that the 70 climatic trajectories now converge onto three attracting trajectories, which reveals that synchronization can be multistable [Pikovsky et al., 2001, p.348], [Balanov et al., 2009, p.94]. This phenomenon is described in the literature as mode-locking [Svensson and Coombes, 2009]. Note that because of the quasiperiodicity of the insolation forcing, these attracting trajectories are no longer periodic nor time-shifted versions of each other. The number of attracting trajectories depends on many factors including the dynamics of the unforced system, the nature of the forcing F(t), and the amplitude γ of the forcing. We will study this in more details in Sec. 4.

Detection of synchronization by the way of the largest Lyapunov exponent (LLE or λ max )
Local or linear stability of an attracting climatic trajectory can be quantified with the largest Lyapunov exponent (LLE) denoted here as λ max [Benettin et al., 1980]. The quantity λ max is a measure of the (average) exponential rate of divergence (λ max > 0) or convergence (λ max < 0) of nearby climatic trajectories. Therefore, a negative value of λ max indicates a locally attracting climatic trajectory or generalized synchronization [Pikovsky et al., 2001, Wieczorek, 2009. A transition from λ max < 0 to λ max = 0 indicates a bifurcation where the attracting climatic trajectory disappears and generalized synchronization is lost. Null and positive values of λ max indicate lack of synchrony (positive λ max indicates chaos but this regime is not encountered here). In the case of periodic forcing, computations of λ max can be easily validated with more precise and reliable numerical bifurcation continuation techniques (see § '41 kyr periodic forcing' below).  The largest Lyapunov exponent λ max is mathematically defined 4 as [Ott, 2002] :

Long-term
where δ Z = [δ x, δ y] are vanishing perturbations about x and y, respectively, governed by the linearization of system Eq. (4). Whereas this classical λ max is defined in long term limit (t → ∞), one can also define [Abarbanel et al., 1991] a short-term version, λ H max , by considering a finite time interval H (H = 50 kyr will be considered in this article): While λ max gives the average or long-term stability information, λ H max can tell us about the behaviour of nearby trajectories within a short time interval H. For example, λ max < 0 does not necessarily imply λ H max < 0 for some suitably chosen H. The definition (6) will be useful in studying the robustness of generalized synchronization in Sec. 6.
For computing λ max , as the system Eq. (4) and the Jacobian have an analytical form, tangent space methods [Kantz and Schreiber, 2004] can be used; technical details are given in the Appendix C.

Influence of the parameters γ and T ULC
The two particular types of synchronization illustrated in Fig. 5(c) and 5(e) have been obtained for a fixed value of the amplitude γ of the external forcing and of the natural period of the unforced paleoclimatic system T ULC . Now, we are equipped to achieve a much broader view of the dynamics by performing a parametric study on these two parameters. The quantity plotted in Figs. 6(a) and 6(b) is the largest Lyapunov Exponent at 3 Myr, far from the transient behaviour so that λ H=3 Myr max may already been considered as a good approximation of λ max . 41 kyr periodic forcing 4 Even if differential versions of the LLE have sometimes been developed mainly for computational efficiency purposes, we however preferred within this article to stick on the original definition of the LLE, because it is more standard and there is no insistent need for lowering computation time in the present framework, as the number of degrees of freedom of the system is reduced. Fig. 6(a) corresponds to the case of the 41 kyr periodic forcing (T F = 41 kyr). The synchronization region (λ max < 0) is composed of several V-shape regions, called Arnol'd tongues (phase-or frequency-locking), originating at 1, 2, 3, etc. times the forcing period T F . These regions correspond to 1:1, 2:1, 3:1 frequencylocking zones (3:2 and 5:2 can also be guessed). Periodic solutions are found within these regions which originate generally speaking at T ULC = (m/n) T F . No synchronization is possible when γ is zero but synchronization may occur already for infinitesimally small γ. Then, for increasing γ, the synchronization region widens and synchronization becomes more stable up to an optimum value γ * of the forcing. When γ > γ * , the synchronization becomes less and less effective, because at large γ the system is too much steered away from its natural dynamics; it may even be driven into chaos at yet higher forcing amplitude [Mettin et al., 1993] but this case is beyond our focus.
In order to perform an accurate validation of the synchronization region given by the LLE (λ max < 0) method, we computed the main Arnol'd tongues boundaries with the more accurate numerical continuation methods such as AUTO [Doedel et al., 2009]. The case of periodic forcing [F(t) = sin(ωt)] with β = 0 has been already extensively studied in the literature, analytically assuming some approximations [Guckenheimer and Holmes, 1983, p.70-75], and using numerical algorithms for pseudo-arc length continuation [Mettin et al., 1993]. The usual approach extends the original nonautonomous system by additional differential equations for the forcing so that the system becomes autonomous, and then explores the (ω, γ) parameter space. Note that the asymmetry introduced here with the parameter β adds slightly more complexity and induces additional features to the diagrams documented in these papers. In this way, we computed Arnol'd tongue boundaries as saddle-node of limit cycle bifurcations for the extended system with α = 11.11 and β = 0.25.
Superposition of LLE calculations and bifurcation boundaries in Fig. 6(a) shows that the synchronization regions obtained with two different techniques match perfectly. This is a confirmation that the method based on the LLE works fine and we will be able to use it for the case of the quasiperiodic insolation forcing. Note that bifurcation boundaries are also drawn in Fig. 6(c) in order to stress the correspondence with yet another method of detecting synchronization that will be discussed in Sec. 4.

Astronomical quasiperiodic forcing
For the case of the quasiperiodic insolation forcing ( Fig. 6(b)), the region of synchronization appears to be in one single piece with some indications of wellseparated tongues (mode locking) at small γ. In other words, whatever the value of the natural period T ULC of the paleoclimatic system, it has a higher probability of being synchronized onto the insolation forcing. Of course, for very low values of γ, there is still no synchronization. Note that γ has been downscaled by σ insol ≈ 5 (compared to σ sine = √ 2/2 for the 41 kyr periodic forcing), in order to keep a realistic range comparison.

Non uniqueness: multistability and basins of attraction
The detection of synchronization using the LLE (λ max < 0) gives only an Yes/No-type of information, without making any distinction between different tongues as this would require information about multistability. For example, Fig. 6(b) indicates synchronization for the parameter settings marked with the symbol '×" but gives no information about the number of attracting trajectories (we know that there are three different attracting trajectories from Fig. 5(e)). To explore the problem of multistable synchronization, we propose a clustering method that not only allows us to detect synchronization, but additionally provides information about the number of attracting trajectories denoted here with N.
Multistability Analysis: numerical estimate of the number of attracting trajectories N by a clustering technique Consider the case of the quasiperiodic insolation forcing with the three attracting trajectories, i.e., N = 3 (Fig. 5(e)). Although N can often be easily assessed visually, we want to automatically detect and count the number of AT s. As a matter of fact, N can be easily estimated in the following way. Fix a time t that defines a two-dimensional (x, y)-section in the (x, y,t) phase space. Then start with a grid of initial conditions at some time t 0 < t and take t − t 0 sufficiently large so that all the initial conditions converge to the attracting trajectories at t. Since each AT is represented by a point on the (x, y)-section, the problem of counting attract-ing trajectories reduces to a simple clustering problem. We designed a suitable automatic cluster detection algorithm that counts the number of clusters to obtain an estimate of N. For example, Fig. 5(f) shows 5 the (x, y)section of the three-dimensional (x, y,t) phase space at t = 550 kyr, given 70 initial conditions at t 0 = 0. The 70 trajectories converge onto three (highly concentrated) clusters corresponding to the three attracting trajectories.
The idea of using clustering analysis for paleoclimatic dynamics comes from the natural fact that clustering is another way of looking at generalized synchronization where negative LLE makes the trajectories cluster more efficiently. This provides another insightful viewpoint on the problem of identification of number of synchronized solutions of the paleoclimatic system: the more stable the synchronization, the more efficient formation of clusters.
Two important aspects of cluster analysis 6 have to be considered to avoid risks of mis-identification of clusters: • the notion of a cluster is based on the threshold distance d T 7 that has to be carefully chosen. If d T is chosen too large, there will be just one cluster including all points; if it is too small, no clusters will form with more than one point.
• in order to have sufficiently well formed clusters, the time interval t −t 0 must be chosen large enough so that the transient behaviour is gone; an illustration of the convergence is given in Fig. 7.
Depending on the type and amplitude of the forcing γ, we can have potentially a whole range of possible number of attracting trajectories N, ranging from one [Tziperman et al., 2006], to a few (two in the 41 kyr periodic forcing example in Figs. 5(c) and 5(d), or three in the quasiperiodic insolation forcing example in Figs. 5(e) and 5(f)). When no forcing is considered 5 See also Fig. 18 for a detailed view. 6 Cluster analysis or clustering is the assignment of a set of observations into subsets (called clusters) so that observations in the same cluster are similar in some sense. This is a common technique for statistical data analysis used in many fields for countless applications. There exists many types of clustering, along with several methods, among which: hierarchical, partitional, spectral, kernel PCA (principal component analysis), k-means, c-means and QT clustering algorithms. 7 This threshold distance d T appears in any computation related to clustering analysis, or analogically Recurrence Plots (RP) analysis in complex networks, for determining neighbours [Marwan et al., 2009. (c) Numerical estimate of the number of attracting trajectories N for the 41 kyr periodic forcing case. In practice, small positive values correspond to synchronization onto a few attracting trajectories, while high values indicate no synchronization. Now the region inside the synchronization tongues is colored in function of N. For the tongue corresponding to a frequency-locking n : 1, we have n attracting trajectories. The bifurcation boundaries of the Arnol'd tongues obtained with the more accurate numerical continuation method AUTO are superimposed, for validation purposes (black curve), and match perfectly.
(d) Numerical estimate of the number of attracting trajectories N in the case of the quasiperiodic insolation forcing given by Eq. (1): the structure is now much more complex, consisting of intermingled series of Arnol'd tongues. The region with one attracting trajectory, corresponding to unique or monostable generalized synchronization, is the largest. However, there are also parameter sets with N = 2, 3 or even more attracting trajectories. (Figs. 5(a) and 5(b)), or there is forcing but no synchronization occurs, we find no clusters at all. This means that there are as many points in the (x, y)-section at time t as initial conditions at time t 0 . Clearly, it is difficult to numerically distinguish between no synchronization and a large number of attracting trajectories (N 1). Therefore, we restrict ourselves to just six different regions in Fig. 6, where we use white to indicate when there are none or more than five attracting trajectories. Now, we apply the numerical clustering analysis in the case of the periodic forcing ( Fig. 6(c)) and of the quasiperiodic forcing ( Fig. 6(d)). We set t = 0 and consider a grid of 49 initial conditions covering x ∈ [−2.2; 2.2] and y ∈ [−2.2; 2.2] at the initial time t 0 = −40T F for the periodic forcing (T F is the period of the forcing), and t 0 = −1600 kyr for the astronomical forcing. Two points in the (x, y)-section are estimated to belong to a different cluster if their Euclidean distance is greater than 0.1.

kyr periodic forcing
An illustration of the three possible synchronized solutions (N = 3) existing for the 3:1 frequency-locking on a periodic forcing is given in Fig. 8, where the response can be locked on one of the periods of the forcing. More generally, N corresponds to the num- ber of forcing cycles associated with the synchronization regime (N = 1 for 1:1; N = 2 for 2:1; N = 3 for 3:2, 3:1, etc.) 8 . The resulting pattern of different N is, as expected, in agreement with the bifurcation diagram ( Fig. 6(c)). For example, N = 3 in the 3:1 tongue. This method allows one to visualize the 4:3 (N = 4) and even the 5:4 (N = 5) tongues to the left of the 3:2 tongue. It is also seen that N is generally larger where different synchronization regimes co-exist; this is the case between the 2:1 and 3:1 regimes around γ = 4.
Astronomical quasiperiodic forcing Fig. 6(d) shows that synchronization occurs for most parameter configurations.
The region with one attracting trajectory (N = 1), corresponding to unique or monostable generalized synchronization [Rulkov et al., 1995], is the largest. However, there are also parameter sets with N = 2, 3 or even more attracting trajectories. They indicate multistable generalized synchronization where different possible stable relationships (2) between the forcing and the oscillator response coexist.
A closer view in the lower values of γ is given in Fig. 9, which allows an insightful physical interpretation.
Three tongues with N = 1, 2, 3 are rooted at T ULC /T ε 1 = 1, 2 and 3, respectively, suggesting a synchronization on the main obliquity component of the astronomical forcing of the same nature as synchronization on a periodic forcing. A series of other synchronization tongues with N > 1 appear; they correspond to 2:1 (N = 2), 3:1 (N = 3), 4:1 (N = 4) and even 5:1 (N = 5) synchro-nization on the three leading components of precession, denoted p 1 , p 2 and p 3 . Consequently, the richness of the astronomical forcing effectively widens the parameter range for which synchronization occurs, compared to a periodic forcing. The phenomenon may be understood intuitively: just as you are more likely to tune on some radio station if you are surrounded by a dozen of free FM emitters, the system is more likely to synchronize on the rich astronomical forcing than on a periodic forcing. Synchronization with N = 1 or 2, found for larger γ, can be interpreted as a form of combined synchronization on both obliquity and precession.
It is crucial to appreciate that synchronized solutions are not periodic and that, unlike the periodic forcing case, different synchronized solutions for a given set of parameters are not time-shifted versions of each other. The idea that different synchronized solutions coexist is of practical relevance for paleoclimate theory. Namely, the set of parameters used to obtain the fit to the paleoclimatic records shown in Fig. 3 give two distinct solutions at t = 0 when started from a grid of initial conditions at t 0 = −700 kyr. Sensitivity studies show that the choice of t − t 0 is sometimes important for estimating correctly N. However, tests with t − t 0 as large as 200 Myr of astronomical time suggest that several attracting trajectories may co-exist at the asymptotic limit of t 0 → −∞.
A similar numerical clustering analysis plot for β = 0.6 is shown in Fig. 10 in order to give an idea of the effect of this parameter (a more detailed analysis is performed in Sect. 5). The main conclusion about the multistability remains, but the particular values of N change, as the intermingled tongue series are different.

Evolving shape of the basins of attraction
Each AT i (i = 1 . . . N) has its own basin of attraction 9 [Barnes and Grimshaw, 1997], that is defined as the set of all initial conditions in the (x, y,t) phase space that converge to that AT i as time tends to infinity. For our nonautonomous system Eq. (4), we can study basins of attraction in the (x, y)-section for different but fixed values of initial time t 0 , and observe how they vary with t 0 . A given initial condition at time t 0 lies in the basin of attraction of AT i if it approaches AT i as time tends to infinity. Technical details about the computation of the basins of attraction by use of the specific classification algorithm developed (see Fig. 18) are given in the Appendix D. Basins of attraction are of major importance because they provide the information about global or nonlinear stability of synchronization. If we care about predictability, basin boundaries indicate when a change in the attracting climatic history is likely.
The evolving shape of the basins of attraction is shown in Figs. 11 and 12, for the case of 41 kyr periodic forcing and the quasiperiodic astronomical forcing, respectively. The evolution is shown as a comic strip, where each subfigure has an x−axis ranging from -1.5 to 1.5, and a y−axis ranging from -2.5 to 2.5 (the axis labels have been removed for a better readability).
In the case of a periodic forcing (two basins), the pattern repeats itself periodically (compare the t 0 = 0 kyr to t 0 = 40 kyr, and to t 0 = 80 kyr subfigures) in Fig. 11. However, in the case of the quasiperiodic forcing (three basins), the pattern is much more intricate and seems not to repeat itself for the time horizon considered here.
The ratio between the area of a basin of attraction and the considered area of the phase space can be interpreted as a probability to converge to the corresponding attracting trajectory when starting from a randomly chosen initial condition. In the case of the periodic forcing, the two AT s are roughly equally likely for all t 0 as could be guessed from Fig. 5(c). However, this is not the case for the quasiperiodic forcing where the probability to reach the same attracting trajectory may vary significantly in time. For example, the yellow basin is rather small at t 0 = 0 kyr but becomes much larger at a later time t 0 = 90 kyr.
In the multistable regime, if an AT i happens to lie sufficiently close to its basin boundary, then small perturbations could make the climate jump to another (coexisting) AT j =i , reducing predictability. This phenomenon is illustrated in Sec. 6.

Influence of the symmetrybreaking parameter β
As the parameter β controls the asymmetry of the glaciation/deglaciation saw-tooth structure (a higher value of β leads to an enhanced asymmetry), it is useful to investigate its effect. We have already indicated in Fig. 10 that multistability depends on β in the case of the quasiperiodic insolation forcing. A more systematic approach encompassing the whole range of β is shown in Figs. 13(a) and 13(b), for the case of 41 kyr periodic forcing and the quasiperiodic astronomical forcing case, respectively. First consider the 41 kyr periodic forcing ( Fig. 13(a)). To understand this Figure, recall that the unforced oscillator (i.e., γ = 0) has a stable fixed point for |β | > 1 and a stable limit cycle for |β | < 1.
The system responds almost linearly to the forcing when |β | is sufficiently large. This explains regions of unique synchronization (N = 1) where only one climate response is possible. The system becomes excitable when |β | is just slightly greater than one. If the forcing is large enough it will excite oscillations. In this case, N is equal to the number of initial conditions if synchronization is lost, or to a smaller number if synchronization occurs.
Consider now the interval −1 < β < 1. For this, keep in mind (i) that the period of the unforced oscillation varies by almost a factor of two within the range 0 < |β | < 1, and (ii) that synchronization requires some relation between the period of the unforced oscillations and the forcing period. Consequently, synchronization on the periodic forcing occurs only for fairly narrow ranges of β that are symmetric around zero. Wieczorek ure reminds us of Arnol'd tongues. The main synchronization regimes detected here correspond to 4:1, 3:1 and 5:2 frequency-locking. Outside these synchronization regimes, the system fails to converge to a sufficiently small set of attracting trajectories, meaning that the forcing is not as an efficient pacemaker. Finally, compare this situation with that obtained with the astronomical forcing ( Fig. 13(b)). Synchronization now occurs in a larger area of the parameter space. Whereas the structure of the periodic forcing is preserved as long as the forcing amplitude is low enough, there is a much more richer and complex pattern of different N for larger γ. This pattern emerges from the interaction with different harmonics and their beatnotes.

Robustness of synchronization
Robustness or reliability of synchronization can be studied in terms of two properties of an attracting climatic trajectory. Local stability analysis based on the short-term LLE (λ H max ) provides information about the short-term local convergence towards the AT . For example, a temporary loss of local stability indicated by λ H max > 0 will cause a temporary loss of synchrony and divergence from the AT even though the trajectory is stable on average (λ max < 0). Global stability analysis based on the geometry of basins of attraction for different AT s provides information about the system response to external perturbations such as random fluctuations. For example, an external perturbation may push a climatic trajectory outside of the basin of attraction of the AT . Robustness and uniqueness of synchronization are closely linked in the sense that the global stability is a factor only when there are coexisting attracting trajectories. Robustness is compromised most when a temporary loss of local stability coalesces with a weakening of the global stability. We will now briefly discuss these two effects that could restrict the prediction horizon for the evolution of climatic trajectories.

Temporary desynchronization via loss of local stability
Some additional experiments made in our paleoclimatic framework reveal another strange behaviour in the system Eq. (4), that can be deduced from a careful inspection of Fig. 14. In the presence of small additive noise, we notice that nearby trajectories could diverge x Figure 14: The temporary divergence of nearby climatic trajectories reveals short-term local instabilities (e.g. around t ≈ 157 kyr). This illustration was obtained by considering a set of 50 random initial conditions at 1], and by adding some noise of small amplitude in the system at every time step in order to trigger the instabilities. The model used is Eq. (4) with α = 11.11, β = 0.25, γ = 0.033, and τ = 33.33.
for some time, like those around t ≈ 157 kyr. Such temporary divergence is similar to desynchronization bursts [Rulkov et al., 1995] and strongly suggests to investigate the evolving sign of the short-term LLE λ H max along the attracting climatic trajectory. We computed λ H=50 kyr max along one of the two attracting trajectories of system Eq. (4), subject to insolation forcing given by Eq. (1). The result is shown in Fig. 15, where the attracting climatic trajectory has been coloured according to the values of λ H=50 kyr max . Although the system is synchronized on a long term (λ max = −0.2 kyr −1 , see Fig. 17), we see here that there exist episodes with positive values of the short-term 10 LLE λ H max , revealing temporary desynchronizations [Wieczorek, 2009]. This explains the divergence of nearby trajectories found in Fig. 14. These results remained unchanged with respect to the most important parameters of the model. For example, our main conclusions about the stability re-(a) (b) Figure 13: Numerical estimate of the number of attracting trajectories N plotted as a function of the asymmetry parameter β and the amplitude of the forcing γ for the ice age model Eq. (4) with α = 11.11, β = 0.25, and τ = 35.09, assuming (left) a 41 kyr periodic forcing and (right) the quasiperiodic astronomical forcing given by Eq.
(1). Synchronization occurs in a larger area of the parameter space in response to the astronomical forcing than in response to the periodic forcing, and the pattern of different N for larger γ is much more richer and complex; this structure emerges from the interaction with different harmonics of the astronomical spectrum and their beatnotes.  main qualitatively valid, even for different values of α (like α = 100), or with a different type of potential (Φ 5 (y) = (y+1.7)(y+1.58)(y+0.8)(y)(y−0.5)), even if the shape and size of the limit cycle and the boundaries of the basins of attraction are of course different. The effect of the insolation function F(t) has also been checked: we compared the attracting trajectories for the insolation given by Eq. (1), compared to those for the insolation given by [Laskar et al., 2004]. As these insolation functions are very similar, the results are also very similar, and no difference was noticed.
At first glance, it may appear that these episodes of temporary divergence are not relevant to the robustness of synchronization because climatic trajectories converge back the the attracting trajectory on a long term. However, other effects may be present that could strongly amplify such temporary divergence. They are identified below. Fig. 12 showing (x, y)-sections with coexisting attracting trajectories in the case of the quasiperiodic insolation, and their basins of attraction for different values of t 0 . Suppose now that the system is subject to additive fluctuations (for example, these may represent volcanic eruptions). Under certain conditions, such external perturbations may cause a displacement of the trajectory to a different basin of attraction, causing a jump 11 to another attracting trajectory.

Consider again
As a further illustration of this idea we show in Fig. 16 two attracting trajectories (in the time series format) that coexist for the same system parameters as those used for the fit of Fig. 3, but with additive fluctuations added to the fast variable (see legend for details). A jump from one trajectory to another at around -475 kyr (arrow) may clearly been identified. This shows that a climatic trajectory is robust against fluctuations if it stays away from the basin boundary but its robustness can weaken significantly due to the weakening of the global stability near the basin boundary.
We conclude that externally triggered jumps between coexisting attracting climatic trajectories seem to be most likely when the temporary desynchronization due to the loss of local stability coalesces with the weakening of the global stability due to the proximity to the basin boundary. 11 In the periodic forcing case the phenomenon of jumping from one attracting trajectory to another in response to a perturbation is called a phase slip [Pikovsky et al., 2001, p.238].  (4) as in Fig. 3 is plotted (black) along with one sample trajectory of the same system (red), but with additive fluctuations added to the fast variable: dy = −τ −1 (α(Φ (y) − x))dt + b dW , with b = 0.5 √ ω ε 1 and W a Wiener-process. The arrow shows the time of the jumping from one to the other climatic attracting trajectory, which reduces the predictability of the timing of the glaciations.

Conclusions
We have identified, illustrated, and provided a systematic study of generalized and multistable synchronization between the climatic glacial/interglacial oscillations and the astronomical forcing. For doing so, a series of appropriate concepts and tools have been developed. A van der Pol-type relaxation oscillator, designed to reproduce the slow-fast dynamics of the paleoclimatic records, has been used for illustration purposes, but the methodology proposed here may of course be applied to other paleoclimatic models.
To study the uniqueness of synchronization, we proposed a convenient concept of the number of attracting trajectories in the phase space of the nonautonomous forced system, each of which corresponds to a synchronized solution. We computed the number of synchronized solutions using a numerical clustering technique, and uncovered that in addition to a unique or monostable synchronization, there are parameter settings where one finds a nonunique or multistable synchronization. At low forcing amplitude we found regions of mode locking where the system synchronizes on the individual components of the astronomical forcing in a way that is similar to frequency-locking on periodic forcing (Arnol'd tongues), giving rise to coexisting synchronized solutions. As the forcing amplitude is increased, the combined effects of precession and obliquity restrict the number of possible synchronized solutions. The emerging stability diagram consists of a large region of monostable synchronization mixed with smaller regions of multistable synchronization. A comparison with periodic forcing shows that the system finds it easier to synchronize to quasiperiodic insolation forcing. It is therefore conceivable that the climate system wandered throughout preferential synchronization regimes on obliquity, precession, or combinations of both, as environmental parameters varied throughout the history of the Pleistocene.
The robustness of generalized synchronization was investigated in terms of the key indicators of stability of synchronized solutions: the long-and short-term largest Lyapunov exponent (local stability), and the geometry of the basins of attraction (global stability). We found that even though the synchronized solutions are locally stable on a long term, there exists episodes where the short-term largest Lyapunov exponent becomes positive, leading to temporary desynchronizations. As a result, climatic trajectories could diverge from the synchronized solution for some period of time (50 kyr typ-Article submitted to Climate Dynamics B. De Saedeleer, M. Crucifix and S. Wieczorek ically). Moreover, we computed the evolving shape of the basins of attraction for the coexisting synchronized solutions, and uncovered that these solutions sometimes approach the basin boundary where they become very susceptible to external perturbations. As a result, a small perturbation could make the climate jump from one synchronized solution to another, reducing predictability. Such jumps seem to be most likely when the temporary loss of the local stability coalesces with the proximity to the basin boundary. In this context, we briefly discussed the effect of stochastic perturbations on the timing of the deglaciations. We also illustrated the difference between the evolving shape of the basins of attraction for periodic and quasiperiodic insolation forcing. In the case of the insolation forcing, we obtained an intricate pattern of basins of attraction that does not appear to repeat itself in time.

Acknowledgements
We are grateful to Guillaume Lenoir for his thorough review of the several versions of the paper. The original idea of using clustering analysis for automatically identifying the number of stable locking states came to the main author, after presentations and dis-   , 1926] is very well known, widely used and deeply studied; here the necessary description which is relevant to the scope of this article is given. For more details the reader is referred to the literature, see e.g. [Strogatz, 1994].
This van der Pol-type oscillator model can be regarded as a special case of the FitzHugh-Nagumo (FHN) model [Kosmidis and Pakdaman, 2003], also known as Bonhoeffer-van der Pol (BVP) model [Barnes and Grimshaw, 1997].
Historically, this model was derived by van der Pol, when he discovered the existence of oscillations in electrical circuits. He found that the oscillation period is determined by τ * = RC (time constant of relaxation) in RC circuits, or by τ * = L/R in RL circuits ; hence he named this oscillation as relaxation oscillations. The interesting characteristics of the relaxation oscillation are the slow asymptotic behavior and the sudden discontinuous jump to another value. This oscillator has therefore been widely used in many fields like in physical and biological sciences, neurology, seismology, lasers, optoelectronics, etc.
This two-dimensional system exhibits the same basic dynamical features (limit cycle and slow-fast dynamics) necessary to fit the paleoclimatic data, except the asymmetry, for which the system will be corrected for, by the introduction of an additional parameter β (see Eq. 4a in Sect. 1).
This van der Pol oscillator is non conservative with a nonlinear damping, governed by the following secondorder differential equation [Balanov et al., 2009] -in the case without forcing : where µ > 0 is a positive constant proportional to the damping. The van der Pol oscillator is a Liénard system [Strogatz, 1994], because f (x) = −µ(1 − x 2 ) is an even function and because g(x) = x is an odd function (see Eq. B.1), in the canonical form : Moreover, since this Liénard system satisfies additionally the Liénard theorem, it has a unique and stable limit cycle 12 in the phase space surrounding the origin. This is a feature which is absolutely required in order to model the glacial-interglacial oscillations.
By way of the Liénard transformation y = x − x 3 /3 − x/µ, the second-order Eq. B.1 can be transformed into an equivalent two-dimensional system of ODE's : which can also be expressed under the equivalent following form (using the convention ε = −1/µ, so that ε → 0 means more and more damping): the associated potential Φ(x) having the shape of a 2well potential. By varying ε, we can adjust the respective time scales of the x and y variables; for ε < 1, the slow-fast variable is x and the slow variable is y. Note that in the main body of this paper, x and y are inverted, so that x will represent the slowly varying ice volume (see Fig. 4), as it is the classical convention for the dynamical theory of paleoclimates. When µ > 0, the system exhibits a stable limit cycle, where energy is conserved. Near the origin x =ẋ = 0, the system is unstable and energy is gained, and far from the origin the system is damped and energy is lost (while when µ = 0, there is no damping, the solution is a pure harmonic signal, and there is conservation of energy everywhere).
There is only one fixed point which is the origin (x,ẋ) = (0, 0).
When x is small, f (x) ≈ −µ (negative damping). Thus, the fixed point (x,ẋ) = (0, 0) is unstable (an unstable focus when 0 < µ < 2, and an unstable node otherwise); see [Hilborn, 2000]. On the other hand, when x is large, f (x) ≈ +µx 2 (positive damping). Therefore, the dynamics of the system is expected to be restricted in some area around the fixed point.
When 0 < µ 1 (small damping), the system can be rewritten in order to avoid division by µ. When µ 1 (large damping), the oscillations become less and less symmetric.
When driven, the van der Pol oscillator can lead to synchronization [Balanov et al., 2009], but also to deterministic chaos [Ruihong et al., 2008], depending of the level of the driving force. Since the paleoclimatic system is driven by the insolation, one could also have synchronization and perhaps chaos -or at least some path on the routes leading to chaos; i.e. Lyapunov exponents becoming positive. This oscillator, or slightly different versions of it (a similar one is the Poincaré oscillator [Glass and Sun, 1994]), has been mathematically largely studied under many aspects: bifurcation structure [Mettin et al., 1993], fixed points and Arnol'd tongues, chaotic dynamics [Chen andChen, 2008, Parlitz andLauterborn, 1987], additive noise [Degli Esposti Boschi et al., 2002], basins of attraction [Barnes and Grimshaw, 1997], etc. In this article, we focus on the synchronization, multistability and predictability properties of a slightly modified version of the classical van der Pol oscillator.
But there is also a big difference of frameworks: usually, a simple periodic forcing is considered, while in the case of the astronomical forcing, the forcing has a much more complex form (quasiperiodic), like the one given in Eq. (1).

Time spent on the slow manifold
In the case of large damping (µ 1, or ε → 0), the oscillations become less and less symmetric, and significant differences becomes to appear between τ slow and τ f ast , as defined in Section 1.
One can more precisely characterize the slow-fast dynamics of the system by considering Eqs. B.5-B.6. We see that we have in any case |ẏ| ∼ ε, which means that y is always slow.
When the system is away from the curve y = Φ (x), we have (|ẋ| ∼ 1/ε) (|ẏ| ∼ ε) : the vector field is mostly horizontal. And the travel speed is v = |ẋ| 2 + |ẏ| 2 ∼ 1/ε 1, so that the system moves quickly in the horizontal direction (y is the slow variable, and x is the slow-fast variable). If y > Φ (x), theṅ x > 0, and the trajectory moves clockwise. See Fig. 4 for an illustration and [Strogatz, 1994] for a discussion about the separation between the two time scales.
When the system enters the region 13 where |y − Φ (x)| ∼ ε 2 , we then have (|ẋ| ∼ ε) ∼ (|ẏ| ∼ ε) : |ẋ| an |ẏ| are about the same order of magnitude, and the travel speed is v ∼ ε → 0 : the system moves slowly along the curve y = Φ (x), and eventually exits from this region, and so on : the system has a stable limit cycle.
Hence we have the following relationships: The period of the oscillations (or of the unperturbed limit cycle), given by is mainly determined by the time during which the system stays around the curve y = Φ (x), which can be roughly estimated to be T ∝ 1/ε, as v ∼ ε → 0. So, the larger the damping (lower ε), the longer the period of the oscillations; that's why an additional time scaling (factor τ) must sometimes be done in order to keep the 100 kyr period for the limit cycle (see Fig 4). Let us mention also that analytical expressions have been derived for the amplitude and period of the limit cycle [D'Acunto, 2006] in both cases (µ small or big), and also for the slow manifold equation [Ginoux and Rossetto, 2006]. There exists also a very fast transition upon variation of a parameter β , from a small amplitude limit cycle to a large amplitude relaxation cycle, explained by the so-called canard phenomenon, cycles, and explosion [Benoît et al., 1981, Guckenheimer et al., 2000, Guckenheimer and Haiduc, 2005.

C Technical details about the calculation of the Lyapunov exponents
We first refer the reader to the seminal papers [Shimada andNagashima, 1979, Benettin et al., 1980]. The methods for computing the Lyapunov (Characteristic) Exponents (LCE) vary depending on the fact that one wishes to achieve the full spectrum of the LCE [Wolf et al., 1985], or only the largest one [Rosenstein et al., 1993]. Analytical derivations of the LCE's of the van der Pol oscillator do also exist [Grasman et al., 2005]. In this research, we computed λ max using the standard method involving a Gram-Schmidt Reorthonormalizaton (GSR) of the 'tangent vectors' [Shimada and Nagashima, 1979, Benettin et al., 1980, Wolf et al., 1985; which is described in the review paper [Ramasubramanian and Sriram, 2000].
We remind here the fundamental principles. Consider an n-dimensional continuous-time dynamical system : where Z and f are n-dimensional vector fields. To determine the LCE's corresponding to some initial condition Z(0), we have to find the long term evolution of the axes of an infinitesimal sphere of states around Z(0). That is to say that we assume [Ott, 2002] δ Z(t) = δ Z(0) e Λt , with Λ = diag(λ 1 , . . . , λ n ), where λ i are the eigenvalues of the system. For this, we consider the linearization of Eq. C.1, given by: where J is the n × n Jacobian matrix defined by J i j = ∂ f i /∂ Z j . Then, starting from a unit vector δ Z(0), the original system given by Eq. C.1 is integrated for Z together with the δ Z tangent system given by Eq. C.2. The evolution of δ Z is such that it tends to align with the most unstable direction (the most rapidly growing one). The choice of the initial vector of the tangent manifold may influence the convergence, but in practice a spinup phase can be performed in order to find the good direction.
The largest Lyapunov exponent λ max is then defined as : It is of course impossible in practice to go to infinity 14 ; so the computation is always truncated to some finite final time, usually of the order of 10 4−5 times the period of the forcing. The convergence can be more rapidly achieved if some transient behaviour is skipped (like illustrated in Fig. 17), i.e. we compute the LCE's 14 The proof of the existence of such a limit has been made by [Oseledec, 1968]. only when we are quite sure to be on the attracting trajectory.
Note that there exists a whole spectrum of Lyapunov exponents λ i (i = 1, 2, ..., n), that can be computed with a unit vector basis, and with a renormalization procedure. Although we are mostly interested in λ max in this article -a positive λ max is associated to a desynchronization -, our subroutine allows to compute all the spectrum of λ i 's of any system. For the sake of flexibility, we used a symbolic software, so that the model could be very easily changed, and all functions (like the Jacobian matrix) are automatically derived once the initial system is given.
One of the standard and popular methods to compute the Lyapunov spectrum of a dynamical system involves a Gram-Schmidt Reorthonormalizaton (GSR) of the 'tangent vectors' [Shimada and Nagashima, 1979, Benettin et al., 1980, Wolf et al., 1985; differential versions of this have also been formulated. Our subroutine includes the GSR procedure, which is required in order to avoid computational overflows, and degeneracy into a single vector. The frequency of reorthonormalization is not critical; as a rule of thumb, GSR is usually performed on the order of once per characteristic period. Here, the normalization time step has been chosen optimally, that is to say the largest possible which preserves the accuracy. Since the GSR never affects the direction of the first vector in a system, this vector tends to seek out the direction in tangent space which is most rapidly growing.
Our subroutine has been validated by comparing our LCE's to those of [Ramasubramanian and Sriram, 2000] for several systems (driven van der Pol, Lorenz '63, etc.); the order of the accuracy achieved is 1 for the Lorenz system 15 .
Coming back to our system of Eqs 4a-4b, for which the Jacobian matrix is: we end up to the following results (Fig. 17): the system Eq. (4) subject to the insolation (Eq. (1)) is synchronized on a long term, since λ max = −0.2 kyr −1 < 0. [kyr ] -1 Figure 17: Convergence of the two LCE's of one of the two attracting trajectories of the system Eq. (4), subject to the insolation (Eq. (1)); some transient is skipped. For τ = 3.33, we have λ max = −0.2 kyr −1 < 0 (for τ = 35.09, the value must be scaled accordingly, which would give -0.019 kyr −1 ), so the system is synchronized on a long term.

Some properties of the LCE's (λ i )
The LCE's are very useful in order to characterize the dynamical behaviour of a system; for example, here are a few interesting properties, valid in the case of the dissipative system Eq. (4) in the autonomous case : • a dynamical system of dimension n has n LCE's and n eigenvectors [Lichtenberg and Lieberman, 1983], • ∑ n i=1 λ i = Tr( J) = div f , both being related to the growth of a volume of dimension n of the phase space [Shimada and Nagashima, 1979], • div( f ) < 0 for a dissipative system (hence the system has at least one negative exponent), • λ i = 0 along a limit cycle (tangent direction of the attractor) [Kantz and Schreiber, 2004], • at least one LCE vanishes if the trajectory of an attractor does not contain a fixed point [Haken, 1983], • the nature of the attractor can be described by analyzing the sign of the LCE's, which gives a qualitative picture of the dynamics [Wolf et al., 1985, Müller, 1995, also function of the dimension (1D to 4D), e.g. an attractor for a dissipative system with one or more λ i > 0 is said to be strange or chaotic; if more than one λ i > 0, there is hyperchaoticity [Rossler, 1979].
• the existence of λ i > 0 is mathematically related to the theory of ergodicity 16 of dynamical systems [Eckmann and Ruelle, 1985], • local bifurcations can be detected by detecting changes of signs of λ i .
The LCE spectrum is closely related to the fractional dimension of the associated strange attractor by the Kaplan-Yorke formula [Kaplan and Yorke, 1979]. There are a number of similar measures of the "strangeness" of strange attractors, like the fractal dimension [Theiler, 1990], information dimension, box-counting dimension, and the correlation dimension [Tsiganis et al., 1999] and exponent [Grassberger and Procaccia, 1983], which allows to distinguish between deterministic chaos and random noise, and is computationally easier.
It is also possible to estimate the LCE of a system by analyzing its time series, but with limited data, or a system subject to non negligible stochastic perturbation, the classical methods may provide incorrect or ambiguous results [Ruelle, 1990], hence require specific methods, like [McCaffrey et al., 1992, Liu et al., 2005.

D Technical details about the computation of the basins of attraction
The practical computation of a basin of attraction is done as follows.
Let us for example come back to the Fig. 5(e), with the three attracting trajectories AT i (x, y,t) due to the insolation forcing. Now, we wonder which initial condition leads ultimately to which of the three AT i . This is the concept of the basin of attraction of a given AT i , classically defined by the set of states that leads to a given AT i . Let us more precisely define the basin of attraction of a given AT i as 16 An attractor is said ergodic (or transitive) if the points fill an entire physically realizable domain; intransitive if there are distinct attractors in several closed subdomains [Saltzman, 2002], and a system is almost intransitive if it resides for long periods of time in one or another domain that is not fully closed off, so that occasional exits from one domain to another can occur.
Article submitted to Climate Dynamics B. De Saedeleer, M. Crucifix and S.
Wieczorek the locus of all points in the (x, y) plane which lead to motion which ultimately converges on that AT i . We initialize many initial conditions on a fine rectangular grid covering the phase space. Each initial condition is then integrated forward to see which AT i its trajectory approached. If the trajectory approached a particular one of the three AT i 's, a dot coloured by the color identifying the AT i is plotted on the grid. For doing this, we need to define a target time at which we will do the classification, and a criteria for the classification.
The classification algorithm is illustrated in Fig. 18, where a cut has been made at a target time t = 550 kyr. The three AT i 's are displayed, together with the location of the trajectories (black circles). To decide if a given trajectory ends up onto a given AT i , we choose a maximum distance from the AT i (dotted circle); taken here to be 1/4 of the minimum distance between two AT i 's. Trajectories falling into that dotted circle are classified as ending into AT i .
We consider a given trajectory starting from t 0 = 0 kyr, integrate it up to t = 550 kyr, and then we examine its position with respect to the AT i at the same time T=550 kyr. If the distance to a given AT k is 'sufficiently small' (see the circles around the AT i ), then this initial condition is coloured in the k th color, associated with the k th attracting trajectories AT i .
Repeating this process for each initial condition on a fine grid of the whole phase space gives the shape of the basins of attraction. As we have three attracting trajectories, we will have three basins of attraction. Such basins of attraction obtained are shown in Fig. 19.
As we are in the case of a non-autonomous system, we then have to repeat this procedure for several starting times t 0 (the position of the AT i 's are constantly evolving, hence the shape of the basins are also varying with time). This has been done to produce the evolving shape of the basins of attraction in Figs. 11 and 12. Note that if t is too close from t 0 , transient behaviours predominate, and not enough time has elapsed in order for the trajectory to be attracted by a given AT i , hence the basins of attraction cannot be defined in that case.
The glacial/interglacial cycles do exist since about 3 million years, but only 8 limit cycles of 100 kyr period have been performed, since the Mid-Pleistocene Transition (MPT). This should however be sufficient to converge onto the attracting trajectories, because the climatic trajectories are rapidly attracted on the limit cycle. So, if the ice age model Eq. (4) is a realistic one, it would be reasonable to state that the transient of climate dynamics has gone, and that we are currently probably Figure 18: Illustration of the classification algorithm, which allows to compute the basins of attractions. The position of the 70 trajectories (black circles) with respect to the three attracting trajectories (symbols) at time t = 550 kyr, starting from t 0 = 20 kyr, are shown. If the trajectory falls into a dotted circle, classification is possible. The 70 trajectories form three highly concentrated clusters, corresponding to the three attracting trajectories.  Figure 19: When integrating the system system Eq. (4), with many initial conditions (grid 201 × 121 = 24 321) starting from time t 0 = 0 kyr, with the insolation forcing given by Eq. (1), we obtain the three basins of attraction above; each basin is coloured with the colour of the corresponding AT i . The function Φ (y) = y 3 /3−y = x, corresponding to the slow manifold, is also shown (dashed curve).

E Instantaneous stability
As we considered the long-term (3 Myr) and short-term (H = 50 kyr) local stabilities, we could wonder why not to consider also, at the other extreme, the 'very shortterm' local stability (H → 0). Even if it has probably less relevance in the context of predictability of glaciations, it may however still be computed and understood as an instantaneous stability. As the horizon of time H tends to zero, it means that we no longer consider any motion in the phase space (no average in time), hence the instantaneous stability becomes also a local property in the phase space. The same formula (Eq. 6) as for the short-term stability is used, but now with an Horizon of H = 1 kyr, representing the very short-term. The result is plotted in Fig. 20, where the trajectory has been coloured with the λ H=1 max values. There are black areas, which means local instantaneous instability. This black zone is always in the same region of the phase space: when y lies within the interval I * = [−1, 1], whatever the initial condition.
We now demonstrate mathematically why the black zone in Fig. 20 left is located within the interval I * = [−1, 1].
Let us derive the theoretical expression of the exactly instantaneous LCE λ inst i . If we call Λ i the eigenvalues of the jacobian matrix J, then the λ inst i are defined as λ inst i = lim t→0 1 t ln Λ i [Drazin and King, 1992], so that we can obtain the largest λ inst i by computing λ inst max = max(ℜ(Λ i )).
As the Jacobian matrix J is a function of y only (and not x), see Eq. C.3, so do λ inst max . The function λ inst max (y) is plotted (Fig. 20, right, black dotted curve), where we clearly see that it is positive on the interval [−1, 1], which is precisely the same interval as I * .
Moreover, the numerical values also match (e.g. the maximum of the function λ inst max has a value of about 0.28, which correspond to the maximum of the color scale).
An interpretation of λ inst max is that a particular trajectory makes a travel through the phase space 'painted' by λ inst max , integrating so the instantaneous stability to achieve the short-term and long-term stabilities.