# Modelling co-evolution of resource feedback and social network dynamics in
human-environmental systems.

Meghdad Saeedian  Department of Mathematics, Imperial College London, London,
SW7 2AZ, United Kingdom Chengyi Tu  School of Economics and Management,
Zhejiang Sci-Tech University; Hangzhou, 310018, China Fabio Menegazzo
Dipartimento di Fisica “G. Galilei”, Università di Padova, Via Marzolo 8,
35131 Padova, Italy Paolo D’Odorico  Department of Environmental Science,
Policy, and Management, University of California, Berkeley, CA, USA Sandro
Azaele  Dipartimento di Fisica “G. Galilei”, Università di Padova, Via Marzolo
8, 35131 Padova, Italy Samir Suweis  Dipartimento di Fisica “G. Galilei”,
Università di Padova, Via Marzolo 8, 35131 Padova, Italy

###### Abstract

Games with environmental feedback have become a crucial area of study across
various scientific domains, modelling the dynamic interplay between human
decisions and environmental changes, and highlighting the consequences of our
choices on natural resources and biodiversity. In this work, we propose a co-
evolutionary model for human-environment systems that incorporates the effects
of knowledge feedback and social interaction on the sustainability of common
pool resources. The model represents consumers as agents who adjust their
resource extraction based on the resource’s state. These agents are connected
through social networks, where links symbolize either affinity or aversion
among them. The interplay between social dynamics and resource dynamics is
explored, with the system’s evolution analyzed across various network
topologies and initial conditions. We find that knowledge feedback can
independently sustain common pool resources. However, the impact of social
interactions on sustainability is dual-faceted: it can either support or
impede sustainability, influenced by the network’s connectivity and
heterogeneity. A notable finding is the identification of a critical network
mean degree, beyond which a depletion/repletion transition parallels an
absorbing/active state transition in social dynamics, i.e., individual agents
and their connections are/are not prone to being frozen in their social
states. Furthermore, the study examines the evolution of the social network,
revealing the emergence of two polarized groups where agents within each
community have the same affinity. Comparative analyses using Monte-Carlo
simulations and rate equations are employed, along with analytical arguments,
to reinforce the study’s findings. The model successfully captures key aspects
of the human-environment system, offering valuable insights to understand how
both the spread of information and social dynamics may impact the
sustainability of common pool resources.

##  1 Introduction

The sustainable use of natural resources that support the metabolism of human
societies is a pressing global challenge [1, 2, 3]. The governance of natural
resources that are used in common by multiple individuals (‘common pool
resources’, or CPRs) is increasingly recognized as a fundamental challenge to
ensure the sustainability of the planet. In the case of open access systems in
which resource harvesting is not regulated, users are (relatively) free to
over-exploit the CPR [4, 5, 6], a phenomenon known as the “Tragedy of the
Commons” [7, 8]. In this case, self-interested, short-term, profit-maximizing
actors, harvest the resource in disregard of its long-term sustainability.
This behavior leads to resource exhaustion, leaving future generations empty-
handed [9, 10]. However, in many cases, CPRs are not open access but the
common property of a community; access is restricted to members of the
community governing the resources according to internal rules and, in many
cases, the sharing of common goals prevents over-exploitation [4, 10].

From a mathematical perspective, a complex and dynamic interaction between
human behavior and CPRs can be modeled as a human-environment system (HES)
[11, 12, 13].

In modeling the HES dynamics, one can simplify the problem by considering only
the interactions between consumers and resources, neglecting the social
structure or resource heterogeneity. In this framework, consumers are treated
as players in a game whose payoffs depend on their extraction rates. Consumer
strategies evolve according to game-theoretic principles, based on the rewards
or costs of their choices relative to other players (or with respect to the
average payoff) [14]. It is well known that in non-repeated games, under the
assumption of rational and purely self-interested agents, players have an
incentive to defect, and defection is the only (Nash) equilibrium. This is a
mathematical analog of the tragedy of the commons [14].

However, to understand the evolution of the HES, it is necessary to consider
not only the dynamics of the consumers [15, 16, 17], but also the dynamics of
the resource [18, 10]. This can be achieved by using one or more differential
equations that describe how the resource changes over time, taking into
account natural growth, trade, and consumption [10]. Moreover, the dynamics of
consumers and of resources are intrinsically coupled. In fact, there exists a
cyclical and responsive relationship between human actions and ecological
changes. Thus, we need a modeling framework in which players’ decisions
directly influence the environmental system, and in turn, the evolving state
of the environment affects future choices and outcomes in the game [19, 20,
21, 10]. Furthermore, in this context, the structure of the social network of
consumers [22] and stochastic fluctuations in resource dynamics [18] may also
have an important effect on the sustainability of CPRs. However, such elements
have not yet been considered in a unified theoretical framework.

In this work, we fill this gap by introducing in the framework of games with
environmental feedback three novel elements: 1) users’ knowledge feedback,
i.e. the behaviour of the player is directly affected by the state of the
resources; 2) The consumer social network structure and dynamics, and in
particular the effect of polarization in the CPR sustainability; 3)The effect
of stochasticity in the resource evolution.

We term knowledge feedback the agents’ behavior of adjusting their extraction
levels according to the resource state. This means, as documented and studied
[23, 24, 19], that consumers defect when the resource is abundant and
cooperate when it is scarce (Fig. 1 (A)).This is rationalized as follows: the
classic replicator equation determines the strategy of a player based on the
average fitness of the other agents, regardless of the level of resources
available. However, we embed the very optimistic assumption that humans will
respond to scarce resources by cooperating more.

In our model, we also introduce ”social relations” among consumers. These
relations are characterized as either positive/negative links, symbolizing
various types of relationships such as friendship or animosity, trust or
distrust, collaboration or intense rivalry, political alignment or opposition,
etc. The intensity of the ”social interactions” is quantified by assessing the
satisfaction level between two interacting individuals, as discussed in [25,
26]. It is important to note that the dynamics of ”social interaction” are not
solely dependent on these social links; the individual states of the agents,
like their opinions or strategies, also play a crucial role. For example, a
friendship between individuals with conflicting opinions may not endure,
leading to dissatisfaction: either they find common ground and align their
views, or their relationship deteriorates. Hence, the ”social interaction”
state is a result of both the individual states of the agents and the nature
of their connection (referred to as pair-connection). In our framework,
enduring or transient social interactions are represented as  positive /
negative interactions, respectively. These interactions are also understood as
satisfying or unsatisfying [25, 26], and in this work, we use these two terms
interchangeably.

One of the main results of our work is that knowledge feedback alone can
sustain CPRs in the mean-field case, where consumers ignore the social
structure and each agent interacts with all others. In other words, we find
that long-term exploitability is achievable if consumers align their behavior
with the resource level. However, when social dynamics are also taken into
account, they tend to split consumers into two groups sharing the same
behavior, leading to a depletion/repletion transition of the CPR, where the
stationary resource level and cooperative fraction among agents depend on the
average network degree. In particular, we find that hubs in the consumer
network favor defectors and cause the depletion phase.

Figure 1: Knowledge Feedback: consumers’ adaptation to environmental
fluctuations. (A) Cooperators (white) are more likely to change their behavior
in resource-rich environments, while defectors (blue) do so in resource-poor
environments. (B) The consumers with different strategies are shown
interacting in binary states: positive (cooperators) and negative (defectors).
The ”social interactions” can have two types of ”social link” states: positive
(solid) or negative (dashed). The combination of two nodes and connecting link
is termed pair-connection, and the sign of each pair-connection is determined
by multiplying the sign of constitute components. Three negative pair-
connection on the left tend to become positive pair-connection on the right to
increase the consumers’ satisfaction.

##  2 Method

###  2.1 Resource Dynamics

We model the CPR dynamics using a logistic function that also considers the
extraction effort of defectors and cooperators as

$\displaystyle\frac{dR}{dt}=TR(1-\frac{R}{K})-R\Big{(}x\hat{e}_{C}+(1-x)\hat{e}_{D}\Big{)},\quad
0\leq\hat{e}_{C}<1\leq\hat{e}_{D}$ (1)

where $T>0$ is the natural growth rate, $K>0$ is the carrying capacity and $x$
is the fraction of cooperators in the systems (that evolve in time, see
below). The extraction rates of a cooperator and a defector are denoted by
$e_{C}$ and $e_{D}$, respectively. If $N$ is the number of agents in the
systems, we define $\hat{e}_{C}=N\cdot e_{C}$, $\hat{e}_{D}=N\cdot e_{D}$ and
assume $\hat{e}_{C}\leq\hat{e}_{D}$, reflecting the fact that cooperators
extract fewer resources than defectors [10].

It is possible to gain a qualitative understanding of the resource dynamics by
rewriting Eq. (1) in terms of a more familiar form of logistic function

$\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle
T^{\prime}R\Big{(}1-\frac{R}{K^{\prime}}\Big{)},$ (2)

where the growth rate, $T^{\prime}$, and carrying capacity, $K^{\prime}$, both
depend on $x(t)$, i.e.,

$\displaystyle
T^{\prime}=\Big{(}1-\frac{x\hat{e}_{C}+(1-x)\hat{e}_{D}}{T}\Big{)}T,\quad
K^{\prime}=\Big{(}1-\frac{x\hat{e}_{C}+(1-x)\hat{e}_{D}}{T}\Big{)}K.$

Depending on $x(t)$ ($0\leq x(t)\leq 1$), the time-dependent growth rate
$T^{\prime}$ can be positive, zero, or negative, but it is always bounded
between an upper and lower limit,
$(T-\hat{e}_{D})<T^{\prime}(t)<(T-\hat{e}_{C})$. For everlasting positive
growth rates, $T^{\prime}(0<t<\infty)>0$, the evolution of the resource
follows a logistic dynamics, while for stable negative growth rates,
$T^{\prime}(0<t<\infty)<0$, the trajectory of $R$ approaches zero. Moreover,
possible scenarios involve oscillations between positive and negative growth
rates, leading to oscillations of $R(t)$ between the aforementioned dynamical
states. Therefore, if $R$ (at stationarity) approaches zero, then the CPR is
in a depleted phase; otherwise it is in a repleted phase.

The rescaled carrying capacity is bound by upper and lower limits,
$(1-\frac{\hat{e}_{D}}{T})K<K^{\prime}(t)<(1-\frac{\hat{e}_{C}}{T})K$. Also,
we introduce a rescaled level of available resources $R^{\prime}$, so we can
interpret it as a probability, i.e.,

$\displaystyle R’=\frac{R}{(1-\frac{\hat{e}_{C}}{T})K},\quad\longrightarrow
0<R’\leq 1$ (3)

Additionally, for simplicity, we set $T=K=1$.

###  2.2 Consumers Strategy Dynamics

The consumers’ strategy dynamics, i.e. how the fraction of cooperators $x$
evolves in time, is governed by the update rules graphically presented in Fig.
2.

We investigate three different scenarios of increasing complexity,
corresponding to different game update rules driving the evolution of the
users’ social networks and of the HES.

First, we study the stochastic dynamics of the HES when only knowledge
feedback is considered and no social structure (NS) is included, as depicted
in Fig. 2 (A).

According to this model, the consumers’ strategies depend only on the
availability of resources: they tend to cooperate when resources are scarce
and to defect when resources are plentiful. To implement this model, we make
two fundamental assumptions: 1) the probability of switching from cooperation
to defection is linearly proportional to the normalized resource level (
$R^{\prime}\in[0,1]$ ); 2) The probability of switching from defection to
cooperation is linearly proportional to the complementary normalized resource
level ( $1-R^{\prime}$ ).

Figure 2: Three models that capture different complex human-environmental
dynamics. The white (blue) circles represent cooperator (defector) and solid
(dash) links represent positive (negative) connections. (A) A unified update
rule is derived from the node dynamics of knowledge feedback without
interaction (NS). (B) A unified update rule is derived from the interplay
between the knowledge feedback and the quenched binary-state interactions,
which we call (QBS). (C) A unified update rule results from the interplay
between the node dynamics of knowledge feedback and the evolving binary-state
social links, which we call (EBS). The pair-connections $a$, $c$ and $e$
denote negative social interactions, while the pairs $b$, $d$ and $f$ indicate
positive social interactions. We illustrate here the social dynamical rules
that transform nanegative interactions into positive interactions through node
or link updates. $R^{\prime}$ is the normalized resource level. In this work,
we have implemented a sequential update scheme.

Then we focus on the essential role of complex social interactions [27]
between consumers. We assume that social connections (link) can be either
positive or negative and that they can be quenched (QBS) or evolve in time
(EBS) [25, 26] as shown in Figs. 2 (B)-(C). Regardless of state of nodes, the
positive/negative connection (link) can be interpreted in different ways,
friendship or animosity, trust or distrust, collaboration or intense rivalry,
political alignment or opposition, etc., depending on the context and the
nature of the problem. We note that many complex systems exhibit binary
connections, similar to attraction/repulsion in physics, friendship/enmity in
opinion formation [25, 26, 28], mutualism/competition in ecology [29, 30], or
excitatory/inhibitory in neuroscience [31].

For the QBS and EBS case, we propose a dynamical update rule based on the
interaction between binary state nodes (defectors / cooperators) and binary
state links (positive / negative), as illustrated in Fig. 2 (B-C). The
negative pair-connections $\\{a,c,e\\}$ correspond to short-lasting social
interactions and tend to evolve into long-lasting ones $\\{b,d,f\\}$ to
maximize the level of satisfaction irreversibly. This condition breaks the
detailed balance condition and drives the dynamics out of equilibrium.

To identify positive / negative pair connections, we can adopt the energy
function approach for signed social networks [32, 33, 34, 35]. In this
approach, the unbiased binary variable ($\pm$) represents the node and link’s
state. The main principle behind this approach is that the negative pair-
connections in social systems are characterized by an odd number of negative
nodes and links, whereas the positive pair-connections are characterized by an
even number of negative nodes and links. To be more precise, the sign
(positive or negative) of a pair-connection is determined by multiplying the
signs of the three components involved in that pair-connection. For example, a
configuration of two defectors with a negative link constitutes a
configuration with three negative components, which results in a negative
pair-connection $c$. These three-component interactions are known as higher-
order interactions.

We present in Appendix. A the interpretations and the respective roles of each
of the six possible pair-connections $\\{a,b,c,d,e,f\\}$ in the update rule,
where pair-connections $\\{a,c,e\\}$ are positive (socially satisfying) and
$\\{b,d,f\\}$ are negative (socially unsatisfying) interactions.

