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

0

a

1.2

1.3

a

c

Empirical studies suggest that) is approximately exponential in temperature, such that it increases by a factor of 2 to 5 for each 10C of warming ( Kirschbaum 1995 ). We choosewith=0.01 (yr) and/10 (). Soil temperature relaxes back to the prescribed atmospheric temperature,, with a time scale dependent on the thermal inertia,=2.5×10(J m), and the soil-to-atmosphere heat transfer coefficient,=5.049×10(J yr). Soil temperature is increased by microbial respiration,), with a constant of proportionality,=3.9×10(J kg), derived from the enthalpy of the respiration reaction ( Thornley 1971 ). This gives the soil temperature equationwith a small parameter ϵ=≈0.064 (kg). Global warming is approximated by an atmospheric temperature rampwith a constant ratein units ofC yr. Equations ( 1.1 )–( 1.3 ) are the. The system ( 1.1 )–( 1.2 ) has a unique asymptotically stable equilibrium for any fixed atmospheric temperature. However, numerical simulations of equations ( 1.1 )–( 1.3 ) suggest that biochemical heat release could destabilize peatland above some critical rate of global warming,. 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.4

a

2. Excitable dynamical systems

whereis the initial atmospheric temperature andis a dimensionless coefficient.

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

2.2

r

whereis the state vector andthe (constant) parameter vector. Ais represented by some prescribed time dependence of the parameter vector,, so the corresponding non-autonomousreadsTo be able to express the phenomenological properties (P1)–(P3) of equations ( 2.1 ) and ( 2.2 ) in more rigorous terms, letdenote a ball of radiusabout)) such that) is the only attractor for system ( 2.1 ) within)], 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.

Figure 1. Examples of phase portraits illustrating two different excitability types in system (2.1) with a stable equilibrium a: (a) type A excitability near a saddle point s and (b) type B excitability near a fold L of critical manifold S=S a ∪L∪S r . Shown are trajectories starting at initial conditions below and above the excitability threshold. The excitable response depends on the return mechanism (dashed) specified by the form of f in (a) and the form of g in (b).

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

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

with slow variable, fast variable, and parameter vector. There is a singular perturbation parameter 03.1 ) satisfies the following conditions:

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

a

r

3.3

a

r

Figure 2. Sketches of phase portraits near locally folded critical manifold S=S a ∪ L∪S r (green) illustrating excitability threshold in (a) the reduced system (3.2) and (b) the singularly perturbed system (3.1) with canard trajectories (blue). denotes the centre manifold of the saddle-node equilibrium L for the one-dimensional layer system (3.3), a denotes the stable equilibrium for system (3.1), S a,ϵ and S r,ϵ (red) denote the attracting and repelling part of the slow manifold for system (3.1), respectively.

that describes the evolution of the slow variableon the critical manifold). In figure 1 ) is partitioned into an attracting part), a repelling part), and the fold point). On the fast time scale/ϵ, one obtains thethat describes the evolution of the fast variablefor fixed, in fast time/ϵ.) and) correspond to the stable and unstable equilibrium, respectively, of the layer system. A fold) of) transverse to the slowdirection corresponds to a saddle-node bifurcation for the layer system. Therefore, in the limitthere is a unique excitability threshold,), whereis the left branch of theof the saddle-node equilibrium) 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 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), 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 2b; 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)/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.4

Figure 3. Response of the excitable system (3.4) to jumps in x of magnitude Δ. We used p 1 =0.5, p 2 =0.1, p 3 =−1 and . (a) Responses for N=3 and different ϵ. (b) Responses for ϵ=10−2 and different N=2,3,4, and (from right to left). Note the logarithmic vertical scale in both panels.

(b) The return mechanism

with Figure 3 shows the response of system ( 3.4 ) to jumps of magnitudein, with dependence on the singular perturbation parameter ϵ and the orderof the specific polynomial function. While the increase in the response remains distinct as an apparent discontinuity for ϵ≤10and≥2, its position on the-axis varies noticeably with ϵ≥10and≤4. The unusually large response max[)/] in figure 3 is shown to scale as ϵowing to the novel return mechanism described in the following section.

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

7

6

5

4

a2

a2

a2

2

a1

a1

a1

1

2

Figure 4. Different return mechanisms ((a) and (c)) give rise to different types of excitable response ((b) and (d)). Shown are trajectories starting at initial conditions below (blue) and above (red) the excitability threshold, and the two nullclines (green). (a), (b) are obtained using equation (3.5) with p 1 =0.5, p 2 =p 3 =−1, ϵ= 10−3 and Q(z,P)=−(z−1.2)3+z−1.2. (c)–(e) are obtained using equation (3.4) with p 1 =0.5, p 2 =0.1, p 3 =−1, ϵ=10−3 and Q(z,P)=z3/3!+z2/2!+z+1. (e) expanded view around the right-most turning point of the (red) excitable trajectory from (c). Note the logarithmic vertical scales in (c)–(e).

