1 Introduction

In the last decades, humanities’ impacts on the environment became increasingly more global in nature, more rapid, extensive, and also threatening for the social and environmental systems themselves [1,2,3]. Scientists refer to this geological era as the Anthropocene due to the impact humanity has and has had on the Earth system [1, 3]. In the Anthropocene, global challenges like climate change and resource scarcity need to be addressed and interventions are required to build a sustainable society which is able to overcome these challenges [4].

To improve the understanding of Earth system dynamics in the Anthropocene, it is crucial to model both the social and the ecological systems as coupled together due to the interconnections and bidirectional feedback processes between humanity and the environment [5,6,7,8,9,10,11], giving rise to the novel field of World-Earth modelling [10]. Examples of such models include network models of individual resource use and social adaptation [12, 13] combined with effects of governance and taxation [14]. Other models investigate land use and social learning [15] or land use and management [16]. There are also models which incorporate not only a social and ecological system but additionally include an economic branch [17].

One important feature of the World-Earth system is that of tipping elements. Tipping elements change rapidly and qualitatively after surpassing a critical threshold, the so-called tipping point [2]. In the past, many natural tipping elements, like the Greenland-ice sheet which can either persist or melt-off completely, due to the melt-elevation feedback, or the Amazon, which does either exist in today’s rainforest state or transition into a Savannah-like state, have been identified and their risk of tipping even within the 2 \(^\circ \)C target of the Paris agreement has been investigated [2, 18].

On the other hand, social tipping and its potential to transform societies, so that the aforementioned (unwanted) natural tipping may be prevented came into research focus recently [19,20,21,22]. Social tipping can for example be observed when a determined minority grows larger than a critical mass and then becomes able to overturn the behaviour or conventions of the majority [23, 24]. Granovetter’s threshold model is an early model formulating and explaining this phenomenon on the example of riots [25]. For the binary decision of not participating or participating in a riot, it is assumed that every individual needs to observe a certain percentage of the population participating in the riot before doing so themselves. This percentage is called their (population-pertaining) threshold. From the distribution of the thresholds, one could then calculate the equilibrium number of participants in the riot. However, real-world threshold distributions are hard to evaluate and reasonable guesses for the distributions lead to everyone or no-one participating in the riot [22]. Two possible extensions to tackle these issues are the introduction of determined minorities and the explanation of an individual’s apparent population-pertaining threshold in terms of its position in the social network and a uniform neighbour-pertaining threshold [22, 26]. This latter type of threshold determines what percentage of a person’s acquaintances need to participate for the person to participate itself. The wide distribution of people’s population-pertaining thresholds then arises, because an emerging cascade of activity reaches different locations in the network at different stages of the riot.

While there is already quite a consolidated understanding of the dynamics within single tipping elements [27, 28], especially from the natural realm [29, 30], it remains an important task of ongoing research to study interactions between different elements, specifically across the climate and the social system [9, 31], since exemplary studies have shown that such interactions can give rise to substantial changes of such systems on the macro-scale [32, 33]. Recent work has investigated also the interaction of tipping elements in natural [34,35,36,37,38], social [19], socio-economic [39], or socio-ecological systems [40]. However, most of the aforementioned studies have either studied tipping elements using identical prototypical normal forms for each element [35,36,37,38, 40] or by focussing on cascading effects and interactions within either only natural or social subsystems [19, 34,35,36,37,38] and thus without explicitly accounting for their potential cross-system interactions.

To bridge this gap, the present work investigates potentially arising dynamics of coupled tipping elements across the natural and social system and proposes a conceptual, low-dimensional model consisting of a stylised environment being described by one macroscopic variable which we call the pollution and a social network of individuals which change their behaviour, i.e., either their contribution to the pollution or their active decision to not contribute to it, according to change probability rates based on three threshold processes:

  1. 1.

    Direct environmental impacts like droughts, heatwaves, or air pollution can threaten human society [41,42,43]. An adaptation to and mitigation of these impacts is necessary in a changing Earth system [44].

  2. 2.

    Social contagion which is recognised as the spread of behavioural patterns among many individuals in an interacting society [45], e.g., social movements [46,47,48] or the spread of social norms [49].

  3. 3.

    Anticipated environmental impacts are rarely discussed in World-Earth models where social system’s are mostly coupled to the environment through precedent or contemporary impacts. However, giving societies the capability of anticipation and foresight might alter the social behavioural patterns significantly. For example, the threat of global climate changing to the worse, like in a “Hothouse Earth” scenario [50], could spark and strengthen climate movements with the goal to prevent consequential environmental impacts [51]. Especially, information and the anticipatory time horizon up to which an individual projects such impacts can have significant influence on the behaviour of individuals and the macroscopic dynamics [11, 52].