###  2.3 Coupled HES Dynamics

We thus propose a co-evolutionary model for HES that incorporates two
interdependent dynamics of human and environment sub-systems, where the
density of cooperators in the human sub-system $x$ and the level of resource
in the environment sub-system $R$ co-evolve.

Putting all the pieces together, the HES dynamics can be modelled by an
undirected network of $N$ individuals, where the nodes represent the consumers
with binary state extraction strategies (cooperator/defector) and the links
represent the binary state social connections. The density of pair-connections
in the network is given by $\rho_{i}=\frac{L_{i}}{L},i\in\\{a,b,c,d,e,f\\}$,
where $L_{i}$ is the number of pair-connections $i$ and $L$ is the total
number of links in the network.

To investigate these co-evolving dynamics we use two methods: microscopic
Monte-Carlo simulations and macroscopic rate equations. We explain these
methods in detail in the following subsections.

###  2.4 Monte-Carlo method

We use Monte-Carlo simulations to study the co-evolutionary dynamics of update
rules (Fig. 2) and resource dynamics (Eq. 1) on complex networks, mainly
Erdős-Rényi networks with size $N$ and mean degree $\mu$, unless otherwise
stated. A Monte Carlo step consists of a sequence of events, each involving
one constituent of the system. The constituents are either nodes (in the case
of knowledge feedback without interaction, Fig. 2 (A)) or pair-connections (in
the case of knowledge feedback with interaction, Fig. 2 (B) and (C)). Thus,
one computer iteration corresponds to $\frac{1}{N}$ or $\frac{1}{L}$ Monte-
Carlo time steps, depending on the case.

  * •

NS model: We use a Monte-Carlo method to simulate the dynamics shown in Fig. 2
(A), where at each iteration, a randomly selected node updates its state
according to a randomly generated number $\xi$ between $0$ and $1$. The node
becomes a cooperator if $R^{\prime}<\xi$, or a defector otherwise. The initial
strategy of the node is irrelevant.

  * •

QBS and EBS models: We use a Monte-Carlo method to implement the dynamics
shown in Fig. 2 (B) and (C), where at each iteration, a randomly selected
pair-connection updates its state according to the corresponding update rules
if it belongs to one of the negative interactions $\\{a,c,e\\}$. Otherwise, no
update occurs. For example, in the QBS dynamics (Fig. 2(B)), if a pair-
connection $e$ is chosen and a random number $\xi$ is generated, the
cooperator switches to defection if $R^{\prime}>\xi$, or the defector switches
to cooperation otherwise. However, in the EBS dynamics (Fig. 2(C)), if a pair-
connection $e$ is chosen, the process works slightly differently. Similarly,
the cooperator switches to defection if $R^{\prime}>\xi$. But for the
complementary case, instead of the defector simply switching to cooperation,
there are two possible outcomes with equal probability: either the defector
switches to cooperation, or the positive link turns into a negative link.

  * •

Resource dynamics: Along with the node and pair-connection updates, the
resource level also changes at each iteration. Using the Euler method and the
updated value of $x$, we can advance the resource dynamics as follows
$R(t+\Delta t)=R(t)+\Delta
t\Bigg{(}R(t)(1-R(t))-R(t)\Big{(}\hat{e}_{C}x(t)+\hat{e}_{D}(1-x(t))\Big{)}\Bigg{)}$,
where for the case without interaction (Fig. 2 (A)), we have $\Delta
t=\frac{1}{N}$, and for the case with interaction (Fig. 2 (B) and (C)), we
have $\Delta t=\frac{1}{L}$.

###  2.5 Macroscopic rate equations method

We can model the temporal evolution of the pair-connection densities in the
social dynamics (Fig. 2 (B) and (C)) as a multi-variable death-birth Markov
process with a master equation (Appendix. C). However, the resource dynamics
only depend on the cooperator density $x$ (Eq. 1), not on the pair-
connections. To address this, we can derive an explicit relation between the
pair-connection densities and the cooperator density for a random homogeneous
network as follows (Appendix. B):

$\displaystyle x$ $\displaystyle=$
$\displaystyle\rho_{a}+\rho_{b}+\frac{\rho_{e}+\rho_{f}}{2}$ (4)
$\displaystyle\phi$ $\displaystyle=$
$\displaystyle\rho_{c}+\rho_{d}+\frac{\rho_{e}+\rho_{f}}{2},$ (5)

where $x$ and $\phi$ are the cooperator and defector densities, respectively,
and $x+\phi=1$. This relation also applies to an Erdős-Rényi network when
$N\gg 1$ (Appendix. B). Using this approach, we can obtain the deterministic
rate equations for the temporal evolution of pair-connection densities for the
evolving dynamicl binary state interaction dynamics (Fig. 2 (C)) as follows

$\displaystyle\ \frac{d\rho_{a}}{dt}$ $\displaystyle=$
$\displaystyle-\rho_{a}+(\mu-1)\Bigg{[}-\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{a}+\frac{(1-R^{\prime})}{2\phi}\Big{(}\rho_{c}+\frac{\rho_{e}}{2}\Big{)}\rho_{f}\Bigg{]}$
(6) $\displaystyle\ \frac{d\rho_{b}}{dt}$ $\displaystyle=$
$\displaystyle(1-R^{\prime})\rho_{a}+\frac{(1-R^{\prime})}{2}\rho_{e}+(\mu-1)\Bigg{[}-\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{b}+\frac{(1-R^{\prime})}{2\phi}\Big{(}\rho_{c}+\frac{\rho_{e}}{2}\Big{)}\rho_{e}\Bigg{]}$
(7) $\displaystyle\ \frac{d\rho_{c}}{dt}$ $\displaystyle=$
$\displaystyle-\rho_{c}+(\mu-1)\Bigg{[}\frac{R^{\prime}}{2x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{f}-\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\frac{\rho_{e}}{2}\Big{)}\rho_{c}\Bigg{]}$
(8) $\displaystyle\ \frac{d\rho_{d}}{dt}$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{c}+R^{\prime}\rho_{e}+(\mu-1)\Bigg{[}\frac{R^{\prime}}{2x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{e}-\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\frac{\rho_{e}}{2}\Big{)}\rho_{d}\Bigg{]}$
(9) $\displaystyle\ \frac{d\rho_{e}}{dt}$ $\displaystyle=$
$\displaystyle-\rho_{e}+(\mu-1)\Bigg{[}\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\Big{(}\rho_{b}-\frac{\rho_{e}}{2}\Big{)}+\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\frac{\rho_{e}}{2}\Big{)}\Big{(}\rho_{d}-\frac{\rho_{e}}{2}\Big{)}\Bigg{]}$
(10) $\displaystyle\ \frac{d\rho_{f}}{dt}$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{a}+(1-R^{\prime})(\rho_{c}+\frac{\rho_{e}}{2})+(\mu-1)\Bigg{[}\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\Big{(}\rho_{a}-\frac{\rho_{f}}{2}\Big{)}+$
(12)
$\displaystyle\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\frac{\rho_{e}}{2}\Big{)}\Big{(}\rho_{c}-\frac{\rho_{f}}{2}\Big{)}\Bigg{]},$

where $\mu$ is the mean degree of the random homogeneous network, and this
coupled differential equation depends on normalized level of resource
$R^{\prime}$. An analytical relationship between $x$ and $R^{\prime}$ is
provided in Eq. D.28 when Eq. 12 reaches the stationary state. Then, an exact
analytical solution for the co-evolution of Eq.1 coupled with Eq. 12 is
obtained by computing the roots of Eq. D.35.

We can also obtain the rate equations for the social dynamics shown in Fig. 2
(B) by using the same method, as shown in Eq. C.32 (Appendix. C). When the
stationary state is reached in Eq. C.32, Eq. D.12 provides an exact analytical
relationship between $x$ and $R^{\prime}$ for an important particular case
where there is no negative interaction. Then, an exact analytical solution is
provided for the co-evolution of Eq.1 coupled with Eq. C.32 presented in Eq.
D.16.

All the above results obtained using this approach agree with those obtained
from the full integration of the differential rate equations (Eq. 12 and Eq.
C.32 coupled with Eq.1).

##  3 Results

We first analyze in detail the dynamics of each of the three model variants
(NS, QBS, EBS) on Erdős-Rényi networks. Then we briefly discuss the main
features of the dynamics on Scale-Free and Small-World networks.

For each realization of the dynamics, we generate a random complex network and
randomly assign the initial states of nodes and links. We denote by $x_{0}$
the initial density of cooperative players, by $l_{0}$ the initial density of
positive links, and by $R^{\prime}_{0}$ the initial value of normalized
resource level.

For each HES dynamics, we compare the Monte-Carlo results with the numerical
solutions of the rate equations, Eq. 12, and - as we will show - we find a
good agreement between them.

We are interested in the conditions that lead to either a depleted or a
repleted phase. A repleted phase is a phase where resources are used
sustainably (${R^{\prime}}^{\text{st}}>0$), while a depleted phase is a phase
where resources are exhausted (${R^{\prime}}^{\text{st}}=0$). This
depletion/repletion transition in resource dynamics also corresponds to a
transition from an absorbing to an active state in the temporal evolution of
pair-connection densities [25, 26]. In an absorbing state, the densities of
negative pair-connections ($\rho_{a}$,$\rho_{c}$ and $\rho_{e}$) decay
exponentially to zero, while in an active state, they reach a non-zero steady
state [25, 26]. The absorbing/active phase can be interpreted in social terms
as an emergent situation where the individuals and their connections are/are
not prone to being frozen in their current social states, respectively. In the
following subsections, we present the results of the co-evolution of Eq. 1 and
the three social dynamics shown in Fig. 2.

###  3.1 No-social interaction (NS).

We consider the HES dynamics with knowledge feedback but without any social
interaction. The temporal evolution of the cooperator density can be written
as

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle-R’x+(1-R’)(1-x).$ (13)

The normalized resource level at steady state in the co-evolution of Eq. 1 and
Eq. 13 is given by
$R^{\prime}(t\rightarrow\infty)={R^{\prime}}^{\text{st}}=\frac{1-\hat{e}_{C}}{1-2\hat{e}_{C}+\hat{e}_{D}}$,
which is always positive for any values of $\hat{e}_{C}$ and $\hat{e}_{D}$.
Finally, we note that $x^{\text{st}}=1-{R^{\prime}}^{\text{st}}$ for any value
of $\hat{e}_{C}$ and $\hat{e}_{D}$, and thus $x^{\text{st}}>0.5$. Therefore,
NS systems always lead to cooperators being the majority in a society. In
other words, for a sustainable cooperative HES to be achieved, the community
of cooperators must always form a majority, $x^{\text{st}}>0.5$. Otherwise,
the system is likely to collapse.

###  3.2 Quenched binary state interaction (QBS).

We now consider the HES dynamics with knowledge feedback and with quenched
binary state interaction, as shown in Fig. 2(B). We simulate the co-evolution
of this social interaction and the resource dynamics on an uncorrelated Erdős-
Rényi network for three different sizes $N=40,400,4000$. We plot the
normalized resource level, $R^{\prime}$, and the cooperator density $x$, at
steady state, as functions of mean degree, $\mu$, in Fig. 3(A) and (C),
respectively. We also show the corresponding curves obtained from numerically
solving the rate equations Eq. C.32. Both Monte-Carlo and the rate equations
simulations indicate that increasing $\mu$, which acts a control parameter,
leads to a transition from a depleted to a repleted state. In fact, in the
limit $\mu\rightarrow N-1$, we qualitatively recover the pure knowledge
feedback case without social structure (see Appendix. H). There is a good
agreement between the results of the rate equations and of the Monte-Carlo
simulation, especially for larger $\mu$ when the Erdős-Rényi network becomes
more homogeneous. Near the critical point, both Monte-Carlo and rate equations
results show qualitative agreement.

Figure 3: Stationary states of the normalized level of resources
${<R^{\prime}>}^{\text{st}}$ and the density of cooperators $<x>^{\text{st}}$
as functions of network mean degree $\mu$. (A) and (C) The QBS dynamics shows
a transition from a depleted state to a repleted state as $\mu$ increases. (B)
and (D) The EBS dynamics shows a similar transition. We compare two dynamics
with networks of different sizes $N=40,400,4000$, at
$(\hat{e}_{C}=0.9,\hat{e}_{D}=1.9)$. We show in all figures the results from
the Monte-Carlo simulation and numerical integration of differential rate
equations, labeled by MC and RE. The initial value of the density of
cooperators, the fraction of positive links, and the normalization level of
the resource are $x_{0}=0.25$ and $l_{0}=0.5$, and
$R^{\prime}_{0}=\frac{2}{3}$, respectively (see Appendix. I). The
${<R^{\prime}>}^{\text{st}}$ (in panel (A) and (B)) is insensitive to
different initial values, as evident in Fig. F.1.

In a finite-size system, for any value of $\mu$, the system may enter a
frustrated jammed state [28, 36, 37], where the quenched binary state
interactions prevent reaching an absorbing state. In this configuration, we
observe a coexistence of negative ($a$,$c$,$e$) and positive ($b$,$d$,$f$)
pair-connections, although the system is not dynamically active, as shown in
Fig. 4 (A)&(B). In other words, applying the update rule on a jammed state
does not change the arrangement of node states in the system, except for some
rare possible node-blinking events. In the structural balance theory, the
probability of reaching a jammed state in fully connected network decreases
rapidly as $N$ increases [28]. We observed a similar behaviour in our system
for large $\mu$. However, examining the probability of the jammed state in
sparse networks is not within the scope of this study. Finally, as evident
from Fig. 3(A) and (C), increasing network connectivity improves the
sustainability of the HES. This is true even in the emergence of jammed state
in the finite size systems.

Figure 4: (A)&(B) Jammed states of the quenched binary state dynamics (QBS)
and (C)&(D) examples of absorbing configurations of the evolving binary state
dynamics (EBS) on a finite system with extraction coefficients
($\hat{e}_{C}=0.9,\hat{e}_{D}=1.9$). The network has $N=40$ nodes and mean
degree $\mu=4$ and $\mu=8$. The initial values are $x_{0}=0.25$, $l_{0}=0.5$
(see Appendix. I), and $R^{\prime}_{0}=\frac{2}{3}$. For both case of QBS and
EBS, the fraction of defectors (blue circles) decreases as the network
connectivity increases.

