Setting up the physical principles of resilience in a model of the Earth System

Resilience is a property of social, ecological, social-ecological and biophysical systems. It describes the capacity of a system to cope with, adapt to and innovate in response to a changing surrounding. Given the current climate change crisis, ensuring conditions for a sustainable future for the habitability on the planet is fundamentally dependent on Earth System (ES) resilience. It is thus particularly relevant to establish a model that captures and frames resilience of the ES, most particularly in physical terms that can be influenced by human policy¹. In this work we propose that resilience can serve as a theoretical foundation when unpacking and describing metastable states of equilibrium and energy dissipation in any dynamic description of the variables that characterise the ES. Since the impact of the human activities can be suitably gauged by the planetary boundaries (PBs) and the planet's temperature is the net result of the multiple PB variables, such as CO₂ concentration and radiative forcing, atmospheric aerosol loading, atmospheric ozone depletion, etc, then resilience features arise once conditions to avoid an ES runaway to a state where the average temperature is much higher than the current one. Our model shows that this runaway can be prevented by the presence of metastable states and dynamic friction built out of the interaction among the PB variables once suitable conditions are satisfied. In this work these conditions are specified. As humanity moves away from Holocene conditions, we argue that resilience features arising from metastable states might be crucial for the ES to follow sustainable trajectories in the Anthropocene that prevent it run into a much hotter potential equilibrium state.

1 Introduction

Over the past decades the human imprint on the Earth System (ES) has been exceptional (Steffen et al., 2015a; Jouffray et al., 2020). While the mass of humans is only about 0.01 % of the total biomass, we have become a dominant force in shaping the face of Earth, including its atmosphere, biosphere, hydrosphere and lithosphere (Vitousek et al., 1997; Crutzen, 2002; Ellis, 2011; Foley, 2011; Nyström et al., 2019), and as of 2020 the global human-made mass surpasses the dry-weight of all living biomass (Elhacham et al., 2020). Thus, humans have become a hyper-keystone species (Worm and Paine, 2016), which rivals geological forces in influencing the trajectory of the ES (Steffen et al., 2018)

A major concern of these changes is the risk of crossing of so-called tipping-points, which refer to the critical threshold at which a small change or event triggers a significant and potentially irreversible (regime) shift in a system (Lenton et al., 2008). Tipping-points have been observed in various systems, such as ecosystems (e.g. food webs, benthic communities), social systems (e.g. norms, policy), economic systems (e.g. market-based economy) and technological systems (e.g. steam engine, smartphone, artificial intelligence) (Scheffer et al., 2001; Scheffer, 2009; Nyborg et al., 2016). Over the past couple of decades there have been raising concerns around the existence of tipping-elements, which are large-scale components (subsystems) of the ES that may transgress a tipping-point (Lenton et al., 2008; Barnosky et al., 2012). Example of such tipping-elements include, the Greenland Ice Sheet, the Atlantic Meridional Overturning Circulation (AMOC), permafrost, monsoon systems, and the Amazon rainforest. Importantly, these tipping-elements interact, which may lead to a cascading behaviour of the entire ES (Wunderling et al., 2024). The consequences of these dynamics for humanity could be colossal (Steffen et al., 2018).

Clearly, knowledge about tipping-points, where they are located, when they are approached and identifying ways to navigate away from them, are key challenges for humanity (Barnosky et al., 2012; Scheffer et al., 2012). Two broad frameworks that could help assist in this regard are planetary boundaries and resilience theory. The two are complementary in the sense that the planetary boundaries provide a quantitative assessment whereas the resilience framework adds a strong theoretical underpinning.

The planetary boundaries (PB) framework (Rockström et al., 2009; Steffen et al., 2015b; Richardson et al., 2023) has been used to define global and regional limits in biophysical processes – “safe operating space” – that must not be crossed if humanity is to stay away from systemic and potentially irreversible shifts in the ES. As such, the planetary boundaries framework serves as a “global dashboard,” tracking humanity's collective impact on key environmental factors that threaten the Earth's ability to sustain human life. More recently, focus has been directed towards exploration of how different boundaries can interact and potentially cascade, thereby shrinking the safe operating space for future human impacts on the ES (Lade et al., 2020). Importantly, while the PB framework highlights the presence of tipping points in biophysical processes, it does not specify their exact thresholds. Instead, it delineates two risk zones: a zone of increasing risk and a high-risk zone. In the former, the further boundary limits are exceeded, the greater the likelihood of causing significant harm–destabilizing critical Earth system processes and disrupting essential life-support functions. In the latter, or high-risk zone, there is a substantial risk of severe and potentially irreversible damage to key planetary functions. In essence, these zones are defined at a precautionary distance from the estimated locations of potential tipping points.

