Chaotic and non-chaotic response to quasiperiodic forcing: limits to predictability of ice ages paced by Milankovitch forcing

It is well known that periodic forcing of a nonlinear system, even of a two-dimensional autonomous system, can produce chaotic responses with sensitive dependence on initial conditions if the forcing induces sufficient stretching and folding of the phase space. Quasiperiodic forcing can similarly produce chaotic responses, where the transition to chaos on changing a parameter can bring the system into regions of strange non-chaotic behaviour. Although it is generally acknowledged that the timings of Pleistocene ice ages are at least partly due to Milankovitch forcing (which may be approximated as quasiperiodic, with energy concentrated near a small number of frequencies), the precise details of what can be inferred about the timings of glaciations and deglaciations from the forcing is still unclear. In this paper, we perform a quantitative comparison of the response of several low-order nonlinear conceptual models for these ice ages to various types of quasiperiodic forcing. By computing largest Lyapunov exponents and mean periods, we demonstrate that many models can have a chaotic response to quasiperiodic forcing for a range of forcing amplitudes, even though some of the simplest conceptual models do not. These results suggest that pacing of ice ages to forcing may have only limited determinism.


Introduction
Periodically forced nonlinear oscillators are textbook examples of nonlinear systems whose attractors can exhibit chaotic behaviour (Ott, 1993). Chaotic attractors can similarly be found for more complex (e.g. quasiperiodic) forcing (Kazuyuki, 1994) along with other structures such as strange nonchaotic attractors (Feudel et al., 2006;Lai et al., 1996).
Since the work of Milankovitch (Milankovitch, 1941), it has been suggested that forcing of the earth's climate by slow, approximately quasiperiodic, changes in the earth's orbit have been responsible for repeated periods of slow growth (glaciation) and rapid retreat (deglaciation) of the Northern Hemisphere ice sheets over the last few million years. These orbital variations, in turn, create oscillations in the northern-hemisphere summer insolation on time scales of obliquity (about 41 kyr) and precession (several periodicities between 19 and 23 kyr, Berger (1978)). The climatic precession cycles are modulated by eccentricity (100 kyr and 400 kyr) (Berger and Loutre, 1991;Laskar et al., 2004). The observed glacial cycles of the Pleistocene transition from rather weak ice ages cycling on a 40 kyr time scale during the early Pleistocene (∼ 2.5 Myr -1 Myr ago) to large-amplitude, asymmetric ice ages cycling on a 100 kyr time scale during the Late Pleistocene (∼ 1 Myr ago -present) (Huybers, 2007;Lisiecki and Raymo, 2005). In an influential paper, Hays, Imbrie and Shackleton (Hays et al., 1976) suggest that the summer insolation forcing is (i) the fundamental cause of the observed oscillations in glacial cycles of the late Pleistocene and (ii) important for pacing these ice ages (in particular, the timing of the deglaciations).
However, since its introduction there have been a number of problems with this hypothesis, particularly for the late-Pleistocene glacial cycles (Imbrie and Imbrie, 1980;Ruddiman et al., 1986;Imbrie et al., 1993;Lisiecki, 2010;Paillard, 2015). (i) During the entire Pleistocene, the astronomical forcing is dominated by variations on the precession (about 23 kyr) and obliquity (about 41 kyr) time scales, but the late Pleistocene response is dominated by variations on the 100 kyr time scale. (ii) The relatively weak forcing by eccentricity on the 100 kyr time scale is decreasing in power throughout the mid and late Pleistocene while the power of the 100 kyr response is growing. (iii) The power in variations of the eccentricity forcing on the 400 kyr time scale is comparable to the power at the 100 kyr time scale, but there is negligible power in the glacial cycle response around 400 kyr. None of results are consistent with the hypothesis that the late-Pleistocene glacial cycles are a direct response to astronomical forcing.
On the other hand, there is ample evidence that the timing of the late-Pleistocene glacial cycles are influenced by the phase of the obliquity and/or eccentricity forcing (Huybers, 2007;Lisiecki, 2010). Feedback processes internal to the climate system, for example affecting the ice mass balance, could then amplify the response to astronomical forcing on specific time scales such that the 100 kyr time scale appears in the late Pleistocene climate response (Abe-Ouchi et al., 2013). Another possibility is that the emergence of the long late-Pleistocene cycles is related to an internal oscillation in the climate system, with approximate period of 100 kyr or longer (Ditlevsen, 2009;Crucifix, 2012;Daruka and Ditlevsen, 2016), excited by the Milankovitch forcing starting around the Mid-Pleistocene transition. Such an internal oscillation has its own periodicity, but if it is paced by the quasi-periodic astronomical forcing the period observed in the final response can be different from the unforced case. In this paper we focus on the modification of an internal oscillation by quasi-periodic forcing.
The word pacing suggests a causal but indirect link between a complex multi-period forcing and a response that is somewhat weaker than synchronization. This relates to the key questions around the predictability of future glaciations and deglaciations. This paper aims to further explore the connection between quasiperiodic forcing of nonlinear oscillators and the predictability of the Pleistocene ice ages. We have been inspired by work of Crucifix and collaborators (Crucifix, 2012;De Saedeleer et al., 2013;Mitsui et al., 2015) to explore several of the oscillator models in the literature under periodic and quasiperiodic forcing. We note that there are however other types of nonlinear resonance, such as response to a resonance system that has no natural oscillations in the absence of a varying input (Marchionne et al., 2016).
There is a variety of low-order models for the late Pleistocene ice ages, ranging from simple (nonlinear) oscillators to more physically motivated models. For example the ubiquitous van der Pol oscillator has been used to describe general features of the ice-age dynamics and has been suggested as a possible minimal model (De Saedeleer et al., 2013). Some of these models such as the van der Pol oscillator are strongly constrained: chaotic attractors are associated with canard behaviour where parts of the attractor are contained in strongly repelling parts of phase space (Itoh and Murakami, 1994). In other words, a forced van der Pol oscillator only gives chaotic attractors for very thin regions in parameter space (Bold et al., 2003). This is however not a universal property of forced nonlinear oscillators -they may have chaotic attractors for much larger regions in parameter space.
For the remainder of Section 1, we briefly review nonlinear oscillator models of the Pleistocene ice ages and possible roles of synchronization and chaos in understanding pacing. In Section 2, we consider the van der Pol oscillator and note that a slightly more complex model (in particular, the van der Pol-Duffing oscillator) can have much wider parameter regimes of stable chaos than seen in (De Saedeleer et al., 2013). Section 3 studies two low order physics-based models of Saltzmann-Maasch 1991, Palliard-Parrenin 2004: we quite readily find attracting chaotic responses of these models to quasiperiodic forcing. Finally, in Section 4, we discuss implications of the results and some open problems. Although our studies do not imply the response of the Pleistocene ice-ages to astronomical forcing is necessarily chaotic, one cannot rule out this possibility. We consider what this would imply about the determinism of deglaciations over long timescales. Some computational details are in Appendix A.

