\section{Introduction} 

Societies were always faced with challenges vastly different in nature. These include external events such as invasions and environmental catastrophes in the case of Mesopotamia \citep{Weiss.1993} or internal challenges like corruption, rebellions and mismanagement in ancient Egypt \citep{Butzer.2012}. Historically, both cases led to societal collapse. Other classic and frequently studied examples include the fall of the Western Roman Empire in the 5th century CE, the collapse of the Maya in 9th century CE or the fall of the Minoan Civilization of Crete in 14th century BCE \citep{Tainter.1988, Middleton.2012}.
Modern societies are facing similarly severe challenges, such as global climate change, pandemics, or financial instabilities. Estimated costs of climate change \citep{Tol.2018} alone will put additional strain on already excessively indebted national states \citep{WorldBank.2020}. While the real estate crisis of 2008 still reverberates in several countries, the Covid-19 pandemic caused major impacts on health systems and potential long-term consequences for economies and socio-cultural systems \citep{Trump.2020}.

The aforementioned examples therefore raise the question whether it is possible to identify underlying principles that determine a society's ability to cope with such large-scale challenges and thus increase their resilience to collapse. As reviewed by \citet{Cumming.2017}, several explanatory models for collapse exist. Explanations focusing on resource limitations such as the classical Malthusian trap \citep{Malthus.1798} and limits to exponential growth and overuse of resources have been widely discussed in the past \citep{Meadows.1972, Rockstrom.2009, Diamond.2011}.

