Chapter 2 Critical Slowing Down: Why Recovery Gets Slower

A valley growing shallow

The folds of the previous chapter have a geometric picture: the stable state a system occupies is like a valley; as pressure rises the valley grows shallower and moves towards the neighbouring ridge, and at the fold the two merge and vanish. The direct consequence of a shallower valley is that, after the same push, the ball rolls back to the bottom more slowly. This is critical slowing down.

In 1984 Wissel turned this into a “universal law”: as a system approaches a threshold, the characteristic time it takes to return to equilibrium after a disturbance grows, and the way it grows can be used to extrapolate the position of the threshold from measurements made far away from it; he tested the method with preliminary harvesting experiments on rotifer populations (Wissel 1984). In six ecological models van Nes and Scheffer found that slower recovery from small perturbations becomes apparent well before the threshold, and argued that the recovery rate (later called “engineering resilience”) can serve as an indicator of the size of the basin of attraction (“ecological resilience”) (Nes and Scheffer 2007).

The square-root law

Calculation

Near an equilibrium $x^*$ , let $y = x – x^*$ ; to first order $\mathrm{d}y/\mathrm{d}t \approx f'(x^*)\,y = -\lambda y$ . Here $\lambda = -f'(x^*)$ is the recovery rate and $T = 1/\lambda$ the characteristic recovery time. At a saddle-node bifurcation $f'(x^*) = 0$ , and near it

$\lambda \propto \sqrt{c_2 – c}, \qquad T \propto (c_2 – c)^{-1/2},$

so $1/T^2$ falls roughly linearly to zero as $c$ increases (Strogatz 2015). Wissel’s extrapolation exploits exactly this: measure the recovery time at several pressures far from the threshold, draw the straight line of $1/T^2$ against pressure, and its intercept with the axis estimates the threshold.

The program computes $\lambda$ , $T$ and $T\sqrt{c_2 – c}$ (roughly constant if the square-root law holds) for the upper stable state of the grazing model; runs two “pulse experiments” (removing 2 % or 10 % of the vegetation at once and timing 90 % recovery); and fits $1/T^2$ over three ranges of pressure to extrapolate the threshold.

Calculation

# Recovery slows near a fold: return time of the upper (vegetated) state of the grazing
# model dx/dt = r x (1 - x/K) - c x^2/(x^2 + a^2) (r = 1, K = 10, a = 1) as c -> c2 = 2.604.
# Linearisation: lam = -f'(x*), return time T = 1/lam. Near a saddle-node lam ~ sqrt(c2 - c),
# so 1/T^2 falls roughly linearly to 0 at c2: the basis of Wissel's (1984) extrapolation.
import math
R, K, A, C2 = 1.0, 10.0, 1.0, 2.604365

def f(x, c):
    return R * x * (1 - x / K) - c * x * x / (x * x + A * A)

def upper_eq(c):
    lo, hi = 4.79, 10.0                    # the upper equilibrium lies above the fold x = 4.78
    for _ in range(80):
        mid = (lo + hi) / 2
        if f(mid, c) > 0: lo = mid
        else: hi = mid
    return lo

def lam(c, h=1e-6):
    x = upper_eq(c)
    return -(f(x + h, c) - f(x - h, c)) / (2 * h)

def pulse_recovery(c, kick, dt=0.001):
    """Push x down by a share kick; time until 90 % of the displacement has decayed."""
    xs = upper_eq(c); x, t = xs * (1 - kick), 0.0
    while xs - x > 0.1 * kick * xs:
        x += f(x, c) * dt; t += dt
        if x < 3.0: return None              # fell into the lower basin: it tipped
    return t

print("   c    x*    lam   T = 1/lam  T*sqrt(c2-c)  pulse 2 %: 90 % back after  pulse 10 %")
for c in (1.0, 1.5, 2.0, 2.3, 2.5, 2.58, 2.60):
    l, big = lam(c), pulse_recovery(c, 0.10)
    print(f"  {c:4.2f}  {upper_eq(c):4.2f}  {l:5.3f}  {1 / l:7.2f}     {math.sqrt(C2 - c) / l:5.2f}"
          f"         {pulse_recovery(c, 0.02):7.2f}           "
          + (f"{big:6.2f}" if big else "TIPPED"))
