<figure><img src="image_1.png" /></figure>

HAL
open science

# Sustainability or societal collapse, dynamics and bifurcations of the handy model

Arnaud Tonnelier

<figure><img src="image_2.png" /></figure>

## To cite this version:

Arnaud Tonnelier. Sustainability or societal collapse, dynamics and bifurcations of the handy model. SIAM Journal on Applied Dynamical Systems, 2023, 22 (3), 10.1137/22m1494336. hal-03659720v2

HAL Id: hal-03659720

https://hal.science/hal-03659720v2

Submitted on 18 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SUSTAINABILITY OR SOCIETAL COLLAPSE, DYNAMICS AND BIFURCATIONS OF THE HANDY MODEL.

## ARNAUD TONNELIER*

**Abstract.** The overexploitation of natural resources by our industrial society questions its long-term sustainability. Recently, a simple nature-society interrelation model, called the HANDY model, has been proposed by Montesharrei et al (2014) to address this concern with a special emphasis on the role of the stratification of the society. In this paper we analyse the dynamics of this model and we explore the influence of two parameters: the nature depletion rate and the inequality factor. We characterize the asymptotic states of the system through a bifurcation analysis and we derive several quantitative predictions. We show that some collapses are irreversible and, depending on the wealth production factor, bistability regimes can be obtained. In particular, a sustainable equilibrium can coexist with cycles of prosperity and collapse. We discuss the possible policies to prevent undesirable fates.

**Key words.** Dynamical system, Bifurcation, Human-Nature Model, Sustainability, Collapse.

**AMS subject classifications.** 37C25, 37C27, 37G15, 37N40

**1. Introduction.** There is a serious concern about the impact of human activities on the environment and several indicators show an increase of natural resources scarcity questioning the long-term viability of our society. The interplay between population and natural resources has been identified as a key feature to understand how a society can reach an equilibrium with a possible overshoot. Agricultural and industrial activities exploit natural resources at a rate faster than its regeneration capacity. In turn, the decrease of natural resources results in a deterioration of human welfare and standards of living with a risk of society instability. This raises the question of a possible societal collapse, a phenomenon that has been observed in numerous civilizations throughout history [4, 11].

Detailed human-nature models have been developed to describe how a society enters a route to collapse but their complexity often prevents any in-depth mathematical analysis [8]. Simplified models based on predator-prey type dynamics provide a simple framework to understand the dynamics of human societies [3]. However, unlike animals, the human population does not interact directly with its environment and the pressure on natural resources is exerted via industrial production and agricultural exploitation. These outputs can be aggregated into a simple variable, the so-called accumulated surplus or Wealth. The wealth is not equally distributed in the population and a small portion of the population, called the Elite, controls the distribution whereas the mass of the population is only allocated a portion of the accumulated surplus. These features have been implemented in the HANDY (Human And Nature Dynamics) model recently proposed by Motesharrei et al. (2014). Since then the HANDY model has attracted great attention and has been largely discussed in different scientific areas.

The HANDY model aims to elucidate the role of the economic stratification on the possible societal collapse. The population is divided into two classes, the Commoners, or workers, who collect (harvest, hunt) and transform the natural resources and the dominant class, the Elites, who control the wealth distribution. The HANDY system is a minimal class-stratification model that can reproduce different historical observations such as societal collapse, smooth transition to a sustainable equilibrium

<footnote>*Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France.</footnote>
or cycles of rise-and-collapse. Numerical simulations [9] reveal two types of collapse that have been named type-L collapse and type-N collapse that correspond to a labor exhaustion with a nature recovering and to a complete exhaustion (population and nature), respectively. These scenario provide a plausible explanation of civilisation collapses with the possible paths that the society could have taken. Collapse is used here as the extreme case, the full collapse, where all the population dies out. A broader description retains a significant decrease of the population size as it is the case when society experiences an overshoot.

The goal of this work is to investigate the dynamics of the HANDY model and to study the different scenarios proposed by Motesharrei et al [9] as illustrated in Fig. 1.1. A special attention is given to the study of the long-term dynamics of the model. We aim to derive mathematical results to support and to understand the numerical simulations already presented on the possible asymptotic states. We focus on the resource exploitation and the economic stratification in order to understand how the human-environment interaction can lead to a population crash. Historical observations and numerical simulations show that very similar environments can lead to radically different society outcomes motivating a bifurcation analysis of the model. We combine qualitative exploration of the dynamics with quantitative results that will be discussed in regards to the typical case of Easter island. Overall, our study shows how an unequal wealth distribution leads to collapse and how the resource exploitation influences the long-term development of a society.

The paper is organized as follows. In section 2 we present the model in details with a rescaling on the original parameters and variables. Preliminary results on dynamics and bifurcations are presented. The next sections are devoted to the different scenarios that describe different societies. In section 5 we summarize our results and discuss their possible implications.

**2. The HANDY model.** The HANDY model describes the interaction between the human population and the natural resources. The formal structure of the model originates from the Lotka-Volterra predator-prey system [7, 13] that has been reformulated by Brander and Taylor (1998) to model the Easter Island's population dynamics. The HANDY model includes two additional features: an economic stratification of the population and a wealth production process. The HANDY model is a four dimensional dynamical system that gives the evolution of the poor and rich populations (Commoners and Elites), noted $x_c$ and $x_e$, respectively, the natural resources or Nature, noted $y$, and the accumulated Wealth, $w$. The model equations are the following:

$$ (2.1) \qquad \begin{cases} \frac{dx_c}{dt} = (\beta_c - f(w, x_c, x_e))x_c, \\ \frac{dx_e}{dt} = (\beta_e - f(\kappa w, x_c, x_e))x_e, \\ \frac{dy}{dt} = \gamma y(\lambda - y) - \delta yx_c, \\ \frac{dw}{dt} = \delta yx_c - g(w, x_c, x_e)(x_c + \kappa x_e). \end{cases} $$

The birth rates of the two populations, $\beta_c$ and $\beta_e$, are assumed to be constant but the death rates depend on wealth and are described with the function:

$$ f(w, x_c, x_e) = \alpha_m + (\alpha_M - \alpha_m) \left(1 - \frac{w}{w_{th}}\right) h(w_{th} - w), $$

where $h$ is the Heaviside step function

$$ h(x) = \begin{cases} 1, & \text{if } x > 0, \\ 0, & \text{if } x \le 0. \end{cases} $$
<figure><img src="image_4.png" /><figcaption>FIG. 1.1. Population evolution of the HANDY model. The total population (in blue) is composed of two classes: Commoner (fulled black line) and Elite (dotted line). We assume that at a given time (that is arbitrarily choosen at t = 800 years) the distribution of wealth changes within the population. When t < 800, the wealth is equally distributed between the two classes (equitable society) and after a period of prosperity a sustainable equilibrium is reached (shaded grey region). At time t = 800, indicated by a red vertical line, the society becomes unequal and a full collapse is observed. Variables and parameters are taken from the original HANDY model [9] where the inequality factor is switched at time t = 800 from 1 to 2.</figcaption></figure>

Alternatively, one can use the sign function such that $h(x) = (1+sgn(x))/2$. Recently, a smooth version of the model has been proposed and studied [10]. The wealth parameter $w_{th}$ is a threshold value below which famine and deprivation start and is defined by

$$w_{th} = \rho(x_c + \kappa x_e),$$

where $\rho > 0$ stands for a consumption per capita. The two constant parameters $\alpha_m$ and $\alpha_M$ ($\alpha_M > \alpha_m > 0$) are the normal and the famine death rates, respectively. The per capita death rates are constant and equal when the accumulated wealth is sufficiently high, i.e. greater than a threshold value that is different for the two populations (weighted by $\kappa \ge 1$ for the Elite)). When the wealth is not sufficient to guarantee the well-being of the population (food is scarce, i.e. $w < w_{th}$), the death rates decrease proportionally to $w/w_{th}$. Note that the death rate of the poor population is always greater than the one of the wealthy part of the society.

The natural resources are aggregated into a single variable, the Nature. The regeneration of the natural resources is described by the logistic function with a saturation at $\lambda > 0$. Parameter $\gamma > 0$ is the regeneration factor. The depletion term $\delta yx_c$ is proportional to both Nature and the number of workers (Commoners) and is modulated by a rate of depletion per worker noted $\delta > 0$. Elites are not involved in the extraction and valorisation of resources.

The gathering and transformation of natural resources are done by the workers and participate to the wealth production. Depletion of natural resources leads to a positive flux of accumulated wealth. The accumulated wealth decreases because of the
<table><caption>TABLE 2.1
Description of the parameters of the HANDY model with their values.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Estimate</th></tr></thead><tbody><tr><td>α<sub>m</sub></td><td>Minimum death rate</td><td>0.01</td></tr><tr><td>α<sub>M</sub></td><td>Maximum death rate</td><td>0.07</td></tr><tr><td>β<sub>c</sub></td><td>Commoner birth rate</td><td>0.03</td></tr><tr><td>β<sub>e</sub></td><td>Elite birth rate</td><td>0.03</td></tr><tr><td>s</td><td>Subsistence salary per capita</td><td>0.0005</td></tr><tr><td>ρ</td><td>Threshold wealth per capita</td><td>0.005</td></tr><tr><td>γ</td><td>Regeneration rate of nature</td><td>0.01</td></tr><tr><td>λ</td><td>Nature carrying capacity</td><td>100</td></tr></tbody></table>