Tainter's {\em theory of collapse of complex societies} emphasises specifically the notion of societal complexity and its tendency to self-amplifying increase in response to stress as the main cause of collapse \citep{Tainter.1988, Tainter.2006}. Within this framework, complexity emerges continuously through problem solving and manifests itself in the differentiation and specialization of social roles, hierarchies and control of behavior, growing population, technical abilities and increased information flow. 
These investments in complexity have an associated cost, such as calories, time, or money, all of which can be reduced to an abstract form of \textit{energy}. Although an increase in complexity is beneficial for societies and adds to well-being, surpluses gained by complexity increments are hypothesized to decrease as the `low-hanging-fruits' start getting less. Thus diminishing marginal returns on investments in complexity (or simply return-on-complexity $ROC$) ultimately drive a society into collapse. In addition, benefits of complexity may be short-natured and quickly consumed by population increases or gains in living standard. In contrast, elevated demands in energy are persistent, leading to new problems, once previously earned energy surpluses are used up. Within the framework proposed in~\citet{Tainter.1988} this process is denoted as the {\em energy-complexity spiral} \citep{Tainter.2006}.

Contemporary evidence for complexity gains in specialization and differentiation of social roles as well as technical abilities are the decrease of labour in the agricultural sector and increase of the finance and communication sectors in the 20th century, while simultaneously efficiency in primary sectors was increased \citep{Feinstein.1999}. Another example for complexity increase is a rising bureaucracy burden in cancer research (formalities, regulations of opening studies). The added complexity diminishes ROIs of research, thus resulting in reduced innovativeness in this field \citep{Steensma.2014}. \citet{Hall.2009} estimates the energy return on investment (EROI) of a sustainable society to be well above 3:1, which can be interpreted as a proxy for complexity. Despite these examples, overall empirical evidence for societal complexity and its change over time is scarce, although it seems ubiquitous in our lives.

Mathematical models and simulations are important tools to study the implications of theories of collapse and resilience and to test competing hypotheses. In the archaeological domain, agent-based models are employed to understand the interaction of water and land-use change \citep{Barton.2010} or to compare trajectories of biocultural evolution \citep{Barton.2012}. Dynamical systems models are frequently employed to model socio-ecological-systems \citep{Kohler.2007, Ostrom.2009, Schlueter.2012}. Examples include a resource economics model of collapse of the Rapa Nui population \citep{Brander.1998} or an exploration of vulnerabilities in forager-resource systems \citep{Freeman.2012}. More recently, \citet{Motesharrei.2014} develop a dynamical systems model to study trajectories to collapse in a coupled hierarchical consumer resource system. Similarly, \citet{Nitzbon.2017} employed a model of coupled carbon, population and capital stocks to identify preconditions of collapse and sustainability.

In the last decade, a new generation of models was conceptualized turning to larger and up to global scales and embracing socio-ecological complexity by coupling concepts from dynamic land-use global models and agent-based models \citep{Rounsevell.2011, Arneth.2014}. Building upon this framework \citet{Brown.2019} identified societal breakdown as an emergent property of large-scale behavioral models of land-use change to different climate-economic scenarios. To strengthen the human side of the equation(s), much attention was payed to increase realism when modelling human decision making in socio-economic and socio-ecological models \citep{Schluter.2017, Schill.2019} as well as modelling collective behaviour and social tipping dynamics \citep{Wiedermann.2020}. Recently, the mentioned efforts precipitated in a comprehensive modelling framework for so-called World-Earth Models \citep{Donges.2020}. 

In this work we explore the theory of the {\em collapse of complex societies} in terms of the mechanisms leading to increasing societal complexity and ultimately to collapse. For this purpose, we implement an agent-based network model consisting of productive and administrative nodes. This implementation successfully recaptures the theory's expected behaviour of collapse due to diminishing marginal returns on complexity. We then go on and study its macroscopic dynamics in detail and identify two escape mechanisms from collapse: increasing the output elasticity of labour, or by allowing social mobility by stochastic transitions of nodes between productive and administrative states. In order to approach the theory of societal collapse from a conceptual angle, we here propose to reduce it to three key elements: 
\begin{enumerate}
    \item The basic currency of any society is energy, since labour and material goods can be viewed as derived from energy.
    \item Problems need to be solved when energy availability is deficient as a consequence of stochastic events or shocks. According to the self-reinforcing process denoted as \textit{energy-complexity spiral} (see above) this process always increases societal complexity.
    \item Gains in complexity are resembled by an increase in administration, since it encapsulates gains in specializations of social roles, hierarchies and control, as well as information flow. Additionally, the size of an administrative body is a very graspable example of complexity.
\end{enumerate} 

The remainder of this paper is organized as follows. We start with a thorough description of our proposed network model in the next section. We first present the basic model and then go on to describe a model variant, which includes a mechanism to counteract collapse. In the results section we first present the stochastic model by discussing exemplary trajectories of the system followed by the derivation of an analytical approximation. Based on this approximation we present a comprehensive ensemble analysis of all crucial model parameters. The paper concludes with a discussion of all relevant results and an outlook to future work.


\section{Model description} 

\subsection{Tainter-inspired network model of a hierarchical society steering into diminishing mar-ginal returns and collapse} 

\begin{figure}[!t]
    \centering
    \includegraphics[width=.65\linewidth]{figures/pub_figure1.png}
    \caption{Sketch of the model setup. Nodes/agents are either attributed as \textit{administrators} $A$ or \textit{labourers} $L$. Labourers that are connected to at least one administrator are denoted as \textit{coordinated labourers} $C$ that produce higher energy outputs due to an increased efficiency. The influence of each administrator is marked with green shading. In case the energy production falls below a critical value, the coordinated labourer that is connected with most other nodes (here marked in orange) becomes an administrator as well.}
    \label{fig:example_network}
\end{figure}

In the following we introduce the model that is used in this work to exemplify the dynamics of collapse that are proposed in~\citet{Tainter.1988}. Its goal is to illustrate the increase in societal complexity as a response to the solving of problems, here the sustaining of
a certain level of energy supply $E$. We define complexity as an increase in the diversification of the society (the creation of a new role) and the corresponding increase in control over the behaviour of others (an increasing administrative body). 
We will show that with only including these two processes, one already observes an increase in an appropriately defined measure of complexity, and, hence, we do not need to specifically account for additional phenomena, such as specialisation or an increase in hierarchies.

In particular, we represent an artificial society by means of a complex network~\citep{newman_structure_2003, castellano_statistical_2009} as complex networks have been successfully used in the past to model social systems with respect to collapse~\citep{heckbert_mayasim:_2013}, sustainability~\citep{wiedermann_macroscopic_2015, barfuss_esd_2017, holstein2020optimization} and their ability to adapt to new conditions~\citep{auer_dynamics_2015, schleussner_clustered_2016}. They are usually comprised of entities of two basic types, i.e., nodes and links. For our case the nodes indicate (representative) small but long-lived entities that supply labor, such as family lines, in a society. Links between the nodes or entities represent some
comparatively stable
social tie, e.g., a professional relationship, between them. $N$ counts the total number of nodes in the network, Fig.~\ref{fig:example_network}.

The sole goal of the society is to produce energy $E$ at a certain level in order to sustain its function. The corresponding per-`capita' (actually, per-node) energy requirement of the society is denoted as $\epsilon$. In order to produce this energy certain nodes (denoted as labourers $L$) harvest an individual energy source $R$. $R$ is subject to external shocks that represent an abstract form of a problem that needs to be solved. In order to increase the energy that is produced from the resource, some labourers $L$ can also be appointed to become administrators $A$. The purpose of the administration is to increase the efficiency of the labourers $L$ that are under its direct management. 
We assume each labourer $L$ who is directly connected with at least one administrator $A$
becomes part of a more efficient body of the labour force which has increasing returns to scale, expressed by an output elasticity of labour of $\phi > 1$ (while the rest of the labour force has an output elasticity of labour of $1$).
We denote this class of more efficient labourers the \textit{coordinated labourers} $C$, Fig.~\ref{fig:example_network}. While they increase the output elasticity of their immediate surrounding, the administrators $A$ themselves do not produce energy any longer (thus indicating a certain cost of maintaining an administrative body). Note that being additionally connected with more than one administrator $A$ does not further increase the output elasticity of the coordinated labourers $C$. Along the lines of ~\citet{Tainter.1988} we measure the resulting societal complexity $S$ in terms of the size of the administration $N_A=|A|$ as it represents the level of behavioral control 
in our model society.

We model the temporal dynamics within our artificial society according to the following rules. At each time step $t$ the available resources per node, $R(t)$, are drawn from a beta distribution with parameters $\alpha = 15$ and $\beta = 1$. The corresponding probability density function is 
\begin{align}
    f(R) = \alpha R(t)^{\alpha-1}. \label{eqn:resource}
