3.1. Characteristic return time

The notion of characteristic return time in the context of resilience for ecological systems already dates back to May [Reference May64]. A very commonly used version is due to Beddington et al. [Reference Beddington, Free and Lawton8] for discrete dynamical systems and then to Pimm and Lawton [Reference Pimm and Lawton80] for the continuous case. It is indistinctly used under different names, for example, return time [Reference Pimm and Lawton80], characteristic return time [Reference Pimm79], engineering resilience [Reference Gunderson35, Reference Holling41, Reference Standish, Hobbs and Mayfield95] and resilience [Reference Neubert and Caswell75]. A qualitative description of the underlying idea is provided by Pimm [Reference Pimm79] as ‘how fast the variables return towards their equilibrium following a perturbation’, or, more specifically, as the ‘time taken for a perturbation to return to $1/e$ of its initial value’. The definition is motivated by the fact that a trajectory starting in the nearby of a locally stable equilibrium $x^*$ will approach it in a time which is proportional to the reciprocal of the eigenvalue with largest real part for the system linearisation at $x^*$ . The presentation below is however given for the more general case of a locally attractive hyperbolic solution (see Definition 2.6).

Definition 3.1 (Characteristic return time). Consider $f\,:\,{\mathbb{R}}^N\to{\mathbb{R}}^N$ satisfying assumption 1. The characteristic return time $T_R$ of the system $\dot x=f(x)$ for the attractor $\mathcal{A}$ is defined as

\begin{equation*} T_R({\mathcal {A}})= \frac {1}{\widehat \alpha }, \end{equation*}

where

\begin{equation*} \widehat \alpha :=\inf \big \{\alpha \gt 0\ \mid \ |Y(t)Y^{-1}(s)|\le Ke^{-\alpha (t-s)},\quad \textit {for all }t\ge s\big \}, \end{equation*}

and $Y\,:\,{\mathbb{R}}\to{\mathbb{R}}^{N\times N}$ is a fundamental matrix solution of $\dot y ={\text{D}}f\big (\widetilde x(t)\big )y$ .

The definition of characteristic return time motivated the introduction of the following indicator of resilience for a stable hyperbolic equilibrium as the rate of decay

\begin{equation*}EV({\mathcal {A}})=\widehat \alpha .\end{equation*}

This indicator has been widely used and studied for both continuous and discrete systems in the case that ${\mathcal{A}}=x^*$ is hyperbolic and attracting (see e.g. Arnoldi et al. [Reference Arnoldi, Loreau and Haegeman4], DeAngelis et al. [Reference DeAngelis, Bartell and Brenkert22, Reference DeAngelis, Mulholland, Palumbo, Steinman, Huston and Elwood23], Harwell and Ragsdale [Reference Harwell and Ragsdale39], Pimm and Lawton [Reference Pimm and Lawton81], Rooney et al. [Reference Rooney, McCann, Gellner and Moore85], Van Nes and Scheffer [Reference Van Nes and Scheffer101], Vincent and Anderson [Reference Vincent and Anderson102]).

Invariance with respect to change of coordinates

The characteristic return time is invariant with respect to change of basis $z=Qy$ , with $Q$ non-singular. Indeed, considering the principal matrix solutions $U(t,s)$ and $V(t,s)$ of $\dot y=A(t)y$ and $\dot z=QA(t)Q^{-1}y$ , respectively, and $y_0\in{\mathbb{R}}^N$ , we have that $Q^{-1}V(t,s)Qy_0=U(t,s)y_0$ . Therefore, if $\dot y=A(t)y$ has an exponential dichotomy on $\mathbb{R}$ with projector the identity and constants $\alpha \gt 0$ and $K\ge 1$ , then we have that

\begin{equation*} |V(t,s)|=|QU(t,s)Q^{-1}|\le \widetilde Ke^{-\alpha (t-s)}. \end{equation*}

For autonomous systems, this fact becomes even more evident because $A$ and $QAQ^{-1}$ have the same set of eigenvalues.

3.2. Reactivity

Neubert and Caswell [Reference Neubert and Caswell75] proposed and studied different indicators with the aim of capturing the transient behaviour of a trajectory starting in the neighbourhood of a stable equilibrium as it may substantially differ from the asymptotic one. Specifically, the reactivity corresponds to the maximum instantaneous rate at which an asymptotically stable linear homogeneous system responds if the initial condition is taken outside the origin. The reactivity of a nonlinear system $\dot x = f(x)$ in the neighbourhood of a stable hyperbolic equilibrium $x^*$ is obtained through its linearisation at $x^*$ .

Definition 3.3 (Reactivity). Let $\dot y= A(t)y$ , with $A\,:\,{\mathbb{R}}\to{\mathbb{R}}^{N\times N}$ locally integrable, be an asymptotically stable linear homogeneous system, and denote by $y(\cdot,t_0,y_0)$ its unique solution satisfying $y(t_0,t_0,y_0)=y_0$ . We shall call the reactivity of the system at time $t_0\in{\mathbb{R}}$ , the quantity

(3.1) \begin{equation} R_{t_0}=\max _{|y_0|\neq 0} \left ( \frac{1}{|y(t,t_0,y_0)|}\frac{d|y(t,t_0,y_0)|}{dt} \right )\bigg |_{t=t_0}. \end{equation}

The system $\dot y= A(t)y$ is called reactive if there is $t_0\in{\mathbb{R}}$ such that $R_{t_0}\gt 0$ and nonreactive otherwise. A nonlinear system $\dot x =f(x)$ satisfying Assumption 1 is called reactive if there exists a neighbourhood of a locally attractive hyperbolic solution $\widetilde x$ such that $\dot y ={\text{D}}f\big (\widetilde x(t)\big )y$ is reactive.

If a system is reactive in a neighbourhood of a locally attractive hyperbolic solution $\widetilde x$ , some trajectories starting in a neighbourhood of $\widetilde x$ may initially move away from $\widetilde x$ , before converging to it. In other words, the finite time Lyapunov exponents for $\widetilde x$ can be positive. An uninformed guess might relate the short-term behaviour of solutions to the real part of the least stable eigenvalue of $A$ . This is true only when $A$ has a set of orthogonal eigenvectors. In such a case, however, a monotonic decay towards zero characterises all the solutions since the eigenvalues of ${\text{D}}f(x^*)$ determine both the asymptotic and the transient behaviour of the system. In other words, a short-time amplification is a possible effect of nonorthogonality of the eigenvector basis – also called non-normality of $A$ . Neubert and Caswell [Reference Neubert and Caswell75] unveil a relation between the reactivity of a linear homogeneous system $\dot y= Ay$ , $A\in{\mathbb{R}}^{N\times N}$ , and the dominant eigenvalue of the symmetric part of the Toeplitz decomposition of $A$ (recall that every real symmetric matrix is Hermitian, and therefore, its eigenvalues are real). The argument works also for nonautonomous linear homogeneous problems. Note that

(3.2) \begin{equation} \begin{split} \frac{d|y(t,t_0,y_0)|}{dt} &=\frac{d}{dt}\sqrt{y(t,t_0,y_0)^\top y(t,t_0,y_0)}\\ &=\frac{\dot y(t,t_0,y_0)^\top y(t,t_0,y_0)+y(t,t_0,y_0)^\top \dot y(t,t_0,y_0)}{2|y(t,t_0,y_0)|} \\ & =\frac{ y(t,t_0,y_0)^\top \big (A(t)^\top +A(t)\big )y(t,t_0,y_0)}{2|y(t,t_0,y_0)|}. \end{split} \end{equation}

Therefore, from (3.1) we obtain that

\begin{equation*} R_{t_0}=\max _{|y_0|\neq 0}\frac { y_0^\top \big (A(t_0)^\top +A(t_0)\big )y_0}{2|y(t,t_0,y_0)|^2}. \end{equation*}

The term on the right-hand side of the previous formula is also known as a Rayleigh quotient, and its maximum is attained at the largest eigenvalue of the matrix $\big (A(t_0)^\top +A(t_0)\big )/2$ [Reference Horn and Johnson42], that is,