# Wissel-style forecast: fit 1/T^2 = lam^2 against c using only c <= cmax, extrapolate to 0
for cmax in (2.0, 2.3, 2.5):
    top = int(cmax * 100)
    pts = [(c / 100, lam(c / 100) ** 2) for c in range(top - 40, top + 1, 5)]
    n = len(pts); mx = sum(p for p, _ in pts) / n; my = sum(q for _, q in pts) / n
    slope = sum((p - mx) * (q - my) for p, q in pts) / sum((p - mx) ** 2 for p, _ in pts)
    print(f"fit on c in [{cmax - 0.4:.1f}, {cmax:.1f}]: predicted threshold"
          f" c = {mx - my / slope:.3f}  (true {C2:.3f})")
#    c    x*    lam   T = 1/lam  T*sqrt(c2-c)  pulse 2 %: 90 % back after  pulse 10 %
#   1.00  8.89  0.781     1.28      1.62            2.98             3.09
#   1.50  8.20  0.645     1.55      1.63            3.61             3.76
#   2.00  7.32  0.473     2.11      1.64            4.92             5.18
#   2.30  6.59  0.333     3.01      1.66            7.03             7.50
#   2.50  5.84  0.192     5.20      1.68           12.25            13.56
#   2.58  5.30  0.092    10.90      1.70           26.21            32.59
#   2.60  5.00  0.038    26.00      1.72           66.09           TIPPED
# fit on c in [1.6, 2.0]: predicted threshold c = 2.584  (true 2.604)
# fit on c in [1.9, 2.3]: predicted threshold c = 2.592  (true 2.604)
# fit on c in [2.1, 2.5]: predicted threshold c = 2.598  (true 2.604)

$T\sqrt{c_2 – c}$ moves only slowly from 1.62 to 1.72: the square-root law holds approximately over a wide range. As the pressure rises from 1.0 to 2.5, the recovery time increases from 1.28 to 5.20, and at 2.60 it is 26.0; the time for 90 % recovery from a 2 % pulse rises from 2.98 to 66.1. One detail deserves attention: at $c = 2.60$ a 10 % pulse already pushes the system past the unstable equilibrium, and the system flips. The closer the threshold, the narrower the basin, and the pulse experiment itself may be the push that carries the system over.

In this noise-free model Wissel’s extrapolation works well: data from 1.6–2.0 give 2.584 and data from 2.1–2.5 give 2.598, against a true value of 2.604. But here $\lambda$ is computed exactly; in real measurements the recovery rate has to be estimated from noisy data.

Slowing down in noise

More often a system is continually hit by small random disturbances, and all we can observe is how it fluctuates.

Calculation

Derivation. The linearised system near an equilibrium plus white noise is the Ornstein–Uhlenbeck process $\mathrm{d}x = -\lambda x\,\mathrm{d}t + \sigma\,\mathrm{d}W$ , whose solution is

$x(t + \Delta t) = e^{-\lambda \Delta t}\,x(t) + \sigma \int_t^{t+\Delta t} e^{-\lambda (t + \Delta t – s)}\,\mathrm{d}W(s).$

The second term is independent of $x(t)$ and has mean zero. Multiplying both sides by $x(t)$ and taking expectations gives the lag-1 autocorrelation at spacing $\Delta t$ , $\mathrm{AC}(1) = e^{-\lambda \Delta t}$ ; by the Itô isometry the stationary variance is

$\operatorname{Var}[x] = \sigma^2 \int_0^\infty e^{-2\lambda s}\,\mathrm{d}s = \frac{\sigma^2}{2\lambda}.$

As $\lambda \to 0$ , $\mathrm{AC}(1) \to 1$ and the variance $\to \infty$ : the system “remembers” disturbances from ever longer ago, and its fluctuations grow (Scheffer et al. 2009). Conversely, $\lambda = -\ln \mathrm{AC}(1)/\Delta t$ . The program simulates the process by the Euler–Maruyama method (step 0.01), sampling once per time unit, 20,000 samples in all.

Calculation

# Critical slowing down in its simplest form: the Ornstein-Uhlenbeck process
#   dx = -lam x dt + sigma dW      (lam = recovery rate, the linearised restoring force)
# Sampled every dt_s: lag-1 autocorrelation AC(1) = exp(-lam dt_s), variance = sigma^2/(2 lam).
# Check both formulas by Euler-Maruyama simulation (step 0.01, sampled every dt_s = 1).
import math, random
rng = random.Random(3)

