The paper studies a novel excitability type where a large excitable response appears when a system’s parameter is varied gradually, or ramped, above some critical rate. This occurs even though there is a (unique) stable quiescent state for any fixed setting of the ramped parameter. We give a necessary and a sufficient condition for the existence of a critical ramping rate in a general class of slow–fast systems with folded slow (critical) manifold. Additionally, we derive an analytical condition for the critical rate by relating the excitability threshold to a canard trajectory through a folded saddle singularity. The general framework is used to explain a potential climate tipping point termed the ‘compost-bomb instability’—an explosive release of soil carbon from peatlands into the atmosphere occurs above some critical rate of global warming even though there is a unique asymptotically stable soil carbon equilibrium for any fixed atmospheric temperature.
An excitable system remains in a quiescent state, if undisturbed, produces a small response to a ‘small’ stimulus but fires a large transient response when the small stimulus exceeds a certain threshold. The notion of excitability was first introduced in biology and physiology in an attempt to understand the spiking behaviour of neurons (Hodgkin & Huxley 1952; FitzHugh 1955, 1961) and their electronic simulators (Nagumo et al. 1962). Soon after this work, excitability was found in chemical reactions (Zaikin & Zhabotinsky 1970; Ruoff 1982). More recently, there have been a number of theoretical and experimental demonstrations of excitability in liquid crystals (Coullet et al. 1994), optical systems including lasers (Dubbeldam et al. 1999; Yacomotti et al. 1999; Tredicce 2000; Wieczorek et al. 2002; Wünsche et al. 2002; Krauskopf et al. 2003; Goulding et al. 2007) and photonic crystals (Yacomotti et al. 2006). A hallmark of excitability is a genuine or apparent discontinuity in the system’s response versus the stimulus strength (FitzHugh 1955; Hoppensteadt & Izhikevich 1997; Doi et al. 1999; Izhikevich 2006). It has become clear that this strongly nonlinear phenomenon is mathematically intriguing and relevant to different fields of science. In this paper, we demonstrate a new form of excitability that is relevant to the stability of peatlands under global warming.
From the original notion of excitability, which is purely phenomenological, one can conclude that an excitable system has to have the following properties:
(P1) A quiescent state.
(P2) An excitability threshold above which the system initially evolves away from the quiescent state giving rise to an excitable response.
(P3) A return mechanism that specifies the type of excitable response and, when the stimulus is off, brings the system back to the quiescent state so it can be excited again.
A typical candidate for a stimulus that perturbs the system from the quiescent state to above the excitability threshold is a fast and large enough change in one of the system parameters, a short impulse, for example (Wünsche et al. 2002). Other possibilities include stochastic fluctuations in the system variable(s) (Lindner et al. 2004). This work studies excitability owing to parameter ramping—a steady, slow and monotonic change in one of the system parameters referred to as the ramped parameter. In the problem considered here, a (unique) quiescent state exists for any fixed setting of the ramped parameter. However, a very large excitable response appears when the parameter is ramped sufficiently fast from one setting to another.
(a) Summary of main results
In §2, we define the excitability phenomena in rigorous terms of dynamical systems, give a brief survey on classification of excitability, and describe two different types of excitability, namely type A and type B. This work focuses on an in-depth analysis of type B excitable models, and §3 describes excitability properties in these models with state jumps.
In §4, we study a general class of type B excitable models with parameter ramping, of which the compost-bomb problem is a representative. We show that if a suitable parameter is ramped in a slow–fast system with an asymptotically stable equilibrium and locally folded critical (slow) manifold, then there may be a critical value of the ramping rate above which an excitable response appears. This result is obtained in two steps using concepts from singular perturbation theory. In the first step, we show that a necessary and sufficient condition for the existence of a critical ramping rate is a folded saddle singularity in the corresponding desingularized system. In the second step, we derive an analytical condition for calculating the critical ramping rate.
The general analysis in §4 is motivated by a need to understand the response of peatlands to global warming (or atmospheric temperature ramping), which represents a potential tipping point in the response of the climate system to anthropogenic forcing (Lenton et al. 2008). It is estimated that peatland soils contain 400–1000 billion tonnes of carbon, which is of the same order of magnitude as the carbon content of the atmosphere. Peat carbon is increased by plant litter and reduced by microbial decomposition in the soil. Peat decomposition is expected to accelerate under global warming, leading to concerns that carbon could be released to the atmosphere, accelerating the rate of carbon dioxide increase and providing a positive feedback to global warming (Cox et al. 2000; Khvorostyanov et al. 2008b). There is also strong empirical evidence of peat instability in response to warm and dry climate anomalies, such as the peatland fires in Russia in summer 2010.
Although many global climate-carbon cycle models predict loss of soil carbon under global warming (Friedlingstein et al. 2006), few properly deal with organic soils such as peats, and all ignore the effects of biochemical heat release associated with microbial decomposition (Khvorostyanov et al. 2008a). In their recent work, Luke & Cox (in press) define the peatland soil system in terms of its vertically integrated soil carbon content, C (kg m−2) and soil temperature, T (°C). Soil carbon is increased by litter fall from plants, Π=1.055 (kg m−2 yr−1), and reduced by microbial decomposition, which depends on the store of carbon, C, and also the temperature sensitivity of the decomposition rate per unit carbon, r(T), 1.1Empirical studies suggest that r(T) is approximately exponential in temperature, such that it increases by a factor of 2 to 5 for each 10°C of warming (Kirschbaum 1995). We choose with r0=0.01 (yr−1) and /10 (°C−1). Soil temperature relaxes back to the prescribed atmospheric temperature, Ta, with a time scale dependent on the thermal inertia, μ=2.5×106 (J m−2 °C−1), and the soil-to-atmosphere heat transfer coefficient, λ=5.049×106 (J yr−1 m−2 °C−1). Soil temperature is increased by microbial respiration, Cr(T), with a constant of proportionality, A=3.9×107 (J kg−1), derived from the enthalpy of the respiration reaction (Thornley 1971). This gives the soil temperature equation 1.2with a small parameter ϵ=μ/A≈0.064 (kg °C−1 m−2). Global warming is approximated by an atmospheric temperature ramp 1.3with a constant rate v in units of °C yr−1. Equations (1.1)–(1.3) are the climate-carbon cycle model with global warming. The system (1.1)–(1.2) has a unique asymptotically stable equilibrium for any fixed atmospheric temperature Ta. However, numerical simulations of equations (1.1)–(1.3) suggest that biochemical heat release could destabilize peatland above some critical rate of global warming, vc. The result is an explosive increase in the soil temperature accompanied by a catastrophic release of peatland soil carbon into the atmosphere, which has been termed the ‘compost-bomb instability’ (Luke & Cox in press). Until now, the reason for the rate dependence of the compost-bomb instability has not been rigorously understood, but the present paper shows that this arises from a novel form of excitability.
In §5, we use the results of §4 to show that for an initial condition at the unique stable equilibrium of equations (1.1)–(1.2) and ϵ small, the climate-carbon cycle model (1.1)–(1.3) undergoes a very large excitable response (the compost-bomb instability) if the global warming rate exceeds the critical value 1.4where T0a is the initial atmospheric temperature and is a dimensionless coefficient.
2. Excitable dynamical systems
To define the excitability phenomenon in the language of dynamical systems, we consider an n-dimensional autonomous dynamical system referred to as a stimulus-free problem 2.1where is the state vector and the (constant) parameter vector. A stimulus is represented by some prescribed time dependence of the parameter vector, , so the corresponding non-autonomous stimulus problem reads 2.2To be able to express the phenomenological properties (P1)–(P3) of equations (2.1) and (2.2) in more rigorous terms, let denote a ball of radius r about a(P(t)) such that a(P) is the only attractor for system (2.1) within Br[a(P)], and define:
A quiescent state is an asymptotically stable state for the stimulus-free problem (2.1).
An excitable response (with respect to suitably chosen δ and σ such that 0<δ≤σ) is a trajectory that starts within Bδ[a(P(0))], leaves Bσ[a(P(t))] for some time, before entering into Bδ[a(P(t))]. As , X(t) may remain within Bδ[a(P(t))] or leave and return to Bδ[a(P(t))] repeatedly.
Given P(t), an excitability threshold is a boundary in the phase space separating initial conditions X(0) that give excitable responses from those that give no excitable responses.
The stimulus-free problem (2.1) is called excitable if it has properties (P1)–(P3).
A typical candidate for a quiescent state is a stable equilibrium but other possibilities include small-amplitude periodic orbits or even ‘small’ chaotic attractors (Marino et al. 2004, 2007; Al-Naimee et al. 2009). An excitable response in definition 2.2 is a trajectory that (repeatedly) moves far enough from a(P(t)) before coming close to a(P(t)), where ‘far enough’ is quantified by the choice of σ. In definition 2.3, if the time-dependent component of P(t), namely , is a solution to an autonomous differential equation, pr can be treated as an additional state variable. Then, the excitability threshold is a boundary in the phase space separating initial conditions (X(0),pr(0)) that give excitable responses from those that give no excitable responses. Unlike a quiescent state, an excitability threshold can be defined in both the stimulus free (2.1) and the stimulus (2.2) problems. However, the threshold may change in the presence of a stimulus and it may depend on the form of P(t).
(a) Classification of excitability
Here, in the spirit of the original work by FitzHugh (1955), we classify excitability by the type of threshold and the resulting increase in the system’s response versus the stimulus strength. Specifically, we distinguish two excitability types. In type A excitability, there is a unique threshold given by the stable manifold of an unstable state, and an increase in the excitable response is discontinuous (related to the singular-point threshold phenomenon of FitzHugh (1955)). In type B excitability, the threshold is not unique, meaning that it weakly depends on the choice of σ in definition 2.2, and an increase in the response is abrupt but continuous (related to the quasi-threshold phenomenon of FitzHugh (1955)). Figure 1 shows examples of type A and B systems that produce excitable responses following a state jump from below to above the excitability threshold. Both examples have the time-scale separation resulting in slow–fast dynamics near a. However, a different shape of the critical manifold (green in figure 1) approximating the slow dynamics, and the fact that the type A system has an unstable (saddle) state near a while the type B system does not, give rise to different types of excitability threshold.
In an example of the type A excitability from figure 1a, the unique excitability threshold for the stimulus-free problem is the stable invariant manifold Ws(s) of the saddle s. For state jumps of magnitude Δ, the plot of the maximum (Euclidean) distance between X(t) and a(P(t)) following a jump versus Δ has a discontinuity defining a unique threshold value for Δ. The return mechanism specified by f(x,p) is less obvious. Previous studies focused on systems near to a homoclinic bifurcation (Kuznetsov 1995). After a single and suitably chosen jump from a to above Ws(s), the system follows the upper branch of Wu(s) and produces a single-pulse response near to a 1-homoclinic bifurcation (Plaza et al. 1997; Ventura et al. 2002; Krauskopf et al. 2003), a n-pulse response near to a n-homoclinic bifurcation (Wieczorek et al. 2002), or an irregular and unpredictable response near to the so-called chaotic Shil’nikov case. An example of the type B excitable system from figure 1b does not have a unique excitability threshold, except in the singular limit , as explained in §3. This system is the main focus of this work and its in-depth analysis is presented in §3 for state jumps and in §4 for parameter ramping.
Here, we consider an excitability problem where the stable state a(P) for system (2.1) exists for any fixed value of the ramped parameter P. More recently, various classifications of excitability have emerged in the context of a different but related problem of excitable bursting where the stable state a(P) for system (2.1) loses stability or disappears at a certain setting of the ramped parameter P (Rinzel & Ermentrout 1989). For example, Hoppensteadt & Izhikevich (1997) and Izhikevich (2006) distinguish between class 1 and class 2 neuronal excitability depending on the bifurcation of the quiescent state and the associated continuous or discontinuous increase in the frequency of bursting. Similarly, Gerstner & Kistler (2002) distinguish between type I and type II excitabilities that are the same as Hoppensteadt & Izhikevich’s class 1 and 2 excitabilities, respectively. More generally, Golubitsky et al. (2001) classify excitable bursters by the codimension of the bifurcation in whose unfolding they first appear. Clearly, class 1 and 2, type I and II as well as Golubitsky’s excitabilities cannot be directly compared to our type A and B excitabilities because they refer to a different dynamical problem. Nonetheless, different excitability problems can be related in the following sense. For example, type A excitability occurs near a saddle-node-homoclinic bifurcation (Kuznetsov 1995) that gives rise to class 1/type I excitability, and type B excitability occurs near a singular Hopf bifurcation followed by a canard explosion (Benoît et al. 1981) that gives rise to class 2/type II excitability (Gerstner & Kistler 2002; Izhikevich 2006).
3. Type B excitability with state jumps
In this section, we fix P and discuss excitability in the (autonomous) stimulus-free problem (2.1) with instantaneous jumps of magnitude Δ in the state vector X. This is a convenient framework for discussing excitability thresholds and the novel return mechanism found in the climate-carbon cycle model (1.1)–(1.2). It is also a good starting point for the analysis of the more complicated problem (2.2) with a time-dependent parameter in §4.
We focus on an example of the type B excitable behaviour shown in figure 1b in a singularly perturbed system of the form 3.1with slow variable , fast variable , and parameter vector . There is a singular perturbation parameter 0<ϵ≪1, meaning that the system evolves on fast (t/ϵ) and slow (t) time scales. It is useful to introduce a critical manifold that approximates the slow dynamics and is the dz/dt=0 isocline when g does not depend on ϵ. Assume that system (3.1) satisfies the following conditions:
(A1) There is an asymptotically stable equilibrium
(A2) ∂g(x,z,P,0)/∂x≠0 for any (x,z)∈S(P) so that the critical manifold can be written as x=h(z,P). Furthermore, near to a(P,ϵ), the critical manifold has a fold transverse to the (slow) x direction, L(P), where ∂h/∂z=0 and ∂2h/∂z2≠0. Differentiating the critical manifold condition, g(h(z,P),z,P,0)=0, with respect to z shows that ∂h/∂z=0 is the same as (∂g(x,z,P,0)/∂z)|x=h(z,P)=0, and ∂2h/∂z2≠0 requires ∂2g(x,z,P,0)/∂z2|x=h(z,P)≠0. Therefore,
(A3) Initial conditions from within in the (x,z) phase space converge to a(P,ϵ) as .
System (3.1) satisfying assumptions (A1)–(A3) has properties (P1)–(P3) so we call this system excitable. In particular, fold L(P) provides an excitability threshold as shown in figure 1b, and g(x,z,P,ϵ) specifies the return mechanism, that is the shape, magnitude and duration of excitable responses. We will now describe the excitability threshold owing to a locally folded critical manifold, and the new return mechanism found in the climate-carbon cycle model (1.1)–(1.2).
(a) Excitability threshold
In contrast to the example from figure 1a, the excitability threshold in the singularly perturbed system (3.1) from figure 1b is not unique. This property is best understood by looking at the limiting problem on different time scales (Arnol’d 1994; Szmolyan & Wechselberger 2001). On the slow time scale t, one obtains the reduced system: 3.2that describes the evolution of the slow variable x on the critical manifold S(P). In figure 1b, S(P) is partitioned into an attracting part Sa(P), a repelling part Sr(P), and the fold point L(P). On the fast time scale τ=t/ϵ, one obtains the layer system: 3.3that describes the evolution of the fast variable z for fixed x, in fast time τ=t/ϵ. Sa(P) and Sr(P) correspond to the stable and unstable equilibrium, respectively, of the layer system. A fold L(P) of S(P) transverse to the slow x direction corresponds to a saddle-node bifurcation for the layer system. Therefore, in the limit there is a unique excitability threshold, (figure 2a), where is the left branch of the centre manifold of the saddle-node equilibrium L(P) for the one-dimensional layer system (3.3).
The geometric singular perturbation theory describes how to paste together the dynamical behaviour of the fast and slow limiting problems to obtain dynamics for the singularly perturbed problem with sufficiently small but non-zero ϵ (Fenichel 1979; Jones 1995). Away from L(P), the normally hyperbolic manifolds Sa(P) and Sr(P) of hyperbolic equilibria for the layer system persist for sufficiently small ϵ as invariant (and still normally hyperbolic) attracting Sa,ϵ(P,ϵ) and repelling Sr,ϵ(P,ϵ) slow manifolds (shown in red in figure 2b) for the singularly perturbed system (3.1). The smaller the ϵ, the closer to L(P) is this persistence valid. In figure 2, the non-hyperbolic equilibrium L(P) for the layer system is the only place where normal hyperbolicity fails on S(P). Near to L(P), Sa,ϵ(P,ϵ) and Sr,ϵ(P,ϵ) typically split, allowing for canard trajectories or ducks that follow the repelling slow manifold for a considerable amount of time before rapidly moving away from Sr,ϵ(P,ϵ) in the fast z direction (figure 2b; Benoît et al. 1981; Benoît 1983; Arnol’d 1994; Krupa & Szmolyan 2001; Szmolyan & Wechselberger 2001). Owing to the ϵ-dependent splitting between Sa,ϵ and Sr,ϵ and the associated continuum of canard trajectories, the unique excitability threshold of the limiting problem does not persist in the singularly perturbed system (3.1). Rather, for different choices of σ in definition 2.2, the excitability threshold is given by a different canard trajectory. One may want to choose σ such that the excitability threshold is the canard trajectory that follows Sr,ϵ(P,ϵ) for the longest time. In contrast to type A excitability, the plot of max[z(t)/za] versus the jump magnitude Δ is continuous giving no unique threshold value for Δ, where za denotes the z component of a(P,ϵ).
We illustrate these effects for the excitable system 3.4with Figure 3 shows the response of system (3.4) to jumps of magnitude Δ in x, with dependence on the singular perturbation parameter ϵ and the order N of the specific polynomial function . While the increase in the response remains distinct as an apparent discontinuity for ϵ≤10−2 and N≥2, its position on the Δ-axis varies noticeably with ϵ≥10−3 and N≤4. The unusually large response max[z(t)/za] in figure 3a is shown to scale as ϵ−1 owing to the novel return mechanism described in the following section.
(b) The return mechanism
The most commonly studied return mechanism is found in variants of the excitable van der Pol and FitzHugh–Nagumo equations (van der Pol 1920; FitzHugh 1961; Nagumo et al. 1962; Doi et al. 1999; Izhikevich 2006) and can be modelled using excitable systems of the form: 3.5Specifically, for Q(z,P)=p7z3+p6z2+p5z+p4 and , the return mechanism is given by a cubic or an N-shaped critical manifold x=−Q(z,P). This is illustrated in figure 4a for Q(z,P)=−(z−1.2)3+z−1.2. Upon perturbing such a system from the stable equilibrium a to above the excitability threshold where , the trajectory rapidly moves away from a towards larger z at a speed ∼O(ϵ−1) until it slows down by the attracting part Sa2 of critical manifold S (figure 4a). Then, the trajectory crosses Sa2 towards smaller x and moves along Sa2 at a speed ∼O(1) until it reaches fold L2. Next, the trajectory rapidly moves towards smaller z, slows down by Sa1, crosses Sa1 towards larger x and approaches a along Sa1. Such an excursion in the phase space results in a square-like excitable response (figure 4b) whose magnitude is given approximately by the horizontal distance between L1 and L2, and whose duration is ∼O(1).
The same type of excitability threshold with a distinctively different return mechanism is found in the climate-carbon cycle model (1.1)–(1.2) and system (3.4). In system (3.4), the critical manifold x=−(p2+p3z)/Q(z,P) grows monotonically with z, then folds and decays monotonically to zero if p4,p6>0, p5≥0 and pj≥0 for j=7,8,…. This is illustrated in figure 4c for Q(z,P)=z3/3!+z2/2!+z+1. Upon perturbing such a system from a to above the excitability threshold where , the trajectory rapidly moves away from a towards larger z but it does not encounter any slow manifold that would abruptly stop its predominantly horizontal motion. Rather, for small enough ϵ, the trajectory converges rapidly to a line x=K−ϵz, where K>0, and moves along this line until it arrives at approximately , where . To see this, assume that the trajectory starts at (z0,x0) at time t0, and P is such that Q(z,P)>Kzn for some K>0 and n≥3 to the right of z0. Then, for fixed x≠0, the term (−1,ϵ−1)xQ(z,P) in the vector field for (dx/dt,dz/dt) will dominate for large enough z, to the extent that away from x=0, trajectories will run parallel to a line x=K−ϵz. As the trajectory approaches x∼0, it sharply turns around by crossing the two nullclines dz/dt=0 and dx/dt=0 (figure 4e). Next, the trajectory moves towards smaller z, slows down by Sa, crosses Sa towards larger x, and approaches a at a speed ∼O(1) along Sa. Such an excursion in the phase space results in an excitable response that is much shorter and has a much larger amplitude when compared with the response in systems (3.5). As , the response amplitude becomes where K depends linearly on x0. This result explains the ϵ−1 scaling of the amplitude of excitable responses plotted in figure 3a. The time at which maximum is reached can be estimated using the time to blow-up for the singular system and will be .
4. Type B excitability with parameter ramping
This section discusses excitability in the (non-autonomous) stimulus problem (2.2), where one component of the parameter vector P(t), namely , is ramped at a rate v, 4.1and all the remaining components of P(t), including and a small parameter , do not vary in time. Henceforth, we refer to pr(t) as the ramped parameter, to the stimulus-free problem (2.1) as the unramped system, and to the stimulus problem (2.2) as the ramped system. In addition to assumptions (A1)–(A3), we make two extra assumptions:
(A4) v is constant and, without any loss of generality, v>0.
(A5) In the unramped system (2.1), given by setting v=0, a stable equilibrium a(P) and excitability threshold, e.g. due to a fold L(P) in type B systems from figure 1b, exist for any fixed setting of the ramped parameter pr(t).
Whether the ramped system (2.2) produces an excitable response depends on the initial condition X(0) and on pr(t). This brings us to an interesting question: what is the excitability threshold in ramped systems that preserve a quiescent state and, given X(0) and pr(0), what is the critical rate, vc, required to trigger an excitable response? In the context of the climate-carbon cycle problem with global warming (1.1)–(1.3), this question has relevance to the notion of dangerous rates of climate change.
We focus on type B excitable systems from figure 1b with parameter ramping. If v is small compared with the fast time scale of system (3.1), that is v≪ϵ−1, then pr(t) becomes an additional slow variable and one has to consider a singularly perturbed problem 4.2with two slow variables and one fast variable . System (4.2) has no equilibrium points when v≠0. Rather, a one-dimensional manifold indicates the position of the stable equilibrium of the unramped system (3.1) for different fixed settings of pr(t) in the (x,z,pr) phase space.
In §3, we showed that the excitability threshold in the unramped system (3.1) is unique only in the limiting case . Guided by this result we seek answers to the questions about the excitability threshold and critical rate in the ramped system (4.2) by analysing the ‘slow’ reduced system given by 4.3With q denoting f or g and k denoting x, z or pr, we simplify notation by introducing In the (x,z,pr) phase space, the reduced system evolves on a folded two-dimensional critical manifold Near to a(p,ϵ), S(p) is partitioned into the attracting part Sa(p), repelling part Sr(p), and a one-dimensional fold transverse to the (slow) x direction which follows from assumptions (A2) and (A5). Furthermore, S(p) is a graph over pr and z, 4.4which follows from in assumption (A2) and the implicit function theorem.
Since excitable responses appear in the fast variable z (figure 4c–d), we need to study the evolution of z on S(p) in slow time t, and this is obtained by differentiating the algebraic constraint in equations (4.3) with respect to t, 4.5Then, the reduced system (4.3) becomes 4.6Furthermore, it is useful to project (4.6) onto the (z,pr) plane using equation (4.4) to get the projected reduced system 4.7Note that system (4.7) is singular at a fold, L(p), where . As a result, z(t) blows up (diverges off to infinity in finite time t) when trajectories reach typical points on L(p). However, trajectories that approach special points on L(p) where in a (hyperbolic) direction so that dz/dt remains finite may cross a fold with finite speed (Benoît 1983). Such special points are called folded singularities (Takens 1976; Arnol’d 1994), pseudo singularities (Argemi 1978) or canard points (Szmolyan & Wechselberger 2001). In graphical terms, folded singularities are intersection points between one-dimensional manifolds L(p) and in the (x,z,pr) phase space. To eliminate certain degenerate cases, we assume that
(A6) is a graph over pr and can be expressed as z=u(pr).
The associated trajectories that cross from Sa(p) to Sr(p) via F(p,v) are called singular canards. Folded singularities are absolutely essential to the understanding of the excitability threshold in slowly ramped systems that preserve stable equilibrium and, to make progress, we need to analyse the flow near folded singularities of the projected reduced system (4.7). The analysis is greatly facilitated by the following scaling also known as desingularization (Dumortier & Roussarie 1996; Krupa & Szmolyan 2001), 4.8which preserves the direction of time on Sa(p) but reverses it on Sr(p). The scaling gives the desingularized system, 4.9Clearly, a folded singularity of the projected reduced system (4.7) is a regular equilibrium of the desingularized system (4.9). This means that we can study the phase portrait of the desingularized system (4.9) in a neighbourhood of F(p,v) using standard techniques. Then, we obtain the phase portrait of the projected reduced system (4.7) simply by changing the direction of time on Sr(p) in the phase portrait of (4.9). In general, folded singularities are classified as folded nodes, folded foci, folded saddles and folded saddle-nodes, based on their classification as equilibria of the desingularized system. The relevant problem of possible types of folded singularities and singular canards in with a folded two-dimensional critical manifold was addressed by Szmolyan & Wechselberger (2001).
In the phase portrait of the projected reduced system (4.7), one should be able to identify up to three different types of initial condition within Sa(p). The first type consists of initial conditions for which trajectories approach L(p) and then blow-up along a fast direction. They represent initial states of the ramped system that lead to excitable responses. The second type consists of initial conditions for which trajectories remain within Sa(p) for all time. If these trajectories additionally approach and remain near a(p,ϵ) as , we say that the ramped system adiabatically follows the stable but changing equilibrium of the unramped system. The third type is initial conditions that give rise to canard trajectories crossing from Sa(p) to Sr(p) via the folded singularity F(p,v). According to definition 2.3, a boundary separating regions of the first and second type is the excitability threshold within Sa(p).
(a) Necessary and sufficient condition for critical rate
Given assumptions (A1)–(A6), a phase portrait of the projected reduced system (4.7) has the following properties:
a(p,ϵ)∩L(p)=∅ or a(p,ϵ)⊂Sa(p) (an intersection between a(p,ϵ) and L(p) would correspond to a Hopf bifurcation (Benoît et al. 1981) of the stable equilibrium for the unramped system (3.1) and violate (A5)).
contains no folds transverse to the x direction nor self-intersections by assumption (A6).
A folded singularity F(p,v) corresponds to an intersection between L(p) and ; this follows from equations (4.7) and definitions of L(p) and .
The vector field points in the direction of increasing pr for all (z,pr)∈S(p)\(L(p)\F(p,v)); this follows from equations (4.7) and (A7).
If v≠0, given in (A2) and definitions of a(p,ϵ), S(p) and , it follows that intersects a(p,ϵ) only at special points (z,pr), such that . Also, in the limit (figure 5a).
One can check that properties (i)–(v) allow for two qualitatively different phase portraits of the projected reduced system (4.7) without a folded singularity, shown in figure 5b,c, one of the eight different portraits with a single-folded singularity shown in figure 5d–k, or combinations of these.
An excitable portrait requires an excitability threshold defining a critical value of v. In the two simple cases, either remains within Sa(p) and can intersect a(p,ϵ) if (figure 5b), or lies entirely within Sr(p) (figure 5c). Neither case involves an excitability threshold within Sa(p) because all initial conditions from Sa(p) remain in Sa(p) for all time (figure 5b) or reach L(p) and blow-up (figure 5c). More interesting are transition cases between figure 5b and c, because these involve a folded singularity F(p,v) (figure 5d–k). A necessary condition for a folded singularity 4.10requires the ramped parameter pr(t) in the fast dz/dt component of the vector field; see appendix A(a) for details. Guided by the grey-shaded regions indicating subsets of Sa(p) that remain in Sa(p) for all time, we identify only two cases that involve an excitability threshold within Sa(p). These are a folded saddle in phase portrait (figure 5d) and a folded saddle-node in figure 5e. Since a folded saddle-node is not structurally stable (Kuznetsov 1995), it will unfold under typical changes in the parameter vector p into either phase portrait in figure 5b or a pair of folded singularities including folded saddle from figure 5d and folded node from figure 5j. This leaves only one phase portrait involving an excitability threshold that is structurally stable.
Hence, under assumptions (A1)–(A6), the ramped system (4.2) has an excitability threshold defining a critical ramping rate, vc, if the desingularized system (4.9) has an equilibrium F(p,v), and 4.11meaning that F(p,v) is a saddle; see appendix A(b) for a derivation of equation (4.11). Furthermore, in the limit the excitability threshold is unique. Within Sa(p), this threshold is typically given by the singular canard trajectory through a folded saddle (shown in dark blue in figure 5d). In the three-dimensional (x,z,pr) phase space, the excitability threshold consists of three components: (I) the singular canard trajectory within Sa(p) and all trajectories that converge to it, (II) the part of the fold curve L(p) past the folded saddle F(p) and all trajectories that converge to it, and (III) the subset of Sr(p) dividing initial conditions that end up on the grey-shaded part of Sa(p) from those that blow-up (figure 6a). For sufficiently small but non-zero ϵ and away from L(p), the (normally hyperbolic) attracting Sa(p) and repelling Sr(p) parts of the critical manifold S(p) perturb to nearby invariant (and still normally hyperbolic) slow manifolds Sa,ϵ(p,ϵ) and Sr,ϵ(p,ϵ), respectively (Fenichel 1979; Jones 1995). At the fold curve, however, Sa,ϵ(p,ϵ) and Sr,ϵ(p,ϵ) typically split except for where they intersect along a special trajectory that is called a maximal canard (figure 6b). Specifically, Szmolyan & Wechselberger (2001) prove that a singular canard through a folded saddle for the limiting problem perturbs to a maximal canard for sufficiently small ϵ. As the excitability threshold is no longer unique in analogy to the two-dimensional problem in figure 2, one may want to choose σ in definition 2.2 such that the threshold component (I) in the singularly perturbed system is the maximal canard trajectory and all trajectories that converge to it.
(b) A framework for calculating the critical rate
To calculate the critical rate of ramping, vc, consider a phase portrait as in figure 5d with folded saddle at and assume the system to be at at t=0. Then, the critical rate is the value of v at which the (v-dependent) excitability threshold crosses . In most cases, the threshold can be computed only numerically as the stable invariant manifold Ws(F) of the saddle equilibrium F(p,v) for the desingularized system (4.9). This means that vc too has to be obtained numerically. However, if is sufficiently close to F(p,v), the threshold can be approximated by its linearization at F(p,v), that is a straight line through F(p,v) in the direction of the eigenvector w(p,v)=(w1(p,v),w2(p,v))T corresponding to the negative eigenvalue of saddle equilibrium F(p,v). Then, vc can be calculated from the condition for the threshold line to cross : 4.12In some cases, this implicit general condition simplifies to an explicit formula for the critical rate vc. For example, later equation (5.5) gives such a formula for the climate-carbon cycle model (1.1)–(1.3).
5. The compost-bomb instability
In this section, we explain the compost-bomb instability reported in Luke & Cox (in press). Specifically, we derive the critical speed of global warming above which the mechanism of type B excitability causes an abrupt increase in peatland soil temperature that is accompanied by a potentially catastrophic release of peatland carbon into the atmosphere.
To use the general framework developed in §4, we recognize that the climate-carbon cycle model with global warming (1.1)–(1.3) is a singularly perturbed system with two slow variables , one fast variable and singular perturbation parameter ϵ. The slow dynamics of equations (1.1)–(1.3) can be approximated by the reduced flow, obtained by setting ϵ=0, evolving on a two-dimensional critical manifold 5.1that has a unique fold transverse to the (slow) C direction The ramped system (1.1)–(1.3) has no equilibrium points but the unramped system, obtained by setting v=0 in equations (1.1)–(1.3), has a unique equilibrium at We find from linearizing (1.1)–(1.2) that (Ceq,Teq) is asymptotically stable if that is for all . In the phase space of the ramped system (1.1)–(1.3), indicates the position of the unique asymptotically stable equilibrium of the unramped system for different but fixed settings of Ta. It follows that, for the given parameter values, the climate-cycle carbon model (1.1)–(1.3) satisfies requirements (A1)–(A5) for all .
To describe the evolution of the fast variable T on S in slow time t, we set ϵ=0 in equations (1.1)–(1.3) to get the reduced system, differentiate the resulting algebraic equation with respect to t, and rewrite the reduced system as 5.2Then, we use the condition for S in equation (5.1) to obtain the projection of the reduced system onto the (T,Ta) plane 5.3that is singular at the fold L where α(T−Ta)−1=0. This singularity is removed by the following scaling which gives the desingularized system 5.4In the desingularized system (5.4), the isocline intersects the isocline L at a point For the given parameter values, condition (4.11) is satisfied so F is a saddle equilibrium of the desingularized system (5.4) (and, hence, a folded saddle of the reduced system (5.3)) with eigenvalues and the eigenvector w=(w1,w2)T corresponding to the negative eigenvalue s− is given by where
We now have all the necessary ingredients to extract the critical rate of global warming, vc, from the general condition (4.12). Given an initial condition (C0,T0,T0a)∈Sa sufficiently close to F, using equation (4.12), and as far as , the climate-carbon cycle model (1.1)–(1.3) will exhibit a very large excitable response (the compost-bomb instability) if the rate of global warming exceeds the critical value 5.5where Elin accounts for the deviation of the actual stable manifold of F from its linearization at the initial condition (C0,T0,T0a) used in the derivation of equation (4.12), and Eϵ is a correction for non-zero ϵ. In the special but realistic case with the initial condition at the unique stable equilibrium of the unramped system, that is , the critical rate formula (5.5) simplifies to 5.6where ‘≈’ arises from neglecting Elin and Eϵ. Furthermore, if 2B2≫1, equation (5.6) simplifies to 5.7that is independent of the initial conditions C0 and T0. A comparison in figure 7a between numerically calculated vc (solid curve) and its approximation given by equations (5.6) (dashed line) and (5.7) (dotted line) reveals a very good agreement for sufficiently small ϵ. A small discrepancy between the solid curve and the dashed line is expected even as because of Elin. Given a value of vc, it is interesting to calculate the critical time, tc, defined as the time to reach the temperature of 100 (°C). Results in figure 7b suggest that about 20 years of global warming at a constant rate vc≈0.09 (°C yr−1) may already cause spontaneous combustion of peatlands.
In figure 8, we fix the initial condition (red dot) and explain three qualitatively different responses of the climate-carbon cycle model (1.1)–(1.3) for different rates of global warming v. To avoid numerical problems encountered in Luke & Cox (in press) owing to the model stiffness illustrated in figure 4e, we used rather than . For v=0.075 (°C yr−1)<vc (figure 8a), the initial condition is below the excitability threshold approximated by the singular canard trajectory (blue), so the system does not produce any excitable responses. Rather, the (red) trajectory moves towards a straight away, meaning that the ramped system quickly approaches and then adiabatically follows the stable but changing equilibrium of the unramped system. However, as the rate of global warming is increased to v=0.09 (°C yr−1) >vc (figure 8b), the folded singularity F and the excitability threshold shift their position such that the same initial condition is now above the excitability threshold. As a result, the ramped system reaches the fold of the slow manifold and exhibits an explosive increase in the soil temperature, T, associated with a catastrophic release of soil carbon into the atmosphere. After the excitable spike, the system returns to the attracting part of the slow manifold approximated by Sa, finds itself below the excitability threshold, approaches and then adiabatically follows the stable but changing equilibrium a of the unramped system. Clearly, for large enough v, the system can exhibit more than one excitable spike. In figure 8c for v=0.3 (°C yr−1)>vc, when the system returns to the attracting part of the slow manifold following the first excitable spike, it still finds itself above the excitability threshold and produces the second spike only after which the trajectory approaches a. Note that the phenomenon of multi-pulse excitability was first described by Wieczorek et al. (2002) in connection with n-homoclinic orbits to a saddle focus. The multiple-spike excitable response described here manifests in a similar way but has a rather different underlying dynamical mechanism.
This work identifies and analyses a novel excitability type in systems with ramping—a steady, slow and monotonic change in one of the parameters, pr, called the ramped parameter. The unramped system considered here is of slow–fast nature and its (unique) asymptotically stable equilibrium exists for any fixed setting of pr. When pr is ramped sufficiently slowly from one setting to another, the system can adiabatically follow the stable but changing equilibrium. However, very large excitable responses appear when pr is ramped above some critical rate of ramping.
We showed that such an excitable system forms a singularly perturbed problem with at least two slow variables and focused on the case with locally folded critical (slow) manifold. Using concepts from singular perturbation theory, we studied the system dynamics in terms of folded singularities and associated canard trajectories corresponding to the intersection of an attracting and repelling slow manifold. In this way, we uncovered possible phase portraits of the ramped system, identified those that give rise to excitability, and gave a necessary and sufficient condition for the existence of a critical ramping rate. Furthermore, we identified a novel type of excitability threshold related to a singular canard trajectory through a folded saddle, and explained the very large amplitude of excitable responses by revealing a novel slow–fast return mechanism that allows for single- or multiple-spike excitable responses. Finally, we derived a general condition for calculating the critical rate of ramping.
The general analysis was motivated by the recently reported ‘compost-bomb instability’ (Luke & Cox in press)—a potentially catastrophic explosive release of peatland soil carbon into the atmosphere as the greenhouse gas carbon dioxide, which could significantly accelerate anthropogenic global warming (Khvorostyanov et al. 2008b). Such apparent discontinuities in the response of the climate system to forcing are commonly termed ‘climate tipping points’. Previous studies of climate tipping points have tended to focus on bifurcations or the levels of equilibrium global warming beyond which each tipping point is likely to be excited (Lenton et al. 2008; Thompson & Sieber in press). In contrast, we have shown here that there is a general class of dynamical systems, including the climate-carbon cycle model (1.1)–(1.3), which define a dangerous rate rather than a dangerous level per se. We suspect that such rate-dependent tipping points are much more common in the climate system than is typically assumed, and suggest that deriving the associated critical rates of global warming, as we have done here for the ‘compost-bomb instability’, would provide valuable guidance for climate change policy.
S.W. thanks Mark Holland for useful remarks. P.A. thanks Kiyoyuki Tchizawa for helpful discussions about canards in singular systems.
- Received September 17, 2010.
- Accepted October 29, 2010.
- This journal is © 2010 The Royal Society
(a) Folded singularity condition
A point (pr,z)∈S is a folded singularity for system (4.7) iff (pr,z) is an equilibrium for system (4.9) or . By definition of L, (pr,z) is an equilibrium for system (4.9) if (pr,z)∈L. Because fS(z,pr,p)=0 iff (pr,z)∈a by definition of a, and a∩L=∅ by property (i), we have fS(z,pr,p)≠0 if (pr,z)∈L. Also, is required by assumption (A2). Hence, (pr,z) is an equilibrium for system (4.9) iff (pr,z)∈L and .
(b) Folded saddle condition
Linearizing system (4.9) at F gives a 2×2 Jacobian matrix with entries , , , and . If J11J22−J12J21<0, the Jacobian matrix has two real eigenvalues of opposite sign so F is a saddle equilibrium of the desingularized system (4.9) and a folded saddle of the projected reduced system (4.7).