Nonlinear oscillator models of the Pleistocene ice ages
Many different mechanisms have been suggested, both in terms of mathematical structure and physical basis, for the remarkable late Pleistocene ice ages. For example, Hagelberg et al. (1994); Imbrie et al. (1993) suggest a (linear) resonance with the eccentricity forcing could explain the 100 kyr periodicity of the late Pleistocene ice ages and at the same time the absence of a 400 kyr periodicity in the records. However, other analyses of proxy records and models seem to suggest they are governed by a nonlinear process that is in some sense phase-locked to astronomical forcing (Tziperman et al., 2006). This mechanism relies on the precession and obliquity time scales of Milankovitch-forcing and requires a nonlinear oscillator whose period varies with amplitude. A nonlinear resonance condition determines the glacial period, which according to the phase-locking mechanism should be 'quantized' (meaning that the ratio of the oscillation period and one forcing period should be a ratio of two integers), although this depends on the definition of a 'phase' which is not obvious for the physical system or the data.
As discussed in Crucifix (2012);De Saedeleer et al. (2013); Mitsui et al. (2015), conceptual models play an important role in terms of providing terminology and mechanisms that might be understood as pacing. It is presumably a sign of usefulness of these ideas that many different nonlinear oscillator models have been found that can produce time series that match the ice volume record. However, it is not even clear whether the glacial cycles should be self-sustaining (existing without forcing) or have stable states that are excited by the forcing. Tziperman et al. (2006) suggest the ice-age problem can be usefully be split into two sub-problems, where one involves the understanding of the pacing and timing of glacial cycles while the other should seek for a physical mechanism giving rise to the glacial cycles. The nonlinear phase locking mechanism provides a potential explanation for the first subproblem, but makes the second more difficult because models with different mechanisms can equally well match the proxy records -this paper concentrates on the first sub-problem.
The models we consider in this paper have the form where some finite-dimensional internal dynamics on x ∈ R d governed by f (x, λ) represents processes that model the ice-age dynamics for a fixed insolation anomaly λ, and Λ(t) represents the time-dependent insolation anomaly induced by astronomical forcing. We measure t in units of kyr. The models we consider all have periodic dynamics for the unforced case (Λ = 0 is constant) and show a nonlinear response to forcing. We consider several types of forcing, the first being periodic Λ 0 (t) = k 1 sin(ω 1 t), where ω 1 represents pure obliquity forcing, ω 1 = 2π 41 . Secondly we consider two-frequency quasiperiodic (QP2) forcing that represent components of the obliquity and precession forcing for (ω 1 , ω 2 ) = ( 2π 41 , 2π 23 ), where k 1,2 gives the relative amplitudes of these frequencies (strictly speaking this is quasiperiodic only if ω 1 /ω 2 is irrational). We also consider more realistic multi-frequency quasiperiodicity (QPn) forcing where the values of s k , c k , ω k are listed in (De Saedeleer et al., 2013, Appendix 1) and we define to provide a comparable normalization of k 1 and k 2 to (3). This is a good model for the last few million years of insolation, and the peak powers of Λ 2,o and Λ 2,p are around 41kyr and 23kyr respectively. Note that for the latitude 65 o N, the precessional components average out over the year, whilst the obliquity components do not. Hence we suggest it is reasonable to vary these components independently as a model for the effects of seasonal integration, even though many authors (e.g. De Saedeleer et al. (2013)) consider only the case k 1 = kµ 1 and k 2 = kµ 2 . In summary, the QPn forcing (5) is a more realistic version of the simplified QP2 forcing (4). We mainly consider four models in our numerical exploration: the 2D models of van der Pol and van der Pol-Duffing, and the 3D models of Saltzman and Maasch (1991) and Paillard and Parrenin (2004). All have been proposed as possible models for the Pleistocene ice ages but we highlight the van der Pol model stands out as being dynamically simpler than the others, in that chaos is confined to very narrow regions.

