Lévy noise versus Gaussian-noise-induced transitions in the Ghil–Sellers energy balance model

We study the impact of applying stochastic forcing to the Ghil–Sellers energy balance climate model in the form of a fluctuating solar irradiance. Through numerical simulations, we explore the noise-induced transitions between the competing warm and snowball climate states. We consider multiplicative stochastic forcing driven by Gaussian and α-stable Lévy – $\mathit{\alpha }\in (\mathrm{0},\mathrm{2})$ – noise laws, examine the statistics of transition times, and estimate the most probable transition paths. While the Gaussian noise case – used here as a reference – has been carefully studied in a plethora of investigations on metastable systems, much less is known about the Lévy case, both in terms of mathematical theory and heuristics, especially in the case of high- and infinite-dimensional systems. In the weak noise limit, the expected residence time in each metastable state scales in a fundamentally different way in the Gaussian vs. Lévy noise case with respect to the intensity of the noise. In the former case, the classical Kramers-like exponential law is recovered. In the latter case, power laws are found, with the exponent equal to −α, in apparent agreement with rigorous results obtained for additive noise in a related – yet different – reaction–diffusion equation and in simpler models. This can be better understood by treating the Lévy noise as a compound Poisson process. The transition paths are studied in a projection of the state space, and remarkable differences are observed between the two different types of noise. The snowball-to-warm and the warm-to-snowball most probable transition paths cross at the single unstable edge state on the basin boundary. In the case of Lévy noise, the most probable transition paths in the two directions are wholly separated, as transitions apparently take place via the closest basin boundary region to the outgoing attractor. This property can be better elucidated by considering singular perturbations to the solar irradiance.

1 Introduction

1.1 Multistability of the Earth's climate

The climate system comprises the following five interacting subdomains: the atmosphere, the hydrosphere (water in liquid form), the upper layer of the lithosphere, the cryosphere (water in solid form), and the biosphere (ecosystems and living organisms). The climate is driven by the inhomogeneous absorption of incoming solar radiation, which sets up nonequilibrium conditions. The system reaches an approximate steady state, where macroscopic fluxes of energy, momentum, and mass are present throughout its domain, and entropy is continuously generated and expelled into the outer space. The climate features variability on a vast range of spatial and temporal scales as a result of the interplay of forcing, dissipation, feedbacks, mixing, transport, chemical reactions, phase changes, and exchange processes between the subdomains (see Peixoto and Oort, 1992, Lucarini et al., 2014a, Ghil, 2015, and Ghil and Lucarini, 2020).

In the late 1960s Budyko (1969) and Sellers (1969) independently proposed that, in the current astronomical and astrophysical configuration, the Earth could support two distinct climates, namely the present-day warm (W) state and a competing one characterised by global glaciation, usually referred to as the snowball (SB) state. Their analysis was performed using one-dimensional energy balance models (EBMs), which, despite their simplicity, were able to capture the essential physical mechanisms in action, i.e. the interplay between two key feedbacks. The Boltzmann feedback is associated with the fact that warmer bodies emit more radiation, and it is a negative, stabilising one. Instead, the instability of the system is due to the presence of the so-called ice–albedo feedback, whereby an increase in the ice-covered fraction of the surface leads to a further temperature reduction for the planet because ice efficiently reflects the incoming solar radiation. These mechanisms are active at all spatial scales, including the planetary one (see Budyko, 1969, and Sellers, 1969). Such pioneering investigations of the multistability of the Earth's climate were later extended by Ghil (1976) (see also the later analysis by Ghil and Childress, 1987), who provided a comprehensive mathematical framework for the problem, based on the study of the bifurcations of the system. The main control parameter defining the stability properties is the solar irradiance S^∗. Below the critical value S_W→SB, only the SB state is permitted, whereas, above the critical value S_SB→W, only the W state is permitted. Such critical values, which determine the region of bistability, are defined by bifurcations that emerge when, roughly speaking, the strength of the positive, destabilising feedbacks becomes as strong as the negative, stabilising feedbacks. Many variants of the models proposed by Budyko and Sellers have been discussed in the literature, all featuring by and large rather similar qualitative and quantitative features (Ghil, 1981; North et al., 1981; North, 1990; North and Stevens, 2006). Furthermore, these models have long been receiving a great deal of attention from the mathematical community regarding the possibility of proving the existence of solutions and evaluating their multiplicity (Hetzer, 1990; Díaz et al., 1997; Kaper and Engler, 2013; Bensid and Díaz, 2019).

Only later were these predictions confirmed by actual data. Indeed, geological and paleomagnetic evidence suggests that, during the Neoproterozoic era, between 630×10⁶ and 715×10⁶ years ago, the Earth went, at least twice, into major long-lasting global glaciations that can be associated with the SB state (see Pierrehumbert et al., 2011, and Hoffman et al., 1998). Multicellular life emerged in our planet shortly after the final deglaciation from the last SB state (Gould, 1989). The robustness and importance of the competition between the Boltzmann feedback and the ice–albedo feedback in defining the global stability properties of the climate has been confirmed by investigations performed using higher complexity models (Lucarini et al., 2010; Pierrehumbert et al., 2011), including fully coupled climate models (Voigt and Marotzke, 2010). While the mechanisms described above are pretty robust, the concentration of greenhouse gases and the boundary conditions defined by the extent and position of the continents have an impact on the values of S_W→SB and S_SB→W, as well as on the properties of the competing states. The presence of multistability has a key importance in terms of determining habitability conditions for Earth-like exoplanets (see Lucarini et al., 2013, and Linsenmeier et al., 2015).

Additionally, several results indicate that the phase space of the climate system might well be more complex than the scenario of bistability described above. Various studies (Lewis et al., 2007; Abbot et al., 2011; Lucarini and Bódai, 2017; Margazoglou et al., 2021) performed with highly nontrivial climate models report the possible existence of additional competing states, up to a total of five (Brunetti et al., 2019; Ragon et al., 2022). In Margazoglou et al. (2021), it is argued that, in fact, one can see the climate as a multistable system where multistability is realised at different hierarchical levels. As an example, the tipping points (Lenton et al., 2008; Steffen et al., 2018) that characterise the current (W) climate state can be seen as a manifestation of a hierarchically lower multistability with respect to the one defining the dichotomy between the W and SB states.

1.2 Transitions between competing metastable states: Gaussian vs. Lévy noise

Clearly, in the case of autonomous systems where the phase space is partitioned in more than one basin of attraction of the corresponding attractors and the basin boundaries, the asymptotic state of the system is determined by its initial conditions. Things change dramatically when one includes time-dependent forcing which allows for transitions between competing metastable states (Ashwin et al., 2012). In particular, following the viewpoint originally proposed by Hasselmann (1976), whereby the fast variables of the climate system act as stochastic forcings for the slow ones (Imkeller and von Storch, 2001), the relevance of studying noise-induced transitions between competing states has become apparent (Hänggi, 1986; Freidlin and Wentzell, 1984). This viewpoint, where the noise is usually assumed to be Gaussian distributed, has provided very valuable insight on the multiscale nature of climatic time series (Saltzman, 2001) and is related to the discovery of phenomena like stochastic resonance (Benzi et al., 1981; Nicolis, 1982).

Metastability is ubiquitous in nature, and advancing its understanding is a key challenge in complex system science at large (Feudel et al., 2018). In general, the transitions between competing metastable states in stochastically perturbed multistable systems take place, in the weak noise limit, through special regions of the basin boundaries, which are named edge states. The edge states are saddles, and the trajectories initialised in the basin boundaries are attracted to them, but there is an extra direction of instability, so that a small perturbation sends an orbit towards one of the competing metastable states with a probability of one (Grebogi et al., 1983; Ott, 2002; Kraut and Feudel, 2002; Skufca et al., 2006; Vollmer et al., 2009). In the case the edge state supports chaotic dynamics, we refer to it as the melancholia (M) state (Lucarini and Bódai, 2017). In previous papers, we have shown that it is possible to construct M states in high-dimensional climate models (Lucarini and Bódai, 2017) and to prove that the nonequilibrium quasi-potential formalism introduced by Graham (1987) and Graham et al. (1991) provides a powerful framework for explaining the population of each metastable state and the statistics of the noise-induced transitions. In the weak noise limit, edge states act as gateways for noise-induced transitions between the metastable states (Lucarini and Bódai, 2019, 2020; Margazoglou et al., 2021); see also a recent study on a nontrivial metastable prey–predator model (Garain and Sarathi Mandal, 2022). The local minima and the saddles of the quasi-potential Φ, which generalises the classical energy landscape for non-gradient systems, correspond to competing metastable states and to edge states, respectively. In our investigation, the climate system is forced by adding a random – Gaussian-distributed – component to the solar irradiance, which impacts, in the form of multiplicative noise, only a small subset of the degrees of freedom of the system. We remark that such a choice of the stochastic forcing does not fully reflect the physical realism, as the variability of the solar irradiance has a more complex behaviour (Solanki et al., 2013). Instead, noise acts as a tool for exploring the global stability properties of the system, and injecting noise as fluctuation of the solar irradiance has the merit of impacting the Lorenz energy cycle, thus effecting all degrees of freedom of the system (Lucarini and Bódai, 2020). See also the recent detailed mathematical analysis of the stochastically perturbed one-dimensional EBMs presented in Díaz and Díaz (2021).

A major limitation of this mathematical framework is the need to rigidly consider Gaussian noise laws, even if considerable freedom is left as to the choice of the spatial correlation properties of the noise. It seems natural to attempt a generalisation by considering the whole class of α-stable Lévy noise laws. Lévy processes (Applebaum, 2009; Duan, 2015), described in detail in Appendix A, are fundamentally characterised by the stability parameter $\mathit{\alpha }\in (\mathrm{0},\mathrm{2}]$, where the α=2 case corresponds to the Gaussian case (which is, indeed, a special Lévy process). In what follows, when we discuss Lévy noise laws, we refer to $\mathit{\alpha }\in (\mathrm{0},\mathrm{2})$.

Note that α-stable Lévy processes have played an important role in geophysics as they have provided the starting point for defining the multiplicative cascades also referred to as universal multifractal. This framework has been proposed as way to analyse and simulate, at climate scales, the ubiquitous intermittency and heavy-tailed statistics of clouds (Schertzer and Lovejoy, 1988), rain reflectivity (Tessier et al., 1993; Schertzer and Lovejoy, 1997), atmospheric turbulence (Schmitt et al., 1996), and soil moisture (Millàn et al., 2016). On longer timescales, multiplicative cascades have been used to interpret temperature records in the summit ice core Schmitt et al. (1995, see Lovejoy and Schertzer, 2013, for a summary of this viewpoint). Mathematicians, on the other hand, have defined a Lévy multiplicative chaos (Fan, 1997; Rhodes et al., 2014) as a more mathematically tractable alternative to the universal multifractal. Finally, we remark that fractional Fokker–Planck equations have been proposed by Schertzer et al. (2001) to investigate the properties of nonlinear Langevin-type equations forced by an α-stable Lévy noise with the goal of analysing and simulating anomalous diffusion.

Following Ditlevsen (1999), it has become apparent that more general classes of α-stable Lévy noise laws might be useful for modelling noise-induced transitions in the climate system like Dansgaard–Oeschger events, which are sequences of periods of abrupt warming followed by slower cooling that occurred during the last glacial period (Barker et al., 2011). The viewpoint by Ditlevsen (1999) was particularly effective in stimulating mathematical investigations into noise-induced escapes from attractors, where, as stochastic forcing, one chooses a Lévy, rather than Gaussian, noise (Imkeller and Pavlyukevich, 2006a, b; Chechkin et al., 2007; Debussche et al., 2013). Such analyses have clarified that a fundamental dichotomy exists with the classical Freidlin and Wentzell scenario mentioned above, even if phenomena like stochastic resonance can also be recovered in this case (Dybiec and Gudowska-Nowak, 2009; Kuhwald and Pavlyukevich, 2016). Whereas, in the Gaussian case, transitions between competing attractors occur as a result of the very unlikely combination of many steps all going in the right direction, in the Lévy case, transitions result from individual, very large and very rare jumps. Recently, Duan and collaborators have made fundamental progress in achieving a variational formulation of the Lévy noise-perturbed dynamical systems (Hu and Duan, 2020) and in developing corresponding methods for data assimilation (Gao et al., 2016) and data analysis (Lu and Duan, 2020). In terms of applications, Lévy noise is becoming a more and more a popular concept and tool for studying and interpreting complex systems (Grigoriu and Samorodnitsky, 2003; Penland and Sardeshmukh, 2012; Zheng et al., 2016; Wu et al., 2017; Serdukova et al., 2017; Cai et al., 2017; Singla and Parthasarathy, 2020; Gottwald, 2021).