###  3.3 Evolving dynamical binary state interaction (EBS)

We now study the phase structure of the co-evolution of the complex social
interaction shown in Fig. 2(C) with resource dynamics given by Eq. 1. We plot
the normalized resource level, $R^{\prime}$, and the cooperator density $x$,
at steady state, as functions of mean degree, $\mu$, for three different sizes
$N=40,400,4000$ in Fig. 3 (B) and (D), respectively. Again, in both the Monte-
Carlo simulations and rate equations of the system, for any (large enough)
sizes, we show the existence of a critical mean-degree $\mu_{c}$ below which a
depleted state occurs. In Monte Carlo simulations, for large sizes
($N=400,4000$), above $\mu_{c}$, the curves overshoot (reach a maximum value,
and for larger $\mu$ they approach a plateau), as seen in Fig. 3(B) and (D).
On the other hand, in the curves obtained from the rate equations the maximum
occurs at $\mu_{c}$ discontinuously, as shown by black-triangle markers.
Interestingly, for larger $\mu$ when the Erdős-Rényi network becomes more
homogeneous, rate equations and Monte-Carlo results are very consistent. For
small size ($N=40$), we observe a sigmoid-like behavior for $R^{\prime}$, from
a depleted phase to a fully repleted phase, as indicated by black-star marker
in Fig. 3(B).

The depletion/repletion phase transition observed in Fig. 3 (B) and (D) is
accompanied by a non-equilibrium absorbing-state transition [25, 26]. By
definition, an absorbing state is a configuration that the model can reach but
from which it cannot escape. In our case, the absorbing state is identified by
a solution when $\rho_{a}=\rho_{c}=\rho_{e}=0$. In general, in the absorbing
phase, the absorbing state is a stable solution against small perturbations,
while in the active phase, the absorbing state becomes unstable (see equation
121 of [38]). To characterize this transition, we monitored the trajectories
of $R^{\prime}$, $x$, and $\rho_{i}$ for three sets of parameters leading to
either an absorbing or an active phase. Fig. 5(A) and (D) represent the
dynamics in an absorbing phase, while (C) and (F) represent an active phase.
The middle panels represent the trajectories of the system for a set of
parameters leading to the maximum value of $R^{\prime}$ (close to the critical
point), as seen in panel (B) of Fig. 3. As shown, in an absorbing phase, the
resource level $R^{\prime}$ and the densities of negative pair-connections
($a,c,e$) decay exponentially to zero, while in an active phase they reach
nonzero steady states. However, in the upper middle panel, although the
dynamics indicate an absorbing phase, $R^{\prime}$ reaches a non-zero plateau.
We speculate that this is due to a strong fluctuating behavior close to
criticality, and related finite size effects. The top panels are obtained from
Monte-Carlo simulations, while the bottom panels are obtained from rate
equations with the same set of parameters. As shown, except for the middle
panels, the Monte-Carlo results agree well with the corresponding rate
equations results.

Figure 5: The temporal evolution of normalized resource level, $R^{\prime}$,
cooperator density, $x$, and pair-connection densities $\rho_{i}$, for
evolving dynamical binary interaction (EBS). The top panels show the results
of the Monte-Carlo simulation of a single run on an Erdős-Rényi network with
$N=4000$. The bottom panels show the results of the corresponding rate
equations. Panels (A) and (D) represent the parameters leading to an absorbing
phase, while panels (C) and (F) represent the parameters leading to an active
phase. Panels (B) and (E) represent the parameters leading to the critical
region where the value of ${R^{\prime}}^{st}$ in panel (B) of Fig. 3 reaches
the maximum value. The initial values for all dynamics are $x_{0}=0.25$,
$l_{0}=0.5$ (see Appendix. I), and $R^{\prime}_{0}=\frac{2}{3}$.

A snapshot of an absorbing configuration for the evolving binary state
dynamics (EBS) on a finite system is shown in Fig. 4 (C)&(D). As seen in this
figure, the fraction of defectors decreases with increasing network
connectivity (Fig. 3(D)), revealing the underlying mechanism behind the
depletion/repletion phase transition in finite systems. In the absorbing
phase, the average characteristic time to reach an absorbing state grows
logarithmically with system size ($<\tau>\sim\log{N}$), while in the active
phase, it grows exponentially ($<\tau>\sim\exp{N}$), as shown in Appendix. E
and Fig. E.1. Finally, as evident from panels (B) and (D) of Fig. 3,
increasing the network connectivity enhances the sustainability of the HES.

###  3.4 Robustness of Results.

We have shown that the co-evolution dynamics on Erdős-Rényi networks, with
relatively symmetric initial conditions, can lead to a transition from
depletion to repletion of CPRs, depending on the network connectivity. In
particular, we have shown that the higher mean degree of the Erdős-Rényi
system promotes CPR sustainability.

Our results are robust by varying initial conditions, the network topology,
and the social model (the similarity of the main results under both the EBS
and QBS models).

Fig. F.1 shows the co-evolution dynamics of ${<R^{\prime}>}^{st}$ and
$<x>^{st}$ as a function of mean degree for the model in Fig.2(C). The Monte-
Carlo and rate equations simulations were initialized in a wide range of
combinations of $R^{\prime}_{0}$ and $x_{0}$ and run on Erdős-Rényi networks.
Regardless of the initial conditions, ${R^{\prime}}^{st}$ evolves toward the
same stationary state dictated by the mean degree. This suggests the co-
evolution rules drive the CPR to a robust equilibrium that depends only on the
mean degree, erasing memory of initial seeds $R^{\prime}_{0}$ and $x_{0}$.
Universality across different initial configurations indicates that the
stationary state is an attractor for the resource dynamics on Erdős-Rényi
networks. Further analyses of the phase space are presented in (Appendix.
F&Appendix. G).

In Fig. F.2, we also plot ${<R^{\prime}>}^{st}$ and $<x>^{st}$ as functions of
network mean degree for four initial cooperator densities
($x_{0}=0.05,0.25,0.75,0.95$), and three network topologies, Erdős-Rényi,
Barabasi-Albert [39] and Watts-Strogatz [40]. We find that the
depletion/repletion transition is qualitatively similar for different networks
heterogeneities. However, the Barabasi-Albert network, which has higher
heterogeneity, exhibits a shift of the critical point for higher average
degree. This implies that the hubs in the network facilitate the depletion
phase. Another notable observation is that, in the depletion phase, the
resource level reaches a depleted state of ${<R^{\prime}>}^{st}=0$, regardless
of the initial condition of the dynamics (absorbing phase in Fig.F.2 in panels
(D), (E), and (F)) or the network topology (absorbing phase in Fig.F.2 panels
(A), (B), and (C)). This indicates that the depletion/repletion transition on
the complex networks are robust to initial conditions.

To vary the social model, we already study the EBS, presented in Fig. 2(C),
and its quenched version, QBS, presented in Fig. 2(B). As shown in Fig. 3 (A)
and (B), the depletion/repletion emerges in both the EBS and QBS models. This
suggests that the characteristic signature of the transition is robust and
does not rely directly on the specific dynamics of the link. However, the
types of transitions in EBS and QBS are different (see Appendix. D). As shown
in Fig. 3 (A), Monte-Carlo and rate equations simulations for the QBS model
both indicate the depletion/repletion transition is continuous (see also the
analytical arguments in Appendix. D.1), while, describing the type of
transition in EBS (Fig. 3 (B)) is a challenging task that we leave for future
work. In Appendix. D.2, we provide some analytical arguments that seem to
indicate that the rate equations, differently from the Monte-Carlo
simulations, show a discontinuous phase transition. Also in (Appendix. G) we
provide some remarks about the critical line of the depletion/repletion
transition in the phase space of ($\hat{e}_{C}$,$\hat{e}_{D}$).

Studying the QBS offers additional insights beyond EBS alone. The QBS retain
most of the key characteristics of the full socio-environmental dynamics.
However, they exhibit a more well-behaved continuous phase transition between
the active and absorbing states, rather than the unusual punctuated transition
observed in the EBS (compare panel (A) and (B) in Fig. 3). Analyzing this
continuous phase transition in the QBS may lend itself more easily to adapting
standard physics frameworks to understand phase transitions. Thus, the QBS
analysis complements the EBS picture while facilitating connection to
established theory.

Finally, it is worth mentioning that the dynamics of QBS and EBS in the
network with a high mean degree behave similarly to NS dynamics. The time
evolution of the density of cooperators $x$ in NS, QBS, and EBS can be
represented mathematically, respectively, as follows:

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle-R’x+(1-R’)(1-x),$ (14)
$\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle\mu\Big{(}-R^{\prime}(\frac{\rho_{a}+\rho_{e}}{2})+(1-R^{\prime})(\frac{\rho_{c}+\rho_{e}}{2})\Big{)},$
(15) $\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle\mu\Big{(}-R^{\prime}(\frac{\rho_{a}+\rho_{e}}{2})+(1-R^{\prime})(\frac{\rho_{c}+\frac{\rho_{e}}{2}}{2})\Big{)}.$
(16)

The derivation of Eq.15 and Eq.16 are presented in Appendix. D in Eq.D.2 and
Eq.D.18, respectively. The Eq.14 on one side and the Eqs.15-16 on the other
side have some similarity. As is evident from Eq.5, $\rho_{a}$ and $\rho_{e}$
are constituent components of $x$, while $\rho_{c}$ and $\rho_{e}$ are
constituent components of $1-x$. As fully clarified in Appendix. H, roughly
speaking, for larger $\mu$ we expect a higher activity of the system, which
results in a higher value of negative pair-connections ($\rho_{a}$,
$\rho_{c}$, and $\rho_{e}$) and, in turn, makes Eqs.15-16 more similar to
Eq.14.

###  3.5 Polarization of the Social Network Dynamics

Quantifying the degree of social polarization is a key issue [41, 42, 25, 26]
we have tackled in this work. A major challenge here is to devise a suitable
measure for polarization. We propose one that is aligned with the aim of this
study. We contend that polarization between two groups should depend not only
on their different strategies but also on their antagonism towards each other.
Polarization arises only when the inter-group connection are negative due to
the different strategies (which is not a trivial condition). For instance, if
$x=0.5$, but all the links between the groups are positive, then there is no
meaningful polarization (which means that the groups cooperate despite having
different opinions).

Hence, we define the polarization index, $r$, as the product of the
polarization due to the node-attribute (Eq. LABEL:x_pair_rel),
$2(\rho_{f}+\rho_{e})$ (which is equal to equation 3 in [41]), and the
polarization due to the link-attribute ( $\frac{\rho_{f}}{\rho_{f}+\rho_{e}}$
). Therefore, the polarization index has the simple form $r=2\rho_{f}$. This
is interesting because $\rho_{f}$ depends on initial conditions in the
depletion phase, but not in the sustainable phase.

The network polarization is shown in Fig. F.2(G), (H) and (I) for the phase
space of ( $\hat{e}_{C}=0.9,\hat{e}_{D}=1.9,\mu$ ).We observe that the
depletion phases have higher polarization than the repletion phases,
indicating that lower network connectivity leads to higher polarization. This
agrees with our intuition. However, we also find that lower extraction
coefficients ( $\hat{e}_{C},\hat{e}_{D}$ ) result in higher polarization (as
evident in Fig. G.1(C)), which is counter-intuitive (See Appendix. G).

In all model configurations (NS, QBS, and EBS), and in both the repleted and
depleted phases, the majority of the consumers deploy a cooperation strategy
(to be quantitative, $x^{\text{st}}>0.5$). In fact, the dynamics of NS always
leads to $x^{\text{st}}>0.5$. On the other hand, the diagram of
$x^{\text{st}}$ in terms of $\mu$ for QBS, and EBS are presented in Fig. 3 (C)
and (D) respectively. As evident, in the depleted phases, for $\mu>2$, the
cooperators build the majority, and in the repleted phase, both systems must
achieve a majority in excess of 90 percent cooperation (the overwhelming
majority). In particular, in the EBS case, Fig. G.1 (B)&(E) shows
$x^{\text{st}}$ as a function of ($\hat{e}_{C},\hat{e}_{D}$). It is
interesting to note that depleted phases tend to have the highest cooperator
density, $x^{\text{st}}>0.7$ (the depleted phase is identified by the black
region in panel (A)&(D)). In other words, the HES system may still not be
sustainable even if the majority of extractors are cooperators.

A likely explanation for the asymmetric behavior of $x^{st}$ is the inherent
asymmetry in the extraction terms in Eq. 1, where the cooperators’ extraction
capacity $\hat{e}_{C}$ is less than the defectors’ extraction capacity
$\hat{e}_{D}$. Specifically, the form of the extraction terms $\hat{e}_{C}x$
and $\hat{e}_{D}(1-x)$ suggests that they contribute comparably to the total
resource dynamics despite the inequality $\hat{e}_{C}<\hat{e}_{D}$. To
compensate, some underlying driver in the system skews the dynamics toward
higher steady-state cooperator densities $x^{st}>(1-x^{st})$. Through this
skew, the lower extraction capacity of cooperators is counterbalanced by their
increased representation, leading $\hat{e}_{C}x$ to be approximately equal to
$\hat{e}_{D}(1-x)$ in equilibrium. The system appears to rebalance the
asymmetry in the extraction terms by shifting the stable ratio of cooperators
to defectors.

##  4 Discussion and Conclusions

In this work, we propose a co-evolutionary model for human-environment systems
(HESs) that incorporates the effects of knowledge feedback and users’ social
interaction dynamics on the sustainability of common pool resources (CPRs).

The proposed framework allows us to capture the complex and dynamic interplay
between human behaviors and resource availability and to explore how different
scenarios affect the sustainability of CPRs. In particular, we have found that
knowledge feedback alone, when fluctuations can be neglected, allows for the
sustainability of the CPRs: the consumers self-organize to reduce the
extraction rate when the resource is low, and this is a sufficient mechanism
to prevent the tragedy of the commons. If stochastic fluctuations on the
resource dynamics are considered, then we might have CPRs collapse.

On the other hand, social interactions among players can generate rich and
diverse patterns of social organization and polarization, which can have
significant implications for the governance of CPR. In particular, both QBS
and EBS can either enhance or hinder sustainability, depending on the
connectivity and heterogeneity of the network. We have shown that there is a
critical network mean degree at which a depletion/repletion phase transition
occurs. This transition is accompanied by an absorbing/active state like of
transition in the social dynamics.