Synchronization and phase locking
In this paper, we try to distinguish between different types of nonlinear models in terms of the response to quasi-periodic forcing. There is a trade-off between simplicity of models for ice age pacing and realism in terms of physical processes and resolution included. While there have been attempts to integrate earth system models of intermediate complexity including ocean, ice sheet and carbon cycle dynamics over several glacial cycles (Abe-Ouchi et al., 2013;Ganopolski and Brovkin, 2017), these models are by far too complex to analyse and understand the process of pacing of the glacial cycles by the astronomical forcing. These questions about pacing can be more adequately answered by studying less complex, simplified models that have been used in the past for the glacial-interglacial cycles. Some of those are physically motivated by including only a limited number of physical processes and restrict the spatial resolution typically to a few "boxes" (e.g. Gildor and Tziperman (2001)), or only one global signal (e.g. Saltzman and Maasch (1991); Paillard (1998)), while others include more mathematical concepts of self-sustained oscillators Crucifix (2012).
The strongest form of pacing one can imagine consists of generalized synchronization to the signal, i.e. there is a functional relationship between the instantaneous value of the forcing to that of the signal -this is clearly not a realistic model for signals such as the proxy of Lisiecki and Raymo (2005) so we should consider weaker versions. This includes, for example, synchronization to a filtered component of the forcing, or a more loose phase locking such that some phase extracted from the forcing remains close to a phase extracted from the response.
The nonlinear phase-locking mechanism of Tziperman et al. (2006) argues that any nonlinear model that becomes appropriately phase-locked would be able to describe the observed ice-volume record. Consequences of nonlinear phase locking suggested by Tziperman et al. (2006) are: (i) Even in the presence of abundant noise in the climate system the phase locking can be effective. (ii) The glacial period should be multiples of the precession and obliquity forcing, although this ratio can change over time. Eccentricity forcing does play a more indirect role by modulating the precession and obliquity forcing amplitudes. (iii) The quasiperiodic nature of the Milankovitch forcing allows for varying glacial periods; however, it was suggested that the timing of glacial terminations could be uniquely determined (even though the phase of the Milankovitch-forcing is not the same during all glacial terminations but spreads around 1/4 of a cycle in the model used by Tziperman et al. (2006)). This last consequence implies the model is insensitive to small changes in initial conditions, because eventually all time series become phase-locked to the same forcing. Similarly, Tzedakis et al.
(2017) use a statistical model to predict the series of past deglaciations suggesting there exists only a small set of possible realisations for the deglaciations during the late Pleistocene. Note that, as highlighted by Marchionne et al. (2016), complex patterns of nonlinear resonance can appear even in the absence of limit cycle oscillations for the unforced system: we do not consider such cases here.
There is no unique way to extract a continuously varying phase from a complex signal unless it is close to periodic: spectral methods have the disadvantage that they are global in time, assume stationarity of the signal and may have significant end effects for finite time intervals. There are many ways to do this, for example by measuring the rotation around a point for a projection of phase space onto a plane. We use a well-studied method of extracting the phase by using the Hilbert transform of the oscillatory signal x(t) (Pikovsky et al., 2001;Le Van Quyen et al., 2001). Let us denote the mean of x(t) as x. This phase is where is such that x + iy is analytic: the integral is understood as a Cauchy principal value. Using this we can compute R(T ), the mean number of rotations of the phase of x(t) over the interval [0, T ], as 2πT and the mean period is 1/R(T ). This works well for signals that have one crossing of the mean x per period but is not so reliable for signals with several crossings per period. For a quasiperiodic signal that has many competing frequencies, different linearly filtered signals will show completely different phases and so conclusions from this (as with any) method of phase extraction need to be treated with care.

