Principles for guiding future research on resilience and tipping points

In the phrase 'resilience of what to what' (Carpenter et al 2001), the first 'what' refers to the some set of outcomes related to the system state under study, such as the state of a lake or a city, while the second 'what' pertains to the external disturbance affecting that system. A system can be conceptualized as any entity that exhibits interrelated components and interactions within a defined boundary. Three types of systems can be distinguished based on the nature of mass and energy exchanges with the environment: (1) an isolated system, which exchanges neither mass nor energy with its surroundings; (2) a closed system, which exchanges energy but not mass with its surroundings; and (3) an open system, which exchanges both mass and energy with its surroundings. Most natural systems are open systems, providing the necessary conditions to maintain complexity and high order in a low-entropy state (Nicolis and Prigogine 1977).

To conduct resilience research, it is essential to identify a state variable (or set of state variables), typically represented on the y-axis of the S-shaped curve in figure 1. The choice of this variable depends on the specific change one aims to address. Let $\overline {x}$ be a state variable, which can represent critical system aspects, such as vegetation density in a grazing system (Noy-Meir 1975) or average earth system temperature (Nicolis and Nicolis 1981). For simplicity, we will focus on a single state variable, though complex systems may involve multiple state variables, as discussed by Boers et al (2022) and Anderies et al (2023). Various ecosystem state variables are reviewed and listed in Yi and Jackson (2021).

The state variable is subject to the following dynamical equation,

$$\begin{align}\frac{{{\text{d}}\overline {x}}}{{{\text{d}}t}} = f\left( {\overline {x}, \alpha } \right),\end{align} \tag{ 1 }$$

where $\alpha$ represents an environmental control parameter, and $f$ is in general a nonlinear function. This equation can be built, for example, based on fundamental principles such as (1) Mass Conservation The change in mass is equal to the incoming mass minus the outgoing mass from the system; or (2) Energy Conservation The change in energy (calculated as temperature multiplied by specific heat) is equal to the incoming energy to the system minus the outgoing energy from the system. For any physical system being modeled, the nonlinearity of $f\left( {\overline {x}{\text{ }}\alpha } \right)$ primarily arises from the complex interplay of positive and negative feedback interactions.

The steady state refers to a condition where the value of the state variable remains constant over time. This occurs when the system's dynamics are balanced, meaning that the processes driving change are exactly counteracted by those opposing it. Mathematically, it is defined as,

$$\begin{align}\frac{{\text{d}{{\overline x}_s}}}{{\text{d}t}} = 0 = f\left( {{{\overline x}_s},{\text{ }}\alpha } \right).\end{align} \tag{ 2 }$$

Due to the system's nonlinearity, solutions of steady state ${\overline {x}_s}$ can be multiple and their stability properties depend on the different values of the parameter $\alpha$.

The solutions determined by equation (2) are also called fixed points or singular points in dynamical system theory. Depending on their stability properties, these fixed points can be further classified as nodes, foci, saddle points, and limit cycles, etc (Strogatz 1994).

In resilience science and ecology, the term 'equilibrium' is often used instead of 'steady state.' The concept of equilibrium originates from thermodynamics and has a precise physical meaning. It characterizes a state of maximum entropy, where there are no gradients or irreversible fluxes (Yi and Jackson 2021).

However, in living systems, maintaining a high level of order under nonequilibrium conditions requires a continuous flow of negative entropy through energy and mass exchanges. As Prigogine (1980) demonstrated, life systems cannot maintain themselves near equilibrium, as the entropy production of such nonequilibria is minimal. Instead, Prigogine clarified that nonequilibrium can be a source of order in a life system that is far from equilibrium. While equilibrium represents one of many possible steady states for a system—where no irreversible processes occur, and entropy is maximized—the steady state of a life system is inherently far from equilibrium.

Resilience refers to a system's ability to handle disturbances by resisting change and recovering from it (Helfgott 2015). Although related, a disturbance is distinct from a perturbation (Rykiel 1985, Yi and Jackson 2021). A disturbance acts as the trigger, while a perturbation is the resulting change in the system's state variables caused by the disturbance. Disturbances are characterized by their type, frequency, and intensity, whereas perturbations reflect the extent of change in an ecosystem. Different ecosystems may experience varying levels of perturbation from the same disturbance. When faced with the same disturbance, an ecosystem exhibiting a smaller perturbation demonstrates greater resistance, whereas a larger perturbation indicates weaker resistance.

The stability of steady states ${\overline {x}_s}{\text{ }}$ can be examined by a linear stability analysis approach (Nicolis and Prigogine 1977). The perturbed steady state can be written as,

$$\begin{align}\overline {x} = {\overline {x}_s} + x^{\prime},\end{align} \tag{ 3 }$$

where $x^{\prime}$ is a deviation from that steady state, called a perturbation or fluctuation. Fluctuations are inherent variations within a system, while perturbations are fluctuations specifically caused by external disturbances. Both involve deviations from the steady state, with their magnitudes depending on the stability of the steady state. Near a stable steady state, both tend to decay, while near a tipping point, they can grow exponentially (figure 1).

Assuming that the perturbation $x^{\prime}$ is sufficiently small, we can linearize equation (1) using a first-order Taylor expansion. Applying the steady state condition (2) we derive the perturbation equation,

$$\begin{align}\frac{{{\text{d}}x^{\prime}}}{{{\text{d}}t}} = \frac{{{\text{d}}f}}{{{\text{d}}\overline {x}}}{\Bigg|_{{{\overline x}_s}}}x^{\prime} = \lambda x^{\prime}.\end{align} \tag{ 4 }$$

