L\'evy-noise versus Gaussian-noise-induced Transitions in the Ghil-Sellers Energy Balance Model

We study the impact of applying stochastic forcing to the Ghil-Sellers energy balance climate model in the form of a fluctuating solar irradiance. Through numerical simulations, we explore the noise-induced transitions between the competing warm and snowball climate states. We consider multiplicative stochastic forcing driven by Gaussian and $\alpha$-stable L\'evy - $\alpha\in(0,2)$ - noise laws, examine the statistics of transition times, and estimate most probable transition paths. While the Gaussian noise case has been carefully studied in a plethora of investigations on metastable systems, much less is known about the L\'evy case, especially in the case of high- and infinite-dimensional systems. In the weak noise limit, the expected residence time in each metastable state scales in a fundamentally different way in the Gaussian vs. L\'evy noise case with respect to the intensity of the noise. In the former case, the classical Kramers-like exponential law is recovered. In the latter case, power laws are found, with the exponent equal to $-\alpha$, in apparent agreement with rigorous results obtained for additive noise in a related - yet different - reaction-diffusion equation as well as in simpler models. This can be better understood by treating the L\'evy noise as a compound Poisson process. The transition paths are studied in a projection of the state space and remarkable differences are observed between the two different types of noise. The snowball-to-warm and the warm-to-snowball most probable transition path cross at the single unstable edge state on the basin boundary. In the case of L\'evy 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.


Multistability of the Earth's Climate
The climate system comprises 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 throughtout 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, the present day Warm (W) state, and a competing one characterised by global glaciation, usually referred to as the Snowball (SB) state. Their analysis was performed using one-dimensional energy balance models (EBM)s, which, despite their simplicity, were able to capture the essential physical mechanism 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 is a negative, stabilizing one. Instead, the instability of the system is due to the presence of the so-called icealbedo feedback: an increase in the ice-covered fraction of the surface leads to further temperature reduction of the planet because ice reflects efficiently 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 snowball state is permitted, whereas above the critical value S SB→W , only the warm 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, stabilizing 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 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 these predictions were confirmed by actual data. Indeed, geological and paleomagnetic evidence suggests that during the Neoproterozoic era, between 630 and 715 million years ago, the Earth went at least twice into major long-lasting global glaciations that can be associated with the SB state; see Pierrehumbert et al. (2011) and Hoffman et al. (1998). Multicellular life emerged in our planet shortly after the final deglaciation from the last SB state (Gould, 1989). The robustness and importance of the competition between the Bolzmann 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 as well as 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 the properties of the competing states. The presence of multistability has a key importance in terms of determining habitability conditions for Earth-like exoplanets; see Lucarini et al. (2013) and Linsenmeier et al. (2015).
Additionally, several results indicate that the phase space of the climate system might well be more complex than the scenario of bistability described above. Various studies (Lewis et al., 2007;Abbot et al., 2011;Lucarini and Bódai, 2017;Margazoglou et al., 2021) performed with highly nontrivial climate models report the possible existence of additional competing states, up to a total of five (Brunetti et al., 2019;Ragon et al., 2021). 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.