Pacing, attractors and the role of chaos
At the weaker end of pacing, one could argue there is some sort of synchronization to forcing if the response loses all information about its initial condition -for example if a neighbourhood of initial conditions converge to the same response trajectories. These response trajectories are pullback attractors for the forced non-autonomous system (De Saedeleer et al., 2013), and there may be several possible local pullback attractors for a given forcing (Kloeden and Rasmussen, 2011).
Given a solution x(t) ∈ R n of (1), we compute the attracting behaviour for some generic initial condition. For stationary forcing the Lyapunov (characteristic) exponents are given by For a typical initial (x, v), the Lyapunov exponent λ can take one of up to n possible values λ 1 ≥ λ 2 ≥ · · · ≥ λ n that represent the possible exponential rates of attraction or repulsion of nearby trajectories: see for example Ruelle and Eckmann (1985). The largest Lyapunov exponent (LLE) λ 1 enables a characterization of the attracting dynamics into chaotic (λ 1 > 0) or non-chaotic (λ 1 ≤ 0). More details of how the LLE is computed is described in Appendix A. Pullback attractors can be characterised using LLEs of the forced system: they are isolated trajectories in the case of negative LLE while a positive LLE implies the pullback attractor is a chaotic set. Even in the case of QPn forcing and a negative LLE one may find sensitive dependence on phases -this corresponds to strange nonchaotic attractor (SNA) (Lai et al., 1996) for the system (1) extended to describe the forcing Λ = L(θ 1 , . . . , θ n ) where (θ 1 , . . . , θ n ) ∈ T n and d dt For SNAs, there is a unique point pullback attractor for any given realisation of the forcing, but the location of that point may vary wildly on changing one of the phases θ i while fixing the others. SNAs have been extensively studied and found in a variety of QP2 forced systems (Feudel et al., 2006), including phase oscillator models for Pleistocene ice age oscillations (Mitsui et al., 2015).
There is a complex relationship between chaotic behaviour and phase locking -as we will see, QP forcing for many nonlinear oscillators with two or more variables readily leads to chaos in the sense of a positive LLE. The presence of chaos in the response does not however mean there is no phase locking -the chaos may be small scale and cause a jitter of the precise phase with no loss of phase locking. It does imply a form of non-determinism in that the future of the trajectory will depend sensitively on initial condition at all points in the future. However for the models we investigate here, we find chaos that is associated with loss of phase locking.

Simple oscillator responses to forcing
We first consider responses of some simple oscillators to periodic and quasiperiodic forcing before moving on to more realistic oscillators in Section 3.

The van der Pol oscillator and generalizations
The modified van der Pol oscillator was suggested as minimal model for the ice age cycles in Crucifix (2012);De Saedeleer et al. (2013). This model has mostly negative LLEs for both periodic and quasiperiodic (two periods) forcing indicating that it unlikely shows a chaotic response to forcing anywhere in the forcing-parameter space.
The van der Pol model of Crucifix (2012); De Saedeleer et al. (2013) involves two dependent variables: x is a slow variable representing deviation of ice volume from some reference and y is a fast variable representing some feedback mechanism with hysteresis.
where the forcing term Λ(t) represents astronomical forcing and can be either periodic (e.g. only obliquity component) or quasiperiodic. The parameters we consider in Table 1 are from De Saedeleer et al. (2013), and for convenience and comparability with their results we use two parameters τ and κ to normalise the period around τ = 1. System (8) can be written more explicitly as the van der Pol oscillator with a linear restoring force (the nonlinearity is only in the damping): A generalization of this can be obtained by replacing the linear restoring force with a nonlinear g(y) as in the following.  (8) and van der Pol Duffing (9) models: all variables are dimensionless, except for time which is in kyr. Note that a nonlinear response has been chosen that give relaxation-like oscillations and κ is chosen to ensure the unforced oscillations have period approximately 100 kyr for time scaling τ = 1.
Note that the system (9) is equivalent to a generalised van der Pol-Duffing oscillator: with nonlinear restoring force g(y). Figure 1 illustrates typical unforced solutions of the two systems in time series and phase portrait. The van der Pol-Duffing oscillator does show regions of potentially chaotic response to periodic and quasiperiodic forcing depending on the amplitude of the forcing.

