Chenciner bubbles and torus break-up in a periodically forced delay differential equation
We study a generic model for the interaction of negative delayed feedback and periodic forcing that was first introduced by Ghil et al. in the context of the El Niño Southern Oscillation (ENSO) climate system. This model takes the form of a delay differential equation and has been shown in previous work to be capable of producing complicated dynamics, which is organised by resonances between the external forcing and dynamics induced by feedback. For certain parameter values, we observe in simulations the sudden disappearance of (two-frequency dynamics on) tori. This can be explained by the folding of invariant tori and their associated resonance tongues. It is known that two smooth tori cannot simply meet and merge; they must actually break up in complicated bifurcation scenarios that are organised within so-called resonance bubbles first studied by Chenciner.
We identify and analyse such a Chenciner bubble in order to understand the dynamics at folds of tori. We conduct a bifurcation analysis of the Chenciner bubble by means of continuation software and dedicated simulations, whereby some bifurcations involve tori and are detected in appropriate two-dimensional projections associated with Poincaré sections. We find close agreement between the observed bifurcation structure in the Chenciner bubble and a previously suggested theoretical picture. As far as we are aware, this is the first time the bifurcation structure associated with a Chenciner bubble has been analysed in a delay differential equation and, in fact, for a flow rather than an explicit map. Following our analysis, we briefly discuss the possible role of folding tori and Chenciner bubbles in the context of tipping.
Compared to fold or saddle-node bifurcations of equilibria or periodic solutions, the phenomenon of folding tori (also referred to as a saddle-node bifurcation of tori or a quasiperiodic fold bifurcation) is technically more challenging. In fact, the term “bifurcation” is used loosely here, because a fold bifurcation of two smooth invariant tori does not actually exist as such in a generic system. Namely, two smooth invariant tori cannot merge at a single bifurcation point. As the two tori approach the fold locus, they lose normal hyperbolicity (smoothness) and then break up in a bifurcation scenario that involves complicated dynamics [3, 8, 54]. The bifurcation scenarios that one encounters near the fold are organised in parameter space into so-called resonance bubbles, which were first studied by Chenciner [10, 11, 12] and are nowadays known as Chenciner bubbles.
The loss of smoothness of tori and their subsequent break-up have been observed in various systems; for example, an electric circuit , cardiac fibrillation data , the column pendulum  and a particle accelerator model . However, understanding the complicated dynamics involved in torus break-up is more challenging. In this paper, we analyse a delay differential equation (DDE), for which a fold bifurcation of tori was observed in . More specifically, we determine precisely how the torus break-up occurs by numerically calculating and identifying the bifurcation structure inside a Chenciner bubble.
The DDE studied in  was introduced by Ghil, Zaliapin and Thompson  as a conceptual model for the El Niño Southern Oscillation (ENSO) system. ENSO is the climatic driver behind the infamous El Niño events, which are characterised by exceptionally high sea-surface temperatures in the eastern equatorial Pacific Ocean that occur sporadically about every 3–7 years. The DDE model, which we refer to here as the GZT model, focusses on the interactions of a negative delayed feedback due to ocean-atmosphere coupling with periodic forcing representing the effect of the seasons. It is a special case of a model introduced by Tziperman, Stone, Cane and Jarosh  and takes the form:
The dependent variable represents the thermocline depth at the eastern boundary of the Pacific Ocean (more specifically, its deviation from the long-term annual mean), which depends on time measured in years. The parameters and of (1) are amplification factors of the negative delayed feedback and periodic forcing terms, respectively. The delay time of the negative feedback mechanism represents the time required for intermediate ocean-atmosphere coupling processes to close the respective feedback loop. In this model the delay is assumed to be constant. The specific form of the coupling function with coupling strength is justified in . Throughout, we set and in order to observe the fold bifurcations of tori seen in . Further details about the climate processes involved can be found in ; also see [29, Sec. 2].
Since folding tori and Chenciner bubbles are generic or typical phenomena, we expect that the case study presented here will be of interest for the bifurcation analysis of models arising in various applications. Because of the simple form that (1) takes, the GZT model could be viewed more generally as a model for the interaction between negative feedback and periodic forcing; very similar models arise, for example, in human motion control , network dynamics  and laser systems . Indeed, folding tori could be present whenever the dynamics involves two or more competing frequencies, or put simply, where tori exist; this includes models in the form of ordinary differential equations with a phase space of dimension at least three.
It is the dimensionality of DDEs that allows this seemingly simple model to generate complicated dynamics. Formally, the GZT model is of the form
where , consists of real parameters and
Note that (2) is nonautonomous, because depends on the time explicitly. For the special case equation (2) is a nonautonomous ordinary differential equation (ODE) and its phase space is , which we also refer to as the physical space. However, as a DDE with , equation (2) has the infinite-dimensional space as its phase space; here is the space of continuous functions over the delay interval with values in . This means that a whole function segment with values in over the time interval is required as part of the initial condition at time ; this is known as an initial history. Details on the general theory of DDEs can be found in, for example, [16, 24, 48].
In the GZT model (1) time enters in the form of a periodic forcing and, hence, the simplest invariant objects are periodic solutions. The stability of a periodic solution is given (as for ODEs) by its Floquet multipliers, which are the eigenvalues of the linearisation along . For a DDE with a finite number of fixed delays, there exists an infinite number of Floquet multipliers, but there are always at most a finite number of unstable ones [16, 24]. If all the Floquet multipliers (except for the trivial Floquet multiplier 1) lie within the unit circle of the complex plane, then is stable; on the other hand, if one or more Floquet multipliers lie outside the unit circle, then is of saddle type; together with the trivial Floquet multiplier 1, the corresponding eigendirections span the (finite-dimensional) linear unstable manifold of . It is convenient for the later discussion to distinguish saddle periodic solutions with one and two unstable Floquet multipliers, and we refer to them as 1-saddle and 2-saddle periodic solutions, respectively.
Simulations of (1) can readily be performed with a fixed-step numerical integration method, where the stepsize is given by uniform discretisation of the history segment. Simulations are only able to reveal stable solutions, and saddle solutions must and can be found with continuation software. We use the continuation package DDE-Biftool [17, 46], which runs in Matlab and is specifically designed for DDEs. Due to the periodic forcing of the GZT model, the simplest invariant solutions are periodic orbits. Once a periodic solution to (1) is found, we utilise DDE-Biftool to numerically continue (or track) the periodic orbit, regardless of its stability, while varying a system parameter. This is achieved by representing the periodic orbit as the solution of an appropriately discretised boundary value problem. DDE-Biftool determines the stability properties of the periodic orbit by calculating its Floquet multipliers, which are used to identify bifurcations. Codimension-one bifurcations — namely, saddle-node (or fold), period-doubling and torus (or Neimarck-Sacker) bifurcations — can in turn be continued in a two-parameter plane by placing appropriate constraints on relevant Floquet multipliers. For background information on continuation methods of DDEs see, for example,  and for implementation details specific to the GZT model, see [29, Sec. 1].
In our earlier work , we conducted a detailed bifurcation analysis of the GZT model (1) with DDE-Biftool. Generally, the dynamics is driven by two independent mechanisms that induce self-sustaining oscillations: the seasonal forcing and the negative delayed feedback. Their interplay can give rise to two-frequency dynamics on an invariant torus, which may be either unlocked or locked. If the frequencies have an irrational ratio, the trajectory on the torus will never close (i.e. it is unlocked) and the solution is called quasiperiodic. If the frequencies have a rational ratio, the trajectory is locked and a pair of periodic solutions (usually, a stable and a 1-saddle periodic solution) exist on the torus. In parameter space, the regions where the dynamics is locked to some fixed period are known as resonance tongues (also called Arnold tongues); these are bounded by saddle-node bifurcations of the respective locked periodic solution.
It was observed in  that, when certain parameter are changed, the invariant torus loses stability and disappears in a fashion that is analogous to a fold or saddle-node bifurcation of periodic solutions. More specifically, we showed that for certain ranges of values, as the parameter of (1) is varied, the observed stable torus approaches a co-existing torus of saddle-type in what appears to be a “fold” of tori.
Our approach to investigating folding tori in the GZT model is to study the locked solutions, organised into resonance tongues, that exist on the tori as they fold. Clearly, as the tori fold, so too do the associated resonance tongues. Inside each resonance tongue, where it folds, exists a Chenciner bubble : within such a bubble an invariant torus undergoes a transition from stable to saddle-type via a loss of normal hyperbolicity and subsequent break up. Therefore, in order to understand the phenomenon of folding tori it is essential to analyse and understand the dynamics in and around such Chenciner bubbles.
Investigating Chenciner bubbles in the context of climate systems could provide important knowledge towards understanding folding tori as a novel type of climate tipping mechanism. This type of tipping should be expected as a feature of systems with more than one frequency due to external forcing and/or internal feedbacks. In fact, oscillations involving multiple frequencies across a wide range are common in global climate systems . Many climate models, as well as observational data, reveal quasiperiodic dynamics (i.e. dynamics on tori); for example, [18, 25, 44, 50, 57]. Nonetheless, as far as we are aware, “multi-frequency climate tipping” has not been considered in the literature yet.
In this paper we consider Chenciner bubbles of folding resonance tongues in the -plane of the GZT model. We determine their structure by bifurcation analysis and demonstrate the resulting complicated dynamics. To this end, we focus on a resonance tongue. It is the largest one in the region of the -plane we consider and will therefore form the largest Chenciner bubble. We employ numerical techniques including the continuation software DDE-Biftool to detail the bifurcation structure of the Chenciner bubble. To our knowledge, the investigation in this paper is also the first time that Chenciner bubbles have been analysed in a DDE. Because it has an infinite-dimensional phase space, in order to interpret the dynamics of the GZT model we consider appropriate projections of the phase space. As part of the bifurcation structure we uncover in the Chenciner bubble, we detect curves of torus bifurcations, neutral saddles, homoclinic bifurcations of periodic solutions, heteroclinic bifurcations of periodic solutions and fold bifurcations of tori in the GZT DDE model. We find that this structure agrees well with the theoretical structure suggested in . One difference we find is that the torus bifurcation inside the Chenciner bubble we study is subcritical, instead of supercritical as in . We give detailed examples of how we calculate two-dimensional projections of the phase space of the DDE model, related to taking a Poincaré section, in order to detect the bifurcations that involve tori. In particular, we compute projections of unstable manifolds to find and illustrate locked tori. Finally, given the additional bifurcation structure present within Chenciner bubbles of folding tori, we briefly discuss the potential role of folding tori as a climate tipping mechanism.
The paper is organised as follows. We set the scene in section 2 by describing the dynamics within a region of the -plane where we expect Chenciner bubbles to exist. In section 3 we review some of the theory behind Chenciner bubbles, including the bifurcation structure as suggested in . We then present and discuss in section 4 a Chenciner bubble of the GZT model. In sections 4.3 and 4.4 we present detailed examples of transitions through certain bifurcations in the Chenciner bubble of the GZT DDE model. Finally, in section 5 we discuss how our results relate to climate tipping and point to some interesting directions of future work.
2 Background on the GZT model
We now present information on the GZT model that demonstrates the existence of the phenomenon of folding tori; see also . Figure 1 is a one-parameter bifurcation diagram in for fixed . The curves labelled ‘torus’ represent maxima of stable tori found by numerical integration, while increasing and decreasing the parameter slowly enough to ensure that the simulation stays on a particular branch of solutions until stability is lost. As a result, scanning up and down in reveals bistabilities. Stable, 1-saddle and 2-saddle periodic solutions are represented by their maxima, as computed with the continuation software DDE-Biftool. The point T denotes a torus bifurcations where the periodic solution loses stability while parameter is decreased. The points SN denote saddle-node bifurcations where stable and 1-saddle periodic solutions of the resonance tongue meet. Because the , , , , and resonance tongues are very small, they are represented by circles in figure 1; open circles represent resonance tongue with both stable and 1-saddle periodic solutions, while filled circles represent resonance tongues with 1-saddle and 2-saddle periodic solutions. In other words, the filled circles represent a branch of solutions on a torus of saddle-type that bridge the two branches of solutions on stable tori. Figure 1 convincingly demonstrates the existence of two loci of folding tori, labelled SNT. These may appear to be just like fold bifurcations of periodic solutions, but as mentioned above, theory predicts more complicated dynamics and bifurcations as the points SNT are approached.
Figure 2 shows three co-existing solutions for and , as indicated by the vertical line in figure 1. The solutions in rows (a) and (c) of figure 2 are stable and are found by numerical integration. The locked solution shown in row (b), on the other hand, is a saddle and calculated with DDE-Biftool. The solutions are displayed (after transients have died down), respectively, as a time series, a projection in the -plane, a stroboscopic trace in the -plane and a logarithmic power spectrum generated by Fourier transform from the time series over 1000 years. The stroboscopic trace is constructed by plotting after each forcing period the first point of the associated function segment, called the headpoint, in projection onto the -plane. Since the forcing period of the GZT model is 1, this means that in column 3 of figure 2 we plot whenever . See [9, 31] for a discussion of stroboscopic and Poincaré traces of DDEs.
The solution in row (a) of figure 2 belongs to the upper branch of stable tori in figure 1 and it is quasiperiodic (or of very high period). The time series in panel (a1) of figure 2 shows a slight modulation of peak heights and the trace of the headpoints in panel (a3) forms a closed curve. As such, the peaks in the power spectrum in panel (a4) are incommensurate with the frequency 1 of the seasonal forcing. The solution in row (b) is a 1-saddle periodic solution from the resonance tongue seen in figure 1. The phase space projection in panel (b2) of figure 2 is a closed curve and the stroboscopic trace in panel (b3) consists of 17 isolated headpoints. The solution in row (c) belongs to the lower branch of stable tori in figure 1 and is also quasiperiodic (or of very high period), which is evidenced by the closed loop in the stroboscopic trace in panel (c3) of figure 2. Because the peaks in panel (a1) are about twice the magnitude of those in panel (c1), the two stable solutions in rows (a) and (c) can be interpreted as two very different climatic states.
We now illustrate folding tori in the -plane. Figure 3 shows maximum maps of the GZT model, which plot the maximum value of attractors found by numerical integration for a range of fixed values and for gradually increasing in panel (a) and decreasing in panel (b), as indicated by the arrows. Maximum maps were calculated in [21, 56], albeit with a constant fixed initial history for each simulation; in contrast, we scan up and down in parameter in order to reveal multistabilities, as in figure 1. We remark that the vertical stripes inside some resonance tongues are not numerical artefacts; they are due to a symmetry property of (1), which manifests itself as a bistability in resonance tongues with even or ; see  for details. The maximum maps in figure 3 are overlaid with curves SN of saddle-node bifurcations of periodic solutions and a curve T of torus bifurcations, computed with DDE-Biftool. The vertical line indicates the parameter slice considered in figure 1.
The maximum maps in figure 3 are divided into two regions of small maxima (green to yellow) and large maxima (orange to red), respectively. The sharp interface in each maximum map, where the maximum of the observed solution changes rapidly from large to small or vice-versa as is increased or decreased, represents where solutions on tori lose stability and disappear. As demonstrated in figure 1, this happens where tori fold.
The folding tori are also evident in the bifurcation set of figure 3. Each resonance tongue is bound by two curves SN of saddle-node bifurcations of periodic solutions, emerging from a resonance point along the curve T of torus bifurcations. We observe that the resonance tongues, which are very narrow near the curve T, fold twice in the -direction. Moreover, the foldings of resonance tongues coincide with the sharp interfaces of the maximum maps. Specifically, the two envelopes of folding resonance tongues form the two boundaries where tori fold.
In the parts of the resonance tongue where the boundaries are drawn in blue in figure 3, one finds pairs of stable and 1-saddle periodic solutions; that is, these periodic solutions lie on a stable torus. In the parts of the resonance tongues where the boundaries are drawn in green, on the other hand, the locked periodic solutions have an additional unstable Floquet multiplier (as was checked with DDE-Biftool). Hence, one finds pairs of 1-saddle and 2-saddle periodic solutions; in other words, these locked periodic solutions lie on a torus of saddle-type, with one unstable direction. At these green segments of resonance tongues correspond to the filled circles in figure 1. It is where the resonance tongues fold in figure 3 and their boundaries change from blue to green that theory, which we briefly review in the next section, predicts the existence of Chenciner bubbles with complicated dynamics.
3 Theory of bifurcating tori and Chenciner bubbles
Bifurcations of tori, also referred to as quasiperiodic bifurcations, have been a subject of renewed interest in recent years; for example, see [5, 28, 32]. More specifically, quasiperiodic fold bifurcations are studied in, for example, [7, 54], where Chenciner bubbles are identified as regions that could potentially contain interesting dynamics.
To set the scene for our brief review of relevant theory, consider a family of autonomous vector fields in three dimensions with two real parameters, and , and a periodic orbit . For simplicity and to connect with the case considered in this paper, we assume that the vector field is a periodically forced two-dimensional oscillator with a constant forcing period of 1. Hence, the three-dimensional phase space is invariant under translation by 1 of the -axis. Without loss of generality we assume further that the origin is an equilibrium of the oscillator for all and , which generates the periodic orbit of period 1 under the forcing. To study the dynamics of the vector field one typically considers the dynamics of the Poincaré map , given as the stroboscopic map of the forcing frequency 1. The origin is a fixed point of and the (nontrivial) Floquet multipliers of are the eigenvalues of the linearisation of the map at the origin. The origin (and hence ) is stable if all the eigenvalues of are strictly inside the unit circle. Assume now that the origin is stable for and loses stability at , when a pair of complex conjugate eigenvalues of moves through the unit circle at .
While the bifurcation at is of codimension one, it is important to realise that the frequency ratio parameter also plays a role. In fact, it is best to consider the situation in the two-parameter -plane. Either is irrational, that is, , or is rational (or resonant), that is, . In either case an smooth (or normally hyperbolic) invariant circle of around the origin is born for (sufficiently small) in what is known as a Neimark-Sacker bifurcation, provided the genericity condition when is rational. In the latter case one speaks of a weak resonance. The cases are called strong resonances and they involve other, more complicated bifurcations and unfoldings [1, 51]. Resonance tongues are found in the -plane, rooted at weak resonance points along the line . For weak resonances, they are bounded by saddle-node bifurcations of a pair of -periodic locked orbits, one of which is an attractor and the other a saddle of . The one-dimensional unstable manifolds of the saddle -periodic orbit meet at the attracting -periodic orbit to form the smooth invariant curve. In between the resonance tongues in the -plane there are smooth curves, rooted at the line at points where is irrational, along which the dynamics on the invariant curve is quasiperiodic, that is, iterates of fill the invariant curve densely. See, for example, [2, 33, 51] for further details.
The above results, including the analysis of the strong resonance, have been obtained by a normal form procedure, developed independently by Arnold  and Takens , that results in a -equivariant planar vector field whose time-one map is an approximation of the -th iterate of the Poincaré map . This reduces the problem to the analysis of a vector field in two dimensions, instead of three: -periodic points of are equilibria of , and invariant curves of are periodic orbits of .
The road towards a theory of folding tori begins with work by Chenciner in [10, 11, 12] and, together with Gasull and Llibre, in  on the unfolding of a degenerate Neimark-Sacker bifurcation of codimension two, where the Neimark-Sacker bifurcation changes from being supercritical to subcritical. This means that the bifurcating torus changes from being attracting to being of saddle type. This bifurcation is now often called a Chenciner bifurcation; it necessarily involves a mechanism of merging and annihilation of these two tori and folding resonance tongues.
The set of further bifurcations associated with a folding resonance tongue, now referred to as a Chenciner bubble, was considered in [12, 13]. Derived is a planar vector field of degree two on a cylinder, which corresponds to a fundamental sector of the -equivariant planar normal vector field mentioned above. Effectively, the rotational -symmetry has been divided out to obtain a planar vector field on a cylinder, that is, with translational symmetry in the -coordinate on the covering space . A two-parameter bifurcation diagram is presented in [12, 13]. It consists of a pair of saddle-node bifurcations (forming the boundary of the resonance tongue), which are connected by curves of homoclinic bifurcations and of saddle-node bifurcations of periodic orbits. These additional curves meet at a special point on a line of symmetry where the vector field is Hamiltonian. This bifurcation diagram therefore represents a minimal set of dynamics but not an overall unfolding.
3.1 Suggested minimal bifurcation diagram of a Chenciner bubble
In  Baesens and MacKay assume an arbitrary weak locking on the torus and that the associated resonance tongue folds in an appropriate parameter plane. They study a more general planar normal form vector field on the cylinder, which includes higher-order terms compared to the one studied by Chenciner. The authors suggest an associated minimal bifurcation diagram for parameter values near the Chenciner bubble, derived with topological argumentations, where the line of Hamiltonian dynamics is unfolded.
Figure 4 reproduces from  this two-parameter bifurcation diagram of the planar vector fields on the cylinder, consisting of a number bifurcations curves and representative phase portraits in open regions and along the bifurcation curves. This figure provides a topological description of the dynamics one can to expect to find inside a Chenciner bubble, and will serve as a useful guide for what one can expect to find in the GZT model in section 4. The phase portraits in figure 4 illustrate a fundamental domain on the covering space: the left and right equilibria need to be identified in order to obtain the phase portrait on the cylinder. Hence, the saddles on the left and right sides are in fact one and the same point on the cylinder, or each others image under the -symmetry. In the centre of the phase portraits are equilibria, which may be stable or unstable. The thicker ovals are contractible periodic orbits and the thicker curves that go across some phase portraits represent periodic orbits around the cylinder.
The top and bottom curves sne in figure 4 are curves of saddle-node bifurcations of equilibria, where a saddle and node are created. They are folded with respect to the vertical direction and represent the boundaries of the resonance tongue. The two curves sne are bridged by several other bifurcation curves. The curve Hopf connects two Takens-Bogdanov bifurcation points, denoted B, one on each of the curves sne. Along Hopf the node at the centre of the phase portrait changes stability and gives birth to a contractible periodic orbit. The dashed curve ns of neutral saddles also connects the two points B; along it the saddle quantity, that is, the sum of the two real eigenvalues at the saddle is zero. The curve ns does not represent a bifurcation, but it still plays a role in how the global bifurcations are organised. Emerging from the Takens-Bogdanov bifurcation points B, are curves chc of contractible homoclinic bifurcations, where the periodic orbit bifurcating from the curve Hopf disappears. The curves chc meet at a point denoted N, where two curves denoted rhc intersect.
Along the curves rhc one finds rotational homoclinic bifurcations, that is, a homoclinic connections that are non-contractible and extend around the cylindrical phase space. Each of the curves rhc starts and ends on opposing curves sne at points labelled Z; these are codimension-two non-central saddle-node homoclinic bifurcations, which mark the change from a saddle-node bifurcation taking place on a periodic orbit (also referred to as SNIPER or SNIC bifurcation), or off it. The intersection point N of the curves chc is a codimension-two point where both non-contractible homoclinic connections occur simultaneously. The resulting homoclinic cycle can be perturbed to contractible homoclinic connections, which explains why N is also the end point of the two curves chc.
The curves rhc are intersected by the curve ns of neutral saddles at points labelled K, which are codimension-two bifurcations where the criticality of the non-contractible homoclinic connections changes; that is, the bifurcating non-contractible periodic orbit changes from attracting to repelling. From the points K emanate curves snp of saddle-node bifurcations of periodic orbits, where a pair of attracting and repelling non-contractible periodic orbits meet and disappear. Note that one of the curves snp also intersects the other curve rhc at the point X, so that three periodic orbits co-exist inside the shaded triangle formed by the points N, K and X.
The bifurcation set in figure 4 represents a consistent bifurcation diagram of the dynamics of a Chenciner bubble, which is robust or structurally stable in the space of two-dimensional vector fields on a cylinder. Since figure 4 takes into account all of the known bifurcations, one also speaks of a minimal bifurcation diagram of a Chenciner bubble. Since the versal unfolding of this phenomenon on the level of planar vector fields is yet unknown, it cannot be ruled out that the curves are arranged differently in the parameter plane and/or further bifurcation curves exist. This is why figure 4 represents one of several possible cases of robust bifurcation diagrams that may exist. For example, in  an extended bifurcation diagram is depicted for the case that the two Takens-Bogdanov bifurcations are oriented differently so that the curves Hopf and ns intersect once to form a figure-eight shape.
3.2 Interpretation in a three-dimensional phase space
The bifurcation set and phase portraits in figure 4 can be interpreted in terms of the two-dimensional Poincaré map and the original vector field in three-dimensions. This constitutes the inverse of the reduction to a two-dimensional approximating vector field and is a standard procedure in bifurcation analysis. The underlying idea is that the map is a perturbation of the composition of the time-one map of the normal form vector field with the action of the -symmetry; for example, see [23, 30]. To begin with, equilibria of the phase portraits in figure 4 are periodic orbits of the same stability in the three-dimensional phase space of , and the curves sne correspond to the curves of saddle-node bifurcations of periodic orbits that form the boundaries of the resonance tongue under consideration. The curve Hopf corresponds to a curve of Neimark-Sacker or torus bifurcations of , and the bifurcating contractible periodic orbit corresponds to an invariant torus around a periodic orbit. As was discussed earlier, the dynamics on the torus may be quasiperiodic or locked, and this is organised by resonance tongues emanating from this torus bifurcation curves at rational values of the frequency ratio; these secondary resonances are not described by the planar normal form .
The curves chc and rhc of homoclinic bifurcations in figure 4, where the stable and unstable manifolds of the saddle equilibrium coincide, translate to a more complicated bifurcation scenario in the three-dimensional vector field . Namely, the two-dimensional stable and unstable manifolds of a periodic orbit of in do not coincide generically. In fact, any curve of homoclinic bifurcation of the planar normal form gives rise (under arbitrarily small generic perturbations) to a pair of curves of first and last homoclinic tangencies between these manifolds. These curves bound a wedged region in parameter space where the two manifolds intersect transversely in a homoclinic tangle. Entering this region corresponds to the break-up of the corresponding attracting or repelling torus. Here the curve rhc corresponds to the break-up of a main torus, while the curve chc corresponds to the break-up of a secondary torus. The respective two curves of homoclinic tangencies have exponential contact at their end points Z and B on the boundary curves of the resonance tongue under consideration.
Finally, the curves snp translate to folding tori in the three-dimensional vector field . As explained above, these “bifurcations” of tori do not actually exist because the tori involved must break up. Thus, in the bifurcation analysis of a single Chenciner bubble one finds necessarily more folding tori with further Chenciner bubbles. More specifically, under perturbation of the time-one map of the planar normal form, the curve snp transforms into a string of Chenciner bubbles associated with secondary (and hence weak) resonances. The result is a hierarchy of Chenciner bubbles occurring on finer and finer scales.
With this interpretation of its relevance for the full dynamics, and in spite of its lack of definiteness, figure 4 provides the most comprehensive picture of the dynamics within a Chenciner bubble in a map of dimension at least two or a vector field of dimension at least three. For example, the investigation in  of a three-dimensional map, describing an adaptively controlled system, identified the (equivalent of) the bifurcation curves sne and Hopf connecting the boundary curves of folding resonance tongues; more specifically, the continuation software MatContM was used to find in a two-parameter plane curves of Neimark-Sacker bifurcations of periodic points that connect the respective curves of saddle-node bifurcations of periodic points bounding the respective resonance tongue. Very recently, the results in  were extended in , again with MatContM, to include the remaining bifurcation curves suggested in .
4 Chenciner bubbles in the GZT model
We now study the bifurcation diagram in a Chenciner bubble of the GZT model (1). To this end, we focus on the resonance tongue in figure 3, because it the lowest-order (weak) resonance tongue with the most pronounced difference between the two saddle-node bifurcation curves that bound it. This helps the exposition since the respective region of the resonance tongue, which we refer to as the Chenciner bubble, is largest; moreover, working with a small period of the locked periodic solutions makes the numerical calculations more manageable. We remark that there are codimension-two Chenciner bifurcations along the curve T, where it switches between subcritical and supercritical. These points are found outside the range shown in figure 3, for example, at ; see .
4.1 Bifurcation diagram in the -plane
Figure 5 shows the bifurcation structure inside the resonance tongue in the -plane near where it folds with respect to , that is, the vertical axis as it does in figure 4. On the bounding curves SN of saddle-node bifurcations of periodic solutions in figure 5 we find Bogdanov-Takens points B from which the curves T of torus bifurcations and NP of neutral saddle periodic orbits emerge; the latter is defined by , where are the two leading Floquet multipliers of the saddle periodic solution. The curves T and NP are calculated by continuation of periodic solutions with DDE-Biftool.
A curve HoC of homoclinic transitions emerges from each of the two points B; the two curves HoC meet at the point N. This point N is the intersection point of two curves HeC of heteroclinic transitions, which each connect to the two opposite curves SN at points Z; note that the two points Z on the left lie outside the range shown in figure 5. From the points K where the curves HeC intersect the curve NP emerge two curves SNT of folding tori. One of the curves SNT crosses T and, hence, intersects the upper curve HoC at the point X. Because the two curves HeC and SNT are very close together, the two points N and X cannot be distinguished on the scale of figure 5.
In this context it is important to recall that, as was discussed in section 3.2, the curves HoC, HeC and SNT figure 5 do not represent smooth bifurcation curves. This means that they cannot be continued with DDE-Biftool. We speak of homoclinic and heteroclinic transitions along HoC and HeC, respectively, because these represent loci of first and last tangencies, with a region of existence of a tangle in between. It turns out that the curves of first and last tangencies are so extremely close together over the range of figure 5 that they cannot be distinguished for all practical purposes. The topological changes when crossing the loci HoC and HeC, on the other hand, can be detected by dedicated simulations of (1). In this way, the locations of HoC and HeC in the -plane can be determined. Specifically, we found the curves HoC by tracking the phase portrait in for a sufficient number of fixed values of , and the curves HeC by tracking the phase portrait in for a sufficient number of fixed values of . Similarly, the locus SNT of where tori fold is not a smooth curve, but rather the envelop of higher-order resonance tongues. It is very impractical to find these resonance tongues to determine SNT. Instead, we again employed the tracking of phase portraits in for a sufficient number of fixed values of to determine the two loci SNT; as one may expect for resonances of very high order, the loci SNT in figure 5 are effectively one-dimensional. We remark that determining the loci HoC, HeC and SNT in the -plane is a considerable effort.
Overall, we find very good agreement between the computed bifurcation diagram figure 5 of the GZT model (1) and the suggested bifurcation diagram figure 4 from . Indeed, the different curves of bifurcations and transitions agree, taking into account the interpretation of the planar normal form for the dynamics of periodic orbits and tori of the DDE. Examples of phase portraits and the respective transitions will be shown and discussed in section 4.3–4.4, where we demonstrate how the dynamics of the GZT model (1) changes when the bifurcation diagram in figure 5 is crossed.
4.2 Criticality of the curve T
We now discuss briefly a feature of the bifurcation diagram figure 5 that differs from the bifurcation diagram in figure 4, namely the criticality of the curve T. In  it is assumed that the Hopf bifurcation is supercritical throughout, so that a family of stable limit cycles exists between the curve Hopf and the curve chc. However, simulations of the GZT model (1) reveal no evidence of a stable torus being born along curve T, suggesting that the torus bifurcation is subcritical throughout. To determine its criticality we now compute resonance tongues of locked periodic solutions that emerge from curve T.
Figure 6 is an enlargement of the bifurcation diagram figure 5, where only curves directly related to the criticality of the curve T are displayed. Again, the two curves SN at the very top and bottom of figure 6 are the boundaries of the original resonance tongue. The additional blue curves represent secondary resonance tongues that are rooted on the curve T, at points of , and resonance with respect to the periodic solution. Hence, the locked periodic solutions in the resonance tongues have periods , and with respect to the one year forcing cycle, respectively. Because they are extremely narrow, as we determined numerically, we compute these resonance tongues via the continuation of a single curves of the respective locked periodic solution. The secondary resonance tongues are in the region where the periodic solution is stable. Moreover, the secondary locked periodic solutions have one unstable Floquet multiplier, which means that they exist on a torus of saddle-type. This confirms that T is indeed a curve of subcritical torus bifurcations. As we have found and will discuss next, the tori of saddle-type that are created at the curve T terminate at curve HoC. This is in contrast to the suggested bifurcation diagram in figure 4, where the curves Hopf and chc of contractible homoclinic bifurcations both concern stable periodic orbits.
4.3 Transition through the curves T and HoC
Figure 7 shows how the phase portrait of the GZT model (1) changes, for and decreasing , during the transition through the curves T and HoC. This figure also serves to illustrate how the curve HoC was identified in the DDE. More specifically, we compute periodic solutions and tori and present them in different ways: in projection onto the -plane and as stroboscopic trace in the -plane, where triangles, crosses and squares represent stable, 1-saddle and 2-saddle periodic solutions, respectively. An enlargement of the stroboscopic trace also shows the one-dimensional trace of the unstable manifolds of the 1-saddle periodic solutions; as in  these curves were calculated via integration of initial conditions very close to the 1-saddle periodic solutions along their unstable eigendirections — technically, within a region referred to as a fundamental domain . Throughout, light and dark shades of the same colours distinguish between the two symmetry-related periodic solutions (whose existence was discussed in section 1); see  for more details. To facilitate the comparison, we provide the corresponding sketch of the normal form on the cylinder in the style of figure 4.
For parameters to the right of curve T in figure 5, as in row (a) of figure 7, there exists two symmetry-related 1-saddle and two symmetrically-related 2-saddle periodic solutions; see panel (a1). In the stroboscopic trace in the -plane of panels (a2) and (a3) they appear as two sets each of period-seven points. The traces of the unstable manifolds of the 1-saddles surround the 2-saddles and converge to a small-amplitude stable torus (not shown in figure 7 for clarity and corresponding to the green/yellow region of the maximum maps in figure 3). When the curve T in figure 5 is crossed for decreasing , the 2-saddle periodic solutions become stable and a torus is created. This is illustrated in row (b) of figure 7. Notice how the torus, which we find as a high-period locked periodic solution, surrounds the respective stable periodic solution; see, in particular, panel (b3) and notice that the unstable manifolds of the 1-saddles effectively remain unchanged.
For , shown in row (c) of figure 7, there is evidence that the unstable manifolds form a homoclinic tangle. More specifically, the unstable manifolds in one direction still converge to the smaller torus. However, the unstable manifolds in the other direction have become more complicated. After encircling the attracting periodic points, some of the integrated points from the fundamental domain converge to the smaller torus, some converge towards the attracting periodic points and some accumulate at the 1-saddle points. In fact, this distribution of integrated points corresponds to an unstable manifold that is subject to infinite folding and stretching, as is typical for a homoclinic tangle. Because of the drastic folding and stretching (of the trace) of the unstable manifold, we do not represent the unstable manifold as curves; this would require considerably more sophisticated calculations . For the practical purpose of detecting the locus HoC in the bifurcation diagram in figure 5, the representation of the homoclinic tangle in panel (c3) of figure 7 provides a clear identification of the homoclinic transition. Finally, row (d) shows the situation to the left of HoC, when the unstable manifolds of the 1-saddles have undergone the homoclinic transition HoC and now spiral onto the attracting periodic solution. Note that we do not find evidence specifically for the first and last homoclinic tangencies and conclude that they occur within a very small parameter range on the scale of the Chenciner bubble.
Figure 7 clearly shows that crossing the curves T and HoC in figure 5 produces qualitative changes of the dynamics that correspond to crossing the curves Hopf and chc in figure 4, when taking into account that the curve T is now subcritical; compare panels (a4)–(c4) of figure 7 with the respective phase portraits of figure 4.
4.4 Transition through the curves HeC and SNT
Figure 8 shows similarly how the phase portrait of the GZT model (1) changes during the transition through the curves HeC and SNT for and increasing . Row (a) is for parameter above both curves HeC and SNT in figure 5, where there exists a pair of symmetry-related stable periodic solutions and a pair of symmetry-related 1-saddle periodic solutions; see figure 8(a1) and (a2). As panel (a3) shows, the two branches of the unstable manifolds of the 1-saddles spiral into neighbouring stable points in the trace. In particular, the spiralling means that the attracting torus is not normal hyperbolic. The moment of the heteroclinic transition HeC is illustrated in row (b) of figure 8: one branch of the unstable manifold of each 1-saddle connects to the respective other symmetry-related 1-saddle. The black curve that illustrates these connections in panel (b3) is actually a stable torus; see also panels (b1) and (b2). Since this invariant curve passes exceptionally close to the periodic points, this situation can be considered the moment of the HeC transition for the practical detection of the curve HeC in the bifurcation diagram. In contrast to the HoC transition in section 4.3, we found no evidence of first and last tangencies or tangles. In between the curves HeC and SNT, as in row (c) of figure 8 the torus is smooth and further away from the 1-saddles; one branch of their unstable manifolds accumulates on the torus, while the other branch goes to the nearby attracting periodic solution. At the curve SNT the stable torus disappears by colliding with a saddle torus. Even though the locus SNT is associated with secondary and higher-order Chenciner bubbles, numerically, this appears to happen instantly and can be detected as such to find SNT as a curve in figure 5. The phase portrait is then as shown in row (d) of figure 8. Again, the sketches of the corresponding phase portraits of the planar normal form clearly show that the transition through HeC and SNT in figure 5 is described well by that through the curves rhc and snp in figure 4.
5 Discussion and interpretation for tipping
This paper took a detailed look at folding tori in a DDE model, which represents the interaction between negative delayed feedback and periodic forcing. In particular, this phenomenon involves Chenciner bubbles, where the tori break up and the dynamics becomes complicated. We analysed the dynamics inside the Chenciner bubble of the GZT model by computing curves of bifurcations of periodic solutions and tori by means of continuation and numerical integration; for some curves this required the computation of unstable manifolds of saddle periodic solutions. The bifurcation set we found inside the Chenciner bubble compares very well with the suggested bifurcation diagram from  of a planar normal form.
Folding tori are of particular interest in the context of the phenomenon of (climate or other) tipping. One speaks of a tipping event when a slight parameter variation of the system creates a qualitative or drastic response . Such sudden changes in observed behaviour can be associated with certain bifurcations  . For example, past data and various mathematical models, such as those in [35, 36, 37, 41, 42], provide evidence of bistability with a fold bifurcation in the North Atlantic thermohaline circulation (THC) — the circulation mechanism that drives the Gulf Stream. In past studies only simple bifurcations of equilibria or periodic solutions are considered. Indeed, the standard case of tipping is that of an attracting equilibrium or periodic orbit approaching a fold bifurcation, in which it then disappears. We propose that the passage of a stable torus through the locus of folding tori is also a type of tipping. In light of the fact that at least two frequencies need to be involved, we refer to this kind of tipping as multi-frequency tipping.
The additional complexity of the dynamics we encountered in folding tori may have important implications for tipping events in systems that are driven by and/or feature multiple frequencies. This is certainly the case for climate systems more generally, which are driven by effects including orbital forcing at various time scales, as well as different types of feedback mechanisms; for example, see  and references therein. A critical and intriguing question is whether multi-frequency tipping can be identified in real-world data or models of at least intermediate complexity. In fact, there are already some possible candidates for multi-frequency tipping events. For example, the above mentioned North Atlantic THC is known to exhibit fold bifurcations. In  a time series analysis is used to predict a tipping event induced by an unspecified bifurcation in a THC model of intermediate complexity. Due to the many degrees of freedom in the model, the detected bifurcation may, in fact, be a quasiperiodic bifurcation involving tori. As a further example, time series analysis is used in  to identify, in data taken from tropical Pacific sediment cores, an ancient greenhouse to icehouse tipping event — described as “fold-like” in . Given the many degrees of freedom of real-world climate systems, this tipping event may be most accurately described by folding tori as well.
A particular aspect of multi-frequency tipping is that Chenciner bubbles are crossed which, depending on the frequencies involved, may be sufficiently large to be of relevance. This raises a number of interesting questions: What influence does the additional bifurcation structure inside the Chenciner bubble have on the observed dynamics when tipping is approached? Can this knowledge provide useful precursors to warn of an imminent multi-frequency tipping event? While answering these question will require considerable follow-on work, we now provide a brief example of multi-frequency tipping that highlights the role of the relative time scales involved.
Figure 9 shows two time series of the GZT model (1) for fixed and with drifting linearly from 2.96 to 3.01, that is, across the Chenciner bubble shown in figure 5. The time series in panels (a) and (b) of figure 9 drift across the Chenciner bubble in 2000 and in 100 years, respectively. While doing so, the curves HeC, T and SN are crossed, and this is indicated by the labelled vertical lines in the time series. This corresponds to the transition from multi-frequency dynamics on a large torus to dynamics on a considerably smaller torus.
In figure 9(a), for a very slow drift in , one can clearly associate the observed behaviour in the time series with the quasi-static passage through the different bifurcations. Before HeC is reached, the dynamics is characterised by an amplitude modulation of about 40 years. After passing through HeC, at , the time series appear to settle on the stable periodic solution, which results in a much shorter period of the amplitude modulation. When T is crossed, at , the periodic solution loses stability in a torus bifurcation, but there is no obvious change in the time series because the instability is only very weak. Finally, the boundary SN of the folding resonance tongue is crossed at . As a result, the periodic solution ceases to exist and the time series settles onto quasiperiodic oscillations of a much smaller amplitude. One could argue that, in this scenario, the heteroclinic transition HeC acts as a precursor to the more dramatic transition at SN. In figure 9(b) the drift occurs over only 100 years. There is a notable change of the frequency of the observed oscillation at HeC, which may again be interpreted as a precursor for the transition to much smaller amplitude oscillations when SN is crossed. However, the changes in the behaviour of the time series appear to be more subtle for this faster drift.
A wider understanding of multi-frequency tipping may be of interest in the contexts of both long-term (palaeoclimatology) and short-term (anthropogenic or human-induced) climate change, and likely also in other fields, such as ecology. In ongoing and future work, it will be very interesting to study specific candidate systems with folding tori to try and identify qualitative changes in their time series data that can be associated with Chenciner bubbles.
Asides from the effect of different time scales, or rather the size of the relevant Chenciner bubble in parameter space, there are other factors that could impede the identification of torus break-up in time series data. Parameter uncertainty could be an issue. In the brief case study presented here, we have analysed only a single Chenciner bubble that exists for a small range of parameters. However, the Chenciner bubble presented here is only one of an infinite number of Chenciner bubbles, of varying sizes, that exist along the locus of folding tori. For example, in figure 3 one would expect a string of Chenciner bubbles along the sharp interface between large and small maxima. Moreover, there are of course many different parameter paths through a Chenciner bubble that a quasi-static solution could take, and this will affect the specific dynamics that are encountered.
A further issue in this context is that time series data taken from observations or complicated models will be subject to noise, for example, reflecting high-frequency effects omitted in conceptual models such as the GZT model. Understanding the extent to which noise can conceal the effects of the bifurcation structure of a Chenciner bubble in data is ongoing work. Finally, more complicated models, compared to the GZT model considered here, will generally dependent not only on time but also on space. Instead of analysing time series data for an entire mesh of points in space, the data is generally converted into an observable (for example, mean sea-surface temperature across the Pacific Ocean) or reduced to a small number of critical components by, for example, a principal component analysis in order to identify strong trends in the data. It is an interesting and practical challenge to determine to what extent such averaging or data reductions preserve signatures of torus break-up.
We thank Claire Postlethwaite and Anup Purewal for helpful discussions and Jan Sieber for his assistance with DDE-Biftool, especially regarding the code for the continuation of neutral saddles. We acknowledge the New Zealand eScience Infrastructure for use of their high-performance computing cluster.
-  V. I. Arnold, Loss of stability of self-oscillations close to resonance and versal deformations of equivariant vector fields, Functional Analysis and its Applications, 11 (1977), pp. 85–92.
-  , Geometrical methods in the theory of ordinary differential equations, vol. 250, Springer Science & Business Media, 2012.
-  D. G. Aronson, M. A. Chory, G. R. Hall, and R. P. McGehee, Bifurcations from an invariant circle for two-parameter families of maps of the plane: a computer-assisted study, Communications in Mathematical Physics, 83 (1982), pp. 303–354.
-  C. Baesens and R. S. MacKay, Resonances for weak coupling of the unfolding of a saddle-node periodic orbit with an oscillator, Nonlinearity, 20 (2007), p. 1283.
-  T. Bakri and F. Verhulst, Bifurcations of quasi-periodic dynamics: torus breakdown, Zeitschrift für angewandte Mathematik und Physik, 65 (2014), pp. 1053–1076.
-  M. d. S. Baptista and I. L. Caldas, Dynamics of the two-frequency torus breakdown in the driven double scroll circuit, Physical Review E, 58 (1998), p. 4413.
-  H. W. Broer, Quasi-periodicity in dissipative and conservative systems, Proceedings Symposium Henri Poincaré, Université Libre de Bruxelles, Solvay Institutes, 5 (2004), p. 2.
-  H. W. Broer, C. Simó, and J. C. Tatjer, Towards global models near homoclinic tangencies of dissipative diffeomorphisms, Nonlinearity, 11 (1998), p. 667.
-  R. C. Calleja, A. Humphries, and B. Krauskopf, Resonance phenomena in a scalar delay differential equation with two state-dependent delays, SIAM Journal on Applied Dynamical Systems, 16 (2017), pp. 1474–1513.
-  A. Chenciner, Bifurcations de points fixes elliptiques-i. courbes invariantes, Publications Mathématiques de l’IHÉS, 61 (1985), pp. 67–127.
-  , Bifurcations de points fixes elliptiques. ii. orbites périodiques et ensembles de cantor invariants, Inventiones mathematicae, 80 (1985), pp. 81–106.
-  , Bifurcations de points fixes elliptiques-iii. orbites périodiques de «petites» périodes et élimination résonnante des couples de courbes invariantes, Publications Mathématiques de l’IHÉS, 66 (1987), pp. 5–91.
-  A. Chenciner, A. Gasull, and J. Llibre, Une description complète du portrait de phase d’un modèle d’élimination résonante, Comptes rendus de l’Académie des sciences. Série 1, Mathématique, 305 (1987), pp. 623–626.
-  A. J. Clarke, An introduction to the dynamics of El Niño & the Southern Oscillation, Academic press, 2008.
-  V. Dakos, M. Scheffer, E. H. van Nes, V. Brovkin, V. Petoukhov, and H. Held, Slowing down as an early warning signal for abrupt climate change, Proceedings of the National Academy of Sciences, 105 (2008), pp. 14308–14312.
-  R. D. Driver, Ordinary and delay differential equations, Springer-Verlag New York Inc., 1977.
-  K. Engelborghs, T. Luzyanina, and G. Samaey, DDE-BIFTOOL: a Matlab package for bifurcation analysis of delay differential equations, TW Report, 305 (2000).
-  S. R. Gamiz-Fortis and R. T. Sutton, Quasi-periodic fluctuations in the Greenland–Iceland–Norwegian Seas region in a coupled climate model, Ocean Dynamics, 57 (2007), pp. 541–557.
-  A. Garfinkel, P.-S. Chen, D. O. Walter, H. S. Karagueuzian, B. Kogan, S. J. Evans, M. Karpoukhin, C. Hwang, T. Uchida, M. Gotoh, O. Nwasokwa, P. Sager, and J. N. Weiss, Quasiperiodicity and chaos in cardiac fibrillation., The Journal of clinical investigation, 99 (1997), pp. 305–314.
-  M. Ghil, Quaternary glaciations-theory and observations, in The Sun in Time, 1991, pp. 511–542.
-  M. Ghil, I. Zaliapin, and S. Thompson, A delay differential model of ENSO variability: parametric instability and the distribution of extremes, Nonlinear Processes in Geophysics, 15 (2008), pp. 417–433.
-  W. Govaerts, Y. A. Kuznetsov, H. G. E. Meijer, and N. Neirynck, A study of resonance tongues near a Chenciner bifurcation using MatcontM, in 7th European Nonlinear Dynamics Conference (ENOC 2011), Sapienza Universita di Roma, 2011.
-  J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, vol. 42, Springer Science & Business Media, 2013.
-  J. K. Hale and S. M. V. Lunel, Introduction to Functional Differential Equations, Springer-Verlag New York Inc., 1993.
-  H. Hartmann, S. Becker, and L. King, Quasi-periodicities in Chinese precipitation time series, Theoretical and Applied Climatology, 92 (2008), pp. 155–163.
-  H. Held and T. Kleinen, Detection of climate system bifurcations by degenerate fingerprinting, Geophysical Research Letters, 31 (2004).
-  T. Insperger, Stick balancing with reflex delay in case of parametric forcing, Communications in Nonlinear Science and Numerical Simulation, 16 (2011), pp. 2160–2168.
-  K. Kamiyama, M. Komuro, T. Endo, and K. Aihara, Classification of bifurcations of quasi-periodic solutions using lyapunov bundles, International Journal of Bifurcation and Chaos, 24 (2014), p. 1430034.
-  A. Keane, B. Krauskopf, and C. M. Postlethwaite, Delayed feedback versus seasonal forcing: resonance phenomena in an El Niño Southern Oscillation model, SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 1229–1257.
-  B. Krauskopf, Strong resonances and takens’s utrecht preprint, in Global Analysis of Dynamical Systems: Festschrift dedicated to Floris Takens for his 60th birthday, Bristol: Institute of Physics Publishing, 2001, pp. 89–111.
-  B. Krauskopf and K. Green, Computing unstable manifolds of periodic orbits in delay differential equations, Journal of Computational Physics, 186 (2003), pp. 230–249.
-  A. P. Kuznetsov and Y. V. Sedova, The simplest map with three-frequency quasi-periodicity and quasi-periodic bifurcations, International Journal of Bifurcation and Chaos, 26 (2016), p. 1630019.
-  Y. A. Kuznetsov, Elements of applied bifurcation theory, vol. 112, Springer Science & Business Media, 2013.
-  T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf, and H. J. Schellnhuber, Tipping elements in the Earth’s climate system, Proceedings of the National Academy of Sciences, 105 (2008), pp. 1786–1793.
-  T. M. Lenton, R. J. Myerscough, R. Marsh, V. N. Livina, A. R. Price, and S. J. Cox, Using genie to study a tipping point in the climate system, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 367 (2009), pp. 871–884.
-  S. Manabe and R. Stouffer, Two stable equilibria of a coupled ocean-atmosphere model, Journal of Climate, 1 (1988), pp. 841–866.
-  R. Marsh, A. Yool, T. Lenton, M. Gulamali, N. Edwards, J. Shepherd, M. Krznaric, S. Newhouse, and S. Cox, Bistability of the thermohaline circulation identified through comprehensive 2-parameter sweeps of an efficient climate model, Climate Dynamics, 23 (2004), pp. 761–777.
-  M. Münnich, M. A. Cane, and S. E. Zebiak, A study of self-excited oscillations of the tropical ocean-atmosphere system. Part II: Nonlinear cases, Journal of the Atmospheric Sciences, 48 (1991), pp. 1238–1248.
-  G. Mustafa and A. Ertas, Experimental evidence of quasiperiodicity and its breakdown in the column-pendulum oscillator, Journal of dynamic systems, measurement, and control, 117 (1995), pp. 218–225.
-  N. Neirynck, W. Govaerts, Y. A. Kuznetsov, and H. G. Meijer, Numerical bifurcation analysis of homoclinic orbits embedded in one-dimensional manifolds of maps, ACM Transactions on Mathematical Software (TOMS), 44 (2018), p. 25.
-  S. Rahmstorf, Bifurcations of the atlantic thermohaline circulation in response to changes in the hydrological cycle, Nature, 378 (1995), p. 145.
-  S. Rahmstorf, M. Crucifix, A. Ganopolski, H. Goosse, I. Kamenkovich, R. Knutti, G. Lohmann, R. Marsh, L. A. Mysak, Z. Wang, et al., Thermohaline circulation hysteresis: a model intercomparison, Geophysical Research Letters, 32 (2005).
-  D. Roose and R. Szalai, Continuation and bifurcation analysis of delay differential equations, in Krauskopf, B., Osinga H.M. and Galán-Vioque J. (Eds.), Numerical continuation methods for dynamical systems, Springer, 2007, pp. 359–399.
-  R. Saha, Periodic fluctuations in deep water formation due to sea ice, arXiv preprint arXiv:1502.06054, (2015).
-  V. Semenov, A. Zakharova, Y. Maistrenko, and E. Schöll, Delayed-feedback chimera states: Forced multiclusters and stochastic resonance, EPL (Europhysics Letters), 115 (2016), p. 10005.
-  J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, and D. Roose, DDE-BIFTOOL Manual — Bifurcation analysis of delay differential equations, arXiv preprint arXiv:1406.7144, (2014).
-  T. Sorrentino, C. Quintero-Quiroz, A. Aragoneses, M. Torrent, and C. Masoller, Effects of periodic forcing on the temporally correlated spikes of a semiconductor laser with feedback, Optics express, 23 (2015), pp. 5571–5581.
-  G. Stépán, Retarded dynamical systems: stability and characteristic functions, Longman Scientific & Technical, 1989.
-  T. Stocker, G. Clarke, H. L. Treut, R. Lindzen, V. Meleshko, R. Mugara, T. Palmer, R. Pierrehumbert, P. Sellers, K. Trenberth, and J. Willebrand, Physical climate processes and feedbacks, in Climate Change 2001: The Scientific Basis. Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change, J. Houghton, Y. Ding, D. Griggs, M. Noguer, P. van der Linden, X. Dai, K. Maskell, and C. Johnson, eds., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2001, ch. 2, pp. 417–470.
-  C. Sun, J. Li, and F.-F. Jin, A delayed oscillator model for the quasi-periodic multidecadal variability of the NAO, Climate Dynamics, 45 (2015), pp. 2083–2099.
-  F. Takens, Forced oscillations and bifurcations, in Global Analysis of Dynamical Systems: Festschrift dedicated to Floris Takens for his 60th birthday, Bristol: Institute of Physics Publishing, 2001, pp. 1–62.
-  J. M. T. Thompson and J. Sieber, Predicting climate tipping as a noisy bifurcation: a review, International Journal of Bifurcation and Chaos, 21 (2011), pp. 399–423.
-  E. Tziperman, L. Stone, M. A. Cane, and H. Jarosh, El Niño chaos: Overlapping of resonances between the seasonal cycle and the Pacific ocean-atmosphere oscillator, Science-AAAS-Weekly Paper Edition-including Guide to Scientific Information, 264 (1994), pp. 72–73.
-  R. Vitolo, H. W. Broer, and C. Simó, Quasi-periodic bifurcations of invariant circles in low-dimensional dissipative dynamical systems, Regular and Chaotic Dynamics, 16 (2011), pp. 154–184.
-  M. Vrahatis, H. Isliker, and T. Bountis, Structure and breakdown of invariant tori in a 4-d mapping model of accelerator dynamics, International Journal of Bifurcation and Chaos, 7 (1997), pp. 2707–2722.
-  I. Zaliapin and M. Ghil, A delay differential model of ENSO variability —ï¿½ï¿½ part 2: Phase locking, multiple solutions and dynamics of extrema, Nonlinear Processes in Geophysics, 17 (2010), pp. 123–135.
-  C. Zühlsdorff, T. Hanebuth, and R. Henrich, Persistent quasi-periodic turbidite activity off Saharan Africa and its comparability to orbital and climate cyclicities, Geo-Marine Letters, 28 (2008), pp. 87–95.