Transitions between Competing Metastable States: Gaussian vs Lévy Noise
Clearly, in the case of autonomous systems where the phase space is partitioned in more than one basin of attraction of the corresponding attractors and the basin boundaries, the asymptotic state of the system is determined by its initial conditions. Things change dramatically when one includes time-dependent forcing which allows for transitions between competing metastable states (Ashwin et al., 2012). In particular, following the viewpoint originally proposed by Hasselmann (1976), whereby the fast variables of the climate system act as stochastic forcings for the slow ones (Imkeller and von Storch, 2001), the relevance of studying noise-induced transitions between competing states has become apparent (Hänggi, 1986;Freidlin and Wentzell, 1984). This viewpoint, where the noise is usually assumed to be Gaussian distributed, has provided very fruitful 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, named edge states. The edge states are saddles: 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 probability 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 Melancholia (M) state (Lucarini and Bódai, 2017). In previous papers, we have shown that it is possible to construct M states in high-dimensional climate models (Lucarini and Bódai, 2017) and to prove that the nonequilibrium quasi-potential formalism introduced by Graham (1987) and Graham et al. (1991) provides a powerful framework for explaining the population of each metastable state and the statistics of the noise-induced transitions. In the weak-noise limit, edge states act as gateways for noise-induced transitions between the metastable states (Lucarini and Bódai, 2019;Lucarini and Bódai, 2020;Margazoglou et al., 2021); see also a recent study on a nontrivial metastable prey-predator model (Garain and Sarathi Mandal, 2022). The local minima and the saddles of the quasi-potential Φ, which generalises the classical energy landscape for non-gradient systems, correspond to competing metastable states and to edge states, respectively.
In our investigation, the climate system is forced by adding a random -Gaussian distributed -component to the solar irradiance, which impacts, in the form of multiplicative noise, only a small subset of the degrees of freedom of the system. We remark that such choice of the stochastic forcing does not fully reflect physical realism, as the variability of the solar irradiance has a more complex behaviour (Solanki et al., 2013). Instead, noise acts as a tool for exploring the global stability properties of the system, and injecting noise as fluctuation of the solar irradiance has the merit of impacting the Lorenz energy cycle, thus effecting all degrees of freedom of the system (Lucarini and Bódai, 2020). See also the recent detailed mathematical analysis of the stochastically perturbed one-dimensional EBMs presented in Díaz and Díaz (2021).
A major limitation of this mathematical framework is the need to rigidly consider Gaussian noise laws, even if considerable freedom is left as to the choice of the spatial correlation properties of the noise. It seems natural to attempt a generalization by considering the whole class of α-stable Lévy noise laws. Lévy processes (Applebaum, 2009;Duan, 2015), described in detail below in Appendix A are fundamentally characterised by the stability parameter α ∈ (0, 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 α ∈ (0, 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 analyze and simulate at climate scales the ubiquitous intermittency and heavy-tailed statistics of clouds (Schertzer and Lovejoy, 1988), rain reflectivity (Tessier et al., 1993;Schertzer and Lovejoy, 1997), atmospheric turbulence (Schmitt et al., 1996), and soil moisture (Millàn et al., 2016). On longer time scales, 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 view point.
Mathematicians, on the other hand, have defined a Lévy multiplicative chaos (Fan, 1997;Rhodes et al., 2014) as a more mathematically tractable alternative to the universal multifractal. Finally, we remark that fractional Fokker-Planck equations have been proposed by Schertzer et al. (2001) to investigate the properties of nonlinear Langevin-type equations forced by a α− stable Lévy noise with the goal of analysing and simulating anomalous diffusion. Following Ditlevsen (1999), it has become apparent that more general classes of α-stable Lévy noise laws might be useful for modelling noise-induced transitions in the climate system like Dansgaard-Oeschger events. The viewpoint by Ditlevsen was particularly effective in stimulating mathematical investigations into noise-induced escapes from attractors where as stochastic forcing one chooses a Lévy, rather than Gaussian, noise (Imkeller and Pavlyukevich, 2006a, b;Chechkin et al., 2007;Debussche et al., 2013). Such analyses have clarified that a fundamental dichotomy exists with the classical Freidlin and Wentzell scenario mentioned above, even if phenomena like stochastic resonance can be recovered also in this case (Dybiec and Gudowska-Nowak, 2009;Kuhwald and Pavlyukevich, 2016). Whereas in the Gaussian case transitions between competing attractors occur as a result of the very unlikely combination of many steps all going in the right direction, in the Lévy case, transitions result from individual, very large and very rare jumps. Recently, Duan and collaborators have made fundamental progresses in achieving a variational formulation of Lévy noise-perturbed dynamical systems (Hu and Duan, 2020) as well as in developing corresponding methods for data assimilation (Gao et al., 2016) and data analysis . In terms of applications, Lévy noise is becoming a more and more a popular concept and tool for studying and interpreting complex systems (Grigoriu and Samorodnitsky, 2003;Penland and Sardeshmukh, 2012;Zheng et al., 2016;Wu et al., 2017;Serdukova et al., 2017;Cai et al., 2017;Singla and Parthasarathy, 2020;Gottwald, 2021).
The contribution by Gottwald (2021) is especially worth recapitulating because of its methodological clarity. There, the idea is, following Ditlevsen (1999), to provide a conceptual deterministic climate model able to generate a Lévy-noise-like signal to describe, at least qualitatively, abrupt climate changes similar to Dansgaard-Oeschger events, which are sequences of periods of abrupt warming followed by slower cooling that occurred during the last glacial period (Barker et al., 2011). A key building block is the idea proposed in Thompson et al. (2017) that a Lévy noise can be produced by integrating the so-called correlated additive and multiplicative (CAM) noise processes, which are defined 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 homogeneization 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 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; the sudden collapse of the West Antarctic ice sheet).

Outline of the Paper and Main Results
We consider here the Ghil-Sellers Earth's EBM Ghil (1976), a diffusive one-dimensional energy balance system, governed by a nonlinear reaction-diffusion parabolic partial differential equation. We stochastically perturb the system by adding random fluctuations to the solar irradiance, therefore the noise is introduced in 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 characterize noise-induced transitions between the competing climate basins and quantify the effect of noise parameters on them by estimating the statistics of escape times and empirically constructing mean transition pathways called instantons.
The results obtained confirm that, in the weak noise limit ε → 0, the mean residence time in each metastable state driven by Gaussian vs. Lévy noise has a fundamentally different dependence on ε. Indeed, as expected, in the Gaussian case the residence time grows exponentially with ε −2 , thus in basic agreement with the well-known Kramers (1940) lsaw and the previous studies performed on climate models (Lucarini and Bódai, 2019;Lucarini and Bódai, 2020). Instead, in the case of α-stable noise laws, the residence time increases with ε −α . We perform simulations for α = {0.5, 1.0, 1.5}. The obtained scaling can be explained by treating effectively the Lévy noise as compound Poisson process and is in agreement with what is found for low dimensional dynamics (Imkeller and Pavlyukevich, 2006a, b), as well as with the infinite dimensional stochastic Chafee-Infante reactiondiffusion equation (Debussche et al., 2013) in the case of additive noise. This might indicate that such scaling laws are more general than what typically considered.
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 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 organized as follows. In Section 2 we present the Ghil-Sellers EBM and summarize its most important dynamical aspects, as well as the steady-state solutions and their stability. The stochastic partial differential equation obtained by randomly perturbing the solar irradiance in the EBM is given in Subsection 3.1, where we also clarify the mathematical meaning of the solution of the stochastic partial differential equation. Subsection 3.2 introduces the mean residence time and most probable transition path between the competing climate states. The numerical methods are also briefly presented. In Section 4 we discuss our main results. In Section 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 behavior 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 energy balance climate model
The Ghil-Sellers EBM (Ghil, 1976) is described by an one-dimensional nonlinear, parabolic, reaction-diffusion partial differential equation (PDE) (1) involving the surface temperature T field and the transformed space variable where φ ∈ [−π/2, π/2] is the latitude. The model describes the processes of energy input, output, and diffusion across the domain and can be written as: 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 right hand side of (1) can be written as and describes the convergence of meridional heat transport performed by the geophysical fluids. The function K(x, T ) is a combined diffusion coefficient, given by 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 -of (1) describes the energy input associated with the absorption of incoming solar radiation can be written as where Q(x) is the incoming solar radiation and α a (x, T ) is the surface reflectivity (albedo), which is expressed as where the subscript {·} c denotes a cutoff for a generic quantity h defined as The term c 2 z(x) in (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 parametrization 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 -of (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 where σ is the Stefan-Boltzmann constant and the emissivity coefficient is expressed as 1 − m tanh(c 3 T 6 ). Such term describes, in a simple yet effective way, the greenhouse effect by reducing infrared radiation losses. The values of the empirical functions at discrete latitudes and empirical parameters c 1 , c 2 , c 3 , c 4 , c 5 , σ , 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 parametrizations 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 -co-exist, see 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, shown in Figure 1a. In Ghil (1976) the three stationary solutions were obtained by equating T t to 0, and it was shown, through linear stability analysis, that the stationary solutions T W and T SB are stable, while T M is unstable. In Bódai et al. (2015) the unstable solution T M was constructed using a modified version of the edge tracking algorithm (Skufca et al., 2006).
Following previous studies (Bódai et al., 2015;Lucarini and Bódai, 2019;Lucarini and Bódai, 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 T and the averaged Equator minus Poles temperature difference ∆T , defined as 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 T , while the large-scale energy transport performed by the geophysical fluids is controlled by ∆T . The boundary between high and low latitude in (11) is established at x = ±1/3, i.e. at 30 • N/S. Additionally, in some visualizations, we consider as a third coordinate the fraction of the surface with a below-freezing temperature (therefore we expect 1 for global glaciation and 0 for no ice). We refer to this variable as I, and it is an attempt to extract an observable that resembles the sea-ice percentage of the Earth's surface. Thus, the In order to analyze the influence of random perturbations on the deterministic dynamics of the climate model described in Section 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 boundary and initial conditions defined by Eq.
(2) apply to the stochastic temperature field T . 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 → SB direction. Instead, a strongly skewed process would have made it very hard to explore the full phase space, because 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 α ∈ (0, 2), so that we consider a jump process. We recall that the jumps become more frequent and less intense as α increases.
We First let us define the concept of mild solution in this context. Let (Ω, F, P) be a given complete probability space and H( · , ·, · ) a separable Hilbert space with norm · and inner product ·, · . Equation (13) can be rewritten in the more general form as follows T Under certain assumptions (Yagi, 2009), the problem (14) is formulated as a Cauchy problem whose local mild solution, a progressively measurable process T (t), for all t ∈ [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 (Green's function in physical terms) having the generalized semigroup property for the family of sector operators with the bounded inverses, and ϒ(T ) = T + F(x, T ), T ∈ H is a nonlinear operator, which we assume to be Lipschitz continuous. Following the abstract theory presented in Yagi (2009), under certain structural assumptions for the operators Ψ and ϒ and for the functional space, one can prove that the solution (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 defineL α=2 (t) =Ẇ (t), where (W (t) t≥0 ) is a Wiener process. We then defineẆ(t) = G(x, T )Ẇ (t) as the generalised derivative of a Wiener process in a suitably defined functional space.

Noise-induced Transitions: Mean Escape Times
By incorporating stochastic forcing into the system, its long-time dynamics change significantly, allowing transitions between the competing basins. This dynamical behaviour is called metastability, and is graphically captured by Figure 2, where in plots (a-b) a typical spatio-temporal evolution of the temperature field is shown, for stability parameters α = 0.5 and α = 1.5, respectively. Instead, in plots (c-d) the temporal noise-induced evolution of global temperature T and the averaged Equator and Poles temperature difference ∆T (as defined in Eqs. (10)-(11)) is shown for the same α. In what follows, we investigate the time statistics and the paths of the transitions between such basins.
In a complete probability space (Ω, F, P) we define the first exit time τ x of a cádlág mild solution T (·; x) of (13) starting at x ∈ D W /SB domain of warm/snowball climate stable state as The mean escape time is then expressed by E[τ x (ω)]. In the case of the infinite dimensional multistable reaction-diffusion system described by Chafee-Infante equation under the influence of additive infinite-dimensional α-stable Lévy noiseα ∈ (0, 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 to a finite state Markov chain with values in the set of stable states. Details of this method are given in Appendix B. Similar results have been obtained for bistable one-dimensional 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| −1−α , as discussed in Appendix A. Let's 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 have 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 −α crit . A similar argument applies when considering transitions triggered by negative fluctuations of the stochastic variable. Small-size events, which occur frequently and correspond to the non-occurrence of jumps, do not actually play any relevant role in determining the transitions, while they are responsible for the variability within each basin of attraction.
We now consider the case α = 2. While the corresponding finite dimensional problem is thoroughly documented in the literature (Freidlin and Wentzell, 1984), and has been applied in a similar context by some of the authors (Lucarini and Bódai, 2019;Lucarini and Bódai, 2020;Ghil and Lucarini, 2020;Margazoglou et al., 2021), the treatment of infinite dimensional SDEs driven by an infinite dimensional Wiener process via the Freidlin-Wentzell theory requires further extension. In the present context, we refer to Budhiraja and Dupuis (2000) and Budhiraja et al. (2008) and references therein where the problem of an infinite dimensional reaction-diffusion equation driven by an infinite dimensional Wiener process has been addressed.
We assume that steady state conditions and ergodicity are met, and we also assume that the analysing system is bistable and a unique edge state is present at the basin boundary, as in the case studied here. In the case of Gaussian noise, transitions between the competing basins of attraction are not determined by a single event as in the 0 < α < 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 has that the mean escape time from either basin of attraction decreases exponentially with increasing noise intensity ε and is given by a generalized Kramers' law where is the height of the quasi-potential barrier in the W attractor, and, correspondingly is the height of the quasi-potential barrier in the SB attractor, and Φ is the Graham's quasipotential mentioned above (Graham, 1987;Graham et al., 1991).
In our case, the theory indicates that if the stochastic forcing is Gaussian, under rather general hypothesis, the instanton will connect the attractor W/SB with the edge state M, which then acts as gateway for noise-induced transitions. Once the quasipotential barrier is overcome, a free fall "relaxation" trajectory links M with the competing attractor SB/W. For equilibrium systems, (e.g. for gradient flows) where detailed balance is achieved, the relaxation and instantonic trajectories within the same basin of attraction are identical. On the contrary, for non-equilibrium systems, the relaxation and instantonic trajectories will differ, and will only meet at the attractor. See a detailed discussion of this aspect and of the dynamical interpretation of the quasi-potential Φ in Lucarini and Bódai (2020) and Margazoglou et al. (2021). Instead, if the noise is of Lévy type, the theory formulated for simpler equations suggests that the instanton will connect the attractor with a region on the basin boundary that is the nearest, in the phase space, to the attractor, as the concept of quasi-potential is immaterial (Imkeller and Pavlyukevich, 2006a, b).
In general, the maximum likelihood transition trajectory T M (t) can be defined Lu and Duan, 2020) as a set of system states at each time moment t ∈ [0,t f ] that maximizes the conditional probability density function p( . | . ; . ) of passage from the origin stable state φ W /SB to the destination stable state φ SB/W and is expressed as where x 0 (x f ) belongs to the basin of attraction D W /SB (D SB/W ) and p ( . | . ) is the probability density function evolving according to the Fokker-Planck equation (Risken, 1996). This method is applicable either if efficient numerical algorithms are available to solve the Fokker-Planck equation associated to the studied stochastically driven system, or, empirically, when considering a large ensemble of simulations. Note that this is not an asymptotic approach, i.e. it does not require a weak noise limit ε → 0 for its application and is applicable for systems with either Gaussian or non-Gaussian noise. Yet, in the weak-noise limit, the definition (18) leads to constructing the optimal transition paths described above.
In the following section, for practical purposes, we construct such optimal transition path in the coarse grained 2D phase space (T , ∆T ) and 3D phase space (T , ∆T , I) of the variables defined in Sect. 2 by averaging the ensemble of transitions connecting the two competing states in the weak noise limit.

Numerical Methods
We solve Eq. (13) through the Matlab pdepe function, which is well suited for solving 1D parabolic and elliptic PDEs. We discretize the 1D space with a regular grid of 201 gridpoints, following Bódai et al. (2015).
The time span of integration t ∈ [0, T f ], varies for different cases, with T f ∈ (10 5 , 15 · 10 5 ) years, with time stepping of one 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 α (t 1 ), ..., L α (t N ) at each moment t j , j ∈ N, are obtained via 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 Figure 2  For the numerical simulations discussed below, we consider α = (0.5, 1.0, 1.5, 2) and ε ∈ (0.0001, 0.3). We select ε in such a way that the noise intensity is strong enough to induce at least 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-organized. Our simulations are performed 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

Results and Discussion
In what follows we aim at addressing three main questions: 1. What is 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,  (Lucarini et al., 2022) containing the raw data produced in this study as well as some illustrative animations portraying noise-induced transitions between the two competing metastable states.

Mean Escape Time
Our analysis confirms that there is fundamental dichotomy in the statistics of mean escape times between Lévy noise and Gaussian noise-induced transitions. Figure 3a shows the dependence of the mean escape time from either attractor on ε and α for the Lévy case. The red circles (blue squares) correspond to escapes from the W (SB) basin; see Lucarini et al. (2022) for additional details. The scaling ∝ ε −α 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 is 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 specifically when multiplicative Lévy noise is considered. The stochastically perturbed trajectories forced by Lévy noise consist of jumps, and the probability of occurrence of a high jump, which can trigger the escape from the reference basin of attraction, is polynomially small in noise intensity ε.
The Gaussian case -where no jumps are present -is portrayed in Fig. 3b. We show in semi-logarithmic scale the mean residence time versus 1/ε 2 . We perform a successful linear fit of the logarithm of the mean residence time in either attractor versus 1/ε 2 , and using Eq. (17), we obtain an estimate of the local quasi-potential barrier ∆Φ W/SB→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.

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

Gaussian noise
We estimate the transition paths by averaging among the escape plus relaxation trajectories using the run performed with the weakest noise, see Table D1. We first perform our analysis in the 2D-projected state space defined by (T , ∆T ). We prescribe two small circular-shaped regions enclosing the two deterministic attractors and search the timeseries of the portions of the whole trajectory that leave one of such 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 T the peak of histogram in the ∆T direction. The distributions are very peaked, and almost indistinguishable estimates for the instantonic and relaxation trajectories are obtained when computing the average of ∆T according to the 2D histogram conditional on the value of T .
In the background of Fig. 4 (a)   Let's provide some physical interpretation of how the transitions occur. Looking at the SB → W most probable path, the escape includes a simultaneous increase in T and ∆T . 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 µ, i.e. µ > µ. While the planet is warming globally, the Equator is warming faster than the Poles, resulting in a positive rate∆ T > 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 of the albedo in this latitude. Once the system crosses the M state, and supposing that persistent µ < µ 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 of sign in the rate of∆ T , and a subsequent decrease of ∆T 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. 4 (a). This implies the presence of less extreme non-equilibrium conditions compared to what observed in Margazoglou et al. (2021), where the W → SB and SB → W transitions occurred through fundamentally different paths; see discussion therein, especially regarding the role of the hydrological cycle. Figure 4 (b) presents the optimal transition paths W → SB and SB → W in a three-dimensional projection where we add as third coordinate the variable I, which indicates the fraction of the surface that has subfreezing temperatures (T < 273.15 K).
On the sides of the figure, two two dimensional projections on the (T , ∆I) and on the (∆T , I) planes are shown. Here, darker brown shadings indicate 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. 4 (a) could have been a simple effect of 2D projection. Instead, we see here that the SB → W and W → SB most probable transition paths also cross in the 3D projection in a well-defined region, which indeed corresponds to the M state (pink square).

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 systems 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 Figure 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 color scheme) in a small region around the deterministic attractors, as most typically the Lévy noise fluctuations of µ 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 a jump process. 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 characterized by the simultaneous decrease of both T and ∆T . 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 the supplementary material.

Lévy noise -Singular Perturbations
Based on what is discussed in Sec. 3.2, we expect that the transitions occur through the nearest region to the outgoing attractor in the basin boundary. We now try to clarify the properties of the most probable escape paths in the Lévy noise case by considering an additional set of simulations, taking inspiration from the edge tracking algorithm (Skufca et al., 2006). The idea is to exploit the fact that large jumps drive the transitions in the Lévy noise case. Starting from the deterministic SB state, we apply in Eq. 6 singular perturbations of the form µ → µ + κδ (t) and bracket the critical value κ SB→W The projections on the 2D phase space spanned by (T , ∆T ) of the supercritical and subcritical paths corresponding to the SB → W (W → SB) transition are shown in Fig. 6a (Fig. 6b)  referring to Lévy noise simulations performed using α = 1.0. Nonetheless, what 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 Lévy noise and the case of singularly perturbed trajectories what matters depends on the discontinuous nature of the Lévy noise.
In panel 6a (6b) the supercritical and subcritical paths are superimposed on the ensemble of 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, 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 ∆T 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. a) b) Figure 7. Same as panels a) and b) of Fig. 6, respectively, using the same 3D projection, lines, color coding as in Fig. 4. The subcritical and supercritical paths are depicted as thin and thick solid lines, respectively. The corresponding dashed lines shown in Fig. 6 are not reported here for simplicity. Here α = 1.0 and ε = 0.004.
Further support to the viewpoint proposed here is given by Fig. 7a and Fig. 7b, which are constructed along the lines of and vanishes elsewhere. The results are virtually unchanged if one considers τ < 1 y because the main dynamical processes of re-equilibration of the system act on longer time scales. The effect of the negative feedbacks of the system start 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 analyses in this direction dealing with stability of the large scale ocean circulation (Stocker and Schmittner, 1997;Lucarini et al., 2005Lucarini et al., , 2007Alkhayuon et al., 2019). Nonetheless, also in this case, the agreement with the results presented in Figs. 6-7 is considerable. We remark that considering longer lasting perturbations allows one to observe W → SB transitions without using at any time (unphysical) negative values for the solar irradiance. This is reassuring in terms of robustness and of physical sense of our results. Further details on the impact of the impact of considering different values of τ are reported in Appendix C.

Conclusions
It is a well-known that, as a result of the competition between the Boltzmann stabilizing feedback and the ice-albedo destabilizing feedback, under current astronomical and astrophysical conditions the climate system is multistable, as at least two competing and distinct climates are present, 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 of 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 modeling the observed climatic phenomena. In this regard, it is important to achieve a deeper understanding of the possible noise-induced transitions between competing stable climate states under α-stable Lévy perturbations and compare them with the Gaussian case.
As a starting point in this direction, we have studied the influence of different noise laws on the metastability properties of the randomly forced Ghil-Sellers EBM, which is governed by a nonlinear, parabolic, reaction-diffusion PDE. In the deterministic version of the model, we have three steady-state solutions: two stable, attractive climate states and one unstable saddle, corresponding to the edge state. The stable states correspond to the well-known W and SB climates. There is a fundamental dichotomy in the properties of the noise-induced transitions determined by whether we consider a stochastic forcing of intensity ε with a Gaussian versus an α-stable Lévy noise law. Note that, instead, the spatial structure of the noise is unchanged. This indicates that the phenomenology associated with the metastable behaviour depends critically on the choice of the noise law.
Not many studies have investigated, numerically or through mathematical theory, the properties of transitions in metastable systems driven by multiplicative Lévy noise, as done here.
First, in the weak noise limit ε → 0, the mean residence times inside either competing basin of attraction for diffusions driven by Gaussian vs. Lévy noise have a fundamentally different dependence on ε. Our results show that the logarithm of the mean residence time for Gaussian diffusions scales with ε −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.
Secondly, 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, yet qualitatively similar. In turn, Lévy diffusions leave the basin through the boundaries 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, 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 Lévy diffusion correspond to supercritical paths associated with Dirac's delta-like singular perturbations to the solar irradiance.
This viewpoint seems of general relevance in 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's 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 of Fig. 5 as resulting from the dynamics of a dynamical system perturbed by Gaussian noise, we would have to conclude that the unperturbed deterministic system possesses at least two edge states on the basin boundary separating the competing basins of attraction; see Margazoglou et al. (2021) for a case where this situation applies. Hence, we would infer fundamentally wrong properties on the geometry of the phase space. Additionally, we would infer fundamentally different properties for the drift term.
Recent developments in data-driven methods based on the formalism of the Kramers-Moyal equation allow to test accurately whether data are compatible with the hypothesis that stochasticity in the dynamics enters as a result of Gaussian noise or more general form of random forcing (Rydin Gorjão et al., 2021b;Li and Duan, 2022). Indeed, we point the reader to the recent contribution by Rydin Gorjão et al. (2021a), which shows that the analysis of proxy climatic datasets indicates the need to go beyond Langevin equation-based modelling, as they discover by 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.
Data availability. All the data used to produce the figures contained in this paper are publicly available as supplementary material in Lucarini et al. (2022) through the data repository https://doi.org/10.6084/m9.figshare.16802503.
Video supplement. Illustrative animations portraying noise-induced transitions can be found in the supplementary material. We further uploaded the relevant material in the YouTube platform with links: Lévy SB→W, Lévy W→SB, Gaussian SB→W, Gaussian W→SB.
Author contributions. All authors contributed equally to this work.
Competing interests. No competing interests are present 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). VL acknowledges the support provided by the EPSRC project EP/T018178/1. This paper is dedicated to K. Hasselmann.
Appendix A: Stochastic perturbations of Lévy type.
In this section we revise the basic properties of a symmetric α-stable Lévy process in a Hilbert space in which the solutions to SPDE (13) are defined. We also repeat some properties in R 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 behavior of the stochastic climate system (13). Let (Ω, F, P) be a given complete probability space and H( · , ·, · ) a separable Hilbert space with norm · and inner product ·, · . A stochastic process (L α (t) t≥0 ) is a symmetric α-stable Lévy process in H if it satisfies: 1) L α (0) = 0, a.s..
2) Independent increments: for any n ∈ N and 0 t 1 < t 2 < · · · < t n−1 < t n the vector is a family of independent random vectors in H.