The solution to this perturbation equation is

$$\begin{align}x^{\prime} \propto {{\text{e}}^{\lambda t}}.\end{align} \tag{ 5 }$$

Thus ${\overline {x}_s}{\text{ is stable}}$ when $\lambda < 0$ and unstable when $\lambda \unicode{x2A7E} 0$. The resilience of the system can be quantified by the recovery rate ${R_\tau }$ (Neubert and Caswell 1997, Arnoldi et al 2018),

$$\begin{align}{R_\tau } = - \frac{1}{{x^{\prime}}}\frac{{{\text{d}}x^{\prime}}}{{{\text{d}}t}} = - \frac{{{\text{dln}}x^{\prime}}}{{{\text{d}}t}} = - \lambda .\end{align} \tag{ 6 }$$

For a stable steady state ($\lambda < 0$), the perturbation decays, and the system returns to the original state. These stable steady states are known as attractors (solid lines in the S-shaped curve in figure 1). The absolute value of $\left| \lambda \right|$ is termed the recovery rate: a larger (more negative) $\left| \lambda \right|$ indicates a faster recovery rate, meaning the system's ability to absorb perturbations—and thus its resilience—is stronger.

Conversely, $\lambda > 0$, the perturbation increases exponentially, causing the system to diverge from the original state. These unstable steady states are noted by dotted line in the S-shaped curve in figure 1. As $\lambda$ becomes less negative, the system's resilience weakens, and the recovery speed slows. At the tipping point where $\lambda = 0$, the system loses its resilience entirely, resulting in an infinite recovery time (i.e. no recovery).

Dakos and Kéfi (2022) describe this form of resilience as engineering resilience due to small perturbations around stable steady-states. A system with greater stability (indicated by a more negative λ) returns more quickly to stable steady-states.

Linear stability analysis is not suitable for handling large perturbations. To address this, the concept of resilience potential $U\left( x \right)$ has been widely used to assess global stability, often illustrated through the ball-and-cup–diagram (figure 1). The potential function can be visualized as a landscape where valleys represent stable states and hills represent unstable states. The deeper the valley, the greater the system's global resilience to disturbances around the associated stable state. A shallow valley indicates that the system can easily shift to a different state, reflecting lower global resilience.

The resilience potential is defined as the integration of the nonlinear function $f\left( {\overline {x}\alpha } \right)$,

$$\begin{align}U\left( {\overline {x}\alpha } \right) = - \mathop {\mathop \int \nolimits }\limits_{\overline x}^{{{\overline x}_s}} f\left( {\xi ,\alpha } \right){\text{d}}\xi .\end{align} \tag{ 7 }$$

Thus, the governing equation (1) can be written as,

$$\begin{align}\frac{{{\text{d}}\overline {x}}}{{{\text{d}}t}} = f\left( {\overline {x} \alpha } \right) = - \frac{{{\text{d}}U}}{{{\text{d}}\overline {x}}}.\end{align} \tag{ 8 }$$

The steady state equation becomes

$$\begin{align}f\left( {{{\overline x}_s},{\text{ }}\alpha } \right) = - \frac{{\text{d}U}}{{\text{d}\overline {x}}}{\Bigg|_{{{\overline x}_s}}} = 0.\end{align} \tag{ 9 }$$

The extrema of $U$ correspond to the steady states ${\overline {x}_s}$. We can consider the linear response near a steady state. From equation (4),

$$\begin{align}\lambda = \frac{{{\text{d}}f}}{{{\text{d}}\overline {x}}} = - \frac{{{{\text{d}}^2}U}}{{{\text{d}}{{\overline x}^2}}}.\end{align} \tag{ 10 }$$

For a bistable system, two stable steady states (${\overline {x}_{s1}},{\text{ }}{\overline {x}_{s3}}$) are separated by an unstable one (${\overline {x}_{s2}}$). The engineering resilience of each stable steady state is defined by the curvature of the resilience potential,

$$\begin{align}{\lambda _{s1}} = - \frac{{{{\text{d}}^2}U}}{{{\text{d}}{{\overline x}^2}}}{\Bigg|_{{{\overline x}_{s1}}}},\,{\lambda _{s3}} = - \frac{{{{\text{d}}^2}U}}{{{\text{d}}{{\overline x}^2}}}{\Bigg|_{{{\overline x}_{s3}}}}.\end{align} \tag{ 11 }$$

Ecological resilience is measured by potential depth with respect to the unstable state (Dakos and Kéfi 2022)

$$\begin{align}\Delta {U_{s1,3}} = U\left( {{{\overline x}_{s2}}} \right) - U\left( {{{\overline x}_{s1,3}}} \right),\end{align} \tag{ 12 }$$

which represents the capacity of an ecosystem to absorb disturbances (figure 1).

It is important to note that the resilience potential approach described above is most applicable in systems that can be approximated as gradient systems. In gradient systems, the dynamics can be described as the downhill motion of a ball in a potential landscape, and this allows for a meaningful definition of the resilience potential. In higher-dimensional systems, the assumption of a gradient structure may not hold, and the potential depth as defined in equation (12) may not adequately capture the dynamics of the system. Therefore, while the ball-and-cup analogy is useful for illustrating basic concepts, it is an approximation that may not capture the full complexity of higher-dimensional or non-gradient systems. In such cases, more sophisticated approaches may be needed to characterize resilience.