consumption of the population and each Elite always takes $κ$ times as much wealth as a Commoner. The rate of consumption is weighted by $g(w)$ given by

$$g(w, x_c, x_e) = s + s \left( \frac{w}{w_{th}} - 1 \right) h(w_{th} - w),$$

where $s > 0$ is a subsistence salary per capita. The consumption rate has two distinct regimes: it is constant when the wealth is sufficiently high and is reduced proportionally to $w/w_{th}$ when the wealth becomes too small.

The dimensionless parameter $κ ≥ 1$ is a weight parameter that measures the degree of inequality of the society and determines the balance of class power between elites and masses. This factor takes part in the consumption of accumulated wealth for both populations and reduces the death rate of the Elite population in comparison to Commoner.

**2.1. Rescaling and Model parameters.** Let us first discuss a possible generalization of the HANDY model where the Nature dynamics is given by

$$\frac{dy}{dt} = \gamma y(\lambda - y) - \mu y x_c,$$

with $μ > 0$, the nature depletion factor, that can be different from the wealth production factor $δ$ (in the Wealth equation). The following change of variables

$$X_c = \frac{\mu}{\delta} x_c, \quad X_e = \frac{\mu}{\delta} x_e, \quad W = \frac{\mu}{\delta} w$$

allows to recover model (2.2). The use of a depletion factor different from the production factor does not qualitatively affect the dynamics. However if we increase, respectively decrease, the depletion factor of the HANDY model (i.e. we have $μ > δ$, respectively $μ < δ$), the population and the produced wealth are decreased, respectively increased, in comparison to the original HANDY model (2.2). As an alternative to the resource degradation factor that appears in the Nature dynamics, the wealth production factor can be used to study the HANDY dynamic. A rescaling of variable could be used to eliminate $δ$ in one equation but we have chosen not do this so that the role of $δ$ could be more easily interpreted.

In what follows, we use the parameter values that were proposed in the original paper [9] and are recalled in table 2.1. A change of variables (see Appendix A) allows to
consider the equivalent HANDY model

$$ (2.2) \quad \begin{cases} \frac{dx_c}{dt} = \left( 1 + 3\left(\frac{w}{x_c+x_e}-1\right)h(x_c+x_e-w) \right)x_c, \\ \frac{dx_e}{dt} = \left( 1 + 3\left(\frac{\kappa w}{x_c+x_e}-1\right)h(x_c+x_e-\kappa w) \right)x_e, \\ \frac{dy}{dt} = 0.5y(100-y) - \delta yx_c, \\ \frac{dw}{dt} = \delta yx_c - 5(x_c+x_e) + 5(x_c+x_e-w)h(x_c+x_e-w) \end{cases} $$

where parameters $\kappa$ and $\delta$ are varied to study different scenarios while other parameters are fixed. Setting parameters to numerical values simplifies the analytical treatment of the model and allows a direct comparison with the analysis made by Montesarei et al [9]. Numerical simulations (not shown) suggest that the two parameters ($\kappa, \delta$) capture the main aspects of the dynamics while other parameters mainly act on the respective time scales of the variables. However a comprehensive analysis of the model remains to be done.

The HANDY model has been used to understand the collapse of multiple types of pre-industrial societies and has been proposed to explain Polynesian, Mayan or Greek population crash. It can provide a guide for present society that faces to resource constraints and population growth. The HANDY model is consistent with several historical records but its quantitative predictions were mainly applied to the historical case of Easter Island. We will discuss the possible consequences of our results in this case. In our rescaled version of the model (2.2), the time is expressed in units of standard variation rate (($\beta_c - \alpha_m$)$^{-1}$ = 50) so that one unit of time in the rescaled model (2.2) corresponds to 50 years. The Commoner population, $x_c$, and the equivalent Elite population, $\kappa x_e$, are recovered multiplying by a factor $\rho^{-1}$ = 200 the values of the corresponding variables in the rescaled model. The original $\delta$ parameter is recovered multiplying by a factor $10^{-4}$. Without loss of generality we set $\lambda$ = 100 (using the change of parameters $\tilde{\delta}$ = $\frac{\lambda\delta}{\rho(\beta-\alpha_m)}$ and $\tilde{\gamma}$ = $\lambda\gamma$). As previously mentioned, it is also possible to set $\delta$ = 1 in the Nature equation. We have decided to take the nature carrying capacity at 100 and to keep the depletion factor in the Nature equation mainly for numerical reasons and to facilitate the interpretation of the results. The nature variable is not rescaled and a direct comparison to the simulations previously published can be done.

The vector field is continuous but its derivative has a jump discontinuity on the two planes $w = x_c + x_e$ and $\kappa w = x_c + x_e$. Both switching planes are crossed transversally and they do not generate sliding solutions [1]. The accuracy of the numerical scheme is slightly affected at the crossing point but the global order of the scheme is preserved (simulations not shown). In this paper, simulations are done with an adaptive fourth order Runge-Kutta scheme.

## 2.2. Overview of the model dynamics and bifurcations.
We examine the dynamics of the HANDY model using $\kappa$ and $\delta$ as parameters. Changing parameter $\kappa$ allows to explore three distinct types of society. The egalitarian case where we formally take $\kappa$ = 0 in the equations describes a society without Elites (or equivalently $x_e(0) = 0$). When $\kappa$ = 1, we obtain the equitable society that is composed of workers and non-workers and the case $\kappa$ > 1 models the unequal society (with Elites and Commoners). Within each type of society, we study the influence on the dynamics of the parameter $\delta$. The goal is to elucidate the long-term dynamics of the model. For this purpose, a special attention is given to the properties of the two trivial fixed points, the nature state $P_n$ = (0, 0, 100, $w$) and the desert state $P_d$ = (0, 0, 0, $w$) (the value of $w$ is not relevant and will be discussed hereafter). The convergence to $P_n$
<figure><img src="image_5.png" /><figcaption>FIG. 2.1. Bifurcation diagram that shows the fixed points or periodic orbits (maxima and minima) of the population as a function of $\delta$. Dotted lines indicate unstable states. The dark curve is for the fixed point, the green curve indicates small amplitude limit cycles and the blue curve is for large amplitude limit cycles. The points, labeled H, Flc1 and Flc2 indicate bifurcation points. The carrying capacity of the system is reached at $(\delta^*, x^*)$. The unstable branches of the limit cycle have been computed using AUTO (with its interface using XPP [5]).</figcaption></figure>

has been reported as a type-L collapse and is related to the disappearance of Labor whereas the case $y = 0$ is categorized as a type-N collapse with the exhaustion of nature. The stability of the corresponding fixed point indicates if the collapse is reversible or not.

Our investigations reveal that a rich dynamic can be obtained for the Egalitarian society whereas the Unequal society shows a stereotyped long term evolution. The different scenarios are examined in separate sections and the structure of the paper reflects this asymmetry with a great attention given to the first scenario (where $x_e = 0$). The solutions of the Egalitarian system are studied in terms of bifurcations using the depletion factor $\delta$ as a bifurcation parameter. Results overview on the dynamics of the models provided by the co-dimension one bifurcation diagram shown in Fig. 2.1. Let us briefly discuss the bifurcations obtained. A transcritical bifurcation occurs at $\delta = \delta_{Tr}$ where the two fixed points, the nature state $P_n$ and the sustainable state $P_s$, interchange their stability. For $\delta < \delta_{Tr}$ the fixed point $P_s$ does not exist (condition $x_s \ge 0$).

At $\delta = \delta_H$ the sustainable equilibrium loses its stability through a supercritical Hopf bifurcation. There exists stable limit cycle for $\delta > \delta_H$ that disappears via a fold of limit cycles at $\delta_{Flc2}$. Another stable limit cycle associated with large amplitude oscillations occurs at a lower value of $\delta$, noted $\delta_{Flc1}$, associated with a fold of limit cycles.
**3. The Egalitarian society**, $\kappa = 0$, $x_e = 0$. In the egalitarian society there is no Elites, $x_e$ is set to 0, and the HANDY model becomes:

$$ (3.1) \qquad \begin{cases} \frac{dx}{dt} = x + 3(w - x)h(x - w), \\ \frac{dy}{dt} = 0.5y(100 - y) - \delta yx, \\ \frac{dw}{dt} = \delta yx - 5x + 5(x - w)h(x - w) \end{cases} $$

where $x = x_c$ is the total population that is composed of workers only.

**3.1. Fixed point and stability.** The egalitarian society has two trivial steady states, $P_d = (0,0,0)$ and $P_n = (0,100,0)$, that exist for every values of the depletion factor $\delta$. The state $P_d$ is the desert state where both population and nature vanish whereas in the nature state $P_n$ the natural resources are fully restored to its maximum value.

