# Harvesting a logistic population, and critical slowing down
A pond holds a population $P(t)$ that, left alone, grows logistically:
$
\frac{dP}{dt} = rP\left(1-\frac{P}{N}\right), \qquad r, N > 0 .
$
Now take fish out. There are two natural ways to write that down, and the two are separated by about four seconds of modeling decision:
$
\textbf{Model A (constant quota):}\quad \frac{dP}{dt} = f_A(P) = rP\left(1-\frac{P}{N}\right) - H,
$
$
\textbf{Model B (constant effort):}\quad \frac{dP}{dt} = f_B(P) = rP\left(1-\frac{P}{N}\right) - HP .
$
In Model A the catch is $H$ fish per unit time, whatever the stock: a quota. In Model B the catch is $HP$, proportional to how many fish there are: fixed effort, variable haul. The letter $H$ has different units in the two models — fish per time in A, inverse time in B — and, as it turns out, the two models fail in entirely different ways.
## Model A: a quota, and a fold
Equilibria are roots of $f_A$:
$
rP - \frac{r}{N}P^{2} - H = 0 \quad\Longrightarrow\quad P^{2} - NP + \frac{NH}{r} = 0 \quad\Longrightarrow\quad P_{\pm} = \frac{N}{2}\left(1 \pm \sqrt{1 - \frac{4H}{rN}}\right).
$
Everything hinges on the discriminant. For $H < rN/4$ there are two equilibria, $P_- < N/2 < P_+$. At $H = rN/4$ they coincide at $N/2$. For $H > rN/4$ there are none. Since $f_A'(P) = r\left(1 - 2P/N\right)$ changes sign at $N/2$,
$
f_A'(P_-) = +r\sqrt{1 - \tfrac{4H}{rN}} > 0 \;\;(\text{source}), \qquad f_A'(P_+) = -r\sqrt{1 - \tfrac{4H}{rN}} < 0 \;\;(\text{sink}).
$
The threshold quota is
$
H^{*} = \frac{rN}{4}, \qquad\text{reached at}\qquad P = \frac{N}{2}.
$
Raise $H$ slowly from zero. The sink $P_+$ slides down from $N$ toward $N/2$; the source $P_-$ rises from $0$ to meet it. At $H^*$ they collide and annihilate, and one hair past $H^*$ the pond has no equilibrium at all: $f_A(P) < 0$ for every $P$, and the stock runs to zero. In the $(H, P)$ plane the equilibria trace a parabola lying on its side with its nose at $(H^*, N/2)$ — a **saddle-node**, or **fold**, bifurcation.
Three features of Model A deserve names, because Model B will have none of them.
> [!warning] What goes wrong with a quota
> **No warning.** Just below $H^*$ the pond looks healthy: it holds steady at $P_+ \ge N/2$, half its carrying capacity or better. The population level itself says nothing about how close the threshold is. The warning is in the shrinking gap between $P_+$ and $P_-$, and in something subtler that the second half of this page is about.
>
> **Finite-time extinction.** $P = 0$ is not an equilibrium of Model A: $f_A(0) = -H < 0$. A stock that falls below $P_-$ does not decay toward zero — it *hits* zero at a definite time, because you are still harvesting a pond with almost nothing in it.
>
> **No recovery.** Once the stock is at zero, lowering the quota back below $H^*$ restores the two equilibria but does nothing for the pond: it sits at $0$, below the source $P_-$, and $f_A(0) < 0$ still. The system remembers its history. That memory is called **hysteresis**.
## Model B: constant effort, and a crossing
Every term of $f_B$ carries a factor of $P$, so factor rather than reaching for the quadratic formula:
$
f_B(P) = P\left[\,r\left(1-\frac{P}{N}\right) - H\,\right] = P\left[(r - H) - \frac{rP}{N}\right],
$
$
P = 0 \qquad\text{and}\qquad P^{*} = N\left(1 - \frac{H}{r}\right).
$
With $f_B'(P) = (r - H) - 2rP/N$,
$
f_B'(0) = r - H, \qquad f_B'(P^{*}) = (r - H) - 2(r - H) = -(r - H).
$
The two derivatives are negatives of each other, so the two equilibria always have opposite stability. For $H < r$ the origin is a source and $P^*>0$ is a sink; at $H = r$ the nonzero equilibrium arrives at the origin; for $H > r$ it has passed *through* into $P < 0$, which is not a population, and the origin has become the sink. Nothing is destroyed. The two branches cross and **trade stability** — a **transcritical** bifurcation at
$
H_{c} = r .
$
> [!note] Why Model B never has Model A's problems
> The whole difference is one structural fact: **is $P = 0$ an equilibrium for every value of $H$?**
>
> In Model B, $f_B(0) = 0$ identically. The extinction state is pinned as an equilibrium no matter what you do. It can never be destroyed and it can never be created, only change stability — which rules out a fold and forces a transcritical crossing. The consequences follow one by one:
>
> - **The decline is visible.** The stable stock is $P^{*} = N(1 - H/r)$, a straight line falling from $N$ at $H = 0$ to $0$ at $H = r$. A manager watching the population sees the trouble for the entire approach.
> - **Extinction is asymptotic, not finite-time.** With $P = 0$ an equilibrium, a trajectory cannot cross it. Above $H_c$ the stock decays like $e^{-(H - r)t}$ — toward zero, never *to* zero.
> - **The failure is reversible.** Cut effort back below $r$ and the origin becomes a source again; any surviving fish grow the stock back up the same straight line to the new $P^*$. No hysteresis. (The model's remnant is a positive real number; a real pond needs an integer, and a count of zero is zero. The model's promise of recovery is only as good as the survivors.)
>
> In Model A the constant $-H$ breaks the pin: $f_A(0) = -H \ne 0$, so nothing anchors any equilibrium, and equilibria are free to meet and vanish. Same fish, same pond. The difference in outcome — reversible decline versus irreversible collapse — is entirely a consequence of how the harvesting term was written down.
<div class="applet" data-applet="harvest-models"></div>
**What "bifurcation" means here.** A bifurcation is a qualitative change in the number or stability of equilibria as a parameter is varied. The sliders above vary $H$; watch the phase lines. In Model A the open and filled circles approach each other, merge into the half-filled marker at $H = H^*$, and disappear, taking the entire green (growing) region with them. In Model B the two circles approach, pass through each other at $H = r$, and swap which is open and which is filled. Things worth trying:
- Set Model A to $H = 0.95H^*$ and look at the trajectories. Most of them still settle at $P_+$, but the trapdoor $P_-$ has risen to meet them: a modest shock now sends the pond to zero in finite time (the red × marks the moment).
- Press **ramp Model A**. One population follows the stable branch down the parabola as $H$ rises, falls off the nose, and then stays at zero all the way back down: the loop does not close. Press **ramp Model B**: the stock slides down the line and slides back up it. The loop closes.
- Set both models to their thresholds ($H = H^*$, $H = r$). In each case one equilibrium has $f' = 0$. Hold that thought.
## Critical slowing down
### What linear stability analysis promises
Near an equilibrium $P_e$ write $P(t) = P_e + \varepsilon(t)$ and Taylor-expand:
$
\dot{\varepsilon} = f(P_e + \varepsilon) = \underbrace{f(P_e)}_{=\,0} + f'(P_e)\,\varepsilon + \tfrac{1}{2}f''(P_e)\,\varepsilon^{2} + \cdots
$
Linear stability analysis keeps the first surviving term and drops the rest:
$
\dot{\varepsilon} \approx f'(P_e)\,\varepsilon \quad\Longrightarrow\quad \varepsilon(t) \approx \varepsilon_0\, e^{\,f'(P_e)\,t}.
$
The sign of $f'(P_e)$ classifies the equilibrium, and the *size* of $f'(P_e)$ is a rate:
$
\tau = \frac{1}{|f'(P_e)|}
$
is the **relaxation time** — the time for a displacement from equilibrium to shrink (or grow) by a factor of $e$. For the sink of Model A,
$
\tau_A = \frac{1}{r\sqrt{1 - 4H/(rN)}} ,
$
which is finite for every $H < H^*$ and **diverges as $H \to H^*$**. This is the subtler warning promised above. As the quota approaches the threshold, the pond does not just sit at a lower level; it recovers from every disturbance more and more slowly. The phenomenon is called **critical slowing down**, and it shows up before the collapse, in the pond's *response* to perturbation rather than in its level.
### Where the promise fails: at the fold
Exactly at $H = H^*$ the two equilibria have merged at $P_e = N/2$ and $f'(N/2) = 0$. The linearization says $\dot{\varepsilon} \approx 0$: nothing moves. That is not "neutrally stable"; it is the linearization declining to answer, because the behaviour is decided by the term it threw away. So keep it. Since $f_A$ is a quadratic in $P$, its Taylor expansion about any point terminates after the second-order term, and the expansion is exact:
$
\dot{\varepsilon} = \tfrac{1}{2} f_A''(N/2)\,\varepsilon^{2} = -\frac{r}{N}\,\varepsilon^{2} .
$
(Equivalently: at $H = H^*$ the right-hand side is a perfect square, $f_A(P) = -\frac{r}{N}\left(P - \frac{N}{2}\right)^2$.) This is separable:
$
\int \frac{d\varepsilon}{\varepsilon^{2}} = -\frac{r}{N}\int dt \quad\Longrightarrow\quad -\frac{1}{\varepsilon} = -\frac{r}{N}t - \frac{1}{\varepsilon_0} \quad\Longrightarrow\quad \boxed{\;\varepsilon(t) = \frac{\varepsilon_0}{1 + \dfrac{r}{N}\varepsilon_0\, t}\;}
$
Compare the two answers:
| | equation | solution | time to go from $\varepsilon_0$ to $\varepsilon$ |
|---|---|---|---|
| off the fold | $\dot\varepsilon = -\lambda\varepsilon$ | $\varepsilon_0 e^{-\lambda t}$ | $\dfrac{1}{\lambda}\ln\dfrac{\varepsilon_0}{\varepsilon}$ — **logarithmic** in the ratio |
| at the fold | $\dot\varepsilon = -c\,\varepsilon^{2}$ | $\dfrac{\varepsilon_0}{1 + c\,\varepsilon_0 t}$ | $\dfrac{1}{c}\left(\dfrac{1}{\varepsilon} - \dfrac{1}{\varepsilon_0}\right)$ — **inverse** in the target |
(Here $\lambda = |f'(P_e)|$ and $c = r/N$.) Exponential decay reaches $1\%$ of anything in $\tau\ln 100 \approx 4.6\tau$. Algebraic decay $\varepsilon \sim \frac{N}{r t}$ takes a time proportional to $1/\varepsilon$: to get ten times closer takes ten times longer. And the formula says something else. For $\varepsilon_0 > 0$ the solution decays; for $\varepsilon_0 < 0$ the denominator reaches zero at $t = N/(r|\varepsilon_0|)$, and $\varepsilon \to -\infty$ in finite time. The merged equilibrium is **half-stable**: attracting from above, repelling below, which is exactly what the half-filled marker on the phase line is telling you.
### A concrete case
Take $r = N = 2$, so $H^* = 1$, and write $a = H^* - H$ for the distance below the threshold. Then $f_A(P) = 2P - P^2 - H = a - (P - 1)^2$, and with $y$ for the population the model becomes
$
\frac{dy}{dt} = a - (y-1)^{2}, \qquad y^{*} = 1 + \sqrt{a}, \qquad f'(y^{*}) = -2\sqrt{a}.
$
This is the same pond, arranged so that the pictures are clean: the sink sits at $1 + \sqrt a$, the source at $1 - \sqrt a$, and $a = 0$ is the fold. Put $\varepsilon = y - y^*$ and the equation for $\varepsilon$ is *exact* and has both terms in it:
$
\dot\varepsilon = a - (\sqrt{a} + \varepsilon)^{2} = -\,2\sqrt{a}\,\varepsilon - \varepsilon^{2}.
$
The linear term has coefficient $\lambda = 2\sqrt a$; the quadratic term is always there with coefficient $1$. Which one wins depends on whether $\varepsilon$ is bigger or smaller than $\lambda$: far from equilibrium the quadratic dominates, close in the linear term does — *unless* $\lambda = 0$, when the quadratic term is all there is. This is a Bernoulli equation ($w = 1/\varepsilon$ makes it linear), and its solution for $\lambda > 0$ is
$
\varepsilon(t) = \frac{\lambda\,\varepsilon_0\, e^{-\lambda t}}{\lambda + \varepsilon_0\left(1 - e^{-\lambda t}\right)},
$
which reduces to $\varepsilon_0 e^{-\lambda t}$ when $\varepsilon_0 \ll \lambda$ and to $\varepsilon_0/(1 + \varepsilon_0 t)$ as $\lambda \to 0$. The time to bring $|\varepsilon|$ from $\varepsilon_0$ down to $\varepsilon$ is
$
t = \frac{1}{\lambda}\ln\!\left[\frac{\lambda + \varepsilon}{\lambda + \varepsilon_0}\cdot\frac{\varepsilon_0}{\varepsilon}\right] \;\xrightarrow{\;\lambda\to 0\;}\; \frac{1}{\varepsilon} - \frac{1}{\varepsilon_0}.
$
**Numbers.** Start $50\%$ above equilibrium, $y(0) = 1.5\,y^*$, and ask how long it takes to come within $1\%$ of $y^*$.
- $a = 1$: $y^* = 2$, $\lambda = 2$, $\tau = 0.5$. Time to $1\%$: $t \approx 1.76$. The linearization alone predicts $\tau \ln 50 \approx 1.96$ — close, and slightly pessimistic, because the $-\varepsilon^2$ term helps at the start.
- $a = 0.01$: $y^* = 1.1$, $\lambda = 0.2$, $\tau = 5$. Time to $1\%$: $t \approx 13.2$.
- $a = 0$: $y^* = 1$, $\lambda = 0$, $\tau = \infty$. Time to $1\%$: $t = 1/0.01 - 1/0.5 = 98$. The linearization predicts that nothing happens at all.
From $a = 1$ to $a = 0$ the return to equilibrium becomes about **fifty-five times slower** — the difference between a second and a minute — while the equilibrium itself moved from $2$ to $1$ and, to anyone watching only the level, looked perfectly steady the whole way.
<div class="applet" data-applet="critical-slowing"></div>
Things worth trying:
- Start at $a = 1$ with the error axes on **semilog**. The actual solution is nearly a straight line: exponential decay, at the slope the linearization predicts. Now walk $a$ down to $0.001$. The straight line tilts toward horizontal and bends: the early part is governed by $-\varepsilon^2$, the late part by $-\lambda\varepsilon$, and the two regimes are visible as two slopes.
- Press **a = 0** and switch to **log–log**. Now the actual solution is the straight line (slope $-1$: $\varepsilon \sim 1/t$), it coincides with the fold prediction, and the linear prediction is a horizontal line at the initial error.
- Still at $a = 0$, drag $y(0)/y^*$ below $1$. The solution runs away from $y^* = 1$ and reaches zero in finite time. Half-stable, seen directly.
- Set $a = 0.01$ and drag $y(0)/y^*$ to $0.8$, which puts $y(0) = 0.88$ below the source at $1 - \sqrt{a} = 0.9$. Extinction in finite time at a quota that is under the threshold.
### The recovery time diverges
Collect the relaxation times as a function of $a$:
$
\tau = \frac{1}{2\sqrt{a}} \propto a^{-1/2}.
$
On log–log axes that is a line of slope $-\tfrac12$. It has no ceiling: halving the distance to the fold multiplies the recovery time by $\sqrt2$, and this continues without bound. The actual time to come within $1\%$ tracks $\tau$ for $a$ not too small, then peels off toward the algebraic limit of $98$ as $a \to 0$, because the quadratic term takes over the last stretch.
<div class="applet" data-applet="slowing-time"></div>
> [!note] Why this is an early-warning signal
> A real pond is never exactly at equilibrium — weather, disease, and luck keep knocking it around. Each knock is a small $\varepsilon_0$, and the pond relaxes back at rate $\lambda = 2\sqrt a$. As the quota is turned up and $a$ shrinks, the knocks take longer and longer to fade. In a time series of the population this shows up as *rising autocorrelation* (this year's deviation predicts next year's, because it has not decayed yet) and *rising variance* (deviations pile up before they can relax). The level of the stock gives no warning; its *fluctuations* do. This is the idea behind the early-warning-signal literature for tipping points in ecosystems, climate, and finance (see M. Scheffer et al., "Early-warning signals for critical transitions," *Nature* 461, 53–59 (2009)), and it is where the second half of this discussion picks up: complex systems, and what slowing down looks like when you can only observe the noise.
---
## Summary
| | Model A: quota $-H$ | Model B: effort $-HP$ |
|---|---|---|
| $P = 0$ an equilibrium? | no: $f(0) = -H$ | yes, always |
| equilibria | $P_\pm = \frac{N}{2}\left(1 \pm \sqrt{1 - 4H/rN}\right)$ | $0$ and $N(1 - H/r)$ |
| bifurcation | saddle-node (fold) at $H^* = rN/4$ | transcritical at $H_c = r$ |
| what happens at threshold | sink and source collide and vanish | sink and source cross and trade stability |
| decline before threshold | stock holds at $\ge N/2$; no visible warning | stock falls linearly to zero; visible |
| extinction | finite time, from any $P_0 < P_-$ | asymptotic, only for $H > r$ |
| reversible? | no (hysteresis) | yes |
| $f'$ at threshold | $0$: linearization silent, decay is algebraic | $0$ at both equilibria as they cross |
The linearization is a first-order tool. It answers the stability question correctly whenever $f'(P_e) \ne 0$, and it also tells you *how fast* — $\tau = 1/|f'(P_e)|$. Watching $\tau$ grow as a parameter is pushed toward a fold is the earliest available warning that the fold is coming. Exactly at the fold the first-order tool fails, and the second-order term takes over with a qualitatively different answer: not exponential, but algebraic, and only from one side.