A bifurcation is a fundamental concept in dynamical systems where a small change in system parameters can cause a qualitative change in the system's behavior. A bifurcation process can be described by phase space language. Phase space is a multidimensional space in which all possible states of a system are represented, with each dimension corresponding to one of the system's variables. In this space, each point represents a unique state of the system at a given time. The solution $\overline {x}\left( t \right)$ of equation (1) with initial condition ${\overline {x}_0}$ traces out a path known as a trajectory in phase space. Fixed points (steady states) are constant over time. These fixed points can be: (1) stable fixed points (attractors): where trajectories tend to converge, indicating the system's ability to return to steady state after a disturbance, and (2) unstable fixed points (repellers): where trajectories move away, indicating that small perturbations can lead to significant changes in the system's state.

When the control parameter $\alpha$ is fixed, the number of fixed points in the phase space and their stability are determined, which in turn determines the flow structure (topological structure) of the phase space trajectories. If a small disturbance occurs to a parameter, the flow structure of the phase space trajectory, along with the number of fixed points and their stabilities, generally remains unchanged.

However, when the control parameter $\alpha$ reaches a critical value ${\alpha _c}$, even a small change in $\alpha$ can fundamentally alter the flow structure of the phase space trajectory, including the number and stability of the fixed points (Glendinning 1994). This phenomenon is known as a bifurcation, with ${\alpha _c}$ referred to as the bifurcation point. The bifurcation point represents the theoretical threshold in system dynamics where changes in $\alpha$ lead to significant modifications in the system's behavior and the stability of fixed points (Glendinning 1994). In the context of resilience theory, this bifurcation point is often called a tipping point. In practical terms, the location of the tipping point may differ from the bifurcation point, as it depends on the rate at which external forcing is applied. A faster rate of change in external forcing can cause the system to cross the tipping point earlier due to delayed responses, also known as rate-induced tipping (Ashwin et al 2012).

At the bifurcation point, several scenarios can occur: a stable steady state may become unstable, with new stable branches (steady states) bifurcating from the original state; two steady states, typically one stable and one unstable, can exchange stability; or a pair of stable and unstable fixed points may collide and annihilate each other or emerge together at the bifurcation point. Additionally, a stable or unstable oscillating solution can emerge as a steady state becomes unstable, as seen in Hopf bifurcations (both supercritical and subcritical).

There are many different types of bifurcation, depending on nonlinearity of $f\left( {{{\overline x}_s},{\text{ }}\alpha } \right)$, as well as the number of steady state solutions and the nature of their stability (Strogatz 1994). Four typical bifurcations are: saddle-node bifurcation, where two fixed points (one stable and one unstable) collide and annihilate each other, leading to a loss of stability; pitchfork bifurcation, where a fixed point transitions from being stable to unstable, resulting in the emergence of two stable branches; transcritical bifurcation, where a stable and an unstable fixed point collide, exchange stability, and continue as two separate branches; and Hopf bifurcation, where a stable fixed point loses stability, giving rise to a limit cycle that leads to periodic oscillations. Different tipping types correspond to these different bifurcation types.

In a deterministic model, such as equation (1), the system settles into a stable state without further fluctuations. However, real ecosystems are subject to random perturbations, much like how Brownian particles experience fluctuating molecular forces, denoted as $F\left( t \right){\text{ }}$ (Risken 1996). To account for this stochasticity, we extend the deterministic model into a stochastic differential equation,

$$\begin{align}\frac{{{\text{d}}x}}{{{\text{d}}t}} = f\left( {x,\alpha } \right) + F\left( t \right) = - \frac{{{\text{d}}U}}{{{\text{d}}x}} + F\left( t \right),\end{align} \tag{ 13 }$$

where $x$ is now a random variable, and $F\left( t \right)$ represents random forces assumed to be white noise, characterized by,

$$\begin{equation}\begin{gathered} \left\langle {F\left( t \right)} \right\rangle = 0, \hfill \\ \left\langle {F\left( t \right)F\left( {t^{\prime}} \right)} \right\rangle = {\sigma ^2}\delta \left( {t - t^{\prime}} \right), \hfill \\ \end{gathered} \end{equation} \tag{ 14 }$$

where $< >$ denotes an expectation value and ${\sigma ^2}$ is the variance of the fluctuating force $F\left( t \right)$ (Nicolis and Nicolis 1981). Equation (13) is known as the ecological Langevin equation (see appendix B for details).

While the deterministic model describes how the system evolves toward a stable state, the Fokker–Planck equation provides a complementary framework by describing how the probability density function $P\left( {x,t} \right)$ evolves over time. This equation quantifies the likelihood of transitions between coexisting stable states, offering insights into ecosystem stability and resilience under stochastic conditions (appendix B). The evolution of the probability density function is governed by:

$$\begin{align}\frac{{\partial P\left( {x,t} \right)}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {\frac{{{\text{d}}U\left( {x,\alpha } \right)}}{{{\text{d}}x}}P\left( {x,t} \right)} \right) + \frac{{{\sigma ^2}}}{2}\frac{{{\partial ^2}P\left( {x,t} \right)}}{{\partial {x^2}}}.\end{align} \tag{ 15 }$$

For a steady state, $\frac{{\partial {P_s}}}{{\partial t}} = 0,$ we obtain,

$$\begin{align}\frac{{{\text{d}}U\left( {x,\alpha } \right)}}{{{\text{d}}x}}{P_s}\left( x \right) + \frac{{{\sigma ^2}}}{2}\frac{{\partial {P_s}\left( x \right)}}{{\partial x}} = {\text{constant}}.\end{align} \tag{ 16 }$$