Responses to periodic and quasiperiodic forcing
The relative lack of chaotic behaviour under periodic, QP2 and QPn forcing for the van der Pol oscillator (8) was noted in De Saedeleer et al. (2013). However this is somewhat specialadding nonlinearity in the restoring force can lead to more prevalence of chaos. Figure 2(a,b) contrasts two examples of timeseries of these oscillators subject to periodic forcing while (c,d) gives a similar contrast for QP2 forcing. For the van der Pol oscillator (Figure 2(a,c)), four randomly chosen initial conditions synchronize to the same response trajectory; while in the case of the van der Pol-Duffing system (Figure 2(b,d)), the synchronization is only intermittent. Scans of the parameter space spanned by the forcing amplitude k and the time scaling τ in Figure 3(a) show the LLE for the van der Pol model (8) with forcing (2), while Figure 3(b) show the LLE for a similar scale with the van der Pol-Duffing generalization (9). Note that τ varies the unforced period with τ = 1 corresponding to the expected 100 kyr in both cases: compare (a,b) to (De Saedeleer et al., 2013, Fig 6). As expected, for small forcing amplitudes k 1 , weak forcing of a limit cycle oscillator gives frequency locking on an attracting invariant torus. These Arnold tongue regions of p:q frequency locking appear (Pikovsky et al., 2001) where the ratio of forcing to natural frequency is close to p : q. The tongues where p and q are small grow preferentially with k 1 , meaning that at for larger k 1 in (a,b) we see most of parameter space is taken up by low order frequency locking with ratios 1:1, 2:1 and 3:1.  Table 1. In (b,d) the solid red lines show the y-nullcline and the dashed green lines shows the x-nullcline: note that there are two additional x-nullclines for the latter system. Observe that both show relaxation oscillations with approximately the same period. Moreover, the critical manifold (y-nullcline) is identical in both cases.  (9) with typical parameters as in Table 1 forced by (3) with four random initial conditions started at time t = 0. For (a,c) there is periodic forcing (3) with k 1 = 1, k 2 = 0 while for (b,d) there is QP2 forcing (3) with k 1 = 1, k 2 = 1. Note the strong synchronization to a single trajectory for the van der Pol system (a,c) and the presence of intermittent fluctuations in the transitions for the the van der Pol-Duffing system (b,d) typical for a chaotic response.
For larger k 1 these regions start to overlap, associated with break-up of the invariant torus. This torus break up means the attractor can explore a higher dimensional region in phase space and potentially become chaotic. However, chaos is notable by its absence in Figure 3(a). In fact, chaotic attractors are forced to live in 'cracks' between frequency locking. This sparsity of chaos for the forced van der Pol oscillator has been studied in the literature (Bold et al., 2003;Itoh and Murakami, 1994), and is associated with the fact that chaotic attractors involves canard trajectories which confine the chaos to thin strips in parameter space for this model.
By contrast, the forced van der Pol-Duffing oscillator has chaos associated with recurrent approaches to a saddle equilibrium which makes it more robust in parameter space, as illustrated in Figure 3 (9) to the same forcing. One can numerically verify that the LLE for (c) is negative while that for (d) is positive and chaotic. This is illustrated in Figure 4 which shows (a,c) LLEs and (b,d) mean periods as a function of k 1 and k 2 for QP2 forcing (3) of (a,b) the van der Pol (8) and (c,d) van der Pol-Duffing (9) systems. Observe the presence of significant regions of chaotic attracting responses for the van der Pol-Duffing system (Figure 4(c)).

Physics-based conceptual model responses to forcing
We consider two different physics-based models of Pleistocene ice-age oscillations. The first of these is the Saltzmann and Maasch model (SM91) (Saltzman and Maasch, 1991) which is a three-variable nonlinear model where all three dynamic variables have similar time scales. The second is the Palliard and Parrenin model (PP04) (Paillard and Parrenin, 2004) which is a relaxation oscillator (with three dynamic variables) where the switch between glacial and interglacial state is induced by a threshold in Antarctic bottom water formation efficiency. We will see that both SM91 and PP04 models show significant areas in parameter space where a chaotic response to periodic or quasiperiodic forcing (either QP2 or QPn) is possible.