Another fixed point, noted $P_s = (x_s, y_s, w_s)$ with $x_s \neq 0$, reported as the sustainable equilibrium, may exist in the phase space region $x > w$. We have

$$ \begin{cases} 3w_s = 2x_s, \\ 0.5(100 - y_s) = \delta x_s, \\ 5w_s = \delta x_s y_s, \end{cases} $$

that defines the following one-parameter family of fixed points

$$ \begin{cases} x_s = \frac{5}{\delta}(10 - \frac{1}{3\delta}), \\ y_s = \frac{10}{3\delta}, \\ w_s = \frac{10}{3\delta}(10 - \frac{1}{3\delta}), \end{cases} $$

The Wealth is directly related to the population ($w_s = 2/3x_s$) and, as expected, the Nature steady state strictly decreases as $\delta$ increases. Fixed points exist only for $x_s \ge 0$ that gives the necessary and sufficient condition of existence

$$ (3.2) \qquad \delta \ge \delta_{Tr} = \frac{1}{30}. $$

Expression (3.2) shows that a sustainable equilibrium is possible only if the depletion rate is sufficiently high. This condition indicates the minimal resource exploitation that can guarantee a subsistence for the population. The steady state $x_s$ is termed the *carrying capacity* for the population and represents the population level that the resources can sustain over the long term. As $\delta$ increases, the population $x_s$ grows and a maximum is reached when the so called *optimal value of $\delta$*, noted $\delta^*$, is reached. As $\delta$ further increases, the population declines and collapses. We have

$$ (3.3) \qquad \delta^* = \frac{1}{15}, $$

and the maximum carrying capacity is thus given by

$$ x_s^* = 375, $$

that gives 75000 units of people in the original (not rescaled) variables of the HANDY model. The carrying capacity is maximized when $y_s = 50$ that is the point of nature's maximal regeneration rate. A generalization of the two expressions (3.2) and (3.3) for arbitrary model parameters is proposed in appendix B.

The two fixed point $P_d$ and $P_n$ are related to the collapse of the society whereas $P_e$
corresponds to a sustainable equilibrium. The stability indicates if these states are
reversible or not.

Near the desert state $P_d$, the Nature dynamics is described by the linearized equation $\dot{y} = 50y$ showing that the state is unstable so that the nature always regenerates. Therefore the type-N collapse in the HANDY model is a reversible collapse and after a while, the nature starts to recover. The analysis of the system at small $y$ values (see Appendix C) reveals that the wealth decreases first and is followed by a population decline. When the population is sufficiently small, the nature regenerates.

**3.1.1. Stability of the nature state.** The fixed point $P_n$ is singular because it lies on the discontinuity plane $x = w$. Its stability is determined by the piecewise linear system

$$ (3.4) \qquad \begin{cases} \frac{dx}{dt} = x + 3(w-x)h(x-w), \\ \frac{dy}{dt} = -50(y-100)-100\delta x, \\ \frac{dw}{dt} = 100\delta x - 5x + 5(x-w)h(x-w). \end{cases} $$

Let us focus on the dynamic in the $(x, w)$ plane that is given by

$$ (3.5) \qquad \begin{cases} \frac{dx}{dt} = 3w - 2x, \\ \frac{dw}{dt} = -5w + 100\delta x, \end{cases} $$

for $x \ge w$, and

$$ (3.6) \qquad \begin{cases} \frac{dx}{dt} = x, \\ \frac{dw}{dt} = -5x + 100\delta x, \end{cases} $$

otherwise. We distinguish between three different configurations that are given by the respective location of the two nullclines $w = 2/3x$, $w = 20\delta x$, and the threshold line $w = x$. The geometrical analysis is presented in Appendix C. The asymptotic stability of the nature state holds for

$$ (3.7) \qquad \delta \le \frac{1}{30}. $$

From the eigenvectors analysis, we find the following asymptotic relations between the population and the Wealth:

$$ (3.8) \qquad w \sim \frac{-3 + \sqrt{9 + 1200\delta}}{6}x, $$

when $\delta \le 0.06$, and

$$ (3.9) \qquad w \sim (-5 + 100\delta)x, $$

otherwise, that shows the dependence between the depletion factor, the Wealth and the population at small $y$ values.

**3.1.2. Stability of the sustainable equilibrium.** The Jacobian matrix of the HANDY system near $P_s$ is given by

$$ \begin{pmatrix} -2 & 0 & 3 \\ -\frac{10}{3} & -\mu & 0 \\ \frac{10}{3} & 50-\mu & -5 \end{pmatrix} $$
where

$$ \mu = \frac{5}{3\delta}. $$

The associated characteristic polynomial is

$$ (3.10) \qquad P(x) = x^3 + (\mu + 7)x^2 + 7\mu x - 10\mu + 500. $$

Condition (3.2) implies $\mu \le 50$. The polynomial coefficients are positive and from the Routh-Hurwitz criterion (the polynomial $a_0 + a_1x + a_2x^2 + a_3x^3$ is stable if and only if $a_1a_2 > a_0a_3$) the characteristic polynomial (3.10) is stable if and only if

$$ 7\mu(\mu + 7) > 500 - 10\mu, $$

that gives

$$ (3.11) \qquad 7\mu^2 + 59\mu - 500 > 0. $$

Condition (3.11) reads $\mu > \mu_+$ where $\mu_+$ is the largest root of the polynomial in the lefthand side of (3.11). We obtain the following necessary and sufficient condition

$$ (3.12) \qquad \delta < \delta_H := \frac{70}{3\sqrt{17481} - 177} \approx 0.32 $$

for $P_s$ to be asymptotically stable. In order to guarantee the existence of the sustainable state, it is necessary to have a depletion factor sufficiently large, above the critical value $\delta_{Tr}$. However to ensure the long-term stability of the state, this factor has to be lower than $\delta_H$ and therefore there exists a range of depletion values for which the system converges towards a sustainable equilibrium. This interval has an amplitude of order ten, i.e. we have $\delta_H \approx 10\delta_{Tr}$.

Beyond stability, an important feature of the qualitative behavior of trajectories is the existence of overshoots in the population dynamic. Overshoots generate a maximum in the population and is followed by a decay to the resting value that is interpreted as a *partial* collapse of the society. To determine under which conditions the approach towards the equilibrium corresponds to a soft landing (monotone convergence) rather than cyclical overshoots (oscillatory convergence), we examine the roots of the characteristic polynomial (3.10). This is done in appendix D and we find that an oscillatory approach to the sustainable equilibrium can occur for $\delta > 0.07$. If the damping is too important the oscillatory part of the trajectory is not observed and the adjustment path is almost monotonic. Numerical simulations reveal that an overdamped population trajectory is obtained for $\delta < 0.078$ and several overshoots are observed for $\delta > 0.12$.

**3.2. Limit cycle - Qualitative analysis.** The HANDY model exhibits sustained periodic oscillations with successive cycles of population growth and decline. This is shown in Fig. 3.1 where rescaled variables are plotted. The qualitative behavior of these solutions can be described as follows. Exploitation of natural resources by the population improves wealth per capita that in turn participates to population growth. As the population increases, stress on the nature becomes important and wealth shrinks due to nature scarcity. When the wealth is low, the population starts to decline that subsequently decreases the pressure on nature. At low population level, natural resources can regenerate and the cycling activity restarts. The different phases of the cycle are the following:
<figure><img src="image_6.png" /><figcaption>FIG. 3.1. The different phases (A-D) of the cyclic activity of the egalitarian HANDY model. Variables are scaled to unity so that only the variation of the variables over time is meaningful. Simulation is done for $\delta = 0.6$ and for an initial condition $(x_0, y_0, w_0) = (1, 100, 0)$ at time $t = 0$. The transitory regime is not shown. Note that phase D is almost pointwise.</figcaption></figure>

* Phase A: Prosperity. This phase is characterized by a population growth, a wealth increase, and a stress on nature. Let us consider that the system is near the nature state $P_n = (1, 100, 0)$. This situation represents the settlement of a small population in a preserved environment. From the analysis done near the nature state for $\delta > 0.06$, the system is locally attracted by the hyperplane (3.8) and there exists a so called hitting time at which the system reaches the switching plane $w = x$ and subsequently we have $w > x$. The Wealth dynamic is approximated by $\dot{w} = (-5 + 100\delta)x$ and the population evolution is described by $\dot{x} = x$ so that both population and accumulated Wealth have an exponential growth (we have $-5 + 100\delta > 1$).

Population growth generates a stress on the nature. The decline of natural resources is described by $\dot{y} = y/2(100 - y - \delta x)$. For $y$ large enough and $x$ small, the Nature variable evolves much faster than the Population and Wealth. Considering $y$ as a fast variable, the instantaneous adjustment of Nature to the critical manifold gives $y(t) = 100 - 2\delta x(t)$. We have $\dot{w} = x(-5 + \delta y)$ and thus the accumulated Wealth is increasing until the Nature reaches the critical value $y = 5/\delta$. At this time the system enters a new mode.