In systems with physical boundaries, it is often assumed that ${P_s}\left( x \right) = 0$ at the boundaries of the state variable x. This assumption allows us to set the constant on the right-hand side of equation (16) to zero. However, for systems where the state space is not bounded, the steady-state probability should be defined across all possible states. Therefore, to accommodate such cases, the constant may not necessarily be zero. Nevertheless, for simplicity, we proceed with the assumption that ${P_s}\left( x \right)$ decays sufficiently fast at large $\left| x \right|$, making the constant effectively zero. Thus, an analytical solution of the steady-state probability density is obtained,

$$\begin{align}{P_s}\left( x \right) = {Z^{ - 1}}{{\text{e}}^{ - \frac{2}{{{\sigma ^2}}}U\left( x \right)}},\end{align} \tag{ 17 }$$

where $Z$ is a normalization factor, determined by,

$$\begin{align}\mathop {\mathop \int \nolimits }\limits_\infty ^{ - \infty } {P_s}\left( x \right){\text{d}}x = 1.\end{align} \tag{ 18 }$$

Thus,

$$\begin{align}Z = \mathop {\mathop \int \nolimits }\limits_\infty ^{ - \infty } {{\text{e}}^{ - \frac{{2U}}{{{\sigma ^2}}}}}{\text{d}}x.\end{align} \tag{ 19 }$$

The resilience potential $U\left( x \right)$ plays a crucial role in expressing the stability of a system, particularly in the presence of stochastic perturbations. The steady-state probability density ${P_s}$ is directly related to the resilience potential and is an important measure of the likelihood of the system occupying a particular state. Specifically, ${P_s}$ reaches its maximum when the resilience potential $U\left( x \right)$ is at a minimum, corresponding to a stable steady state $\left( {{{\overline x}_{s1}},{{\overline x}_{s3}}} \right)$, and ${P_s}$ reaches its minimum when $U\left( x \right)$ is at a maximum, corresponding to an unstable steady state $\left( {{{\overline x}_{s2}}} \right)$.

This relationship is mathematically expressed as:

$$\begin{align}\frac{{{{\text{d}}^2}U}}{{{\text{d}}{x^2}}}{\Bigg|_{{{\overline x}_s}}}\left\{ {\begin{array}{*{20}{l}} {\kern-3pt< 0,{\text{ maximum }}U,{\text{ minimum }}{P_s},{{\overline x}_s}{\text{ is unstable}}} \\ {\kern-3pt> 0,{\text{minimum }}U,{\text{maximum }}{P_s},{{\overline x}_s}{\text{ is stable }}} .\end{array}} \right.\end{align} \tag{ 20 }$$

The stability of these steady states is governed by the curvature of the resilience potential $U\left( x \right)$. In particular, the second derivative $\frac{{{{\text{d}}^2}U}}{{{\text{d}}{x^2}}} = U^{{\prime} {\prime}}$ determines whether a steady state is stable or unstable. When $U^{{\prime} {\prime}} > 0$, the system is in a stable steady state, where small perturbations decay, and the system tends to return to the stable state. Conversely, when $U^{{\prime} {\prime}} < 0$, the system is in an unstable steady state, where small perturbations grow, leading the system away from the original state.

A critical condition occurs when $U^{{\prime} {\prime}} = 0$. This marks the bifurcation point, where the system undergoes a qualitative change in behavior. At this point, the system's stability changes, leading to a tipping point where a small perturbation can cause a significant shift in the system's state.

For a bistable system, the steady-state probability density ${P_s}\left( x \right)$ can be calculated under three conditions (figure 2): (i) $U\left( {{{\overline x}_{s1}}} \right) < U\left( {{{\overline x}_{s3}}} \right)$, (ii)${\text{ }}U\left( {{{\overline x}_{s1}}} \right) > U\left( {{{\overline x}_{s3}}} \right)$, (iii) $U\left( {{{\overline x}_{s1}}} \right) \simeq U\left( {{{\overline x}_{s3}}} \right)$, which is related to the normalization factor $Z$. Although obtaining an exact solution of equation (19) is challenging for arbitrary values of $\sigma$, for small $\sigma$, the steepest descent method (appendix A) can provide an approximate solution.

Zoom In Zoom Out Reset image size

The steady-state probability density can be expressed as:

$$\begin{align}{P_s}\left( x \right) & = \left\{ \sigma \sqrt {\frac{\pi }{{U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right)}}} {{\text{e}}^{ - \frac{2}{{{\sigma ^2}}}U\left( {{{\overline x}_{s1}}} \right)}} \right.\nonumber\\ &\quad \left.+ \sigma \sqrt {\frac{\pi }{{U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s3}}} \right)}}} {{\text{e}}^{ - \frac{2}{{{\sigma ^2}}}U\left( {{{\overline x}_{s3}}} \right)}} \right\}^{ - 1}{{\text{e}}^{ - \frac{2}{{{\sigma ^2}}}U\left( x \right)}},\end{align} \tag{ 21 }$$

which can be simplified under the three conditions:

• (i)  
  If $U\left( {{{\overline x}_{s1}}} \right) < U\left( {{{\overline x}_{s3}}} \right)$, due to the inverse small factor $\frac{1}{{{\sigma ^2}}}$ in the exponent, the difference between two minima $U\left( {{{\overline x}_{s1}}} \right){\text{ and }}U\left( {{{\overline x}_{s3}}} \right)$ will be amplified enormously. The term containing $U\left( {{{\overline x}_{s3}}} \right)$ will become vanishingly small in (21) (Nicolis and Nicolis 1981). By expanding $U\left( x \right)$ into a Taylor series around the stable steady state ${\overline {x}_{s1}}$, and under the condition $2\Delta U/{\sigma ^2} \gg 1$, ${P_s}\left( x \right)$ around ${\overline {x}_{s1}}$ simplifies to,

  $$\begin{align}{P_s}\left( x \right) \approx \frac{1}{\sigma }\sqrt {\frac{{U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right)}}{\pi }} {{\text{e}}^{ - \frac{1}{{{\sigma ^2}}}U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right){{\left( {x - {{\overline x}_{s1}}} \right)}^2}}}.\end{align} \tag{ 22 }$$

  This is a Gaussian distribution centered on ${\overline {x}_{s1}},$ dominates and attracts all initial conditions in the limit of longtime (light blue in figure 2).