The contribution by Gottwald (2021) is especially worth recapitulating because of its methodological clarity. There, the idea is, following Ditlevsen (1999), to provide a conceptual deterministic climate model able to generate a Lévy-noise-like signal to describe, at least qualitatively, abrupt climate changes similar to Dansgaard–Oeschger events. A key building block is the idea proposed in Thompson et al. (2017) that a Lévy noise can be produced by integrating the so-called correlated additive and multiplicative (CAM) noise processes, which are defined by starting from standard Gaussian processes. The other key ingredient is to consider the atmosphere as the fast component in the multiscale model and deduce, using homogenisation theory (Pavliotis and Stuart, 2008; Gottwald and Melbourne, 2013), that its influence on the slower climate components can be closely represented as a Gaussian forcing. Finally, the temperature signal is cast as the integral of a CAM process.

We remark that Gaussian and Lévy noise can be associated with stochastic forcings of a fundamentally different nature. One might think of Gaussian noise as being associated to the impact of very rapid unresolved scales of motion on the resolved ones Pavliotis and Stuart (2008). Instead, one might interpret α-stable Lévy noise as describing, succinctly, the impact of what in the insurance sector are called acts of God (e.g. an asteroid hitting the Earth, a massive volcanic eruption, or the sudden collapse of the West Antarctic ice sheet).

1.3 Outline of the paper and main results

We consider here the Ghil–Sellers Earth's EBM (Ghil, 1976), a diffusive, one-dimensional energy balance system, governed by a nonlinear reaction–diffusion parabolic partial differential equation. We stochastically perturb the system by adding random fluctuations to the solar irradiance; therefore, the noise is introduced in a multiplicative form. We study the transitions between the two competing metastable climate states and carry out a comparison of the effect of considering Lévy vs. Gaussian noise laws of weak intensity ε.

The main challenges of the problem are (a) the fact that we are considering dynamical processes occurring in infinite dimensions (Doering, 1987; Duan and Wang, 2014; Alharbi, 2021) and (b) the consideration of multiplicative Lévy noise laws (Peszat and Zabczyk, 2007; Debussche et al., 2013). We characterise noise-induced transitions between the competing climate basins and quantify the effect of noise parameters on them by estimating the statistics of escape times and empirically constructing mean transition pathways called instantons.

The results obtained confirm that, in the weak noise limit ε→0, the mean residence time in each metastable state driven by Gaussian vs. Lévy noise has a fundamentally different dependence on ε. Indeed, as expected, in the Gaussian case, the residence time grows exponentially with ε⁻², which is, thus, in basic agreement with the well-known Kramers (1940) law and the previous studies performed on climate models (Lucarini and Bódai, 2019, 2020). Instead, in the case of α-stable noise laws, the residence time increases with ε^−α. We perform simulations for $\mathit{\alpha }=\mathit{\{}\mathrm{0.5},\mathrm{1.0},\mathrm{1.5}\mathit{\}}.$ The obtained scaling can be explained by effectively treating the Lévy noise as a compound Poisson process, and this is in agreement with what is found for low-dimensional dynamics (Imkeller and Pavlyukevich, 2006a, b) and for the infinite dimensional stochastic Chafee–Infante reaction–diffusion equation (Debussche et al., 2013) in the case of additive noise. This might indicate that such scaling laws are more general than what has been so far assumed.

Furthermore, we find clear confirmation that, in the case of Gaussian noise in the weak noise limit, the escape from either attractor's basin takes place through the edge state. Indeed, the most probable paths for both thawing and freezing processes meet at the edge state and have distinct instantonic and relaxation sections. In turn, for Lévy noise in the weak noise limit, the escapes from a given basin of attraction occur through the boundary region closest to the outgoing attractor. Hence, the paths are very different from the Gaussian case (especially so for the freezing transition) and, somewhat surprisingly, are identical regardless of the value of α considered. These properties can be better understood by studying the impact of including singular perturbations to the value of the solar irradiance.

The rest of the paper is organised as follows. In Sect. 2, we present the Ghil–Sellers EBM and summarise its most important dynamical aspects and the steady-state solutions and their stability. The stochastic partial differential equation obtained by randomly perturbing the solar irradiance in the EBM is given in Sect. 3, where we also clarify the mathematical meaning of the solution of the stochastic partial differential equation. Section 3 also introduces the mean residence time and most probable transition path between the competing climate states. The numerical methods are also briefly presented. In Sect. 4, we discuss our main results. In Sect. 5, we present our conclusions and perspectives for future investigations. Finally, Appendix A presents a succinct description of α-stable Lévy processes, Appendix B sketches the derivation of the scaling laws for mean residence times presented in Debussche et al. (2013), Appendix C explores the behaviour and dynamics of singular Lévy perturbations of different duration, and Appendix D presents a tabular summary of the statistics of the problem.

2 The Ghil–Sellers energy balance climate model

The Ghil–Sellers EBM (Ghil, 1976) is described by a one-dimensional nonlinear, parabolic, reaction–diffusion partial differential equation (PDE) involving the surface temperature T field and the transformed space variable $x=\mathrm{2}\mathit{\varphi }/\mathit{\pi }\in [-\mathrm{1},\mathrm{1}]$, where $\mathit{\varphi }\in [-\mathit{\pi }/\mathrm{2},\mathit{\pi }/\mathrm{2}]$ is the latitude. The model describes the processes of energy input, output, and diffusion across the domain and can be written as follows:

$$\begin{array}{}\text{(1)}& C\left(x\right)T_t=D_I(x,T,T_x,T_{xx})+D_{\text{II}}(x,T)-D_{\text{III}}\left(T\right),\end{array}$$

where C(x) is the effective heat capacity, and $T=T(x,t)$ has boundary and initial conditions, as follows:

$$\begin{array}{}\text{(2)}& T_x(-\mathrm{1},t)=T_x(\mathrm{1},t)=\mathrm{0},T(x,\mathrm{0})=T_{\mathrm{0}}\left(x\right).\end{array}$$

The equation does not depend explicitly on the time t. The subscripts _t and _x refer to partial differentiation. The first term – D_I – on the right-hand side of Eq. (1) can be written as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& D_I(x,T,T_x,T_{xx})\\ & ={\displaystyle \frac{\mathrm{4}}{{\mathit{\pi }}^{\mathrm{2}}\cos (\mathit{\pi }x/\mathrm{2})}}\left[\cos \right(\mathit{\pi }x/\mathrm{2}\left)K\right(x,T)T_x{]}_x,\end{array}\end{array}$$

and describes the convergence of meridional heat transport performed by the geophysical fluids. The function K(x,T) is a combined diffusion coefficient, expressed as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }K(x,T)=k_{\mathrm{1}}\left(x\right)+k_{\mathrm{2}}\left(x\right)g\left(T\right),\text{with}\text{(5)}& {\displaystyle }g\left(T\right)={\displaystyle \frac{c_{\mathrm{4}}}{T^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{c_{\mathrm{5}}}{T}}\right).\end{array}$$

The empirical functions k₁(x) and k₂(x) are eddy diffusivities for sensible and latent heat, respectively, and g(T) is associated with the Clausius–Clapeyron relation, which describes the relationship between temperature and saturation water vapour content of the atmosphere.

The second term – D_II – on the right-hand side of Eq. (1) describes the energy input associated with the absorption of incoming solar radiation and can be written as follows:

$$\begin{array}{}\text{(6)}& D_{\text{II}}(x,T)=\mathit{\mu }\mathcal{Q}\left(x\right)\left[\mathrm{1}-{\mathit{\alpha }}_{\mathrm{a}}(x,T)\right],\end{array}$$

where 𝒬(x) is the incoming solar radiation, and α_a(x,T) is the surface reflectivity (albedo), which is expressed as follows:

$$\begin{array}{}\text{(7)}& {\mathit{\alpha }}_{\mathrm{a}}(x,T)={\left\{b\left(x\right)-c_{\mathrm{1}}\left(T_m+min\left[T-c_{\mathrm{2}}z\left(x\right)-T_m,\mathrm{0}\right]\right)\right\}}_{\mathrm{c}},\end{array}$$