The role of hubs in facilitating defection and depletion, or the effect of
environmental feedback in promoting cooperation and repletion, are two main
results of our work.

To study the robustness and generality of our results, we have investigated
different types of game-theoretic dynamics to capture the complex and adaptive
behaviors of consumers who exploit CPRs, whose dynamics are modeled through a
logistic equation. We have also considered different network topologies and
initial conditions.

We have also measured the degree of polarization in the network, which depends
on both the node and link attributes, and we have found that lower network
connectivity and extraction coefficients lead to higher polarization. In this
regard, we have observed that even if the majority of society cooperates, the
human-environmental system can be set in depleted phases.

From a methodological point of view, we have compared the results obtained
with Monte-Carlo simulations with the analytical and numerical ones from the
rate equations. Although the two methodologies have given the same qualitative
results, the nature of the sustainable-depleted transition was different in
these two approaches. In this regard, more theoretical analyses and
investigations are needed. Other limitations remain and may define the
directions for future research. One limitation is that we have assumed that
consumers are homogeneous in their extraction coefficients and knowledge
feedback strategies, which may not be realistic in some cases. A possible
extension is to introduce heterogeneity among consumers in these aspects and
study how it affects the co-evolutionary dynamics and outcomes. Another
limitation is that we have neglected the role of governance systems or
institutions in regulating CPRs, which are known to play an important role in
enhancing or undermining sustainability. A possible extension could thus
incorporate governance systems or institutions as another sub-system in our
model, and study how they interact with consumers and resources. A third
limitation is that we have used a simple logistic equation to model the
resource dynamics, which may not capture some important features of CPRs such
as non-linearity, threshold effects, or resilience. A possible extension of
our model is thus developing more realistic models of resource dynamics and
study how they affect the co-evolutionary dynamics and outcomes.

## Acknowledgments

M.S. acknowledges the fellowship from the Department of Physics And Astronomy
‘G. Galilei’, University of Padova.

## Code availability

The GitHub link which include some part of the repository of the codes that we
used for reproducing the numerical results:
https://github.com/MeghdadSaeedian/human-environmental-system

## Appendices

##  Appendix A Interpretation of positive and negative social interactions

Here, we provide the interpretations and respective roles for each of the six
possible pair-connections $\\{a,b,c,d,e,f\\}$ in the update rule. The pair-
connections $\\{a,c,e\\}$ represent positive interactions (socially
satisfying), while $\\{b,d,f\\}$ represent negative interactions (socially
unsatisfying).

  * •

Pair-connection $a$: Two cooperators are connected with negative link, such as
intense rivalry, or distrust. We regard this configuration as a unsatisfying
pair-connection because they share a common noble goal (self-organization to
achieve sustainable HES) and a negative link would hinder this goal.
Therefore, if they act rationally, this configuration would eventually change.
Either, for the sake of a sustainable HES, the negative link turns into a
positive one ($a$ to $b$), or, one of the individuals succumbs to temptation
and switches its strategy from cooperation to defection ($a$ to $f$).
According to the update rules, the former transition is more likely to occur
in resource scarcity conditions and the latter in resource abundance.

  * •

Pair-connection $c$: Two defectors are connected with negative link. The
rationale for classifying this configuration as a unsatisfying pair-connection
is related to the concept of knowledge feedback. In resource abundance,
negative link is pointless, whereas a positive one ($c$ to $d$) can make
resource exploitation more efficient, so it is rational to interact
positively. However, in the case of resource scarcity, according to knowledge
feedback, one of the individuals may decide to change its strategy from
defection to cooperation ($c$ to $f$).

  * •

Pair-connection $e$: Two consumers with different strategies are connected
with positive link, such as collaboration or friendship. This is a conflicting
situation because a key condition for any collaboration is to have the same
strategy and approach, otherwise they cannot work together coherently. In this
case, there are three ways to achieve a positive pair-connection. According to
knowledge feedback, in resource abundance people are more likely to become
greedy ($e$ to $d$), and in resource scarcity, they tend to be more committed
to self-organization. The latter effort can be realized through the
complementary transitions ($e$ to $f$) and ($e$ to $b$) equally. In simple
terms, ($e$ to $f$) can be interpreted as a situation in which the cooperator
boycotts the defector. On the other hand, ($e$ to $b$) can be interpreted as a
situation in which the cooperator convinces the defector to become a
cooperator.

  * •

Pair-connection $b$: Two cooperators are connected with positive link. Having
a positive bond between two players with the same strategies prevents any
tendency to change the node’s or link’s state.

  * •

Pair-connection $d$: Two defectors are connected with positive link. Having a
positive bond between two players with the same strategies prevents any
tendency to change the node’s or link’s state.

  * •

Pair-connection $f$: Two consumers with different strategies are connected
with negative link, such as opposition or mutual boycott. In fact, when two
people have conflicting opinions or opposing strategies, the tendency is to
view your opposition negatively rather than forge a positive connection built
on understanding [43].

All the above-mentioned statements are summarized in the update rule shown in
Fig. 2 (C) and Appendix. B.

##  Appendix B Relationship between the density of cooperators x and the
density of pair-connections $\boldsymbol{\rho}$

Since the links among nodes in homogeneous random networks are generated
randomly with probability $p=\frac{\mu}{N}$ (see equation (3.2) [44]), we can
relate the pair-connection densities, $\rho_{i}$, to the cooperator density
$x(t)$. For example, in an Erdős-Rényi network with $N\gg 1$, the
$\rho_{a}+\rho_{b}$ is equal to the fraction of the number of pair-connections
within the cooperators over the total number of pair-connections in the
network. The former is $\frac{n(n-1)p}{2}$ and the latter is
$\frac{N(N-1)p}{2}$, where $n$ is the number of cooperators. The density of
pair-connections within the defectors, $\rho_{c}+\rho_{d}$, is obtained in a
similar way using the number of defectors $(N-n)$. The density of interface
pair-connections between the cooperators and defectors, $\rho_{e}+\rho_{f}$,
can be obtained by subtracting $\rho_{a}+\rho_{b}+\rho_{c}+\rho_{d}$ from one.
All these relations are as follows [25]:

$\displaystyle\rho_{a}+\rho_{b}$ $\displaystyle\simeq$
$\displaystyle\frac{n(n-1)}{N(N-1)}\xrightarrow{n,N\gg 1}x^{2}$
$\displaystyle\rho_{c}+\rho_{d}$ $\displaystyle\simeq$
$\displaystyle\frac{(N-n)(N-n-1)}{N(N-1)}\xrightarrow{n,N\gg 1}(1-x)^{2}$
$\displaystyle\rho_{e}+\rho_{f}$ $\displaystyle\simeq$
$\displaystyle\frac{N(N-1)-(N-n)(N-n-1)-n(n-1)}{N(N-1)}\xrightarrow{n,N\gg
1}2x(1-x)$

where $x=\frac{n}{N}$. Since the rate equations mimics the system in the
thermodynamic limit, the approximate equality sign ($\simeq$) becomes an exact
equality ($=$). Therefore, the final relation can be written as follows

$\displaystyle x$ $\displaystyle=$
$\displaystyle\rho_{a}+\rho_{b}+\frac{\rho_{e}+\rho_{f}}{2}=1-\Big{(}\rho_{c}+\rho_{d}+\frac{\rho_{e}+\rho_{f}}{2}\Big{)}.$
(B.2)

##  Appendix C Master equation and rate equations

We can use formalisms like the master equation to describe the internal
fluctuations, which are generated by the dynamics itself as a result of their
interactions and are thus intrinsic to the dynamics [45]. This description
suits the nature of the evolution of the pair-connections shown in Fig. 2.
Therefore, we attempt to write down a master equation for the time-evolution
of pair-connections. This dynamics is equivalent to a hypothetical multi-
variable birth-death Markov process. For example, in the Fig. 2 (C) for the
Monte-Carlo process, when a pair-connection $e$ is selected for updating,
technically, one pair-connection $e$ vanishes (disappears) and simultaneously
either of pair-connections $b$, $d$, or $f$ with associated probabilities
emerges. This conversion occurs through two types of updating: link-update and
node-update. Link updates only affect one pair-connection at a time (e.g. Fig.
C.1). However, in the node-updates, the status of all pair-connections
connected to the updated node are affected (e.g. Figs. C.2 and C.3). When a
node-update occurs, one pair-connection is converted due to the direct
execution of the update-rule. In addition, $\mu-1$ pair-connections connected
to the updated node are converted due to the indirect execution of the update-
rule (e.g. Figs. C.2 and C.3). For the mean-field modeling of the latter case,
at any moment in the dynamics, we assume that the pair-connections on the
network are uniformly distributed. We observed that this approximation works
well. In this appendix, we provide details on how to write down the rate
equations for the update rule shown in Fig. 2(C). A similar approach can be
used for the update rule Fig. 2(B).

The first step to write down the master equation is to define the system: the
system of interest here is the vector of the number of pair-connections on
homogeneous networks. In homogeneous random networks and large Erdős-Rényi
networks (thermodynamic limit), the total number of links is $L=N\mu/2$ ($\mu$
is the mean degree of the network) such that

$\displaystyle L$ $\displaystyle=$ $\displaystyle
L_{a}+L_{b}+L_{c}+L_{d}+L_{e}+L_{f}$ (C.1)

where $L_{i},\quad i\in\\{a,b,c,d,e,f\\}$ is the number of pair-connections.
Therefore, we can define a multi-variable system using the $L_{i}$’s as follow

$\displaystyle\mathbf{L}$ $\displaystyle=$
$\displaystyle\\{L_{a},L_{b},L_{c},L_{d},L_{e},L_{f}\\},$ (C.2)

and, a single jumps out (in) of (to) the system $\mathbf{L}$, through the
variable $L_{i}$ can be define by $\mathbf{L}-\mathbf{r}^{\pm}_{i}$, where

$\displaystyle{\mathbf{r}}_{i}^{+}$ $\displaystyle=$
$\displaystyle+\\{\delta_{a,i},\delta_{b,i},\delta_{c,i},\delta_{d,i},\delta_{e,i},\delta_{f,i}\\}$
(C.3) $\displaystyle{\mathbf{r}}_{i}^{-}$ $\displaystyle=$
$\displaystyle-\\{\delta_{a,i},\delta_{b,i},\delta_{c,i},\delta_{d,i},\delta_{e,i},\delta_{f,i}\\}\quad
i\in\\{a,b,c,d,e,f\\}$ (C.4)

$\displaystyle\delta_{j,i}=\begin{cases}1&i=j\\\
0&\text{otherwise.}\end{cases}$ (C.5)

The ${\mathbf{r}}_{i}^{+}$’s are the unit basis vectors. The $\mathbf{L}$ and
${\mathbf{r}}_{i}^{\pm}$ are extensive quantities. The corresponding intensive
quantities are defined here by $\boldsymbol{\rho}$ and
$\boldsymbol{\epsilon}^{\pm}_{i}$, where

$\displaystyle\boldsymbol{\rho}$ $\displaystyle=$
$\displaystyle\\{\rho_{a},\rho_{a},\rho_{a},\rho_{c},\rho_{d},\rho_{e},\rho_{f}\\},\quad\rho_{i}=\frac{L_{i}}{L}$
(C.6)

and

$\displaystyle\boldsymbol{\epsilon}_{i}^{+}$ $\displaystyle=$
$\displaystyle+\epsilon\\{\delta_{a,i},\delta_{b,i},\delta_{c,i},\delta_{d,i},\delta_{e,i},\delta_{f,i}\\}$
(C.7) $\displaystyle\boldsymbol{\epsilon}_{i}^{-}$ $\displaystyle=$
$\displaystyle-\epsilon\\{\delta_{a,i},\delta_{b,i},\delta_{c,i},\delta_{d,i},\delta_{e,i},\delta_{f,i}\\}\quad
i\in\\{a,b,c,d,e,f\\}$ (C.8)

where $\epsilon=\frac{1}{L}$ and note $\boldsymbol{\epsilon}_{i}^{\pm}$’s are
along the unit base vectors.

Therefore, the time evolution of the probability to find the system at state
$\mathbf{L}$ can be written as

$\displaystyle\frac{\partial P(\mathbf{L},t)}{\partial t}=$
$\displaystyle\sum_{i}\Bigg{(}\sum_{{\mathbf{r}}_{i}^{\pm}}W(\mathbf{L}-{\mathbf{r}}_{i}^{\pm},{\mathbf{r}}_{i}^{\pm},t)P(\mathbf{L}-{\mathbf{r}}_{i}^{\pm},t)$
(C.9)
$\displaystyle-\sum_{{\mathbf{r}}_{i}^{\pm}}W(\mathbf{L},{\mathbf{r}}_{i}^{\pm},t)P(\mathbf{L},t)\Bigg{)}$

where $W(\mathbf{L}-{\mathbf{r}}_{i}^{\pm},{\mathbf{r}}_{i}^{\pm},t)$ is the
transition probability per unit time to pass state
$\mathbf{L}-{\mathbf{r}}_{i}^{\pm}$ to $\mathbf{L}$ (Consequently,
$W(\mathbf{L},{\mathbf{r}}_{i}^{\pm},t)$ is to pass state $\mathbf{L}$ to
$\mathbf{L}+{\mathbf{r}}_{i}^{\pm}$). As seen, the superscript $+$ ($-$)
always indicates an increasing (a decreasing) event in density of pair
connection $i$. Therefore, the rate $W(\dots)$ holding ${\mathbf{r}}$ with
superscript + (-) implies an increasing (decreasing) rate. The transition
probability per unit time is often expressed as a scaling law of the form
below based on a dimensionless parameter $L$ [46] (see also equation II.3 of
ref. [45])