• (ii)  
  If $U\left( {{{\overline x}_{s1}}} \right) > U\left( {{{\overline x}_{s3}}} \right)$, similarly, ${P_s}\left( x \right)$ around ${\overline {x}_{s3}}$ simplifies to,

  $$\begin{align}{P_s}\left( x \right) \approx \frac{1}{\sigma }\sqrt {\frac{{U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s3}}} \right)}}{\pi }} {{\text{e}}^{ - \frac{1}{{{\sigma ^2}}}U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s3}}} \right){{\left( {x - {{\overline x}_{s3}}} \right)}^2}}}.\end{align} \tag{ 23 }$$

  This is another Gaussian distribution centered on ${\overline {x}_{s3}}$.
• (iii)  
  If $U\left( {{{\overline x}_{s1}}} \right) \simeq U\left( {{{\overline x}_{s3}}} \right)$, both terms in (21) must be retained with $U\left( {{{\overline x}_{s1}}} \right) = U\left( {{{\overline x}_{s3}}} \right) = {U_m}$. This results in,

  $$\begin{align}{P_s}\left( x \right) = \frac{1}{\sigma }{\left\{ {\sqrt {\frac{\pi }{{U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right)}}} + \sqrt {\frac{\pi }{{U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right)}}} } \right\}^{ - 1}}{{\text{e}}^{ - \frac{2}{{{\sigma ^2}}}\left\{ {U\left( x \right) - {U_m}} \right\}}}.\end{align} \tag{ 24 }$$

Resulting in a two-hump distribution with equal-height maxima (red dashed lines in figure 2).

The stability of a system at a steady state ${\overline {x}_{s1}}$ is closely tied to the curvature of its resilience potential, $U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right)$, at that state. A steeper potential well, indicated by a larger $U^{{\prime} {\prime}}\left( {{{\overline x}_{s1}}} \right)$, signifies strong local stability, as small perturbations will result in only minor deviations, and the system will quickly return to its stable state. However, this increased local stability may correspond to a narrower basin of attraction, meaning the system is less tolerant to larger perturbations. If the system is pushed outside this narrow basin, it may transition to a different state, indicating lower global resilience against significant disturbances.

In formulas (22) and (23), the term $x - {\overline {x}_{s1}}$ and $x - {\overline {x}_{s3}}$ represent the distance from the steady states ${\overline {x}_{s1}}$ and ${\overline {x}_{s3}}$, respectively. This distance defines the basin of attraction for each steady state. The basin of attraction refers to the horizontal extent of the potential well around a stable steady state, encompassing all initial states from which the system will return to that steady state after a perturbation. The width of these basins, determined by the curvature $U^{{\prime} {\prime}}\left( {{{\overline x}_{s1}}} \right)$ or $U^{{\prime} {\prime}}\left( {{{\overline x}_{s3}}} \right)$, plays a crucial role in the system's stability and resilience. A steeper curvature (larger $U^{{\prime} {\prime}}$) means a narrower basin of attraction, reflecting high local stability but lower global resilience. Conversely, a shallower curvature (smaller $U^{{\prime} {\prime}}$) indicates a broader basin, meaning the system can absorb larger disturbances, demonstrating higher global resilience.

The interplay between stability and resilience is evident in the shape and size of the basin of attraction (Livina et al 2011). A steep potential well suggests strong stability but a potentially smaller basin, while a broader basin reflects higher resilience, allowing the system to tolerate larger disturbances. The figure 2 illustrates these relationships by showing how a wider basin of attraction correlates with higher resilience, giving the system state (represented by the ball) more room to move without losing stability.

The fundamental characteristics of the steady-state solution (21) of the Fokker–Planck equation depend on the resilience potential. Specifically, the depth of the bistable potential well plays a crucial role in determining the likelihood of the system being in a particular state. Variations in the depth of the potential well are linked to the model parameters (see figure 2). As these parameters change, the ecosystem system may transition from one potential well to another, overcoming potential barriers in the process. The duration required for this transition is referred to as the tipping time (Dakos and Kéfi 2022), denoted by τ.

To estimate the tipping time or understand the transition processes between alternative stable states, one must solve the time-dependent solution of the Fokker–Planck equation. Various methods have been developed for these solutions since Kramers' pioneering work in 1940, with different researchers contributing, such as Caroli et al (1979) and Tomita et al (1976). For instance, Tomita et al (1976) initially transformed the Fokker–Planck equation into a Schrödinger form and subsequently employed traditional quantum mechanics methods to solve for the eigenvalues of the Schrödinger equation. The tipping time is estimated below (details provided in Tomita et al (1976)),