The resilience concept describes the extent to which a system can resist and develop (e.g. ecosystems or the the entire ES) with change by absorbing recurrent perturbations, deal with uncertainty and risk, and still sustain its key properties (Holling, 2001; Folke, 2006). This conception of resilience is based on the understanding that humans and Nature are deeply interconnected through feedbacks between social and ecological components, which together influence overall behavior and dynamics (Biggs et al., 2012). This interdependence defines a social-ecological system (Berkes and Folke, 1998) in which human well-being and prosperity rely on the stability and functioning of the Earth system (Folke et al., 2011). Multiple states (regimes), tipping-points and self-reinforcing feedback mechanisms (hysteresis) are a central feature of resilience (Holling, 2001). In cases where resilience is high, a powerful shock – such as, storms, large wildfires, pest outbreaks in ecosystems, or armed conflicts, trade wars, and supply chain disruptions in social systems – is required to push the system beyond a tipping-point and into another state. However, gradual (creeping) change – such as, loss biodiversity, habitat fragmentation and pesticide resistance in ecosystems, or growing inequality and changing social norms in society – erodes resilience of the current state. This makes the system vulnerable even to smaller perturbations. Once the system finds itself in this new state it can be difficult, or even impossible to reverse due to self-reinforcing feedback mechanisms (Scheffer et al., 2001; Scheffer, 2009; Nyström et al., 2019). Within the context of PB variables, species extinction (i.e. biodiversity loss PB) represents an irreversible process. Resilience has also been suggested as a conceptual framework that could assist in developing paths towards sustainability (Folke et al., 2016). Hence, it can serve as a theoretical and practical foundation for the planetary boundaries framework. An important point to bear in mind, however is that resilience is a property of a system and is neither “good” nor “bad” per se. It can help maintain the current state of a system no matter whether it is deemed desirable or undesirable. The Holocene epoch has allowed development of agriculture, permanent settlements, and the emergence of complex human societies, so maintaining Holocene-like conditions can be deemed desirable, and safeguarding of resilience that support these conditions of critical importance for humanity (Steffen et al., 2018).

Bearing in mind the resilience concept and its importance we aim in this work to specify, in the context of a thermodynamical model of the ES, what are the physical properties that manifest themselves collectively as resilience features of the ES. Our starting point is a thermodynamical model of the ES from Holocene state conditions to other potentially stable states, which can be regarded as phase transitions and admit a description through the Landau-Ginzburg Theory (LGT) (Bertolami and Francisco, 2018, 2019; Barbosa et al., 2020). The LGT is a theoretical framework used in physics to describe phase transitions, such as when a material changes from a solid to a liquid state or a magnetic material loses its magnetism. Here we use the LGT to describe the transitions the ES has gone throughout the history of Earth.

As we shall review in the next section, this framework allows for determining the equilibrium states of the ES in terms of the planet's biophysical subsystems or processes that are, due to the impact of the human activities, the driving forces that dominate its evolution. In the Anthropocene, human activities are here collectively denoted by H. In the phase-transition model discussed in Bertolami and Francisco (2018, 2019) and Barbosa et al. (2020), H was considered an external field, however, in the present work, we admit that through policies and actions, the dynamic features of the ES can be altered so to modify the topographic landscape of possible Anthropocene trajectories. Way to do so include, mitigation strategies, such as halting deforestation and changing agricultural practices that contribute to CO₂ emission; transformation strategies, such as shifting from fossil fuel-based economies to ones based or renewable energy, and; restoration strategies, such as restoration of degraded ecosystems and CO₂ capture technologies.

As previously discussed, the proposed Landau-Ginzburg model allows for getting the evolution equation of the ES, the so-called Anthropocene equation, and to associate the sharp rise of the physical parameters that characterise the ES to the great acceleration of the human activities (Bertolami and Francisco, 2018), which became conspicuous from the second half of the 20th century and onwards (Steffen et al., 2015a).