3) Stationary increments
This law in R n is a symmetric α-stable distribution L(.) = S α ((t − l) 1 α , 0, 0), i.e., zero skewness and shift parameters, with a stability parameter α ∈ (0, 2] and a scaling parameter (t − l) 1 α . The stable distribution by the generalised central limit theorem (Schertzer and Lovejoy (1997)) is a limit in distribution as n → ∞ of the normalized sum Y n = 1 B n n ∑ i=1 (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 (Kuske and Keller (2000); Burnecki et al. (2015)) with 0 < α ≤ 2, c 1 and c 2 positive constants, r 1 (x) → 0 as x → −∞ and r 2 (x) → 0 as x → +∞. When this condition holds and α = 2 we can set B n = h(n), where h(n) satisfies h 2 = n ln h, and the stable distribution is just the Gaussian law. 4) Sample paths are continuous in probability, i.e. for any t 0 and η > 0 For α ∈ (0, 2) the symmetric α-stable Lévy process in R n has characteristic function of the form where C(α) = π −1/2 Γ((1+α)/2)Γ(n/2) Γ((n+α)/2) , and Γ(.) is the Gamma function. In the case where α = 2 we set C(2) = 1/2 and (A5) becomes the characteristic function of a standard Brownian motion. However, 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. Firstly, the α-stable Lévy process is a discontinuous, pure jump process, while the Brownian motion has continuous paths. Secondly, the Brownian motion has moments of all orders, whereas E |L α (t)| γ < ∞ iff γ < α. It can also be proved that the tails of L α (t) are heavy, i.e. P (L α (t) > u) ∼ u −α , u → ∞ , quite the opposite of the exponentially light Gaussian tails. Moreover for α ∈ (0, 1), the path variation of L α (t) is bounded on finite time intervals, and unbounded for α ∈ [1, 2).
Although neither the incremental nor the marginal distribution 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: 1) (Lévy-Khintchine formula) Its characteristic function is here 1 S is the indicator function for a set S, taking 1 on S and 0 otherwise, and ν is a Borel measure (so called the Lévy jump measure) in H for which H (1 ∧ y 2 )ν(dy) < ∞ with 1 ∧ y 2 = min{1, y 2 }. A Borel measure, as well, can be defined as the expected value of the number of jumps of specified size Q in the unit time interval, i.e. ν(Q) = EN(1, Q)(ω), ω ∈ Ω.
2) (Lévy-Itô decomposition) For any sequence of positive radii r n → 0 and O n = {y ∈ H | r n+1 < ||y|| r n } there exist a sequence of independent compensated compound Poisson processes (L n (t)) t 0 , n 0 in H with jump measures ν n (B) = ν(B ∩ O n ) for B ∈ B(H) the Borel σ -algebra in H and n 1, which satisfy P-almost surely for all t 0 L n (t) = L n (t) − t H yν n (dy), n 1.
3) Its Lévy jump measure ν is symmetric in the sense that ν(−Q) = ν(Q) for Q ∈ B(H) and has the specific geometry where r = y and s = y/ y and σ : B(∂ B 1 (0)) → [0, ∞) 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 R n is defined by with the intensity constant c(n, α) = αΓ((n+α)/2) 2 1−α π n/2 Γ(1−α/2) where Γ(.) is the Gamma function; see Duan (2015); Applebaum (2009). One can come to a more intuitive interpretation of the stability parameter α ∈ (0, 2) variation: for smaller values of α, the process is characterized by higher jumps with a lower frequency. As α increases, jumps decrease in height and the frequency of their occurrence increases.
Appendix B: Probabilistic theory for the Lévy noise-induced escape.
We briefly recapitulate here the main ideas behind the proof given in Debussche et al. (2013) of how the mean residence time in the competing metastable states of stochastically perturbed Chafee-Infante reaction-diffusion PDE scales with the intensity ε of the additive L(t) α−stable Lévy noise that acts as stochastic forcing.
One proceeds by considering the decomposition of the driving Lévy process with regularly varying the jump measure ν into small ξ ε and large η ε jump components. Let ∆ t L = L(t) − L(t−) denote the jump increment of L at time t 0, and 1 ε ρ for ε, ρ ∈ (0, 1) the jump height threshold of L. The process η ε is a compound Poisson process consisting of all jumps of height ∆ t L > ε −ρ with intensity and the jump probability measure outside the ball 1 ε ρ B 1 (0) by where B 1 (0) is a ball of unit radius in H centered at the origin. The occurrence time of a k-th large jump is defined recursively by The waiting times between sucessive η ε t jumps have an exponential distribution Z k − Z k−1 ∼ Exp(β ε ).
Small jump processes ξ ε = L − η ε due to the symmetry of Lévy measure ν is 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 = Z k − Z k−1 , the solution of (13) following the deterministic path (1) returns to a small vicinity of the stable equilibria φ W /SB When a first large jump occurs, the solution process moves to the neighboring domain of attraction with probability This is the probability that at time t i there will be a jump increment ∆ t i L that exceeds the distance between the attractor and its domain of attraction boundary, expressed by the jump probability measure (B2). In the zero-noise limit the mean residence time in a basin of attraction is given by i.e. by the sum of all the mean values of large jump occurrence time times the probability that the minimum of a sample of size i of jump increments is sufficiently large to get into the neighboring 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 behavior of the random dynamical system is known as a metastability.
In Debussche et al. (2013) it was proved that in the time scale λ (ε) = ν( 1 ε B c 1 (0)), ε > 0 the metastable shifting of the diffusion process between neighborhoods of the two attractors represents a continuous time Markov chain in state space {φ SB , φ W } with a transition rate matrix Q where µ(·) is the limit measure of ν.

Appendix C: Transitions induced by singular perturbations
In Sect. 4.2.3, and in particular Figs. 6-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 µ → µ + κδ (t) (where µ = µ 0 = 1.05 throughout) we have been able to bracket the critical values κ W ↔SB,s crit allowing for transitions between the two attractors. The expression κδ (t) is approximated as κ∆ τ (t) where ∆ τ (t) = 1/τ if 0 < t < τ and vanishes elsewhere. In Sect. 4.2.3 the results are shown for τ = 1y.
We performed additional simulations to locate the supercritical and subcritical values of κ for τ = 1 month (m), 6 m, 2 y, and 4 y. The corresponding supercritical values of κ SB→W,u crit and κ W →SB,u crit are shown in Table C1. In Fig. C1 we plot the corresponding supercritical transition trajectories for the values of Tab. C1 at different duration. Notice that now we use colored solid lines for the supercritical cases. To estimate the basin boundary we record the final point of when the forcing was active, in the (T , ∆T ) 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 has folds than cannot be captured in the simpled 2-dimensional projection used in Fig. C1.
Finally, as stated earlier, from the third column of Table C1, we remark that, when considering forcings with duration of e.g.
2 y 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 µ + κ/τ.

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