where the subscript {⋅}_c denotes a cutoff for a generic quantity h, defined as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{h}_{\mathrm{c}}=\left\{\begin{array}{ll}{h}_{min}& h\le h_{min},\\ h& h_{min}<h<h_{max},\\ h_{max}& h_{max}\le h.\end{array}\right.\end{array}$$

The term c₂z(x) in Eq. (7) indicates the difference between the sea level and surface level temperatures, and b(x) is a temperature independent empirical function of the albedo. The parameterisation given in Eqs. (7)–(8) encodes the positive ice–albedo feedback. The relative intensity of the solar radiation in the model can be controlled by the parameter μ.

The last term – D_III – on the right-hand side of Eq. (1) describes the energy loss to space by outgoing thermal planetary radiation and is responsible for the negative Boltzmann feedback. It can be written as follows:

$$\begin{array}{}\text{(9)}& D_{\text{III}}\left(T\right)=\mathit{\sigma }{T}^{\mathrm{4}}\left[\mathrm{1}-m\tanh \left(c_{\mathrm{3}}{T}^{\mathrm{6}}\right)\right],\end{array}$$

where σ is the Stefan–Boltzmann constant, and the emissivity coefficient is expressed as 1−mtanh (c₃T⁶). Such a term describes, in a simple yet effective way, the greenhouse effect by reducing infrared radiation losses. The values of the empirical functions $C\left(x\right),Q\left(x\right),b\left(x\right),z\left(x\right),k_{\mathrm{1}}\left(x\right),k_{\mathrm{2}}\left(x\right)$ at discrete latitudes and empirical parameters $c_{\mathrm{1}},c_{\mathrm{2}},c_{\mathrm{3}},c_{\mathrm{4}},c_{\mathrm{5}},\mathit{\sigma },m,T_m$ are taken from Ghil (1976), as modified in Bódai et al. (2015). The choice of the empirical functions and parameters are extensively discussed in Ghil (1976). Of course, one might reasonably wonder about the robustness of our modelling strategy. Indeed, a plethora of EBMs analogous to the one described here have been presented in the literature, where slightly different parameterisations for the diffusion operator, for the albedo, and for the greenhouse effect are introduced. Such models are in fundamental agreement, both in terms of physical (Ghil, 1981; North et al., 1981; North, 1990; North and Stevens, 2006) and mathematical properties (Hetzer, 1990; Díaz et al., 1997; Kaper and Engler, 2013; Bensid and Díaz, 2019).

Figure 1 (a) Stationary solutions T_W(x), T_SB(x), and T_M(x) in Kelvins (K) of the zonally averaged energy balance model Eq. (1). (b) Bifurcation diagram of the average temperature $\stackrel{\mathrm{‾}}{T}$ as a function of control parameter μ.

In this study, we consider μ=1.05. For this value of μ, two stable asymptotic states – the W and the SB states – co-exist (see Fig. 1b). Indeed, a codimension of one manifold separates the basins of attraction of the W and SB states. We refer to D^W (D^SB) as the basin of attraction of the W (SB state). We refer to B as the basin boundary, which includes a single edge state M. Therefore, the system has three stationary solutions T_W(x), T_SB(x), and T_M(x) for the W, SB, and M state, respectively, as shown in Fig. 1a. In Ghil (1976), the three stationary solutions were obtained by equating T_t to 0, and it was shown, through linear stability analysis, that the stationary solutions T_W and T_SB are stable, while T_M is unstable. In Bódai et al. (2015) the unstable solution T_M was constructed using a modified version of the edge-tracking algorithm (Skufca et al., 2006).

Following previous studies (Bódai et al., 2015; Lucarini and Bódai, 2019, 2020; Margazoglou et al., 2021), when visualising our results, we apply a coarse graining to the phase space of the model. In what follows, we perform a projection on the plane spanned by the spatially averaged temperature $\stackrel{\mathrm{‾}}{T}$ and the averaged Equator minus the poles' temperature difference ΔT, which is defined as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }\stackrel{\mathrm{‾}}{T}=\left[T\right(x,t){]}_{\mathrm{0}}^{\mathrm{1}},\text{(11)}& {\displaystyle }\mathrm{\Delta }T=\left[T\right(x,t){]}_{\mathrm{0}}^{\mathrm{1}/\mathrm{3}}-[T(x,t){]}_{\mathrm{1}/\mathrm{3}}^{\mathrm{1}},\text{ where}\text{(12)}& {\displaystyle }\left[T\right(x,t){]}_{x_l}^{x_h}={\displaystyle \frac{{\int }_{x_l}^{x_h}\cos (\mathit{\pi }x/\mathrm{2})T(x,t)\mathrm{d}x}{{\int }_{x_l}^{x_h}\cos (\mathit{\pi }x/\mathrm{2})\mathrm{d}x}}.\end{array}$$

Such a representation allows for a minimal yet still physically relevant description of the system. Indeed, changes in the energy budget of the system (warming versus cooling) are, to a first approximation, related to variations in $\stackrel{\mathrm{‾}}{T}$, while the large-scale energy transport performed by the geophysical fluids is controlled by ΔT. The boundary between high and low latitude in Eq. (11) is established at $x=\pm \mathrm{1}/\mathrm{3}$, i.e. at 30^∘ N/S. Additionally, in some visualisations, we consider, as a third coordinate, the fraction of the surface with a below-freezing temperature (therefore, we expect 1 for global glaciation and 0 for no ice). We refer to this variable as I, and it is an attempt to extract an observable that resembles the sea ice percentage of the Earth's surface. Thus, the stationary solutions T_W(x), T_SB(x), and T_M(x), in terms of ΔT and $\stackrel{\mathrm{‾}}{T}$, correspond to ΔT_W = 16 K, ΔT_SB = 8.3 K, ΔT_M = 17.5 K, ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{W}}$ = 297.7 K, ${\stackrel{\mathrm{‾}}{T}}_{\text{SB}}$ = 235.1 K, ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{M}}$ = 258 K, I_W=0.2, I_SB=1, and I_M=1.

3 Background and methodology

3.1 Stochastic energy balance model

In order to analyse the influence of random perturbations on the deterministic dynamics of the climate model described in Sect. 2, we perturb the relative intensity μ of the solar irradiance by including a symmetric α-stable Lévy process and rewrite Eq. (1) in the form of the following stochastic partial differential equation (SPDE):

$$\begin{array}{}\text{(13)}& \begin{array}{rl}C\left(x\right){\mathcal{T}}_t=& D_I\left(x,\mathcal{T},{\mathcal{T}}_x,{\mathcal{T}}_{xx}\right)\\ & +D_{\text{II}}(x,\mathcal{T})\left(\mathrm{1}+\mathit{\epsilon }/\mathit{\mu }{\dot{L}}^{\mathit{\alpha }}\left(t\right)\right)-D_{\text{III}}\left(\mathcal{T}\right),\end{array}\end{array}$$

where the boundary and initial conditions defined by Eq. (2) apply to the stochastic temperature field 𝒯. Here the parameter ε>0 controls the noise intensity, and (L^α(t)_t≥0) is a symmetric α-stable process defined in Appendix A. We consider symmetric processes because we want to have a simple mathematical model allowing for transitions in both the SB→W and the W direction →SB direction. Instead, a strongly skewed process would have made it very hard to explore the full phase space because a lack of symmetry would invariably favour one of the two transitions. As mentioned before, we refer to the Lévy case if the stability parameter $\mathit{\alpha }\in (\mathrm{0},\mathrm{2})$, so that we consider a jump process. We recall that the jumps become more frequent and less intense as α increases.

We define $\dot{\mathcal{L}}\left(t\right)=\mathcal{Q}\left(x\right)[\mathrm{1}-{\mathit{\alpha }}_{\mathrm{a}}(x,\mathcal{T}\left)\right]{\dot{L}}^{\mathit{\alpha }}\left(t\right)$, as the generalised derivative of a stochastic process in a suitably defined functional space. Equation (13) features multiplicative noise. The research interest on this type of SPDE (Doering, 1987; Peszat and Zabczyk, 2007; Duan and Wang, 2014; Alharbi, 2021) is mainly focused on defining weak, strong, mild, and martingale solutions, in specifying under which conditions these solutions exist and are unique, and in constructing numerical approximation schemes for the solutions (Davie and Gaines, 2000; Cialenco et al., 2012; Burrage and Lythe, 2014; Jentzen and Kloeden, 2009; Kloeden and Shott, 2001), among other aspects.

First, let us define the concept of a mild solution in this context. Let $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ be a given complete probability space and $H(\Vert \cdot \Vert ,\langle \cdot ,\cdot \rangle )$ a separable Hilbert space with a norm $\Vert \cdot \Vert$ and inner product $\langle \cdot ,\cdot \rangle$. Equation (13) can be rewritten in the more general form, as follows:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\mathcal{T}}_t=A\left(x\right)\left[E\right(x,\mathcal{T}){\mathcal{T}}_x{]}_x+F(x,\mathcal{T})+\mathit{\epsilon }G(x,\mathcal{T}\left){\dot{L}}^{\mathit{\alpha }}\right(t),\\ & {\mathcal{T}}_x(-\mathrm{1},t)={\mathcal{T}}_x(\mathrm{1},t)=\mathrm{0},\\ & \mathcal{T}(x,\mathrm{0})=T_{\mathrm{0}}\left(x\right),\end{array}\end{array}$$

where $A,E,F,G$ are Lipschitz functions defined on $[-\mathrm{1},\mathrm{1}]\times H$ and $G(x,\mathcal{T}){\dot{L}}^{\mathit{\alpha }}\left(t\right)=\dot{\mathcal{L}}\left(t\right)$. Under certain assumptions (Yagi, 2010), the problem (Eq. 14) is formulated as a Cauchy problem, for which a local mild solution, a progressively measurable process 𝒯(t), for all $t\in [\mathrm{0},t_F]$ and T₀∈H has the following integral representation:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathcal{T}\left(t\right)=& \mathrm{\Psi }\left(t\right)T_{\mathrm{0}}+\underset{\mathrm{0}}{\overset{t}{\int }}\mathrm{\Psi }(t-s)\mathrm{{\rm Y}}\left(\mathcal{T}\right(s\left)\right)\mathrm{d}s\\ & +\mathit{\epsilon }\underset{\mathrm{0}}{\overset{t}{\int }}\mathrm{\Psi }(t-s)G\left(\mathcal{T}\right(s\left)\right)\mathrm{d}\mathit{\beta }\\ & +\mathit{\epsilon }\underset{\mathrm{0}}{\overset{t}{\int }}\mathrm{\Psi }(t-s)G\left(\mathcal{T}\right(s\left)\right)\mathrm{d}\mathit{\gamma },\end{array}\end{array}$$

where the dependence on x is kept implicit, and β (γ) is the Poisson random measure (compensated Poisson random measure) defined through Lévy–Itô decomposition. Instead, Ψ(t) with t⩾0 is the evolution operator having the generalised semigroup property for the family of sector operators with the bounded inverses, and $\mathrm{{\rm Y}}\left(\mathcal{T}\right)=\mathcal{T}+F(x,\mathcal{T})$, 𝒯∈H is a nonlinear operator, which we assume to be Lipschitz continuous. Following the abstract theory presented in Yagi (2010), under certain structural assumptions for the operators Ψ and Υ and for the functional space, one can prove that the solution Eq. (15) is the unique local mild solution of Eq. (14).

As mentioned above, things are radically different for the special case α=2, which corresponds to Gaussian noise. In this case, we revisit Eq. (14), and we define ${\dot{L}}^{\mathit{\alpha }=\mathrm{2}}\left(t\right)=\dot{W}\left(t\right)$, where (W(t)_t≥0) is a Wiener process. We then define $\dot{\mathcal{W}}\left(t\right)=G(x,\mathcal{T})\dot{W}\left(t\right)$ as the generalised derivative of a Wiener process in a suitably defined functional space.

3.2 Noise-induced transitions: mean escape times

By incorporating stochastic forcing into the system, its long-time dynamics change significantly, allowing transitions between the competing basins. This dynamical behaviour is called metastability and is graphically captured by Fig. 2, where, in Fig. 2a and b, a typical spatiotemporal evolution of the temperature field is shown for the stability parameters α=0.5 and α=1.5, respectively. Instead, in Fig. 2c and d, the temporal evolution of the global temperature $\stackrel{\mathrm{‾}}{\mathcal{T}}$ and of the averaged Equator and poles' temperature difference Δ𝒯 (as defined in Eqs. 10–11) is correspondingly shown. In what follows, we investigate the time statistics and the paths of the transitions between such basins.

In a complete probability space $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ we define the first exit time τ_x of a cádlág mild solution $\mathcal{T}(\cdot ;x)$ of Eq. (13), starting at the $x\in D^{\mathrm{W}/\text{SB}}$ domain of a W/SB climate stable state as follows:

$$\begin{array}{}\text{(16)}& {\mathit{\tau }}_x\left(\mathit{\omega }\right)=inf\left\{t>\mathrm{0}\left|{\mathcal{T}}_t\right(\mathit{\omega },x)\notin D^{\mathrm{W}/\text{SB}}\right\},\mathit{\omega }\in \mathrm{\Omega },x\in H.\end{array}$$

The mean escape time is then expressed by 𝔼[τ_x(ω)]. In the case of the infinite dimensional multistable reaction–diffusion system described by Chafee–Infante equation under the influence of additive infinite-dimensional α-stable Lévy noise – $\mathit{\alpha }\in (\mathrm{0},\mathrm{2})$ – it was shown (Debussche et al., 2013) that, in the weak noise limit, ε→0 the mean escape time from one of the competing basins of attraction increases as ε^−α. In such a limit, the jump diffusion system reduces the Markov chain to a finite state, with values in the set of stable states. Details of this method are given in Appendix B. Similar results have been obtained for bistable one-dimensional stochastic differential equations (SDEs; Imkeller and Pavlyukevich, 2006a, b). The basic reason behind this result is that, in order to study the transitions between the competing basins of attraction, one can treat the Lévy noise as a compound Poisson process, where jumps arrive randomly, according to a Poisson process, and the size of the jumps x is given by a stochastic process that obeys a specified probability distribution. For a symmetric α-stable Lévy process, such a distribution asymptotically decreases as $|x{|}^{-\mathrm{1}-\mathit{\alpha }}$, as discussed in Appendix A. Let us assume that positive values of x bring the state of the system closer to the basin boundary (as in the case of positive fluctuations of the solar irradiance when studying escapes from the SB state). Assuming a simple geometry for the basin boundary, we find that a transition takes place when an event larger than a critical value x_crit>0 is realised. The probability of such an event scales with $x_{\text{crit}}^{-\mathit{\alpha }}$. A similar argument applies when considering transitions triggered by negative fluctuations of the stochastic variable. Small-size events, which occur frequently and correspond to the non-occurrence of jumps, do not actually play any relevant role in determining the transitions, while they are responsible for the variability within each basin of attraction.

We now consider the case α=2. While the corresponding finite dimensional problem is thoroughly documented in the literature (Freidlin and Wentzell, 1984) and has been applied in a similar context by some of the authors (Lucarini and Bódai, 2019, 2020; Ghil and Lucarini, 2020; Margazoglou et al., 2021), the treatment of infinite dimensional SDEs driven by an infinite dimensional Wiener process via the Freidlin–Wentzell theory requires further extension. In the present context, we refer to Budhiraja and Dupuis (2000) and Budhiraja et al. (2008) and references therein, where the problem of an infinite dimensional reaction–diffusion equation driven by an infinite dimensional Wiener process has been addressed.