def simulate(lam, sigma, n=20000, h=0.01, every=100):
    x, out, s = 0.0, [], sigma * math.sqrt(h)
    for _ in range(2000):                      # burn-in
        x += -lam * x * h + s * rng.gauss(0, 1)
    for i in range(n * every):
        x += -lam * x * h + s * rng.gauss(0, 1)
        if i % every == every - 1:
            out.append(x)
    return out

def ac1_var(y):
    m = sum(y) / len(y)
    v = sum((u - m) ** 2 for u in y) / len(y)
    c = sum((y[i] - m) * (y[i + 1] - m) for i in range(len(y) - 1)) / (len(y) - 1)
    return c / v, v

print(" lam   sigma   AC(1) theory  simulated   Var theory  simulated   lam = -ln AC(1)")
for lam, sigma in ((1.0, 0.5), (0.5, 0.5), (0.2, 0.5), (0.05, 0.5), (0.5, 1.0)):
    a, v = ac1_var(simulate(lam, sigma))
    print(f" {lam:4.2f}   {sigma:4.2f}     {math.exp(-lam):.3f}       {a:.3f}"
          f"      {sigma ** 2 / (2 * lam):7.3f}    {v:7.3f}       {-math.log(a):.3f}")
print("halving lam: AC(1) rises and Var doubles; doubling sigma: Var x4, AC(1) unchanged")
#  lam   sigma   AC(1) theory  simulated   Var theory  simulated   lam = -ln AC(1)
#  1.00   0.50     0.368       0.351        0.125      0.125       1.046
#  0.50   0.50     0.607       0.609        0.250      0.252       0.495
#  0.20   0.50     0.819       0.820        0.625      0.631       0.199
#  0.05   0.50     0.951       0.954        2.500      2.636       0.048
#  0.50   1.00     0.607       0.612        1.000      1.004       0.491
# halving lam: AC(1) rises and Var doubles; doubling sigma: Var x4, AC(1) unchanged

Formula and simulation agree: as $\lambda$ falls from 1 to 0.05, the autocorrelation rises from 0.37 to 0.95 and the variance from 0.125 to 2.5. At $\lambda = 0.05$ the simulated variance is about 5 % above theory (2.636 against 2.500): the slower the recovery, the more correlated neighbouring samples are, and only a few hundred of the 20,000 samples are effectively independent. The closer the threshold, the larger the indicators, and the more data are needed to estimate them. The last row shows something else: doubling the noise multiplies the variance by four (0.252 → 1.004) while the autocorrelation barely changes (0.609 → 0.612). Rising variance alone may only mean that outside disturbances have grown. On this basis Ditlevsen and Johnsen stressed that an early warning requires variance and autocorrelation to rise together, and that conclusions drawn from only one of the two signals are invalid (Ditlevsen and Johnsen 2010).

A trajectory approaching the threshold

Calculation

Add noise to the grazing model ( $\sigma = 0.15$ ) and let the pressure rise from 1.0 to 2.8 over 1,000 time units, recording once per unit. Find the moment of collapse ( $x < 3$ ), drop the last 10 points, detrend with a Gaussian kernel (bandwidth 25), compute the lag-1 autocorrelation and variance of the residuals in a moving window half as long as the record, and measure their upward trend with Kendall’s $\tau$ (following Dakos and colleagues 2008, 2012).

Calculation

# Early-warning indicators before a fold: dx = [x (1 - x/10) - c x^2/(x^2+1)] dt + 0.15 dW,
# c rising from 1.0 to 2.8 over 1,000 time units (fold at 2.604), one sample per unit.
# Before the shift: Gaussian-kernel detrending (bandwidth 25); window = half the record;
# AC(1) and variance of the residuals; trend by Kendall's tau (after Dakos et al. 2008, 2012).
import math, random

def simulate(seed, T=1000, sigma=0.15, c0=1.0, c1=2.8, dt=0.05):
    rng = random.Random(seed); x, xs, cs = 8.9, [], []
    for t in range(T):
        c = c0 + (c1 - c0) * t / T
        for _ in range(int(1 / dt)):
            x += ((x * (1 - x / 10) - c * x * x / (x * x + 1)) * dt
                  + sigma * math.sqrt(dt) * rng.gauss(0, 1))
        xs.append(x); cs.append(c)
    return xs, cs