However, as will be seen below, the original model did not exhibit explicit features that resemble resilience. This is the main purpose of the present work. As the model is based on thermodynamical arguments, one must seek for physical properties that would lead to a more resilient behaviour of the ES. In the context of the model, resilience is regarded as the resistance the ES shows in changing from one equilibrium state to another. At the present transient period, the Anthropocene, it has been hypothesised that the ES is moving away from the Holocene equilibrium state to a new state, potentially a Hothouse Earth state (Steffen et al., 2018) (Fig. 1). As we shall see, our results show that resilience is associated to the existence of metastable states and explicit dissipation of energy that prevent the ES to runway towards the Hothouse Earth state.

Figure 1 A schematic illustration of the evolution of the Earth System with a start from the Neolithic revolution (∼12 000 years ago). Leading up to its current state (i.e. “warm Holocene Earth state”) 7 of 9 planetary boundaries have been transgressed. A continuation on this pathway suggests that the Earth system may end up in a Hothouse Earth state (Steffen et al., 2018) (left pathway). However, explicit dissipation of energy, policies and actions geared at building resilience of a metastable “Holocene-like Earth state” (see also Fig. 2) could provide an opportunity to build a trajectory toward a future “cooling Earth state” (right pathway).

A pleasing feature of the proposed description is that it allows for drawing trajectories of the ES in the phase space of model's variables. By considering that the PBs and the ensued temperature display dynamics that are affected by PBs self-interactions which are shown to be different from zero (Barbosa et al., 2020), two well defined and distinct sets of trajectories were identified upon assumptions about the evolution of the PB: a linear growth of the human activities, H(T)=H₀t, where H₀ is an arbitrary constant, from which follows that all ES trajectories starting at the Holocene are led to Hothouse Earth state (Steffen et al., 2018) (Fig. 1) with a necessarily higher temperature than the Holocene average temperature (Bertolami and Francisco, 2019); if instead, the increase of the human activities impact on the ES obey a discrete logistic map (May, 1976; Jakobson, 1981; Kingsland, 1995), trajectories can display bifurcations or chaotic behaviour (Bernardini et al., 2025). Of course, as human activities are bounded by the finiteness of resources, the logistic map might be a more accurate description of its behaviour, although it is not quite clear what is the time span elapsed between successive steps of the logistic map. In any case, it is relevant to keep in mind that a too fast increase might give origin to trajectory bifurcations or even chaotic behaviour, which, of course, precludes predictions and control measures on the evolution of the ES.

In this work we extend the previous studies of the ES model carried out in Bertolami and Francisco (2018, 2019), Barbosa et al. (2020), and Bernardini et al. (2025) on various aspects. Previously, we aimed to show the inevitability of the Hothouse Earth state given the disestablishing nature of the human activities and the interplay among the PBs. Here, we consider the dynamic features arising from the self-interactions of the 9 identified PBs, here generically denoted as h_i, $i=\mathrm{1},\mathrm{\dots },\mathrm{9}$, and show the specific conditions to implement resilience in the eleven dimensional space $(\mathit{\psi },h_i,F(\mathit{\psi },h_i\left)\right)$, where (ψ) denotes the temperature and F(ψ,h_i) the Helmoholtz free energy (cf. below). Resilience can be regarded as a set of measures that prevent or delay the evolution of the ES towards a Hothouse Earth state and ensuring that this state is as close as possible to the Holocene state². This can be implemented by creating metastable states to avoid a runaway situation due to a barrier that arises as higher-order terms into the Helmholtz free energy are introduced (cf. discussion below). A further requirement is dynamic friction, that is friction introduced via a kinetic energy-type term, to restrict the change of state in the phase space. This is a fairly natural condition as any realistic system dissipates energy. The specific conditions for the ES to acquire effective resilience features will be discussed below. Trajectories of the ES without and with resilience are depicted in Figs. 1 and 2 respectively (cf. a detailed discussion below).

Figure 2 Free energy in function of the temperature, planetary boundaries (H) and resilience features (matastable state).

This paper is organised as follows: in Sect. 2 we review the cardinal aspects of the LGT of the ES and discuss the most relevant features of the dynamical system emerging from the model; in Sect. 3, we discuss the implementation of the resilience features in the model and connect them to properties that any model of the ES should have. Finally, in Sect. 4 we present our conclusions and discuss how our work can be extended to address several issues concerning features and transformation of the global social-ecological system.