We assume that steady-state conditions and ergodicity are met, and we also assume that the analysing system is bistable and a unique edge state is present at the basin boundary, as in the case studied here. In the case of Gaussian noise, transitions between the competing basins of attraction are not determined by a single event as in the $\mathrm{0}<\mathit{\alpha }<\mathrm{2}$ case but, instead, occur as a result of very unlikely combinations of subsequent realisations of the stochastic variable acting as a forcing. In the weak noise limit, the transitions occur according to the least unlikely (yet very unlikely) chain of events, whose probability is described using a large deviation law (Varadhan et al., 1985). One finds that the mean escape time from either basin of attraction decreases exponentially with increasing noise intensity ε and is given by a generalised Kramers' law, as follows:

$$\begin{array}{}\text{(17)}& \mathbb{E}\left[{\mathit{\tau }}_{\mathrm{W}/\text{SB}}\right(\mathit{\epsilon }\left)\right]\approx \exp \left({\displaystyle \frac{\mathrm{2}\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{W}\to \mathrm{M}/\text{SB}\to \mathrm{M}}\left(\mathcal{T}\right)}{{\mathit{\epsilon }}^{\mathrm{2}}}}\right),\end{array}$$

where $\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{W}\to \mathrm{M}}={\mathrm{\Phi }}_{\mathrm{M}}\left(\mathcal{T}\right)-{\mathrm{\Phi }}_{\mathrm{W}}\left(\mathcal{T}\right)$ is the height of the quasi-potential barrier in the W attractor; correspondingly, $\mathrm{\Delta }{\mathrm{\Phi }}_{\text{SB}\to \mathrm{M}}\left(\mathcal{T}\right)={\mathrm{\Phi }}_{\mathrm{M}}\left(\mathcal{T}\right)-{\mathrm{\Phi }}_{\text{SB}}\left(\mathcal{T}\right)$ is the height of the quasi-potential barrier in the SB attractor, and Φ is the Graham quasi-potential mentioned above (Graham, 1987; Graham et al., 1991).

3.3 Noise-induced transitions: most probable transition paths

In the weak noise limit, the most probable path to escape an attractor is defined by a class of trajectories named instantons (Grafke et al., 2015, 2017; Bouchet et al., 2016; Grafke and Vanden-Eijnden, 2019) or maximum likelihood escape paths (Lu and Duan, 2020; Dai et al., 2020; Hu and Duan, 2020; Zheng et al., 2020). However, note that different noise laws result into possibly radically different instantonic trajectories (Dai et al., 2020; Zheng et al., 2020).

In our case, the theory indicates that, if the stochastic forcing is Gaussian, under a rather general hypothesis, the instanton will connect the attractor W/SB with the edge state M, which then acts as gateway for noise-induced transitions. Once the quasi-potential barrier is overcome, a free-fall relaxation trajectory links M with the competing attractor SB/W. For equilibrium systems, (e.g. for gradient flows), where a detailed balance is achieved, the relaxation and instantonic trajectories within the same basin of attraction are identical. On the contrary, for non-equilibrium systems, the relaxation and instantonic trajectories will differ and will only meet at the attractor. (See a detailed discussion of this aspect and of the dynamical interpretation of the quasi-potential Φ in Lucarini and Bódai, 2020, and Margazoglou et al., 2021). Instead, if the noise is of Lévy type, the theory formulated for simpler equations suggests that the instanton will connect the attractor with a region on the basin boundary that is the nearest, in the phase space, to the attractor, as the concept of the quasi-potential is immaterial (Imkeller and Pavlyukevich, 2006a, b).

In general, the maximum likelihood transition trajectory 𝒯_M(t) can be defined (Zheng et al., 2020; Lu and Duan, 2020) as a set of system states at each time moment $t\in [\mathrm{0},t_f]$ that maximises the conditional probability density function $p(.|.;.)$ of the passage from the origin stable state ${\mathit{\varphi }}^{\mathrm{W}/\text{SB}}$ to the destination stable state ${\mathit{\varphi }}^{\text{SB}/\mathrm{W}}$ and is expressed as follows:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\mathcal{T}}_{\mathrm{M}}\left(t\right)\\ & =\arg \underset{x}{max}\left[p\left(\mathcal{T}\right(t)=x|\mathcal{T}\left(\mathrm{0}\right)=x_{\mathrm{0}};\mathcal{T}\left(t_f\right)=x_f)\right]\\ & ={\displaystyle \frac{p\left(\mathcal{T}\left(t_f\right)=x_f\left|\mathcal{T}\right(t)=x\right)\cdot p\left(\mathcal{T}\left(t\right)=x\left|\mathcal{T}\right(\mathrm{0})=x_{\mathrm{0}}\right)}{p\left(\mathcal{T}\right(t_f)=x_f|\mathcal{T}\left(\mathrm{0}\right)=x_{\mathrm{0}})}},\end{array}\end{array}$$

where x₀ (x_f) belongs to the basin of attraction $D^{\mathrm{W}/\text{SB}}$ ($D^{\text{SB}/\mathrm{W}}$), and $p(.|.)$ is the probability density function evolving according to the Fokker–Planck equation (Risken, 1996). This method is applicable either if efficient numerical algorithms are available to solve the Fokker–Planck equation associated to the studied stochastically driven system, or, empirically, when considering a large ensemble of simulations. Note that this is not an asymptotic approach, i.e. it does not require a weak noise limit ε→0 and is applicable for systems with either Gaussian or non-Gaussian noise. Yet, in the weak noise limit, the definition (Eq. 18) leads to constructing the optimal transition paths described above.

In the following section, for practical purposes, we construct such optimal transition path in the coarse-grained 2D phase space ($\stackrel{\mathrm{‾}}{\mathcal{T}},\mathrm{\Delta }\mathcal{T}$) and 3D phase space ($\stackrel{\mathrm{‾}}{\mathcal{T}},\mathrm{\Delta }\mathcal{T},\mathcal{I})$ of the variables defined in Sect. 2 by averaging the ensemble of transitions connecting the two competing states in the weak noise limit.

Figure 2 The metastable behaviour of the solution path of a stochastic energy balance model (Eq. 13), for ε=0.04, T₀ = 300 K, and $t\in (\mathrm{0},\mathrm{300})$ years, for (a) α=0.5, (b) α=1.5, and the (c) temperature average $\stackrel{\mathrm{‾}}{\mathcal{T}}$ and (d) temperature contrast at low and high latitudes. Δ𝒯, for ε=0.01, $\mathit{\alpha }=\mathrm{0.5},$andα=1.5. Red, green, or blue dash-dotted lines portray the stationary climate states of ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{W}}$/${\stackrel{\mathrm{‾}}{T}}_{\mathrm{M}}$/${\stackrel{\mathrm{‾}}{T}}_{\text{SB}}$, respectively.

3.4 Numerical methods

We solve Eq. (13) through the MATLAB pdepe function, which is well suited for solving 1D parabolic and elliptic PDEs. We discretise the 1D space with a regular grid of 201 grid points, following Bódai et al. (2015).

The time span of integration $t\in [\mathrm{0},T_f]$, varies for different cases, with $T_f\in ({\mathrm{10}}^{\mathrm{5}},\mathrm{15}\times {\mathrm{10}}^{\mathrm{5}})$ years, with a time stepping of 1 year. Each year, we consider a different value for the relative solar irradiance by extracting a random number Z_j (see Eq. 19). To simulate the stochastic noise term εL^α(t), which is added in the parameter μ in Eq. (13), we use the recursive algorithm from Duan (2015). The process values $L^{\mathit{\alpha }}\left(t_{\mathrm{1}}\right),\mathrm{\dots },L^{\mathit{\alpha }}\left(t_N\right)$ at each moment t_j, j∈ℕ are obtained via the following:

$$\begin{array}{}\text{(19)}& L^{\mathit{\alpha }}\left(t_j\right)=L^{\mathit{\alpha }}\left(t_{j-\mathrm{1}}\right)+(t_j-t_{j-\mathrm{1}}{)}^{\frac{\mathrm{1}}{\mathit{\alpha }}}{Z}_j,j=\mathrm{1},\mathrm{\dots },N,\end{array}$$

where the second term is an independent increment, and Z_j are the independent standard symmetric α-stable random numbers generated by an algorithm in Weron and Weron (1995, see also Grafke et al., 2015, for a detailed explanation of the steps above). For illustrative reasons, some sample solutions of Eq. (13) for different values of α are shown in Fig. 2a and b.

For the numerical simulations discussed below, we consider $\mathit{\alpha }=(\mathrm{0.5},\mathrm{1.0},\mathrm{1.5},\mathrm{2})$ and $\mathit{\epsilon }\in (\mathrm{0.0001},\mathrm{0.3})$. We select ε in such a way that the noise intensity is strong enough to induce at least an order of 10 transition, given our constraints in the time length of the simulations, and weak enough that we are not far from the weak noise limit, where the scaling laws discussed above apply and transitions paths are well organised. Our simulations are performed by taking the Itô interpretation for the stochastic equations.

We remark that, when we consider Lévy noise, it does happen that, for some years, the solar irradiance has negative values. Of course these conditions bear no physical relevance, and are a necessary result of considering unbounded noise. Nonetheless, we have allowed for this to occur in our simulations in order to be able to stick to the desired mathematical framework. We remind the reader that this study does not aim at capturing, with any high degree of realism, the description of the actual evolution of climate. At any rate, as can be understood from the discussion below in Sect. 4.2.2 and from what is reported in Appendix C, were we to consider longer-lasting (e.g. 2 years vs. 1 year) fluctuations of the solar irradiance, a satisfactory exploration of the transitions between the competing W and SB states would be possible with a greatly reduced occurrence of such unphysical events, and the basic reason for this being the presence of a larger factor $(t_j-t_{j-\mathrm{1}}{)}^{\frac{\mathrm{1}}{\mathit{\alpha }}}$ in Eq. (19).

4 Results and discussion

In what follows, we aim at addressing three main questions: (1) what are the temporal statistics of the SB→W and W→SB transitions? (2) What are the typical transition pathways? (3) What are the fundamental differences between transitions caused by Gaussian vs. Lévy noise? A summary of the results of the numerical simulations is given in Table D1 in Appendix D, including the sample size, i.e. number of transitions, point estimates for mean escape time, and their 0.95 confidence intervals for exits from both the W and SB basins. See the data availability section for information on how to access the supplement (Lucarini et al., 2022), which contains the raw data produced in this study and some illustrative animations portraying noise-induced transitions between the two competing metastable states.

4.1 Mean escape time

Our analysis confirms that there is a fundamental dichotomy in the statistics of mean escape times between Lévy noise and Gaussian noise-induced transitions.

Figure 3a shows the dependence of the mean escape time from either attractor on ε and α for the Lévy case. The red circles (blue squares) correspond to escapes from the W (SB) basin (see Lucarini et al., 2022, for additional details). The scaling $\propto {\mathit{\epsilon }}^{-\mathit{\alpha }}$ presented in Eq. (B6) is shown by the dotted black line for each value of α. We also portray the best power law fit of the mean residence time with respect to ε for each value of α; the confidence intervals of the exponent are shown in Table 1. Our empirical results seem to indicate, at least in this case, an agreement with the ε^−α scaling presented and discussed earlier in the paper. This points at the possibility that the ε^−α scaling might apply in more general conditions than what has been, as of yet, rigorously proven, and this is specifically the case when multiplicative Lévy noise is considered. The stochastically perturbed trajectories forced by Lévy noise consist of jumps, and the probability of occurrence of a high jump, which can trigger the escape from the reference basin of attraction, is polynomially small in noise intensity ε.