$$\begin{align}\tau &= \frac{1}{{{\lambda _1}}} = \frac{{2\pi }}{{\sqrt {\left| {U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s2}}} \right)} \right|} }}{{\text{e}}^{\frac{{2U\left( {{{\overline x}_{s2}}} \right)}}{{{\sigma ^2}}}}}\left\{ \sqrt {U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right)} {{\text{e}}^{\frac{{2U\left( {{{\overline x}_{s1}}} \right)}}{{{\sigma ^2}}}}}\right.\nonumber\\ &\quad \left. + \sqrt {U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s3}}} \right)} {{\text{e}}^{\frac{{2U\left( {{{\overline x}_{s3}}} \right)}}{{{\sigma ^2}}}}} \right\}^{ - 1}.\end{align} \tag{ 25 }$$

The tipping time depends on the initial position of the system with a specific potential well:

If the system initially resides in the left well of the potential and $U\left( {{{\overline x}_{s1}}} \right) > U\left( {{{\overline x}_{s3}}} \right)$, then (25) reduces to

$$\begin{equation}\tau = \frac{{2\pi }}{{\sqrt {\left| {U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s2}}} \right)} \right|U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s1}}} \right)} }}{{\text{e}}^{\frac{{2\left( {U\left( {{{\overline x}_{s2}}} \right) - U\left( {{{\overline x}_{s1}}} \right)} \right)}}{{{\sigma ^2}}}}}.\end{equation} \tag{ 26 }$$

Conversely, if the system starts in the right well of the potential, and $U\left( {{{\overline x}_{s1}}} \right) < U\left( {{{\overline x}_{s3}}} \right)$, equation (25) reduces to

$$\begin{equation}\tau = \frac{{2\pi }}{{\sqrt {\left| {U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s2}}} \right)} \right|U{^{{\prime} {\prime}}}\left( {{{\overline x}_{s3}}} \right)} }}{{\text{e}}^{\frac{{2\left( {U\left( {{{\overline x}_{s2}}} \right) - U\left( {{{\overline x}_{s3}}} \right)} \right)}}{{{\sigma ^2}}}}}.\end{equation} \tag{ 27 }$$

The tipping time depends exponentially on the ratio of the barrier height $\Delta U$ to the variance ${\sigma ^2}$/2, making it highly sensitive to changes in the relative magnitude of these two quantities. When the ratio $2\Delta U/{\sigma ^2} \gg 1,$ the tipping time is amplified to an exceptionally long duration. This can be attributed to a combination of two possible factors reasons: either the variance ${\sigma ^2}$ is small and the random disturbance is insufficient to drive the system to a new stable state, or the potential well is deep.

The term 'resilience' was initially introduced into ecological literature by Holling (1973) as a 'measure of the persistence of systems and their ability to absorb change and disturbance while still maintaining the same relationships between populations or state variables.'

Holling's definition adopts a broader ecological perspective, recognizing resilience as the system's capacity to recover from disturbances while preserving its fundamental structures and functions. Current classifications distinguish between engineering resilience and ecological resilience, each with distinct characteristics and applications.

In this context, ecological resilience refers to the ability of a system to recover from a perturbed state to one of its stable steady states, without altering the self-organized processes and structures that define it. By contrast, engineering resilience is typically quantified through the recovery rate (λ, a dominant eigenvalue) or recovery time (the reciprocal of the recovery rate), determined using traditional linear stability analysis. This approach characterizes the decay of a perturbed state near a stable steady state, which is the essence of engineering resilience (Dakos and Kéfi 2022). However, engineering resilience applies only to systems with a single stable steady state and small perturbations.

In reality, perturbations might be large enough to drive ecosystems to alternative stable steady states, where resilience is measured by the resilience potential depth (as in equation (12)), a concept termed ecological resilience (Holling 1973, Peterson et al 1998, Gunderson 2000, Dakos and Kéfi 2022). The resilience potential is valuable not only for understanding global stability—allowing the system to recover from all possible perturbations—but also for systems capable of transitioning from one stable steady state to another, including through tipping points (at the threshold value of the control parameter or bifurcation point).

However, ecological resilience diverges from the original meaning, which emphasized that a perturbed system must return to its stable steady state, maintaining unchanged self-organized processes and structures. Each of the multiple stable steady states may have its unique self-organized structures (Peterson et al 1998). From a deterministic perspective, when one stable steady state becomes unstable, the stable and unstable steady states coincide, forcing the system to tip toward an alternative stable state (see figure 3 in Yi and Jackson (2021)). Resilience, therefore, is not defined for unstable systems (Pimm 1984), although some researchers refer to the long-term rate of return near the tipping point as asymptotic resilience (Arnoldi et al 2018).

From a stochastic perspective, the system can transition from a more stable state with a deeper potential well to a less stable state with a shallower potential well in response to extreme events such as droughts, flooding, and wildfires. However, the more stable steady state, having the highest probability density, attracts all initial conditions in the long run (figure 2). The potential function incorporates the results of linear stability analysis (where the recovery rate corresponds to the curvature of the potential), extending beyond resilience. While resilience lacks a clear definition at the tipping point, the potential function provides a more evident and comprehensive representation.

In practice, ecosystem resilience is often estimated using a resilience index based on observational data. However, there is no universal resilience index. Typically, a resilience index is defined as a measure of the damage and recovery of system state variables relative to unperturbed states (Yi and Jackson 2021). The steady states before and after a perturbation can be approximated by a multi-year average.

Identifying the ideal state variable for a resilience index is challenging due to its synthetic nature. Nonetheless, various observational data, such as the tree-ring width index (RWI), normalized difference vegetation index (NDVI), enhanced vegetation index, and leaf area index (LAI), can serve as proxies for system state variables. The choice of a specific proxy depends on the investigator's aims and may vary from case to case (see appendix C "Ecosystem Resilience Index).