\end{align}

The stochasticity of $R(t)$ here represents the occurrence of shocks or problems to be solved by society. Smaller values of $R(t)$ represent larger shocks or problems. We have chosen $f(R(t))$ such that drawing $R(t)$ mostly results in $R(t)$ close to $1$ (i.e., no or small problems) with a rare chance of $R(t) \ll 1$ (i.e., large problems). Depending on the number of labourers $N_L(t)$ (uncoordinated) and $N_C(t)$ (coordinated) the total produced energy $E(t)$ is computed as
\begin{align}
E = R (N_L + N_C^{\phi}).\label{eqn:energy}
\end{align}

If energy per capita falls below a minimum value representing the society's vital needs, i.e., if $E(t)/N < \epsilon$ for some parameter $\epsilon$, we assume the society tries to solve this problem by appointing exactly one (additional) administrator such that the output elasticity of some labourers is increased. The new administrator is selected as follows:
 \begin{itemize}
    \item[(a)] If currently no administrators exist in the society (i.e., if $N_A=0$) the labourer $L$ with the most network connections becomes the sole administrator.
    \item[(b)] If at least one administrator currently exists (i.e., if $N_A>0$), the coordinated labourer who has the most network connections becomes an additional administrator.
\end{itemize}
The model then proceeds to the next time step and the availability of $R(t)$ as well as the produced energy $E(t)$ is again computed according to Eqn.~\eqref{eqn:resource} and  Eqn.~\eqref{eqn:energy}.

Following~\citep{Tainter.2006} we define the potential onset of societal collapse, i.e., the vicinity to a critical state, once the system approaches diminishing returns on complexity $ROC$. We operationalize $ROC$ into our model by computing the differences in energy $E$ and complexity $S$ after the recruitment of one additional administrator $A$ at time $t$, 
\begin{align}
    ROC(t) = \frac{E(t) - E(t-1)}{S(t) - S(t-1)}.
\end{align}

Note that $ROC(t)$ depends on the specific structure of the underlying social network, in particular on who the newly appointed administrator $A$ is connected to and how many labourers turn from $L$ to $C$ as a consequence. As outlined above we measure societal complexity $S$ simply in terms of the size of the administration $N_A$, i.e., $S(t)=N_A(t)$. Since the size of the administration is increased by one due to the recruitment of a single additional administrator at time $t$ we obtain $N_A(t) = N_A(t-1)+1$ such that $S(t) - S(t-1) = 1$. Therefore the return on complexity $ROC$ reduces to
\begin{align}
    ROC(t) = E(t) - E(t-1).
\end{align}

Once $ROC(t)$ becomes negative our artificial society is expected to begin its decline into collapse since the energy production $E$ becomes smaller and approaches zero with increasing time $t$ and a correspondingly ever increasing size of the administration $N_A$. In the following section we therefore directly display results for the energy production $E(t)$ over time and interpret a positive slope as increasing returns to complexity and a negative slope as decreasing returns to complexity.

For the initial model setup we consider an Erd\H{o}s-R\'enyi random network~\citep{erdos_evolution_1960}. It consists of $N$ ($N=400$ for our case) nodes that are all initially assumed to be labourers ($N_L(0)=N$, $N_A(0)=N_C(0)=0$). Additionally, we place a link between each pair of nodes with a fixed probability $\rho$ (which we will vary between $\rho=0$ and $\rho=0.1$ for our analyses). Thus, $\rho$ gives the expected density of links in the resulting network, while $N(N-1)/2$ is the maximum number of possible links. Starting from this setup, we simulate the system for a maximum of $t_{max}=10,000$ time steps following the model logic that is given above. In the case when all nodes becomes administrators (i.e., if $N_A(t)=N$) we stop the simulation even though $t_{max}$ may not have been reached. In that case no more energy is produced ($E=0$ due to the lack of labourers $L$) and we consider the society to be entirely collapsed.

\subsection{Social mobility as a possible countermeasure to collapse}

The basic model that we introduced in the previous section represents closely the assumptions put forward in~\citet{Tainter.1988}. Here a society may get caught in the energy-complexity spiral as the increased or even flat demand for energy increases its complexity (measured as the size of the administrative body $N_A(t)$) when being faced with ever new problems.

One obvious downside of the above model setup is that, once selected, an administrator $A$ does not change its role back to become a labourer $L$ (or $C$) again. This decision is motivated by findings that social mobility is predominantly stable or upward rather than downward within  \citep{Stawarz.2018} and across generations~\citep[e.g.,][]{adermon_intergenerational_2018}. Similarly, it is only possible to become an administrator $A$ if one is selected from the set of coordinated labourers $C$, meaning that a spontaneous jump in hierarchy is permitted within the logic of the model.