$\displaystyle W(\mathbf{L}-{\mathbf{r}}_{i}^{\pm},{\mathbf{r}}_{i}^{\pm},t)$
$\displaystyle=$ $\displaystyle
f(L)L\Bigg{[}\psi_{0}\Big{(}\boldsymbol{\rho}-\boldsymbol{\epsilon}^{\pm}_{i},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\epsilon\psi_{1}\Big{(}\boldsymbol{\rho}-\boldsymbol{\epsilon}^{\pm}_{i},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\dots\Bigg{]}$

Lets $f(L)=1$, and consider the fact that $P(\rho,t)=P(\mathbf{L},t)$ (see
equation II.4 of Ref. [45]). The master equation Eq. C.9 for the intensive
macro-variable $\boldsymbol{\rho}$ can be written as

$\displaystyle\epsilon\frac{\partial P(\boldsymbol{\rho},t)}{\partial t}$
$\displaystyle=$
$\displaystyle\sum_{i}\Bigg{(}\sum_{\boldsymbol{\epsilon}^{\pm}_{i}}\Big{[}\psi_{0}\Big{(}\boldsymbol{\rho}-\boldsymbol{\epsilon}^{\pm}_{i},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\epsilon\psi_{1}\Big{(}\boldsymbol{\rho}-\boldsymbol{\epsilon}^{\pm}_{i},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\dots\Big{]}P(\boldsymbol{\rho}-\boldsymbol{\epsilon}^{\pm}_{i},t)$
(C.11)
$\displaystyle-\sum_{\boldsymbol{\epsilon}^{\pm}_{i}}\Big{[}\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\epsilon\psi_{1}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\dots\Big{]}P(\boldsymbol{\rho},t)\Bigg{)}$

Assuming that each of the function
$\psi_{0}\Big{(}\boldsymbol{\rho}-\boldsymbol{\epsilon}^{\pm}_{i},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}$
has Taylor expansion with respect to the first argument (So far, we follow the
derivation in section III of ref. [46]). Hence we obtain the Kramers-Moyal
expansion (Equation 4.18 of [47]) for multivariate system Eq. C.11 as

$\displaystyle\frac{\partial P(\boldsymbol{\rho},t)}{\partial t}$
$\displaystyle=$
$\displaystyle\sum_{m=1}^{\infty}\frac{(-1)^{m}}{m!}\epsilon^{m-1}\sum_{{j_{1}}\dots
j_{m}}\frac{\partial}{\partial\rho_{j_{1}}\dots\partial\rho_{j_{m}}}A^{(m)}_{{j_{1}}\dots
j_{m}}(\boldsymbol{\rho},t)P(\boldsymbol{\rho},t)$ (C.12)

where $j_{\nu}\in\\{a,b,c,d,e,f\\}$ and $\nu=1\dots m$

$\displaystyle A^{(m)}_{{j_{1}}\dots j_{m}}(\boldsymbol{\rho},t)$
$\displaystyle=$
$\displaystyle\sum_{i}\Bigg{(}\sum_{\boldsymbol{\epsilon}^{\pm}_{i}}(r_{j_{1}}\dots
r_{j_{m}})\Big{[}\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\epsilon\psi_{1}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}+\dots\Big{]}\Bigg{)}$

where $r_{j_{\nu}}$ ($\nu=1\dots m$) is the $j^{th}$ components of the vector
${\mathbf{r}}_{i}^{\pm}$
($\boldsymbol{\epsilon}^{\pm}_{i}=\epsilon{\mathbf{r}}_{i}^{\pm}$). Since
${\mathbf{r}}_{i}^{+}$’s are unit basis vectors, the multiplication
$(r_{j_{1}}\dots r_{j_{m}})$ is always zero except the case that $j_{\nu}$’s
are the same. This leads to

$\displaystyle\sum_{{j_{1}}\dots
j_{m}}\frac{\partial}{\partial\rho_{j_{1}}\dots\partial\rho_{j_{m}}}\sum_{i}\sum_{\boldsymbol{\epsilon}^{\pm}_{i}}(r_{j_{1}}\dots
r_{j_{m}})\longrightarrow\sum_{j}\Big{(}\frac{\partial}{\partial\rho_{j}}\Big{)}^{m}\sum_{\boldsymbol{\epsilon}^{\pm}_{j}}({\pm
1})^{m}.$ (C.13)

Note, index $j$ appears in
$\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{\pm}_{i},t\Big{)}$
and its higher orders. Therefore, Eq. C.12 is reduced to

$\displaystyle\frac{\partial P(\boldsymbol{\rho},t)}{\partial t}$
$\displaystyle=$
$\displaystyle\sum_{m=1}^{\infty}\frac{(-1)^{m}}{m!}\epsilon^{m-1}\sum_{j}\Big{(}\frac{\partial}{\partial\rho_{j}}\Big{)}^{m}A^{(m)}_{j}(\boldsymbol{\rho},t)P(\boldsymbol{\rho},t)$
(C.14)

where

$\displaystyle A^{(m)}_{j}(\boldsymbol{\rho},t)$ $\displaystyle=$
$\displaystyle\sum_{\boldsymbol{\epsilon}^{\pm}_{j}}(\pm
1)^{m}\Big{[}\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{\pm}_{j},t\Big{)}+\epsilon\psi_{1}\Big{(}\boldsymbol{\rho},\epsilon^{\pm}_{j},t\Big{)}+\dots\Big{]}$
(C.15)

Therefore, truncating the terms for $m>2$, the Fokker-Plank equation for
multi-variable system of $\boldsymbol{\rho}$ can be obtained as follow

$\displaystyle\frac{\partial P(\boldsymbol{\rho},t)}{\partial t}$
$\displaystyle=$
$\displaystyle-\sum_{j}\frac{\partial}{\partial\rho_{j}}A^{(1)}_{j}(\boldsymbol{\rho},t)P(\boldsymbol{\rho},t)+$
(C.17)
$\displaystyle\frac{\epsilon}{2}\sum_{j}\Big{(}\frac{\partial}{\partial\rho_{j}}\Big{)}^{2}A^{(2)}_{j}(\boldsymbol{\rho},t)P(\boldsymbol{\rho},t)$

where the drift vector, $A^{(1)}_{j}(\boldsymbol{\rho},t)$, and the diagonal
component of diffusion matrix, $A^{(2)}_{j}(\boldsymbol{\rho},t)$, are as
follow

$\displaystyle A^{(1)}_{j}(\boldsymbol{\rho},t)$ $\displaystyle=$
$\displaystyle\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{+}_{j},t\Big{)}-\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{-}_{j},t\Big{)}$
(C.18) $\displaystyle A^{(2)}_{j}(\boldsymbol{\rho},t)$ $\displaystyle=$
$\displaystyle\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{+}_{j},t\Big{)}+\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{-}_{j},t\Big{)}$

Since the diffusion matrix has only diagonal elements, we can write the Ito
stochastic differential equation for Eq. C.17 using the formulas in Gardiner
[48], equations 4.3.41 and 4.3.42.

$\displaystyle\frac{d\rho_{j}(t)}{dt}=A_{j}^{(1)}(\boldsymbol{\rho},t)+\sqrt{\epsilon
A_{j}^{(2)}(\boldsymbol{\rho},t)}\zeta$ (C.19)

where $\zeta$ is a multi-variable Wiener process. As evident, for a large size
system $\epsilon\rightarrow 0$, the second term in the right side is
negligible. Therefore, we arrive at a well defined rate equations as (see also
equation 3.5.19 of Gardiner [48])

$\displaystyle\frac{d\rho_{j}(t)}{dt}=\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{+}_{j},t\Big{)}-\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{-}_{j},t\Big{)}$
(C.20)

The transition probability per unit time
$\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{+(-)}_{j},t\Big{)}$
is the total rate by which the density of variable $\rho_{j}$ is increasing
(decreasing). The total rate is calculated by adding up all the rates by which
jumping events occur per unit of time, as follow

$\displaystyle\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{+}_{j},t\Big{)}=\sum_{\nu}w_{\nu\rightarrow
j}+\sum_{\nu\omega i}w_{\nu\rightarrow\omega}^{i\rightarrow
j},\quad\nu\in\\{a,c,e\\},\quad\omega\in\\{b,d,f\\},\quad
i\in\\{a,b,c,d,e,f\\}$
$\displaystyle\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{-}_{j},t\Big{)}=\sum_{\omega}w_{j\rightarrow\omega}+\sum_{\nu\omega
i}w_{\nu\rightarrow\omega}^{j\rightarrow
i},\quad\nu\in\\{a,c,e\\},\quad\omega\in\\{b,d,f\\}.\quad
i\in\\{a,b,c,d,e,f\\}$

where the $w_{\nu\rightarrow j}$ ($w_{j\rightarrow\omega}$) is the rate at
which jumping events are made into (out of) the variable $\rho_{j}$ directly
due to execution of the update rule from pair-connection $\nu$ to $j$ ($j$ to
$\omega$). Examples of these rates can be seen in Fig. C.1. The
$w_{\nu\rightarrow\omega}^{i\rightarrow j}$
($w_{\nu\rightarrow\omega}^{j\rightarrow i}$) is the rate at which jumping
events are made into (out of) the variable $\rho_{j}$ indirectly due to the
node-update from $\nu$ to $\omega$ where the pair-connection $i$ is connected
$\nu$ and $j$ is connected $\omega$ ($j$ is connected $\nu$ and $i$ is
connected $\omega$). Example of these transition probability rates can be seen
in Figs. C.2 and C.3. Finally, by definition, the rates
$w_{\nu\rightarrow\omega}^{i\rightarrow j}$ and $w_{i\rightarrow j}$ are
obtained by multiplying the probability of the occurrence of jumping event,
$\mathpzc{p}$, and the frequency of the occurrence of that event per unit of
time, $\mathpzc{f}$. In other word $w=\mathpzc{f}.\mathpzc{p}$. Note, this
approach relies on an approximation that performs well when correlations are
not strong. Also, some of these rates are zero, for instance $w_{e\rightarrow
d}^{b\rightarrow c}=0$, or $w_{c\rightarrow b}=0$.

Figure C.1: The conversion rates of the pair-connections according to the
update rule applied in the Monte-Carlo simulation. For illustration, we show
$w_{e\rightarrow\omega}$ where $\omega\in(b,d,f)$. We explain the first
transition ($w_{e\rightarrow f}$) in detail. Given that a pair-connection $e$
with probability $\rho_{e}$ is randomly chosen for updating, it changes to
pair-connection $f$ with probability $\frac{1-R^{\prime}}{2}$. Thus, the
transition probability is $\mathpzc{p}=\frac{1-R^{\prime}}{2}\rho_{e}$.
Furthermore, since this conversion affects only one pair-connection, the event
frequency per unit time is $\mathpzc{f}=1$. Hence, $w_{e\rightarrow
f}=\frac{1-R^{\prime}}{2}\rho_{e}$. The top panel indicate a link-update,
while the middle and bottom panels capture the node-updates. Figure C.2: The
conversion rates of the pair-connections according to the update rule induced
in the Monte-Carlo simulation. For illustration, we show
$w_{\nu\rightarrow\omega}^{i\rightarrow e}$ where $i\in(a,b,c,d,e,f)$. We
explain the first transition ( $w_{a\rightarrow f}^{b\rightarrow e}$ ) in
detail. Given that a pair-connection $a$ with probability $\rho_{a}$ is
randomly chosen for updating, it changes to pair-connection $f$ with
probability $R^{\prime}$. However, the pair-connection $b$ can be linked to
the pair $a$ through the updated node. In the homogeneous network, the
probability that a pair-connection $b$ is connected to a white node
(cooperator) is $\frac{\rho_{b}}{x}$, as shown in Eq. 5. This is because, the
only pair-connections that can be linked to a white node (cooperator) are
$\\{a,b,e,f\\}$, but the contributions of these connections for $e$ and $f$
are half of $a$ and $b$. Therefore,
$\mathpzc{p}=R^{\prime}\rho_{a}\frac{\rho_{b}}{x}$. Furthermore, since each
node has a $\mu$ neighbor, node-update on a pair-connection (e.g. Fig. (C.1)
middle and bottom) will indirectly affect all other $\mu-1$ pair-connections
linked to the updated node. Thus, the frequency of occurrence of this event
per unit of time is $\mathpzc{f}=\mu-1$. Finally, $w_{a\rightarrow
f}^{b\rightarrow e}=(\mu-1)R^{\prime}\frac{\rho_{a}\rho_{b}}{x}$. Figure C.3:
The conversion rates of the pair connections according to the update rule
induced in the Monte-Carlo simulation. For illustration, we show
$w_{\nu\rightarrow\omega}^{e\rightarrow i}$ where $i\in(a,b,c,d,e,f)$. We
explain the second transition from top ( $w_{c\rightarrow f}^{e\rightarrow b}$
) in detail. Given that a pair-connection $c$ with probability $\rho_{c}$ is
randomly chosen for updating, it changes to pair-connection $f$ with
probability $1-R^{\prime}$. However, the pair-connection $e$ is linked to the
pair $c$ through a blue node. In the homogeneous network, the probability that
a pair-connection $e$ is connected to a blue node (defector) is
$\frac{\rho_{e}}{2\phi}$, as shown in Eq. 5. This is because, the only pair-
connections that can be linked to a blue node (defector) are $\\{c,d,e,f\\}$,
but the contributions of these connections of type $e$ and type $f$ are half
of the $d$ and $c$ types. Therefore,
$\mathpzc{p}=(1-R^{\prime})\rho_{c}\frac{\rho_{e}}{2\phi}$. Furthermore, since
each node has a $\mu$ neighbor, node-update on a pair-connection (e.g. Fig.
(C.1) middle and bottom) will indirectly affect all other $\mu-1$ pair-
connections linked to the updated node. Thus, the event frequency per unit
time is $\mathpzc{f}=\mu-1$. Finally, $w_{c\rightarrow f}^{e\rightarrow
b}=(\mu-1)(1-R^{\prime})\frac{\rho_{c}\rho_{e}}{2\phi}$.

For illustration, we calculate the rate equations Eq. C.20 for the pair
connection $e$ for which all the direct and indirect transition probability
rate are shown in the Fig. C.1 and Figs. C.2 and C.3, respectively.

$\displaystyle\frac{d\rho_{e}}{dt}$ $\displaystyle=$
$\displaystyle\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{+}_{e},t\Big{)}-\psi_{0}\Big{(}\boldsymbol{\rho},\boldsymbol{\epsilon}^{-}_{e},t\Big{)}$
$\displaystyle\frac{d\rho_{e}}{dt}$ $\displaystyle=$
$\displaystyle\Bigg{[}\Big{(}w^{b\rightarrow e}_{a\rightarrow
f}+w^{d\rightarrow e}_{c\rightarrow f}+w^{d\rightarrow e}_{e\rightarrow
b}+w^{b\rightarrow e}_{e\rightarrow d}\Big{)}\Bigg{]}-$
$\displaystyle\Bigg{[}\Big{(}w_{e\rightarrow f}+w_{e\rightarrow
b}+w_{e\rightarrow d}\Big{)}+\Big{(}w^{e\rightarrow d}_{a\rightarrow
f}+w^{e\rightarrow b}_{c\rightarrow f}+w^{e\rightarrow b}_{e\rightarrow
b}+w^{e\rightarrow d}_{e\rightarrow d}\Big{)}\Bigg{]}$