The Gaussian case – where no jumps are present – is portrayed in Fig. 3b. We show, in semi-logarithmic scale, the mean residence time versus $\mathrm{1}/{\mathit{\epsilon }}^{\mathrm{2}}$. We perform a successful linear fit of the logarithm of the mean residence time in either attractor versus $\mathrm{1}/{\mathit{\epsilon }}^{\mathrm{2}}$, and using Eq. (17), we obtain an estimate of the local quasi-potential barrier $\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{W}/\text{SB}\to \mathrm{M}}$, which is half of the slope of the corresponding straight lines of the linear fit (see the last column of Table 1). We conclude that, for μ=1.05, the local minimum of Φ corresponding to the W state is deeper than the one corresponding to the SB state.

Figure 3 Estimates of the mean escape time 𝔼(τ) (in years) from the W (blue circles) and SB (red squares) states, as a function of the noise intensity ε. (a) Lévy noise for α=0.5, 1.0, and 1.5, with the dotted line being the corresponding prediction from Eq. (B6), while the straight green (SB) and dashed black (W) lines are the fittings of Eq. (B6) of the relevant dataset. (b) Gaussian noise, with straight green (SB) and dashed black (W) lines being the fit of Eq. (17), is shown.

Table 1 Estimates of the exponent α via the fitting of Eq. (B6) for the Lévy case (three first columns) and of the energy barrier $\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{W}/\text{SB}\to \mathrm{M}}$ via the fitting of Eq. (17) for the Gaussian case (last column). In parenthesis, the estimated error of the last digit is shown.

4.2 Escape paths for the noise-induced transitions

We now explore the geometry of the transition paths associated with the metastable behaviour of the system. We first discuss the case of Gaussian noise because it is indeed more familiar and more extensively studied.

4.2.1 Gaussian noise

We estimate the transition paths by averaging among the escape plus relaxation trajectories using the run performed with the weakest noise (see Table D1). We first perform our analysis in the 2D-projected state space defined by $(\stackrel{\mathrm{‾}}{\mathcal{T}},\mathrm{\Delta }\mathcal{T})$. We prescribe two small circular-shaped regions enclosing the two deterministic attractors and search the time series of the portions of the whole trajectory that leave one of these regions and reach the other one. This creates two subsets of our full dataset from which we build a 2D histogram for each of the SB→W and W→SB transitions in the projected space. We then estimate the most probable transition paths by finding, for each bin value of $\stackrel{\mathrm{‾}}{\mathcal{T}}$, the peak of the histogram in the Δ𝒯 direction. The distributions are very peaked, and almost indistinguishable estimates for the instantonic and relaxation trajectories are obtained when computing the average of Δ𝒯 according to the 2D histogram conditional on the value of $\stackrel{\mathrm{‾}}{\mathcal{T}}$.

Figure 4 (a) Invariant measure in the 2D-projected state space defined by $\left(\stackrel{\mathrm{‾}}{\mathcal{T}},\mathrm{\Delta }\mathcal{T}\right)$. The coloured points indicate the deterministic attractors of the SB (blue), M (green), and W (red) states, and the blue (red) line is the stochastically averaged transition paths for the W→SB (SB→W) transitions. Dashed (solid) lines are the relaxation (instantonic) trajectories. The arrows show the direction of transitions. The top left inset shows that the dark blue (red) contours portray the ensembles of the transition paths between W→SB (SB→W). Here, the system is driven by Gaussian noise with ε=0.14. (b) The invariant measure and most probable transition paths (W→SB in blue and SB→W in red) in the 3D-projected state space are defined by $\left(\stackrel{\mathrm{‾}}{\mathcal{T}},\mathrm{\Delta }\mathcal{T},\mathcal{I}\right)$. The darker brown shading indicates a higher probability density for the corresponding isosurface. The 2D projections on the $\left(\stackrel{\mathrm{‾}}{\mathcal{T}},\mathcal{I}\right)$, and (Δ𝒯,ℐ) planes are shown. The location of the M state is given by a pink square. Here, the system is driven by Gaussian noise, with ε=0.16.

In the background of Fig. 4a, we show the empirical estimate of the invariant measure in the 2D-projected state space defined by $(\stackrel{\mathrm{‾}}{\mathcal{T}},\mathrm{\Delta }\mathcal{T})$. Additionally, we indicate the position of the deterministic attractors, where the blue (red) circle corresponds to the SB (W) state and of the M state (green square). In the inset of Fig. 4a, we present the ensemble of W→SB (SB→W) transitions as deep blue (red) contours. The most probable transition paths are shown in blue for the W→SB and in red for the SB→W. The instantonic portion of the blue (red) line is the one connecting the W (SB) attractor to the M state and is portrayed as a solid line, while the relaxation portion, connecting the M state with the SB (W) attractor, is portrayed as a dashed line. Within each basin of attraction, the instantonic and relaxation trajectories do not coincide, and, instead, only meet at the corresponding attractor and at the M state. This is particularly clear for the W state. The presence of such a loop, proving the existence of non-vanishing probability currents and the breakdown of detailed balance, is a signature of non-equilibrium dynamics, which was also observed in Margazoglou et al. (2021) and has, instead, gone undetected in Lucarini and Bódai (2019) and Lucarini and Bódai (2020). See Lucarini et al. (2022) for some illustrative simulations of the transitions.

Let us provide some physical interpretation of how the transitions occur. Looking at the SB→W most probable path, the escape includes a simultaneous increase in $\stackrel{\mathrm{‾}}{\mathcal{T}}$ and Δ𝒯. In practice, a SB→W transition takes place when, starting at the SB state, one has a (rare) sequence of positive anomalies in the fluctuating solar irradiance $\stackrel{\mathrm{̃}}{\mathit{\mu }}$, i.e. $\stackrel{\mathrm{̃}}{\mathit{\mu }}>\mathit{\mu }$. While the planet is warming globally, the Equator is warming faster than the poles, resulting in a positive rate $\dot{\mathrm{\Delta }\mathcal{T}}>\mathrm{0}$, because it receives, in relative and absolute terms, more incoming solar radiation. Considering that the Equator also in the SB state is warmer than the poles, the melting of the ice conducive to the transition occurs first at the Equator, with a subsequent decrease in the albedo in low latitudes. Once the system crosses the M state, and supposing that persistent $\stackrel{\mathrm{̃}}{\mathit{\mu }}<\mathit{\mu }$ do not appear at this stage, the system will relax towards the W state. The relaxation includes a consistent global warming of the planet but with a change in sign in the rate of $\dot{\mathrm{\Delta }\mathcal{T}}$, and a subsequent decrease in Δ𝒯, implying that as soon as the temperature at the Equator has risen enough, the poles will then warm at a faster pace because the ice–albedo effect kicks in.

The global freezing of the planet associated with the W→SB transition is qualitatively similar but not identical to the reverse SB→W process. Notice a considerable overlap of the transition paths ensembles in both basins of attraction, shown as red and blue contours in the inset of Fig. 4a. This implies the presence of less extreme non-equilibrium conditions compared to what was observed in Margazoglou et al. (2021), where the W→SB and SB→W transitions occurred through fundamentally different paths (see the discussion therein, especially regarding the role of the hydrological cycle).

Figure 4b presents the optimal transition paths W→SB and SB→W in a three-dimensional projection, where we add, as a third coordinate, the variable ℐ, which indicates the fraction of the surface that has subfreezing temperatures (𝒯 < 273.15 K). On the sides of the figure, two two-dimensional projections on the (𝒯,Δℐ) and on the (Δ𝒯,ℐ) planes are shown. Here, darker brown shadings indicate the higher density of points and the red and blue dots sample the highest probability for the SB→W and W→SB transitions paths, respectively. One could argue that the presence of an intersection between the SB→W and W→SB highest probability transition paths in Fig. 4a could have been a simple effect of 2D projection. Instead, we see here that the SB→W and W→SB most probable transition paths also cross in the 3D projection in a well-defined region, which indeed corresponds to the M state (pink square).

4.2.2 Lévy noise

There is scarcity of rigorous mathematical results regarding the weak noise limit of the transition paths between competing states in metastable stochastic systems forced by multiplicative Lévy noise. Indeed, the derivation of analytical results for this type of system largely remains an open problem. Recently, for stochastic partial differential equations with additive Lévy and Gaussian noise, the Onsager–Machlup action functional has been derived in Hu and Duan (2020), leading to a precise formulation of the most probable transition paths. Hence, we do not have solid mathematical results to interpret what we describe below, where, instead, we need to use heuristic arguments. As far as we know, this is the first attempt to estimate the most probable transition pathway between the metastable states in infinite stochastic systems with multiplicative pure Lévy process.

Figure 5 A two-dimensional projection of the invariant measure on $\left(\stackrel{\mathrm{‾}}{\mathcal{T}},\mathrm{\Delta }\mathcal{T}\right)$ for different choices of α for Lévy noise. (a) α=0.5 and ε=0.0001, (b) α=1 and ε=0.004, and (c) α=1.5 and ε=0.01. The blue (red) line corresponds to the W→SB (SB→W) most probable transition path, and the arrows show the direction of the transitions. The coloured points indicate the deterministic attractors of the SB (blue), M (green), and W (red) states. In the inset of the left top corner plot in panels (a–c), the dark blue (red) contours are the ensembles of the transition paths between W→SB (SB→W).

A striking feature in Fig. 5 is that the invariant measure and the structure of the most probable transition paths (SB→W and W→SB), in the weak noise limit, are fundamentally different between the Lévy case and the Gaussian one. The invariant measure is highly peaked (dark red in the colour scheme) in a small region around the deterministic attractors, as, most typically, the Lévy noise fluctuations of $\stackrel{\mathrm{̃}}{\mathit{\mu }}$ are very small. Additionally, the most probable transition paths depend very weakly on the chosen value for the stability parameter α. This suggests that the geometry of most probable path of transitions does not depend on the frequency and height of the Lévy diffusion jumps but rather on the qualitative fact that we are considering jump processes. Note that each panel of Fig. 5 is computed using data coming from the weakest noise considered for the corresponding value of α (see Table D1).

The W→SB most probable transition path is characterised by the simultaneous decrease in both $\stackrel{\mathrm{‾}}{\mathcal{T}}$ and Δ𝒯. This implies that the jump leads to a rapid and direct freezing of the whole planet. The stochastically averaged path crosses the basin boundary far from the M state. The most probable SB→W transition follows, instead, a path that is somewhat similar to the one found for the Gaussian case. We then argue that the closest region in the basin boundary to SB attractor is not too far from the M state. Further insight on difference between the Gaussian and Lévy case can be found by looking at the animations included in Lucarini et al. (2022).

Figure 6 Comparison between the supercritical and subcritical paths that are constricted, using singular perturbations and the ensemble of trajectories corresponding to Lévy noise-induced transitions. (a) SB→W transitions. (b) W→SB transitions. Singular perturbations are where the subcritical paths are depicted as thin solid and dashed lines. The supercritical paths, leading to transition, are depicted as thick solid and dashed lines. See the text for further details. Here, α=1.0 and ϵ=0.004.

Figure 7 Same as Fig. 6, respectively, using the same 3D projection, lines, and colour coding as in Fig. 4. The subcritical and supercritical paths are depicted as thin and thick solid lines, respectively. The corresponding dashed lines shown in Fig. 6 are not reported here, for simplicity. Here, α=1.0 and ϵ=0.004.

4.2.3 Lévy noise – singular perturbations