def detrend(y, bw=25):
    out = []
    for i in range(len(y)):
        js = range(max(0, i - 3 * bw), min(len(y), i + 3 * bw + 1))
        w = [math.exp(-0.5 * ((j - i) / bw) ** 2) for j in js]
        out.append(y[i] - sum(a * y[j] for a, j in zip(w, js)) / sum(w))
    return out

def stats(w):                                  # lag-1 autocorrelation, variance
    m = sum(w) / len(w); v = sum((u - m) ** 2 for u in w)
    return sum((w[i] - m) * (w[i + 1] - m) for i in range(len(w) - 1)) / v, v / len(w)

def kendall(y):
    s = sum((y[j] > y[i]) - (y[j] < y[i]) for i in range(len(y)) for j in range(i + 1, len(y)))
    return s / (len(y) * (len(y) - 1) / 2)

def indicators(pre, step=5):
    res = detrend(pre); win = len(res) // 2
    idx = list(range(win, len(res) + 1, step))
    a, v = zip(*(stats(res[i - win:i]) for i in idx))
    return idx, a, v

if __name__ == "__main__":
    xs, cs = simulate(11)
    shift = next(t for t, x in enumerate(xs) if x < 3.0)
    print(f"shift (x < 3) at t = {shift}, c = {cs[shift]:.3f}; deterministic fold at c = 2.604")
    idx, a, v = indicators(xs[:shift - 10])     # drop the last 10 points: the fall itself
    print(" window ends at t   c     AC(1)   variance")
    for k in sorted({round(i * (len(idx) - 1) / 6) for i in range(7)}):
        print(f"      {idx[k]:4d}       {cs[idx[k] - 1]:.2f}   {a[k]:.3f}   {v[k]:.4f}")
    print(f"Kendall tau: AC(1) {kendall(a):+.2f}, variance {kendall(v):+.2f}")
# shift (x < 3) at t = 913, c = 2.643; deterministic fold at c = 2.604
#  window ends at t   c     AC(1)   variance
#        451       1.81   0.479   0.0179
#        526       1.94   0.469   0.0185
#        601       2.08   0.513   0.0197
#        676       2.21   0.556   0.0212
#        751       2.35   0.571   0.0225
#        826       2.48   0.596   0.0232
#        901       2.62   0.737   0.0352
# Kendall tau: AC(1) +0.91, variance +0.92

A simulated trajectory approaching a fold (top) and the indicators in a moving window (bottom). The vegetation collapses at t = 913, when c = 2.64, slightly beyond the deterministic fold at 2.60. Before the collapse, the lag-1 autocorrelation rises from about 0.48 to about 0.74 and the variance from about 0.018 to about 0.035. — A simulated trajectory approaching a fold (top) and the indicators in a moving window (bottom). The vegetation collapses at $t = 913$ , when $c = 2.64$ , slightly beyond the deterministic fold at 2.60. Before the collapse, the lag-1 autocorrelation rises from about 0.48 to about 0.74 and the variance from about 0.018 to about 0.035.

The autocorrelation rises from 0.48 to 0.74 and the variance from 0.018 to 0.035, with Kendall’s $\tau$ of +0.91 and +0.92; the rise accelerates in the second half, just as the square-root law predicts. When the last window ends, $c = 2.62$ , already past the deterministic fold, yet the system is still in the upper state: the “bottleneck” left behind by the fold makes the system take a while to leave (Volume 13, Chapter 4). This is an idealised single demonstration: the system really is approaching a threshold slowly, the parameters are known, the noise is constant and the record is long. The next chapter assesses systematically when such warnings can be trusted.

Evidence

Reviews and palaeoclimate. The 2009 review by Scheffer and colleagues gathered slower recovery, rising autocorrelation and variance, changing skewness and flickering into “early-warning signals” that may apply generally to ecosystems, climate and financial markets (Scheffer et al. 2009); Carpenter and Brock showed in a model of lake eutrophication that rising variance can be a leading indicator (Carpenter and Brock 2006). Dakos and colleagues analysed eight reconstructed abrupt climate shifts (the end of the greenhouse Earth about 34 million years ago, the end of the Younger Dryas, the Bølling–Allerød transition, the desertification of North Africa about 5,000 years ago and the ends of four glaciations); autocorrelation rose before all eight, with a combined probability of $P < 0.003$ in surrogate-data tests. But the authors also reported that bootstrap tests suggested that about half of the positive trends could arise by chance, that the three strongest records had more than 450 points each while the others had fewer than 150, so that results from the shorter series are less reliable, that the detrending method affects the results, and that slowing down appears only when a system approaches a threshold gradually (Dakos et al. 2008).