However, in order to increase social mobility, we also aim to investigate a version of our model that allows for the nodes to randomly change their state with a fixed small probability. Such dynamics can be interpreted as the loss or gain in socio-economic status within or across generations. Within the logic of our model such an additional process is implemented as follows. At the beginning of each time step $t$ every node changes its state with probability $p_e$ either from $A$ to $L$ or $C$ (depending on whether it is connected to another administrator), or from $L$ or $C$ to $A$. We expect that this allows to regulate the complexity of the society and possibly avoid its unlimited growth (as for this setting the complexity may also decrease as the problems/shocks become smaller again). Hence, such process should cause a sustainable ratio between the size of the administration $N_A$ and the body of labourers $N_L+N_C$. Note that if we set $p_e=0$, we again obtain the dynamics of the model setup that is outlined in the model description.


\section{Results and Discussion}

\subsection{Example trajectories} 

Figure \ref{fig:fig2_exampleTrajectories}A shows the implementation of the original approach to Tainter's theory of societal collapse with dynamics described above in the model description.
Initially the model society benefits from its response to external shocks. The reason is that initially each new administrator $A$ has a much higher number of connections to labourers $L$, which are thus turned into coordinated labourers $C$ producing energy output at increasing returns to scale (here $\phi = 1.05$). As the amount of labourers $L$ decreases, the marginal returns on complexity ($ROC$)  of each new administrator $A$ decreases until it becomes negative. Figure~\ref{fig:fig2_exampleTrajectories}B displays the energy output as a function of administration which resembles well the assumed underlying principle of diminishing marginal returns on complexity in \citet{Tainter.2006}. In our first approximation to the theory of collapse, the model society can only react to shocks by choosing the best possible administrator each time the energy requirements are not met. We interpret the formation of an administrative body as an increase in societal complexity. After a period of increased energy output, marginal returns decrease while correspondingly the share of administration rises, caught in an energy-complexity-spiral until the society collapses (here $E = 0$), Fig.~\ref{fig:fig2_exampleTrajectories}A.

\begin{figure}[!t]
    \centering
    \includegraphics[width=\linewidth]{figures/pub_figure2.png}
    \caption{Exemplary network simulations of a Tainter-like model society ($N = 400$) according to the model description. The panels show only the first 5000 timesteps to focus on the initial dynamic. Blue curves show the share of administration in the network (light blue: Administration as result of decreased resource availability, dark blue: Administration resulting from exploration with rate $p_e$). Orange curves show the average energy produced per node. A) shows the typical development of a network reacting to shocks by changing one node from coordinated labourer $C$ to administrator $A$. B) displays a moving average of the energy measured in A) against the size of the administration. The black curve indicates stylized parabola-shaped diminishing marginal returns as in ~\citep{Tainter.2006}. C) displays the network development at an intermediate exploration probability ($p_e = 0.00275)$. D) shows the development of a network high probability of exploration ($p_e = 0.02$).}
    \label{fig:fig2_exampleTrajectories}
\end{figure}

Next, we study whether small adaptation mechanism of the model society enable the model to overcome the fast collapse shown in Fig.~\ref{fig:fig2_exampleTrajectories}A. For this purpose we allow for a random exploration of node states 
enabling the model agents to change their state ($A \longrightarrow L/C$ or vice versa) with a low probability ($p_e$ = 0.00275). For a network of $N=400$ individuals this amounts to an average of 1.1 nodes changing their status at random per time step. 


As figure \ref{fig:fig2_exampleTrajectories}C shows, an exploration rate of this low magnitude can already be insufficient to overcome the collapse that follows from the event-triggered selection of administrators indicated by a higher survival time of the society, compare Fig.~\ref{fig:fig2_exampleTrajectories}A/C. However, energy output $E$ is still very low at $p_e = 0.00275$, due to the high degree of administration in the network over sustained periods of time.

We also study the case of larger rates of exploration, i.e, $p_e=0.02$, Fig.~\ref{fig:fig2_exampleTrajectories}D. We find an initially sharp increase in administration, facilitated by the high exploration probability. The trajectory then converges to a stable administrator share slightly above $N_A/N=0.5$. This fraction is positively offset by additional effect based promotion of nodes, because energy production per capita is always below the threshold ($\epsilon = 1.0$). After reaching this (meta-)stable state, the model society seems to survive for a long time (at least longer than $t=5000$) while faced with a decreased mean energy production at around $E\approx 0.5$.

It is noteworthy that in each instance of the simulation (Fig.~\ref{fig:fig2_exampleTrajectories}A, C, D) the society reaches an energy production above the energy requirements, which is then however quickly surpassed until the energy level either tends towards zero for zero or low exploration rates in Fig.~\ref{fig:fig2_exampleTrajectories}A,C or converges to a stable regime for large exploration rates in Fig.~\ref{fig:fig2_exampleTrajectories}D.

\subsection{Deterministic macroscopic approximation of the stochastic micro-model} 

We now derive a macroscopic approximation of the above stochastic micro-model in terms of an ordinary differential equation describing the average time evolution of an aggregate quantity (here the total number $N_A$ of administrators). This time evolution is governed by average transition rates between the three groups. The rates are based on approximations of the probabilities with which individual agents switch between the groups and the assumption that $N$ is large so that the law of large numbers applies.