2 A Thermodynamical Model for the Earth System

We first review the main features of the proposed model for the ES (Bertolami and Francisco, 2018) and discuss in the next section the conditions to extend it in order to explicitly exhibit resilient properties.

The proposal of Bertolami and Francisco (2018) is to regard transitions of the ES as phase transitions which can be described by the LGT through an order parameter, ψ, and natural parameters (astronomical, geophysical, internal). In the Anthropocene, the natural forces average out to zero and the system is driven by the strength of the human activities, collectively denoted by H. In this approach, the thermodynamic description of the system is obtained through the Helmholtz free energy, F, which can be written as an analytic function of an order parameter, ψ, which is chosen to be the reduced temperature relative to Holocene average temperature, 〈T_H〉, $\mathit{\psi }:=(T-\langle T_{\mathrm{H}}\rangle )/\langle T_{\mathrm{H}}\rangle$. Thus, in the Anthropocene, disregarding the spatial variation of ψ, one can write (Bertolami and Francisco, 2018, 2019):

$$\begin{array}{}\text{(1)}& F(\mathit{\psi },H)=F_{\mathrm{0}}+a{\mathit{\psi }}^{\mathrm{2}}+b{\mathit{\psi }}^{\mathrm{4}}-\mathit{\gamma }H\mathit{\psi },\end{array}$$

where F₀, a, b and γ are constants. The linear term in ψ corresponds to the human activities, which at the Anthropocene can match the quadratic and quartic contributions due to natural causes (astronomic, geological internal).

The strength of the human activities are probed by their impact via the PBs (Steffen et al., 2015b; Rockström et al., 2009), h_i, $i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{9}$ with respect to their Holocene values. Given that the PB can interact among themselves, the most general expression for H is given by (Bertolami and Francisco, 2019):

$$\begin{array}{}\text{(2)}& H=\sum _{i=\mathrm{1}}^{\mathrm{9}}{h}_i+\sum _{i,j=\mathrm{1}}^{\mathrm{9}}{g}_{ij}{h}_i{h}_j+\sum _{i,j,k=\mathrm{1}}^{\mathrm{9}}{\mathit{\alpha }}_{ijk}{h}_i{h}_j{h}_k+\mathrm{\dots },\end{array}$$

where [g_ij] is a non-degenerate, det[g_ij]≠0 9×9 matrix. Similar conditions should be imposed on the coefficients α_ijk and β_ijkl of the higher-order interaction terms. In principle, these interactions terms are sub-dominating, however, their importance has to be established empirically. As pointed out in Bertolami and Francisco (2019), the interaction terms may lead to new equilibrium states and suggest some mitigation strategies depending on their sign and strength in the matrix entries (Bertolami and Francisco, 2019). This will be explicitly discussed in the next section. In Barbosa et al. (2020), it was shown that the interaction term between the climate change variable, CO₂ concentration, say, h₁, and the oceans acidity, say, h₂, was non-vanishing and contributed to about 10 % of the value of the individual contributions themselves.

In order to introduce resilience features into the model, that is, resistance to change from one equilibrium state into another, we have to consider, contrary to previous works (Bertolami and Francisco, 2018, 2019; Barbosa et al., 2020; Bernardini et al., 2025), that the PBs are dynamical variables that are not only passively changed due to human activities, but that can be actively altered so to boost the resilience features of the ES. This allows us to project how the ES would behave depending on its initial state and subsequent trajectory in the phase space of the model, specified through the variables $(\mathit{\psi },\dot{\mathit{\psi }},h_i,\dot{h_i})$. Thus, for a given set of initial conditions, corresponding to a state $\left(\mathit{\psi }\right(\mathrm{0}),\dot{\mathit{\psi }}(\mathrm{0}),h_i(\mathrm{0}),\dot{h_i}(\mathrm{0}\left)\right)$ in the phase space, one can, in principle, obtain the trajectories, orbits, in the phase space after solving the initial value problem through the evolution equations of the system. The equations of motion are obtained through the Lagrangian or equivalently through the Hamiltonian formalism. The latter, yielding to first order differential equations, is more suitable to establish a dynamical system in its canonical form.