Finally, following the same approach for other pair-connections we arrive at

$\displaystyle\ \frac{d\rho_{a}}{dt}$ $\displaystyle=$
$\displaystyle-\rho_{a}+(\mu-1)\Bigg{[}R^{\prime}\Big{(}-\frac{\rho_{a}^{2}}{x}-\frac{\rho_{e}\rho_{a}}{x}\Big{)}+(1-R^{\prime})\Big{(}+\frac{\rho_{c}\rho_{f}}{2\phi}+\frac{1}{2}\frac{\rho_{e}\rho_{f}}{2\phi}\Big{)}\Bigg{]}$
(C.21) $\displaystyle\ \frac{d\rho_{b}}{dt}$ $\displaystyle=$
$\displaystyle(1-R^{\prime})\rho_{a}+\frac{(1-R^{\prime})}{2}\rho_{e}+(\mu-1)\Bigg{[}R^{\prime}\Big{(}-\frac{\rho_{a}\rho_{b}}{x}-\frac{\rho_{e}\rho_{b}}{x}\Big{)}+(1-R^{\prime})\Big{(}\frac{\rho_{c}\rho_{e}}{2\phi}+\frac{1}{2}\frac{\rho_{e}^{2}}{2\phi}\Big{)}\Bigg{]}$
(C.22) $\displaystyle\ \frac{d\rho_{c}}{dt}$ $\displaystyle=$
$\displaystyle-\rho_{c}+(\mu-1)\Bigg{[}R^{\prime}\Big{(}\frac{\rho_{a}\rho_{f}}{2x}+\frac{\rho_{e}\rho_{f}}{2x}\Big{)}+(1-R^{\prime})\Big{(}-\frac{\rho_{c}^{2}}{\phi}-\frac{1}{2}\frac{\rho_{e}\rho_{c}}{\phi}\Big{)}\Bigg{]}$
(C.23) $\displaystyle\ \frac{d\rho_{d}}{dt}$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{c}+R^{\prime}\rho_{e}+(\mu-1)\Bigg{[}R^{\prime}\Big{(}\frac{\rho_{a}\rho_{e}}{2x}+\frac{\rho_{e}^{2}}{2x}\Big{)}+(1-R^{\prime})\Big{(}-\frac{\rho_{c}\rho_{d}}{\phi}-\frac{1}{2}\frac{\rho_{e}\rho_{d}}{\phi}\Big{)}\Bigg{]}$
(C.24) $\displaystyle\ \frac{d\rho_{e}}{dt}$ $\displaystyle=$
$\displaystyle-\rho_{e}+(\mu-1)\Bigg{[}R^{\prime}\Big{(}-\frac{\rho_{a}\rho_{e}}{2x}+\frac{\rho_{a}\rho_{b}}{x}-\frac{\rho_{e}^{2}}{2x}+\frac{\rho_{e}\rho_{b}}{x}\Big{)}$
(C.25)
$\displaystyle+(1-R^{\prime})\Big{(}\frac{\rho_{c}\rho_{d}}{\phi}-\frac{\rho_{c}\rho_{e}}{2\phi}+\frac{1}{2}\frac{\rho_{e}\rho_{d}}{\phi}-\frac{1}{2}\frac{\rho_{e}^{2}}{2\phi}\Big{)}\Bigg{]}$
$\displaystyle\ \frac{d\rho_{f}}{dt}$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{a}+(1-R^{\prime})(\rho_{c}+\frac{\rho_{e}}{2})+(\mu-1)\Bigg{[}R^{\prime}\Big{(}-\frac{\rho_{a}\rho_{f}}{2x}+\frac{\rho_{a}^{2}}{x}-\frac{\rho_{e}\rho_{f}}{2x}+\frac{\rho_{e}\rho_{a}}{x}\Big{)}$
(C.26)
$\displaystyle+(1-R^{\prime})\Big{(}\frac{\rho_{c}^{2}}{\phi}-\frac{\rho_{c}\rho_{f}}{2\phi}+\frac{1}{2}\frac{\rho_{e}\rho_{c}}{\phi}-\frac{1}{2}\frac{\rho_{e}\rho_{f}}{2\phi}\Big{)}\Bigg{]}$

For the quenched binary state dynamics (Fig. 2(B)), the rate equations can be
derived as we did for Eq. C.26. Therefore, the rate equations for the quenched
dynamics can be written as

$\displaystyle\ \frac{d\rho_{a}}{dt}$ $\displaystyle=$
$\displaystyle-R^{\prime}\rho_{a}+(\mu-1)\Bigg{[}-\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{a}+\frac{(1-R^{\prime})}{2\phi}\Big{(}\rho_{c}+\rho_{e}\Big{)}\rho_{f}\Bigg{]}$
(C.27) $\displaystyle\ \frac{d\rho_{b}}{dt}$ $\displaystyle=$
$\displaystyle(1-R^{\prime})\rho_{e}+(\mu-1)\Bigg{[}-\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{b}+\frac{(1-R^{\prime})}{2\phi}\Big{(}\rho_{c}+\rho_{e}\Big{)}\rho_{e}\Bigg{]}$
(C.28) $\displaystyle\ \frac{d\rho_{c}}{dt}$ $\displaystyle=$
$\displaystyle-(1-R^{\prime})\rho_{c}+(\mu-1)\Bigg{[}\frac{R^{\prime}}{2x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{f}-\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\rho_{e}\Big{)}\rho_{c}\Bigg{]}$
(C.29) $\displaystyle\ \frac{d\rho_{d}}{dt}$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{e}+(\mu-1)\Bigg{[}\frac{R^{\prime}}{2x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\rho_{e}-\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\rho_{e}\Big{)}\rho_{d}\Bigg{]}$
(C.30) $\displaystyle\ \frac{d\rho_{e}}{dt}$ $\displaystyle=$
$\displaystyle-\rho_{e}+(\mu-1)\Bigg{[}\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\Big{(}\rho_{b}-\frac{\rho_{e}}{2}\Big{)}+\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\rho_{e}\Big{)}\Big{(}\rho_{d}-\frac{\rho_{e}}{2}\Big{)}\Bigg{]}$
(C.31) $\displaystyle\ \frac{d\rho_{f}}{dt}$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{a}+(1-R^{\prime})\rho_{c}+$ (C.32)
$\displaystyle(\mu-1)\Bigg{[}\frac{R^{\prime}}{x}\Big{(}\rho_{a}+\rho_{e}\Big{)}\Big{(}\rho_{a}-\frac{\rho_{f}}{2}\Big{)}+\frac{(1-R^{\prime})}{\phi}\Big{(}\rho_{c}+\rho_{e}\Big{)}\Big{(}\rho_{c}-\frac{\rho_{f}}{2}\Big{)}\Bigg{]}$

##  Appendix D Analytical Solution for the Rate Equations

Analytical solutions are presented here for the coevolution of resource
dynamics (Eq. 1) and rate equations of QBS (Eq. C.32) and EBS (Eq. 12). These
solutions enable us to analytically study the type of transitions that occur.
For QBS, we found an exact analytical solution for an important particular
case of Eq. C.32 where all the quenched interactions are negative. This
solution was then coupled in an exact way with the dynamics given by Eq. 1.
For EBS we found an exact analytical solution for Eq. 12. This exact solution
was then coupled with the dynamics provided by Eq. 1 and solved numerically.

In the case of QBS, the analytical results shown in Fig. D.1(A) describe a
particular system that retains the key property of a phase transition that is
present in the general case. Thus, this particular case can be considered a
close approximation of the general QBS dynamics, and can be studied in place
of the more complex general system. The similarity is evident through
comparison of Fig. D.1(A) and Fig. 3(A), which both exhibit the continuous
phase transition characteristic of the general QBS case. In case EBS, the
analytical results given by Fig. D.1(B) are quantitatively in excellent
agreement with the numerical integration of the differential equations
presented by black triangle marker in Fig. 3(B).

Figure D.1: Comparison of analytical solution and numerical integration of
differential equations for normalized resource levels as a function of mean
degree. Panel (A) shows results for QBS (Eq. D.16) and panel (B) shows results
for EBS (Eq. D.35). In both panels, there is excellent agreement between the
analytical solution and numerical integration. The dynamics describing the
system in panel (A) represent a particular case of general QBS dynamics which
retains the key feature of a continuous phase transition, as shown in Fig.
3(A). Thus panel (A) serves as an important approximation of the general case.
In Panel (B), the plateau amplitude derived analytically demonstrates
quantitative agreement with the results from numerical integration, as
evidenced by the black triangle marker in Fig. 3(B). Finally, the discrepancy
between the critical points obtained from the analytical and numerical
integration methods is an area requiring further investigation. One potential
explanation is that the non-zero analytical solutions (lines) in the regions
that numerical integration (circles) yields zero ($R^{\prime st}=0$) is
unstable.

###  D.1 A solution for co-evolution of recourse dynamics (Eq. 1) and rate
equation of QBS (Eq. C.32)

We begin by expressing the derivative of the cooperator density, $\dot{x}$,
with respect to the density of negative pair-connections. Differentiating both
sides of Eq. 5 allows us to obtain the derivative of the cooperator density
as:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$
$\displaystyle\frac{d\rho_{a}}{dt}+\frac{d\rho_{b}}{dt}+\frac{1}{2}\Big{(}\frac{d\rho_{e}}{dt}+\frac{d\rho_{f}}{dt}\Big{)}$
(D.1) $\displaystyle\ \dot{x}$ $\displaystyle=$
$\displaystyle\mu\Big{(}-\frac{R^{\prime}}{2}(\rho_{a}+\rho_{e})+\frac{1-R^{\prime}}{2}(\rho_{c}+\rho_{e})\Big{)}.$
(D.2)

Throughout this section, all variables are at their stationary states,
however, for simplicity, we omit the superscript of stationary. Moreover, the
case of $R=\rho_{a}=\rho_{c}=\rho_{e}=0$ is a trivial solution. In the
following we focus on the solutions simultaneously positive
$\rho_{a}+\rho_{c}+\rho_{e}>0$ and $R>0$ (replete and active phase).

At the stationary state, Eq. D.2 gives a useful relation between the negative
pair-connections as

$\displaystyle\rho_{c}+\rho_{e}=\frac{R^{\prime}}{1-R^{\prime}}\Big{(}\rho_{a}+\rho_{e}\Big{)}$
(D.3)

In the second step, using Eq. D.3, we try to find independent solution for
$\rho_{i}$’s and $x$ in terms of $\mu$ and $R$. The system of equations
involving in this step can be sorted as

$\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-R^{\prime}\rho_{a}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}-\frac{\rho_{a}}{x}+\frac{\rho_{f}}{2\phi}\Bigg{]}$
(D.4) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle(1-R^{\prime})\rho_{e}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}-\frac{\rho_{b}}{x}+\frac{\rho_{e}}{2\phi}\Bigg{]}$
(D.5) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-(1-R^{\prime})\rho_{c}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{\rho_{f}}{2x}-\frac{\rho_{c}}{\phi}\Bigg{]}$
(D.6) $\displaystyle\ 0$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{e}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{\rho_{e}}{2x}-\frac{\rho_{d}}{\phi}\Bigg{]}$
(D.7) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-\rho_{e}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{1}{x}(\rho_{b}-\frac{\rho_{e}}{2})+\frac{1}{\phi}(\rho_{d}-\frac{\rho_{e}}{2})\Bigg{]}$
(D.8) $\displaystyle\ 0$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{a}+(1-R^{\prime})\rho_{c}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{1}{x}(\rho_{a}-\frac{\rho_{f}}{2})+\frac{1}{\phi}(\rho_{c}-\frac{\rho_{e}}{2})\Bigg{]}$
(D.9) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle\mu\Bigg{[}-\frac{R^{\prime}}{2}(\rho_{a}+\rho_{e})+\frac{1-R^{\prime}}{2}(\rho_{c}+\rho_{e})\Bigg{]}$
(D.10) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-1+\rho_{a}+\rho_{b}+\rho_{c}+\rho_{d}+\rho_{e}+\rho_{f}$ (D.11)

While the analytical solution to Eq.D.11 is cumbersome, we can study a
simplified scenario that still captures key system behaviors. Particularly, if
we assume no positive connections exist ($\rho_{b}=\rho_{e}=\rho_{d}=0$), the
system remains active while preserving symmetries in interaction types between
the within-cooperators and within-defectors interactions. The analytical
solution of this important particular case yields results very similar to the
general case of QBS demonstrated in Fig.3(A). Focusing on this particular case
allows analytical tractability while eliciting the core system dynamics.

Substituting $\rho_{e}=0$ in Eq.D.11 and employing a standard software of
mathematical analysis (Wolfram Mathematica, see section Code Availability), an
independent analytical solution for $x$ is obtained as

$\displaystyle
x(R^{\prime},\mu)=\frac{g_{0}(R^{\prime},\mu)+\sqrt{g_{1}(R^{\prime},\mu)}}{g_{2}(R^{\prime},\mu)}$
(D.12)

where

$\displaystyle\ g_{0}(R^{\prime},\mu)$ $\displaystyle=$
$\displaystyle-1-2q+2{R^{\prime}}+2q{R^{\prime}}$ (D.13) $\displaystyle\
g_{1}(R^{\prime},\mu)$ $\displaystyle=$ $\displaystyle
1-4R+4q^{2}{R^{\prime}}+4{R^{\prime}}^{2}-4q^{2}{R^{\prime}}^{2}$ (D.14)
$\displaystyle\ g_{2}(R^{\prime},\mu)$ $\displaystyle=$ $\displaystyle
2(-1-q+2{R^{\prime}}+2q{R^{\prime}})$ (D.15)

and q=$\mu-1$. The Eq. D.12 is a analytical solution for the rate equations
Eq. C.32. In order to find a solution for co-evolution of QBS (Eq. C.32) and
resource dynamics (Eq. 1) we can substitute Eq. D.12 in the resource dynamics
at the stationary state when $R>0$ as follow

$\displaystyle
0=1-(1-\hat{e}_{C})R^{\prime}-x(R^{\prime},\mu)(\hat{e}_{C}-\hat{e}_{D})-\hat{e}_{D}.$
(D.16)