There are two processes which make the number of administrators change, exploration and targeted recruiting in response to energy demands 
Due to exploration, from the $N_A$ many administrators, on average $p_e N_A$ many leave the administration per unit time, and from the $N - N_A$ many non-administrators, on average $p_e (N - N_A)$ many are hired as administrators per unit time, making a balance of $p_e (N - 2 N_A)$ additional administrators on average per unit time.

For targeted recruiting, we make the simplifying assumption that the density of links inside and between the three groups remains approximately constant and can thus be estimated by the overall link density $\rho$ of the Erd\"os-Renyi network. Then, if $N_A$ is the number of administrators, the probability that a non-administrator is not linked to any of the $N_A$ administrators (and is thus un-coordinated) is approximately
$(1-\rho)^{N_A}$, hence the numbers of un-coordinated and coordinated labourers are approximately
$L \approx (N - N_A) (1-\rho)^{N_A}$ and
$C \approx (N - N_A) (1 - (1-\rho)^{N_A})$,
hence we have
$E = R e$ with $e = L + C^\phi \approx (N - N_A) (1-\rho)^{N_A} + [(N - N_A) (1 - (1-\rho)^{N_A})]^\phi$.

An additional administrator is recruited iff $E/N < \epsilon$, i.e., iff $R < N \epsilon / e$. The probability that $E/N < \epsilon$ is then $P = F(\epsilon / e, \beta, \alpha)$, where $F(\cdot, \beta, \alpha)$ is the cumulative distribution function of the beta distribution.
The expected number of administrators hired additionally due to shocks per time unit is then also equal to this probability $P$.

In all, we get the approximation
\begin{align*}\label{eqn:approx}
    \frac{dN_A}{dt} &\approx p_e(N - 2 N_A) + F\left(\frac{N \epsilon}{(N - N_A) (1 - \rho)^{N_A} + [(N - N_A) (1 - (1 - \rho)^{N_A})]^\phi}, \beta, \alpha\right), \\
    C &\approx (N - N_A) (1 - (1 - \rho)^{N_A}), \\
    L &\approx (N - N_A) (1 - \rho)^{N_A}
\end{align*}
as long as $N_A<N$.

\begin{figure}[!t]
    \centering
    \includegraphics[width = \linewidth]{figures/pub_figure3.png}
    \caption{Analytic approximation of the model dynamics with the share of administration $N_A/N$ (A) and corresponding energy production per capita (B). Transparent solid lines show the results of an ensemble of 100 simulations of the microscopic network model for different values of $p_e$. Dashed lines display the respective analytical approximation. (C) Histograms of collapse frequencies of the model, i.e., frequency of times at which energy production approaches zeros, $E=0$. Bars located at time $t \geq 5000$ show right-censored cases of the simulation.}
    \label{fig:fig3_macro_approx}
\end{figure}

Simulation results are well predicted by the approximation, Fig.~\ref{fig:fig3_macro_approx}A,B. Particularly at higher values of $p_e$, the approximation yields excellent results both for the energy production and the administration share. Only at exploration probability of $p_e=0$ (dashed line and light colors in Fig.~\ref{fig:fig3_macro_approx}), the macroscopic approximation slightly overestimates the velocity of collapse. This can be explained by one simplifying assumption of the macroscopic approximation, namely the network degree or local link density to be homogeneous, i.e., approach $\rho$, across the entire network. In contrast, due to the links being put at random when setting up the network, this distribution in a node's number of connections becomes heterogeneous with some nodes showing an above-average number of connections. Because those nodes are preferred over those with few connections in the process of administration selection, initially selected administrators will result in an above average energy return. This produces the prolonged slow-growth phase of the administration in the numerical simulation. In contrast, the macroscopic approximation assumes that at a given point in time any new administrator would have the same average effect. 

To illustrate the probability of collapse under different exploration assumptions, histograms of survival times from the stochastic model are computed, Fig~\ref{fig:fig3_macro_approx}C. At high rates of social mobility $p_e$ the model converges to a stable state where the society persists for infinite times due to its ability absorb shocks. For lower $p_e \approx 0$ we also find isolated cases of infinite survival even though the society mostly collapses comparatively early, Fig.~\ref{fig:fig3_macro_approx}C. Such observations can be explained by a fragmentation of the network into smaller worker networks that are mutually disconnected. The scattered survival times along the time axis for intermediate exploration rates illustrate stochastic survival of the model society due to a random sequences of shocks. 

\subsection{Influence of model parameters on survival time} 

In the third step of our analysis we use the approximation proposed in the previous section to estimate model outcomes for a broad range of parameter values $\rho$, $\phi$ and $p_e$. Note, that the additional parameters of resource availability ($\alpha$, $\beta$) and threshold ($\epsilon$) have a major influence on the outcome of the model as well. Specifically, a high probability of low resource availability leads to much shorter survival times and vice versa while a low threshold $\epsilon$ to appoint further administrators considerably increases survival times. However, since these effects did not reveal any additionally remarkable results (not shown), we focus our attention on the major drivers of survival time and energy output, i.e, the exploration probability $p_e$, link density $\rho$ and output elasticity of labour $\phi$. 