The Lagrangian function must include, besides the potential, which is given by the free energy, a set of kinetic energy terms for the canonical coordinates. The simplest possible kinetic term is a quadratic term proportional to the squared first derivative of each coordinate. Thus, we can write the following Lagrangian:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathcal{L}(q,\dot{q})=& {\displaystyle \frac{\mathit{\mu }}{\mathrm{2}}}{\dot{\mathit{\psi }}}^{\mathrm{2}}+{\displaystyle \frac{\mathit{\nu }}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{\mathrm{9}}{\dot{h_i}}^{\mathrm{2}}\\ & -F_{\mathrm{0}}-a{\mathit{\psi }}^{\mathrm{2}}-b{\mathit{\psi }}^{\mathrm{4}}+\mathit{\gamma }H\mathit{\psi },\end{array}\end{array}$$

where μ and ν are arbitrary constants and the dots stand for time derivatives. The constant ν is assumed to be the same for all PB variables.

Aiming to get the Hamiltonian function, we evince the relevant canonical conjugate momenta associated to ψ and to a generic PB variable, h_i:

$$\begin{array}{}\text{(4)}& {\displaystyle }{p}_{\mathit{\psi }}={\displaystyle \frac{\partial \mathcal{L}}{\partial \dot{\mathit{\psi }}}}=\mathit{\mu }\dot{\mathit{\psi }},\text{(5)}& {\displaystyle }{p}_{h_i}={\displaystyle \frac{\partial \mathcal{L}}{\partial \dot{h_i}}}=\mathit{\nu }\dot{h_i},\end{array}$$

from which follows the Hamiltonian function

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathcal{H}(\mathit{\psi },p)=& {\displaystyle \frac{p_{\mathit{\psi }}^{\mathrm{2}}}{\mathrm{2}\mathit{\mu }}}+\sum _{i=\mathrm{1}}^{\mathrm{9}}{\displaystyle \frac{p_{h_i}^{\mathrm{2}}}{\mathrm{2}\mathit{\nu }}}\\ & +F_{\mathrm{0}}+a{\mathit{\psi }}^{\mathrm{2}}+b{\mathit{\psi }}^{\mathrm{4}}-\mathit{\gamma }H\mathit{\psi },\end{array}\end{array}$$

and Hamilton's equations,

$$\begin{array}{}\text{(7)}& {\displaystyle }\dot{\mathit{\psi }}={\displaystyle \frac{\partial \mathcal{H}}{\partial p_{\mathit{\psi }}}},\dot{p_{\mathit{\psi }}}=-{\displaystyle \frac{\partial \mathcal{H}}{\partial \mathit{\psi }}},\text{(8)}& {\displaystyle }\dot{h_{\mathrm{1}}}={\displaystyle \frac{\partial \mathcal{H}}{\partial p_{h_i}}},\dot{p_{h_i}}=-{\displaystyle \frac{\partial \mathcal{H}}{\partial h_i}}.\end{array}$$

The equations of motion read, considering for while just the contribution from the lowest order terms in Eq. (2):

$$\begin{array}{}\text{(9)}& \mathit{\mu }\ddot{\mathit{\psi }}=-\mathrm{2}a\mathit{\psi }-\mathrm{4}b{\mathit{\psi }}^{\mathrm{3}}+\mathit{\gamma }H\end{array}$$

and

$$\begin{array}{}\text{(10)}& \mathit{\nu }\ddot{h_i}=\mathit{\gamma }\mathit{\psi }.\end{array}$$

To exemplify the behaviour of variables ψ and h_i, let us obtain the resulting solutions for the simple case considered in Bertolami and Francisco (2019). For b≃0, we can neglect the cubic term in the equation of motion for ψ to get the equation of an harmonic oscillator under the action of an external force, H(t). This yields for the simple case of an initial linear time evolution,

$$\begin{array}{}\text{(11)}& H\left(t\right)=H_{\mathrm{0}}t,\end{array}$$

for an equilibrium initial state, $\dot{\mathit{\psi }}\left(\mathrm{0}\right)=\mathrm{0}$, the analytical solution:

$$\begin{array}{}\text{(12)}& \mathit{\psi }\left(t\right)={\mathit{\psi }}_{\mathrm{0}}\cos \left(\mathit{\omega }t\right)+\mathit{\alpha }t,\end{array}$$

where $\mathit{\omega }=\sqrt{\mathrm{2}a/\mathit{\mu }}$ is an angular frequency, $\mathit{\alpha }=\mathit{\gamma }{H}_{\mathrm{0}}/\mathrm{2}a$ and ψ₀ is an arbitrary constant fixed by the initial conditions.

The solution for the impact on the PB, h_i(t), which initially behaves collectively as Eq. (11), that is ${\sum }_{i=\mathrm{1}}^{\mathrm{9}}{h}_i(t\simeq \mathrm{0})=H_{\mathrm{0}}$, quickly evolves to a cubic growth in time:

$$\begin{array}{}\text{(13)}& h_i\left(t\right)=A\cos \left(\mathit{\omega }t\right)+B{t}^{\mathrm{3}}+{\mathit{\alpha }}_{i}t,\end{array}$$

where $A=-\mathit{\gamma }{\mathit{\psi }}_{\mathrm{0}}/\mathit{\nu }{\mathit{\omega }}^{\mathrm{2}}$, $B=\mathit{\alpha }\mathit{\gamma }/\mathrm{6}\mathit{\nu }$, for an arbitrary α_i.

These solutions show that if the temperature ψ grows from an initial linear collective behaviour of the PBs, H=H₀t, then it quickly drives the h_is towards a cubic growth. Clearly, this model shows no resilience features as depicted in Fig. 1, where one clearly sees that from the Holocene, Anthropocene trajectories inevitably evolve towards a Hothouse Earth state.

In what follows we shall consider the introduction into the free energy function of a cubic term for ψ and higher than linear order terms for the PBs as these will allow for metastable states to arise, thus leading to bounded solutions for ψ and the PBs. Metastable states correspond to potential intermediate energy states between the Holocene state and the Hothouse Earth least energy state. In the LGT, metastable states can be considered and studied through cubic terms in the Helmholtz free energy. The conditions for the appearance of metastable states were already discussed in a completely different context, namely in a proposal to classify rocky planets (Bertolami and Francisco, 2022), using the ideias developed in Bertolami and Francisco (2018, 2019) and Barbosa et al. (2020) to describe the ES. In concrete terms, cubic terms might arise from PB interactions that have a strong dependence on the temperature.

Before concluding this discussion it is worth stressing once again that the behaviour of the ES depends crucially on the assumptions about the evolution of the PB. Indeed, as pointed out in the introduction, the supposition that human activities grow linearly as in Eq. (11) implies, as exemplified above, that ES trajectories lead to a potential “Hothouse Earth” state (Bertolami and Francisco, 2019) as discussed by Steffen et al. (2018). However, if the human activities impact on the ES behaves as a discrete logistic map³, as suggested in Bernardini et al. (2025), then evolution will depend the rate of growth of human activities as solutions admit regular trajectories as well as trajectories that present bifurcations and even chaotic behaviour. In the next section we shall consider the features that must be introduced in the Helmholtz free energy and the conditions they must satisfy in order to avoid the ES evolves towards the Hothouse Earth state.

3 Setting up the physical principles of resilience

As mentioned above, in this model resilience features are associated to bounded trajectories in the Anthropocene and these ask for the existence of metastable states. In the LGT the metastable states arise by intruding cubic terms on the free energy. As pointed out in Bertolami and Francisco (2022), the introduction of a cubic term allows for a richer variety of equilibrium states. Indeed, consider the free energy:

$$\begin{array}{}\text{(14)}& F(\mathit{\psi },H)=F_{\mathrm{0}}+a{\mathit{\psi }}^{\mathrm{2}}-c|\mathit{\psi }{|}^{\mathrm{3}}+b{\mathit{\psi }}^{\mathrm{4}}-\mathit{\gamma }H\mathit{\psi },\end{array}$$

where we assume that constants b, c and γ are positive, while constant a can be negative.

The existence of extrema is given by two conditions. The first one reads:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial F(\mathit{\psi },H)}{\partial \mathit{\psi }}}=\mathrm{0}=\mathrm{2}a\mathit{\psi }-\mathrm{3}c{\mathit{\psi }}^{\mathrm{2}}+\mathrm{4}b{\mathit{\psi }}^{\mathrm{3}}-\mathit{\gamma }H.\end{array}$$