\begin{equation*}R_{t_0}=\lambda _{\text{dom}}\big (H(A(t_0)\big )\in {\mathbb {R}},\qquad \text {where } H(A)=\frac {A(t_0)^\top +A(t_0) }{2}.\end{equation*}

In the context of autonomous non-normal linear operators, reactivity is also known as the numerical abscissa of $A$ [Reference Embree and Trefethen27].

Invariance with respect to change of coordinates

In general, the reactivity is not preserved under a change of basis $z=Qy$ , with $Q$ non-singular. If $Q$ is orthogonal, that is ( $Q^{-1}=Q^\top$ ), the eigenvalues of the symmetric parts of $A(t_0)$ and $QA(t_0)Q^{-1}$ are the same and reactivity is conserved.

Example 3.4. In this example, we compare three planar asymptotically stable linear systems $\dot{x}=A_i x,$ where $i=\{1,2,3\}:$

\begin{equation*} A_1= \begin{pmatrix} -2 & \quad 0\\[3pt] 5 & \quad -1 \end{pmatrix},\ \ A_2= \begin{pmatrix} -2 & \quad 1\\[4pt] \sqrt{26} + \sqrt{50} & \quad - \sqrt{26} - \sqrt{50} - 1 \end{pmatrix},\ \ A_3= \begin{pmatrix} -2 & \quad 0\\[4pt] 1 & \quad -1 \end{pmatrix}. \end{equation*}

All three systems have the same value of the dominant eigenvalue and consequently the same characteristic return time given as $T_R(0)=1$ . The systems $A_1$ and $A_2$ are reactive with a positive value of reactivity given by $R_0=(\sqrt{26} - 3)/2\approx 1.05.$ The system $A_3$ is nonreactive with $R_0=(\sqrt{2} - 3)/2\approx -0.79.$ In the first row of Figure 3 , the plots show the time evolution of the magnitudes of trajectories starting from different perturbed initial conditions around the origin.

Figure 3. The transient behaviour of three systems $A_1,A_2,A_3$ from Example 3.4 is investigated. In the first row, the time evolution of the magnitude of four trajectories starting from different initial conditions around the origin is depicted. In the second row, we see the respective amplification envelopes of each system. Even though the systems have the same characteristic return time, and additionally, the systems $A_1$ and $A_2$ have also the same value of reactivity, the indicator amplification envelope and its characteristics $\rho _{\text{max}}$ and $t_{\text{max}}$ are able to capture different transient behaviour.

Although the trajectories of the reactive systems $A_1$ and $A_2$ can exhibit a transient growth in their magnitude, this is not possible for the nonreactive system $A_3$ . Furthermore, this example shows that despite having the same characteristic return time and reactivity indicator values, the transient behaviour of the trajectories in systems $A_1$ and $A_2$ may vary substantially. This motivates the introduction of the so-called amplification envelope (see Subsection 3.3 ). The amplification envelopes of three considered systems are plotted in the second row of Figure 3 .

3.3. Maximal amplification and maximal amplification time

Example 3.4 shows that reactive systems with same reactivity and characteristic return time do not necessarily share the same transient behaviour. This fact led Neubert and Caswell [Reference Neubert and Caswell75] to the definition of amplification envelope, which records the maximal deviation from the attractor for trajectories starting at perturbed initial conditions of fixed norm.

Definition 3.5 (Amplification envelope). The amplification envelope for an asymptotically stable linear system $\dot y= A(t)y$ , with $A\,:\,{\mathbb{R}}\to{\mathbb{R}}^{N\times N}$ locally integrable, is defined as the continuous function

(3.3) \begin{equation} \begin{aligned} \rho \,:\,{\mathbb{R}}^+\times{\mathbb{R}}\to{\mathbb{R}}^+,\quad (t, t_0)\mapsto \rho (t,t_0)=\sup _{ |y_0|\neq 0} \frac{|y(t+t_0,t_0,y_0)|}{|y_0|}, \end{aligned} \end{equation}

where $y(\cdot,t_0,y_0)$ is the unique solution of $\dot y= A(t)y$ satisfying $y(t_0,t_0,y_0)=y_0$ . Therefore, $\rho (t,t_0)$ is in fact equal to the operator norm (induced by the Euclidean norm) of the principal matrix solution $Y(t+t_0,t_0)$ of $\dot y =A(t)y$ for the initial time $t_0$ , that is

\begin{equation*} \rho (t,t_0)=\big \|Y(t+t_0,t_0)\big \|_{op}. \end{equation*}

For a nonlinear system $\dot x =f(x)$ such that $f$ satisfies Assumption 1, the amplification envelope in a neighbourhood of a hyperbolic solution $\widetilde x$ is defined as the amplification envelope of the linear system $\dot y ={\text{D}}f(\widetilde x(t))y$ .

Note in particular that if $A(t)=A\in{\mathbb{R}}^{N\times N}$ for all $t\in{\mathbb{R}}$ , the amplification envelope depends only on the variable $t\in{\mathbb{R}}^+$ , that is

\begin{equation*} \rho (t,t_0)=\rho (t)=\big \|e^{At}\big \|_{op}. \end{equation*}

In such a case, a constant $\mathcal{K}(A)$ , known as the Kreiss constant of $A$ and defined as $\mathcal{K}(A)=\sup _{ Re(z)\gt 0} Re(z)\|(zI-A)^{-1}\|_{op}$ , can be used to identify upper and lower estimates for the maximum amplification over time, using Kreiss Matrix Theorem [Reference Trefethen99]. In particular,

\begin{equation*} \mathcal {K}(A)\le \max _{t\gt 0}\rho (t)\le eN\mathcal {K}(A). \end{equation*}

The exact computation of the Kreiss constant is itself the object of study [Reference Apkarian and Noll3]. A method to address the time-dependent case on a finite time interval requires the use of an augmented Lagrangian whose first variations are set to zero to obtain a set of algebraic-differential equations that can be solved forward and backward in time iteratively [Reference Schmid90].

From the amplification envelope, we can derive the two following indicators of resilience. The maximal amplification corresponds to the maximal magnification over time for trajectories starting in the nearby of an attractor, whereas the maximal amplification time records the occurrence of the maximal amplification.

Definition 3.6 (Maximal amplification and maximal amplification time). Consider the setting of Definition 3.5 . The maximal amplification and maximal amplification time are respectively defined as:

\begin{equation*}\rho _{\text{max}}(t_0)=\max _{t \geq 0}\rho (t,t_0),\qquad t_{\text{max}}(t_0)=\text {argmax}_{t \geq 0}\rho (t,t_0).\end{equation*}

3.4. Intrinsic stochastic and deterministic invariability

The amplification envelope provides an effective tool to describe the local transient behaviour of a system close to an attractor if it is affected by an ‘isolated and impulsive’ perturbation. Arnoldi et al. [Reference Arnoldi, Loreau and Haegeman4] investigate the same transient behaviour when the system is subjected to a time-dependent (random or deterministic) forcing. In the spirit of this paper, we generalise the presentation of [Reference Arnoldi, Loreau and Haegeman4] to hyperbolic solutions (see Definition 2.6) up to the point where this makes sense.

Consider a continuous and bounded function $A\,:\,{\mathbb{R}}\to{\mathbb{R}}^{N\times N}$ such that the linear homogeneous system

(3.4) \begin{equation} \dot y= A(t)y \end{equation}

has an exponential dichotomy on $\mathbb{R}$ with projector the identity (see Definition 2.5). This means that, denoted by $U(t,s)$ the principal matrix solution of (3.4) at time $s\in{\mathbb{R}}$ , one has that there are constants $\alpha \gt 0$ and $K\ge 1$ such that

(3.5) \begin{equation} |U(t,s)|\le Ke^{-\alpha (t-s)},\quad \text{for all }t\ge s. \end{equation}

Firstly, consider the stochastic differential problem with additive white noise

(3.6) \begin{equation} dy=A(t)y\,dt+S\,dW(t) \end{equation}

where $S\in{\mathbb{R}}^{N\times M}$ , and $ W= (W_1,\dots W_M)^\top$ is a vector of $M$ independent Brownian motions. As $t\to \infty$ , the distribution of each strong solution of (3.6) converges to a stationary Gaussian distribution centred at the origin. Moreover, the covariance matrix of the system at time $t\in{\mathbb{R}}$ and initial data at $t_0\in{\mathbb{R}}$ is given by

(3.7) \begin{equation} C(t,t_0)=\int _{t_0}^{t}U(t,s)\Sigma U(t,s)^\top \,ds, \end{equation}

where $\Sigma :=SS^\top$ . Notice also that the differentiation of $C(t,t_0)$ with respect to $t$ shows that it solves the matrix differential equation

(3.8) \begin{equation} \frac{d}{dt}C= A(t)C+CA(t)^\top +\Sigma. \end{equation}

In particular, $C(t,-\infty )$ is the only bounded solution of (3.8) over the whole real line, and (3.5) and (3.7) together imply that (3.8) has a local attractor (see e.g. [[Reference Kloeden and Rasmussen50], Theorem 1.23]) that must thus coincide with $C_*(\Sigma )=\{(t,C(t,-\infty )\mid t\in{\mathbb{R}}\}$ . $C_*(\Sigma )$ shall be called the stationary covariance of the system.

In order to construct an indicator of intrinsic resilience against stochastic perturbations, one looks at the largest stationary response among the possible perturbations of given ‘magnitude’.

Definition 3.8 (Intrinsic stochastic invariability). Consider a continuous dynamical system induced by an ordinary differential equation $\dot x =f(x)$ , $x\in{\mathbb{R}}^N$ , satisfying 1 and the assumptions in Section 2. Moreover, assume that $\widetilde x$ is a locally attractive hyperbolic solution and consider the local attractor ${\mathcal{A}}=\{\widetilde x(t)\mid t\in{\mathbb{R}}\}$ . Fixed a norm $\|\cdot \|$ on the matrix space ${\mathbb{R}}^{N\times N}$ (e.g. the operatorial norm, the Frobenius norm), the stochastic variability $\mathcal{V}_S({\mathcal{A}})$ and the intrinsic stochastic invariability of the attractor $\mathcal{A}$ with respect to $\|\cdot \|$ are respectively defined as

\begin{equation*} \mathcal {V}_S({\mathcal {A}})=\sup _{\substack {\Sigma \ge 0,\,\|\Sigma \|=1\\ t\in {\mathbb {R}}}}\|C(t,-\infty )\|_{op}\qquad and \qquad \mathcal {I}_S({\mathcal {A}})=\frac {1}{2\mathcal {V}_S({\mathcal {A}})}, \end{equation*}

where $C_*(\Sigma )$ is the stationary covariance of (3.6) with $A(t)={\text{D}}f\big (\widetilde x(t)\big )$ .

Remark 3.9. If ${\mathcal{A}}=x^*$ , then $A={\text{D}}f(x^*)$ , and (3.8) does not depend on time. In this case, one can show that the linear operator $\widehat A(C):=AC+CA^\top$ has $N^2$ eigenvalues of the form $\lambda _j+\lambda _i$ , where $\lambda _i,\lambda _j$ are eigenvalues of $A$ , and the stationary covariance matrix $C_*(\Sigma )$ is given by the unique solution of the Lyapunov equation $AC+CA^\top +\Sigma =0$ [[Reference Berglund and Gentz10], Lemma 5.1.2]. Therefore, one has that

\begin{equation*} \mathcal {V}_S({\mathcal {A}})=\sup _{\Sigma \ge 0,\,\|\Sigma \|=1}\|-\widehat A^{-1}(\Sigma )\|. \end{equation*}

On the other hand, one can consider a deterministic forcing of (3.4), that is,

(3.9) \begin{equation} \dot y=A(t)y+g(t), \end{equation}

where $g\,:\,{\mathbb{R}}\to{\mathbb{R}}^N$ is a bounded function. Denoted by $Y(t)$ a fundamental matrix solution of (3.4), note that

\begin{equation*} x(t,g)=Y(t)\int _{-\infty }^tY^{-1}(s)g(s)\,ds,\quad t\in {\mathbb {R}}, \end{equation*}

is a solution of (3.9), which can be regarded as the stationary system response, and its mean square deviation from the origin (which is the global attractor of the homogeneous problem) is

(3.10) \begin{equation} m(g):=\lim _{T\to \infty }\frac{1}{T}\int _0^T|x(s,g)|^2\,ds. \end{equation}

Arnoldi et al. [Reference Arnoldi, Loreau and Haegeman4] suggest to use the largest mean square deviation of $x(t,g)$ from the origin over the possible perturbations $g$ of given norm, as an indicator of resilience of the attractor.

Definition 3.10 (Intrinsic deterministic invariability). Consider a continuous dynamical system induced by an ordinary differential equation $\dot x =f(x)$ , $x\in{\mathbb{R}}^N$ , satisfying 1 and the assumptions in Section 2. Moreover, assume that $\widetilde x$ is a locally attractive hyperbolic solution and consider the local attractor ${\mathcal{A}}=\{\widetilde x(t)\mid t\in{\mathbb{R}}\}$ . The deterministic variability $\mathcal{V}_D({\mathcal{A}})$ and the intrinsic deterministic invariability of the attractor $\mathcal{A}$ are respectively defined as

(3.11) \begin{equation} \mathcal{V}_D({\mathcal{A}})=\sup _{\|g\|_\infty =1}2\sqrt{m(g)},\qquad and \qquad \mathcal{I}_D({\mathcal{A}})=\frac{1}{\mathcal{V}_D({\mathcal{A}})}, \end{equation}

where $m(g)$ is the mean square deviation of the stationary response of (3.9) with $A(t)={\text{D}}f\big (\widetilde x(t)\big )$ .

Remark 3.11. If ${\mathcal{A}}=x^*$ , then $A={\text{D}}f(x^*)$ , and one can use standard frequency analysis to improve the information given by (3.11). This is particularly true when one limits the possible deterministic perturbation to the so-called ‘wide-sense stationary signals’. Indeed, since any deterministic signal can be developed into a sum of harmonic terms, or Fourier modes, and due to the fact that in the linear approximation, the system response to a general perturbation is equal to the sum of the system response to the single-frequency modes, a convexity argument yields that the perturbation generating the largest system response is a single-frequency mode [Reference Arnoldi, Loreau and Haegeman4]. Therefore, when $g$ is a single-frequency periodic forcing one obtains that

\begin{equation*} \mathcal {V}_D({\mathcal {A}})=\sup _{\omega \in {\mathbb {R}}}\|\text {i}\omega -A^{-1}\|, \end{equation*}

where $\omega$ is the forcing frequency of $g$ , $\text{i}$ is the imaginary unit, and $\|\cdot \|$ is the induced matrix operator norm from the norm on $\mathbb{R}^N$ .

When ${\mathcal{A}}=x^*$ , an important result in [Reference Arnoldi, Loreau and Haegeman4] establishes a chain of order for some of the indicators of local resilience presented in this section.

Proposition 3.12. If ${\mathcal{A}}=x^*$ is a hyperbolic stable equilibrium, the following chain of inequalities holds true:

\begin{equation*} -R_0\le \mathcal {I}_S({\mathcal {A}})\le \mathcal {I}_D({\mathcal {A}})\le EV({\mathcal {A}}). \end{equation*}

4.1. Latitude in width, distance to threshold and precariousness

The notion of ‘width’ of a basin of attraction as an indicator of resilience dates back directly to the seminal work by Holling [Reference Holling40]. Loosely speaking, the width of a basin of attraction corresponds to the length of the ‘minimal’ segment crossing through the attractor and intersecting (at its extrema) the boundary of the basin of attraction, if it applies (see illustrative sketch in Figure 4). Although a precise mathematical formulation seems hard to find in the literature, the geometrical representations used in many works (e.g. Peterson et al. [Reference Peterson, Allen and Holling78] and Walker et al. [Reference Walker, Holling, Carpenter and Kinzig103]) permit to precisely formalise the underlying idea.

Figure 4. Three scenarios of a planar dynamical system are depicted. The grey regions represent the basins of attraction for the different attractors in blue. The corresponding boundaries are depicted as a black dashed line. (a) Latitude in width $L_w$ for this equilibrium is given as the length of the red dashed line. Distance to threshold is given as the length of the green line. (b) System with a limit cycle attractor. $L_w$ is given as the length of the red dashed line and distance to threshold $DT$ as the length of the green line. Note that the lines have different locations. (c) Example of an attractor with a basin of attraction that has infinite latitude in width indicator ${L_w= + \infty }.$ Note the distance to threshold given as the length of the green line is finite.

Definition 4.1 (Latitude in width). Consider the (possibly empty) set

(4.1) \begin{equation} S= \big \{(y,z)\in{\mathbb{R}}^{2N}\,\big |\,y, z \in \partial{\mathcal{B}(\mathcal{A})}, \textit{ and }\alpha y+ (1-\alpha )z \in{\mathcal{A}}\textit{ for some } \alpha \in (0,1) \big \}. \end{equation}

The latitude in width of an attractor ${\mathcal{A}}\subset{\mathbb{R}}^N$ for a continuous flow $\phi \,:\,{\mathbb{R}}\times{\mathbb{R}}^N\to{\mathbb{R}}^N$ is defined as

\begin{equation*} L_w({\mathcal {A}})=\begin {cases} \infty,&\textit {if }S=\emptyset \\ \inf _{(y, z) \in S}|y-z|,&\textit {otherwise}. \end {cases} \end{equation*}

If $L_w({\mathcal{A}})\lt \infty$ , the inferior in the previous formula is in fact a minimum as shown in the next result.

Proposition 4.2. Under the notation and assumptions of Definition 4.1 and if $L_w \lt +\infty$ , we have:

\begin{equation*}L_w({\mathcal {A}}) = \inf _{(y, z) \in S} |y-z| = \min _{(y, z) \in S} |y-z|.\end{equation*}

Proof. By definition $\mathcal{A}$ is a compact subset of ${\mathbb{R}}^N$ . Therefore, there is $\rho \gt L_w({\mathcal{A}})\gt 0$ such that for all $a_1,a_2 \in{\mathcal{A}}$ , $|a_1-a_2|\lt \rho$ . Moreover, notice that $\partial{\mathcal{B}(\mathcal{A})}$ is nonempty since $L_w({\mathcal{A}})\lt \rho \lt \infty$ . By the compactness of $\mathcal{A}$ , there are $a_1,\dots, a_n\in{\mathcal{A}}$ such that

\begin{equation*} {\mathcal {A}}\subset \bigcup _{i=1}^n \overline {B}_{\rho }(a_i)\,=\!:\,P. \end{equation*}

In particular, notice that if $y,z\in S$ and $|y-z|\lt \rho$ , then $(y,z)\in P\times P$ , which is a compact subset of ${\mathbb{R}}^{2N}$ . Therefore, by the continuity of the norm in ${\mathbb{R}}^N$ and Weierstrass Theorem, we immediately obtain the result.

Already in Peterson et al. [Reference Peterson, Allen and Holling78], while still employing the idea of width of the basin of attraction (under the name of ecological resilience), a focus was aimed at understanding the smallest perturbation of initial conditions able to drive the system away from the original attractor (or equivalently the largest perturbation of initial conditions for which the systems returns to the attractor). This approach, qualitatively, yet unequivocally presented by Beisner et al. [Reference Beisner, Dent and Carpenter9], emphasises the importance of the points of minimal distance between an attractor and the boundary of its basin of attraction (in some cases called threshold or separatrix). The definition has the advantage of revealing that an attractor with infinite latitude in width (see Definition 4.1) is not necessarily ‘very resilient’ (see Figure 5). This idea has been used consistently ever since for example in Lundström and Adainpää [Reference Lundström and Aidanpää62], Lundström [Reference Lundström61], Kerswell et al. [Reference Kerswell, Pringle and Willis47], Klinshov et al. [Reference Klinshov, Nekorkin and Kurths48], Mitra et al. [Reference Mitra, Kurths and Donner74] and Fassoni and Yang [Reference Fassoni and Yang29] where it is called precariousness (on the subject see also Definition 4.7 and Proposition 4.8 below).

Figure 5. Representation of the dynamics induced by the planar system $\dot{x}_1=-x_1+10x_2$ , $\dot{x}_2=x_2(10\exp \big(\frac{-x_1^2}{100}\big)-x_2)(x_2-1)$ . The system has three equilibria $X_0$ , $X_s$ and $X_1$ . The stable and unstable manifolds of the saddle-node $X_s$ are depicted in red. The stable manifold of $X_s$ is the separatrix between the basins of attraction of $X_0$ and $X_1$ , respectively painted in green and white. A qualitative behaviour of the system can be deduced via the vector field (blue arrows) and a few trajectories in solid blue. Both stable equilibria have an infinite latitude in width, while their distance to threshold can be made as a small as wished through a suitable scaling. The resilience of this system has been thoroughly analysed by Kerswell et al. [Reference Kerswell, Pringle and Willis47].

Definition 4.3 (Distance to threshold). Consider the (possibly empty) set

(4.2) \begin{equation} S=\{(a,y)\in{\mathbb{R}}^{2N}\mid a \in{\mathcal{A}}, y \in \partial{\mathcal{B}(\mathcal{A})} \}. \end{equation}

The distance to threshold for an attractor ${\mathcal{A}}\subset{\mathbb{R}}^N$ for a continuous flow $\phi \,:\,{\mathbb{R}}\times{\mathbb{R}}^N\to{\mathbb{R}}^N$ is defined as

\begin{equation*} DT({\mathcal {A}}) =\begin {cases} \infty,&\textit {if }S=\emptyset \\ \inf _{(a, y) \in S}|a-y|,&\textit {otherwise}. \end {cases} \end{equation*}

When dealing with a concrete model of a real phenomenon, it is possible to tune the previous definitions according to the available information on the phenomenon. In particular, if one can anticipate the set of possible perturbed initial conditions, the calculation of the latitude in width and the distance to threshold can be restricted to a suitable subset $C$ of the phase space called region of interest. For example, such an adjustment on the distance to threshold has been explicitly presented in [Reference Meyer68] where the set $S$ in (4.2) is substituted by

(4.3) \begin{equation} S_C=\{(a,y)\in{\mathbb{R}}^{2N}\mid a \in{\mathcal{A}}, y \in \partial{\mathcal{B}(\mathcal{A})}\cap C\}, \end{equation}

for some $C\subset{\mathbb{R}}^N$ . Likewise, the latitude in width can be modified with respect to the region of interest. Then, for some fixed $C\subset{\mathbb{R}}^N$ , the set $S$ (4.1) is changed for

\begin{equation*} S_C=\{(y,z)\,:\, (y,z) \in S \text { and } y \in C\}. \end{equation*}

Example 4.5.

(4.4) \begin{equation} \begin{aligned} \dot{x}&= x(x-1/\varepsilon )(x+\varepsilon ) \quad \quad {where }\quad x \in{\mathbb{R}}, \varepsilon \in (0,1] \end{aligned} \end{equation}

The dynamical system induced by (4.4) has a stable equilibrium at $x=0$ , which we consider as the relevant attractor, and two unstable fixed points at $-\varepsilon$ and $1/\varepsilon .$ See Figure 6 for a sketch of the phase space. Notice that when $\varepsilon$ becomes smaller, the distance to threshold $DT(0)=\varepsilon$ decreases, but the latitude in width $L_w(0)=1/\varepsilon +\varepsilon$ increases. If the set of anticipated perturbed initial conditions $C={\mathbb{R}}^+$ is considered, the distance to threshold becomes $DT_{{\mathbb{R}}^+}(0)=1/\varepsilon,$ since we have the set $S_C=\{(0,1/\varepsilon )\}$ (see (4.3)), while the latitude in width $L_{w,{\mathbb{R}}^+}(0)$ remains the same. On the other hand, if we fix $C=(0,\varepsilon ),$ the set $S_C$ becomes empty for both distance to threshold and latitude in width and therefore $DT_{(0,\varepsilon )}(0)=L_{w,(0,\varepsilon )}(0)=+\infty .$

Figure 6. Sketch of the phase space for the dynamical system induced by (4.4) and of the distance to threshold and latitude in width of the attractor ${\mathcal{A}}=\{0\}$ . Additionally, we can see that when we specify the anticipated perturbations as ${\mathbb{R}}^+,$ the distance to threshold changes, whereas, the latitude in width remains the same (for details, see Example 4.5).

Example 4.6. Let us consider the example presented in Remark 2.2 (see also Figure 7 ). Depending on the different choice of the local attractor we may obtain different values of the distance to threshold $DT.$ If we consider only the limit cycle of radius $1$ as the attractor ${\mathcal{A}}_1,$ the unstable equilibrium at the origin is part of the basin’s boundary $\partial \mathcal{B}(\mathcal{A}_1).$ Therefore, the distance to threshold $DT({\mathcal{A}}_1)=1$ . However, if we consider as a local attractor the set of points in the limit cycle and the origin, ${\mathcal{A}}_2={\mathcal{A}}_1\cup \{(0,0)\},$ or the whole closed ball of radius one ${\mathcal{A}}_3=B_1$ , we obtain $DT({\mathcal{A}}_2)=DT({\mathcal{A}}_3)=2.$ It is obvious that in terms of actual resilience of the system’s state against a perturbation the origin plays a negligible role.

Figure 7. Sketch of the phase space for the dynamical system induced by (2.2) and representation of the distance to threshold depending on the choice of the local attractor (see Remark 2.2). For the local attractor ${\mathcal{A}}_1$ consisting of the points in the limit cycle of radius $1$ , $DT({\mathcal{A}}_1)=1$ . For the local attractor ${\mathcal{A}}_2={\mathcal{A}}_1\cup \{(0,0)\}$ , $DT({\mathcal{A}}_2)=2.$ The choice of the attractor ${\mathcal{A}}_2$ admits to consider only distances to the ‘significant’ parts of the boundary.

Instead of focusing on an attractor and its basin, in practice it is convenient to look at the current state of the system within the basin. The idea of studying the minimal distance from the boundary of the basin has qualitatively presented by Walker et al. [Reference Walker, Holling, Carpenter and Kinzig103], where it is described as only one of four components of ecological resilience together with the latitude in width (see Definition 4.1), the resistance of a gradient system (see Definition 5.3) and the cross-scale interaction among the previous components and eventual subsystems. We push the idea in [Reference Walker, Holling, Carpenter and Kinzig103] a bit forward and propose a mathematical definition which probes also the points outside the basin of attraction. This is not only a reasonable generalisation, as it also allows to quantify the difficulty for a point outside the basin to enter it after a given perturbation, but it also allows us to later establish a connection with a recently presented indicator (see Definition 5.5).

Definition 4.7 (Precariousness). The function $P_{\mathcal{A}}\,:\,{\mathbb{R}}^N\to{\mathbb{R}}\cup \{\infty \}$ defined by

\begin{equation*} x_0\mapsto P_{\mathcal {A}}(x_0)=\begin {cases} \infty &\textit {if }\partial {\mathcal{B}(\mathcal{A})}=\varnothing,\\[4pt] \displaystyle \inf _{y\in \partial {\mathcal{B}(\mathcal{A})}}|x_0-y|&\textit {if }\partial {\mathcal{B}(\mathcal{A})}\neq \varnothing \textit { and } x_0\in {\mathcal{B}(\mathcal{A})}\\[4pt] \displaystyle -\inf _{y\in \partial {\mathcal{B}(\mathcal{A})}}|x_0-y|&\textit {if }\partial {\mathcal{B}(\mathcal{A})}\neq \varnothing \textit { and } x_0\notin {\mathcal{B}(\mathcal{A})}\\ \end {cases} \end{equation*}

is called the precariousness of the point $x_0\in{\mathbb{R}}^N$ with respect to the basin of attraction $\mathcal{B}(\mathcal{A})$ .

The following simple proposition establishes a relation between the distance to threshold and the asymptotic value of precariousness.

Proposition 4.8. Consider a continuous dynamical system $\phi$ on ${\mathbb{R}}^N$ . If $\partial{\mathcal{B}(\mathcal{A})}$ is nonempty, then

\begin{equation*} \lim _{t\to \infty }P_{\mathcal {A}}\big (\phi (t,x_0)\big )\ge DT({\mathcal {A}}),\quad \textit {for all }x_0\in {\mathcal{B}(\mathcal{A})}, \end{equation*}

and the equality holds if every trajectory in $\mathcal{A}$ is dense.

Proof. The first inequality is a direct consequence of the definition of basin of attraction (see Definition 2.1) and of distance to threshold (see Definition 4.3) that is for any $x_0\in{\mathcal{B}(\mathcal{A})}$ ,

\begin{align*} \begin{split} \lim _{t\to \infty }P_{\mathcal{A}}\big (\phi (t,x_0)\big )&=\lim _{t\to \infty } \inf _{y\in \partial{\mathcal{B}(\mathcal{A})}}|\phi (t,x_0)-y| =\lim _{t\to \infty } \inf _{\substack{x\in \omega (x_0),\\y\in \partial{\mathcal{B}(\mathcal{A})}}}|\phi (t,x_0)\pm x-y|\\ &\ge \inf _{\substack{x\in \omega (x_0),\\y\in \partial{\mathcal{B}(\mathcal{A})}}}|x-y|\ge \inf _{\substack{a\in{\mathcal{A}},\\y\in \partial{\mathcal{B}(\mathcal{A})}}}|a-y|= DT({\mathcal{A}}). \end{split} \end{align*}

If all trajectories in $\mathcal{A}$ are dense, the equality is a direct consequence of the fact that $\omega (x_0)\subset{\mathcal{A}}$ is invariant for the flow and thus it is dense in $\mathcal{A}$ .

4.2. Latitude in volume and Basin stability

For one-dimensional dynamical systems, the latitude in width of a basin of attraction (see Definition 4.1) coincides, in fact, with its Lebesgue measure $\mu ({\mathcal{B}(\mathcal{A})})$ (recall that a basin of attraction is always an open set and therefore Lebesgue-measurable; see Proposition 2.4). The notion acquires special interest when the study can be restricted to a measurable set $C\subset{\mathbb{R}}^N$ of finite measure since one can interpret the ‘relative volume’ of the basin of attraction $\mathcal{B}(\mathcal{A})$ (or the portion of $\mathcal{B}(\mathcal{A})$ which lies in $C$ ) in terms of the probability that a perturbed initial condition still belongs to $\mathcal{B}(\mathcal{A})$ .

According to Grümm [ [Reference Gruemm34], Page 4], Holling and collaborators are the first to use this approach. Grümm [ [Reference Gruemm34], Page 7] proposes to measure an attractor’s resilience as the volume of its basin of attraction. Wiley et al. [Reference Wiley, Strogatz and Girvan107] and Fassoni and Yang [Reference Fassoni and Yang29] suggest to measure the volume in a relative sense to some bounded subset, as follows.

Definition 4.9 (Latitude in volume). Let $C\subset{\mathbb{R}}^N$ be Lebesgue-measurable and such that $0\lt \mu (C)\lt \infty$ . We define the latitude in volume of the attractor $\mathcal{A}$ with respect to the region of interest C as:

\begin{equation*}L_v({\mathcal {A}},C) = \frac {\mu ({\mathcal{B}(\mathcal{A})}\cap C)}{\mu (C)}.\end{equation*}

$L_v$ returns a nondimensional value in $[0,1]$ . Particularly, $L_v({\mathcal{A}},C)=0$ if and only if all the trajectories starting in $C$ (except at most those starting in a negligible subset of $C$ ) do not converge to $\mathcal{A}$ as $t\to \infty$ . On the other hand, $L_v({\mathcal{A}},C)=1$ if and only if all the trajectories starting in $C$ (except at most those starting in a negligible subset of $C$ ) converge to $\mathcal{A}$ as $t\to \infty$ .

Remark 4.10 (Region of interest). The ‘region of interest’ $C$ implicitly reveals the characteristics of the models to whom the indicator of latitude in volume should be applied. The system at hand should be at most subjected to isolated impulsive perturbations of bounded magnitude $\delta \gt 0$ , where the term isolated means that the interval of time between two occurrences should be longer than the transient required by the system to reach an ‘indistinguishable proximity’ to the attractor $\mathcal{A}$ from any point in the basin of attraction. Then, the subset $C$ is defined as the set of possible perturbed initial conditions from the attractor [Reference Menck, Heitzig, Marwan and Kurths67], that is $B_\delta ({\mathcal{A}})$ . If additional information on the probability distribution of such perturbations is available, the indicator can be refined (see Definition 4.11). On the other hand, this indicator assumes additional meaning in the case of multi-stable systems [Reference Fassoni and Yang29, Reference Wiley, Strogatz and Girvan107]. In this case the basin’s volume can be interpreted as a measure of the ‘attractor’s relevance’ in the phase space, where it competes against other attractors and their basins over the phase space. Thus, the region of interest should include all the (relevant) attractors ${\mathcal{A}}_i$ , $i=1,\dots,\nu$ (assuming that they are included in a set of bounded measure) as well as the possible perturbed initial conditions of magnitude $\delta$ from them, that is $\bigcup _{i=1}^\nu B_\delta ({\mathcal{A}}_i)\subset C$ . Even so, there is no consensus in the literature as how to choose the region $C$ itself. For example, in [Reference Fassoni and Yang29] the smallest $N$ -dimensional interval $I$ containing $\bigcup _{i=1}^\nu B_\delta ({\mathcal{A}}_i)$ is considered (see the conceptual sketch in Figure 8b). This approach privileges the immediate calculation of the measure of $I$ .

Figure 8. Depiction of two different approaches in choosing the region of interest $C$ depending on the available information. The relevant attractor is sketched as a black dot and identified by the symbol $\mathcal{A}$ , whilst its basin of attraction corresponds to the grey area denoted by the symbols $\mathcal{B}(\mathcal{A})$ . The region of interest $C$ is shown as the area encircled by a black dashed line and denoted by $C$ . (a) The latitude in volume of the attractor $\mathcal{A}$ is calculated with respect to the region of interest $C$ corresponding to the set of expected perturbations of the initial conditions (red region). (b) The latitude in volume of the attractor $\mathcal{A}$ is calculated with respect to the $n-$ dimensional interval $C$ containing also the further attractors ${\mathcal{A}}_1,{\mathcal{A}}_2$ as well as the set of each local attractor’s expected perturbed initial conditions (red regions).

Menck et al. [Reference Menck, Heitzig, Kurths and Schellnhuber66, Reference Menck, Heitzig, Marwan and Kurths67] propose a variation of the latitude in volume (see Definition 4.9) where a weight over the phase space is considered, according to an a priori chosen probability density $\rho$ of the expected perturbations of initial conditions from the attractor. This indicator and its numerical calculation have been investigated in several works including Lundström [Reference Lundström61], Lundström and Adainpää [Reference Lundström and Aidanpää62], Evans and Swartz [Reference Evans and Swartz28], and Mitra et al. [Reference Mitra, Choudhary, Sinha, Kurths and Donner73].

Definition 4.11 (Basin stability). Let $\rho \,:\,{\mathbb{R}}^N\to{\mathbb{R}}^+$ be a Lebesgue integrable function such that $\int _{{\mathbb{R}}^N} \rho (x)\, dx =1$ . The basin stability $S_{\mathcal{B}(\mathcal{A})}$ of an attractor $\mathcal{A}$ with respect to $\rho$ is defined as:

(4.5) \begin{equation} S_{\mathcal{B}(\mathcal{A})}(\rho )= \int _{{\mathbb{R}}^N} \chi _{\mathcal{B}(\mathcal{A})}(x)\rho (x) \, dx, \end{equation}

where $\chi _Z(x)$ is the characteristic function of $Z \subset{\mathbb{R}}^N$ . $S_{\mathcal{B}(\mathcal{A})}$ returns a nondimensional value in $[0,1]$ .

Remark 4.12. Given $C\subset{\mathbb{R}}^N$ Lebesgue-measurable and such that $0\lt \mu (C)\lt \infty$ and considered

\begin{equation*} \rho _C(x) = \begin {cases} 1/\mu (C) & \quad \text {for } x \in C, \\[4pt] 0 & \quad \text {else,}\\ \end {cases} \end{equation*}

one immediately has that $S_{\mathcal{B}(\mathcal{A})}(\rho _C)=L_v({\mathcal{A}},C)$ .

5.1. Return time

Between the ’70s and the ’90s of the last century, the notion of return time as an indicator of resilience (firstly introduced by O’Neill [Reference O’Neill76]) gathered attention and stimulated a fair amount of research (e.g. Cottingham and Carpenter [Reference Cottingham and Carpenter19], DeAngelis et al. [Reference DeAngelis, Bartell and Brenkert22, Reference DeAngelis, Mulholland, Palumbo, Steinman, Huston and Elwood23], Harte [Reference Harte38]). The reason lies in the fact that this definition was able to capture the transient behaviour of the system after an initial perturbation in contrast with the more classical characteristic return time (see Definition 3.1). Although there are several variations in the literature (mostly differing only for the constant of normalisation), they all are similar to the definition below [Reference Neubert and Caswell75, Reference O’Neill76].

Definition 5.1 (Return time). Consider a continuous dynamical system $\phi$ on ${\mathbb{R}}^N$ , and assume that $\mathcal{A}$ is an attractor for $\phi$ and $\mathcal{B}(\mathcal{A})$ its basin of attraction. The return time for the system from $x_p \in{\mathcal{B}(\mathcal{A})}\setminus{\mathcal{A}}$ to $\mathcal{A}$ is defined as:

(5.1) \begin{equation} \overline T_R({\mathcal{A}}, x_p) = \frac{1}{\mathrm{dist}(x_p,{\mathcal{A}})}\int _{0}^{\infty } \mathrm{dist}(\phi (t,x_p),{\mathcal{A}})\, dt, \end{equation}

where $\mathrm{dist}(A,B)$ denotes the Hausdorff semi-distance between the sets $A,B\subset{\mathbb{R}}^N$ . We shall use the symbol $\overline T_R$ instead of the more classical $T_R$ in order to avoid confusion with the characteristic return time presented in Subsection 3.1 . Moreover, for a fixed measurable set $C\subseteq{\mathcal{B}(\mathcal{A})}\setminus{\mathcal{A}}$ , the maximal return time and average return time from $C$ to $\mathcal{A}$ are respectively defined as

\begin{equation*} T_R^{\text{sup}}({\mathcal {A}},C)=\sup _{x_p \in C} \overline T_R({\mathcal {A}}, x_p),\quad \textit {and}\quad T_R^{\text{mean}}({\mathcal {A}},C)=\frac {1}{\mu (C)}\int _{C} \overline T_R({\mathcal {A}},x)\ d\mu (x). \end{equation*}

The definition of return time is a valuable theoretical tool, but the indefinite integral in (5.1) is a clear obstacle in practical applications. Therefore, more calculable indicators, still accounting for the transient dynamics (see Neubert and Caswell [Reference Neubert and Caswell75]), have been developed and were frequently preferred in applications.

The following result shows (as expected) that the integral in (5.1) is finite when the attractor consists of a hyperbolic solution (see Definition 2.6).

Proposition 5.2. Consider a continuous dynamical system on ${\mathbb{R}}^N$ induced by an ordinary differential equation $\dot x =f(x)$ , with $f$ continuously differentiable. Furthermore, assume that the attractor $\mathcal{A}$ is made of the points of a locally attractive hyperbolic solution $\widetilde x$ (see Definition 2.6 ). Then,

\begin{equation*} \overline T_R({\mathcal {A}}, x_p)\lt +\infty,\qquad \textit {for all }x_p \in {\mathcal{B}(\mathcal{A})}. \end{equation*}

Proof. Due to the First Approximation Theorem (see [[Reference Hale36], Theorem III.2.4]), there is $\delta \gt 0$ such that if $\overline x\in{\mathcal{B}(\mathcal{A})}$ satisfies $|\overline x-\widetilde x(s)|\le \delta$ for some $s\in{\mathbb{R}}$ , then for all $t\ge s$ , $|\phi (t,\overline x)-\widetilde x(t)|\lt Ke^{-\beta (t-s)}$ for some $K\gt 1$ and $\beta \gt 0$ . On the other hand, by definition, for every $x_p\in{\mathcal{B}(\mathcal{A})}\setminus{\mathcal{A}}$ there is $s\gt 0$ such that $\mathrm{dist}(\phi (s,x_p),{\mathcal{A}})\lt \delta$ . Then, the result is a consequence of the linearity of the integral in (5.1).

5.2. Resistance of a gradient system

Stability landscapes are a commonly used tool in mathematical ecology as they provide immediate intuitive information on the dynamical features of equilibria (when they exist) of low-dimensional systems. Limitations and warnings of cautious employment of this instrument have been pointed out by many (see [Reference Meyer68] e.g.). In mathematical terms, a stability landscape can be seen as the potential $V$ of a gradient system (and related to the notion of Lyapunov function [Reference Conley17]). This means that the phenomenon can be represented via an ordinary differential equation of the type $\dot x= - \nabla V(x)$ , where $V$ is a twice continuously differentiable real valued function. Then, one can consider the minimal work required to bring the system from an attractor $\mathcal{A}$ to a perturbed initial condition outside the basin of attraction $\mathcal{B}(\mathcal{A})$ as an indicator of resilience (see e.g. Lundström [Reference Lundström61]).

Definition 5.3 (Resistance of a gradient system). Consider a continuous dynamical system induced by an ordinary differential equation $\dot x= - \nabla V(x)$ , where $V\,:\,{\mathbb{R}}^N\to{\mathbb{R}}$ is a twice continuously differentiable. The resistance (potential) of a gradient system with respect to the attractor $\mathcal{A}$ is

\begin{equation*} W({\mathcal {A}})=\inf _{x_0\notin {\mathcal{B}(\mathcal{A})},\,a\in {\mathcal {A}}}\big (V(x_0)-V(a)\big ). \end{equation*}

Remark 5.4. Walker et al. [Reference Walker, Holling, Carpenter and Kinzig103] call resistance of an attractor $\mathcal{A}$ the quantity defined by

\begin{equation*} R_W({\mathcal {A}})=\frac {W({\mathcal {A}})}{L_w({\mathcal {A}})}, \end{equation*}

where $L_w({\mathcal{A}})$ is the latitude in width of $\mathcal{A}$ . As pointed out also by Meyer [Reference Meyer68], care should be taken in the interpretation of this indicator. In particular, if the timescales of disturbance and internal dynamics are considerably different, the slope of the potential function at the state of the system does not necessarily relate to the capability of the system to recover.

5.3. Resilience boundary

Meyer et al. [Reference Meyer, Hoyer-Leitzel and Iams70] suggest to measure the resilience of an attractor by introducing and applying the so-called ‘flow-kick framework’, a sequence of repeated ‘impulsive’ disturbances. Given a continuous dynamical system induced by an ordinary differential equation $\dot x =f(x)$ , $x\in{\mathbb{R}}^N$ , assume that solutions are defined for all $t\in{\mathbb{R}}^+$ . For fixed $\tau \gt 0$ , $\kappa \in{\mathbb{R}}^N$ and $a\in{\mathcal{A}}$ , a flow-kick trajectory $\widetilde x_{\tau,\kappa }(\cdot, a)$ starting in $a\in{\mathbb{R}}^N$ is constructed by considering a partition of the positive real line in intervals $[j\tau,(j+1)\tau )$ , $j\in{\mathbb{N}}$ , and a sequence of functions $\widetilde x_j\,:\,[j\tau,(j+1)\tau ]\to{\mathbb{R}}^N$ such that $\widetilde x_j$ solves the Cauchy problem

\begin{equation*} \dot x=f(x),\qquad \widetilde x_j(j\tau )=\widetilde x_{j-1}(j\tau )+\kappa,\quad \text { and }\quad \widetilde x_0(0)=a\in {\mathcal {A}}. \end{equation*}

One may interpret $\kappa$ as a kick, while $\tau$ represents the flow time, and a pair $(\tau,\kappa )\in (0,\infty )\times{\mathbb{R}}^N$ is called a disturbance pattern. Notice that except for the trivial case where $\kappa$ is the origin in ${\mathbb{R}}^N$ , a flow-kick trajectory is not a solution of the original problem.

Definition 5.5 (Flow-kick resilience set and resilience boundary). Under the previously stated assumptions and notation, a disturbance pattern $(\tau,\kappa )\in (0,\infty )\times{\mathbb{R}}^N$ is said to be in the flow-kick resilience set if

\begin{equation*} \inf _{t\gt 0}P_{\mathcal {A}}\big (\widetilde x_{\tau,\kappa }(t, a)\big )\gt 0\quad \textit {for all }a\in {\mathcal {A}}, \end{equation*}

where $P_{\mathcal{A}}(y)$ is the precariousness of $y\in{\mathbb{R}}^N$ (see Definition 4.7 ). Moreover, the boundary of the flow-kick resilience set in $(0,\infty )\times{\mathbb{R}}^N$ is called the (flow-kick) resilience boundary for the basin of attraction $\mathcal{B}(\mathcal{A})$ .

If $N\gt 1$ , the resilience boundary should be considered as a precautionary threshold. Indeed, it allows to identify the flow-kick trajectories which leave $\mathcal{B}(\mathcal{A})$ at some point, but does not allow to distinguish them from those who eventually return to $\mathcal{B}(\mathcal{A})$ [Reference Meyer, Hoyer-Leitzel and Iams70]. On the other hand, if $N=1$ and $\widetilde x_{\tau,\kappa }(\overline t, a)\notin{\mathcal{B}(\mathcal{A})}$ for some $a\in{\mathcal{A}}$ and disturbance pattern $(\tau,\kappa )$ , the same is true for all $t\ge \overline t$ . In fact, under appropriate assumptions, the scalar case allows to associate a quantitative indicator based on the area of the region in between the flow-kick resilience set and the line $|\kappa |=DT({\mathcal{A}})$ ; see also Figure 11, where $DT({\mathcal{A}})$ denotes the distance to threshold for $\mathcal{A}$ . This indicator has the merit of establishing a relation with distance to threshold (Definition 4.3) and distance to bifurcation (Definition 6.1).

Figure 11. Example of flow-kick resilience and resilience boundary for scalar models in population dynamics. The upper left-hand side panel portrays two populations. Population 1 modelled by $\dot{x}=x(1-\frac{x}{100kt} )(\frac{x}{20kt}-1 ),$ (solid blue line) where $kt$ stands for kilo-tonnes, and population 2 by $\dot{x}=x(1-\frac{x}{100kt})(\frac{x}{20 kt}-1)(0.0002x^2-0.024x+1.4)$ (dashed red line). Both systems have an attracting fixed point at ${\mathcal{A}}=\{100kt\}$ whose basin of attraction is marked in grey. The two systems have same distance to threshold and characteristic return time. Their resilience boundaries are shown in the upper right-hand side panel where also three kick flow patterns are highlighted. In the lower panels, two pairs of flow-kick trajectories are shown for population 1 on the left-hand side and population 2 on the right-hand side. In black the flow-kick trajectory relative to the disturbance pattern $(6\text{ months},-24kt)$ , in red the one relative to $(12\text{ months},-48kt)$ and in green the one relative to $(18\text{ months},-30kt)$ . The disturbance pattern $(6\text{ months},-24kt)$ lies in the resilience set for population 1, but outside the resilience set of population 2. We refer the reader to [Reference Meyer, Hoyer-Leitzel and Iams70] for further details.

5.4. Intensity of attraction

The intensity of attraction aims to measure the strength of the attraction of an attractor. It was first developed by McGehee [Reference McGehee65] in the context of continuous maps acting on a compact metric space, and it has been recently extended by Meyer [Reference Meyer69] and Meyer and McGehee [Reference Meyer and McGehee71] to the continuous-time dynamics induced from an ordinary differential equation on ${\mathbb{R}}^N$ . For the sake of consistency with the setup of this work, we shall present the latter.

Consider a continuous dynamical system induced by an ordinary differential equation $\dot x =f(x)$ , with $f\,:\,U\subset{\mathbb{R}}^N\to{\mathbb{R}}^N$ , $U$ open, and $f$ bounded and globally Lipschitz continuous on $U$ . Moreover, consider an interval $I\subset{\mathbb{R}}$ containing zero and an essentially bounded function $g\in L^\infty ([0,T],{\mathbb{R}}^N)$ , that is, $g\,:\,[0,T] \to{\mathbb{R}}^n$ such that $\inf \{c\geq 0\,:\, \mu (\{t \in [0,T]\,:\,|g(t)|\gt c\})=0\}\lt \infty$ and the perturbed system

(5.2) \begin{equation} \dot{x}=f(x)+g(t). \end{equation}

Equation (5.2) is well-defined (it is in fact a special type of so-called Carathéodory differential equations), and it satisfies the conditions of existence and uniqueness for the associated Cauchy problems [Reference Coddington and Levinson16]. In fact, it is also assumed that for any $x_0\in{\mathbb{R}}^N$ the solution $x_g(\cdot,x_0)$ of (5.2) with $x(0)=x_0$ can be continuously extended to $[0,T]$ . Notice also that the function

\begin{equation*} \phi \,:\,[0,T]\times {\mathbb {R}}^N\times L^\infty \to {\mathbb {R}}^N,\quad (t,x_0,g)\mapsto \phi (t,x_0,g)=x_g(t,x_0) \end{equation*}

vary continuously in $[0,T]\times{\mathbb{R}}^N\times L^\infty$ [Reference Meyer and McGehee71]. As a matter of fact, under the given assumptions, the problem also induces a continuous skew-product flow [Reference Longo56–Reference Longo, Novo and Obaya58]. The fundamental idea behind the concept of intensity is to conceive each bounded set of perturbing functions $g\in L^\infty$ as a class of controls, and to investigate which parts of the phase space can be reached through them. In this way one is able to relate the transient dynamics of a system to the minimal class of controls, which are able to drive the points of an attractor away from it for all future time. In order to state the mathematical definition, let us recall some notation. Given $S \subset{\mathbb{R}}^n$ , we denote by $ \phi (t,S,g)$ the set $\bigcup _{x_0 \in S}\phi (t,x_0,g)$ and call the reachable set from $S$ via controls of norm $r\gt 0$ within the interval of time $I$ , the set

(5.3) \begin{equation} P_r(S,[0,T])=\bigcup _{t \in [0,T]} \bigcup _{g \in F_r([0,T])} \phi (t,S,g), \end{equation}

where $F_r([0,T]) =\{g\in L^\infty ([0,T],{\mathbb{R}}^n)\,:\, \left \lVert{g}\right \rVert _\infty \lt r\}$ .

Definition 5.6 (Intensity of attraction). The intensity of attraction $\mathcal{I}({\mathcal{A}},[0,T])$ of an attractor $\mathcal{A}$ is defined as:

(5.4) \begin{equation} \mathcal{I}({\mathcal{A}},[0,T]) = \sup \{r \geq 0 \mid P_r({\mathcal{A}},[0,T])\subset K \subset{\mathcal{B}(\mathcal{A})} \textit{ for some compact } K\subset{\mathbb{R}}^N \}. \end{equation}

Interestingly, the intensity of attraction provides information on the persistence of an attractor not only against a time-dependent perturbation. Meyer et al. prove that for any vector field, which differs in norm from $f$ less than its intensity of attraction, a compact positively invariant subset of ${\mathbb{R}}^N$ exists which contains the attractors of both systems (see [[Reference Meyer and McGehee71], Theorem 6.12]).

5.5. Expected escape times

It is a well-known fact that the presence of a stochastic perturbation can have a non-negligible impact on the dynamics of a system and possibly trigger a critical transition [Reference Kuehn51]. Dennis et al. [Reference Dennis, Assas, Elaydi, Kwessi and Livadiotis24] who, in turn, refer to Gardiner [Reference Gardiner31] for the main idea suggest to judge the resilience of an attractor for a system subjected to white noise via the notion of mean first passage time, also known as expected escape time.

Definition 5.7 (Expected escape times). Consider a scalar stochastic differential equation

(5.5) \begin{equation} dX =f(X)\,dt+\sqrt{\nu (X)}\,dW(t), \end{equation}

where $\nu \,:\,{\mathbb{R}}\to{\mathbb{R}}^+$ is continuous, and $ W$ is a standard Brownian motion. If (5.5) has a stationary distribution with probability density $p(x)$ , then the area under $p(x)$ between $x_0$ and $x_1$ gives the long-term proportion of time that the process $X$ will spend in the interval $(x_0,x_1)$ . The expected escape times correspond to the mean first passage times for a process starting at $x_0$ to increase (resp. decrease) to $x$ ,

\begin{align*} \begin{split} \tau _1(x)&=2\int _{x_0}^x\frac{\int _0^yp(z)\,dz}{\nu (x)p(y)}\,dy\quad \text{for }x_0\lt x, \\[4pt] \tau _2(x)&=2\int ^{x_0}_x\frac{\int _y^\infty p(z)\,dz}{\nu (x)p(y)}\,dy\quad \text{for }x_0\gt x. \end{split} \end{align*}

The expected escape times represent the mean amount of time necessary for the process $X$ to attain the value $x$ starting from $x_0\lt x$ and $x_0\gt x$ , respectively. As explained in [Reference Dennis, Assas, Elaydi, Kwessi and Livadiotis24], one or both these integrals can be infinite. For example if $x=0$ acts like an absorbing value and the system lies close to it, there is a positive probability that the first escape time takes an infinite value and so $p(x)$ is not integrable near $0$ . An analogous reasoning holds for the second integral if $\infty$ acts as an absorbing boundary. In this case, the right tail of $p$ is not integrable.

Invariance with respect to change of coordinates

The expected escape times are presented for scalar stochastic equations. A linear change of coordinates in this context would mean a scaling of the only available dependent variable. For example, consider $Y=aX$ , with $a\gt 0$ . Then, one obtains the stochastic differential equation $dY =af(Y/a)\,dt+\sqrt{a^2\nu (Y/a)}\,dW(t)$ . Note that if $p(z)$ is the stationary distribution for $X$ , then $\widetilde p(z)=p(z/a)/a$ is the stationary distribution for $Y$ . Then, the expected escape times for the process $Y$ starting at $y_0$ are

\begin{align*} \begin{split} \widetilde \tau _1(y)&=2\int _{y_0}^y\frac{\int _0^\mu \widetilde p(z)\,dz}{a^2\nu (y/a)\widetilde p(\mu )}\,d\mu \quad \text{for }y_0\lt y, \\[5pt] \widetilde \tau _2(y)&=2\int ^{y_0}_y\frac{\int _\mu ^\infty \widetilde p(z)\,dz}{a^2\nu (y/a)\widetilde p(\mu )}\,d\mu \quad \text{for }y_0\gt y. \end{split} \end{align*}

For example, let us consider the first integral; then; we have that for $y_0\lt y$ and using $\widetilde p(z)=p(z/a)/a$ , and suitable changes of variables,

\begin{align*} \begin{split} \widetilde \tau _1(y)&=2\int _{y_0}^y\frac{\frac{1}{a}\int _0^\mu p(z/a)\,dz}{a\nu (y/a) p(\mu/a)}\,d\mu = 2\int _{y_0/a}^{y/a}\frac{\int _0^{a\eta } p(z/a)\,dz}{a\nu (y/a) p(\eta )}\,d\eta =2\int _{y_0/a}^{y/a}\frac{\int _0^{\eta } p(\zeta )\,d\zeta }{\nu (y/a) p(\eta )}\,d\eta \\ \end{split} \end{align*}

which is equal to $\tau _1(y/a)$ for the process $X$ starting at $y_0/a$ . An analogous result holds for $\tau _2$ . General diffeomorphic reparametrisations of $\mathbb{R}$ do not in general preserve the value of the indicators.

6.1. Distance to bifurcation

Bifurcation theory is a classical branch of nonlinear dynamical systems which studies the lack of topological equivalence in phase space upon (smooth) change of one or more of the system’s parameters. The classic theory for autonomous dynamical systems encompasses the bifurcations of local steady states, periodic orbits, as well as homoclinic and heteroclinic orbits [Reference Kuznetsov55]. For its characteristics, a bifurcation point entails a profound change in the dynamics of the system, which may include the change of stability or the disappearance of a considered attractor. In practical terms, a bifurcation point can sometimes be undesirable, especially when the system exhibits hysteresis or irreversibility [Reference Carpenter, Ludwig and Brock13]. It is therefore natural that the distance to a bifurcation point has been suggested as a measure of resilience [Reference Dobson25, Reference Ives and Carpenter43].

Definition 6.1 (Distance to bifurcation). Let $\mathcal{A}$ be an attractor for the dynamical system induced by (6.1) $_{\lambda _0}$ and denote by $\Lambda _{\text{bif}}({\mathcal{A}})\subset{\mathbb{R}}^M$ the set of bifurcation points for $\mathcal{A}$ . The distance to bifurcation for $\mathcal{A}$ is defined as

\begin{equation*} D_{\rm {bif}}({\mathcal {A}})=\inf _{\lambda ^* \in \Lambda _{\text{bif}}} |\lambda _0-\lambda ^*|. \end{equation*}

Notably, distance to bifurcation can be interpreted as a distance to threshold (see Definition 4.3) in parameter space.

6.2. Resistance, elasticity and persistence

Harrison [Reference Harrison37] examines the transient behaviour of a system which undergoes an instantaneous change of parameters and returns to the original values, just as instantly, after an interval of time. The system’s response during and after the perturbation are suggested as indicators of the system’s capacity to withstand a change and to recover. Harrison names these features respectively resistance and resilience (although, to avoid confusion, we will use Orians’s nomenclature [Reference Orians77] and call Harrison’s resilience elasticity). Using Webster’s words [Reference Webster, Waide and Patten105], the first describes ‘the ability to resist displacement’, whereas the latter ‘the rate of recovery to the original condition’. With the notation set at the beginning of the section, we have the following definitions:

Definition 6.2 (Resistance). Given an attractor $\mathcal{A}$ for (6.1) $_{\lambda _0}$ , and fixed $T\gt 0$ such that for all $\lambda \in \Lambda$ , and $a\in{\mathcal{A}}$ , $\phi _\lambda (\cdot,a)$ is defined at least in $[0,T]$ , the resistance of system (6.1) $_{\lambda _0}$ at $\mathcal{A}$ with respect to $\Lambda$ and the stress period $[0,T]$ is defined as:

\begin{equation*}R({\mathcal {A}},\Lambda, T)=\sup _{\substack {t\in [0,T],\, a\in {\mathcal {A}}_{\lambda _0}\\ \lambda \in \Lambda }}|\phi _{\lambda }(t,a)-\phi _{\lambda _0}(t,a)|.\end{equation*}

Definition 6.3 (Elasticity). Given an attractor $\mathcal{A}$ for (6.1) $_{\lambda _0}$ , and fixed $T\gt 0$ such that $\phi _{\lambda }(T,{\mathcal{A}})\subset \mathcal{B}({\mathcal{A}})$ for all $\lambda \in \Lambda$ , consider the continuous function $\Phi _{\lambda _0}\,:\,[T,\infty )\times{\mathcal{A}}\times \Lambda \to{\mathbb{R}}^N$ defined by

\begin{equation*} \Phi _{\lambda _0}(t,a, \lambda ):=|\phi _{\lambda _0}\big (t-T,\phi _{\lambda }(T,a)\big )-\phi _{\lambda _0}(t,a)|. \end{equation*}

The elasticity of the system (6.1) $_{\lambda _0}$ at $\mathcal{A}$ with respect to $\Lambda$ and the stress period $[0,T]$ is defined as:

\begin{equation*}E({\mathcal {A}},\Lambda, T)=\sup _{\substack {t\in [T,\infty ),\, a\in {\mathcal {A}}_{\lambda _0}\\ \lambda \in \Lambda }}\left (\frac {1}{\Phi _{\lambda _0}(t,a, \lambda )}\frac {d\, \Phi _{\lambda _0}(t,a, \lambda )}{dt}\right ). \end{equation*}

In Figure 12, it is possible to appreciate the complementary information provided by resistance and elasticity in a simple scalar differential model from population dynamics.

Figure 12. Effect of a stress period of Harrison’s type on two populations following the logistic model $\dot{x}=rx(1-x/K)$ with carrying capacity $K=1$ , and growth rates, respectively, equal to $r_1=1$ and to $r_2=2$ . A stress period $[0,1]$ is applied to both models during which the carrying capacity is diminished of the $20\%$ (i.e. $\Lambda =\{0.8\})$ . The resistance of the attractor $x^*=1$ is lower for the first species than the second. Conversely, the elasticity of the attractor $x^*=1$ is higher for the first species than the second.

Remark 6.4. The definitions of resistance and elasticity presented above admit a natural generalisation. Instead of using the original trajectory inside the attractor as a reference, in practice it may be sufficient to record the displacement from the attractor or the rate of convergence to the attractor. This would lead to the following weaker definitions

\begin{equation*} R^w({\mathcal {A}},\Lambda, T)=\sup _{\substack {t\in [0,T],\, a,a'\in {\mathcal {A}}_{\lambda _0}\\ \lambda \in \Lambda }}|\phi _{\lambda }(t,a)-a'|, \end{equation*}

and

\begin{equation*} \Phi ^w_{\lambda _0}(t,a, \lambda ):=\inf _{a'\in {\mathcal {A}}}|\phi _{\lambda _0}\big (t-T,\phi _{\lambda }(T,a)\big )-a'|, \end{equation*}

\begin{equation*}E^w({\mathcal {A}},\Lambda, T)=\sup _{\substack {t\in [T,\infty ),\, a\in {\mathcal {A}}_{\lambda _0}\\ \lambda \in \Lambda }}\left (\frac {1}{\Phi ^w_{\lambda _0}(t,a, \lambda )}\frac {d\, \Phi ^w_{\lambda _0}(t,a, \lambda )}{dt}\right ). \end{equation*}

Ratajczak et al. [Reference Ratajczak, D’Odorico, Collins, Bestelmeyer, Isbell and Nippert83] use a setup similar to Harrison’s to address the problem of the interaction between the duration and the ‘intensity’ of a stress period in determining a tipping of the system from the original basin of attraction to a different one (if it applies). In accordance with Harrison’s nomenclature, we shall call this property persistence.

Definition 6.5 (Persistence). The persistence of an attractor $\mathcal{A}$ for (6.1) $_{\lambda _0}$ with respect to a stress intensity $\lambda \in \Lambda$ is given by

\begin{equation*} P_\lambda ({\mathcal {A}})=\sup \{T\gt 0\mid \phi _\lambda (t,{\mathcal {A}})\subset \mathcal {B}({\mathcal {A}})\textit { for all }0\le t\le T\}. \end{equation*}

The persistence of an attractor $\mathcal{A}$ for (6.1) $_{\lambda _0}$ with respect to a stress duration $T\gt 0$ is given by

\begin{equation*} P_T({\mathcal {A}})=\sup \{\rho \gt 0\mid \phi _\lambda (t,{\mathcal {A}})\subset \mathcal {B}({\mathcal {A}})\textit { for all }0\le t\le T\textit { and } |\lambda -\lambda _0|\le \rho \}. \end{equation*}

Remark 6.6. It is clear that Harrison’s formulations of resistance and elasticity as well as of persistence generalise to any perturbation on a finite time interval rather than just a piece-wise constant change of parameters. Consider (6.1) $_{\lambda _0}$ for some $\lambda _0\in{\mathbb{R}}^M$ , $T\gt 0$ , and a set $\Lambda$ of functions from ${\mathbb{R}}\times{\mathbb{R}}^N$ into ${\mathbb{R}}^M$ such that if $\lambda \in \Lambda$ , then for every compact set $K\subset{\mathbb{R}}^N$ , $\sup _{x\in K}\lambda (t,x)$ is locally integrable and $\lambda (t,x)=\lambda _0$ for all $x\in{\mathbb{R}}^N$ and $t\notin [0,T]$ . Hence, one can consider the problem

\begin{equation*} \dot x =f\big (x,\lambda (t,x)\big ), \end{equation*}

and still define resistance, elasticity and persistence as above.

6.3. Resilience against rate-induced tipping

Sizable changes in the output of a system upon a negligible variation of the input are referred to as critical transitions or tipping points. Motivated by current challenges in nature and society [Reference Gladwell32, Reference Scheffer88], the study of the several mechanisms leading to a critical transition has experienced significant mathematical interest as well [Reference Kuehn51, Reference Kuehn52]. In recent years for example, the phenomenon of rate-induced tipping [Reference Wieczorek, Ashwin, Luke and Cox106] has been proposed as an alternative mechanism for a critical transition, with respect to the more classical autonomous bifurcations and noise-induced tipping [Reference Ashwin, Perryman and Wieczorek5]. Rate-induced tipping can be seen as a special type of nonautonomous bifurcation for the differential equation

(6.2) \begin{equation} \dot x =f\big (x,\gamma (rt)\big ),\qquad x\in{\mathbb{R}}^N,t\in{\mathbb{R}},r\gt 0 \end{equation}

obtained from (6.1) through a time-dependent variation of parameters $\gamma \,:\,{\mathbb{R}}\to{\mathbb{R}}^M$ which is assumed to be bounded, continuous and asymptotically constant: $\lim _{t\to -\infty }\gamma (t)=\lambda _0,\ \lim _{t\to \infty }\gamma (t)=\lambda _\infty$ . In particular, the parameter $r\gt 0$ represents the rate at which the time-dependent variation takes place. We shall maintain a concise presentation and leave the detailed assumptions and technical aspects for Appendix A.

Typically, the study of rate-induced tipping is carried out under the assumption that the time-dependent variation of the parameters does not cross any point of autonomous bifurcation; that is we assume that for all $\lambda$ in the image of $\gamma$ , the dynamical systems induced by (6.1) $_\lambda$ , are all topologically equivalent. One might expect that for any fixed attractor ${\mathcal{A}}_0$ of the past limit-problem, and denoted by ${\mathcal{A}}_f$ the ‘corresponding’ attractor of the future limit-problem, for all $r\gt 0$ there is an attractor of (6.2) $_r$ which approaches ${\mathcal{A}}_0$ as $t\to -\infty$ and ${\mathcal{A}}_f$ as $t\to \infty$ . However, this is not always the case and, depending on $r$ , a local attractor of (6.2) $_r$ limiting at ${\mathcal{A}}_0$ as $t\to -\infty$ may fail to track ${\mathcal{A}}_f$ as $t\to \infty$ . This phenomenon is called rate-induced tipping, and the specific value $r^*_{\gamma }({\mathcal{A}}_0)$ of the parameter where it occurs is called critical rate (see Definition A.1). It seems therefore natural that the critical rate $r^*_{\gamma }({\mathcal{A}}_0)$ determines a threshold of resilience for (6.1) $_{\lambda _0}$ against the time-dependent change of parameters $\gamma$ , in the sense that for all $r\lt r_\gamma ^{*}$ , where $r_\gamma ^{*}$ is the inferior over all the attractors of (6.1) $_{\lambda _0}$ , the tracking of attractors from the past to the future is always verified.

This statement assumes further significance if we look at rate-induced tipping as of a nonautonomous bifurcation. Then, the bifurcation point $r_\gamma ^{*}$ can be used to construct a distance to bifurcation analogous to the one in Definition 6.1.

Definition 6.7 (Distance to rate-induced tipping). Let ${\mathcal{A}}_0$ be an attractor of the past limit-problem, and ${\mathcal{A}}_f$ the ‘corresponding’ attractor of the future limit-problem. The distance to rate-induced tipping for (6.2) and for ${\mathcal{A}}_0$ is defined as

\begin{equation*} D_{\rm {rate}}({\mathcal {A}}_0)=|r-r^*_\gamma ({\mathcal {A}}_0)|. \end{equation*}

It should be noted, that in contrast to autonomous bifurcations, estimating the value of a critical rate can be hard if not impossible [Reference Kuehn and Longo53]. There is, furthermore, a substantial difference from Definition 6.1. The tipping point $r_\gamma ^{*}$ depends on the function $\gamma$ that realises the variation of parameters. Therefore, it is beneficial to investigate the continuous dependence of $r^*_\gamma$ with respect to the variation of $\gamma$ . We hereby propose a first result in this direction, when all the attractors of (6.1) $_\lambda$ are hyperbolic equilibria. The result guarantees the lower semi-continuity (and if it applies the continuity) of $r^*_\gamma$ with respect to $\gamma$ in a specific set of functions: fixed $t_0\in{\mathbb{R}}$ , we consider the set $\Omega (\gamma, t_0)$ of twice continuously differentiable functions $\omega \,:\,{\mathbb{R}}\to \Lambda$ such that $\omega (t)=\gamma (t)$ for all $t\le t_0$ and $\lim _{t\to \infty }\omega (t)=\lambda \in \Lambda$ . For a detailed explanation of the symbols appearing in the statement and for the proof, please consult Appendix A.

The model

We shall study the parametric differential problem

(7.1) \begin{equation} \dot{x}=f(x,r,L)=rx\left (1-\frac{x}{K}\right )\left (\frac{x}{L}-1\right ), \end{equation}

where $x\in{\mathbb{R}}^+$ represents the population size, $r\gt 0$ the intrinsic growth rate, $K\gt 0$ the carrying capacity of the environment, and $0\lt L\le K$ an Allee threshold. For the sake of convenience, we will set the carrying capacity to $K=1$ . For $0\lt L\lt 1$ , the problem has three equilibria, at $x_0=0$ , $x_1=1,$ and $x_L=L$ . The first two are stable and the latter unstable. Therefore, $L$ determines a threshold below which the population cannot persist.

7.1. Resilience of the attractor $x_1=1$

In this subsection, we apply the presented indicators to the attractor $x_1=1$ of the considered model $\dot x=f(x,r,L)$ . Note that when $L=1$ , the system undergoes a transcritical bifurcation and the equilibria $x_L$ and $x_1$ collide. When the calculations cannot be performed explicitly, the software Matlab is used for numerical computation and simulation. In particular, we use the function ode45 with the options on the relative and absolute tolerance respectively set to RelTol = 1e-12 and AbsTol = 1e-12.

Nonlinear transient dynamics

The return time (Definition 5.1) does not admit an explicit closed form. Therefore, we employ a routine to estimate it numerically. In Figure 13 we show an approximation of the mean return time from the set $C=(L,1)$ estimated via a Monte Carlo simulation (uniform sampling on the set $(L+10^{-7},1)$ with $10000$ points) depending on the parameters $r\in [0.01,0.5]$ and $L\in [0.05,0.95]$ . For all points, the return time is approximated via a finite time integration ending when the considered trajectory achieves $1-10^{-10}$ . Note that the considered problem can also be interpreted as a gradient system. Thus, we are able to derive its resistance $W$ (Definition 5.3) as

\begin{equation*} W(x_1)=\int _L^{1} f(x,r,L)\, dx=-\frac {(L-1)^3(L+1)r}{12L}. \end{equation*}

The intensity of attraction $\mathcal{I}$ (Definition 5.6) is obtained by calculating the maximum of $f$ on the interval $[L,1].$ We show the behaviour of $\mathcal{I}(x_1)$ depending on $r\in [0.01,0.5]$ and $L\in [0.05,0.95]$ in Figure 14.

Figure 13. Comparison of the indicators for the attractor $x_1=1$ of the model (7.1) upon the variation of the parameters $r\in [0.01,0.5]$ and $L\in [0.05,0.95]$ . From top to bottom: opposite of the real part of the dominant eigenvalue $EV=\frac{1}{T_R}$ , distance to threshold $DT$ and the reciprocal of the mean return time $T_R^{\text{mean}}(1,(L+10^{-7},1))$ . On the left, we see the scaled indicator values for the given parametric subspace. In the middle a heat map for the parametric subspace, where yellow indicates the highest and dark blue the lowest values of the indicators (see the logarithmic colour scale on the right). Black lines correspond to the contour lines, which mark the parameter combinations with the same value of the indicators. If the environmental changes drive the parameters along the contour curves, all indicators, apart from $DT$ , are not able to capture the approaching bifurcation independently. On the right, we see the dependence of the indicator values on the parameter $L.$ The three pictures on each row showcase the same surface viewed by different viewpoints.

Figure 14. Continuation of the comparison of the indicators for different parameter choices of the model (7.1): resistance (potential) $W,$ intensity of attraction $\mathcal{I},$ resistance $R$ and elasticity $E$ (please note that in order to have higher resilience corresponding to a higher numerical value, we worked with the reciprocal of $R$ and $E$ ). We consider perturbation of the environmental capacity to $K_p=0.9,$ for time $t=[0,10].$ On the left, we see the scaled indicator values for a given parametric subspace: $[0.01,0.5]= r \times L=[0.05,0.95]$ for resistance (potential) and intensity and $[0.01,0.5]= r \times L=[0.05,0.89]$ for resistance and elasticity. In the middle, there is a heat map for the parametric subspace, where yellow indicates the highest and blue the lowest values of the indicators (see the logarithmic colour scale on the right). Black lines correspond to the contour lines, which mark the parameter combinations with the same values of the indicator. On the right, we see the dependence of the indicator values on the parameter $L.$ The three pictures on each row showcase the same surface viewed by different viewpoints. We see that resistance and elasticity behave in an opposite manner, elasticity being higher for lower values of $\{r,L\}$ while resistance being lower.

7.1.1. Parametric study

We perform a parametric analysis of (7.1), with the aim of portraying the behaviour of some of the considered indicators as the transcritical bifurcation point $L=1$ is approached (results are summarised in Figures 13 and 14).

We focus on the parameter subspace given by $r\in [0.01,0.5]$ , $L\in [0.05,0.95]$ and calculate seven indicators (see previous Subsection 7.1 for details): reciprocal of characteristic return time $EV=\frac{1}{T_R}$ , distance to threshold $DT$ , mean return time $T_R^{\text{mean}}(1,(L+10^{-7},1)),$ resistance $W,$ intensity of attraction $\mathcal{I}$ , resistance $R$ and elasticity $E$ . The last two indicators are calculated for the perturbation of the environmental capacity $K_p=0.9$ for a stress period $T=10.$

For a better comparison of the indicators, the estimates shown in Figures 13 and 14 are scaled to lie in $[0,1]$ by $(x-x_{\min })/(x_{\max }-x_{\min }),$ where $x$ represents the indicator, and minimum $x_{\min }$ and maximum $x_{\max }$ corresponds to the local minimum and maximum of the indicator on the considered subset of parameters. For the same reason, instead of $T_R^{\text{mean}},R,E$ , we represent their reciprocal so that the value $1$ always corresponds to the highest estimate of resilience. Additionally, resistance and elasticity we further restrict the parameter space so to guarantee that these indicators are well-defined.

It is possible to appreciate that upon fixing a value $r\gt 0$ and increasing $L$ , all the considered indicators capture a loss of resilience for the system as the bifurcation point approaches. The only exception is given by the resistance which increases since the rate of exponential decay towards the attractor decreases in the nearby of the bifurcation point. This phenomenon is also directly related to the so-called critical slowing-down effect [Reference Dakos, Carpenter, van Nes and Scheffer21, Reference Scheffer, Bascompte and Brock89, Reference Strogatz96]. On the other hand, there are also parameter combinations, indicated by the black contour lines, where the indicator stays constant.

Figure 15. Five population strategies with different choices of parameters in equation (7.1). The table on the left-hand side contains the parameters for each species. The plot on the right-hand side showcases the graphs of $f(x,r,L)$ for $x\in [0,1]$ depending on the chosen values of parameters $r$ and $L$ .

Figure 16. Each row represents one indicator and each column one species. The chromatic scale applies row-wise as the values of the indicators have not been normalised. Darker tones of blue correspond to higher values of resilience. With the exception of species 4, all the other species are considered as the most resilient for at least one of the represented indicators.

7.1.2. Comparison of five species

Finally, we aim to showcase the behaviour of the indicators with respect to five different species growth strategies for the Allee effect model (7.1). The different values of the parameters as well as the respective graphs of $f(x,r,L)$ for $x\in [0,1]$ are shown in Figure 15. Given the previous considerations, we shall focus on the following indicators of resilience: the opposite of the dominant eigenvalue $EV=1/T_R,$ distance to threshold $DT,$ mean return time $T_R^{\text{mean}}(1,(L+10^{-7},1)),$ resistance (potential) of a gradient system $W,$ intensity of attraction $\mathcal{I},$ resistance $R$ and elasticity $E.$

The numerical values for these indicators are shown in Figure 16. In all cases, $x_1=1$ is the reference attractor. Each row corresponds to a single indicator, each column to a species. The chromatic scale applies row-wise as the values of the indicators have not been normalised. Darker tones of blue correspond to higher values of resilience. It is apparent how answering the question ‘Which species is the most resilient?’ is not trivial. The ranking of the resilience of each species varies depending on the considered indicator.

In particular, it is interesting to notice that, not only there is no overall agreement on the resilience ranking, but some indicators ( $EV$ and $DT$ e.g.) provide antipodal responses. It is possible to single out three main clusters: a) $EV$ and $1/T_R^{\text{mean}}$ , b) $1/E$ , $I$ and $W$ , c) $1/R$ and $DT$ . With the exception of species ‘4’, all the other species are considered as the most resilient for at least one of the represented indicators.

Finally, in Figure 17 we can see the resilience boundaries for the considered five species. As studied in [Reference Meyer, Hoyer-Leitzel and Iams70], it is possible to extrapolate a measure of resilience for the attractor by calculating the area between the resilience boundary and the line $\kappa = DT(x_1)=1-L$ (see Subsection 5.3). However, since the distances to threshold for the five species are different it seems natural to ask if the obtained values are representative and comparable. The right panel of Figure 17 shows the numerical values obtained for the indicator before and after dividing them by the respective distances to threshold.

Figure 17. On the left-hand side, the resilience boundaries of five species (see Figure 15) modelled through (7.1). On the right-hand side, first row, the numerical values of the areas between the resilience boundary and the respective distance to threshold calculated for each species. The second row contains the same values now divided by the respective distances to threshold $DT(x_1)=1-L$ . The colour gradient should be intended row-wise. Darker shades of blue correspond to a higher resilience.

7.2. Resilience of the attractor $x_0=0$

In this subsection, we perform a (shorter) analysis of resilience for the attractor $x_0=0$ . In order to solve the singularity arising from posing $L=0$ in (7.1), we shall study the problem

(7.2) \begin{equation} \dot{x}=g(x,r,L)=rx\left (1-x\right )\left (x-L\right ), \quad x\ge 0, \end{equation}

which is equivalent to (7.1) when $L\gt 0$ (and $K=1$ ), via the variation of time scale $t\mapsto Lt$ . When $L$ tends to zero, the continuous variation of the solutions can be deduced using Tikhonov’s theorem [Reference Tikhonov98]. Note that when $L=0$ , (7.2) undergoes a transcritical bifurcation and the equilibria $x_L$ and $x_0$ collide.

Nonlinear transient dynamics

Same as for the attractor $x_1=1$ , we use a numerical routine to approximate the return time (Definition 5.1). The results are now shown in Figure 18. We have used the same setting and stopping criteria for the Monte Carlo simulation as before (uniform sampling on the set $(0, L-10^{-7})$ with $10000$ points) depending on the parameters $r\in [0.01,0.5]$ and $L\in [0.05,0.95]$ . The resistance $W$ (Definition 5.3) is obtained as

\begin{equation*} W(x_0)=\int _L^{0} g(x,r,L)\, dx=\frac {L^3(2-L)r}{12}. \end{equation*}

The intensity of attraction $\mathcal{I}$ (Definition 5.6) is obtained by calculating the maximum of $g$ on the interval $[0,L].$ We show the behaviour of $\mathcal{I}(x_0)$ depending on $r\in [0.01,0.5]$ and $L\in [0.05,0.95]$ in Figure 18.

Figure 18. Comparison of the indicators for the attractor $x_0=0$ of the model (7.2) upon the variation of the parameters $r\in [0.01,0.5]$ and $L\in [0.05,0.95]$ . From top to bottom: $EV=\frac{1}{T_R},$ $DT,$ $T_R^{\text{mean}}(0,(0, L-10^{-7})),$ $W(0)$ and $I(0)$ . On the left, we see the scaled indicator values for the given parametric subspace. In the middle a heat map for the parametric subspace, where yellow indicates the highest and dark blue the lowest values of the indicators (see the logarithmic colour scale on the right). Black lines correspond to the contour lines, which mark the parameter combinations with the same value of the indicators. If the environmental changes drive the parameters along the contour curves, all indicators, apart from $DT$ , are not able to capture the approaching bifurcation independently. On the right, we see the dependence of the indicator values on the parameter $L.$ The three pictures on each row showcase the same surface viewed by different viewpoints.

7.3. Synthesis, comparisons and remarks

In this section, we carried out the analysis of a scalar population model with Allee effect (7.1) in terms of the resilience of its attractors at $x_1=1$ and $x_0=0$ . The different outcomes given by some of the indicators of resilience presented in this paper have been tested against the variation of the growth rate and of the Allee threshold (both treated as parameters). For the attractor $x_1$ , the indicators were also compared across five different species with fixed pairs of parameters (see part 7.1.1 and 7.1.2). The resulting absence of a non-ambiguous answer to the question ‘which species is more resilient?’ may seem contradictory at first, but it is a known feature of resilience. The measured resilience depends on both the specific dynamic feature under analysis and the class of acting perturbations. As pointed out by Carpenter et al. [Reference Carpenter, Walker, Anderies and Abel14] and Tamberg et al. [Reference Tamberg, Heitzig and Donges97], one needs to carefully define the context and questions before coming to conclusions about resilience.

For instance, if we are interested in the survival of a species, that is, the resilience of the equilibrium $x_1$ , the following conclusions may be draft. If we are solely interested in the largest one-time shock that the species can withstand (e.g. a sudden ecological catastrophe), but we are not concerned with the dynamics of the recovery, the most suitable indicator is distance to threshold $DT$ (or other basin shape indicators – for their comparison see Proposition 4.13). As explained in subsection 7.1 and graphically shown in Figure 13, $DT$ is exclusively dependent on the value of the Allee threshold $L$ and is not changed by the magnitude of the growth rate $r$ . From a biological point of view, this might be interpreted as the number of individuals in the species that have to suddenly disappear for the species to go extinct, independently of its growth rate. If the Allee threshold is crossed, the growth rate will not help. The lower $L$ is, the weaker the Allee effect, which means the species population can sustain bigger shocks from which they can eventually recover. Thus, we observe the highest resilience for values of $L$ close to $0$ and the lowest for high Allee thresholds around $1$ (see Figure 13).

Conversely, when the focus is on the speed of recovery after a single ‘small’ perturbation of the number of individuals in the species, a local indicator such as $EV$ may better fulfil the task. Indeed, $EV$ corresponds to the exponential rate of convergence towards a hyperbolic attracting solution. In our model, this indicator depends on both the intrinsic growth rate $r$ and threshold $L$ (see subsection 7.1 and Figure 13). A higher growth rate $r$ results in a faster recovery, whereas the higher the Allee threshold $L$ , the slower the recovery. As extensively explained in Section 3, all the local indicators are reliable only when ‘small’ perturbations are taken into account. Whenever a larger perturbation is considered, the local indicators such as $EV$ might not be informative enough and we have to turn to other transient dynamics indicators. One straightforward choice is the global version of the local indicator $EV\,:\,$ the mean return time $T_R^{\text{mean}}.$ In our model, $T_R^{\text{mean}}$ of the attractor $x_1=1$ shows overall similar behaviour as the local indicator $EV$ . This is mostly due to the extreme simplicity of the model at hand. In general, however, these indicators may provide substantially different results – examples of this can be found in part 7.1.2.

Next, we focus on continuous and repeated perturbations. In biological terms, we can consider a perturbation caused by a permanent outflow of species (imagine e.g. fishing, hunting, or illness). This can be modelled as $\dot{x}=f(x)+g(t),$ where $f(x)$ is the unperturbed dynamics and $g(t)$ is the perturbation term (e.g. fishing). Then, the indicator intensity $\mathcal{I}$ will find the critical magnitude of the outflow term $g(t)$ after which the species can still recover. The overall trend seen in Figure 14 is similar to that of $EV$ – the higher the reproductive speed in terms of growth rate $r$ , the better can the species compensate for the outflow, although there are still differences with other indicators (see part 7.1.2).

Another type of perturbation, different from those already mentioned, is recurring stress on the species numbers. This is where periods of no disturbance are followed by strong and rapid stress on the species population. Think, for instance, of a seasonal illness or periodic fishing. This can be modelled in terms of periodic, discrete shocks of some strength. In this case, the resilience boundary is the most obvious choice of indicator. It captures both the magnitude of shocks that a species can sustain to not fall below the critical population numbers given by the Allee threshold, as well as, the recovery capabilities of the species during the periods without perturbation. High values of growth rate $r$ help the recovery rate, and low values of Allee threshold $L$ permit severe shocks to the species number to still lie in the recoverable numbers (see part 7.1.2 for comparison).

A different interpretation for the transient dynamics is given by the indicator of resistance (potential) which quantify how much work needs to be put into eradicating the species at hand (see Figure 14). In our context, this is possible due to the specific structure of the equation that can be written as a gradient problem. It is not surprising that this indicator provides estimates in accordance with the other indicators of resilience for transient dynamics – the higher the growth rate, the harder it is to eliminate the species. Likewise, the weaker the Allee effect, the harder it is to eradicate species below the threshold $L.$

On the other hand, the changing environment is another type of stress species factor. In the previous subsection, we discussed the effect of lowering the species’ carrying capacity. This could be caused, for example, by deteriorating environmental conditions. In terms of resilience, we might ask the following questions: once the environmental stress factor has ended, what is the impact on the species’ population numbers? If the species can recover (which is not always possible in general) how fast the recovery is? The former can be answered by measuring the resistance (Harrison) $R$ , the latter by elasticity $E$ (see Figure 14). Resistance (Harrison) quantifies how the population numbers decrease after the stress ends. What we can notice is that populations with higher growth rates are more sensitive to the decline of carrying capacity. The same is true for low values of the Allee threshold. These species decline faster due to the way the model is defined where the overall tendency of the vector field to have a higher magnitude also leads to a faster decrease in the population numbers. Elasticity, on the other hand, determines the speed of recovery after the stress ends by measuring the species population’s ability to flourish again. Notice that for each species it finds the slowest recovery speed, thus, the worst-case scenario. The starting position of recovery is described by resistance; therefore, they are related. We notice that species with higher growth rates $r$ , are better equipped to recover, by quickly reproducing their numbers, also for low Allee thresholds. Apart from resistance, this trend is consistent with other indicators mentioned. We can see that the species can have different strategies to fight temporal environmental stresses: typically, resistance and elasticity are opposite aspects of resilience. Finally, we like to stress that similar comments can be made for the equilibrium $x_0=0$ , which represents the extinction state. The interest of the extinctions state lies in the opportunity to control infesting species or diseases. Through the parametric study carried out in Subsection 7.2, and the numerical evidence showcased in Figure 18, we illustrated, once again, the different outcomes of resilience indicators. Particular emphasis is posed on the mean return time which displays a new phenomenon: the long persistence of an infesting species (which is eventually destined to extinction) when the Alee threshold approaches the carrying capacity. In practical terms, this case seems quite rare in reality where the Alee threshold is usually considered in $(0,0.5)$ .