\begin{figure}[!t]
    \centering
    \includegraphics[width = \linewidth]{figures/pub_figure4.png}
    \caption{Survival time analysis of the parameters exploration probability ($p_e$), link density ($\rho$) and efficiency ($\phi$).
    (A) The effect of $p_e$ between zero and $\frac{1}{N}$ for some select values of $\phi$ on median survival time over all computed values of $\rho$. The vertical dashed line defines the exploration threshold of $\frac{1}{N}$.
    Lower panels ($B_{1-4}$) display the relation between $\phi$ and $\rho$ for different exploration probabilities. Grey areas in the lower panels indicate a survival time $\geq 10000$, potentially going to infinity.}
    \label{fig:parameters}
\end{figure}

Panels \ref{fig:parameters}$\text{B}_{1-4}$ show the survival time as a function of link density $\rho$ and output elasticity of labour $\phi$. Figure~\ref{fig:parameters}$\text{B}_1$ displays a parameter analysis for the case of no exploration, $p_e=0$ and can hence be interpreted as showing minimum requirements of a successful administration. The most important requirement being a sufficiently large connectivity of the social network (indicated by the increasingly grey area with increasing $\rho$). Also, the extreme case of $\rho = 0.0$, i.e, a network with no connections, must be noted since in this case nodes cannot be converted  at all to administrators and thus the model remains in its initial state. As $p_e$ approaches $\frac{1}{N}$  (Fig.~\ref{fig:parameters}$\text{B}_{2-3}$) the relevance of $\rho$ diminishes and the survival time becomes mainly a function of the output elasticity $\phi$. 

Figure \ref{fig:parameters}A demonstrates the effect of exploration for selected values of $\phi$. 
We have chosen a range of values for $\phi$ that roughly correspond to the order of magnitude of returns to scale reported in
\citet{balevzentis2015trends,diewert2011returns}.
For exploration probabilities $p_e<1/N$, i.e., when less than one individual per time unit switches their type, collapse is likely, which is indicated by reduced median survival times, Fig.~\ref{fig:parameters}A. 
In fact, our macroscopic approximation shows that $N_A=N$ is a stable fixed point of the approximate dynamics for $p_e < 1/N$ since the RHS of Eqn.\ \ref{eqn:approx} remains strictly positive for $N_A\to N$. This means that the labour force will vanish in finite time, thereby also resulting in a finite survival time.
Additionally, in this ``low exploration regime'', a larger output elasticity of labour $\phi$ tends to increase the median survival time. 
Indeed, the term $F(\cdot)$ in Eqn.\ \ref{eqn:approx} is a decreasing function of $\phi$ for large $N_A$, which means a larger $\phi$ slows down the collapse towards $N_A=N$.
The effect of a $p_e$ on survival time is more difficult to understand because it is twofold, as can again be seen in Eqn.\ \ref{eqn:approx}. Initially, when still $N_A<N/2$, larger $p_e$ increases $dN_A/dt$ and speeds up the process of recruiting administrators, as seen on the left of Fig.~\ref{fig:fig3_macro_approx}A. As soon as $N_A>N/2$, the reverse happens and a larger $p_e$ slows down that process. For small $p_e$ below some turning point value $p_e^0<1/N$, the first effect dominates, so that survival time decreases. For $p_e$ between $p_e^0$ and $1/N$, however, the second effect dominates, so that survival time increases again.
For some values of $\phi$, the two effects are clearly distinguishable in different phases of the evolution, as can be seen in the realization depicted in Fig.~\ref{fig:two_phase_admin_growth}, where the system stays a long time close to $N_A=N/2$ before reaching $N_A=A$ and thereby collapses. 
As $p_e$ approaches the critical value of $1/N$, this effect becomes ever stronger and expected survival times approach infinity.

\begin{figure}[!t]
    \centering
    \includegraphics[width = \linewidth]{figures/pub_figure3_b.png}
    \caption{Trajectories of administration share for $\phi = 1.18$ and $\rho = 0.2$ for two choices of $p_e$ below the critical value $p_e=\frac{1}{N}$ above which collapse of the system is mitigated. The trajectories show how the administration share is stabilized by exploration. This mechanism effectively slows down the total change in administration size until a critical administration size is reached at which administration growth enters the event triggered phase.}
    \label{fig:two_phase_admin_growth}
\end{figure}

Of course, since the microscopic network model is of stochastic nature, individual survival times still vary quite a bit between realizations for the same set of parameters. As survival time also depends on other parameters, we display only the central tendency of survival time in Fig.~\ref{fig:fig3_macro_approx}A. Since we terminate our simulations after $t=10000$ time steps, we cannot distinguish higher survival times, which means we cannot use the arithmetic mean of survival times as the measure of central tendency. This is why we display the median instead, which has the additional advantage of being a more robust statistic than the mean.

 
For larger values of $\phi$, we see that some realizations survive for very long even though $p_e<1/N$. This is because for large enough $\phi$, the approximate dynamics Eqn.\ \ref{eqn:approx} have a second stable fixed point despite $N_A=N$. Indeed, for $\phi\to\infty$, the $F(\cdot)$ term vanishes and $dN_A/dt=0$ for some value close to $N_A=N/2$. In that case, realizations of the stochastic model starting with small $N_A$ are likely to stay close to $N_A=N/2$ for very long (which is thus a ``metastable'' state) before eventually escaping into the basin of attraction of the ``collapse'' fixed point $N_A=N$ due to a large enough shock. 