Evidence

Experiments. In Peter Lake, Wisconsin, Carpenter and colleagues added largemouth bass in stages in 2008–2009 to change the food web gradually, with neighbouring Paul Lake as a reference; more than a year before the food-web transition was complete, the variance, autocorrelation and other indicators of chlorophyll in Peter Lake were already anomalous, while those of the reference lake were not. This was a controlled whole-ecosystem experiment, but with only one manipulated lake there was no true replication, and the precise nonlinear mechanism is not known (Carpenter et al. 2011). Dai and colleagues worked with budding yeast growing on sucrose: cooperative breakdown of the sucrose creates a strong Allee effect, and a daily dilution factor acts as a death rate; between dilution factors of 500 and about 1,600 the system has two stable states, survival and extinction. Near the fold the coefficient of variation rose more than fourfold and the autocorrelation time from under 2 days to over a week; skewness showed no systematic change (Dai et al. 2012).

This positive evidence has one feature in common: either the system was pushed towards the threshold deliberately and slowly (Peter Lake, yeast), or the analysis was done after the transition had happened (palaeoclimate). The latter carries an easily overlooked statistical problem, taken up in the next chapter.

Chapter summary

Near a fold the recovery rate falls to zero as a square root; in the grazing model the recovery time measured far from the threshold extrapolates the threshold fairly accurately. Under random disturbance, slower recovery appears as rising lag-1 autocorrelation $e^{-\lambda\Delta t}$ and variance $\sigma^2/(2\lambda)$ ; rising variance alone may only mean more noise. Similar signals have been seen in an experimental lake, in yeast and in palaeoclimate records. But the closer the threshold, the more dangerous the measurement and the more data it takes.

References

Carpenter, S. R., and W. A. Brock. 2006. “Rising Variance: A Leading Indicator of Ecological Transition.” Ecology Letters 9 (3): 311–18.

Carpenter, S. R., J. J. Cole, M. L. Pace, et al. 2011. “Early Warnings of Regime Shifts: A Whole-Ecosystem Experiment.” Science 332 (6033): 1079–82. https://doi.org/10.1126/science.1203672.

Dai, Lei, Daan Vorselen, Kirill S. Korolev, and Jeff Gore. 2012. “Generic Indicators for Loss of Resilience Before a Tipping Point Leading to Population Collapse.” Science 336 (6085): 1175–77. https://doi.org/10.1126/science.1219805.

Dakos, Vasilis, Marten Scheffer, Egbert H. van Nes, Victor Brovkin, Vladimir Petoukhov, and Hermann Held. 2008. “Slowing down as an Early Warning Signal for Abrupt Climate Change.” Proceedings of the National Academy of Sciences 105 (38): 14308–12. https://doi.org/10.1073/pnas.0802430105.

Ditlevsen, Peter D., and Sigfus J. Johnsen. 2010. “Tipping Points: Early Warning and Wishful Thinking.” Geophysical Research Letters 37 (19): L19703. https://doi.org/10.1029/2010GL044486.

Nes, Egbert H. van, and Marten Scheffer. 2007. “Slow Recovery from Perturbations as a Generic Indicator of a Nearby Catastrophic Shift.” The American Naturalist 169 (6): 738–47. https://doi.org/10.1086/516845.

Scheffer, Marten, Jordi Bascompte, William A. Brock, et al. 2009. “Early-Warning Signals for Critical Transitions.” Nature 461 (7260): 53–59. https://doi.org/10.1038/nature08227.

Strogatz, Steven H. 2015. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. 2nd ed. Westview Press.

Wissel, C. 1984. “A Universal Law of the Characteristic Return Time Near Thresholds.” Oecologia 65 (1): 101–7. https://doi.org/10.1007/BF00384470.

Fan, Y. (2026). Phase Transitions, Critical Slowing Down and Resilience in Complex Dynamical Systems: From Saddle-Node Bifurcations and Early-Warning Signals to Cascading Failures, Read Beside Classical Chinese Thought on Incipient Change [复杂系统相变、临界慢化与生态韧性的非线性动力学] (Version 1.0.0). 5D Library. https://doi.org/10.5281/zenodo.23770082 · Licensed under Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0).