In this context, Xu et al (2022) employed the NDVI ratio as a resistance index to investigate the sensitivity of forest canopy height to drought, revealing that tall trees are more vulnerable under such conditions. Similarly, both Rocha (2022) and Zampieri et al (2021) utilized gross primary productivity (GPP) as a proxy to construct resilience indices. The advantage of this method lies in its ability to leverage Earth System Model-generated databases for a global-scale examination of GPP resilience under various emissions scenarios.

The tipping point is defined by the critical value of a control parameter at which the system undergoes an abrupt shift from a self-organized structure (loses stability) to a new one. A slight change in environmental conditions near the tipping point can lead to a catastrophic transformation in the system's structure and function. Therefore, identifying where and when a tipping point is likely to occur is crucial for developing adaptation strategies and preparedness policies (Valdes 2011).

The characteristics exhibited by a system as it approaches a tipping point, known as tipping behavior or tipping dynamics, have generated significant attention in recent research (Bastiaansen et al 2022, Yi et al 2024). Tipping behavior includes phenomena such as critical slowing down (Prettyman et al 2022), early warning signals (Scheffer et al 2009), late warning signals (aka delayed transition or ghost) (Vidiella et al 2018), Turing bifurcation (Rietkerk et al 2021), and tipping cascades (Rocha et al 2018, Klose et al 2021).

Critical slowing down refers to the diminishing recovery rate or resilience as the system nears the tipping point, with recovery time increasing. During this process, the system's short-term memory increases, leading to rising variance and autocorrelation (Held and Kleinen 2004, Livina and Lenton 2007, Dakos et al 2015, Prettyman et al 2018). Consequently, time series methods based on power spectrum analysis have been widely used to develop indicators of critical slowing down (Bastiaansen et al 2022, Prettyman et al 2022, Clarke et al 2023). However, predicting future tipping points using these indicators remains extremely challenging. Dakos et al (2015) noted that these indicators sometimes fail to accurately signal an impending transition.

In contrast, the tipping points associated with past abrupt climate transitions between glaciations and interglacials are more clearly evidenced in paleoclimate data. For instance, Ramadhin et al (2021) found that increased temperature variance preceded glacial termination. The rising variance near tipping points is considered a useful indicator for threshold tipping in various ecological and social systems, as this method does not require specific model parameters or detailed dynamical information (Carpenter and Brock 2006).

The so-called climate tipping elements (subsystems of the earth system that may exhibit tipping behavior), such as the Greenland Ice Sheet, the Atlantic Meridional Overturning Circulation, permafrost, monsoon systems, and the Amazon rainforest, may become unstable at their individual tipping points, leading to shifts in their self-organized structures and functions. The feedback interactions among these tipping elements, driven by rising temperatures, mean that the tipping behavior of one element can influence others. This interplay may result in simultaneous or cascading tipping events (Klose et al 2021, Franzke et al 2022). Understanding cascading tipping dynamics is currently at a conceptual stage (Rocha et al 2018, Klose et al 2021), and empirical evidence of interactions among tipping elements remains an open question (Scheffer and Van Nes 2018). Recently, Liu et al (2023) employed a climate network approach to explore the propagation pathway of the teleconnection between Tibetan Plateau snow cover instability and Amazon rainforest dieback.

Translating real-world ecological, climate, or social system data into the ball-and-cup framework presents significant challenges. To address the resilience of 'what to what,' both a state variable and a control parameter must be identified, as discussed in section 2. Directly mapping data onto the S-shaped curve is complex, but we can evaluate whether the system fits one of two cases: (1) stable but fluctuating, or (2) undergoing a regime shift.

For case (1), where fluctuations occur around a stable state, multi-year averages can be used to approximate the steady state, while larger fluctuations help define perturbations. From there, we can assess resilience by examining the system's recovery capacity and resistance (Yi and Jackson 2021). For case (2), where regime shifts may be occurring, sliding window averages of time-series data can reveal transitions between stable states. Examples include the collapse of Saharan vegetation (Scheffer and Carpenter 2003) and regime shifts in the Pacific Ocean ecosystem (Hare and Mantua 2000).

In addition to identifying regime shifts from time-series data, an important approach involves detecting the characteristics of tipping behavior, as discussed in section 6.3. Key features include critical slowing down (Prettyman et al 2022), early warning signals (Scheffer et al 2009), delayed transitions or 'ghost' states (Vidiella et al 2018), Turing bifurcation (Rietkerk et al 2021), and tipping cascades. Techniques such as spectral analysis, time-delay embedding, and phase-space reconstruction provide valuable methods for translating raw data into a dynamical systems framework that aligns with the ball-and-cup model.

High-dimensional systems, characterized by numerous interacting variables, present significant challenges for analysis and interpretation. Simplifying these complex models is crucial to understanding emergent behaviors like tipping points and regime shifts. A key strategy in model reduction is identifying the dominant variables or modes that capture the system's essential dynamics while discarding less critical components. By isolating these core drivers, the system can be projected onto a lower-dimensional space where its fundamental behaviors become more tractable.

Dimensionality reduction techniques, such as principal component analysis (Jolliffe 2002) and manifold learning (Tenenbaum et al 2000), offer powerful tools for this task. These methods enable the identification of key patterns and correlations within the data, which can then be used to create reduced-order models that retain the essential non-linearities and bifurcations responsible for bistable or multistable behavior.

Additionally, non-dimensionalization (Barenblatt 1996) and timescale separation (Haken 1983) are effective in simplifying models by identifying slow and fast processes, allowing researchers to focus on the critical long-term dynamics while approximating or ignoring faster, transient effects. Through techniques like center manifold reduction (Carr 1981) or adiabatic elimination (Guckenheimer and Holmes 1983), complex models can be reduced to a manageable form while preserving their ability to capture tipping point dynamics.