When $p_e$ crosses the critical value $1/N$, the collapse fixed point $N_A=N$ of the approximate dynamics becomes unstable, so that collapse in the can only occur due to a sequence of large shocks which becomes increasingly unlikely. The median survival time is then very large.
For large $p_e$, one can also see from Eqn.\ \ref{eqn:approx} at which levels of $N_A$ occur most likely over the course of the simulation. Since the $F(\dots)$ term is between 0 and 1, $dN_A/dt=0$ implies $N/2 \le  N_A \le N/2 + 1/2p$. 
In other words, when either $p_e$ or $\phi$ are large enough, one can expect that after some transient phase, there will be only slightly more administrators than labourers on average for a long time, before eventually a large shock eventually causes the system collapse after all. 

\section{Conclusion and Outlook}

We have developed an agent-based network model to illustrate emergent dynamics from the {\em theory of societal collapse} as postulated in~\citet{Tainter.1988}. In our model, problem solving increases complexity and ultimately leads to collapse of the system by diminishing marginal returns. Problems were abstractly represented by random shocks in energy production, which were counteracted by adding complexity to the networked social system in form of an increasing administration. We were able to derive a well-performing macroscopic approximation of the proposed model, thus providing a very simple mathematical description of the theory in form of an ordinary differential equation. 
Additionally, we carried out an analysis of the effect of random status exploration where agents can change their status from a productive to an administrative state (and vice versa). We found that a small exploration rate of at least $p_e=1/N$ is already sufficient to counterbalance the progression towards a network consisting only of administrators under the influence of shocks and reach a stable state where survival times are very long. However, smaller exploration rates have adverse effects on survival time, as they only facilitate the spread of administration over the network. 

In addition, we find that a minimal link density $\rho$ in the network is essential for long survival times, but has no further effect beyond a threshold value. Also, large output elasticity of labour $\phi$, i.e., the effectiveness of administrators to increase productivity of labour, can increase the survival time. When output elasticity of labour is large enough, networks may survive for very long until a large enough shock destroys them, regardless of the exploration rate. Finally, we note that with increasing time $t$ and correspondingly increasing share of administration in the network, the potential to change and react to external influences was reduced because marginal returns on complexity diminished. By a similar argument \citet{Allison.2004} show how, for instance, the Western Australian Agricultural Region could end up in a `Lock-in Trap', which is characterized by low ability to adapt to changing external circumstances. 

In summary, within the scope of our model two main strategies allow mitigation of collapse: (i) an increase in the output elasticity of labour or (ii) a random transition of people between productive and administrative roles at a sufficiently large rate.

With regard to open research topics posed by \citet{Cumming.2017}, such as moving the debate from `whether' to `why' collapse occurred and possible strategies to ``avoid, slow or hasten collapse'', we wish to contribute our findings to the discourse of the drivers of socio-ecological resilience and collapse. Growing complexity has recently been identified as a source of increased societal risk due to, e.g., high degrees of socio(-ecological) interconnections that can lead to cascading failures \citep{Helbing.2013, Rocha.2018}. In this regard complexity is considered mainly a facilitator of collapse by shaping systems to become more vulnerable. In our research we show that complexity may be even more directly linked to collapse since societies may reach a point where associated costs outweigh the benefits, thereby directly causing their breakdown. This dimension is, to our understanding, an underrepresented view of socio-ecological collapse. 

Similar to \citet{Motesharrei.2014} we find that an unequal society in which the elite consumes more (in their case) or produces less (in our case) than the commoners is likely to run into collapse. In addition, we found that the absence of an elite (i.e.\ egalitarian society), which is comparable to the random state exploration mechanism in our model, produced more pathways to avoid collapse than a hierarchical (complex) society. In comparison, previously reported simulations of competing pre-modern societies showed that collapse fades once societies become complex and can develop intensive agriculture at sufficiently large scale \citep{Guzman.2018}. This finding contrasts our result for two specific reasons: First, in \citet{Guzman.2018} complexity is a binary state and cannot increase further once a society is denoted as being {\em complex}. Second and more importantly, \citet{Guzman.2018} model a warrior class comparable to administrators in our model, which is directly coupled to available food (energy). Hence complexity in terms of hierarchies cannot exceed the level of existing energy and thereby grow to an unsustainable size. A more detailed investigation into the coupling of complexity and economic development in future work would certainly help in understanding the importance of complexity in collapse.

\citet{Butzer.2012} points out that collapse itself should be modelled by sophisticated socio-ecological models, which do not draw on simplifying assumptions. Nevertheless, \citeauthor{Butzer.2012} concludes that preconditions for breakdown are typically economic or fiscal decline. We agree and emphasize that the our model of the theory of {\em collapse of complex societies} offers a potential path leading to exactly such preconditions.

Like in \citet{Watson.1983} or \citet{Granovetter.1978}, the objective of our research was not to predict any guidelines for resilient societies, but to illustrate the mechanism of collapse due to diminishing marginal returns of an increasingly complex society. We acknowledge that our approach to modelling collapse is highly conceptual, nevertheless our results underline the relevance of understanding the emergence and consequences of complexity in the resilience discourse. 