* Phase B: Population peak and nature exhaustion. This phase shows a decline in accumulated wealth, a population peak and a nature exhaustion. When the natural resources become lower than $5/\delta$, the accumulated wealth starts to decrease whereas the population is still increasing and thus there exists a hitting time at which the hyperplane $w = x$ is reached (and subsequently $x(t) > w(t)$). Right after the hyperplane crossing, the population continues to grow (we have $\dot{x} = 3w - 2x$) until the population nullcline is crossed (when $x = 3/2w$) and the population peak is reached. The natural resources are still declining and are almost totally consumed at the end of this period.

* Phase C: Full type-N collapse. The resource scarcity does not allow to restore a sufficient level of wealth and population is declining that in turn affects the wealth production. Thus, both population and wealth are decreasing and the
system approaches the desert state $P_d = (0, 0, 0)$, with $x(t) > w(t)$.

* Phase D: Fast Nature regeneration. Near the desert state, the Nature dynamics is approximated by $\dot{y} = 50y$ that implies a fast nature regeneration. The second order term in the Nature equation can no longer be neglected and the Nature dynamics is governed by the logistic growth $\dot{y} = y/2(100 - y)$. The Nature dynamics can be approximated by an instantaneous jump from $y \sim 0$ to $y \sim 100$; the system switches from the desert state to the nature state and the cycle can restart.

The trajectory during the oscillations successively lies on the two parts of the Nature nullcline, one part is associated with a fast adjustment of the Nature and is described by the surface $y = 100 - 2\delta x$, while the other is related to the slow variation of the Nature near the branch $y = 0$. Both surfaces are unstable, the former because of the positive interactions between wealth and population, and the latter due to the fast nature regeneration.

Critical questions that are not addressed by the qualitative analysis are: what are the characteristic values of the system variables in each mode?, what is the period of the oscillations?, what are the time durations of the different phases?. To tackle these issues, a quantitative analysis is presented in the next section.

**3.3. Sustained oscillations - Analytical results.** The quantitative analysis of the model exploits the piecewise linear structure of the $(x, w)$ vector field and several approximations of the nonlinear Nature dynamics.

Let us consider an initial condition $(x_0, 100, w_0)$ in a region of the state space where $x_0 \ll 1$ and $w_0 \ll 1$. The quasi static approximation of the Nature dynamics gives

$$ (3.13) \qquad y(t) = 100 - 2\delta x(t) $$

that is valid as long as $y$ is large enough. If we neglect the initiation period where $x(t) > w(t)$, i.e. we assume that the system quickly reaches the state space region $x < w$ (see Appendix C), the population dynamics is thus given by

$$ (3.14) \qquad \frac{dx}{dt} = x, $$

and the accumulated Wealth is

$$ \frac{dw}{dt} = -5x + \delta xy = x(-5 + 100\delta - 2\delta^2 x) $$

that gives, using (3.14),

$$ (3.15) \qquad w(t) = 100\delta x(t) - 5x(t) - \delta^2 x(t)^2, $$

where we neglect the terms in $x_0$ and $w_0$. In the regime where natural resources are available, the population growth exponentially and the time evolution of the accumulated wealth can be expressed as the transformation of the nature capacity minus two terms: one reflecting the nature depletion and the other representing the consumption of the population.

From (3.15), the accumulated Wealth reaches a maximum for $x = 50/\delta$ given by

$$ (3.16) \qquad w_{max} = 2500 - \frac{250}{\delta}, $$

where terms in $1/\delta^2$ have been neglected. After peak wealth, the nature is almost totally transformed and the production term can be neglected in the Wealth dynamics
so that the decrease of wealth is governed by $\dot{w} = -5x$. If we shift the time so that $t = 0$ at the wealth peak, the leading order approximation of the Wealth is given by

$$ (3.17) \qquad w(t) = 2500 - 5x(t) $$

that shows a decline of the accumulated wealth proportional to population size. The population continues to grow according to (3.14) until the threshold plane $w = x$ is reached. From (3.17), the crossing point is given by

$$ (3.18) \qquad w = x = \frac{1250}{3} \sim 416.67. $$

After this switching point, the accumulated wealth is not sufficient to ensure population growth and population begins to experience famine. The time evolution of the population is described by $\dot{x} = 3w - 2x$. If we take the switching point as the time origin, the population and the accumulated Wealth are given by

$$ (3.19) \qquad x(t) = \frac{1250}{3} (2e^{-2t} - e^{-5t}), \qquad w(t) = \frac{1250}{3} e^{-5t}, $$

respectively. The population reaches a maximum given by

$$ (3.20) \qquad x_{max} = 500 \left(\frac{4}{5}\right)^{\frac{2}{3}} \sim 430.89, $$

and the corresponding Wealth value is

$$ (3.21) \qquad w = \frac{1250}{3} \left(\frac{4}{5}\right)^{\frac{5}{3}}. $$

In comparison with the values at the threshold point, there is only a small increase of the population size (around 3%) whereas a wealth fall is observed (31%). It is interesting to observe that the maximum reaches by the population in the sustainable regime (the maximum carrying capacity) is 13% lower than the maximum reaches during the cycling activity of the system. After population peak, both wealth and population decay exponentially and we have $w < x$.

The different approximations are illustrated in Fig 3.2. In panel (a), we plot the population trajectory that is used to compute the approximations of the Nature and Wealth dynamics. The two insets show that the initiation period ($x > w$ in phase A) is very narrow and that the threshold crossing in phase B is very close to the population peak. Numerically we have $x = w \approx 452.4$ at the threshold point and we have $max_t(x(t)) = 467.9$ that should be compared to the analytical results (3.18) and (3.20), respectively. Panel (b) shows the Nature trajectory and its approximation (3.13) in phase A, the two curves almost coincide. In panel (c), we plot the accumulated Wealth and the three different approximations (3.15), (3.17) and (3.19) in the phases A, B and C, respectively. The quasi static approximation (3.13) yields an under-estimation of the Nature that implies an under-approximation of the wealth production. The discrepancy is illustrated in Fig. 3.2 (c) with the inset near the peak wealth where we have $max_t(w(t)) = 2292$ that has to be compared with $w_{max} = 2143$ given by (3.16) for $\delta = 0.7$.

The time evolutions of both the population and the accumulated Wealth illustrate the co-called Seneca effect: the decline is faster than the growth. The population
<figure><img src="image_7.png" /><figcaption>FIG. 3.2. Trajectories of (a) the population, (b) the Nature and (c) the accumulated Wealth and the corresponding approximations in the different phases of the oscillations. Dotted lines represent the computed trajectories and the associated approximations (panels (b) and (c)) are shown in dotted lines. In (a), the trajectory of the Wealth variable is also shown (red line) near the beginning of the cycle and near the peak population (see the two insets). In panel (b), the Nature and its approximation are almost indistinguishable (see the zoom in phase A). In phase B, a more accurate approximation starting at $y = 5/\delta$ (dotted horizontal line) is plotted and can be viewed more precisely in the lower inset. In (c), the Wealth trajectory and its approximation are shown in two distinct zooms, one near the transitions from phase A to B and the other near the transition from B to C. The dotted red line indicates the time at which the threshold $w = x$ is crossed (in region B). In (d), the normalized trajectories are presented with the different notations for the characteristic times of the dynamics. In (a)-(c), simulations are done for $\delta = 0.7$ and with an initial condition $(x_0, y_0, w_0) = (1, 100, 0)$. In (d) we use $\delta = 0.5$. The initial transient has been removed and the origin of time ($t = 0$) is set at the point $y = 5/\delta$ when the nature is regenerating (that graphically appears as an instantaneous jump from 0 to 100).</figcaption></figure>