Specifically, for)=and, the return mechanism is given by a cubic or an-shaped critical manifold=−). This is illustrated in figure 4 for)=−(−1.2)−1.2. Upon perturbing such a system from the stable equilibriumto above the excitability threshold where, the trajectory rapidly moves away fromtowards largerat a speed ∼(ϵ) until it slows down by the attracting partof critical manifold). Then, the trajectory crossestowards smallerand moves alongat a speed ∼(1) until it reaches fold. Next, the trajectory rapidly moves towards smaller, slows down by, crossestowards largerand approachesalong. Such an excursion in the phase space results in a square-like excitable response ( figure 4 ) whose magnitude is given approximately by the horizontal distance betweenandand whose duration is ∼(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 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 t 0 , 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 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

4. Type B excitability with parameter ramping

wheredepends linearly on. This result explains the ϵscaling of the amplitude of excitable responses plotted in figure 3 . The time at which maximum is reached can be estimated using the time to blow-up for the singular system and will be

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

r

(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 p r (t).

and all the remaining components of), includingand a small parameter, do not vary in time. Henceforth, we refer to) as the, to the stimulus-free problem ( 2.1 ) as the, and to the stimulus problem ( 2.2 ) as the. In addition to assumptions (A1)–(A3), we make two extra assumptions:

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 1b 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.2

r

r

with two slow variablesand one fast variable. System ( 4.2 ) has no equilibrium points when≠0. Rather, a one-dimensional manifoldindicates the position of the stable equilibrium of the unramped system ( 3.1 ) for different fixed settings of) in the () 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.3

r

r

a

r

r

4.4

Withdenotingoranddenotingor, we simplify notation by introducingIn the () phase space, the reduced system evolves on a folded two-dimensional critical manifoldNear to,ϵ),) is partitioned into the attracting part), repelling part), and a one-dimensional fold transverse to the (slow)directionwhich follows from assumptions (A2) and (A5). Furthermore,) is a graph overandwhich follows fromin 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.5

4.6

r

4.7

r

(A6) is a graph over p r and can be expressed as z=u(p r ).

Then, the reduced system ( 4.3 ) becomesFurthermore, it is useful to project ( 4.6 ) onto the () plane using equation ( 4.4 ) to get the projected reduced systemNote that system ( 4.7 ) is singular at a fold,), where. As a result,) blows up (diverges off to infinity in finite time) when trajectories reach typical points on). However, trajectories that approach special points on) wherein a (hyperbolic) direction so that d/dremains finite may cross a fold with finite speed ( Benoît 1983 ). Such special pointsare called Arnol’d 1994 ), Argemi 1978 ) or Szmolyan & Wechselberger 2001 ). In graphical terms, folded singularities are intersection points between one-dimensional manifolds) andin the () phase space. To eliminate certain degenerate cases, we assume that

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

a

r

4.9

r

which preserves the direction of time on) but reverses it on). The scaling gives theClearly, 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) using standard techniques. Then, we obtain the phase portrait of the projected reduced system ( 4.7 ) simply by changing the direction of time on) 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 inwith 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 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.

Figure 5. Phase portraits of the projected reduced system system (4.7) satisfying assumptions (A1)–(A6) corresponding to the projections of critical manifold S onto the (z,p r ) plane in the reduced system (4.6). S a ,L and S r denote the attracting part, fold and the repelling part of S, respectively. Initial conditions from S a that remain in S a for all time are shaded in grey. a denotes the stable equilibrium of the unramped system (dashed curve), denotes the isocline of desingularized system (4.9) (green curve) and F denotes folded singularities (black dots). Canard trajectories that cross from S a to S r through F are shaded in light blue and unique singular canards tangent to the strong stable manifold of the equilibrium for the desingularized system are in dark blue.

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 5b), or lies entirely within S r (p) (figure 5c). 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 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.10

r

a

a

a

requires the ramped parameter) in the fast d/dcomponent of the vector field; see appendix A() for details. Guided by the grey-shaded regions indicating subsets of) that remain in) for all time, we identify only two cases that involve an excitability threshold within). These are a folded saddle in phase portrait ( figure 5 ) and a folded saddle-node in figure 5 . Since a folded saddle-node is not structurally stable ( Kuznetsov 1995 ), it will unfold under typical changes in the parameter vectorinto either phase portrait in figure 5 or a pair of folded singularities including folded saddle from figure 5 and folded node from figure 5 . 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.11

a

r

a

r

a

a

r

a,ϵ

r,ϵ

a,ϵ

r,ϵ

Figure 6. A sketch of the phase portrait of the ramped system (4.2) near a folded saddle F. (a) Three components of the excitability threshold (I)–(III) in the limit . Perturbations that take trajectories from the shaded region across (I), (II) or (III) will cause an excitable response. The trajectory through F is the singular canard (dark blue). S a , L and S r denote the attracting part, fold and repelling part of critical manifold S, respectively (cf. figure 5d). (b) Phase portrait for 0<ϵ≪1. The slow attracting manifold, S a,ϵ , and the slow repelling manifold, S r,ϵ , intersect along the maximal canard trajectory (dark blue).

