## Abstract

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.

## 1. Introduction

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.* 2008*b*). 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.* 2008*a*). 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 *r*_{0}=0.01 (yr^{−1}) and /10 (^{°}C^{−1}). Soil temperature relaxes back to the prescribed atmospheric temperature, *T*_{a}, with a time scale dependent on the thermal inertia, *μ*=2.5×10^{6} (J m^{−2} ^{°}C^{−1}), and the soil-to-atmosphere heat transfer coefficient, *λ*=5.049×10^{6} (J yr^{−1} m^{−2} ^{°}C^{−1}). Soil temperature is increased by microbial respiration, *Cr*(*T*), with a constant of proportionality, *A*=3.9×10^{7} (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 *T*_{a}. However, numerical simulations of equations (1.1)–(1.3) suggest that biochemical heat release could destabilize peatland above some critical rate of global warming, *v*_{c}. 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 *T*^{0}_{a} 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 *B*_{r}[*a*(*P*)], and define:

### Definition 2.1

*A quiescent state* is an asymptotically stable state for the stimulus-free problem (2.1).

### Definition 2.2

*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.

### Definition 2.3

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.

### Definition 2.4

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, *p*_{r} can be treated as an additional state variable. Then, the *excitability threshold* is a boundary in the phase space separating initial conditions (*X*(0),*p*_{r}(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 1*a*, the unique excitability threshold for the stimulus-free problem is the stable invariant manifold *W*^{s}(*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 *W*^{s}(*s*), the system follows the upper branch of *W*^{u}(*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 1*b* 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 1*b* 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 d*z*/d*t*=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 ∂^{2}*h*/∂*z*^{2}≠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 ∂^{2}*h*/∂*z*^{2}≠0 requires ∂^{2}*g*(*x*,*z*,*P*,0)/∂*z*^{2}|_{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 1*b*, 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 1*a*, the excitability threshold in the singularly perturbed system (3.1) from figure 1*b* 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 1*b*, *S*(*P*) is partitioned into an attracting part *S*_{a}(*P*), a repelling part *S*_{r}(*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*/ϵ. *S*_{a}(*P*) and *S*_{r}(*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 2*a*), 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 *S*_{a}(*P*) and *S*_{r}(*P*) of hyperbolic equilibria for the layer system persist for sufficiently small ϵ as invariant (and still normally hyperbolic) attracting *S*_{a,ϵ}(*P*,ϵ) and repelling *S*_{r,ϵ}(*P*,ϵ) *slow manifolds* (shown in red in figure 2*b*) 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*), *S*_{a,ϵ}(*P*,ϵ) and *S*_{r,ϵ}(*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 *S*_{r,ϵ}(*P*,ϵ) in the fast *z* direction (figure 2*b*; Benoît *et al.* 1981; Benoît 1983; Arnol’d 1994; Krupa & Szmolyan 2001; Szmolyan & Wechselberger 2001). Owing to the ϵ-dependent splitting between *S*_{a,ϵ} and *S*_{r,ϵ} 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 *S*_{r,ϵ}(*P*,ϵ) for the longest time. In contrast to type A excitability, the plot of max[*z*(*t*)/*z*^{a}] versus the jump magnitude **Δ** is continuous giving no unique threshold value for **Δ**, where *z*^{a} 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*)/*z*^{a}] in figure 3*a* 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*)=*p*_{7}*z*^{3}+*p*_{6}*z*^{2}+*p*_{5}*z*+*p*_{4} and , the return mechanism is given by a cubic or an *N*-shaped critical manifold *x*=−*Q*(*z*,*P*). This is illustrated in figure 4*a* 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 *S*_{a2} of critical manifold *S* (figure 4*a*). Then, the trajectory crosses *S*_{a2} towards smaller *x* and moves along *S*_{a2} at a speed ∼*O*(1) until it reaches fold *L*_{2}. Next, the trajectory rapidly moves towards smaller *z*, slows down by *S*_{a1}, crosses *S*_{a1} towards larger *x* and approaches *a* along *S*_{a1}. Such an excursion in the phase space results in a square-like excitable response (figure 4*b*) whose magnitude is given approximately by the horizontal distance between *L*_{1} and *L*_{2},
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*=−(*p*_{2}+*p*_{3}*z*)/*Q*(*z*,*P*) grows monotonically with *z*, then folds and decays monotonically to zero if *p*_{4},*p*_{6}>0, *p*_{5}≥0 and *p*_{j}≥0 for *j*=7,8,…. This is illustrated in figure 4*c* for *Q*(*z*,*P*)=*z*^{3}/3!+*z*^{2}/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 (*z*^{0},*x*^{0}) at time *t*_{0}, and *P* is such that *Q*(*z*,*P*)>*Kz*^{n} for some *K*>0 and *n*≥3 to the right of *z*^{0}. Then, for fixed *x*≠0, the term (−1,ϵ^{−1})*xQ*(*z*,*P*) in the vector field for (*dx*/d*t*,d*z*/d*t*) 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 d*z*/d*t*=0 and *dx*/d*t*=0 (figure 4*e*). Next, the trajectory moves towards smaller *z*, slows down by *S*_{a}, crosses *S*_{a} towards larger *x*, and approaches *a* at a speed ∼*O*(1) along *S*_{a}. 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 *x*^{0}. This result explains the ϵ^{−1} scaling of the amplitude of excitable responses plotted in figure 3*a*. 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 *p*_{r}(*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 1*b*, exist for any fixed setting of the ramped parameter*p*_{r}(*t*).

Whether the ramped system (2.2) produces an excitable response depends on the initial condition *X*(0) and on *p*_{r}(*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 *p*_{r}(0), what is the *critical rate*, *v*_{c}, 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 1*b* with parameter ramping. If *v* is small compared with the fast time scale of system (3.1), that is *v*≪ϵ^{−1}, then *p*_{r}(*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 *p*_{r}(*t*) in the (*x*,*z*,*p*_{r}) 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 *p*_{r}, we simplify notation by introducing
In the (*x*,*z*,*p*_{r}) 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 *S*_{a}(*p*), repelling part *S*_{r}(*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 *p*_{r} and *z*,
4.4which follows from in assumption (A2) and the implicit function theorem.

Since excitable responses appear in the fast variable *z* (figure 4*c*–*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*,*p*_{r}) 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 d*z*/d*t* 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*,*p*_{r}) phase space. To eliminate certain degenerate cases, we assume that

(A6) is a graph over

*p*_{r}and can be expressed as*z*=*u*(*p*_{r}).

The associated trajectories that cross from *S*_{a}(*p*) to *S*_{r}(*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 *S*_{a}(*p*) but reverses it on *S*_{r}(*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 *S*_{r}(*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 *S*_{a}(*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 *S*_{a}(*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 *S*_{a}(*p*) to *S*_{r}(*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 *S*_{a}(*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*,ϵ)⊂*S*_{a}(*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

*p*_{r}for all (*z*,*p*_{r})∈*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*,*p*_{r}), such that . Also, in the limit (figure 5*a*).

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 5*b*,*c*, one of the eight different portraits with a single-folded singularity shown in figure 5*d*–*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 *S*_{a}(*p*) and can intersect *a*(*p*,ϵ) if (figure 5*b*), or lies entirely within *S*_{r}(*p*) (figure 5*c*). Neither case involves an excitability threshold within *S*_{a}(*p*) because all initial conditions from *S*_{a}(*p*) remain in *S*_{a}(*p*) for all time (figure 5*b*) or reach *L*(*p*) and blow-up (figure 5*c*). More interesting are transition cases between figure 5*b* and *c*, because these involve a folded singularity *F*(*p*,*v*) (figure 5*d*–*k*). A necessary condition for a folded singularity
4.10requires the ramped parameter *p*_{r}(*t*) in the fast d*z*/d*t* component of the vector field; see appendix A(*a*) for details. Guided by the grey-shaded regions indicating subsets of *S*_{a}(*p*) that remain in *S*_{a}(*p*) for all time, we identify only two cases that involve an excitability threshold within *S*_{a}(*p*). These are a folded saddle in phase portrait (figure 5*d*) and a folded saddle-node in figure 5*e*. 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 5*b* or a pair of folded singularities including folded saddle from figure 5*d* and folded node from figure 5*j*. 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, *v*_{c}, 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 *S*_{a}(*p*), this threshold is typically given by the singular canard trajectory through a folded saddle (shown in dark blue in figure 5*d*). In the three-dimensional (*x*,*z*,*p*_{r}) phase space, the excitability threshold consists of three components: (I) the singular canard trajectory within *S*_{a}(*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 *S*_{r}(*p*) dividing initial conditions that end up on the grey-shaded part of *S*_{a}(*p*) from those that blow-up (figure 6*a*). For sufficiently small but non-zero ϵ and away from *L*(*p*), the (normally hyperbolic) attracting *S*_{a}(*p*) and repelling *S*_{r}(*p*) parts of the critical manifold *S*(*p*) perturb to nearby invariant (and still normally hyperbolic) slow manifolds *S*_{a,ϵ}(*p*,ϵ) and *S*_{r,ϵ}(*p*,ϵ), respectively (Fenichel 1979; Jones 1995). At the fold curve, however, *S*_{a,ϵ}(*p*,ϵ) and *S*_{r,ϵ}(*p*,ϵ) typically split except for where they intersect along a special trajectory that is called a *maximal canard* (figure 6*b*). 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, *v*_{c}, consider a phase portrait as in figure 5*d* 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 *W*^{s}(*F*) of the saddle equilibrium *F*(*p*,*v*) for the desingularized system (4.9). This means that *v*_{c} 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*)=(*w*_{1}(*p*,*v*),*w*_{2}(*p*,*v*))^{T} corresponding to the negative eigenvalue of saddle equilibrium *F*(*p*,*v*). Then, *v*_{c} 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 *v*_{c}. 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 (*C*^{eq},*T*^{eq}) 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 *T*_{a}. 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*,*T*_{a}) plane
5.3that is singular at the fold *L* where *α*(*T*−*T*_{a})−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*=(*w*_{1},*w*_{2})^{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, *v*_{c}, from the general condition (4.12). Given an initial condition (*C*^{0},*T*^{0},*T*^{0}_{a})∈*S*_{a} 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 *E*_{lin} accounts for the deviation of the actual stable manifold of *F* from its linearization at the initial condition (*C*^{0},*T*^{0},*T*^{0}_{a}) 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 *E*_{lin} and *E*_{ϵ}. Furthermore, if 2*B*^{2}≫1, equation (5.6) simplifies to
5.7that is independent of the initial conditions *C*^{0} and *T*^{0}. A comparison in figure 7*a* between numerically calculated *v*_{c} (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 *E*_{lin}. Given a value of *v*_{c}, it is interesting to calculate the critical time, *t*_{c}, defined as the time to reach the temperature of 100 (^{°}C). Results in figure 7*b* suggest that about 20 years of global warming at a constant rate *v*_{c}≈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 4*e*, we used rather than . For *v*=0.075 (^{°}C yr^{−1})<*v*_{c} (figure 8*a*), 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}) >*v*_{c} (figure 8*b*), 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 *S*_{a}, 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 8*c* for *v*=0.3 (^{°}C yr^{−1})>*v*_{c}, 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.

## 6. Conclusions

This work identifies and analyses a novel excitability type in systems with ramping—a steady, slow and monotonic change in one of the parameters, *p*_{r}, 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 *p*_{r}. When *p*_{r} 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 *p*_{r} 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.* 2008*b*). 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.

## Acknowledgements

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

## Appendix A

#### (a) Folded singularity condition

A point (*p*_{r},*z*)∈*S* is a folded singularity for system (4.7) iff (*p*_{r},*z*) is an equilibrium for system (4.9) or . By definition of *L*, (*p*_{r},*z*) is an equilibrium for system (4.9) if (*p*_{r},*z*)∈*L*. Because *f*^{S}(*z*,*p*_{r},*p*)=0 iff (*p*_{r},*z*)∈*a* by definition of *a*, and *a*∩*L*=∅ by property (i), we have *f*^{S}(*z*,*p*_{r},*p*)≠0 if (*p*_{r},*z*)∈*L*. Also, is required by assumption (A2). Hence, (*p*_{r},*z*) is an equilibrium for system (4.9) iff (*p*_{r},*z*)∈*L* and .

#### (b) Folded saddle condition

Linearizing system (4.9) at *F* gives a 2×2 Jacobian matrix with entries , , , and . If *J*_{11}*J*_{22}−*J*_{12}*J*_{21}<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).