Future research should prioritize the development of reduction methods that not only simplify models but also retain the non-linear feedbacks and threshold effects that are critical for understanding resilience and tipping points. Such reduced models can serve as a valuable tool for exploring how high-dimensional systems transition between stable states, providing insights into critical transitions in ecological or climate systems.

Feedback loops are fundamental in shaping system resilience. Positive feedbacks can destabilize a system, driving it toward tipping points, while negative feedbacks work to stabilize the system and restore balance. Without feedback mechanisms, the S-shaped curve that is commonly used to represent resilience dynamics would not exist; the relationship between the steady-state and the control parameter would instead be linear. Introducing feedback processes to generate the characteristic S-shaped curve is both crucial and challenging, as it requires identifying the specific mechanisms that amplify or dampen system responses.

Managing feedback networks requires identifying processes that can be influenced to maintain system resilience. Anderies et al (2002) illustrate this in a fire-driven rangeland system, where positive feedback between grazing pressure and shrub dominance drives the system toward a less desirable, shrub-dominated state, reducing resilience. Conversely, negative feedback from controlled fire suppresses shrub growth, allowing the system to recover to a grass-dominated state and expanding the basin of attraction of the desirable high-biomass steady states. By adjusting grazing intensity and fire frequency, managers can disrupt destabilizing feedbacks, reduce perturbations, and prevent the system from tipping into a degraded low-biomass steady states.

Random fluctuations, or 'noise,' play a critical role in the dynamics of systems, potentially pushing them towards tipping points or delaying their approach. As discussed in sections 5.2, 5.3 and 6.2, noise influences the shape of the potential well (section 5.3), the tipping time and potential depth (section 6.2), and the overall stochastic nature of resilience (section 5.2). Understanding these effects requires careful analysis of how external variability interacts with system stability, particularly near critical thresholds.

A classic example of noise-induced transitions is provided by Nicolis and Nicolis (1981), who demonstrated how stochastic perturbations in nonlinear dynamical systems could induce shifts between stable states. In bistable systems, noise reduces the effective depth of the potential well, increasing the likelihood of transitions to alternative states. This phenomenon is particularly relevant for ecological systems, where random environmental factors, such as climate variability or disturbances like wildfires, can destabilize ecosystems and drive them toward regime shifts.

Future research could delve deeper into stochastic differential equations and noise-induced transitions to explore how random variability affects system behavior in both ecological and coupled human-natural systems. Specifically, studies could examine the role of climate variability, economic shocks, or other external drivers in modulating resilience, either reinforcing stability or pushing the system toward a critical transition. By quantifying how noise interacts with control parameters, researchers can better predict how random perturbations affect the likelihood and timing of tipping points, shedding light on the stochastic aspects of resilience and stability.

Tipping points often occur at localized or smaller scales but can propagate to trigger large-scale system changes. For instance, local desertification might spread to regional ecosystems, leading to broader climate shifts. Understanding the interactions between scales is essential for predicting large-scale transitions. Local positive feedbacks such as between plant density and rain water infiltration could be intrinsically linked to regional positive feedbacks such as between albedo and rainfall, thereby influencing the large-scale climate (Rietkerk et al 2011). Future research could focus on multiscale modeling approaches that link micro-level processes (e.g. local species extinction) to macro-level outcomes (e.g. global biodiversity loss). Network theory and spatial modeling could provide frameworks to understand how smaller tipping elements interact to produce emergent, large-scale system changes. An important area of future research is how stabilizing and destabilizing interactions between potential earth tipping elements in the climate system (Lenton et al 2008) could increase or decrease the probability of tipping and invoke or counteract tipping cascades (Wunderling et al 2024). Although helpful as a first analysis, simple conceptual models cannot fully capture the intersystem and large-scale emerging stability of such complex systems. Thus, it is still unknown how earth system elements interact in more complex Earth System Models (Drijfhout et al 2015) and in the real global Earth system.

Many systems exhibit complex behaviors across multiple scales, from local subsystems to global networks. Developing multi-scale models that effectively bridge these scales is a significant challenge. Such models should account for how isolated bistable dynamics (e.g. within plant communities, water resources, or microbial populations) couple and interact to influence larger-scale phenomena (e.g. regional climate systems or ecosystem shifts).

Future work could explore the dynamics of individual subsystems in isolation and how their coupling drives emergent properties at broader scales. Tools like hierarchical modeling (Cressie and Wikle 2011), scale-transition theory (Chesson 2018), and cross-scale feedback analysis (Levin 1998) are promising approaches for linking local processes to global trends. For example, understanding how tipping events in localized bistable systems aggregate and propagate through networks could reveal new insights into system-wide resilience and critical transitions.

By explicitly considering both the isolated behavior of subsystems and their interactions, multi-scale models can provide a more comprehensive framework for studying the interplay between small-scale and large-scale dynamics in complex systems.

Turing bifurcations before tipping can lead to spatial pattern formation within so-called Busse balloons, preventing tipping. This is because spatial pattern formation may persists within such Busse balloons well beyond tipping points (Rietkerk et al 2021). Such persistence may also be the case for so-called coexistence patterns, another kind of spatial pattern formation, where two stable states coexist in the same spatial domain. Crucial questions are how these types of spatial pattern formation can adapt in the face of ongoing environmental change and how this is affected by environmental spatial heterogeneity, noise and rate of environmental change.