Solving Eq.D.16 analytically derives an independent relationship between
$R^{\prime}$ and the control parameters ($\mu$, $\hat{e}_{C}$, $\hat{e}_{D}$).
Since this analytical solution is lengthy, it is presented in the GitHub
repository Mathematica file detailing the full solution to Eq.D.16. This
separate file allows the interested reader to examine the full analytical
details, while avoiding cluttering the current discussion. The key insight is
that an independent functional connection between $R^{\prime}$ and ($\mu$,
$\hat{e}_{C}$, $\hat{e}_{D}$) emerges from analytically solving the dynamics
of the co-evolving quenched system.

The solution to Eq.D.16 is plotted as a function of $\mu$ in Fig.D.1(A). There
is evident qualitative agreement between the analytical solution of this
particular case and the general continuous phase transition seen in Fig. 3(A).

###  D.2 A solution for co-evolution resource dynamics (Eq. 1) and rate
equation of EBS (Eq. 12)

We first derive an exact analytical solution to Eq.12. This establishes the
baseline social evolutionary dynamics in isolation. We then numerically solve
for the co-evolution of Eq.12 coupled to the resource dynamics captured by Eq.
1.

We start by expressing the derivative of the cooperator density in terms of
the negative pair connections as

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$
$\displaystyle\frac{d\rho_{a}}{dt}+\frac{d\rho_{b}}{dt}+\frac{1}{2}\Big{(}\frac{d\rho_{e}}{dt}+\frac{d\rho_{f}}{dt}\Big{)}$
(D.17) $\displaystyle\ \dot{x}$ $\displaystyle=$
$\displaystyle\mu\Big{(}-\frac{R^{\prime}}{2}(\rho_{a}+\rho_{e})+\frac{1-R^{\prime}}{2}(\rho_{c}+\frac{\rho_{e}}{2})\Big{)}.$
(D.18)

Next, we rewrite Eq. D.18 at the stationary state and obtain a useful relation
between the negative pair-connections as

$\displaystyle\rho_{c}+\frac{\rho_{e}}{2}=\frac{R^{\prime}}{1-R^{\prime}}\Big{(}\rho_{a}+\rho_{e}\Big{)}$
(D.19)

In the second step, using Eq. D.19, we try to find independent solution for
$\rho_{i}$’s and $x$ in terms of $\mu$ and $R$. The system of equations
required to find the independent solutions are as follow

$\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-\rho_{a}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}-\frac{\rho_{a}}{x}+\frac{\rho_{f}}{2\phi}\Bigg{]}$
(D.20) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle(1-R^{\prime})\rho_{a}+\frac{(1-R^{\prime})}{2}\rho_{e}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}-\frac{\rho_{b}}{x}+\frac{\rho_{e}}{2\phi}\Bigg{]}$
(D.21) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-\rho_{c}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{\rho_{f}}{2x}-\frac{\rho_{c}}{\phi}\Bigg{]}$
(D.22) $\displaystyle\ 0$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{c}+R^{\prime}\rho_{e}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{\rho_{e}}{2x}-\frac{\rho_{d}}{\phi}\Bigg{]}$
(D.23) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-\rho_{e}+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{1}{x}(\rho_{b}-\frac{\rho_{e}}{2})+\frac{1}{\phi}(\rho_{d}-\frac{\rho_{e}}{2})\Bigg{]}$
(D.24) $\displaystyle\ 0$ $\displaystyle=$ $\displaystyle
R^{\prime}\rho_{a}+(1-R^{\prime})(\rho_{c}+\frac{\rho_{e}}{2})+(\mu-1)R^{\prime}(\rho_{a}+\rho_{e})\Bigg{[}\frac{1}{x}(\rho_{a}-\frac{\rho_{f}}{2})+\frac{1}{\phi}(\rho_{c}-\frac{\rho_{f}}{2})\Bigg{]}$
(D.25) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle\mu\Bigg{[}-\frac{R^{\prime}}{2}(\rho_{a}+\rho_{e})+\frac{1-R^{\prime}}{2}(\rho_{c}+\frac{\rho_{e}}{2})\Bigg{]}$
(D.26) $\displaystyle\ 0$ $\displaystyle=$
$\displaystyle-1+\rho_{a}+\rho_{b}+\rho_{c}+\rho_{d}+\rho_{e}+\rho_{f}$ (D.27)

Employing a standard software of mathematical analysis (Wolfram Mathematica,
see section Code Availability), an independent analytical solution for $x$ is
obtained as

$\displaystyle
x(R^{\prime},\mu)=\frac{g_{0}(R^{\prime},\mu)-\sqrt{g_{1}(R^{\prime},\mu)}}{g_{2}(R^{\prime},\mu)}$
(D.28)

where

$\displaystyle\ g_{0}(R^{\prime},\mu)$ $\displaystyle=$
$\displaystyle{R^{\prime}}+3qR^{\prime}-7R^{2}-8q{R^{\prime}}^{2}+13{R^{\prime}}^{3}+q{R^{\prime}}^{3}-4{R^{\prime}}^{4}+8q{R^{\prime}}^{4}-4qR^{5}$
(D.29) $\displaystyle\ g_{1}(R^{\prime},\mu)$ $\displaystyle=$ $\displaystyle
2-22{R^{\prime}}+97{R^{\prime}}^{2}-9q^{2}{R^{\prime}}^{2}-214{R^{\prime}}^{3}+84q^{2}{R^{\prime}}^{3}+235{R^{\prime}}^{4}-280q^{2}{R^{\prime}}^{4}$
(D.30)
$\displaystyle-98{R^{\prime}}^{5}+398q^{2}{R^{\prime}}^{5}-17{R^{\prime}}^{6}-131q^{2}{R^{\prime}}^{6}+22{R^{\prime}}^{7}-298q^{2}{R^{\prime}}^{7}-4{R^{\prime}}^{8}$
(D.31)
$\displaystyle+388q^{2}{R^{\prime}}^{8}-184q^{2}{R^{\prime}}^{9}+32q^{2}{R^{\prime}}^{10}$
(D.32) $\displaystyle\ g_{2}(R^{\prime},\mu)$ $\displaystyle=$
$\displaystyle-2+14{R^{\prime}}+6q{R^{\prime}}-36{R^{\prime}}^{2}-28q{R^{\prime}}^{2}+38{R^{\prime}}^{3}+40q{R^{\prime}}^{3}-10{R^{\prime}}^{4}$
(D.33) $\displaystyle-22q{R^{\prime}}^{4}+4q{R^{\prime}}^{5}$ (D.34)

and q=$\mu-1$. The Eq. D.28 is an exact analytical solution for the rate
equations Eq. 12. In order to find a solution for co-evolution of EBS (Eq. 12)
and resource dynamics (Eq. 1) we can substitute Eq. D.28 in the resource
dynamics at the stationary state when $R>0$ as follow

$\displaystyle
0=1-(1-\hat{e}_{C})R^{\prime}-x(R^{\prime},\mu)(\hat{e}_{C}-\hat{e}_{D})-\hat{e}_{D}.$
(D.35)

Solving Eq.D.35 derives an independent relationship between $R^{\prime}$ and
the control parameters ($\mu$, $\hat{e}_{C}$, $\hat{e}_{D}$). The details to
solving the Eq.D.35 presented in GitHub (see Code Availability). Note that by
truncating on the lower orders in Eq. D.34, an approximate analytical solution
can be obtained, but here we tried to keep all terms.

The solution to Eq. D.35, shown in Fig. D.1(B), displays discontinuity above
the critical point. The solution also approaches a nonzero plateau whose
amplitude quantitatively matches numerical integration (also see the black
triangle marker in Fig. 3(B)). However, further investigation is needed to
clarify the discrepancy in critical points between the analytical and
numerical integration approaches.

##  Appendix E Fluctuation-driven absorbing states in finite size systems

The existence of absorbing states implies that a finite system will eventually
fluctuate into absorbing state. However, the key is determining how the
characteristic absorption time, $\tau$, depends on system size $N$. Fig. E.1
show the average of the characteristic time as a function of system size. The
mean time $\mu$, plotted in panel A, representing relatively smaller mean
degree ($<\mu>=4$), show a logarithmic increase with system size
$<\tau>\sim\text{log}(\beta N)$, while in panels B and C ($<\tau>=12,16$),
grows exponentially with $N$ as $<\tau>\sim e^{\alpha N}$. For sufficiently
large N, an exponential form reflects the rarity of absorption events. This
exponential divergence signifies that decay to the absorbing state requires
enormously long times for large systems. From a statistical mechanics
perspective, the active phase persists indefinitely in the thermodynamic limit
as $N\rightarrow\infty$. Fluctuation-driven absorption vanishes, and the
active state is a true macroscopic phase, stabilized exponentially with size.

Figure E.1: Average characteristic time to reach an absorbing state as a
function of system size for two sets of parameters leading to either an
absorbing phase in panel (A) or an active phase in panels (B) and (C). As
shown, the curves in the absorbing (active) phase are approximately
logarithmic (exponential). All the computed values are shown with their error
bars.

##  Appendix F Effect of initial condition and network topology in the
depletion/repletion transition

The dynamic behavior of complex systems can depend sensitively on both the
initial conditions and the structure of the underlying interaction network.
Here we provide the phase transition diagram for the co-evolution EBS dynamics
given in Fig. 2 (C) for a wide range of initial condition in Fig. F.1 and
three underlying network topology in Fig. F.2.

Figure F.1: Effects of initial conditions on ${<R^{\prime}>}^{\text{st}}$,
$<x>^{\text{st}}$, and $<\rho_{f}>^{\text{st}}$. Using the  EBS dynamics model
outlined in Fig. 2(C), we simulated across a broad spectrum of initial
conditions on Erdős-Rényi networks. Specifically, we systematically varied the
initial settings of $x_{0}\in\\{0.1,0.5,0.9\\}$ and
$R^{\prime}_{0}\in\\{0.1,0.5,0.9\\}$, covering the nine combinations of these
parameterized starting states. Based on these input configurations, we
executed Monte-Carlo simulations alongside numerical integration of coupled
differential rate equations governing the dynamics. The consistent legend is
only presented in panel (B). In panel (A), we indicate the depletion/repletion
transition in ${<R^{\prime}>}^{\text{st}}$. As is evident, in Monte Carlo
simulation, all curves collapse in the same steady state, except those
starting from $R^{\prime}_{0}=0.1$, where the trajectories fluctuate. Even
these exceptional cases exhibit attraction to the common steady state
statistics when examined under an expanded lens, by aggregating many
stochastic realizations and simulating larger system sizes. In panel (B) and
(C), we plot the ${<x>}^{\text{st}}$ and ${<\rho_{f}>}^{\text{st}}$,
respectively. As evident, in the absorbing phase (here, $\mu\lessapprox 5$),
the density of cooperators in (B) and the level of polarization in (C) depend
on the initial condition $R^{\prime}_{0}$ (for open access code see Code
Availability section).

Figure F.2: Effects of network topology on ${<R^{\prime}>}^{\text{st}}$ and
$<x>^{\text{st}}$. (A) and (D) Results for the Erdős-Rényi network. (B) and
(E) Results for the Barabasi-Albert network. (C) and (F) Results for the
Watts-Strogatz network. We show how these quantities vary with network mean
degree $\mu$ for four initial cooperator densities (
$x_{0}=0.05,0.25,0.75,0.95$ ) and three network topologies: Erdős-Rényi,
Barabasi-Albert and Watts-Strogatz (with rewiring probability 0.05). The
initial values are $l_{0}=0.5$ and $R^{\prime}_{0}=\frac{2}{3}$ for all cases.
The network topology does not qualitatively affect the deplete/replete
transition, but there is a slight correlation between the critical values and
the network heterogeneity. The Barabasi-Albert network, which is the most
heterogeneous, has the largest critical value, while the Watts-Strogatz
network, which is the least heterogeneous, has the smallest critical value. In
the absorbing phase, ${<R^{\prime}>}^{\text{st}}$ and $<x>^{\text{st}}$ depend
on the initial conditions, but in the active phase they do not. The Monte-
Carlo and rate equations results agree well. (G), (H) and (I)
${<\rho_{f}>}^{\text{st}}$ as a function of $\mu$, which measures the level of
polarization, $r=2\rho_{f}$. The ${<\rho_{f}>}^{\text{st}}$ depends on the
initial conditions in the depleted phase, but not in the repleted phase. The
depleted phase also shows a higher degree of polarization than the repleted
phase.

##  Appendix G Depletion/Repletion transition in ($\hat{e}_{C},\hat{e}_{D}$)
phase space

The critical line of depletion/repletion transition is extended in phase space
of ($\hat{e}_{C},\hat{e}_{D}$) as presented in Fig. G.1(A)&(D). The essence of
this transition seems to be discontinuous due to the abrupt change in the
normalized level of resource at the stationary state, $<R^{\prime}>^{st}$. In
order to make an accurate claim on the type of transition, extensive finite
size analysis is required. This is beyond the scope of this manuscript, so we
leave it here for future research.

In Fig. G.1(B)&(E), the average of density of cooperators is presented in
color-code in ($\hat{e}_{C},\hat{e}_{D}$) phase space. As evident, even in the
depleted phase, the density of cooperators almost always builds the majority
of the population ($x^{st}\gtrsim 0.5$). In other words, as also shown in the
main text, even if the majority of the population cooporates, the HER system
is not guaranteed to be sustainable.

The case shown in Fig. G.1(B) demonstrates a counterintuitive result - the
density of cooperators in the depleted phase (characterized by high
$\hat{e}_{C}$ and $\hat{e}_{D}$) is actually higher than the density of
cooperators in the repleted phase. This seems paradoxical at first glance.
However, if we plot the steady-state density of cooperators ($x^{st}$) as a
function of the network mean degree ($\mu$) for any given combination of
$\hat{e}_{C}$ and $\hat{e}_{D}$ values, we see that $x^{st}$ is uniformly
lower in the absorbing phase than in the active phase. This trend is
illustrated in Fig. G.2(B). This illustrates that the density of cooperators
decreases when transitioning from the active state to the absorbing state.
However, the cooperator density remains higher than that of defectors in both
phases. Consider the average values, for instance.

In Fig. G.1(C)&(F), we present the $<\rho_{f}>^{st}$ as index of polarization
(see section 3.5). As seen, the lower extraction coefficients lead to higher
polarization.