In contrast to more economically spirited models, we do not include in our simple model any form of competition, economic externality, or resource depletion, which in the real world can be alternative or additional drivers of collapse. An additional important limitation of our model is that we incorporated no mechanism for reduction in bureaucracy beyond random state exploration in contrast to the automated increase in complexity due to shocks. However, one major argument in ~\citet{Tainter.1988} is that endogenous reduction of complexity as problem solution strategy is not observed in historical cases, at least it is not a dominant process. Instead we assume that the model's administration possesses characteristics of persistent power structures \citep{Acemoglu.2008}, which come into existence much more easily than they are removed. Further we compare the administration to expanding elite bureaucracies according to \citet{Parkinson.19.11.1955} and \citet{Weber.1922}. As an extension, it would be interesting to consider more complex dynamics of intra-generational and inter-generational mobility of social status and wealth and their respective change over time. 

Finally, we acknowledge that our proposed model only represents one feature of complexity, i.e., the emergence of one additional hierarchy layer, the administration. Of course, modelling efforts exist which consider much more intricate mechanisms, like human behaviour, shared resource systems, interacting societies, etc. \citep{Arneth.2014, Rounsevell.2011, Schill.2019, Schluter.2017, Nitzbon.2017, Motesharrei.2014}. But keeping Occam's razor in mind, we believe that this very simple model (without population increase, additional hierarchy layers, dynamic networks, etc.) is sufficient to make a case for the direct influence of complexity with respect to social-ecological collapse. 

Our model should be extended in further studies to analyze how complexity relates to collapse in potentially more realistic scenarios that could, e.g., be represented by a further social stratification or and more complex network topologies. In this context, complexity costs of social connectivity should be considered. Additionally, impacts of innovations and energy substitutes (coal, oil, etc.) on collapse could be studied to identify additional strategies to slow down or avoid collapse. Furthermore, the interrelation between energy and complexity could also be studied in more elaborate World-Earth and Human-Earth system models to study under which circumstances the results obtained in this work could be transferred to less abstract applications. For this purpose, appropriate measures of energy and complexity need to be defined depending on the specific model and research question at hand.

While increasing complexity has brought vast benefits to modern industrialized societies, our findings also raise the question to which degree resilient societies can increase complexity without creating large-scale risks of collapse due to an ever increasing demand for energy. Along those lines, it is of great interest to further investigate whether the rise of complexity always comes with an increased risk or whether certain forms of complexity exist that do not necessarily increase the risk of societies for collapse. Such an assessment would allow to develop design principles for resilient infrastructures and social structures that do not facilitate complexity in the understanding of~\citet{Tainter.1988} thereby allowing modern societies to escape the energy-complexity spiral to path their the trajectory of the socio-ecological Earth System toward a sustainable trajectory.


\section*{Acknowledgements}

The idea for this work was developed during a {\em Studienkolleg} working group on {\em The Earth as a complex system} funded by the German National Academic Foundation (Studienstiftung des deutschen Volkes). The work was carried out in the scope of the COPAN collaboration at the Potsdam Institute for Climate Impact Research. M.W.\, J.H.\ and J.F.D.\ were supported by the Leibniz Association (project DominoES) and the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (ERC grant agreement No.\ 743080 ERA). F.S. was supported by the German National Academic Foundation. The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for providing resources on the high-performance computer system at PIK. For fruitful discussions we sincerely thank Jürgen Kurths, Lea Tamberg, Adrian Lison, Sascha Haupt, Vivian Grudde, Bastian Ott and other members of the Studienkolleg's working group 4. We thank Wolfgang Lucht for inspiring discussions that triggered our interest in Tainter's theory.

\section*{Appendix A: Optimal number of administrators}

From the analytical approximation, we know that the expected energy input $E$ is proportional to
\begin{align}
    E_0 &:= (N - N_A) (1 - \rho)^{N_A} + [(N - N_A) (1 - (1 - \rho)^{N_A})]^\phi,
\end{align}
which has a maximum in $N_A$ where
\begin{align}
    0 = \partial_{N_A} E_0 &=
    - [1 - (N - N_A) \ln(1 - \rho)](1 - \rho)^{N_A} \\
    &+ \phi[(N - N_A) (1 - (1 - \rho)^{N_A})]^{\phi - 1}
    [(1 - \rho)^{N_A} [1 - (N - N_A) \ln(1 - \rho)] - 1].
\end{align}
E.g., for $N=400$, $\rho=0.05$, and $\phi=1.05$, the optimal number of administrators is $N_A\approx 34$ or $\approx 9\%$, leading to an $E\approx 460$ or $E/N\approx 1.15$.

\section*{Model Documentation}

The model code is available at CoMSES OpenABM \href{https://www.comses.net/codebase-release/e03dc306-e6c5-4c8d-bccb-34ce60a7b0fc/}{https://www.comses.net/codebase-release/e03dc306-e6c5-4c8d-bccb-34ce60a7b0fc/} version 1.0.0
 