(b) A framework for calculating the critical rate

meaning that) is a saddle; see appendix A() for a derivation of equation ( 4.11 ). Furthermore, in the limitthe excitability threshold is unique. Within), this threshold is typically given by the singular canard trajectory through a folded saddle (shown in dark blue in figure 5 ). In the three-dimensional () phase space, the excitability threshold consists of three components: (I) the singular canard trajectory within) and all trajectories that converge to it, (II) the part of the fold curve) past the folded saddle) and all trajectories that converge to it, and (III) the subset of) dividing initial conditions that end up on the grey-shaded part of) from those that blow-up ( figure 6 ). For sufficiently small but non-zero ϵ and away from), the (normally hyperbolic) attracting) and repelling) parts of the critical manifold) perturb to nearby invariant (and still normally hyperbolic) slow manifolds,ϵ) and,ϵ), respectively ( Fenichel 1979 Jones 1995 ). At the fold curve, however,,ϵ) and,ϵ) typically split except for where they intersect along a special trajectory that is called a). 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 choosein 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.

To calculate the critical rate of ramping, v c , 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 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.12

c

5. The compost-bomb instability

In some cases, this implicit general condition simplifies to an explicit formula for the critical rate. For example, later equation ( 5.5 ) gives such a formula for the climate-carbon cycle model ( 1.1 )–( 1.3 ).

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

a

that has a unique fold transverse to the (slow)directionThe ramped system ( 1.1 )–( 1.3 ) has no equilibrium points but the unramped system, obtained by setting=0 in equations ( 1.1 )–( 1.3 ), has a unique equilibrium atWe find from linearizing ( 1.1 )–( 1.2 ) that () is asymptotically stable ifthat 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. 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.2

a

5.3

a

5.4

1

2

−

Then, we use the condition forin equation ( 5.1 ) to obtain the projection of the reduced system onto the () planethat is singular at the foldwhere)−1=0. This singularity is removed by the following scalingwhich gives the desingularized systemIn the desingularized system ( 5.4 ), theisoclineintersects theisoclineat a pointFor the given parameter values, condition ( 4.11 ) is satisfied sois a saddle equilibrium of the desingularized system ( 5.4 ) (and, hence, a folded saddle of the reduced system ( 5.3 )) with eigenvaluesand the eigenvector=(corresponding to the negative eigenvalueis given bywhere

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 (C0,T0,T0 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.5

lin

a

ϵ

5.6

lin

ϵ

5.7

c

lin

c

c

c

Figure 7. (a) Critical rate of global warming, v c , as a function of ϵ=μ/A with fixed A=3.9×107 (J kg−1) and varied μ, obtained using (solid curve) numerical solutions to equations (1.1)–(1.3) as well as analytical formulae (dashed line) (5.6) and (dotted line) (5.7) derived for . (b) Critical time to reach 100°C, t c , calculated for v=v c +0.001 (°C yr−1). The initial condition is the stable equilibrium of the unramped system with T0 a =0 (°C). The dots indicate ϵ≈0.064 (kg °C−1 m−2) for μ=2.5×106 (J m−2 °C−1).

whereaccounts for the deviation of the actual stable manifold offrom its linearization at the initial condition () used in the derivation of equation ( 4.12 ), andis 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 towhere ‘≈’ arises from neglectingand. Furthermore, if 2≫1, equation ( 5.6 ) simplifies tothat is independent of the initial conditionsand. A comparison in figure 7 between numerically calculated(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 asbecause of. Given a value of, it is interesting to calculate the critical time,, defined as the time to reach the temperature of 100 (C). Results in figure 7 suggest that about 20 years of global warming at a constant rate≈0.09 (C yr) 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)<v c (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) >v c (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 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 8c 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.

Figure 8. Spiky excitable responses in the climate-carbon cycle model with global warming (1.1)–(1.3) with ϵ≈0.064 (kg °C−1 m−2) and . The solution starting at the equilibrium of the unramped system (C,T,T a )=(50,8.15,0) (red dot) is shown as (left column) a trajectory in the (C,T,T a ) phase space and (right column) a response of the soil temperature T in time t. Different rows show (a) sub-threshold response for v= 0.075 (°C yr−1)<v c , (b) single-spike excitable response for v= 0.09 (°C yr−1)>v c and (c) double-spike excitable response for v=0.3 (°C yr−1) >v c . For reference, included are: the two-dimensional critical manifold S=S a ∪L∪S r (grey surface), the unique fold L of S and the unique asymptotically stable equilibrium a of the unramped system (black curves), the folded saddle F (black dot) and the singular canard trajectory (blue) indicating the excitability threshold within S a (for ). Note the logarithmic T-scale.

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

Acknowledgements S.W. thanks Mark Holland for useful remarks. P.A. thanks Kiyoyuki Tchizawa for helpful discussions about canards in singular systems.

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 fS(z,p r ,p)=0 iff (p r ,z)∈a by definition of a, and a∩L=∅ by property (i), we have fS(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).

Footnotes