Based on what is discussed in Sect. 3.2, we expect that the transitions occur through the nearest region to the outgoing attractor in the basin boundary. We now try to clarify the properties of the most probable escape paths in the Lévy noise case by considering an additional set of simulations, taking inspiration from the edge-tracking algorithm (Skufca et al., 2006). The idea is to exploit the fact that large jumps drive the transitions in the Lévy noise case. Starting from the deterministic SB state, we apply, in Eq. (6), singular perturbations of the form $\mathit{\mu }\to \mathit{\mu }+\mathit{\kappa }\mathit{\delta }\left(t\right)$ and bracket the critical value ${\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W}}$, leading to a transition to the W state as ${\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W}}\in [{\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W},\mathrm{s}},{\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W},\mathrm{u}}]$, where the simulation performed with the supercritical (subcritical) value of $\mathit{\kappa }={\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W},\mathrm{u}}$ ($\mathit{\kappa }={\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W},\mathrm{s}}$) asymptotically reaches the competing (comes back to the initial) steady state. We obtain ${\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W},\mathrm{s}}\approx \mathrm{1.149}$  y and ${\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W},\mathrm{u}}\approx \mathrm{1.1492}$  y. Starting from the W state, we follow a similar procedure and find ${\mathit{\kappa }}_{\text{crit}}^{\mathrm{W}\to \text{SB}}\in [{\mathit{\kappa }}_{\text{crit}}^{\mathrm{W}\to \text{SB},\mathrm{u}},{\mathit{\kappa }}_{\text{crit}}^{\mathrm{W}\to \text{SB},\mathrm{s}}]$. We obtain ${\mathit{\kappa }}_{\text{crit}}^{\mathrm{W}\to \text{SB},\mathrm{s}}\approx -\mathrm{1.3458}$  y and ${\mathit{\kappa }}_{\text{crit}}^{\mathrm{W}\to \text{SB},\mathrm{u}}\approx -\mathrm{1.346}$  y. In both cases, the value of κ of the supercritical and subcritical paths differs by δκ≈0.0002  y.

The projections on the 2D phase space spanned by (𝒯,Δ𝒯) of the supercritical and subcritical paths corresponding to the SB→W (W→SB) transition are shown in Fig. 6a (Fig. 6b), using the thick and thin black dashed lines, respectively. The basin boundary is indicated by a cyan (magenta) line for the SB→W (W→SB) transition. The steps to estimate the basin boundary are presented in Appendix C. Note that we are using, as background, the invariant measure and the subset of transitions referring to Lévy noise simulations performed using α=1.0. Nonetheless, what is discussed below would apply equally well had we chosen to consider as background, instead, data coming from the simulation performed with α=0.5 or α=1.5. Indeed, for the link we propose between transitions due to the Lévy noise and the case of singularly perturbed trajectories, what matters depends on the discontinuous nature of the Lévy noise.

In Fig. 6a (Fig. 6b) the supercritical and subcritical paths are superimposed on the ensemble of the trajectories of the SB→W (W→SB) transitions due to Lévy noise. The lines are better visible in the insets. By construction, after the perturbation is applied, the supercritical and subcritical orbits are close to the basin boundary. Hence, they are attracted towards the M state before being repelled towards the final asymptotic state. For comparison, we also portray, for both the SB→W and the W→SB case, an additional pair of supercritical and subcritical paths that are constructed using values of κ that differ by δκ=0.3  y, which is a much larger difference than the one mentioned previously (for the dashed lines). The paths depicted as thick (thin) solid lines cross (do not cross) the basin boundary. When looking at the W→SB transitions due to Lévy noise, we understand that, if the perturbation sends the orbit near the basin boundary, then the subsequent evolution of the system follows the supercritical paths. Of course, the Lévy perturbation often overshoots the basin boundary. In this case, after the transition, the orbit is not necessarily attracted towards the M state, whereas it converges directly to the final SB state. The signature of the attracting influence of the M state persists in the stochastically averaged transition trajectory. Note the bending towards higher values of Δ𝒯 before the eventual convergence to the SB state. In the case of the SB→W transitions, a similarly good correspondence between the supercritical path and the stochastically averaged transition trajectory is found.

Further support to the viewpoint proposed here is given by Fig. 7a and b, which are constructed along the lines of Fig. 4b and portray the supercritical and subcritical paths and the stochastically averaged transition trajectory realised via Lévy noise, with α=1.0 for the SB→W and W→SB case, respectively.

The transitions shown in Figs. 6 and 7 have been obtained by considering a discrete approximation of Dirac's δ, where the forcing acts at a constant value for 1  y. Specifically, Dirac's δ(t) is approximated as Δ_τ(t), where ${\mathrm{\Delta }}_{\mathit{\tau }}\left(t\right)=\mathrm{1}/\mathit{\tau }$, if $\mathrm{0}<t<\mathit{\tau }$, and vanishes elsewhere. The results are virtually unchanged if one considers τ<1  y because the main dynamical processes of the re-equilibration of the system act on longer timescales. The effect of the negative feedbacks of the system starts to become apparent when considering slower perturbations, lasting 2 or more years. Indeed, the resilience of the system to transitions is reduced when, ceteris paribus, faster perturbations are considered (see the analyses in this direction dealing with stability of the large-scale ocean circulation; Stocker and Schmittner, 1997; Lucarini et al., 2005, 2007; Alkhayuon et al., 2019). Nonetheless, also in this case, the agreement with the results presented in Figs. 6 and 7 is considerable. We remark that considering longer-lasting perturbations allows one to observe W→SB transitions without using, at any time, (unphysical) negative values for the solar irradiance. This is reassuring in terms of the robustness and of the physical sense of our results. Further details on the impact of considering different values of τ are reported in Appendix C.

5 Conclusions

It is a well-known that, as a result of the competition between the Boltzmann stabilising feedback and the ice–albedo destabilising feedback, under current astronomical and astrophysical conditions, the climate system is multistable, as at least two competing and distinct climates are present, i.e. the W and the SB. More recent investigations indicate that the partition of the phase space of the climate system might be more complex, as more than two asymptotic states might be present, some of them, possibly, associated with small basins of attraction.

For deterministic multistable systems, the asymptotic state of an orbit depends uniquely on the initial condition, and, specifically, on which basin of attraction it belongs to. The presence of stochastic forcing allows for transitions to occur between competing basins, thus giving rise to the phenomenon of metastability. Gaussian noise as a source of stochastic perturbations has been widely studied by the scientific community in recent years and provided very fruitful insight into the multiscale nature of the climatic time series. However, it has become apparent that more general classes of α-stable Lévy noise laws might also be suitable for modelling the observed climatic phenomena. In this regard, it is important to achieve a deeper understanding of the possible noise-induced transitions between competing stable climate states under α-stable Lévy perturbations and compare them with the Gaussian case.

As a starting point in this direction, we have studied the influence of different noise laws on the metastability properties of the randomly forced Ghil–Sellers EBM, which is governed by a nonlinear, parabolic, reaction–diffusion PDE. In the deterministic version of the model, we have three steady-state solutions, i.e. two stable, attractive climate states and one unstable saddle, corresponding to the edge state. The stable states correspond to the well-known W and SB climates. There is a fundamental dichotomy in the properties of the noise-induced transitions determined by whether we consider a stochastic forcing of intensity ε with a Gaussian versus an α-stable Lévy noise law. Note that, instead, the spatial structure of the noise is unchanged. This indicates that the phenomenology associated with the metastable behaviour depends critically on the choice of the noise law. Not many studies have investigated, numerically or through mathematical theory, the properties of transitions in metastable systems driven by multiplicative Lévy noise, as done here.

First, in the weak noise limit ε→0, the mean residence times inside either competing basin of attraction for diffusions driven by Gaussian vs. Lévy noise have a fundamentally different dependence on ε. Our results show that the logarithm of the mean residence time for Gaussian diffusions scales with ε⁻², while, instead, a much weaker dependence is found for the Lévy case. Indeed, we find that the mean residence time is proportional to ε^−α, where α is the stability parameter of the noise law. This result is in agreement with what has been proven in some special cases for additive Lévy noise and might indicate that these scaling properties are more general than usually assumed. We propose a simple argument based on approximating the Lévy noise as a composed Poisson process to support the applicability of the result in general circumstances, but, clearly, detailed mathematical investigations in this direction are needed.

Second, the results obtained for the most probable transition paths confirm that, in the weak noise limit, escapes from basins of attraction driven by Gaussian noise take place through the edge state. Additionally, instantonic and relaxation portions within each basin of attraction are clearly distinct, indicating nonequilibrium conditions that are, yet, qualitatively similar. In turn, Lévy diffusions leave the basin through the boundary region closest to the outgoing attractor, which seems to be the vicinity of the edge state when the thawing transition is considered. The freezing transition, instead, proceeds along a path that is fundamentally different. Finally, the most probable transition paths for the Lévy case appear to depend very weakly on the value of the stability parameter α, but seem, instead, to be determined by the nature of the Lévy noise of being a jump process. Indeed, we suggest that these properties can be better understood by considering that, to a first approximation, the transitions due to the Lévy diffusion correspond to supercritical paths associated with Dirac's δ-like singular perturbations to the solar irradiance. This viewpoint seems of general relevance to other problems where Lévy noise is responsible for exciting transitions between competing metastable states.

Our findings provide strong evidence that choosing noise laws other than Gaussian leads to fundamental changes in the metastability properties of a system, both in terms of statistics of the transitions between competing basins of attraction and most probable paths for such transitions. Leaving the door open for general noise laws might be relevant, both for interpreting observational data and for performing modelling exercises for the climate system and complex systems in general.

Let us give an example of the impact of making a wrong assumption on the nature of the acting stochastic forcing. Were we to naively interpret one of the panels in Fig. 5 as resulting from the dynamics of a dynamical system perturbed by Gaussian noise, then we would have to conclude that the unperturbed deterministic system possesses at least two edge states on the basin boundary separating the competing basins of attraction (see Margazoglou et al., 2021, for a case where this situation applies). Hence, we would infer fundamentally wrong properties on the geometry of the phase space. Additionally, we would infer fundamentally different properties for the drift term.

Recent developments in data-driven methods based on the formalism of the Kramers–Moyal equation allow one to test accurately whether data are compatible with the hypothesis that stochasticity in the dynamics enters as a result of Gaussian noise or more general form of random forcing (Rydin Gorjão et al., 2021b; Li and Duan, 2022). Indeed, we point the reader to the recent contribution by Rydin Gorjão et al. (2021a), which shows that the analysis of proxy climatic datasets indicates the need to go beyond Langevin equation-based modelling, as they discover that it is necessary to treat noise as the sum of continuous and discontinuous processes. This indicates the need to consider, in future modelling exercises, the possibility of investigating the properties of metastable systems where the stochastic forcing comes as the result of simultaneous Gaussian and α-stable Lévy noise perturbations.

Appendix A: Stochastic perturbations of Lévy type

In this section, we provide a summary of the basic properties of a symmetric α-stable Lévy process in a Hilbert space in which the solutions to SPDE (Eq. 13) are defined. We also repeat some properties in the ℝⁿ space that are more familiar to a wide audience of readers. It is pertinent to refer to the distribution law of Lévy increments, its characteristic function, the Lévy–Itô decomposition, and the Lévy jump measure for a deeper study of the metastable behaviour of the stochastic climate system (Eq. 13). Let $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ be a given complete probability space and $H(\Vert \cdot \Vert ,\langle \cdot ,\cdot \rangle )$ a separable Hilbert space with the norm $\Vert \cdot \Vert$ and inner product $\langle \cdot ,\cdot \rangle$. A stochastic process (L^α(t)_t≥0) is a symmetric α-stable Lévy process in H if it satisfies the following:

1. L^α(0)=0, almost surely.

2. Independent increments – for any n∈ℕ and $\mathrm{0}\mathit{⩽}{t}_{\mathrm{1}}<t_{\mathrm{2}}<\mathrm{\cdots }<t_{n-\mathrm{1}}<t_n$, the vector is as follows:

   $$\begin{array}{}\text{(A1)}& \left(L^{\mathit{\alpha }}\right(t_{\mathrm{1}})-L^{\mathit{\alpha }}(t_{\mathrm{0}}),\mathrm{\dots },L^{\mathit{\alpha }}(t_n)-L^{\mathit{\alpha }}(t_{n-\mathrm{1}}\left)\right),\end{array}$$

   where there is a family of independent random vectors in H.