In particular, we incorporate the social system via an agent-based model (ABM) similar to Ref. [22], since ABMs have been proven as a promising tool to model social behaviour and decision-making [53,54,55]. In ABMs, individuals, called agents, adjust their behaviour according to certain deterministic or stochastic rules after assessing their individual situation [54]. One of the strengths of ABMs lies in the microscopic details which arise from the simulation of individuals and can induce emergent macroscopic phenomena [56]. Within ABMs, non-linear features, such as local heterogeneous structures, can be observed and analysed to help understanding the impact of disturbances or causes of unexpected outcomes [7]. However, ABMs come with significant computational costs and are hard to analyse formally. Therefore, we also apply a mean-field approximation of the microscopic behaviour of the network of interacting agents.

The model simulations presented in this work reveal that a direct coupling between a dynamic environment and social dynamics leads to one globally stable attractor. However, additionally accounting explicitly for the effects of anticipation induces metastability where additional environmental states that are not within or close to that stable attractor persist over a long time.

The remainder of this paper is structured as follows. First, in Sect. 2, the micro-level socio-ecological model itself is introduced and the implementation of the three aforementioned processes is explained. Second, a mean-field approximation of the model is given in Sect. 2.3. Third, we give an overview of the model behaviour, first for a static and then a dynamic environment, and compare it to the mean-field approximation in Sect. 3.1. In Sect. 3.2, one can find the main results of the paper, metastability for high anticipation times resulting in a stabilisation of potentially unpolluted environments and thus social tipping. We close this paper with a discussion of the results, possible connections to real-world parameters, and a breif outlook in Sect. 4.

2 Model

To illustrate the influence of the environment on social activation and vice versa, we formulate a conceptual network model which is tuneable from a mainly social dynamic version similar to Granovetter’s threshold model [25] as in Ref. [22] to a coupled socio-ecological agent-based model. A focus in the model is the effect of anticipation. In Fig. 1, we provide a schematic overview of the model’s setup which illustrates the different interaction mechanisms between the environment and the social system, explained in detail below.

Fig. 1
figure 1

Schematics for the model setup. The environment is represented by just one variable, the global pollution level Y, which changes proportional to its own value and the social systems mean state, the share of polluting agents X. The social system consists of individual agents which are either adding to the pollution or not. They change their binary state \(x_i\) depending on the current pollution level (direct impacts), the pollution level they anticipate in the future (anticipated impacts) and the behaviour of their direct neighbours in the network (social contagion)

Full size image

First, we define agent i’s pollution state as:

$$\begin{aligned} x_i = {\left\{ \begin{array}{ll} 1,\quad \text {for }\text {polluting} \text { (i.e., not mitigating pollution) agents}\\ 0,\quad \text {for }\text {non-polluting}\text { (i.e., mitigating pollution) agents} \end{array}\right. }. \end{aligned}$$
(1)

Agents are able to switch from one state to the other, and especially, the direction from polluting to actively non-polluting is of importance if one wants to guide the system towards an unpolluted state. The mean share of polluting agents will be denoted \(X= \frac{1}{N}\sum _i x_i\) where N is the number of agents. We use this macroscopic variable X to couple the social system to the environment which is represented by only one macroscopic variable, the pollution Y.

2.1 Ecological dynamics

As greenhouse gases, especially CO\(_2\), are the major driving force of anthropogenic climate change [57, 58], we motivate our environmental module with their fundamental dynamics. Carbon added to the atmosphere, for example by burning fossil fuel, does not stay in the atmosphere, but is constantly removed through natural carbon sinks like terrestrial ecosystems or the ocean, with the latter being the major carbon sink in the Earth system [59, 60]. For simplicity, the carbon decay is modelled by a linear differential equation with a single average lifetime \(\tau \) [59, 61, 62]. Additionally, pollution is added proportionally to the share of polluting agents X:

$$\begin{aligned} \dot{Y} = \frac{1}{\tau } (X-Y). \end{aligned}$$
(2)

In particular, Y is measured in units of maximal equilibrium pollution, so that if all N agents pollute consistently (\(X=1\)), then \(Y\rightarrow 1\). The average lifetime of pollution, \(\tau > 0\), allows to scale the system to a quasi-static environment by letting \(\tau \rightarrow \infty \). Here, the term quasi-static means that the environmental dynamics become very slow in comparison to the dynamics of the social system for which we will explain all further details below in Sect. 2.2. However, a reasonable estimate for a real-world average lifetime of atmospheric carbon is approximately 50 years [60, 61]. Consequently, for \(\tau =1\), one time unit in our model roughly corresponds to 50 years in a real-world representation of the Earth’s atmosphere. We are aware that this implementation of the environment is very conceptual and that there are also comparatively more advanced models of the carbon cycle [59, 63,64,65]. However, for the qualitative behaviour of our model only some form of natural decay is necessary, so the environmental dynamics from Eq. (2) suffices.

2.2 Social dynamics

The social dynamics of the share of polluting agents X arise from the dynamics of the single agents’ states \(\mathbf{x} = (x_i)_{i=1}^N\). The agents change their states with the change probability rates \(p_i^-(\mathbf{x},Y)\) from polluting to non-polluting (\(x_i=1\rightarrow x_i=0\)) and \(p_i^+(\mathbf{x},Y)\) from non-polluting to polluting (\(x_i=0\rightarrow x_i=1\)):

$$\begin{aligned} p_i^-(\mathbf{x},Y)&= \alpha \cdot p_\text {dir}^-(Y) + \beta \cdot p_{i\text {,soc}}^-(\mathbf{x})\cdot p_\text {ant}^-(X,Y),\end{aligned}$$
(3)
$$\begin{aligned} p_i^+(\mathbf{x},Y{)}&=\alpha \cdot p_\text {dir}^+(Y) + \beta \cdot p_{i\text {,soc}}^+(\mathbf{x})\cdot p_\text {ant}^+(X,Y). \end{aligned}$$
(4)

There are three contributions to the change probability rates which we will motivate and explain in detail in the following paragraphs: those due to direct environmental impacts, \(p_\text {dir}^\pm (Y)\), due to social contagion, \(p_{i,\text {soc}}^\pm (\mathbf{x})\), and due to anticipated environmental impacts, \(p_\text {ant}^\pm (X,Y)\). The first contribution is scaled by a parameter \(\alpha \ge 0\) which we call the vulnerability, since one could interpret it as how susceptible an agent is to immediate changes in their environment. The second and third contributions are coupled together and are scaled by a parameter \(\beta \ge 0\) which we call the farsightedness, since it weights the anticipation effect. This coupling can be understood as social contagion leading to a spread of the action due to anticipated environmental impacts. Note that among these factors, only \(p_{i,\text {soc}}^\pm (\mathbf{x})\) depends on the microstate \(\mathbf{x}\), while the others only depend on the macrostate (X, Y).

2.2.1 Direct environmental impacts

Direct environmental impacts like droughts or storms have triggered sudden societal changes in the past, including the onset of migratory patterns or even the collapse of whole nations [66, 67]. In our model, we link the behaviour of agents to the current pollution level Y through a threshold function. If the pollution is above the pollution threshold parameter \(\gamma \in [0,1]\), a polluting agent has a higher probability to become non-polluting and vice versa:

$$\begin{aligned} p_\text {dir}^-(Y)&= H\left( Y-\gamma \right) , \end{aligned}$$
(5)
$$\begin{aligned} p_\text {dir}^+(Y)&= 1-H\left( Y-\gamma \right) , \end{aligned}$$
(6)

where H is the Heaviside step function.

2.2.2 Social contagion

The embedding of individuals in a social network and the observation of their immediate surroundings has been proven to have a large influence on people’s behaviour [45, 49, 68, 69]. Recently proposed network models allow to incorporate and explain such dynamics in the context of opinion formation [70], behaviour [68, 71] and information spread [72] or social movements and collective action [69]. Granovetter’s threshold model [25] and recent network-based extensions [22, 73,74,75] describe such influence by other individuals or direct neighbours, also known as social contagion [76]. From Refs. [22, 26], we take the idea to define a microscopic neighbour-pertaining threshold \(\chi \in [0,1]\) which influences the likelihood of the agents’ behavioural change. If the share of polluting neighbours in the fixed social network G, \(\frac{X_i}{k_i}\), is lower than the social threshold \(\chi \), a polluting agent becomes more likely to convert to being non-polluting and vice versa. Here, \(X_i = |\{j:(i,j)\in E,~x_j=1\}|\) gives the number of polluting neighbours, \(k_i = |\{j:(i,j)\in E\}|\) is i’s degree (number of neighbours), and E is the fixed set of network edges:

$$\begin{aligned} p_{i,\text {soc}}^-(\mathbf{x})&= 1-H\left( \frac{X_i}{k_i}-\chi \right) , \end{aligned}$$
(7)
$$\begin{aligned} p_{i,\text {soc}}^+(\mathbf{x})&= H\left( \frac{X_i}{k_i}-\chi \right) . \end{aligned}$$
(8)

2.2.3 Anticipated environmental impacts

The anticipation of potential environmental catastrophes might spark, strengthen, and justify climate movements like the recent Fridays for Future-movement [51]. Specifically, the anticipation time, i.e., how far one extrapolates the trajectory of the Earth system into the future, can have a big impact on the social urge to act [11, 52]. We assume the simplest possible form of anticipation, a linear extrapolation, and therefore define the anticipated pollution as:

$$\begin{aligned} Y_\text {ant}=Y+\theta {\dot{Y}}, \end{aligned}$$
(9)

where \(\theta > 0\) is the anticipation time. If the anticipated pollution is above the pollution threshold \(\gamma \), a polluting agent becomes more likely to get non-polluting. Likewise, if it is below, a non-polluting agent becomes more likely to get polluting:

$$\begin{aligned} p_\text {ant}^-(X,Y)&= H\left( Y + \theta {\dot{Y}} - \gamma \right) \nonumber \\&= H\left( Y + \frac{\theta }{\tau } (X-Y) - \gamma \right) , \end{aligned}$$
(10)
$$\begin{aligned} p_\text {ant}^+(X,Y)&= 1 - H\left( Y + \theta {\dot{Y}} - \gamma \right) \nonumber \\&= 1 - H\left( Y + \frac{\theta }{\tau } (X-Y) - \gamma \right) . \end{aligned}$$
(11)

Plugging these factors into Eqs. (3, 4) yields the full change probability rates:

$$\begin{aligned} p^-_i(\mathbf{x},Y)&= \alpha H(Y - \gamma ) + \beta \left[ 1 - H(X_i - \chi k_i)\right] \nonumber \\&\quad \times H\left( Y + \frac{\theta }{\tau } (X-Y) - \gamma \right) , \end{aligned}$$
(12)
$$\begin{aligned} p^+_i(\mathbf{x},Y)&= \alpha \left[ 1 - H(Y - \gamma )\right] + \beta H(X_i - \chi k_i) \nonumber \\&\quad \times \left[ 1 - H\left( Y + \frac{\theta }{\tau } (X-Y) - \gamma \right) \right] , \end{aligned}$$
(13)

which together with Eq. (2) fully define our system for a given network structure.

2.3 Mean-field approximation

Applying a mean-field approximation and therefore eliminating the explicit network dependence allows us to get important insights into the systems dominant dynamics. We therefore assume that the network is large, densely connected, and has a small diameter, such that changes spread fast through the network and thus, most of the time, the share of polluting neighbours is approximately equal for all nodes. Consequently, the latter can be approximated well by the overall share of polluting agents:

$$\begin{aligned} \frac{X_i}{k_i} \approx X. \end{aligned}$$
(14)

Plugging this into Eqs. (7, 8) cancels out the network dependence and all agents have approximately the same change probability rates:

$$\begin{aligned} p^-(\mathbf{x},Y)&\approx {{\tilde{p}}}^-(X,Y) = \alpha H(Y- \gamma ) + \beta \left[ 1 - H(X - \chi )\right] \nonumber \\&\quad \times H\left( Y + \theta \dot{Y} - \gamma \right) , \end{aligned}$$
(15)
$$\begin{aligned} p^+(\mathbf{x},Y)&\approx {{\tilde{p}}}^+(X,Y) = \alpha \left[ 1 - H(Y-\gamma )\right] + \beta H(X - \chi )\nonumber \\&\quad \times \left[ 1 - H\left( Y + \theta \dot{Y} - \gamma \right) \right] . \end{aligned}$$
(16)

Since the network is also assumed to be very large, the actual stochastic change in X can be approximated by its expectation value and we obtain finally the deterministic differential equations:

$$\begin{aligned} \dot{X}&= (1 - X ) {{\tilde{p}}}^+(X,Y) - X {{\tilde{p}}}^-(X,Y),\nonumber \\ \dot{Y}&= \frac{1}{\tau } (X - Y). \end{aligned}$$
(17)

Through this approximation, we reduced the number of equations from \(2N+1\) to only two (however non-smooth) differential equations which can be solved numerically. In addition, a piece-wise analytical solution is possible, since the change rates \(p^\pm \) are constant except when X hits the threshold \(\chi \) or when Y or \(Y_\text {ant}\) hits the threshold \(\gamma \). For constant \(p^\pm \), the differential equations become linear in X and Y and can be solved through an exponential ansatz and the method of varying constants (see S1). However, due to the non-smoothness, it is not possible to find a full solution, but only such a piece-wise solution. For the numerical implementation, we used the piece-wise calculated equations for X and Y and probed through integrating for discrete time steps which particular piece-wise solution has to be used each time t.Footnote 1

3 Results

In the first part of this section, we compare the dynamics of the microscopic model to the dynamics of the above introduced mean-field approximation. First, the case of a static environment is studied, i.e., the social dynamics are significantly faster than the environmental dynamics (\(\alpha \,,\beta =\text {const.},\;\tau \rightarrow \infty \)), to show that our model displays a bistable behaviour as also found in other works on social tipping [22, 74]. Second, we study the case of a dynamic environment, i.e., the social system evolves on a similar timescale as the environmental system, and we investigate the attractors of this feedback system. Subsequently, we use the mean-field approximation and identify multi-stability in the regime of long anticipation times.

3.1 Comparing the microscopic model to the mean-field approximation

Fig. 2
figure 2

Comparison of the microscopic model (a, b, e, f) to the mean-field approximation (c, d, g, h). a–d Static environment with parameters \(\gamma =0.6\), \(\chi = 0.4 \), \(\alpha = 0 \), \(\beta = 10\), \(\tau = 10^6\), \(\theta = 1\), and three different initial conditions (dashed, dotted, and dash-dotted). In (a, c), the distribution of trajectories \(p\left( X(t=20)\right) \) with \(X(t=20)\) is illustrated through a histogram (initial conditions \((X_0,Y_0)\) are taken from the unit interval \(X_0,Y_0\in [0,1]\) with a step size 0.05 each). Generally, two main attractors at \(X=0\) and \(X=1\) can be seen. Panels (b, d) show trajectories with \(Y(t=0)=1\) and three different \(X(t=0)\) for the microscopic simulations and the mean-field approximation, respectively. The black lines indicate the thresholds \(\gamma \) and \(\chi \). e–h Dynamical environment with \(\tau =1\), a direct coupling of the social system to the environment, \(\alpha =0.1\), and other parameters as before. There is only one main attractor at \(X(t=20)\approx \gamma \) for both the microscopic simulations (e, f) and the mean-field approximation (g, h). Deviations can be found in the transient behaviour (\(t<2\)) as seen in the trajectories, as shown in (f, h), while the long-term behaviour is similar. In the insets, oscillations around the pollution threshold \(\gamma \) are shown for one trajectory (\(X(t=0)=0.75,\,Y(t=0)=1\))

Full size image

We compare the mean-field approximation with microscopic model simulations for Erdős–Rényi networks with \(N=200\) nodes and an average degree of \({\bar{k}}=10\). Calculations are robust against larger network size, so that finite-size effects are negligible in the presented results (not shown). In Fig. 2a–d, results for a static environment (\(\tau = 10^6\)) are displayed. Additionally, the social system does not directly respond to the environment (\(\alpha =0\)) but only to the anticipated impacts (\(\beta =10\)). These parameter choices ensure that the dynamics are dominantly driven by social contagion, see Eqs. (3, 4), so that we can compare our model in this boundary case against other studies [22, 74]. The microscopic results are shown in Fig. 2a, b and the mean-field approximation results in Fig. 2c, d. In Fig. 2a, c, the distribution of the long-term share of polluters, \(p(X(t=20))\), is shown as a histogram. In both cases, we observe two prominent maxima at \(X=0\) and \(X=1\) which compares well to previous network-based threshold models that are showing bistability, as well [22, 74]. The exemplary trajectories in Fig. 2b, d show that the trajectories approach \(X=0\) if \(Y\geqslant \gamma \) and \(X<\chi \). In the microscopic case, Fig. 2b, trajectories with \(X(t=0)\) slightly above \(\chi \) can approach the attractor at \(X=0\), even though the social system would be above the social threshold. This can be explained through the local heterogeneity of the network leading to some nodes in an environment with \(\frac{X_i}{k_i}<\chi \), even though \(X\geqslant \chi \). These nodes can then start a cascade eventually leading to all nodes approaching the same state. In the mean-field approximation, Fig. 2d, such a behaviour cannot be observed as \(X<\chi \) and \(Y+\theta \dot{Y}\geqslant \gamma \) results in all change rates becoming zero. To account for these mismatches, future work should explore more advanced approximation techniques, such as heterogeneous mean-field approximations [77] to provide a better representation of the model’s dynamics for those limiting cases as well. Since most other trajectories and the distributions of steady states, however, match well across the numerical experiments and the analytical estimations we consider our approximation to perform sufficiently well for the purpose of the present study.

Fig. 3
figure 3

a A metastability diagram with respect to anticipation time \(\theta \) is shown for model parameters of \(\alpha =0.1\), \(\beta =10\), \(\tau = 1\), \(\gamma =0.6\), and \(\chi = 0.4\). Every column corresponds to a histogram obtained for one \(\theta \) representing the occurrences of share of polluters \(p(X(t=20))\). The black dashed lines indicate the values of \(\theta \) in panels (b–d). b Phase space diagram with two exemplary trajectories (blue and red lines) for \(\theta =5\) and other parameters as before. Line colours darken with time t. D gives the diagonal where \(X=Y\) and S the line where \(Y_\text {ant}=\gamma \). R, L, and H indicate locations where the trajectories show discontinuous transitions. A more-detailed discussion can be found in the Supplement (see S2 and Fig. S.1). c, d Polluting agents \(X(t=20)\) as a function of the initial condition plotted as a colour code for \(\theta =20\) (c) and \(\theta =100\) (d)

Full size image

In Fig. 2e–h, we introduce a dynamic environment (\(\tau =1\)) and couple the social system directly to the environment (\(\alpha =0.1\)). The system again reaches its attractor before \(t=20\), and in Fig. 2e, g, the distribution of \(X(t=20)\) is displayed. In contrast to the static environment case, it peaks around \(X\approx Y\approx \gamma \) for both the microscopic simulations (Fig. 2e) and the mean-field calculations (Fig. 2g). The existence of only one main attractor around \(X\approx Y\approx \gamma \) is emphasized through the trajectories in Fig. 2f, h which in the long-term show the same behaviour for microscopic and mean-field simulations. However, in the transient phase \((t<2)\), the microscopic simulations still differ from the mean-field solution due to the heterogeneity of the underlying network. The coupling to a dynamic environment, Fig. 2e–h, leads to a loss of the multi-stability of the pure social model shown in Fig. 2a–d which can be interpreted as a form of “damping”. Additionally, one can observe that \(|X-Y|\) becomes very small which should be expected, since due to Eq. (2), the pollution always follows the value of X. However, X and Y do not completely converge, but oscillate around the pollution threshold \(\gamma \) which is illustrated in the insets in Fig. 2f, h for one illustrative instance of X(t) (\(X(t=0)=0.75\)). During each period of this oscillation, both Y and \(Y_\text {ant}\) cross the threshold \(\gamma \) twice, so the trajectory is continuous but not smooth. Due to the stochasticity of the microscopic simulations, the oscillations in Fig. 2f are less regular than the oscillations in the mean-field approximation, see the inset in Fig. 2h.

We saw that the mean-field approximation is not accurate in the transient behaviour, see Fig. 2f, h for \(t<2\), but generally provides comparable results in the long-term behaviour. Therefore, it qualifies for a qualitative assessment of the model dynamics which we will use in the next section to analyse the trajectories for different anticipation times \(\theta \). In a similar spirit as above, we suggest to explore the potential for higher order mean-field approximations [77] to further improve the representation of the model’s dynamics in the transient phase as well.

3.2 Metastability for high anticipation time

The introduction of anticipation is one of the key features of our model and therefore needs special evaluation. Thus, we show in Fig. 3a the long-term share of polluters \(p(X(t=20))\) resolved for varying anticipation times and model parameters \(\alpha =0.1\), \(\beta =10\), and \(\tau =1\) as a histogram. For \(\alpha >\beta \), the dynamics are dominated by direct environmental impacts and the system converges towards one attractor at the pollution threshold \(\gamma \) similar to Fig. 2e–h. The choices of \(\alpha =0.1\) and \(\beta =10\) ensure that we can observe social dynamics which are driven dominantly by social contagion processes instead of behavioural changes due to direct environmental impacts. Our choice of \(\tau \) allows for the dynamic evolution of the environment on similar timescales as the social dynamics, compare Fig. 2e–h.

As our main focus is the qualitative model behaviour, we use the mean-field approximation given by Eq. (17) which leads to an overall good agreement to the numerical experiments (compare Fig. 2). In Fig. 3, the results for the dynamics with different anticipation times \(\theta \) are displayed. On the left-hand side of Fig. 3a, the existance of only one attractor at \(X=Y=\gamma =0.6\) is clear. This reproduces for small \(\theta \) the results from the previous section, compare Fig. 2g. In Fig. 3b, a phase space diagram is shown for two exemplary trajectories (red and blue lines) in which the trajectories approach the line where the anticipated pollution is equal to \(Y_\text {ant} =Y+\frac{\theta }{\tau }(X-Y) =\gamma \) (dotted line S in Fig. 3b) and then approach the attractor at \(X=Y=\gamma \) through a large number of small steps. However, between \(\chi \) and \(\gamma \), the trajectories can deviate from S. In the Supplement (see S2), we show that, indeed, \(X=Y=\gamma \) is a stable fixed point of the deterministic macroscopic approximation whenever \(\chi<\gamma <\beta /(\alpha +\beta )\) (see lines L, H, and R in Fig. 3b, respectively) and \(\theta >\tau \).

Despite this unique attractor that determines the long-term fate of the system, two kinds of long metastable behaviour emerge if the anticipation time \(\theta \) is large enough. In the snapshot at time \(t=20\) shown in Fig. 3a as a function of \(\theta \), one can see that first an upper “branch” appears for \(\theta \gtrsim 10\), representing trajectories that at some point visit a region of both high X and Y and then decrease very slowly towards \(\gamma \), like the red one, as shown in Fig. 3b. More interestingly, also a lower branch becomes visible for \(\theta \gtrsim 20\), representing trajectories that at some point visit a region of both low X and Y, then increase slowly towards \(X=\chi \) and \(Y_\text {ant}=\gamma \) before then suddenly moving fast to \(X,Y>\gamma \) and eventually converging to the attractor, like the blue one shown in Fig. 3b. The exact dependence of the system’s position at \(t=20\) on its initial condition is shown in Fig. 3c, d for \(\theta =20\) and \(\theta =100\), red and blue areas representing \(X(t=20)>\gamma \) and \(X(t=20)<\gamma \), respectively.

The metastability in both branches arises because of two effects. First, close to line S, the sign of \(\dot{X}\) alternates fast (as seen in the insets in Figs. 2f, h, 3b), keeping the system close to that line. Second, the larger \(\theta \), the closer line S is to the diagonal D, causing \(|\dot{Y}|\) to become very small when on S. Since in the macroscopic approximation \(\dot{Y}\) has a constant sign throughout such a phase, it slowly converges to either \(X=Y=\gamma \) (red) or \(X=\chi \), \(Y_\text {ant}=\gamma \) (blue). In the latter case, once the system crosses \(X=\chi \) (line L), it moves fast towards the upper right. In addition to these long-term effects that can be well explained by the macroscopic approximation the micro-model exhibits even longer metastable phases when \(\theta \) becomes large. In that case, the lines S and D are close enough, so that the finite-size changes of the micro-model can move the system across D while close to S. As a result, not only the sign of \(\dot{X}\) but also that of \(\dot{Y}\) will alternate. In the limit of \(\theta \rightarrow \infty \), S and D coincide and the system may stay forever in any point on S with \(X<\chi \) or \(\gamma<X<\beta /(\alpha +\beta )\) (see S3 in the Supplement for details).

4 Discussion and conclusion

We have proposed a conceptual co-evolutionary model of social contagion, environmental dynamics, and associated immediate and anticipated impacts that serve as a mediator between the two subsystems. The coupling of the social system to a dynamic environment has led to a dampening of the social dynamics onto one attractor. The influence of anticipation on the dynamics has opened the possibility to guide the system towards a region with low shares of polluting agents X.

In Sect. 3.1, we recovered that the pure social dynamics can be multi-stable if the environment is static. However, introducing a direct coupling to a dynamic environment dampened the multi-stability and one main attractor formed exactly at the pollution threshold \(\gamma \). This can be explained by the general tendency of the agents to lower their environmental impact if the pollution is high but also reduce their mitigation cost if the pollution is low. If the agents act short-sightedly and adjust their behaviour directly with regard to the current pollution level, this leads to a fast convergence towards a trade-off scenario between environmental and “economic” benefits. Thus, if the pollution threshold were to represent the level of pollution the agents really consider “optimal”, the short-sighted behaviour would be perfectly rational, since it reaches and stabilises the optimum as fast as possible. In our more “boundedly rational” interpretation of the model, however, the pollution threshold does not represent the optimal level of pollution. Instead, it is assumed to be much lower, but represents the level at which the problem of pollution is perceived as so urgent by the agent that immediate action seems due. For that case, we considered an additional farsighted driver for behaviour that is, however, modulated by social dynamics: agents may stop polluting already when they anticipate pollution to cross the pollution threshold at a certain future time and enough of their neighbours are non-polluting already. If this additional neighbour-pertaining threshold \(\chi \) is well below the pollution threshold and the anticipation time is high, this has the effect that trajectories starting at low pollution levels become metastable, staying for very long time close to the (lower) neighbour-pertaining threshold. Afterwards, it suddenly increases to a level above the pollution threshold. From here, the trajectories eventually converge to the threshold from above (blue line in Fig. 3b). The downside of this encouraging result is that now trajectories starting at high pollution levels also become metastable above the threshold and decrease towards it only very slowly without an intermediate low pollution phase (red line in Fig. 3b). This metastable regime seems to form rather nonlinearly when \(\theta \) is increased. Even though Fig. 3a reminds of a bifurcation diagram, there is no actual bifurcation: \(X=Y=\gamma \) remains the only attractor of the deterministic macroscopic approximation throughout. The stochastic microscopic model, however, can in principle visit and stay in the region around \(X=Y=\chi \) frequently due to its finite-size stochastic jumps. Thus, we argue to find a process that one may denote as social tipping, even though it is not possible to define only one critical anticipation time above which the metastable behaviour predominantly exists. We conclude that the larger the anticipation time, the longer we expect the system to remain in the metastable regime before potentially (at least in the mean-field approximation) converging back to the stable attractor. The desired time that the system should stay within the transient metastable regime may very well vary across different real-world applications. Hence, our results illustrate that social tipping can be interpreted as more than simply a parameter-induced (saddle-node) bifurcation. Instead, our model results align with recent definitions of social tipping that specifically account for the existence of qualitatively alternative (meta-)stable states whose existence can not solely be explained by a single control parameter [20, 78].

If one wants to decrease the long-term average pollution level of the modelled system by tweaking its parameters or directly influencing its state, one could consider the following strategy. Whenever the system is at \(X>\chi \) (more agents polluting than the neighbour-pertaining threshold), it first has to be moved to \(X<\chi \). In the real world, this could be achieved by temporary measures such as taxation or punishment for pollution or subsidies for environmentally friendly behaviours. Note that a mere reduction of X below \(\gamma \) but not below \(\chi \), or a reduction of Y instead of X would not suffice, since then X and Y would grow again due to the social dynamics. If the anticipation time \(\theta \) can be influenced, then once the system is at \(X<\chi \), \(\theta \) could be raised to increase the length of the resulting metastable phase with \(X<\chi \). How this could be achieved in reality is less clear, since the anticipation time is strongly linked to how much agents care about the future. If the neighbour-pertaining threshold \(\chi \) can be influenced, a trade-off emerges between (i) lowering \(\chi \) to achieve lower pollution levels during the metastable phase, and (ii) keeping \(\chi \) large enough, so that X can always be brought back to \(X<\chi \) by political measures once a metastable phase has ended. Note that such political measures seem not essential during a metastable phase but may of course help to prolong it.

Even though our conceptual model is relatively simple, it provides important and surprising insights into the socio-ecological dynamics with anticipation as a key feature. However, it offers room for adjustments, e.g., the implementation of more realistic environmental dynamics which themselves could be subject to tipping such as previously described for ice dynamics [30, 36].

Studying the effect of logistically growing change probabilities can alter a systems overall behaviour and lead to novel trajectories [79, 80] which would be a promising addition to our model, as well. In general, studying other effects in addition to social activation, such as adaption, which has been shown to play a crucial role for the stability of co-evolutionary socio-ecological systems [12, 14, 81, 82], is a promising step to broaden the application of our model and strengthen its meaningfulness. An additional step in this direction would be the discussion of other potentially more realistic network topologies, such as scale-free or small-world networks, in their effect on the simulated dynamics. Along those lines, advanced approximation techniques that are specifically tailored to such network topologies, e.g., heterogeneous mean-field approximations [77], could thereby further strengthen the overlap between numerical results and analytical calculation and at the same time explain some of the mismatches that we have observed in our present work. In summary, our work provides a valuable first step in investigating potentially emergent dynamics from interactions between individuals in a dynamic and accessible environment.