decreases twice as fast as it increases (compare (3.14) with (3.19) and a factor 5 is obtained for the accumulated Wealth (compare (3.15)) in phase A with the after threshold dynamic in phase C (see (3.19)).

To describe more precisely the characteristic time scales of the dynamics (see Fig.3.2 (d)), we have the following theorem:
THEOREM 3.1. The period of the oscillations is approximated by

$$ (3.22) \qquad T \sim 45\delta $$

where the duration of the prosperity period (phase A) is given by

$$ (3.23) \qquad t_A \sim \frac{2T}{3}, $$

the duration of the population collapse (phase C) is

$$ (3.24) \qquad t_C \sim \frac{T}{3} $$

The asymptotic approximation of the inter-peak interval (phase B) is

$$ (3.25) \qquad t_B \sim \ln \delta + \ln \frac{25}{3} + \frac{1}{3} \ln \frac{5}{4} $$

*Proof.* Let $x_0$ be the value of $x(t)$ at the beginning of the cycle (phase A) and let $\tilde{t}_A$ be the time at which the system reaches the threshold point $x = w$. Using (3.14) and (3.18), we have

$$ (3.26) \qquad x_0 e^{\tilde{t}_A} = \frac{1250}{3}. $$

Let us consider $\tilde{t}_C$ the time interval from the crossing point $x = w$ (in phase B) to the endpoint of the cycle, i.e. $x(\tilde{t}_C) = x_0$. From (3.19) we have

$$ (3.27) \qquad x_0 = \frac{2500}{3} e^{-2\tilde{t}_C} $$

where the term $e^{-5\tilde{t}_C}$ has been neglected, i.e. we assume $\tilde{t}_C \gg 1$. We thus have $\tilde{t}_A = 2\tilde{t}_C$ where the term $\ln 2$ is assumed to be small comparing to $\tilde{t}_C$. By definition, the period of the cycle is $T = \tilde{t}_A + \tilde{t}_C$ and we obtain

$$ (3.28) \qquad T = 3\tilde{t}_C. $$

The time $\tilde{t}_C$ corresponds to the duration of the after threshold dynamics, i.e. the time interval where $w < x$, whereas $t_C$ is associated with the decaying phase of the population. Therefore, the time interval $\tilde{t}_C - t_C$ is the time to peak population when the origin of time $t = 0$ is taken at the threshold point. Using (3.19) we solve $\frac{dx}{dt} = 0$ and we obtain

$$ (3.29) \qquad \tilde{t}_C - t_C = \frac{1}{3}(\ln 5 - \ln 4) \sim 0.07, $$

showing that the peak population is reached just after crossing the threshold, i.e. the population begins to decline shortly after the starvation point. Therefore, $\tilde{t}_C$ is an approximation $t_C$ and using (3.28) we obtain (3.24).

Peak wealth is reached when $y = 5/\delta$. Using (3.13) and (3.14), the time duration of phase A is given by

$$ (3.30) \qquad x_0 e^{t_A} = \frac{-5 + 100\delta}{2\delta^2}. $$
Using (3.27), we find

$$ (3.31) \qquad e^{t_A} = \frac{3}{50\delta} e^{2\tilde{t}_C} $$

where we assume $\delta \gg 0.05$. We assumed that $\tilde{t}_C \gg 1$ so we can ignore the terms outside the exponentials in (3.31) and we find $t_A = 2\tilde{t}_C$. Using (3.28) we obtain (3.23).

From (3.26) and (3.30) we have

$$ (3.32) \qquad e^{\tilde{t}_A - t_A} = \frac{2500\delta^2}{300\delta - 15} \sim \frac{25}{3}\delta $$

By definition the time duration of the phase B is $t_B = \tilde{t}_A - t_A + \tilde{t}_C - t_C$. Using (3.29) and (3.32), we find (3.25).

The estimation of the period $T$ requires a more accurate approximation of the $y$ dynamics. Let us consider system at the time of peak wealth $t_A$ where the Nature is given by $y = \frac{5}{\delta}$. For $\delta$ sufficiently large we have $y \ll 50$. Thus we consider the following approximation of the Nature equation

$$ (3.33) \qquad \frac{dy}{dt} = y(50 - \delta x) $$

For $t \in [t_A, \tilde{t}_A]$, we have $x(t) = x_1 e^t$ where $x_1 = \frac{50}{\delta} - \frac{5}{2\delta^2}$. Integration of (3.33) gives

$$ (3.34) \qquad y(t) = \frac{5}{\delta} e^{50(t-t_A)+\delta x_1(1-e^{t-t_A})}. $$

The corresponding approximation is plotted in Fig. 3.2 (b) (see the inset near the transition from phase A to phase B). Using (3.32), we obtain the following approximation of the nature state at time $\tilde{t}_A$

$$ (3.35) \qquad y(\tilde{t}_A) = y_1 \sim \frac{5}{\delta} e^{-\frac{1250}{3}\delta}, $$

where we keep only the terms of order one in $\delta$ in the exponential. Now consider the dynamics of the Nature in the last time interval of the cycle. Using approximation (3.33) and the population evolution (3.19), we find

$$ (3.36) \qquad y(t) = y_1 e^{50t + \frac{1250\delta}{3}(e^{-2t}-1) - \frac{250\delta}{3}(e^{-5t}-1)}, $$

where $t=0$ is set at time $t = \tilde{t}_A$. Approximation of $\tilde{t}_C$ is obtained when the nature starts to regenerate and we fix this value at $y(\tilde{t}_C) = 5/\delta$ (the exact value does not affect the approximation since the constant terms are neglected in the asymptotic expansion of $\tilde{t}_C$). Using (3.36) we find the leading order approximation

$$ (3.37) \qquad \tilde{t}_C \sim 15\delta $$

that gives (3.22) using (3.28).

The different approximations that have been used in the proof rely on the smallness of constant parameters compared to $\delta$ (more precisely, $10\delta \ll 1$ is sufficient). The assumptions on the characteristic times can be easily checked a posteriori for $\delta$ large enough. □
During the repetitive pattern of collapse-and-regeneration, the system stays successively near the nature state $P_n$ and near the desert state $P_d$. The trajectory is trapped near these two unstable fixed points for a significant amount of time that increase linearly with $\delta$ and with a time ratio of 2/3 and 1/3, respectively. The reason of the apparent stability of the desert state relies on the smallness of the Nature variable $y(t)$ when the system is approaching $P_d$. The regeneration of the nature is fast but it requires a drastic decrease of the population and the system can leave the desert state only when the population size is sufficiently small. Conversely, at the beginning of the cycle, the stress on the nature appears to be very small due to the smallness of the population size $x(t)$. Nature starts to decrease only when the population has sufficiently recovered from its depletion. More precisely, at the end of the cycle, when society collapses, the population size, noted $x_0$, is given by (using (3.27) and (3.37))

$$ (3.38) \qquad x_0 \sim \frac{2500}{3} e^{-30\delta}. $$

The population grows again when a sufficient level of wealth is produced. Since wealth creation needs workers, the nature remains in a preserved state before suffering alteration. The fast transition from $P_d$ to $P_n$ together with the slow dynamics near these two unstable states suggest the possible existence of Canard solutions. Additional work is necessary to clarify this point but is beyond the scope of this paper.

The accuracy of the results presented in theorem 3.1 is assessed by comparing the asymptotic approximations to numerical simulations. Plot reveals that even for small $\delta$ values, the time features of the oscillations are well predicted. The main discrepancy is observed for the duration of the phase A where the approximation of $t_A$ corresponds to the total duration of the two phases A and B. Further calculations (not shown) at 0(1) order gives

$$ (3.39) \qquad T \approx 45\delta - 3\ln\delta - 10 $$

ensuring a better fit of the period. Expression (3.39) should be compared with the least squares regression line that is given by $T = 44.5\delta - 5.7$ (see Fig. 3.3).

**3.4. Consequences for the Easter Island population dynamics.** The parameters used for the simulations of the egalitarian society (see Table 2.1) are compatible with Easter Island records. In this case, as depicted in Fig. 2.1, the cyclic population dynamics is determined by two fold limit cycle bifurcations. One bifurcation occurs for a depletion rate given by

$$ \delta_{Flc2} \approx 0.33 $$

that leads for $\delta > \delta_{Flc2}$ to the oscillatory cycles previously analysed. The other fold bifurcation is at the origin of the large amplitude oscillations and occurs at a lower value of $\delta$, noted $\delta_{Flc1}$, and is approximated by

$$ \delta_{Flc1} \approx 0.23. $$

At the fold bifurcation $F_{lc1}$, oscillations have a minimal period of $T \approx 3.81$ that corresponds to 190 years (238 years is obtained using approximation (3.39)). The large amplitude oscillations give rise to a periodic pattern of prosperity-and-collapse where the maximum reaches by the population during prosperity is almost independent of $\delta$. We have $x_{max} \approx 450$ that is higher than the population capacity previously calculated
<figure><img src="image_8.png" /><figcaption>FIG. 3.3. Period of the oscillations (solid black curve) as a function of $\delta$. The grey regions indicate the duration of the different phases of the cycle: phases A, B and C, from bottom to top, respectively. The dotted black lines are the approximations of the phase durations given in theorem (3.1). The approximated durations of phase A and phase C are calculated using the period obtained from the numerical integration. The approximation for the phase A almost coincides with the total duration of the two phases, A and B, whereas phases B and C are accurately estimated. The dotted blue line is obtained from the approximation of the period given by (3.39). The green line is the least squares regression line and is given by $T = 44.5\delta - 5.7$. The insets show the two simulations (rescaled variables) of the system trajectories associated with the bounds of the delta interval, $\delta = 0.237$ (left) and $\delta = 0.5$ (right). A time simulation of $t = 10$ and $t = 30$ is considered, respectively. To fit in one plot, the population is divided by 3 and the accumulated Wealth by 7.</figcaption></figure>

$(x_s^* = 375)$. The minimum population during the cycle is almost zero indicating a full-collapse even near the bifurcation point. This regime is avoided by a policy that keeps the depletion factor below the critical value $\delta_H$. When the sustainable threshold value $\delta_H$ is exceeded, the society first experiences small-amplitude oscillations (with a fairly high frequency, $T \approx 1.0$, i.e. 50 years) before falling into large periodic variations when the critical value $\delta_{Flc2}$ is exceeded. To return in the sustainable regime, the depletion factor must be reduced but, due to the hysteresis effect, a stronger reduction of $\delta$ is necessary, up to the value $\delta_{Flc1}$, to recover the stable sustainable equilibrium. The system presents two bistability regimes. One occurs in a narrow range of $\delta$ values, $\delta \in (\delta_H, \delta_{Flc2}) \sim (0.32, 0.33)$, and is associated with the coexistence of small and large amplitude oscillations. For $\delta$ between $\delta_{Flc1}$ and $\delta_H$ (i.e. $\delta \in (0.23, 0.32)$), the system is also bistable and the society can evolve in two radically different ways. Depending on the initial condition, the trajectories converge towards a steady state or towards a limit cycle. The boundary between the two basins of attraction is shown in the $(x, y, w)$ space in Figure 3.4. The separatrix is a tunnel-shaped surface that poorly depends on the nature state except at small $y$ values where a strong deformation is observed as $y$ tends towards 0. The abrupt variation of the cross sections occurs in a small range of $y$ values (see Fig. 3.4b).

An external input on the bistable HANDY oscillator can initiate a regime shift with the transition from the limit cycle to the sustainable equilibrium and vice versa. In response to a stimulus, the system may turn from one state to another. Let us consider that the system is oscillating and that a simple impulsive control on the Wealth or on the Nature is applied in order to recover a sustainable equilibrium. As indicated by
<figure><img src="image_9.png" /><figcaption>FIG. 3.4. Bistable HANDY oscillator for $\delta = 0.29$. (a) Coexistence in the (x, y, w) state space of a large amplitude stable limit cycle (blue curve) and a stable fixed point (blue point). The red curves are cross sections of the separatrix surface of the two basins of attraction that are obtained for different values of y. The corresponding level curves are shown in the (x, w) plane in panel (b) and the associated y values are indicated. The fixed point is also shown (we have $y_s \approx 11.5$). In (a), the green lines show parts of the limit cycle where it is possible to switch from the limit cycle to the fixed point using an impulsive control on w (in the flyover zone of the separatrix) or on y (in the collapse phase where y and w are small). Another view on the dynamic is shown in the inset.</figcaption></figure>

the green zone on the limit cycle in Figure 3.4(a), there is a window of opportunity to operate such transition. For a Nature induced transition, a small additional amount of nature when the system is on the attracting manifold of the desert state can initiate a switch from the limit cycle to the fixed point. Alternatively, a strong wealth reduction in the flyover zone of the basin of attraction of the fixed point (region of the state space where the limit cycle passes over the separatrix) is able to make a transition to the sustainable region. This is illustrated in Figure 3.5 where an external control on the Nature (a) or on the Wealth (b) is applied.
<figure><img src="image_10.png" /><figcaption>FIG. 3.5. HANDY oscillator as a bistable switch. The rows show the time evolution of the population, the Nature and the accumulated Wealth, respectively. On the left, an external input is applied on the Nature that shifts $y$ from almost zero to 5. In the right column, a strong wealth reduction is carried out. The external input is applied at the time marked by a vertical dashed line and is indicated by an arrow. In both cases, the system switches from a cycle of prosperity-and-collapse to a sustainable equilibrium. Simulations are done for $\delta = 0.29$ with an initial condition on the cycle near the nature state $(x_0, y_0, w_0) = (4.8, 97.2, 16.7)$).</figcaption></figure>