The Saltzman and Maasch 1991 model
The models of Saltzman et al. (Maasch and Saltzman, 1990;Saltzman, 2002;Saltzman and Maasch, 1988, 1990, 1991Saltzman and Verbitsky, 1993) established the idea that the (late Pleistocene) glacial cycles appear as a limit cycle in the unforced system then synchronised in some way to the orbital forcing. The dynamics (and consequently the specific form of the limit cycle and the bifurcation leading to the limit cycle) varies across the different models. In most of them, it is assumed that the background climate slowly varies throughout the Pliocene-Pleistocene ('tectonically driven' decline in atmospheric CO 2 -concentration) and the model is formulated as an anomaly model to a (slowly evolving) background climate. (c) (d) Figure 3: Scan of (a,c) LLEs and (b,d) mean periods for typical initial conditions of (a,b) van der Pol equation (8) and (c,d) van der Pol Duffing equation (9) with periodic forcing (2). We vary the time scaling τ (where τ = 1 corresponds to an unforced period 100 kyr) on the horizontal axis and the forcing amplitude k = k 1 on the vertical axis. Observe in both cases there is an Arnold tongue structure for small k = k 1 and regions of 1:1, 2:1 and 3:1 locking (left to right) that persist to large k = k 1 for (a,b) but include regions of more complex chaotic dynamics for (c,d) (compare with (De Saedeleer et al., 2013, Figure 6(a,c))).  (8) and (c,d) van der Pol Duffing equation (9) with QP2 forcing (3). The forcing amplitudes k 1 and k 2 of the QP2 forcing are varied on the two axes, while the unforced period is 100 kyr in both cases (τ = 1). Crucifix (2012) analysed the bifurcation structure of two of these models (Saltzman and Maasch, 1990Maasch, , 1991 for the full (non-anomaly) equation with respect to one parameter, the slowly drifting 'tectonic' CO 2 decline F µ (t). In these models, much of the interesting dynamics depends on the specific form of the CO 2 equation, which at the same time is the most problematic to interpret physically. The model SM91 we use here is from Saltzman and Maasch (1991) and includes one specific representation of the carbon cycle dynamics. The model couples the dynamics of the ice volume I, CO 2 concentration µ and deep-ocean temperature θ as in Crucifix (2012): The insolation forcing is Λ(t) and other parameters are as in Table 2: as in Crucifix (2012), we consider the original variables instead of the anomaly variables in that the latter constrained the possible bifurcations and attractors of the system. In Crucifix (2012), this is analysed with respect to variations in the parameter F µ , while the other forcings Λ(t) and F θ are set to zero. We use here fixed values of the 'tectonic forcing' F µ = F θ = 0 (see Crucifix (2012)) and periodic, QP2 or QPn forcing Λ(t). We fix k R = 0.4 to give comparability with other models in terms of response to forcing. Figure 5(a) illustrates time series of a typical response of this system (10) to periodic forcing (2) while (b) shows a corresponding response to QP2 forcing (3). The evolution of four initial conditions are shown on the same axes -note that in both of these cases there is synchronization to the same trajectory, typical of an attractor with a negative LLE. For k 1 = 0 (precession only) there is also locking present for larger k 2 (not shown). Note the presence of a region of chaotic locking to approximately 100 kyr with positive LLEs: Figure 5(c) shows an example of where the four time series repeatedly and intermittently synchronize for periods of time while (d) shows a projected Poincaré section sampled when the phase of obliquity is zero modulo 2π for a longer single trajectory from (c). Figure 6(a,c) illustrates LLEs and (b,d) mean periods for a scan with periodic forcing (2) for this model. Panels (a,b) show scans for varying time scaling τ (which scales the unforced oscillation period) and forcing amplitude κ = k 1 while panels (c,d) show scans through amplitudes of the components for QP2 forcing (3). There are clearly regions of chaotic behaviour in parameter space, mostly concentrated at larger values of the precessional component amplitude k 2 = 1 and k 1 small. For k 2 = 0 (obliquity only) observe that there is a large region of 1:1 locking to the obliquity for larger k 1 , and a region of 3:1 locking for smaller k 1 .  Table 2. Panel (a) shows responses to periodic forcing (2) with k 1 = 1 and (b) shows response to QP2 forcing (3) with k 1 = k 2 = 2. In both cases there is synchronization to a single trajectory. Panel (c) shows QP2 forcing with k 1 = 0.8, k 2 = 0.8 and a response consistent with chaos. The Poincaré section (d) corresponds to one of the trajectories in (c) sampled at zero phase of the obliquity forcing.

The Paillard-Parrenin 2004 model
The model of Paillard and Parrenin (2004) (PP04) gives an explicit qualitative physical explanation for switches between glacial and interglacial regimes. The key assumption here is that deglaciations are induced by changes in atmospheric CO 2 . In this case, the (slow) transition from interglacial to glacial can be explained by Milankovitch theory, where the summer insolation in northern high latitudes play a major role. The faster glacial to interglacial transition (climate escapes from a deep glacial) is caused by the release of CO 2 from the deep ocean into the atmosphere. In this model, this is achieved by an oceanic switch, characterised by a salty bottom water formation efficiency parameter F . During cold periods, the Southern Ocean bottom water formation is strong, the deep ocean is strongly stratified with cold and salty Antarctic bottom water (AABW) and can store a lot of carbon. However, a few thousand years after a glacial maximum, the Antarctic ice sheet reaches its maximum extent making salty AABW formation through brine rejection on the continental shelf more difficult, the Southern Ocean stratification weakens while at the same time atmospheric CO 2 rises rapidly. The model involves three variables, the global ice volume V forced by CO 2 and Northern Hemisphere summer insolation, the Antarctic ice sheet extent A driven by sea level changes (via V ) and the atmospheric CO 2 concentration C connected to the deep water state (glacial Value Units Physical interpretation α 1 1.673915 × 10 16 kg yr −1 Constant ice growth rate. α 2 9.52381 × 10 15 kg yr −1 Effect of CO 2 and θ on ice dynamics. α 3 10 −4 yr −1 Inverse linear response time ice. β 1 0.5118377 ppm yr −1 CO 2 coefficients. β 2 6.258680 × 10 −3 yr −1 CO 2 coefficients. β 3 2.639456 × 10 −5 (ppm yr) −1 CO 2 coefficients. β 4 3.628118 × 10 −8 (ppm 2 yr) −1 CO 2 coefficients. β 5 5.833333 × 10 −3 ppm ( o C yr) −1 CO 2 coefficients. γ 1 1.85125 × 10 −3 o C yr −1 Rate of growth of θ. γ 2 1.125 × 10 −23 o C(kg yr) −1 Effect of ice on θ. γ 3 2.5 × 10 −4 yr −1 Inverse linear response time θ.
Insolation normalization τ 1 Time scaling The oceanic (salty bottom water) switch parameter F is given by: The forcing I 65N = Λ(t) and the parameter values are given in Table 3. We approximate the Heaviside function in (11) using H(x) = (1 + tanh(Kx))/2 with K = 100. It should be noted that the only nonlinearity in this system is the Heaviside function (or its tanh approximation) which provides a state-dependent switch between two linear modes: (i) an ice-accumulation mode where the system relaxes towards a low CO 2 , high ice volume state and (ii) an ice-ablation mode where the systems relaxes towards a high CO 2 , low ice volume state.  Table 3: Chosen values of parameters and their interpretation for the PP04 model (11) of Paillard and Parrenin (2004): these correspond to the parameter values used in Crucifix (2012). As V , A and C are dimensionless quantities, most parameters are dimensionless.
The unforced system exhibits internal oscillations at a 132 kyr period for the parameters used by Paillard and Parrenin (2004). Figure 7(a) illustrates time series of a typical response of this system to periodic forcing -in this case there are several attractors corresponding to mode-locking to different phases. Figure 7(b) show responses of this system to QP2 forcing (4) with k 1 = k 2 = 1 while (c) shows the responses for k 1 = 0.5, k 2 = 0.2 and an example with apparent chaotic response (see also Figure 11 in the Appendix). Computing LLEs for a scan through k 1 the amplitude of the obliquity and k 2 precessional components, Figure 8(a,b) shows the response of PP04 (11) to periodic forcing (2)  Finally, in Figure 9, we illustrate scans of LLE for (a) SM91 and (b) PP04 with Milankovitch forcing (5) with varying obliquity and precessional amplitudes k 1 and k 2 . Observe a similar distribution of regions of locking as in Figures 6 and 8, but, as expected, with a considerably more complex bifurcation structure associated with the QPn forcing. Note that the forcing in De Saedeleer et al. (2013) corresponds to a point on a line (k 1 , k 2 ) = (kµ 1 , kµ 2 ) ≈ k(26, 76).

Discussion
In summary, a lack of chaotic response to forcing is a property of some nonlinear oscillators that have specific structures -in particular, for the van der Pol oscillator in the relaxation oscillation regime, the chaos is associated with canard solutions which only appear as part of the attractor for narrow regions of phase space. Phase oscillators can also avoid chaotic responses; for example, consider the phase oscillator model of Mitsui et al. (2015): where φ ∈ R/2πZ is a periodic phase variable. The one-dimensional nature of the system means that LLEs on attractors are non-positive: chaos is not possible, even though strange non-chaotic attractors (SNAs) can appear in this system Mitsui et al. (2015). We show here that more general oscillator models, including some physically motivated ones, can show larger regions in parameter space where there are chaotic attractors. We find this for the van der Pol-Duffing oscillator as well as for the models of Saltzman and Maasch (1991) and Paillard and Parrenin (2004), and this is the case not only for periodic forcing, but also for two-frequency forcing and a realistic approximation of Milankovitch forcing. We expect to find similar behaviour even in systems that have no stable limit cycle in the unforced case, but if there is non-linear resonance with the forcing (Marchionne et al., 2016).
Certain features need to be present to show prevalent chaotic response, firstly the forcing needs to be of sufficient amplitude and, at least for two dimensional oscillators, there needs to be some shear in phase space near the limit cycle. This has been found to give rise to chaotic response in a variety of systems in that it gives a mechanism for the required stretchand-fold in phase space (Blackbeard et al., 2011). To illustrate, consider the forced shear  (11) with parameters as in Table 3: observe in this case there is no 1:1 locking but rather a subharmonic locking. Panel (a) shows responses to periodic forcing (2) with k 1 = 1 and (b) shows response to QP2 forcing (3) with k 1 = k 2 = 2. Panel (c) shows QP2 forcing with k 1 = 0.5, k 2 = 0.2. Poincaré section (d) corresponding to one of the trajectories in (c) sampled at zero phase of the obliquity forcing. oscillator of Lin and Young (2008): where θ ∈ R/(2πZ) is a phase variable and y ∈ R. If we choose parameters P = 100, σ = 3, λ = 0.2, γ = 0.1 and scan through k 1 and k 2 we find large regions of strongly chaotic responses for QP2 forcing (4); see Figure 10.

Implications for modelling the Pleistocene ice ages
This study does not imply either the presence or absence of a chaotic response to forcing of the physical 'ice-age' oscillator by astronomical forcing. However, it does suggest that a chaotic response cannot be ruled out simply because it is rare for the forced van der Pol oscillator. Indeed, both of the physics-based models show more chaos more readily, but still only in regions of parameter space that may or may not be relevant to the late Pleistocene ice ages. We do note however that in order to get an approximately 100 kyr response to a periodic forcing at approximate 41 kyr period, we need to be in a region of moderate forcing amplitude -too great and we find 1:1 locking, too little and we find no locking. This happens to also be a region where we can find regions of chaotic response for the PP04 and SM91 models. If there is a chaotic response with LLE λ > 0, this would place a fundamental limit to the predictability of future ice ages, and in particular the onset of rapid deglaciations. Roughly speaking, suppose an estimation of the error in initial fraction of ice for an accurate ice age model is 0 < ε < 1. This will grow at approximate rate λ, meaning there is effectively a time horizon T h such that εe λT h = O(1), i.e. T h = − ln ε/λ. Since convergence of Lyapunov exponents is highly non-uniform, a more sophisticated estimate along the lines discussed in Heydt and Ashwin (2017) which includes state dependence of the convergence could imply the existence of longer or shorter time horizons depending on current state. For the examples considered in Figure 11, the implied LLE is approximately 0.004 corresponding to a timescale of e-folding of errors that is λ −1 ≈ 250kyr. As for weather forecasting, the time horizons represent a limit beyond which detailed forecasts are likely to be no better than average state estimation.
We remark on a number of possible interesting directions for further work. One is clearly to work with more realistic and accurate models of the ice-age oscillations. Another is to improve the forcing -the Milankovitch 65 o N summer solstice forcing is clearly only an idealization. However there are still considerable challenges to use data to constrain the models in a meaningful way, in particular to determine whether a model with positive LLE or negative LLE is more appropriate, and if negative whether any response is an SNA or not. Indeed, it would be interesting to explore the appearance of SNAs that may appear during the breakdown to chaos (Feudel et al., 2006) for these systems.   (15). Observe convergence to a LLE which, in this case, is positive.