The resulting cubic equation admits at least one real solution, say, ψ_M, meaning that there are at least two metastable states, ψ_M and −ψ_M. Clearly, ψ_M≠0 as far as H≠0.

However, the unboundedness of the evolution of the variables (ψ,h_i) is due to the unboundedness of the PBs. Recent assessment of the PBs has shown that 7 out of the 9 PBs have gone beyond their Holocene values where they were at equilibrium, a state usually referred to as Safe Operating Space (SOS).

The motion in the eleven-dimensional configuration space, $(\mathit{\psi },h_i,F(\mathit{\psi },h_i\left)\right)$, is quite complex, so in order to simplify the analysis we consider one single generic PB, h_i, and assume that the remaining ones are unchanged⁴. The free energy can be written explicitly in terms of the high order contributions in H depicted in Eq. (2). We consider the essential set of terms in order to carry out the minimisation procedure, that is:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}F(\mathit{\psi },H)=& \widehat{F_{\mathrm{0}}}+a{\mathit{\psi }}^{\mathrm{2}}-c|\mathit{\psi }{|}^{\mathrm{3}}+b{\mathit{\psi }}^{\mathrm{4}}\\ & -\mathit{\gamma }(h_i+g_i{h}_i^{\mathrm{2}}+b_i{h}_i^{\mathrm{3}})\mathit{\psi },\end{array}\end{array}$$