An increased nature recovery in the collapse period ($y$ and $w$ small and $x$ near 200) or a severe wealth contraction near the peak wealth allow to switch the system from the cycle of prosperity-and-collapse to the sustainable equilibrium. To do this, there is a set of possible values for the external control signal that depends on the state of the system, i.e. its position on the cycle, at the time of the signal. In Figure 3.5(a), adding a small amount of nature avoid the complete collapse of the society. The severe wealth reduction shown in Figure 3.5 (b) stops the population growth and stabilizes the system to the fixed point. The control on the Nature is done later on the society cycle in comparison to the one performed on the Wealth. Conversely, it is more easy to leave the sustainable equilibrium. This is suggested by the size and the shape of the basin of attraction of the sustainable equilibrium shown in Figure 3.4 (a). Different distinct scenario destabilize the society through cycles of prosperity and decline: an increase of wealth or an almost total nature exhaustion or the arrival (or departure) of population.
**4. The non-egalitarian society, $\kappa \ge 1$.** In this section, we consider a stratified society with an inequality factor $\kappa \ge 1$. For $\kappa = 1$ the HANDY model describes an equitable society and we will show that the long-term dynamic can be directly derived from the egalitarian society. The case $\kappa > 1$ leads to the destabilization of the sustainable equilibrium and to the global stability of the nature state, i.e. the population dies out.

**4.1. The equitable society, $\kappa = 1$.** For $\kappa = 1$ the society is said to be equitable and the population can be classified into workers and non-workers. The rescaled equations that describe the evolution of the equitable society are

$$ 
\left\{ 
\begin{aligned}
\frac{dx_c}{dt} &= x_c + 3x_c \left( \frac{w}{x_c+x_e} - 1 \right) h(x_c + x_e - w), \\
\frac{dx_e}{dt} &= x_e + 3x_e \left( \frac{w}{x_c+x_e} - 1 \right) h(x_c + x_e - w), \\
\frac{dy}{dt} &= 0.5y(100 - y) - \delta y x_c, \\
\frac{dw}{dt} &= \delta y x_c - 5(x_c + x_e) + 5(x_c + x_e - w)h(x_c + x_e - w).
\end{aligned}
\right.
$$

Both population have the same per capita change rate and we have

$$ \frac{\dot{x}_c}{x_c} = \frac{\dot{x}_e}{x_e}, $$

that implies

$$ x_e(t) = r x_c(t) $$

where the ratio $r$ is determined by the initial condition, $r = x_e(0)/x_c(0)$. For $x_c(0) = 0$, it is easy to show that the system converges to the nature state. If we note

$$ x = x_e + x_c = (1 + r)x_c, $$

the total population and

$$ (4.1) \qquad \Delta = \frac{\delta}{1+r}, $$

the rescaled depletion factor then the variables $(x, y, w)$ satisfy the egalitarian HANDY model (3.1) with a depletion factor given by (4.1). Therefore, the evolution of the equitable society is obtained from the one of the egalitarian society where the depletion rate is increased by a factor $r$. This factor is determined by the initial configuration of the equitable society and increases as the initial Elite population grows. The egalitarian society presents the same asymptotic states than the equitable society but, in contrast with the egalitarian society, the interactions of the various bifurcations gives rise to different multistability regimes: bistability with the coexistence of the nature state and the sustainable equilibrium, tristability with the additional existence of a stable limit cycle and a quadrastable configuration as illustrated in Figure 4.1.

Depending on the initial condition the ratio $r > 0$ changes so that the equivalent egalitarian society, with the rectified depletion (4.1), can be shifted into different regimes. Let us consider $\delta = 0.4$. For an initial condition such that $r$ is small we have $\Delta \sim 0.4$ and the equivalent egalitarian society has a large amplitude limit cycle so that the egalitarian system also presents large sustained oscillations. The initial condition can be tuned so that the narrow regime of small oscillations is reached ($r = 0.25$ gives $\Delta = 0.32$). For an initial condition such that $r = 2$, we have $\Delta = 0.13 < \delta_{Flc}$ and the egalitarian society presents a sustainable equilibrium and thus the equitable society too. For $r = 20$, we have $\Delta \sim 0.019 < \delta_{Tr}$ and the system converges towards the nature state.
<figure><img src="image_11.png" /><figcaption>FIG. 4.1. Quadrastable egalitarian HANDY model. We use $\delta = 0.4$ and four different initial conditions are considered for the Elite population. A zoom near the origin is shown in the inset. The red curve shows a full collapse while the green curve indicates that a sustainable equilibrium is reached. The blue curve represents small oscillations originating from a supercritical Hopf bifurcation of the fixed point. The black curve show repetitive cycles of collapse-and-regeneration.</figcaption></figure>

**4.2. The unequal society.** Model (2.2) with $\kappa > 1$ represents the unequal society. For parameter values used in [9], the HANDY model is *Elite-dominated*, i.e. the rate of change of individuals in the Elite population is always greater than the one of the Commoner population. Akhavan and Yorke (2020) proved that for a large class of Elite-dominated models, the commoner population always collapses, i.e. $x_c(t) \to 0$. The proof is based on three arguments : the existence of a trapping region in the phase space, the non existence of a (strictly positive) equilibrium in this region, and the existence of a Lyapunov function. More precisely, we have

$$ \left(\frac{x_c}{x_e}\right)' = \left(\frac{x'_c}{x_c} - \frac{x'_e}{x_e}\right) \frac{x_c}{x_e}, $$

so that the ratio $x_c(t)/x_e(t)$ is strictly decreasing when the per capita change rate for the Elite population is greater than the per capita change rate of the Commoner population. Models considered in [2] include a modified HANDY system where a leak term is added in the Wealth equation. Numerical simulations suggest that both populations decline to extinction in the original Elite-dominated HANDY model. The economic stratification of the society with an unequal wealth distribution induces the global stability of the type L collapse. The only way to maintain a population over time is to decrease the inequality factor. In the following, we give more details on the dynamics of the unequal society.

Clearly, the desert state $P_d = (0,0,0,0)$ and the nature state $P_n = (0,0,100,0)$ are both trivial fixed points. No other fixed point exists in the region $x_c + x_e \le w$. The desert state $P_d$ is unstable because of the self-regeneration of the nature. The sustainable equilibrium $P_s = (x_{c,s}, 0, y_s, w_s)$ where $(x_{c,s}, y_s, w_s)$ is given by (3.2)) is also a fixed point that lies in the region $x_c + x_e < \kappa w$ of the phase space. This fixed point exists only when $\kappa > 3/2$ (condition $x_{c,s} < \kappa w_s$). Since the Elite dynamic is given by $\dot{x}_e = x_e$, the Elite population is growing and the fixed point is unstable. Therefore, the existence of an Elite population yields to the destabilisation of the
sustainable equilibrium. It is easy to see that there is no equilibrium in the region $x_c + x_e > \kappa w$.

