the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Lévy noise versus Gaussiannoiseinduced transitions in the Ghil–Sellers energy balance model
Valerio Lucarini
Larissa Serdukova
Georgios Margazoglou
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 noiseinduced 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 infinitedimensional 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 Kramerslike 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 snowballtowarm and the warmtosnowball 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.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 presentday warm (W) state and a competing one characterised by global glaciation, usually referred to as the snowball (SB) state. Their analysis was performed using onedimensional 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 socalled ice–albedo feedback, whereby an increase in the icecovered 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^{6} and 715×10^{6} years ago, the Earth went, at least twice, into major longlasting 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 Earthlike 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 timedependent 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 noiseinduced 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 highdimensional climate models (Lucarini and Bódai, 2017) and to prove that the nonequilibrium quasipotential 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 noiseinduced transitions. In the weak noise limit, edge states act as gateways for noiseinduced 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 quasipotential Φ, which generalises the classical energy landscape for nongradient systems, correspond to competing metastable states and to edge states, respectively. In our investigation, the climate system is forced by adding a random – Gaussiandistributed – 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 onedimensional 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 heavytailed 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 Langevintype 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 noiseinduced 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 noiseinduced 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 GudowskaNowak, 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 noiseperturbed 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évynoiselike 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 socalled 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, onedimensional 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 noiseinduced 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 ε^{−2}, which is, thus, in basic agreement with the wellknown 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 lowdimensional 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 steadystate 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.
The Ghil–Sellers EBM (Ghil, 1976) is described by a onedimensional 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:
where C(x) is the effective heat capacity, and $T=T(x,t)$ has boundary and initial conditions, as follows:
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 righthand side of Eq. (1) can be written as follows:
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:
The empirical functions k_{1}(x) and k_{2}(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 righthand side of Eq. (1) describes the energy input associated with the absorption of incoming solar radiation and can be written as follows:
where 𝒬(x) is the incoming solar radiation, and α_{a}(x,T) is the surface reflectivity (albedo), which is expressed as follows:
where the subscript {⋅}_{c} denotes a cutoff for a generic quantity h, defined as follows:
The term c_{2}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 righthand 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:
where σ is the Stefan–Boltzmann constant, and the emissivity coefficient is expressed as 1−mtanh (c_{3}T^{6}). 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).
In this study, we consider μ=1.05. For this value of μ, two stable asymptotic states – the W and the SB states – coexist (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 edgetracking 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{\u203e}}{T}$ and the averaged Equator minus the poles' temperature difference ΔT, which is defined as follows:
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{\u203e}}{T}$, while the largescale 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 belowfreezing 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{\u203e}}{T}$, correspond to ΔT_{W} = 16 K, ΔT_{SB} = 8.3 K, ΔT_{M} = 17.5 K, ${\stackrel{\mathrm{\u203e}}{T}}_{\mathrm{W}}$ = 297.7 K, ${\stackrel{\mathrm{\u203e}}{T}}_{\text{SB}}$ = 235.1 K, ${\stackrel{\mathrm{\u203e}}{T}}_{\mathrm{M}}$ = 258 K, I_{W}=0.2, I_{SB}=1, and I_{M}=1.
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):
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:
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_{0}∈H has the following integral representation:
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 Noiseinduced transitions: mean escape times
By incorporating stochastic forcing into the system, its longtime 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{\u203e}}{\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:
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 infinitedimensional α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 onedimensional 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. Smallsize events, which occur frequently and correspond to the nonoccurrence 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 steadystate 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:
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 quasipotential 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 quasipotential barrier in the SB attractor, and Φ is the Graham quasipotential mentioned above (Graham, 1987; Graham et al., 1991).
3.3 Noiseinduced 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 VandenEijnden, 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 noiseinduced transitions. Once the quasipotential barrier is overcome, a freefall 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 nonequilibrium 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 quasipotential Φ 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 quasipotential 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(\phantom{\rule{0.33em}{0ex}}.\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}.\phantom{\rule{0.33em}{0ex}};\phantom{\rule{0.33em}{0ex}}.\phantom{\rule{0.33em}{0ex}})$ 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:
where x_{0} (x_{f}) belongs to the basin of attraction ${D}^{\mathrm{W}/\text{SB}}$ (${D}^{\text{SB}/\mathrm{W}}$), and $p\phantom{\rule{0.33em}{0ex}}(\phantom{\rule{0.33em}{0ex}}.\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}.\phantom{\rule{0.33em}{0ex}})$ 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 nonGaussian 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 coarsegrained 2D phase space ($\stackrel{\mathrm{\u203e}}{\mathcal{T}},\mathrm{\Delta}\mathcal{T}$) and 3D phase space ($\stackrel{\mathrm{\u203e}}{\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.
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:
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 longerlasting (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).
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 noiseinduced 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 noiseinduced 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 semilogarithmic 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 quasipotential 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.
4.2 Escape paths for the noiseinduced 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 2Dprojected state space defined by $(\stackrel{\mathrm{\u203e}}{\mathcal{T}},\mathrm{\Delta}\mathcal{T})$. We prescribe two small circularshaped 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{\u203e}}{\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{\u203e}}{\mathcal{T}}$.
In the background of Fig. 4a, we show the empirical estimate of the invariant measure in the 2Dprojected state space defined by $(\stackrel{\mathrm{\u203e}}{\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 nonvanishing probability currents and the breakdown of detailed balance, is a signature of nonequilibrium 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{\u203e}}{\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{\u0303}}{\mathit{\mu}}$, i.e. $\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\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 nonequilibrium 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 threedimensional 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 twodimensional 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 welldefined 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.
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{\u0303}}{\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{\u203e}}{\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).
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 edgetracking 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 reequilibration 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 largescale 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 longerlasting 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.
It is a wellknown 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 noiseinduced 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 steadystate solutions, i.e. two stable, attractive climate states and one unstable saddle, corresponding to the edge state. The stable states correspond to the wellknown W and SB climates. There is a fundamental dichotomy in the properties of the noiseinduced 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 ε^{−2}, 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 datadriven 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 equationbased 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.
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 ℝ^{n} 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:

L^{α}(0)=0, almost surely.

Independent increments – for any n∈ℕ and $\mathrm{0}\mathit{\u2a7d}{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.

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(tl\left)\right).\end{array}$$This law in ℝ^{n} is a symmetric αstable distribution $\mathfrak{L}(.)={S}_{\mathit{\alpha}}\left(\right(tl{)}^{\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 $(tl{)}^{\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]\phantom{\rule{0.25em}{0ex}}x{}^{\mathit{\alpha}},x<\mathrm{0},\\ & =\mathrm{1}[{c}_{\mathrm{2}}+{r}_{\mathrm{2}}(x\left)\right]\phantom{\rule{0.25em}{0ex}}{x}^{\mathit{\alpha}},x>\mathrm{0},\end{array}\end{array}$$with $\mathrm{0}<\mathit{\alpha}\le \mathrm{2}$, c_{1}, and c_{2} positive constants, r_{1}(x)→0 as $x\to \mathrm{\infty}$ and r_{2}(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^{2}=nln h, and the stable distribution is just the Gaussian law.

Sample paths are continuous in probability, i.e. for any t ⩾ 0 and η>0, as follows:
$$\begin{array}{}\text{(A4)}& \underset{{\scriptscriptstyle 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 ℝ^{n} 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)\phantom{\rule{0.25em}{0ex}}t\phantom{\rule{0.25em}{0ex}}\left\rightu{}^{\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}}\phantom{\rule{0.25em}{0ex}}\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}\phantom{\rule{0.25em}{0ex}}\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}\phantom{\rule{0.25em}{0ex}}\left({L}^{\mathit{\alpha}}\right(t)>u)\sim {u}^{\mathit{\alpha}},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}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:

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{{\scriptscriptstyle 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{\u2a7d}\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{{\scriptscriptstyle 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)$, ω∈Ω.

In the Lévy–Itô decomposition, for any sequence of positive radii r_{n}→0 and ${\mathcal{O}}_{n}=\mathit{\{}y\in H\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}{r}_{n+\mathrm{1}}<\left\righty\left\right\mathit{\u2a7d}{r}_{n}\mathit{\}}$, there exists a sequence of independent compensated compound Poisson processes $\left({\overline{L}}_{n}\right(t){)}_{t\mathit{\u2a7e}\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{\u2a7e}\mathrm{1}.\end{array}$$If L^{α}(t) is a symmetric αstable Lévy process in ℝ^{n} 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{\u0303}}{N}(t,\mathrm{d}y)+\underset{\parallel y\parallel \mathit{\u2a7e}\mathrm{1}}{\int}yN(t,\mathrm{d}y),\end{array}$$where $\stackrel{\mathrm{\u0303}}{N}(\mathrm{d}t,\mathrm{d}x)=N(\mathrm{d}t,\mathrm{d}x){\mathit{\nu}}_{\mathit{\alpha}}\left(\mathrm{d}x\right)\phantom{\rule{0.25em}{0ex}}\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{\u2a7e}\mathrm{1}$) jumps are controlled by $\stackrel{\mathrm{\u0303}}{N}(t,\mathrm{d}y)$ ( N(t,dy) ).

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{{\scriptscriptstyle Q}}{\int}\mathit{\nu}\left(\mathrm{d}y\right)=\underset{{\scriptscriptstyle 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 ℝ^{n} 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\rightu{}^{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.
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:
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:
where B_{1}(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:
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:
When a first large jump occurs, the solution process moves to the neighbouring domain of attraction with the following probability:
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 zeronoise limit, the mean residence time in a basin of attraction is given by the following:
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),\phantom{\rule{0.25em}{0ex}}\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:
where μ(⋅) is the limit measure of ν.
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.
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 2dimensional 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}$.
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.
All the data used to produce the figures contained in this paper are publicly available in the supplement (Lucarini et al., 2022) through the data repository https://doi.org/10.6084/m9.figshare.16802503.
Illustrative animations portraying noiseinduced transitions can be found in Lucarini et al. (2022, https://doi.org/10.6084/m9.figshare.16802503). We, furthermore, uploaded the relevant material to YouTube at the following: https://youtu.be/2l77YZ9VKlo (last access: 5 May 2022), which is the Lévy SB→W; https://youtu.be/aReg0KYevU (last access: 5 May 2022), which is the Lévy W→SB; https://youtu.be/qrL9A2QNYZs (last access: 5 May 2022), which is the Gaussian SB→W; and https://youtu.be/zIzNX9gkCTo (last access: 5 May 2022), which is the Gaussian W→SB.
VL conceptualized the paper, developed the methodology, validated the results, conducted the data analysis. and took the lead role in writing and revising the paper. LS and GM conceptualized the paper, developed the methodology and software, and helped with the writing and revising of the paper.
At least one of the (co)authors is a member of the editorial board of Nonlinear Processes in Geophysics. The peerreview process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This article is part of the special issue “Centennial issue on nonlinear geophysics: accomplishments of the past, challenges of the future”. It is not associated with a conference.
The authors wish to thank Peter Ashwin, Reyk Börner, Jesus I. Diaz, Jinqiao Duan, Michael, Ghil, Tobias Grafke, Serafim Kalliadasis, Alessandro Laio, XueMei Li, and Greg Pavliotis, for the useful exchanges on various topics covered in this paper. Valerio Lucarini wishes to thank Myles Allen, for suggesting that we look into 3D projections of the phase space when studying transitions paths. The authors acknowledge the support provided by the EU Horizon 2020 project TiPES (grant no. 820970). Valerio Lucarini acknowledges the support provided by the EPSRC project (grant no. EP/T018178/1). This paper is dedicated to K. Hasselmann.
This research has been supported by the Horizon 2020 (TiPES; grant no. 820970) and the Engineering and Physical Sciences Research Council (grant no. EP/T018178/1).
This paper was edited by Daniel Schertzer and reviewed by two anonymous referees.
Abbot, D. S., Voigt, A., and Koll, D.: The Jormungand global climate state and implications for Neoproterozoic glaciations, J. Geophys. Res.Atmos., 116, D18103, https://doi.org/10.1029/2011JD015927, 2011. a
Alharbi, R.: Nonlinear parabolic stochastic partial differential equation with application to finance, Doctoral thesis (PhD), University of Sussex, Brighton, http://sro.sussex.ac.uk/id/eprint/96730 (last access:5 May 2022), 2021. a, b
Alkhayuon, H., Ashwin, P., Jackson, L. C., Quinn, C., and Wood, R. A.: Basin bifurcations, oscillatory instability and rateinduced thresholds for Atlantic meridional overturning circulation in a global oceanic box model, P. Roy. Soc. AMath. Phy., 475, 20190051, https://doi.org/10.1098/rspa.2019.0051, 2019. a
Applebaum, D.: Lévy Processes and Stochastic Calculus, Cambridge Studies in Advanced Mathematics, 2nd edn., Cambridge University Press, https://doi.org/10.1017/CBO9780511809781, 2009. a, b
Ashwin, P., Wieczorek, S., Vitolo, R., and Cox, P.: Tipping points in open systems: bifurcation, noiseinduced and ratedependent examples in the climate system, Philos. T. R. Soc. A, 370, 1166–1184, https://doi.org/10.1098/rsta.2011.0306, 2012. a
Barker, S., Knorr, G., Edwards, R. L., Parrenin, F., Putnam, A. E., Skinner, L. C., Wolff, E., and Ziegler, M.: 800,000 Years of Abrupt Climate Variability, Science, 334, 347–351, https://doi.org/10.1126/science.1203580, 2011. a
Bensid, S. and Díaz, J. I.: On the exact number of monotone solutions of a simplified Budyko climate model and their different stability, Discrete Cont. Dyn.B, 24, 1033–1047, 2019. a, b
Benzi, R., Sutera, A., and Vulpiani, A.: The mechanism of stochastic resonance, J. Phys. AMath. Gen., 14, L453–L457, https://doi.org/10.1088/03054470/14/11/006, 1981. a
Bódai, T., Lucarini, V., Lunkeit, F., and Boschi, R.: Global instability in the Ghil–Sellers model, Clim. Dynam., 44, 3361–3381, 2015. a, b, c, d
Bouchet, F., Gawedzki, K., and Nardini, C.: Perturbative calculation of quasipotential in nonequilibrium diffusions: a meanfield example, J. Stat. Phys., 163, 1157–1210, https://doi.org/10.1007/s1095501615032, 2016. a
Brunetti, M., Kasparian, J., and Vérard, C.: Coexisting climate attractors in a coupled aquaplanet, Clim. Dynam., 53, 6293–6308, https://doi.org/10.1007/s00382019049267, 2019. a
Budhiraja, A. and Dupuis, P.: A variational representation for positive functionals of infinite dimensional Brownian motion, Probability and Mathematical Statistics–Wroclaw University, 20, 39–61, 2000. a
Budhiraja, A., Dupuis, P., and Maroulas, V.: Large deviations for infinite dimensional stochastic dynamical systems, Ann. Probab., 36, 1390–1420, https://doi.org/10.1214/07AOP362, 2008. a
Budyko, M. I.: The effect of solar radiation variations on the climate of the Earth, Tellus, 21, 611–619, https://doi.org/10.3402/tellusa.v21i5.10109, 1969. a, b
Burnecki, K., Wylomanska, A., and Chechkin, A.: Discriminating between Light And heavytailed distributions with limit theorem, PLoS ONE, 10, e0145604, https://doi.org/10.1371/journal.pone.0145604, 2015. a
Burrage, K. and Lythe, G.: Accurate stationary densities with partitioned numerical methods for stochastic partial differential equations, Stochastics and Partial Differential Equations: Analysis and Computations, 2, 262–280, https://doi.org/10.1007/s4007201400328, 2014. a
Cai, R., Chen, X., Duan, J., Kurths, J., and Li, X.: Lévy noiseinduced escape in an excitable system, J. Stat. Mech. Theory E., 2017, 063503, https://doi.org/10.1088/17425468/aa727c, 2017. a
Chechkin, A., Sliusarenko, O., Metzler, R., and Klafter, J.: Barrier crossing driven by Levy noise: Universality and the Role of Noise Intensity, Phys. Rev. E, 75, 041101, https://doi.org/10.1103/PhysRevE.75.041101, 2007. a
Cialenco, I., Fasshauer, G. E., and Ye, Q.: Approximation of stochastic partial differential equations by a kernelbased collocation method, Int. J. Comput. Math., 89, 2543–2561, https://doi.org/10.1080/00207160.2012.688111, 2012. a
Dai, M., Gao, T., Lu, Y., Zheng, Y., and Duan, J.: Detecting the maximum likelihood transition path from data of stochastic dynamical systems, Chaos, 30, 113124, https://doi.org/10.1063/5.0012858, 2020. a, b
Davie, A. M. and Gaines, J. G.: Convergence of numerical schemes for the solution of parabolic stochastic partial differential equations, Math. Comput., 70, 121–134, https://doi.org/10.1090/s0025571800012242, 2000. a
Debussche, A., Högele, M., and Imkeller, P.: The dynamics of nonlinear reactiondiffusion equations with small lévy noise, in: Lecture Notes in Mathematics, Springer, Berlin, https://doi.org/10.1007/9783319008288_1, 2013. a, b, c, d, e, f, g
Díaz, G. and Díaz, J. I.: Stochastic energy balance climate models with Legendre weighted diffusion and an additive cylindrical Wiener process forcing, Discrete Cont. Dyn.S., https://doi.org/10.3934/dcdss.2021165, 2021. a
Díaz, J. I., Hernández, J., and Tello, L.: On the Multiplicity of Equilibrium Solutions to a Nonlinear Diffusion Equation on a Manifold Arising in Climatology, J. Math. Anal. Appl., 216, 593–613, https://doi.org/10.1006/jmaa.1997.5691, 1997. a, b
Ditlevsen, P. D.: Observation of αstable noise induced millennial climate changes from an icecore record, Geophys. Res. Lett., 26, 1441–1444, https://doi.org/10.1029/1999GL900252, 1999. a, b, c
Doering, C. R.: A stochastic partial differential equation with multiplicative noise, Phys. Lett. A, 122, 133–139, https://doi.org/10.1016/03759601(87)907912, 1987. a, b
Duan, J.: An introduction to stochastic dynamics, Cambridge University Press, New York, 2015. a, b, c
Duan, J. and Wang, W.: Effective Dynamics of Stochastic Partial Differential Equations, Elsevier, Boston, https://doi.org/10.1016/C2013015235X, 2014. a, b
Dybiec, B. and GudowskaNowak, E.: Lévy stable noiseinduced transitions: stochastic resonance, resonant activation and dynamic hysteresis, J. Stat. Mech. Theory E., 2009, P05004, https://doi.org/10.1088/17425468/2009/05/p05004, 2009. a
Fan, A. H.: Sur les chaos de Lévy stables d'indice $\mathrm{0}<\mathit{\alpha}<\mathrm{1}$, Ann. Sci. Math. Québec, 1, 53–66, 1997. a
Feudel, U., Pisarchik, A. N., and Showalter, K.: Multistability and tipping: From mathematics and physics to climate and brain–Minireview and preface to the focus issue, Chaos, 28, 033501, https://doi.org/10.1063/1.5027718, 2018. a
Freidlin, M. I. and Wentzell, A. D.: Random perturbations of dynamical systems, Springer, New York, 1984. a, b
Gao, T., Duan, J., Kan, X., and Cheng, Z.: Dynamical inference for transitions in stochastic systems with αstable Lévy noise, J. Phys. AMath. Theor., 49, 294002, https://doi.org/10.1088/17518113/49/29/294002, 2016. a
Garain, K. and Sarathi Mandal, P.: Stochastic sensitivity analysis and early warning signals of critical transitions in a tristable prey–predator system with noise, Chaos, 32, 033115, https://doi.org/10.1063/5.0074242, 2022. a
Ghil, M.: Climate stability for a Sellerstype model, J. Atmos. Sci., 33, 3–20, https://doi.org/10.1175/15200469(1976)033<0003:CSFAST>2.0.CO;2, 1976. a, b, c, d, e, f
Ghil, M.: EnergyBalance Models: An Introduction, in: Climatic Variations and Variability: Facts and Theories: NATO Advanced Study Institute First Course of the International School of Climatology, Ettore Majorana Center for Scientific Culture, Erice, Italy, March 9–21, 1980, edited by: Berger, A., Springer Netherlands, Dordrecht, 461–481, https://doi.org/10.1007/9789400985148_27, 1981. a, b
Ghil, M.: A mathematical theory of climate sensitivity or, How to deal with both anthropogenic forcing and natural variability?, in: Climate Change: Multidecadal and Beyond, edited by Chang, P. C., Ghil, M., Latif, M., and Wallace, J. M., World Scientific/Imperial College Press, 31–51, 2015. a
Ghil, M. and Childress, S.: Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics, SpringerVerlag, Berlin, 1987. a
Ghil, M. and Lucarini, V.: The physics of climate variability and climate change, Rev. Mod. Phys., 92, 035002, https://doi.org/10.1103/RevModPhys.92.035002, 2020. a, b
Gottwald, G.: A model for DansgaardOeschger events and millennialscale abrupt climate change without external forcing, Clim. Dynam., 56, 227–243, https://doi.org/10.1007/s0038202005476z, 2021. a, b
Gottwald, G. A. and Melbourne, I.: Homogenization for deterministic maps and multiplicative noise, P. Roy. Soc. AMath. Phy., 469, 20130201, https://doi.org/10.1098/rspa.2013.0201, 2013. a
Gould, S. J.: Wonderful Life: The Burgess shale and the Nature of History, W.W. Norton, New York, 1989. a
Grafke, T. and VandenEijnden, E.: Numerical computation of rare events via large deviation theory, Chaos, 29, 063118, https://doi.org/10.1063/1.5084025, 2019. a
Grafke, T., Grauer, R., and Schäfer, T.: The instanton method and its numerical implementation in fluid mechanics, J. Phys. AMath. Theor., 48 333001, https://doi.org/10.1088/17518113/48/33/333001, 2015. a, b
Grafke, T., Schäfer, T., and VandenEijnden, E.: Long Term Effects of Small Random Perturbations on Dynamical Systems: Theoretical and Computational Tools, in: Recent Progress and Modern Challenges in Applied Mathematics, Modeling and Computational Science, edited by: Melnik, R., Makarov, R., and Belair, J., Fields Institute Communications, Springer, New York, NY, https://doi.org/10.1007/9781493969692_2, pp. 17–55, 2017. a
Graham, R.: Macroscopic potentials, bifurcations and noise in dissipative systems, in: Fluctuations and Stochastic Phenomena in Condensed Matter, edited by: Garrido, L., Springer Berlin Heidelberg, 1–34, ISBN',9783540474012, 1987. a, b
Graham, R., Hamm, A., and Tél, T.: Nonequilibrium potentials for dynamical systems with fractal attractors or repellers, Phys. Rev. Lett., 66, 3089–3092, https://doi.org/10.1103/PhysRevLett.66.3089, 1991. a, b
Grebogi, C., Ott, E., and Yorke, J. A.: Fractal Basin Boundaries, LongLived Chaotic Transients, and UnstableUnstable Pair Bifurcation, Phys. Rev. Lett., 50, 935–938, https://doi.org/10.1103/PhysRevLett.50.935, 1983. a
Grigoriu, M. and Samorodnitsky, G.: Dynamic Systems Driven by Poisson/Lévy White Noise, in: IUTAM Symposium on Nonlinear Stochastic Dynamics, edited by: Namachchivaya, N. S. and Lin, Y. K., Springer Netherlands, Dordrecht, 319–330, https://doi.org/10.1007/9789401001793_28, 2003. a
Hänggi, P.: Escape from a metastable state, J. Stat. Phys., 42, 105–148, 1986. a
Hasselmann, K.: Stochastic climate models, Part I. Theory, Tellus, 28, 473–485, 1976. a
Hetzer, G.: The structure of the principal component for semilinear diffusion equations from energy balance climate models, Houston J. Math., 16, 203–216, 1990. a, b
Hoffman, P. F., Kaufman, A. J., Halverson, G. P., and Schrag, D. P.: A Neoproterozoic Snowball Earth, Science, 281, 1342–1346, https://doi.org/10.1126/science.281.5381.1342, 1998. a
Hu, J. and Duan, J.: OnsagerMachlup action functional for stochastic partial differential equations with Levy noise, arXiv [preprint], https://doi.org/10.48550/ARXIV.2011.09690, 4 December 2020. a, b, c
Imkeller, P. and Pavlyukevich, I.: First exit times of SDEs driven by stable L'evy processes, Stoch. Proc. Appl., 116, 611–642, https://doi.org/10.1016/j.spa.2005.11.006, 2006a. a, b, c, d
Imkeller, P. and Pavlyukevich, I.: Lévy flights: transitions and metastability, J. Phys. AMath. Gen., 39, L237–L246, https://doi.org/10.1088/03054470/39/15/l01, 2006b. a, b, c, d
Imkeller, P. and von Storch, J. S.: Stochastic Climate Models, Birkhauser, Basel, 2001. a
Jentzen, A. and Kloeden, P. E.: The numerical approximation of stochastic partial differential equations, Milan J. Math., 77, 205–244, https://doi.org/10.1007/s0003200901000, 2009. a
Kaper, H. and Engler, H.: Mathematics and climate, SIAM, Philadelphia, 2013. a, b
Keller, J. and Kuske, R.: Rate of convergence to a stable law, SIAM J. Appl. Math., 61, 1308–1323, https://doi.org/10.1137/s0036139998342715, 2000. a
Kloeden, P. E. and Shott, S.: Linearimplicit strong schemes for itôgalkerin approximations of stochastic PDES, Journal of Applied Mathematics and Stochastic Analysis, 14, 697341, https://doi.org/10.1155/S1048953301000053, 2001. a
Kramers, H.: Brownian motion in a field of force and the diffusion model of chemical reactions, Physica, 7, 284–304, https://doi.org/10.1016/S00318914(40)900982, 1940. a
Kraut, S. and Feudel, U.: Multistability, noise, and attractor hopping: The crucial role of chaotic saddles, Phys. Rev. E, 66, 015207, https://doi.org/10.1103/PhysRevE.66.015207, 2002. a
Kuhwald, I. and Pavlyukevich, I.: Stochastic Resonance in Systems Driven by αStable Lévy Noise, International Conference on Vibration Problems 2015, Procedia Engineer., 144, 1307–1314, https://doi.org/10.1016/j.proeng.2016.05.129, 2016. a
Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping elements in the Earth's climate system, P. Nat. Acad. Sci. USA, 105, 1786–1793, 2008. a
Lewis, J. P., Weaver, A. J., and Eby, M.: Snowball versus slushball Earth: Dynamic versus nondynamic sea ice?, J. Geophys. Res., 112, C11014, https://doi.org/10.1029/2006JC004037, 2007. a
Li, Y. and Duan, J.: Extracting Governing Laws from Sample Path Data of NonGaussian Stochastic Dynamical Systems, J. Stat. Phys., 186, 30, https://doi.org/10.1007/s1095502202873y, 2022. a
Linsenmeier, M., Pascale, S., and Lucarini, V.: Climate of Earthlike planets with high obliquity and eccentric orbits: Implications for habitability conditions, Planet. Space Sci., 105, 43–59, https://doi.org/10.1016/j.pss.2014.11.003, 2015. a
Lovejoy, S. and Schertzer, D.: The Weather and Climate: Emergent Laws and Multifractal Cascades, Cambridge University Press, Cambridge, 2013. a
Lu, Y. and Duan, J.: Discovering transition phenomena from data of stochastic dynamical systems with Lévy noise, Chaos, 30, 093110, https://doi.org/10.1063/5.0004450, 2020. a, b, c
Lucarini, V. and Bódai, T.: Edge states in the climate system: exploring global instabilities and critical transitions, Nonlinearity, 30, R32–R66, https://doi.org/10.1088/13616544/aa6b11, 2017. a, b, c
Lucarini, V. and Bódai, T.: Transitions across Melancholia States in a Climate Model: Reconciling the Deterministic and Stochastic Points of View, Phys. Rev. Lett., 122, 158701, https://doi.org/10.1103/PhysRevLett.122.158701, 2019. a, b, c, d, e
Lucarini, V. and Bódai, T.: Global stability properties of the climate: Melancholia states, invariant measures, and phase transitions, Nonlinearity, 33, R59–R92, https://doi.org/10.1088/13616544/ab86cc, 2020. a, b, c, d, e, f, g
Lucarini, V., Calmanti, S., and Artale, V.: Destabilization of the thermohaline circulation by transient changes in the hydrological cycle, Clim. Dynam., 24, 253–262, https://doi.org/10.1007/s003820040484z, 2005. a
Lucarini, V., Calmanti, S., and Artale, V.: Experimental mathematics: Dependence of the stability properties of a twodimensional model of the Atlantic ocean circulation on the boundary conditions, Russ. J. Math. Phys., 14, 224–231, https://doi.org/10.1134/S1061920807020124, 2007. a
Lucarini, V., Fraedrich, K., and Lunkeit, F.: Thermodynamic analysis of snowball Earth hysteresis experiment: Efficiency, entropy production and irreversibility, Q. J. Roy. Meteor. Soc., 136, 2–11, https://doi.org/10.1002/qj.543, 2010. a
Lucarini, V., Pascale, S., Boschi, R., Kirk, E., and Iro, N.: Habitability and Multistability in Earthlike Planets, Astron. Nachr., 334, 576–588, https://doi.org/10.1002/asna.201311903, 2013. a
Lucarini, V., Blender, R., Herbert, C., Ragone, F., Pascale, S., and Wouters, J.: Mathematical and physical ideas for climate science, Rev. Geophys., 52, 809–859, https://doi.org/10.1002/2013RG000446, 2014a. a
Lucarini, V., Serdukova, L., and Margazoglou, G.: Lévynoise versus Gaussiannoiseinduced Transitions in the GhilSellers Energy Balance Model, figshare, https://doi.org/10.6084/m9.figshare.16802503, 2022. a, b, c, d, e, f
Margazoglou, G., Grafke, T., Laio, A., and Lucarini, V.: Dynamical landscape and multistability of a climate model, P. Roy. Soc. AMath. Phy., 477, 20210019, https://doi.org/10.1098/rspa.2021.0019, 2021. a, b, c, d, e, f, g, h, i
Millàn, H., Cumbrera, R., and Tarquis, A. M.: Multifractal and Levystable statistics of soil surface moisture distribution derived from 2D image analysis, Appl. Math. Model., 40, 2384–2395, https://doi.org/10.1016/j.apm.2015.09.063, 2016. a
Nicolis, C.: Stochastic aspects of climatic transitions – response to a periodic forcing, Tellus, 34, 308–308, https://doi.org/10.3402/tellusa.v34i3.10817, 1982. a
North, G. and Stevens, M.: Energybalance climate models, in: Frontiers of Climate Modeling, edited by: Kiehl, J. T. and Ramanathan, V., Cambridge University Press, Cambridge, 52–72, https://doi.org/10.1017/CBO9780511535857.004, 2006. a, b
North, G. R.: Multiple solutions in energy balance climate models, Global Planet. Change, 2, 225–235, https://doi.org/10.1016/09218181(90)90003U, 1990. a, b
North, G. R., Cahalan, R. F., and Coakley Jr., J. A.: Energy balance climate models, Rev. Geophys., 19, 91–121, https://doi.org/10.1029/RG019i001p00091, 1981. a, b
Ott, E.: Chaos in dynamical systems, 2nd edn., Cambridge University Press, Cambridge, https://doi.org/10.1017/CBO9780511803260, 2002. a
Pavliotis, G. and Stuart, A.: Multiscale methods, Texts in applied mathematics, Springer, New York, NY, 2008. a, b
Peixoto, J. P. and Oort, A. H.: Physics of Climate, AIP Press, New York, New York, 1992. a
Penland, C. and Sardeshmukh, P. D.: Alternative interpretations of powerlaw distributions found in nature, Chaos, 22, 023119, https://doi.org/10.1063/1.4706504, 2012. a
Peszat, S. and Zabczyk, J.: Stochastic Partial Differential Equations with Levy Noise: An Evolution Equation Approach, Cambridge University Press, https://doi.org/10.1017/cbo9780511721373, 2007. a, b
Pierrehumbert, R., Abbot, D., Voigt, A., and Koll, D.: Climate of the Neoproterozoic, Annu. Rev. Earth Pl. Sc., 39, 417–460, https://doi.org/10.1146/annurevearth040809152447, 2011. a, b
Ragon, C., Lembo, V., Lucarini, V., Vérard, C., Kasparian, J., and Brunetti, M.: Robustness of Competing Climatic States, J. Climate, 35, 2769–2784, https://doi.org/10.1175/JCLID210148.1, 2022. a
Rhodes, R., Sohier, J., and Vargas, V.: Levy multiplicative chaos and star scale invariant random measures, Ann. Probab., 42, 689–724, https://doi.org/10.1214/12AOP810, 2014. a
Risken, H.: The Fokker–Planck equation, Springer, Berlin, 1996. a
Rydin Gorjão, L., Riechers, K., Hassanibesheli, F., Witthaut, D., Lind, P. G., and Boers, N.: Changes in stability and jumps in Dansgaard–Oeschger events: a data analysis aided by the Kramers–Moyal equation, Earth Syst. Dynam. Discuss. [preprint], https://doi.org/10.5194/esd202195, in review, 2021a. a
Rydin Gorjão, L., Witthaut, D., Lehnertz, K., and Lind, P. G.: ArbitraryOrder FiniteTime Corrections for the KramersMoyal Operator, Entropy, 23, 517, https://doi.org/10.3390/e23050517, 2021b. a
Saltzman, B.: Dynamical Paleoclimatology: Generalized Theory of Global Climate Change, Academic Press New York, New York, 2001. a
Schertzer, D. and Lovejoy, S.: Multifractal simulations and analysis of clouds by multiplicative processes, Atmos. Res., 21, 337–361, https://doi.org/10.1016/01698095(88)90035X, 1988. a
Schertzer, D. and Lovejoy, S.: Universal multifractals do exist!: Comments on “a statistical analysis of mesoscale rainfall as a random Cascade”, J. Appl. Meteorol., 36, 1296–1303, https://doi.org/10.1175/15200450(1997)036<1296:UMDECO>2.0.CO;2, 1997. a, b
Schertzer, D., Larchevêque, M., Duan, J., Yanovsky, V. V., and Lovejoy, S.: Fractional Fokker–Planck equation for nonlinear stochastic differential equations driven by nonGaussian Lévy stable noises, J. Math. Phys., 42, 200, https://doi.org/10.1063/1.1318734, 2001. a
Schmitt, F., Lovejoy, S., and Schertzer, D.: Multifractal analysis of the Greenland IceCore Project climate data, Geophys. Res. Lett., 22, 1689–1692, https://doi.org/10.1029/95GL01522, 1995. a
Schmitt, F., Schertzer, D., Lovejoy, S., and Brunet, Y.: Multifractal temperature and flux of temperature variance in fully developed turbulence, Europhys. Lett., 34, 195–200, https://doi.org/10.1209/epl/i1996004384, 1996. a
Sellers, W. D.: A global climatic model based on the energy balance of the earthatmosphere system, J. Appl. Meteorol., 8, 392–400, 1969. a, b
Serdukova, L., Zheng, Y., Duan, J., and Kurths, J.: Metastability for discontinuous dynamical systems under Lévy noise: Case study on Amazonian Vegetation, Sci. Rep.UK, 7, 9336, https://doi.org/10.1038/s41598017076868, 2017. a
Singla, R. and Parthasarathy, H.: Quantum robots perturbed by Levy processes: Stochastic analysis and simulations, Commun. Nonlinear Sci., 83, 105142, https://doi.org/10.1016/j.cnsns.2019.105142, 2020. a
Skufca, J. D., Yorke, J. A., and Eckhardt, B.: Edge of Chaos in a Parallel Shear Flow, Phys. Rev. Lett., 96, 174101, https://doi.org/10.1103/PhysRevLett.96.174101, 2006. a, b, c
Solanki, S. K., Krivova, N. A., and Haigh, J. D.: Solar Irradiance Variability and Climate, Annu. Rev. Astron. Astr., 51, 311–351, https://doi.org/10.1146/annurevastro082812141007, 2013. a
Steffen, W., Rockström, J., Richardson, K., Lenton, T. M., Folke, C., Liverman, D., Summerhayes, C. P., Barnosky, A. D., Cornell, S. E., Crucifix, M., Donges, J. F., Fetzer, I., Lade, S. J., Scheffer, M., Winkelmann, R., and Schellnhuber, H. J.: Trajectories of the Earth System in the Anthropocene, P. Nat. Acad. Sci. USA, 115, 8252–8259, https://doi.org/10.1073/pnas.1810141115, 2018. a
Stocker, T. F. and Schmittner, A.: Influence of CO_{2} emission rates on the stability of the thermohaline circulation, Nature, 388, 862–865, https://doi.org/10.1038/42224, 1997. a
Tessier, Y., Lovejoy, S., and Schertzer, D.: Universal multifractals: theory and observations for rain and clouds, J. Appl. Meteorol. Clim., 32, 223–250, https://doi.org/10.1175/15200450(1993)032<0223:UMTAOF>2.0.CO;2, 1993. a
Thompson, W. F., Kuske, R. A., and Monahan, A. H.: Reduced αstable dynamics for multiple time scale systems forced with correlated additive and multiplicative Gaussian white noise, Chaos, 27, 113105, https://doi.org/10.1063/1.4985675, 2017. a
Varadhan, S. R. S.: Large deviations and applications, Society for Industrial and Applied Mathematics Philadelphia, 75 pp., https://doi.org/10.2307/2287939, 1985. a
Voigt, A. and Marotzke, J.: The transition from the presentday climate to a modern Snowball Earth, Clim. Dynam., 35, 887–905, https://doi.org/10.1007/s0038200906335, 2010. a
Vollmer, J., Schneider, T. M., and Eckhardt, B.: Basin boundary, edge of chaos and edge state in a twodimensional model, New J. Phys., 11, 013040, https://doi.org/10.1088/13672630/11/1/013040, 2009. a
Weron, A. and Weron, R.: Computer simulation of Levy alphastable variables and processes, Chaos – The Interplay Between Stochastic and Deterministic Behaviour, edited by: Garbaczewski, P., Wolf, M., and Weron, A., Springer Berlin Heidelberg, Berlin, Heidelberg, 379–392, ISBN 9783540447221, 1995. a
Wu, J., Xu, Y., and Ma, J.: Lévy noise improves the electrical activity in a neuron under electromagnetic radiation, PloS one, 12, e0174330–e0174330, https://doi.org/10.1371/journal.pone.0174330, 2017. a
Yagi, A.: Dynamical Systems, in: Abstract Parabolic Evolution Equations and their Applications, Springer Monographs in Mathematics, Springer Berlin Heidelberg, https://doi.org/10.1007/9783642046315, 2010. a, b
Zheng, Y., Serdukova, L., Duan, J., and Kurths, J.: Transitions in a genetic transcriptional regulatory system under Lévy motion, Sci. Rep.UK, 6, 29274, https://doi.org/10.1038/srep29274, 2016. a
Zheng, Y., Yang, F., Duan, J., Sun, X., Fu, L., and Kurths, J.: The maximum likelihood climate change for global warming under the influence of greenhouse effect and Lévy noise, Chaos, 30, 013132, https://doi.org/10.1063/1.5129003, 2020. a, b, c
 Abstract
 Introduction
 The Ghil–Sellers energy balance climate model
 Background and methodology
 Results and discussion
 Conclusions
 Appendix A: Stochastic perturbations of Lévy type
 Appendix B: Probabilistic theory for the Lévy noiseinduced escape
 Appendix C: Transitions induced by singular perturbations
 Appendix D: Estimates for the mean escape time
 Data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 The Ghil–Sellers energy balance climate model
 Background and methodology
 Results and discussion
 Conclusions
 Appendix A: Stochastic perturbations of Lévy type
 Appendix B: Probabilistic theory for the Lévy noiseinduced escape
 Appendix C: Transitions induced by singular perturbations
 Appendix D: Estimates for the mean escape time
 Data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References