3. Stationary increments – for 0 ⩽ l<t random vectors L^α(t)−L^α(l) and L^α(t−l) have the same law 𝔏(.) in H, as follows:

   $$\begin{array}{}\text{(A2)}& \mathfrak{L}\left(L^{\mathit{\alpha }}\right(t)-L^{\mathit{\alpha }}(l\left)\right)=\mathfrak{L}\left(L^{\mathit{\alpha }}\right(t-l\left)\right).\end{array}$$

   This law in ℝⁿ is a symmetric α-stable distribution $\mathfrak{L}(.)=S_{\mathit{\alpha }}\left(\right(t-l{)}^{\frac{\mathrm{1}}{\mathit{\alpha }}},\mathrm{0},\mathrm{0})$, i.e. zero skewness and shift parameters, with a stability parameter $\mathit{\alpha }\in (\mathrm{0},\mathrm{2}]$ and a scaling parameter $(t-l{)}^{\frac{\mathrm{1}}{\mathit{\alpha }}}$. The stable distribution by the generalised central limit theorem (Schertzer and Lovejoy, 1997) is a limit in the distribution as n→∞ of the normalised sum $Y_n=\frac{\mathrm{1}}{B_n}\sum _{i=\mathrm{1}}^n(X_i-M_n)$ of n independent identically distributed random variables X_i, with a common probability distribution function F(x), which does not necessarily have to have moments of both the first and second order. A necessary and sufficient condition for this is as follows (Keller and Kuske, 2000; Burnecki et al., 2015):

   $$\begin{array}{}\text{(A3)}& \begin{array}{rl}F\left(x\right)& =[c_{\mathrm{1}}+r_{\mathrm{1}}(x\left)\right]|x{|}^{-\mathit{\alpha }},x<\mathrm{0},\\ & =\mathrm{1}-[c_{\mathrm{2}}+r_{\mathrm{2}}(x\left)\right]x^{-\mathit{\alpha }},x>\mathrm{0},\end{array}\end{array}$$

   with $\mathrm{0}<\mathit{\alpha }\le \mathrm{2}$, c₁, and c₂ positive constants, r₁(x)→0 as $x\to -\mathrm{\infty }$ and r₂(x)→0 as $x\to +\mathrm{\infty }$. When this condition holds and α=2, then we can set B_n=h(n), where h(n) satisfies h²=nln h, and the stable distribution is just the Gaussian law.

4. Sample paths are continuous in probability, i.e. for any t ⩾ 0 and η>0, as follows:

   $$\begin{array}{}\text{(A4)}& \underset{l\to t}{lim}\mathbb{P}(\Vert L^{\mathit{\alpha }}(t)-L^{\mathit{\alpha }}(l)\Vert >\mathit{\eta })=\mathrm{0}.\end{array}$$

   For $\mathit{\alpha }\in (\mathrm{0},\mathrm{2})$ the symmetric α-stable Lévy process in ℝⁿ has a characteristic function of the following form:

   $$\begin{array}{}\text{(A5)}& {\displaystyle }\mathbb{E}\left[e^{i\langle u,L^{\mathit{\alpha }}\left(t\right)\rangle }\right]=e^{-C\left(\mathit{\alpha }\right)t\left|\right|u|{|}^{\mathit{\alpha }}},u\in {\mathbb{R}}^n,t\ge \mathrm{0},\end{array}$$

   where $C\left(\mathit{\alpha }\right)={\mathit{\pi }}^{-\mathrm{1}/\mathrm{2}}\frac{\mathrm{\Gamma }\left(\right(\mathrm{1}+\mathit{\alpha })/\mathrm{2})\mathrm{\Gamma }(n/\mathrm{2})}{\mathrm{\Gamma }\left(\right(n+\mathit{\alpha })/\mathrm{2})}$, and Γ(.) is the Gamma function. In the case where α=2, we set $C\left(\mathrm{2}\right)=\mathrm{1}/\mathrm{2}$ and (A5) becomes the characteristic function of a standard Brownian motion. However, the Brownian motion cannot be seen as a weak limit of α-stable Lévy process because of the divergence C(α)→∞ as α→2. The properties of the sample paths of L^α(t) are, in fact, quite different for α=2 and α<2. First, the α-stable Lévy process is a discontinuous, pure jump process, while the Brownian motion has continuous paths. Second, the Brownian motion has moments of all orders, whereas $\mathbb{E}\left|L^{\mathit{\alpha }}\right(t){|}^{\mathit{\gamma }}<\mathrm{\infty }$ if γ<α. It can also be proved that the tails of L^α(t) are heavy, i.e. $\mathbb{P}\left(L^{\mathit{\alpha }}\right(t)>u)\sim u^{-\mathit{\alpha }},u\to \mathrm{\infty }$ , which is quite the opposite of the exponentially light Gaussian tails. Moreover, for $\mathit{\alpha }\in (\mathrm{0},\mathrm{1})$, the path variation in L^α(t) is bounded on finite time intervals and unbounded for $\mathit{\alpha }\in [\mathrm{1},\mathrm{2})$.

Although neither the incremental nor the marginal distributions of a Lévy process in general are representable by the elementary functions, the Lévy motion is completely determined by the Lévy–Khintchine formula, which specifies the characteristic function of the Lévy process.

If L^α(t) is a symmetric α-stable Lévy process in H, then, in the following:

1. The characteristic function of the Lévy–Khintchine formula is as follows:

   $${\mathrm{\Lambda }}_t\left(h\right)=\mathbb{E}\left[e^{i\langle h,L^{\mathit{\alpha }}\left(t\right)\rangle }\right]=e^{t\mathit{\psi }\left(h\right)},h\in H,t\ge \mathrm{0},$$

   where, in the following:

   $$\begin{array}{}\text{(A6)}& \mathit{\psi }\left(h\right)=\underset{H}{\int }\left(e^{i\langle h,y\rangle }-\mathrm{1}-i\langle h,y\rangle {\mathbf{1}}_{\mathit{\{}\mathrm{0}<\Vert y\Vert \mathit{⩽}\mathrm{1}\mathit{\}}}\right)\mathit{\nu }\left(\mathrm{d}y\right),\end{array}$$

   where 1_S is the indicator function for a set S, taking 1 on S and 0, otherwise, and ν is a Borel measure (also called the Lévy jump measure) in H, for which $\underset{H}{\int }(\mathrm{1}\wedge \Vert y{\Vert }^{\mathrm{2}})\mathit{\nu }\left(\mathrm{d}y\right)<\mathrm{\infty }$ with $\mathrm{1}\wedge \Vert y{\Vert }^{\mathrm{2}}=min\mathit{\{}\mathrm{1},\Vert y{\Vert }^{\mathrm{2}}\mathit{\}}$. A Borel measure, in addition, can be defined as the expected value of the number of jumps of specified size Q in the unit time interval, i.e. $\mathit{\nu }\left(Q\right)=\mathbb{E}N(\mathrm{1},Q)\left(\mathit{\omega }\right)$, ω∈Ω.

2. In the Lévy–Itô decomposition, for any sequence of positive radii r_n→0 and ${\mathcal{O}}_n=\mathit{\{}y\in H|r_{n+\mathrm{1}}<\left|\right|y\left|\right|\mathit{⩽}{r}_n\mathit{\}}$, there exists a sequence of independent compensated compound Poisson processes $\left({\overline{L}}_n\right(t){)}_{t\mathit{⩾}\mathrm{0}}$, n⩾0 in H, with jump measures ${\mathit{\nu }}_n\left(B\right)=\mathit{\nu }(B\cap {\mathcal{O}}_n)$ for B∈ℬ(H), the Borel σ-algebra in H, and n⩾1, which satisfy ℙ, almost surely for all t⩾0 as follows:

   $$\begin{array}{}\text{(A7)}& {\displaystyle }& {\displaystyle }L\left(t\right)=\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\overline{L}}_n+L_{\mathrm{0}}\left(t\right),\text{(A8)}& {\displaystyle }& {\displaystyle }{\overline{L}}_n\left(t\right)=L_n\left(t\right)-t\underset{H}{\int }y{\mathit{\nu }}_n\left(\mathrm{d}y\right),n\mathit{⩾}\mathrm{1}.\end{array}$$

   If L^α(t) is a symmetric α-stable Lévy process in ℝⁿ with the generating triplet $(\mathrm{0},\mathrm{0},{\mathit{\nu }}_{\mathit{\alpha }})$, then there exists an independent Poisson random measure N on ${\mathbb{R}}^{+}\times ({\mathbb{R}}^n\setminus \mathit{\{}\mathrm{0}\mathit{\left\}}\right)$ (quantifying the number of jumps of L^α(t)), such that, for each t⩾0, the following applies:

   $$\begin{array}{}\text{(A9)}& L^{\mathit{\alpha }}\left(t\right)=\underset{\parallel y\parallel <\mathrm{1}}{\int }y\stackrel{\mathrm{̃}}{N}(t,\mathrm{d}y)+\underset{\parallel y\parallel \mathit{⩾}\mathrm{1}}{\int }yN(t,\mathrm{d}y),\end{array}$$

   where $\stackrel{\mathrm{̃}}{N}(\mathrm{d}t,\mathrm{d}x)=N(\mathrm{d}t,\mathrm{d}x)-{\mathit{\nu }}_{\mathit{\alpha }}\left(\mathrm{d}x\right)\mathrm{d}t$ is the compensated Poisson random measure, and ν_α(dx) is the jump measure. The small $\parallel y\parallel <\mathrm{1}$ (large $\parallel y\parallel \mathit{⩾}\mathrm{1}$) jumps are controlled by $\stackrel{\mathrm{̃}}{N}(t,\mathrm{d}y)$ ( N(t,dy) ).

3. Its Lévy jump measure ν is symmetric in the sense that $\mathit{\nu }(-Q)=\mathit{\nu }\left(Q\right)$ for Q∈ℬ(H) and has the following specific geometry:

   $$\begin{array}{}\text{(A10)}& \mathit{\nu }\left(Q\right)=\underset{Q}{\int }\mathit{\nu }\left(\mathrm{d}y\right)=\underset{Q}{\int }{\displaystyle \frac{\mathrm{d}r}{r^{\mathrm{1}+\mathit{\alpha }}}}\mathit{\sigma }\left(\mathrm{d}s\right),\end{array}$$

   where $r=\Vert y\Vert$ and $s=y/\Vert y\Vert$ and $\mathit{\sigma }:\mathcal{B}(\partial B_{\mathrm{1}}(\mathrm{0}\left)\right)\to [\mathrm{0},\mathrm{\infty })$ is an arbitrary finite Radon measure on the unit sphere of H. The jump measure for a symmetric α-stable Lévy motion L^α(t) in ℝⁿ is defined by the following:

   $$\begin{array}{}\text{(A11)}& {\mathit{\nu }}_{\mathit{\alpha }}\left(\mathrm{d}u\right)=c(n,\mathit{\alpha }){\displaystyle \frac{\mathrm{d}u}{\left|\right|u|{|}^{n+\mathit{\alpha }}}},\end{array}$$

   with the intensity constant $c(n,\mathit{\alpha })=\frac{\mathit{\alpha }\mathrm{\Gamma }\left(\right(n+\mathit{\alpha })/\mathrm{2})}{{\mathrm{2}}^{\mathrm{1}-\mathit{\alpha }}{\mathit{\pi }}^{n/\mathrm{2}}\mathrm{\Gamma }(\mathrm{1}-\mathit{\alpha }/\mathrm{2})}$, where Γ(.) is the Gamma function (see Duan, 2015, and Applebaum, 2009).

One can come to a more intuitive interpretation of the stability parameter $\mathit{\alpha }\in (\mathrm{0},\mathrm{2})$ variation. For smaller values of α, the process is characterised by higher jumps with a lower frequency. As α increases, jumps decrease in height, and the frequency of their occurrence increases.

Appendix B: Probabilistic theory for the Lévy noise-induced escape

We briefly recapitulate here the main ideas behind the proof given in Debussche et al. (2013) of how the mean residence time in the competing metastable states of stochastically perturbed Chafee–Infante reaction–diffusion PDE scales with the intensity ε of the additive L(t) α-stable Lévy noise that acts as stochastic forcing.