Let us now determine the behavior of the system near the nature state $P_n$. The linearized system near the fixed point is not defined in the region $x_e + x_c > w$. Motivated by the results on Elite-dominated models [2], we assume that the dynamics near the nature state is dominated by the Elite population, i.e. $x_e \gg x_c$ and $x_e > \kappa w$. Thus, we approximate the HANDY model by

$$ \begin{cases} \dot{x}_c = -2x_c, \\ \dot{x}_e = -2x_e + 3\kappa w, \\ \dot{w} = 100\delta x_c - 5w, \end{cases} $$

and it is straightforward to show that $(x_c(t), x_e(t), w(t)) \to (0, 0, 0)$. The size of the Commoner population with respect to the produced wealth is given by

$$ \frac{w(t)}{x_c(t)} \to \frac{100\delta}{3}. $$

The evolution of the ratio between the Elites and the Commoners is described by

$$ \frac{x_e(t)}{x_c(t)} = 100\delta\kappa t + \frac{x_e(0)}{x_c(0)}, $$

showing that the Elite population is expected to dominate with a growing rate causing society to collapse.

Numerical simulations indicate that the nature state is globally stable and, depending on parameter values, long term oscillations can be observed before population dies out. When the initial Commoner's population is sufficiently large, a cyclical activity can be observed with several peaks that eventually leads to the irreversible society collapse (see Fig. 4.2 for $\kappa = 2$, $\delta = 0.3$, and $(x_c(0), x_e(0), y(0), w(0)) = (100, 1, 100, 0)$). This pattern (damping cycles of prosperity and crisis) is easily obtained for $\delta$ sufficiently large and for $\kappa$ close to 1 (almost equitable society).

Let us conclude with few remarks: the Commoner's population collapse occurs independently of the Nature or Wealth dynamics and is inherent to the Elite domination in the model. For $0 < \kappa < 1$, the society can be thought of as a two classes society with workers and non-workers where the non-workers have deteriorated living conditions and we have $x_e/x_c \to 0$.

**5. Discussion.** The rigorous analysis of detailed models of nature-society interactions reveals to be considerably complex and mathematically tractable models are desirable for a better understanding of the dynamics and to increase our ability to identify the crucial elements to avoid society collapse. In their pioneering work, Brander and Taylor (1998) have proposed a simple model that exploits the similarities between human-nature interrelations and prey-predator interactions with a particular emphasis on the role of the regeneration rate (or the intrinsic growth) of renewable resources. In contrast with two species competition models, the interaction in the HANDY model is mediated by an additional variable, the so called accumulated Wealth, and the focus is placed on the effects of the depletion rate on the dynamics. In addition, the HANDY model incorporates another salient feature of the society: the economic stratification.

Similarly to the Brander and Taylor (BT) model, different stable fixed points can be obtained in the HANDY model with either a monotonic adjustment or cyclical
<figure><img src="image_12.png" /><figcaption>FIG. 4.2. Unequal society. (a) Plot of the two classes of population (blue is for the Commoner and black for the Elite), the Nature (in green) and the accumulated Wealth (in red). After a transient where the system presents several peaks, a full collapse of the society is observed. Before oscillatory death, the system shows different cycles with an increasing part of the Elite's population whereas the accumulated Wealth presents a maximum followed by a Commoner's population peak. In (b) the evolution of the ratio between the Commoner and the Elite is shown.</figcaption></figure>

overshootings, a feature that has been proposed as an explanation of the Easter Island declines compared with the evolution of other Polynesian Islands. For all social stratifications, type N collapses, i.e. exhaustion of nature and population extinction, are reversible because of the self-regeneration of the nature inherited from the BT assumption on natural resources. However, unlike the HANDY system, the BT model does not capture a recurrent pattern observed in the history of civilizations: the periodic cycles of rise-and-collapse. In the one-class HANDY model (egalitarian society, i.e. without an Elite population) this pattern repeats indefinitely for a sufficiently strong nature-wealth transformation factor. As the production wealth factor changes, the one-class HANDY model presents different asymptotic states: the nature state, a sustainable equilibrium, small oscillations and large cycles of prosperity and collapse. For moderate values of the production factor, there are two bistability regimes. One
is associated with the coexistence of a sustainable equilibrium and a large amplitude
limit cycle. A perturbation on a variable can cause the system to initiate or terminate
oscillations. In the second bistability regime, the society first experiences small am-
plitude oscillations as a premise of a more dramatic pattern of collapse-and-recovery
that finally occurs when the production continues to grow. The two-classes society
with an equal wealth distribution (equitable society), presents similar dynamics to
the one-class society provided that the production factor is adapted to the size of
the Elite population. In particular the production factor has to be increased in or-
der to maintain the same population size. Depending on the initial size of the Elite
population, the equitable society supports the multistability between the different
asymptotic states found for the egalitarian society. In the unequal society, a sustain-
able equilibrium can not be reached and the only way to avoid an irreversible society
collapse is to return to the egalitarian society.

The subsystems for industrial production and for the agricultural sector are not de-
scribed explicitly in the HANDY model but are aggregated into a single variable,
the accumulated Wealth. However, despite its simplicity, the HANDY model extracts
some generic properties of society collapse. At first, the natural resources play the role
of the motor of growth and population increases over time. However, this trend can-
not be sustained indefinitely due to the depletion of natural resources. The economic
variable, the wealth in the HANDY model, reaches a maximum which poorly depends
on the productivity factor but which is directly related to the complete transforma-
tion of the nature. Subsequently, a population peak is generated. The population and
wealth maxima are greater than those associated with the carrying capacity (15% and
900% greater for the population and wealth, respectively) indicating that the limits
of growth are reached and that the society is on track for a global collapse. Popu-
lation decline follows a simple pattern initiated by the scarcity of natural resources
that induces a wealth reduction that in turn affects the population. If the wealth
is reinterpreted as industrial outputs and agricultural production, the HANDY tra-
jectories are reminiscent to what is observed in more realistic models as the world3
model [8, 12]. In particular, the decreasing phase is more abrupt than the growing
phase. This property does not rely on the different time constants in the dynamics
but is rather due to the nature capacity compared to the small ecological footprint
at the beginning of the prosperity period. In the HANDY model, the prosperity-and-
collapse cycle with a population growth (2/3 of the period) followed by a population
decline (1/3 of the period) illustrates the well known Seneca effect : " fortune is of
sluggish growth, but ruin is rapid ". After crossing the critical wealth threshold, the
living conditions are deteriorating and there is only 15 years of population growth
before the society begins to fall apart.

Exploitation of natural resources generates industrial and agricultural outputs that
improve human welfare and induce population growth. In turn, population increase
creates stress on nature. The balance between positive and negative feedbacks can
become unstable and the crucial question is that of the policy to be adopted to avoid
human population collapse. Despite detailed models are more suitable to address
this issue [8, 6], some general observations can be drawn from the HANDY model
analysis. A high population size requires a significant rate of nature depletion that
automatically leads the society towards a prosperity-and-crises pattern. Therefore,
limiting the nature transformation is fundamental to avoid collapse. The deliberate
choice of reducing industrial outputs has been recognised as essential to reach sustain-
ability. However this policy is not sufficient if the wealth is not equally distributed
and prioritizing an equal wealth distribution or favorizing social mobility [2] is also
essential.

**Appendix A. The HANDY model and parameters.** The general form of the HANDY model is given by