where we have aggregated all contributions to the quadratic and cubic terms in h_i, a generic PB, within the constants g_i and b_i. To ensure boundedness it is necessary that g_i is negative and that b_i is positive.

Thus, from Eq. (16), one gets the condition:

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\partial F(\mathit{\psi },h_i)}{\partial h_i}}=\mathrm{1}+\mathrm{2}{g}_i{h}_i+\mathrm{3}{b}_i{h}_i^{\mathrm{2}}=\mathrm{0},\end{array}$$

which admits real non-vanishing solutions, h_iM. as far as $g_i^{\mathrm{2}}>\mathrm{3}{b}_{\mathrm{1}}$ for b_i≠0 or $h_{iM}=-\frac{\mathrm{1}}{\mathrm{2}{g}_i}$ if b₁=0.

The general conditions to ensure that the extremum (ψ_M,h_iM) corresponds to a minimum and hence to a metastable state are given by:

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}F({\mathit{\psi }}_M,h_{iM})}{\partial {\mathit{\psi }}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}F({\mathit{\psi }}_M,h_{iM})}{\partial h_i^{\mathrm{2}}}}-{\left({\displaystyle \frac{{\partial }^{\mathrm{2}}F({\mathit{\psi }}_M,h_{iM})}{\partial \mathit{\psi }\partial h_i}}\right)}^{\mathrm{2}}>\mathrm{0}.\end{array}$$

and

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}F({\mathit{\psi }}_M,h_{iM})}{\partial {\mathit{\psi }}^{\mathrm{2}}}}>\mathrm{0},\end{array}$$

which yield the following relationships:

$$\begin{array}{}\text{(20)}& g_i<-\mathrm{3}{b}_i{h}_{iM}\end{array}$$

and

$$\begin{array}{}\text{(21)}& \mathrm{2}a-\mathrm{6}c\left|{\mathit{\psi }}_M\right|+\mathrm{12}b{\mathit{\psi }}_M^{\mathrm{2}}>\mathrm{0}.\end{array}$$

Satisfying these conditions imply the ES can settle in a the metastable state, (ψ_M,h_iM), that is, the system shows resilience and does not runaway towards the “Hothouse Earth” state as depicted in Fig. 2 as far as $\mathrm{3}{b}_i<g_i^{\mathrm{2}}<\mathrm{9}{b}_i^{\mathrm{2}}{h}_{iM}$.

Notice that the conditions for the existence of a metastable state can be met if g_i<0 even if coefficients b_i vanish. This is quite welcome as these coefficients are associated to higher-order interaction terms, which from phenomenological considerations, are presumably small. On the other hand, a non-vanishing and negative contribution from the quadratic term $h_i^{\mathrm{2}}$ is absolutely necessary. Actually, the concrete case studied in Barbosa et al. (2020) shows that this is possible. Furthermore, condition Eq. (21) can be satisfied if a<0.