Figure G.1: Depletion/Repletion transition in ($\hat{e}_{C},\hat{e}_{D}$)
phase space. In this figure, we plot the color-coded order parameters ( $R$,
$x$, and $f$ ) as functions of the control parameters (
$\hat{e}_{C}$,$\hat{e}_{D}$, and $\mu=4$ ) in a raster plot. The simulations
run for $N=4000$. The top (bottom) panels show the results from Monte Carlo
(rate equations) simulations. (A) and (D) The deplete/replete transition in
the plane of extraction coefficients ( $\hat{e}_{C}$ and $\hat{e}_{D}$ ). The
deplete (replete) phase is identified where $<R^{\prime}>^{st}=0$
($<R^{\prime}>^{st}>0$). In other words, a distinctive border between black
and cooper regions marks the critical line of the transition. The panel (A)
suggests a discontinuous phase transition due to the abrupt change in color
codes at the transition boundary. Further simulations (not presented here)
reveal the critical line shifts to larger values of $\hat{e}_{C}$ and
$\hat{e}_{D}$ as $\mu$ increases. (B) and (E) It is almost always the case
that the cooperators outnumber the population by more than half. This is scary
because sometimes even above eighty percent of cooperators does not guarantee
the sustainability HES. (C) and (F) Lower extraction coefficients lead to
higher polarization. Figure G.2: (A) shows the distribution of the
characteristic absorption time $\tau$ for network mean degrees of $\mu=4$ and
$\mu=12$, representing the absorbing and active phases respectively. (B)
displays the distribution of the density of cooperators $x$ for the same two
network conditions leading to the absorbing ($\mu=4$) and active ($\mu=12$)
phases. Comparing the two panels illustrates how the absorptive dynamics that
emerge at lower mean degrees impact both the time to reach a frozen state as
well as the eventual composition of cooperators within the frozen states.
Finally note, the distribution of the characteristic absorption time $\tau$
differs markedly between the active phase ($\mu=12$) and the absorbing phase.
In the active phase, $\tau$ has an extended, multi-scaling distribution,
unlike the more limited distribution seen in the absorbing phase. A log-log
plot is therefore used in panel (A) to depict the $\tau$ distributions for
both the active and absorbing phases on the same axes. This allows
visualization of the full range of timescales exhibited in the active phase,
which spans orders of magnitude longer than the absorption times typically
observed in the absorbing phase.

##  Appendix H Individual vs. social dynamics

When the social network is densely connected, the social dynamics of the QBS &
EBS models depicted in Fig. 2(B)&(C) exhibit behavior similar to that of the
individual dynamics NS shown in Fig. 2(A).

By definition the NS is a dynamic solely on the nodes, neglecting the presence
of interactions among users. This mechanism holds when a consumer makes a
decision without the influence of other consumers. On the other hand, the case
of QBS & EBS is defined solely for the links that are conditioned by the
strategy of consumers on their sides. These mechanisms are true when we
attempt to model the dependence of a consumer’s decision on the status of
other consumers (link and node).

The key concept of social interactions in QBS & EBS is the reluctance of
individual to update their strategy stemming from the positive bound with the
like-minded friends (pair-connection $b$ and $d$) or negative connections with
opponent (pair-connection $f$). Therefore, when there is a high density of
$b$, $d$, and $f$ pair-connections within the network, the level of activity
of both nodes and links tends to be low (and vice-versa). The update rules
presented in Fig. 2(B)&(C) attempt to increase the density of positive pair-
connections ($b$,$d$,$f$) by converting the negative pair-connections
($a$,$c$,$e$) to the positive ones ($b$,$d$,$f$). However, the node-updates,
as evidently shown in Fig. C.2 and Fig. C.3, change the status of all pair-
connections that are connected to the updated-nodes, slowing down the increase
of positive interactions ($b$,$d$,$f$). This effect is even amplified when
connectivity in the network increases. Therefore, increasing $\mu$ tends to
stabilize the density of positive pair connections ($b$,$d$,$f$), by
increasing the level of activity of the nodes in the social system. This makes
the QBS & EBS similar to the dynamics of NS where there is no link to prevent
nodes activity.

##  Appendix I Initial condition

To ensure that the initial condition satisfies the relations, we use the
following expression as the initial condition when integrating numerically

$\displaystyle\rho_{a}(0)$ $\displaystyle=$ $\displaystyle
x_{0}^{2}(1-\ell_{0})$ $\displaystyle\rho_{b}(0)$ $\displaystyle=$
$\displaystyle x_{0}^{2}\ell_{0}$ $\displaystyle\rho_{c}(0)$ $\displaystyle=$
$\displaystyle(1-x_{0})^{2}(1-\ell_{0})$ $\displaystyle\rho_{d}(0)$
$\displaystyle=$ $\displaystyle(1-x_{0})^{2}\ell_{0}$
$\displaystyle\rho_{e}(0)$ $\displaystyle=$ $\displaystyle
2x_{0}(1-x_{0})\ell_{0}$ $\displaystyle\rho_{f}(0)$ $\displaystyle=$
$\displaystyle 2x_{0}(1-x_{0})(1-\ell_{0})$ (I.1)

where $\ell_{0}$ is the fraction of friendly links.

## References

  * [1] Thomas Dietz, Elinor Ostrom, and Paul C Stern.  The struggle to govern the commons.  science, 302(5652):1907–1912, 2003.
  * [2] Elinor Ostrom.  The challenge of common-pool resources.  Environment: Science and Policy for Sustainable Development, 50(4):8–21, 2008.
  * [3] Elinor Ostrom.  A general framework for analyzing sustainability of social-ecological systems.  Science, 325(5939):419–422, 2009.
  * [4] Elinor Ostrom.  Governing the commons: The evolution of institutions for collective action.  Cambridge university press, 1990.
  * [5] Elinor Ostrom, Roy Gardner, James Walker, and Jimmy Walker.  Rules, games, and common-pool resources.  University of Michigan press, 1994.
  * [6] Frank Van Laerhoven and Elinor Ostrom.  Traditions and trends in the study of the commons.  International journal of the Commons, 1(1):3–28, 2007.
  * [7] Garrett Hardin.  The tragedy of the commons: the population problem has no technical solution; it requires a fundamental extension in morality.  science, 162(3859):1243–1248, 1968.
  * [8] Garrett Hardin.  Extensions of” the tragedy of the commons”.  Science, 280(5364):682–683, 1998.
  * [9] Oliver P Hauser, David G Rand, Alexander Peysakhovich, and Martin A Nowak.  Cooperating with the future.  Nature, 511(7508):220–223, 2014.
  * [10] Chengyi Tu, Paolo D’Odorico, Zhe Li, and Samir Suweis.  The emergence of cooperation from shared goals in the governance of common-pool resources.  Nature Sustainability, 6(2):139–147, 2023.
  * [11] Billie Lee Turner, Pamela A Matson, James J McCarthy, Robert W Corell, Lindsey Christensen, Noelle Eckley, Grete K Hovelsrud-Broda, Jeanne X Kasperson, Roger E Kasperson, Amy Luers, et al.  Illustrating the coupled human–environment system for vulnerability analysis: three case studies.  Proceedings of the National Academy of Sciences, 100(14):8080–8085, 2003.
  * [12] Richard D Horan, Eli P Fenichel, Kevin LS Drury, and David M Lodge.  Managing ecological thresholds in coupled environmental–human systems.  Proceedings of the National Academy of Sciences, 108(18):7333–7338, 2011.
  * [13] Chris T Bauch, Ram Sigdel, Joe Pharaon, and Madhur Anand.  Early warning signals of regime shifts in coupled human–environment systems.  Proceedings of the National Academy of Sciences, 113(51):14560–14567, 2016.
  * [14] Roger B Myerson et al.  Fundamentals of social choice theory.  Quarterly Journal of Political Science, 8(3):305–337, 2013.
  * [15] Peter Hammerstein and Reinhard Selten.  Game theory and evolutionary biology.  Handbook of game theory with economic applications, 2:929–993, 1994\.
  * [16] Guy Hoffman, Jodi Forlizzi, Shahar Ayal, Aaron Steinfeld, John Antanitis, Guy Hochman, Eric Hochendoner, and Justin Finkenaur.  Robot presence and human honesty: Experimental evidence.  In proceedings of the tenth annual ACM/IEEE international conference on human-robot interaction, pages 181–188, 2015.
  * [17] Timothy N Cason, Daniel Friedman, and Ed Hopkins.  Cycles and instability in a rock–paper–scissors population game: A continuous time experiment.  Review of Economic Studies, 81(1):112–136, 2014.
  * [18] Hendrik Santoso Sugiarto, Lock Yue Chew, Ning Ning Chung, and Choy Heng Lai.  Complex social network interactions in coupled socio-ecological system: Multiple regime shifts and early warning detection.  New Developments in Computational Intelligence and Computer Science (Institute for Natural Sciences and Engineering, Montclair, NJ), 2015\.
  * [19] Joshua S Weitz, Ceyhun Eksin, Keith Paarporn, Sam P Brown, and William C Ratcliff.  An oscillating tragedy of the commons in replicator dynamics with game-environment feedback.  Proceedings of the National Academy of Sciences, 113(47):E7518–E7525, 2016.
  * [20] Yu-Hui Lin and Joshua S Weitz.  Spatial interactions and oscillatory tragedies of the commons.  Physical review letters, 122(14):148102, 2019.
  * [21] Andrew R Tilman, Joshua B Plotkin, and Erol Akçay.  Evolutionary games with environmental feedbacks.  Nature communications, 11(1):915, 2020.
  * [22] NN Chung, LY Chew, and CH Lai.  Influence of network structure on cooperative dynamics in coupled socio-ecological systems.  Europhysics Letters, 104(2):28003, 2013.
  * [23] Maja Schlüter, Alessandro Tavoni, and Simon Levin.  Robustness of norm-driven cooperation in the commons.  Proceedings of the Royal Society B: Biological Sciences, 283(1822):20152431, 2016.
  * [24] David Dubois, Derek D Rucker, and Adam D Galinsky.  Social class, power, and selfishness: when and why upper and lower class individuals behave unethically.  Journal of personality and social psychology, 108(3):436, 2015.
  * [25] Meghdad Saeedian, Maxi San Miguel, and Raul Toral.  Absorbing phase transition in the coupled dynamics of node and link states in random networks.  Scientific reports, 9(1):1–14, 2019.
  * [26] Meghdad Saeedian, Maxi San Miguel, and Raul Toral.  Absorbing-state transition in a coevolution model with node and link states in an adaptive network: Network fragmentation transition at criticality.  New Journal of Physics, 22(11):113001, 2020.
  * [27] Claudio Castellano, Santo Fortunato, and Vittorio Loreto.  Statistical physics of social dynamics.  Reviews of modern physics, 81(2):591, 2009.
  * [28] Tibor Antal, Pavel L Krapivsky, and Sidney Redner.  Dynamics of social balance on networks.  Physical Review E, 72(3):036121, 2005.
  * [29] Stefano Allesina and Si Tang.  Stability criteria for complex ecosystems.  Nature, 483(7388):205–208, 2012.
  * [30] Samir Suweis, Filippo Simini, Jayanth R Banavar, and Amos Maritan.  Emergence of structural and dynamical properties of ecological mutualistic networks.  Nature, 500(7463):449–452, 2013.
  * [31] Ilenia Apicella, Daniel Maria Busiello, Silvia Scarpetta, and Samir Suweis.  Emergence of synchronised and amplified oscillations in neuromorphic networks with long-range interactions.  Neurocomputing, 461:716–726, 2021.
  * [32] Giuseppe Facchetti, Giovanni Iacono, and Claudio Altafini.  Computing global structural balance in large-scale signed social networks.  Proceedings of the National Academy of Sciences, 108(52):20953–20958, 2011.
  * [33] Rajeev Singh, Subinay Dasgupta, and Sitabhra Sinha.  Extreme variability in convergence to structural balance in frustrated dynamical systems.  EPL (Europhysics Letters), 105(1):10003, 2014.
  * [34] A Kargaran and GR Jafari.  Heider and coevolutionary balance: From discrete to continuous phase transition.  Physical Review E, 103(5):052302, 2021.
  * [35] Tuan Minh Pham, Imre Kondor, Rudolf Hanel, and Stefan Thurner.  The effect of social balance on social fragmentation.  Journal of the Royal Society Interface, 17(172):20200752, 2020.
  * [36] M Saeedian, N Azimi-Tafreshi, GR Jafari, and J Kertesz.  Epidemic spreading on evolving signed networks.  Physical Review E, 95(2):022314, 2017.
  * [37] Seth A Marvel, Steven H Strogatz, and Jon M Kleinberg.  Energy landscape of social balance.  Physical review letters, 103(19):198701, 2009.
  * [38] Haye Hinrichsen.  Non-equilibrium critical phenomena and phase transitions into absorbing states.  Advances in physics, 49(7):815–958, 2000.
  * [39] Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  science, 286(5439):509–512, 1999.
  * [40] Duncan J Watts and Steven H Strogatz.  Collective dynamics of ‘small-world’networks.  nature, 393(6684):440–442, 1998.
  * [41] Snehal M Shekatkar and Sukratu Barve.  Importance of initial conditions in the polarization of complex networks.  Europhysics Letters, 122(3):38002, 2018.
  * [42] Marilena Hohmann, Karel Devriendt, and Michele Coscia.  Quantifying ideological polarization on a network using generalized euclidean distance.  Science Advances, 9(9):eabq2044, 2023.
  * [43] Ana-Maria Bliuc, Craig McGarty, Emma F Thomas, Girish Lala, Mariette Berndsen, and RoseAnne Misajon.  Public division about climate change rooted in conflicting socio-political identities.  Nature Climate Change, 5(3):226–229, 2015.
  * [44] Alain Barrat, Marc Barthelemy, and Alessandro Vespignani.  Dynamical processes on complex networks.  Cambridge university press, 2008.
  * [45] Werner Horsthemke and Léon Brenig.  Non-linear fokker-planck equation as an asymptotic representation of the master equation.  Zeitschrift für Physik B Condensed Matter, 27(4):341–348, 1977\.
  * [46] Peter Hänggi.  On derivations and solutions of master equations and asymptotic representations.  Zeitschrift für Physik B Condensed Matter, 30(1):85–95, 1978\.
  * [47] Jose L Garcia-Palacios.  Introduction to the theory of stochastic processes and brownian motion problems.  arXiv preprint cond-mat/0701242, 2007.
  * [48] Crispin W Gardiner et al.  Handbook of stochastic methods, volume 3.  springer Berlin, 1985.