$$ (5.1) \quad \left\{ \begin{aligned} \frac{dx_c}{dt} &= \left[ 1 + a\left(\frac{w}{x_c+x_e} - 1\right)h(x_c+x_e-w) \right] x_c, \\ \frac{dx_e}{dt} &= \left[ b + a\left(\frac{\kappa w}{x_c+x_e} - 1\right)h(x_c+x_e-\kappa w) \right] x_e, \\ \frac{dy}{dt} &= \gamma y(100-y) - cxy, \\ \frac{dw}{dt} &= cxy - d\left((x_c+x_e) - (x_c+x_e-w)h(x_c+x_e)\right), \end{aligned} \right. $$

where parameters (a, b, c, d, e) are related to the original parameters of the model as follows:

$$ a = \frac{\alpha_M - \alpha_m}{\beta_c - \alpha_m}, $$

$$ b = \frac{\beta_e - \alpha_m}{\beta_c - \alpha_m}, $$

$$ c = \frac{\delta}{\rho(\beta - \alpha_m)}, $$

$$ d = \frac{s}{\rho(\beta - \alpha_m)}. $$

The following rescaled variables have been used:

$$ \tilde{x}_c = \rho x_c, $$

$$ \tilde{x}_e = \rho \kappa x_e, $$

$$ \tilde{t} = (\beta_c - \alpha_m)t, $$

where the tilde has been dropped in (5.1) for convenience. Parameter values used in [9] (see Table 2.1) give

$$ a = 3, $$

$$ b = 1. $$

**Appendix B. Transcritical bifurcation and optimal capacity.** One can easily calculate using the original parameters of the HANDY model that the trans-critical bifurcation occurs at

$$ \delta = \frac{s(\beta_c - \alpha_M)}{\lambda(\alpha_m - \alpha_M)} $$

that gives the transcritical bifurcation for the rescaled model at

$$ \delta_{Tr} = \frac{s(\beta_c - \alpha_M)}{\rho\lambda(\alpha_m - \alpha_M)(\beta_c - \alpha_m)} $$

and its numerical value calculated in (3.2) using table 2.1.

The optimal value of $\delta$, noted $\delta^*$, can be expressed in terms of the original parameters as

$$ \delta^* = \frac{2s(\beta_c - \alpha_M)}{\rho\lambda(\alpha_m - \alpha_M)(\beta_c - \alpha_m)} $$

that gives expression (3.3) using the numerical values table 2.1.
<figure><img src="image_13.png" /><figcaption>FIG. 5.1. Vector field of the HANDY model in the $(x, w)$ plane for small $y$ values. The switching line ($x = w$) and the population nullcline are depicted in black and blue, respectively. Two representative trajectories are shown in red. In the grey region, the population grows before decreasing. After threshold crossing, a full collapse is observed where both population and wealth decline.</figcaption></figure>

**Appendix C. Stability of $P_d$ and $P_n$ in the egalitarian society.** We study the stability of the desert state $P_d = (0, 0, 0)$ and the nature state $P_n = (0, 100, 0)$ in the egalitarian society (HANDY model with $x_e = 0$). The vector field is nonsmooth at $x = w$, and classical techniques of local analysis based on the Jacobian matrix cannot directly applied. The discontinuity occurs in the equations for the population and for the accumulated Wealth whereas the Nature equation is smooth. We analyze the vector field in the $(x, w)$ plane using simple geometrical considerations.

**Desert state, $P_d$.** The linearization of the Nature dynamics near the desert state indicates that the complete exhaustion of nature is an unstable state and regeneration will occur. However, as long as $x > \frac{1}{6}(50 - y/2)$ the Nature variable is decreasing and remains small until the population itself becomes sufficiently small. To understand the evolution of the system we analyze the behavior of the population and the accumulated Wealth when $y$ is small. At small $y$ values, population and Wealth equations are given by

$$ \begin{cases} \frac{dx}{dt} = -2x + 3w, \\ \frac{dw}{dt} = -5w, \end{cases} \quad \text{for } x \ge w, \text{ and} \quad \begin{cases} \frac{dx}{dt} = x, \\ \frac{dw}{dt} = -5x, \end{cases} \quad \text{for } x \le w. $$

The corresponding vector field is shown in Fig. 5.1. Simple geometrical considerations allow to show that trajectories in the $(x, w)$ plane are attracted by the invariant line $x = 0$ and converge to the origin. The Wealth variable is always decreasing and is followed by a decrease in population. Depending on the initial condition (with respect to the population nullcline $w = 2/3x$) population shows a phase of growth before declining.

**Nature state, $P_n$.** Near the maximum nature capacity $y \sim 100$, the approximation of the HANDY model is given by

$$ \begin{cases} \frac{dx}{dt} = -2x + 3w, \\ \frac{dw}{dt} = -5w + 100\delta x, \end{cases} \quad \text{for } x \ge w, \text{ and} \quad \begin{cases} \frac{dx}{dt} = x, \\ \frac{dw}{dt} = -5x + 100\delta x, \end{cases} \quad \text{for } x < w. $$

Depending on the position of the Wealth nullcline $w = 20\delta x$ with respect to the switching line $w = x$ and the population nullcline $w = 2/3x$, we distinguish between
<figure><img src="image_14.png" /><figcaption>FIG. 5.2. Vector field of the HANDY model in the $(x, w)$ plane for a state variable $y$ near the nature capacity ($y \sim 100$). The switching line is depicted in black. The blue line represents the population nullcline and the green line is for the Wealth nullcline. The light grey region is the attracting region corresponding to the $\delta$ range indicates below the associated plot. The dotted line is the attracting eigenvector associated with the numerical value used for the plot: (a) $\delta = 0.025$, (b) $\delta = 0.04$, and (c) $\delta = 0.07$.</figcaption></figure>

three different configurations that are illustrated in Fig. 5.2. Simple geometrical considerations allow to show the stability of the nature state for $\delta < 1/30$.

In the critical case where $\delta = 1/30$, the approximated system does not allow to conclude on the stability of the nature state and higher order terms have to be considered in the local approximation of the HANDY system. From the Nature equation, it can be shown that the Nature variable converges inferiorly to 100 and the non linear term $\delta xy$ will be slightly lower than $100\delta x$ when $y$ is close to 100 that can be reinterpreted as $\delta$ slightly smaller than $1/30$. Therefore, this configuration is similar to the case depicted in Fig. 5.2(a) and the system converges towards the nature state.

**Appendix D. Oscillatory approach to the sustainable equilibrium in the egalitarian society.** We determine when the characteristic polynomial (3.10) has imaginary roots. We use the following property: the third order polynomial $a_0 + a_1x + a_2x^2 + x^3$ has imaginary roots if and only if

$$4p^3 + 27q^2 < 0,$$

where

$$p = a_1 - \frac{a_2^2}{3} \quad \text{and} \quad q = a_0 - \frac{a_1a_2}{3} + \frac{2a_2^3}{27}.$$

Using (3.10) we find

$$7436000 - 430720\mu - 17761\mu^2 + 3106\mu^3 - 89\mu^4 < 0$$

that gives the approximation $\mu > 24.31$, or equivalently,

$$\delta > 0.069.$$

**Acknowledgments.** We would like to acknowledge Camille Julliard and Vincent Acary for helpful discussions.

REFERENCES
[1] V. ACARY, O. BONNEFON, AND B. BROGLIATO, *Nonsmooth Modeling and Simulation for Switched Circuits*, Springer Verlag, 69, Lecture Notes in Electrical Engineering (2010).

[2] N. AKHAVAN, AND J. A. YORKE, *Population collapse in Elite-dominated societies: a differential equations model without differential equations*, SIAM J. Appl. Dyn. Syst. 19 (2020), pp. 1736–1757.

[3] J. A. BRANDER, AND M. S. TAYLOR, *The simple economics of Easter Island: a Ricardo Malthus model of renewable resource use*, Am. Econ. Rev., 88 (1998), pp. 119–138.

[4] J. DIAMOND, *Collapse: How Societies Choose to Fail or Succeed*, Penguin, New York, 2005.

[5] B. ERMENTOUT, *Simulating, Analyzing, and Animating Dynamical Systems*, SIAM ed., 2002.

[6] G. HERRINGTON, *Update to limits to growth: Comparing the world3 model with empirical data*. Journal of Industrial Ecology, 25 (2021), pp. 614–626.

[7] A. J. LOTKA, *Elements of Physical Biology*, Williams and Wilkins, 1925.

[8] D. H. MEADOWS, D. L. MEADOWS, J. RANDERS, W. W. BEHRENS III, *The Limits to Growth*, Universe Books, 1972.

[9] S. MOTESHARREI, J. RIVAS, AND E. KALNAY, *Human and nature dynamics (HANDY): Modeling inequality and use of resources in the collapse or sustainability of societies*, Ecol. Econ., 101 (2014), pp. 90–102.

[10] S. ROMAN, S. BULLOCK, AND M. BREDE, *Coupled societies are more robust against collapse: a hypothetical look at Easter Island*, Ecol. Econ., 132 (2017), pp. 264–278.

[11] J. A. TAINTER, *The Collapse of Complex Societies*, Cambridge University Press, 1988.

[12] W. THISSEN, *Investigations into the World3 model: The capital and resource subsystem*, IEEE Trans. Syst., Man, Cybern., 6 (7), (1976), pp. 455–466.

[13] V. VOLTERRA, *Fluctuations in the abundance of a species considered mathematically*. Nature, 118 (1926), pp. 558-560.