Another feature associated to resilience is the “inertia” that the ES shows in changing from a given state to another. This feature can be identified with the ubiquitous dissipation of energy present in any physical system. Most often dynamical dissipation processes can be described through velocity-proportional frictional forces which imply that just part of the free energy of a system is turned into kinetic energy, that is, motion of the system. In the Lagrangian/Hamiltonian formalism for a particle, the effect of these forces can be accounted through the Rayleigh dissipation function, $R=-\mathit{\kappa }{p}^{\mathrm{2}}/m$, where κ is a constant, p is the canonical conjugate momentum and m the mass of the particle.

For the ES, introducing dissipation through the Rayleigh function implies that the left hand side of the equations of motion (9) and (10) acquire extra terms $-{\mathit{\kappa }}_{\mathit{\psi }}\dot{\mathit{\psi }}$ and $-{\mathit{\kappa }}_{h_i}\dot{h_i}$, respectively. The effect of these terms is to reduce the amplitude of the motion of the ES once it goes from one state to another, thus acting as a resistance of the system to the change of its state. This can be clearly associated to resilience.

These considerations are sufficient for setting the physical conditions for the resilience of the ES. As we have seen, a metastable state corresponding to the solution (ψ_M,h_iM) of Eqs. (15), (18), and (19), whose free energy (Eq. 16) coefficients satisfy the conditions Eq. (20) and Eq. (21) together with the unavoidable dynamic friction/energy dissipation that exists in any system are the physical properties that endow the ES for having a resilient behaviour. For sure, further research is needed in order to establish which PBs are more suitable for setting up the conditions obtained above. This means that the PB properties concerning their dependence on the temperature and strength of their self-interaction and with other PBs must be further studied.

4 Conclusions

In this work we have considered the physical principles to ascertain the conditions of resilience in a LGT model of the ES. In order to implement resilience features we have endowed and considered modifications of the free energy so to ensure the existence of metastable states. Furthermore, we have modelled the ES capability to remain in an equilibrium state by arguing that it can be suitably prevented to runway towards a potential Hothouse Earth state by the presence of metastable states whose existence conditions were explicitly shown and the unavoidable dissipation of energy during the evolution of the relevant variables.

Indeed, we have shown that, thanks to the PBs interactions, a metastable state (ψ_M,h_iM) can exist if the conditions, Eqs. (20) and (21), for the coefficients of the free energy, Eq. (16), are satisfied. As pointed out in the above discussion, these conditions can be satisfied even if coefficients b_i vanish as far as g_i<0.

Based on the observational data, it is possible to infer that the metastable state found above might correspond either to an actual state that the ES is close to reach or to a state that can be reached by policy and actions (i.e. mitigation, transformation and restoration strategies) to drive the ES away from the Anthropocene traps it seems to be currently entangled in (see Søgaard Jørgensen et al., 2023, for a description of the 14 major Anthropocene traps).

A recent assessment has shown that 7 out of the 9 PBs have been crossed (Kitzmann et al., 2025) meaning that the evolution of most of the PBs is uncontrolled. Moreover it is unclear if the ES has already reached a point of no return, but it is evident that urgent measures to reverse the current development are needed. In fact, no single set of measures seems to be sufficient to halt the evolution of the PBs beyond the safe operating state. Two of the PBs that deserve particular attention are climate change and biosphere integrity. Both are deemed “core” because their essential role in the ES. The climate system reflects the distribution and balance of energy at the Earth's surface, while the controls material and energy flows, helping to strengthen the systems's resilience against both rapid and long-term changes. This calls for a concerted action involving stewardship measures (Steffen et al., 2015a, 2011; Bertolami, 2022), bringing into the economy (internalising) the workings of the ES (see e.g. Bertolami, 2024) and making them become part of revised economic paradigms (Sureth et al., 2023; Bertolami and Gonçalves, 2024, 2025), mitigation strategies that may include technological carbon sequestration (see e.g. Bertolami, 2025; Bertolami and de Matos, 2024, and refs. therein), and storage as means to curb climate overshoot, to avoid irreversible changes to the ES that will compromise the navigation space for the future generations. Given that tipping in some major ecosystems constituting the Earth System is already being observed, such as in Pacific coral reef systems, and that other tipping elements including the Amazon rainforest and the polar ice sheets may be approaching the threshold of self-reinforcing tipping dynamics (Lenton et al., 2025), a critical question arises as to whether irreversible damage to the Earth System has already occurred, or is imminent. The answer comes only through the understanding of the mechanisms of resilience and how their boosting, through the PB interactions, can be effective. We hope that our work can provide a modest help in this respect.