One proceeds by considering the decomposition in the driving Lévy process by regularly varying the jump measure ν into small ξ^ε and large η^ε jump components. Let ${\mathrm{\Delta }}_{t}L=L\left(t\right)-L(t-)$ denote the jump increment of L at time t⩾0, and $\frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\rho }}}$ for ε, $\mathit{\rho }\in (\mathrm{0},\mathrm{1})$ the jump height threshold of L. The process η^ε is a compound Poisson process consisting of all jumps of height $\Vert {\mathrm{\Delta }}_{t}L\Vert >{\mathit{\epsilon }}^{-\mathit{\rho }}$, with the following intensity:

$$\begin{array}{}\text{(B1)}& {\mathit{\beta }}_{\mathit{\epsilon }}=\mathit{\nu }\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\rho }}}}{B}_{\mathrm{1}}^c\left(\mathrm{0}\right)\right)\approx {\mathit{\epsilon }}^{\mathit{\alpha }\mathit{\rho }},\end{array}$$

and the jump probability measure outside the ball $\frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\rho }}}{B}_{\mathrm{1}}\left(\mathrm{0}\right)$ is as follows:

$$\begin{array}{}\text{(B2)}& \mathit{\nu }\left(\cdot \cap {\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\rho }}}}{B}_{\mathrm{1}}^c\left(\mathrm{0}\right)\right)/{\mathit{\beta }}_{\mathit{\epsilon }},\end{array}$$

where B₁(0) is a ball of unit radius in H centred at the origin. The occurrence time of a kth large jump is defined recursively by the following:

$$\begin{array}{}\text{(B3)}& {\mathcal{Z}}_{\mathrm{0}}=\mathrm{0},{\mathcal{Z}}_k=inf\mathit{\{}t>{\mathcal{Z}}_{k-\mathrm{1}}|\Vert {\mathrm{\Delta }}_{t}L\Vert >{\mathit{\epsilon }}^{-\mathit{\rho }}\mathit{\}},k\mathit{⩾}\mathrm{1}.\end{array}$$

The waiting times between successive ${\mathit{\eta }}_t^{\mathit{\epsilon }}$ jumps have an exponential distribution ${\mathcal{Z}}_k-{\mathcal{Z}}_{k-\mathrm{1}}\sim \text{Exp}\left({\mathit{\beta }}_{\mathit{\epsilon }}\right)$.

Small jump processes ${\mathit{\xi }}^{\mathit{\epsilon }}=L-{\mathit{\eta }}^{\mathit{\epsilon }}$, due to the symmetry of Lévy measure ν, are a mean zero martingale in H with finite exponential moments. Probabilistic events causing small jumps in the stochastic solution of the system are not able to overcome the force of its deterministic stable state, and therefore, do not contribute to the exit from the basin of attraction. Formally, during the time between two large jumps $t_k={\mathcal{Z}}_k-{\mathcal{Z}}_{k-\mathrm{1}}$, the solution of Eq. (13), following the deterministic path (Eq. 1), returns to a small vicinity of the stable equilibria ${\mathit{\varphi }}^{\mathrm{W}/\text{SB}}$, as follows:

$$\begin{array}{}\text{(B4)}& \underset{x\in D^{\mathrm{W}/\text{SB}}}{sup}\underset{{\mathcal{Z}}_{k-\mathrm{1}}\mathit{⩽}t\mathit{⩽}{\mathcal{Z}}_k}{sup}\Vert \mathcal{T}\left(t\right)-T\left(t\right)\Vert \to \mathrm{0}\text{for}\mathit{\epsilon }\to \mathrm{0}.\end{array}$$

When a first large jump occurs, the solution process moves to the neighbouring domain of attraction with the following probability:

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}& \mathbb{P}\left({\mathit{\varphi }}^{\mathrm{W}/\text{SB}}+\mathit{\epsilon }{\mathrm{\Delta }}_{t_i}L\notin D^{\mathrm{W}/\text{SB}}\right)\\ & =\mathbb{P}\left({\mathrm{\Delta }}_{t_i}L\in {\displaystyle \frac{\mathrm{1}}{\mathit{\epsilon }}}\left[{\left(D^{\mathrm{W}/\text{SB}}\right)}^c-{\mathit{\varphi }}^{\mathrm{W}/\text{SB}}\right]\right)\\ & ={\displaystyle \frac{\mathit{\nu }\left(\frac{\mathrm{1}}{\mathit{\epsilon }}\left[(D^{\mathrm{W}/\text{SB}}{)}^c-{\mathit{\varphi }}^{\mathrm{W}/\text{SB}}\right]\cap \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\rho }}}{B}_{\mathrm{1}}^c\left(\mathrm{0}\right)\right)}{\mathit{\nu }\left(\frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\rho }}}{B}_{\mathrm{1}}^c\right(\mathrm{0}\left)\right)}}\approx {\mathit{\epsilon }}^{\mathit{\alpha }(\mathrm{1}-\mathit{\rho })}.\end{array}\end{array}$$

This is the probability that, at time t_i, there will be a jump increment ${\mathrm{\Delta }}_{t_i}L$ that exceeds the distance between the attractor and its domain of attraction boundary, as expressed by the jump probability measure (Eq. B2). In the zero-noise limit, the mean residence time in a basin of attraction is given by the following:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}& \mathbb{E}\left[\mathit{\tau }\right(\mathit{\epsilon }\left)\right]\\ & \approx \sum _{i=\mathrm{1}}^{\mathrm{\infty }}\mathbb{E}\left[{\mathcal{Z}}_i\right]\mathbb{P}\left[inf\left\{j:{\mathit{\varphi }}^{\mathrm{W}/\text{SB}}+\mathit{\epsilon }{\mathrm{\Delta }}_{t_j}L\notin D^{\mathrm{W}/\text{SB}}\right\}=i\right]\\ & \approx \mathbb{E}\left[t_{\mathrm{1}}\right]\mathbb{P}\left({\mathit{\varphi }}^{\mathrm{W}/\text{SB}}+\mathit{\epsilon }{\mathrm{\Delta }}_{t_{\mathrm{1}}}L\notin D^{\mathrm{W}/\text{SB}}\right)\\ & \cdot \sum _{i=\mathrm{1}}^{\mathrm{\infty }}i{\left(\mathrm{1}-\mathbb{P}\left[{\mathit{\varphi }}^{\mathrm{W}/\text{SB}}+\mathit{\epsilon }{\mathrm{\Delta }}_{t_{\mathrm{1}}}L\notin D^{\mathrm{W}/\text{SB}}\right]\right)}^{i-\mathrm{1}}\\ & \approx {\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\alpha }\mathit{\rho }}}}{\mathit{\epsilon }}^{\mathit{\alpha }(\mathrm{1}-\mathit{\rho })}{\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\alpha }(\mathrm{1}-\mathit{\rho })}}}\right)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\alpha }}}},\end{array}\end{array}$$

that is, by the sum of all the mean values of a large jump occurrence time times the probability that the minimum of a sample of size i of jump increments is sufficiently large to reach into the neighbouring domain of attraction. Thus, at the random time instant of large jumps, the solution process transitions, in an abrupt move, from one attractor to another. Such behaviour of the random dynamical system is known as a metastability.

In Debussche et al. (2013), it was proven that, in the timescale $\mathit{\lambda }\left(\mathit{\epsilon }\right)=\mathit{\nu }\left(\frac{\mathrm{1}}{\mathit{\epsilon }}{B}_{\mathrm{1}}^c\left(\mathrm{0}\right)\right),\mathit{\epsilon }>\mathrm{0}$ the metastable shifting of the diffusion process between neighbourhoods of the two attractors represents a continuous time Markov chain in a state space {ϕ^SB,ϕ^W}, with a transition rate matrix 𝔔, as follows:

$$\begin{array}{}\text{(B7)}& \mathfrak{Q}={\displaystyle \frac{\mathrm{1}}{\mathit{\mu }\left(B_{\mathrm{1}}^c\left(\mathrm{0}\right)\right)}}\left(\begin{array}{cc}-\mathit{\mu }\left(\right(D^{\text{SB}}-{\mathit{\varphi }}^{\text{SB}}{)}^c)& \mathit{\mu }\left(\right(D^{\text{SB}}-{\mathit{\varphi }}^{\text{SB}}{)}^c)\\ \mathit{\mu }\left(\right(D^{\mathrm{W}}-{\mathit{\varphi }}^{\mathrm{W}}{)}^c)& -\mathit{\mu }\left(\right(D^{\mathrm{W}}-{\mathit{\varphi }}^{\mathrm{W}}{)}^c)\end{array}\right),\end{array}$$

where μ(⋅) is the limit measure of ν.

Appendix C: Transitions induced by singular perturbations

In Sect. 4.2.3, and in particular Figs. 6 and 7, we have studied the effect of singular perturbations of a Lévy kick. The idea is that transitions in a system perturbed by Lévy noise are primarily driven by rare large jumps. By applying a singular perturbation of the form $\mathit{\mu }\to \mathit{\mu }+\mathit{\kappa }\mathit{\delta }\left(t\right)$ (where $\mathit{\mu }={\mathit{\mu }}_{\mathrm{0}}=\mathrm{1.05}$ throughout), we have been able to bracket the critical values ${\mathit{\kappa }}_{\text{crit}}^{\mathrm{W}\leftrightarrow \text{SB},\mathrm{s}}$ allowing for transitions between the two attractors. The expression κδ(t) is approximated as κΔ_τ(t), where ${\mathrm{\Delta }}_{\mathit{\tau }}\left(t\right)=\mathrm{1}/\mathit{\tau }$ if $\mathrm{0}<t<\mathit{\tau }$ and vanishes elsewhere. In Sect. 4.2.3 the results are shown for τ=1 year.

Table C1 Values of supercritical ${\mathit{\kappa }}_{\text{crit}}^{\mathrm{\dots },\mathrm{u}}$ for SB→W (second column) and W→SB (third column) for different approximations Δ_τ(t) of Dirac's δ(t).

Figure C1 Comparison between the supercritical paths constructed using singular perturbations of different duration τ, where, in the legend, m stands for month and y stands for years. (a) SB→W transitions. (b) W→SB transitions. Different values of τ correspond to different coloured lines. The basin boundaries are now depicted as thick black lines. Here, α=1.0 and ϵ=0.004.

We performed additional simulations to locate the supercritical and subcritical values of κ for τ = 1 month, 6 months, 2 years, and 4 years. The corresponding supercritical values of ${\mathit{\kappa }}_{\text{crit}}^{\text{SB}\to \mathrm{W},\mathrm{u}}$ and ${\mathit{\kappa }}_{\text{crit}}^{\mathrm{W}\to \text{SB},\mathrm{u}}$ are shown in Table C1. In Fig. C1, we plot the corresponding supercritical transition trajectories for the values of Table C1 at different durations. Notice that now we use coloured solid lines for the supercritical cases. To estimate the basin boundary, we record the final point of when the forcing was active, in the (𝒯,Δ𝒯) projected space, for each duration. This point for each case is particularly visible, in the insets of Fig. C1, as a rapid reflection of the trajectory, which then follows closely the basin boundary (depicted as a thick black line). The basin boundary we can explore through this procedure is then estimated by linking the points obtained when considering various values of τ. Notice that the estimated basin boundaries are slightly different when looking at the two SB→W and W→SB transitions, as the basin boundaries have folds than cannot be captured in the sampled 2-dimensional projection used in Fig. C1.

Finally, as stated earlier, from the third column of Table C1, we remark that, when considering forcings with durations of, for example, 2 years and longer, transitions from the W to the SB state can be achieved while retaining, at all times, a positive value for the solar irradiance because, while the forcing is active, its value is $\mathit{\mu }+\mathit{\kappa }/\mathit{\tau }$.

Appendix D: Estimates for the mean escape time

We report in Table D1 a summary of the statistics of the escape times from the W state and from the SB state for various choices of the noise law.

Table D1 Estimates and 95 % confidence intervals for the mean escape time τ from the W state and from the SB state, for Lévy noise with (a) α=0.5, (b) α=1.0, (c) α=1.5, and (d) Gaussian noise. N^o indicates the average number of transitions occurring in 10⁵ years of temporal evolution.