Global Stability Properties of the Climate: Melancholia States, Invariant Measures, and Phase Transitions
GGlobal Stability Properties of the Climate: Melancholia States, InvariantMeasures, and Phase Transitions
Valerio Lucarini , , and Tam´as B´odai , , , Department of Mathematics and Statistics, University of Reading, Reading UK Centre for the Mathematics of Planet Earth, University of Reading, Reading UK CEN, University of Hamburg, Hamburg, Germany Pusan National University, Busan, Republic of Korea Center for Climate Physics, Institute for Basic Science, Busan, Republic of Korea (Dated: March 12, 2020) a r X i v : . [ phy s i c s . a o - ph ] M a r bstract For a wide range of values of the incoming solar radiation, the Earth features at least two attracting states,which correspond to competing climates. The warm climate is analogous to the present one; the snowballclimate features global glaciation and conditions that can hardly support life forms. Paleoclimatic evidencessuggest that in the past our planet flipped between these two states. The main physical mechanism respon-sible for such an instability is the ice-albedo feedback. In a previous work, we defined the Melancholiastates that sit between the two climates. Such states are embedded in the boundaries between the two basinsof attraction and feature extensive glaciation down to relatively low latitudes. Here, we explore the globalstability properties of the system by introducing random perturbations as modulations to the intensity of theincoming solar radiation. We observe noise-induced transitions between the competing basins of attraction.In the weak noise limit, large deviation laws define the invariant measure, the statistics of escape times,and typical escape paths called instantons. By constructing the instantons empirically, we show that theMelancholia states are the gateways for the noise-induced transitions. In the region of multistability, in thezero-noise limit, the measure is supported only on one of the competing attractors. For low (high) valuesof the solar irradiance, the limit measure is the snowball (warm) climate. The changeover between the tworegimes corresponds to a first-order phase transition in the system. The framework we propose seems ofgeneral relevance for the study of complex multistable systems. Finally, we put forward a new method forconstructing Melancholia states from direct numerical simulations, which provides a possible alternativewith respect to the edge-tracking algorithm.
MSC Classification: 76U05 Rotating fluids; 82C26 Dynamic and nonequilibrium phase transitions (gen-eral); 60Gxx Stochastic processes; 37D45 Strange attractors, chaotic dynamics; 85A20 Planetary atmo-spheres; 76E20 Stability and instability of geophysical and astrophysical flows . INTRODUCTION In the late 1960s and in the 1970s, Budyko, Selllers, and Ghil [1–3] proposed the idea thatthe Earth, in the current astrophysical and astronomical configuration, supports two co-existingattractors, the warm (W) state we live in, and the so-called snowball (SB) state, which is charac-terised by global glaciation and globally averaged surface temperature of about 200-220 K. Usingparsimonious yet physically meaningful energy balance models, they indicated that the bistabilityof the climate system is the result of the competition between the positive ice-albedo feedback(a glaciated surface reflects the incoming radiation more effectively) and the negative Boltzmannradiative feedback (a warmer surface emits more radiation to space). The relevance of the theo-retical ansatz became apparent when paleoclimatic data showed that, indeed, our planet has beenflipping in and out of states of global glaciations corresponding to the predicated SB states duringthe Proterozoic, about 650 Mya [4–6]. W According to these energy balance models, the Earth’sclimate is bistable for a substantial range of values of the solar irradiance S ∗ , which include thepresent day value. Below the critical value S ∗ W → SB , only the SB state is permitted, whereas abovethe critical value S ∗ SB → W , only the W state is permitted. Such critical values, which determinethe boundaries of the region in parametric space where bistability is realised, are defined by bi-furcations that occur when, roughly speaking, the strength of the positive, destabilising feedbacksbecomes as strong as the negative, stabilizing feedbacks. Indeed, models of different levels ofcomplexity ranging up to the state-of-the-art Earth System Models currently used for climate pro-jections agree on predicting the existence of multistability in the climate system and point to thefundamental mechanisms described above as responsible for it, as well as providing values for S ∗ W → SB that are in broad agreement with those obtained using simple models [7–9]. We remarkthat both the concentration of greenhouse gases and the position of the continents have an impacton the values of S ∗ W → SB and S ∗ SB → W and on the properties of the W and SB states [10]. Extremelyhigh values of the concentration of CO seem to be needed to deglaciate from a SB state [11].Improving our understanding of the critical transitions associated to such a bistability is a keychallenge of geosciences and has strong implications also in terms of the quest for understand-ing or anticipating planetary habitability. Planets in the habitable zone have astronomical andastrophysical configurations that allow, in principle, the presence of water at surface. Therefore, S ∗ W → SB defines the cold boundary of the habitable zone. Clearly, an exoplanet in the habitablezone can be in the regime of bistability: if the planet is in the SB state, it will have very hard3ime supporting life . Additionally, astronomical parameters such as the obliquity of the planet[12–14], eccentricity [12], or the length of the year [15, 16] can have a dramatic effect on theproperties of multistability of a planet in the habitable zone, up to the point of erasing it altogether.In particular, planets with Earth-like atmospheres and high seasonal variability can have ice-freeareas at much larger distance from the host star than planets without seasonal variability, whichleads to a substantial expansion of the outer edge of the habitable zone. Additionally, one expectsthat tidally locked planets with an active carbon cycle can never be found in a SB state [17]. A. Multistability of the Climate System
Further investigations have proposed the possibility of the existence of alternative cold stateswith respect to the SB one. Such states are characterised by the existence of a thin strip of ice-freeregion near the equator. Clearly, this possibility has key implications in terms of habitability andthe evolution of life on Earth. Different physical mechanisms have been proposed for explainingthe existence of such a state, based either on the role of a dynamical ocean [18] or the specificproperties of the albedo of sea ice [19]. In fact, in a previous work [20] we have found that,indeed, a third co-existing stable state, intermediate between the SB and the W climate, can befound in a climate model featuring a very simplified representation of the oceanic heat transport,so that one might expect that the existence of more than two competing attractors could be a ratherrobust property of the climate system.Indeed, the dynamical landscape of the Earth system might be even more complex than what isusually expected. A recent study [21], performed using a rather sophisticated climate model runusing aquaplanet boundary conditions (without continents), indicates the existence of at least fivecompeting climate states, ranging from a snowball to a very warm state without sea ice.Additionally, the physics and the chemistry of the climate system feature further complexitieswhen one considers even warmer conditions. For sufficiently large values of S ∗ , the W climatestate loses its stability as a result of the dramatic strengthening of the positive feedback associatedto the presence of water vapour in the atmosphere. Indeed, warm conditions favour, through thethermodynamic effect associated with the Clausius-Clayperon relation, the retaining capacity ofwater vapour of the atmosphere. The water vapour is a powerful greenhouse gas, as it is active in The project EDEN (http://project-eden.space/) combines ideas and methods in astrophysics, planetary sciences, andastrobiology to search for and characterize nearby habitable worlds. S ∗ nears the critical value S ∗ W → SB with the system being in the W state, the climatic enginebecomes more efficient, because larger temperature gradients are realised inside the domain. Suchan increased efficiency leads to a stronger atmospheric circulation, which is fuelled by tempera-ture gradients and tends to reduce them by transporting heat from warm to cold regions, actingas a non-trivial diffusion process. Such a non-linear equilibration mechanism acts as a negativefeedback and, broadly speaking, is a macroscopic manifestation of the second law of thermody-namics. One of the results of the heat transport performed by the atmospheric circulation is thestabilization of the ice-line. When S ∗ = S ∗ W → SB , the ice-albedo feedback becomes as strong asthe negative feedbacks of the system, and the system flips to SB state with the ice-line reachingthe equator. Similar mechanisms are in place when the system is in the SB state and S ∗ nears,instead, S ∗ SB → W . When S ∗ = S ∗ SB → W , the ice begins to melt near the equator, leading to a rapidpoleward retreat of the ice line.At the critical transitions the climate system is not anymore able to dampen fluctuations dueto (infinitesimal) external forcings. Using methods borrowed from transfer operator theory [25],the investigation of the behaviour of the same model used in [9] has indeed shown that when S ∗ nears S ∗ W → SB , the spectral gap - defined as the absolute value of the (negative) imaginary part ofthe subdominant Ruelle-Pollicott pole [26, 27] - of the transfer operator constructed in a suitablydefined reduced phase space vanishes. As a result, exponential decay of correlation is lost and thesystem experiences what is often referred to as critical slowing down [28]; see also Refs. [29, 30].Far from critical transitions, it has been shown [31–33] that it is possible to perform climatechange projections resulting from a time-dependent CO forcing using Ruelle response theory[34]. In the case of perturbations not depending explicitly on time, response theory allows one todescribe how the measure of the system changes differentiably with respect to small changes inthe dynamics of the system. In the case of time-dependent perturbations, response theory makesit possible to reconstruct the measure supported on the pullback attractor [35–37] (see also the5losely related concept of snapshot attractor [38, 39]) the non-autonomous system through a per-turbative approach around a reference state, which, in the case of the climate studies referred tohere, corresponds to the pre-industrial conditions. When we are nearing a critical transition, it isreasonable to expect a monotonic decrease of the spectral gap [40, 41]. Hence, in the vicinity ofthe critical transitions the presence of a vanishing spectral gap leads to having a vanishing radiusof expansion for response theory [42]. Indeed, it is expected that the (near) closure of the spectralgap is associated to a strongly enhanced sensitivity of the system’s statistics to perturbations [43].See [44] for a thorough discussion on the various regimes of climatic response to forcings and ofthe relationship between climate change and climate variability across various temporal scales. B. Melancholia States of the Climate System
A key question is what lies in-between the stable co-existing climates within the region of theparameters space where bistability is found. In simple models, it is often possible to identify un-stable solutions sitting in-between the two stable climates. Such unstable solutions are embeddedin the boundary between the two basins of attraction and, are roughly speaking, ice-covered up tothe mid-latitudes. These solutions are saddles because they attract orbits starting from initial con-ditions on the basin boundary but are not stable, as small generic perturbations pushing the orbitoutside the basin boundary with probability one and then lead to the system falling eventually intoone of the competing attractors [3, 45].When studying more comprehensive climate models featuring chaotic dynamics, things are, asdescribed in the next section, considerably more complex from a mathematical point of view, andthe individuation of the unstable saddle solutions is much harder [46–49]. Since these solutionsare unstable, they cannot be found by direct numerical simulation. In a previous investigation[20], we adapted the edge tracking algorithm [50, 51] introduced for constructing the edge states,i.e. the special solutions separating laminar from long-lived turbulent regimes of motion in afluid dynamical setting . We used a recursive technique of bisection on the initial conditions forshadowing trajectories on the basin boundary separating the two co-existing W and SB states, andmanaged to populate the corresponding saddles.The analysis was performed using an intermediate complexity climate model with O (10 ) de-grees of freedom. The saddles we found had the remarkable property of featuring chaos (we Note that in such systems, technically, there are no competing attractors, because the only true attractor is thelaminar state. A different approach based on control theory aims at finding unstable solutions by a feedback loop involvingchanges in the value of the control parameter defining the region of bistability [52, 53]. )b)FIG. 1. Bifurcation diagram for the coupled climate model studied in [20]. Panel a) Globally averagedocean temperature (cid:104) T S (cid:105) vs µ . Bistability is found for a large range of values of µ . Red continuous line:W attractor; blue continuous line: SB attractor; green continuous line: chaotic M state; purple continuousline: mean properties of the symmetry broken chaotic M state: red dashed line: warm side of the symmetrybroken chaotic M state; blue dashed line: cold side of the symmetry broken chaotic M state; green dottedline; transient symmetric chaotic M state; empty squares: warm side (red), cold side (blue) and averageproperties of the third attractor. The W-to-SB tipping point is located at µ ∼ . ; the SB-to-W tippingpoint is located at µ ∼ . ;. See [20]. (cid:104) T S (cid:105) and of the temperature difference be-tween low and high latitudes ∆ T S as a function of the relative solar irradiance µ = S ∗ /S ∗ , where S ∗ is the present day value. By focusing on the M states we have been able to show the existenceof much richer than previously thought dynamical landscape.As discussed in [20], up to µ ∼ . , the M state is characterised by longitudinal symmetryin its statistical properties, just as the boundary conditions of the system, are, indeed, longitudi-nally symmetric. The chaotic dynamics manifests itself as weather variability in a form not toodissimilar from the usual one observed in stable climates. Nonetheless, on long time scales, orbitsinitialised near the M states drift to either the W or the cold SB state, as a result of the dominatingpositive ice-albedo feedback. For µ ∼ . , we observed that the symmetric M state becomestransient, evolving very slowly (on a time scale scale much longer than the other ones typical ofthe system) into a symmetry-broken state, where very cold and very warm conditions co-exist,separated by two regions of very strong longitudinal temperature gradient. The two regions fea-ture rather different dynamical behaviour and the boundary between them rotates very slowly intime. The nontrivial bifurcation associated to such a symmetry break leads to dynamical regimesthat resembles chimera states in extensive systems [54, 55]. The third climate state mentionedbefore exists in a small parametric window near µ ∼ . [20].We remark that the dynamical systems viewpoint clarifies that the critical transitions for theW (snowball) state occurring when S ∗ approaches the critical value S ∗ W → SB ( S ∗ SB → W ) are associ-ated to the collision between the M state and one of the W (snowball) climates, according to thedynamical scenario of boundary crisis [48]. The system’s reduced ability to dampen fluctuationsnear the tipping points and the associated shrinking of the spectral gap described above can beseen, dynamically, as the result of the fact that the attractor attracts orbits less effectively in itsimmediate neighbourhood because of the presence of a nearby M state state [28]. C. This Paper: Goals and Main Results
Studying multistable systems in general, and the climate system in particular, using determin-istic autonomous dynamical systems faces two important issues, both resulting from the fact thatthe phase space is partitioned into disjoint invariant sets - the various basins of attractions and theirrespective boundaries. First, it is not possible to account for transitions between the co-existing8asins of attraction. Instead, transitions between distinct regimes of motion are observed in manysystems of interest. Additionally, one cannot establish an ergodic, physically relevant invariantmeasure, as the co-existing attractors are disjoint. Assigning a weight to each of them is, indeed,a highly arbitrary operation. Hence, one cannot answer the question of how likely it is for thesystem of interest to be found at a given time in a specific regime of motion. Instead, all the sta-tistical properties of an orbit are conditional on which invariant set its initial conditions belong to.Indeed, building upon the preliminary results reported in a short communication [56] and givingthem a much broader scope and setting them a in much more robust mathematical framework, wewant to address these points. For the benefit of the reader, we report below the main goals of ourinvestigation and the main findings presented in this paper.We will introduce a stochastic forcing to the model studied in [20] as a Gaussian perturbationwith variance σ modulating the intensity of the incoming solar radiation. This impacts the radia-tive budget of the planet in a very nontrivial way and, thanks to the ice-albedo feedback, acts asa multiplicative noise. Additionally, we conjecture, supported by physical and qualitative math-ematical arguments, that the combination of noise and nonlinear deterministic dynamics leads toa hypoelliptic diffusion process, i.e. the noise propagates to all degrees of freedom of the sys-tem [57]. The presence of a random forcing allows the system to perform transitions between theneighbourhoods of the deterministic attractors, crossing the basin boundaries that, in the unper-turbed case, are impenetrable; see some classical results in [58–60]. The consideration of a randomforcing will allow us to construct the ergodic invariant measure of the system by performing longintegrations, and investigate the properties of noise-induced transitions.Following [49, 61–63], we show that, in the weak noise limit and under suitable conditions, theinvariant measure can be written as a large deviation law with the following notable properties.The exponent is given by minus the quasi-potential function divided by σ / . We will show howto compute the quasi-potential from the drift and volatility fields, and show that the drift field canbe decomposed in two contributions having radically different dynamical meaning. Additionally,the quasi-potential is a Lyapunov function and provides a clear picture of the evolution of thesystem in the absence of noise. In the region of bistability, the quasi-potential has local minimaassociated to the attractors and has a saddle behavior at the M states of the deterministic system.In the region where only one stable state is realised in the deterministic case, the quasi-potentialhas just one local (and global) minimum, corresponding to the unique attractor.In the case of multistability, the logarithm of the average permanence time in a basin of at-9raction increases linearly with the product of /σ times the difference between the value of thequasi-potential at the M state through which the transitions takes place and its value at the cor-responding attractor. We also show that the stochastically averaged exit trajectories connect theattractor with the M state, which is indeed the most likely exit point from the deterministic basinof attraction. In the weak noise limit, such paths correspond to the instantons [59, 60, 64–66]. Inour simulations, we show that such an identification becomes more accurate as the intensity of thenoise is decreased. Exploiting this property, we also propose a new method for constructing the Mstates via direct numerical integration of stochastic differential equations.We discover that, since generally the local minima of the quasi-potential corresponding to thetwo attractors have different values, in the zero-noise limit only one attractor - the one corre-sponding to the global minimum of the quasi-potential - is populated. Nonetheless, an individualtrajectory may in fact persist near the competing metastable state for a very long time, as the per-manence time in the corresponding basin of attraction also diverges (but at a slower rate than forthe asymptotic state). As indicated by physical intuition, for low values of S ∗ , the noise selectsas limit measure the SB state, while for high values of S ∗ , the limit measure is the W attractor.The changeover takes place for a critical value of the relative solar irradiance µ = µ crit , where µ crit ≈ . . For such a value of µ , the equivalent of a first order phase transition for equilibriumsystems takes place in the system. The order parameter that more obviously captures the transitionis the globally averaged surface temperature.The paper is structured as follows. In Section II we summarise the main mathematical conceptsand ideas we have used to frame our investigation, to perform the data analysis, and to interpretour results. In Section III we report the modelling suite we have used, describe the data we haveproduced, and give information on how they have been processed. In Section IV we present anddiscuss our results. In Section V we propose and test numerically a new method for construct-ing the M states by looking at the intersections of instantons constructed by taking into accountthree different noise laws superimposed on to the same deterministic dynamics. In Section VI wesummarize and interpret our findings, describe the limitations of our work, and propose ideas forfuture investigations. In Appendix A, we speculate on the fact that the framework presented inthe paper is potentially suitable for clarifying the role played in complex historical processes by contingency , as discussed by Gould [67] in the context of evolutionary biology.10 I. MATHEMATICAL BACKGROUNDA. Geometry of the Phase Space: Attractors and M States
The investigation of systems possessing multiple steady states is an extremely active area ofinterdisciplinary research, encompassing mathematics, natural sciences, and social sciences. Therecent review by Feudel et al. [68], introducing a special journal issue on the topic, gives a rathercomplete overview of the state-of-the-art of the ongoing activities on this topic and provides sev-eral interesting examples. In the context of Earth sciences, it is more common to refer to criticaltransitions in multistable systems using the expression tipping points , which has received a greatpopularity following a paper by Lenton et al. [69].A possible way to introduce the mathematical background for multistable systems is the follow-ing. We consider a smooth autonomous continuous-time deterministic dynamical system definedon a smooth finite-dimensional compact manifold M (in what follows, a subset of R N ). We as-sume that the dynamical system is dissipative, so that the phase space continuously contracts andthe Lebesgue measure is not conserved. The orbit evolves from an initial condition x ∈ M attime t = 0 . We define x ( t, x ) = S t ( x ) as the orbit at a generic time t , where S t indicates the evo-lution operator. The corresponding set of ordinary differential equations written componentwiseis: ˙ x i = F i ( x ) , i = 1 , . . . , N, (1)where F ( x ) = d / d τ S τ ( x ) | τ =0 is a smooth N -dimensional vector field. We define such a dynam-ical system as multistable if it features more than one asymptotic states. Specifically, we meanthat there are two or more attractors Ω j , j = 1 , . . . , J , each a corresponding basin of attractionhaving a finite Lebesgue measure. Each of the attractors is an invariant set and is the support of aninvariant measure of the system. The asymptotic state where the trajectory falls into is determinedby its initial conditions, and the phase space is partitioned between the basins of attraction B l ofthe various attractors Ω j and the boundaries ∂B l , l = 1 , . . . , L of such basins.We assume, for simplicity, that within each basin of attraction an ergodic measure can be de-fined as the limit of the empirical measure constructed by averaging over infinitely long forwardtrajectories for Lebesgue almost all initial conditions in the the basin of attraction. In other terms,all the invariant measures of the system can be written as a convex sum of its j = 1 , . . . , J ergodiccomponents, each supported on the corresponding attractor Ω j , j = 1 , . . . , J .11e also assume that orbits initialized on the basin boundaries ∂B l , l = 1 , . . . , L are attractedtowards invariant saddles. In each basin boundary ∂B l there can be one or more of such sad-dles. We assume that each saddle features only one unstable direction. Such instability repelstrajectories initialised near the saddle towards either of the competing attractors.In general, one can be in the situation where the asymptotic dynamics on one or more of thecompeting attractors can be chaotic (meaning here that at least one Lyapunov exponent is positiveand unstable periodic orbits are dense). In some systems, chaotic dynamics can be realised on oneor more of the saddles embedded in the basin boundaries [46–48, 70]. These are what we refer toas M states. M states are non-trivial geometrical objects being the support of a typically nontrivialmeasure [49].A fascinating aspect of multistable systems is the following. If the first Lyapunov exponent of achaotic saddle is larger than the inverse of the life time of the saddle state itself (which measures therate at which orbits initialised near the saddle state are repelled towards the two nearby asymptoticstates), then the basin boundary separating the basin of attraction of the two asymptotic sets hasco-dimension strictly smaller than one. In two-dimensional maps, it has been proved that in thiscase the basin boundary is not a manifold, but rather a rough fractal [46, 49, 70]. The authorshave recently proposed a generalisation of this result to the case of N -dimensional maps [71].In [20], despite high-dimensionality of the system, one could detect that the co-dimension of thebasin boundary is strictly smaller than one. In fact, such co-dimension was found to be very closeto zero, as a result of time scale difference between the most relevant instabilities: the climaticinstability due to the ice-albedo feedback acting near the M state is much slower than the fastweather-like instability associated to baroclinic processes occurring on the M state. This impliesthat near the basin boundary there is basically no predictability of the second kind in the sense ofLorenz [72]. The basin boundary is, de facto , a grey zone, and it is difficult to assess where orbitsinitialised near the boundary will end up. B. Impact of Stochastic Perturbations: Invariant Measure and Noise-induced Transitions
The goal of our investigation here is to analyse the impact of imposing a random forcing onthe deterministic dynamics discussed above. Processes that can be described as a noise-inducedescape from an attractor have long been studied in the natural sciences; see [58–60]. We generaliseEq. 1 by adding a stochastic component. We then consider a stochastic differential equation (SDE)12n Itˆo form written as dx i = F i ( x ) dt + σs ( x ) ij dW j , (2)where x , F ∈ R N , F and s ij ∈ R N × N are smooth, σ ≥ , dW j is the increment of an N − dimensional brownian motion, and C ij ( x ) = s ik ( x ) s jk ( x ) is the noise covariance matrix.We consider the case of a hypoelliptic diffusion process. This amounts to assuming that while thecovariance matrix of noise can be singular, the drift term modified according to the Stratonovichconvention and the columns of the matrix s satisfy the so-called H ¨ormander condition, i.e., theLie algebra generated by them has dimension N everywhere [73]. As a result of this, a smoothinvariant density with respect to Lebesgue is realised because the noise is propagated to all thecoordinates through the drift term; see [57] . We discuss below in Sect. III B why such a mathe-matically important assumption can be heuristically justified in the specific case studied here.Taking inspiration from the Freidlin-Wentzell [61] theory and modifications thereof [49, 62,63], in the weak-noise limit σ → we seek a special functional form for the invariant measure.Indeed, we look for a large deviation law: Π σ ( x ) ∼ exp (cid:18) − x ) σ (cid:19) , (3)where the rate function Φ( x ) is referred to as quasi-potential, and we have neglected the pre-exponential term. Specifically, the symbol ∼ in Eq. 3 implies that Φ( x ) = − / σ → σ log Π σ ( x ) .The function Φ( x ) can be obtained as follows. We solve the stationary Fokker-Planck equationcorresponding to Eq. 2: ∂ j J j ( x ) = 0 , J j ( x ) = − F j ( x )Π σ ( x ) + σ ∂ i ( C ij ( x )Π σ ( x )) , (4)where J is the current density. We then consider the weak noise limit, and use as ansatz theexpression given in Eq. 3. We obtain the following Hamilton-Jacobi equation [74]: F i ( x ) ∂ i Φ( x ) + C ij ( x ) ∂ i Φ( x ) ∂ j Φ( x ) = 0 . (5)This equation allows one to express Φ in terms of the drift and volatility fields introduced in Eq.2. The quasi-potential Φ can also be computed by solving the variational problem associated withthe Freidlin-Wentzell action [75]. Finally, alternative routes for computing Φ have been proposedin [76, 77]. Assuming an elliptic diffusion process is extremely restrictive because it requires all degrees of freedom to be drivenby Gaussian noise. Φ is far from trivial, yet of great interest in many applications; seee.g., [78] for the case of biological systems. Brackston et al. [79] have recently proposed an algo-rithm for estimating Φ in the case that the governing equations are polynomial and involves solvingan optimization over the coefficients of a polynomial function. Instead, Tang et al. [80] proposed avariational method for estimating in the populations corresponding to each deterministic attractorwithout resorting to computing the invariant measure.Following [62, 63], we now describe the dynamical meaning of Φ . Indeed, solving the previousHamilton-Jacobi equation corresponds to the fact that it is possible to write the drift vector field asthe sum of two vector fields: F i ( x ) = R i ( x ) − C ij ( x ) ∂ j Φ( x ) (6)that are mutually orthogonal, so that R i ( x ) ∂ i Φ( x ) = 0 . In the case Eq. 2 describes a thermody-namical system near equilibrium, R defines the time reversible dynamics, while F − R defines theirreversible, dissipative dynamics [81]. One finds that d Φ( x ) /dt = − C ij ( x ) ∂ i Φ( x ) ∂ j Φ( x ) + R i ( x ) ∂ i Φ( x ) = − C ij ( x ) ∂ i Φ( x ) ∂ j Φ( x ) . (7)As a result of this, Φ has the role of a Lyapunov function whose decrease describes the convergenceof an orbit to the attractor. Specifically, Φ( x ) has local minima at the deterministic attractors Ω j , j = 1 , . . . , J , and has a saddle behaviour at the saddles Π l , l = 1 , . . . , L . If an attractor (saddle) ischaotic, Φ has constant value over its support, which can then be a strange set [62, 63].Note that in the standard case of dynamics taking place in an energy landscape defined by a(confining) potential U ( x ) and noise correlation matrix proportional to the identity (obtained bysetting F i ( x ) = − ∂ i U ( x ) and C ij ( x ) = in the previous equations) one has Φ = U . Additionally,one derives ˙ U ( x ) = − ∂ i U ( x ) ∂ i U ( x ) < and U ( x ) is a Lyapunov function, and one recovers anequilibrium state, where detailed balance applies and, by definition, the current vanishes ( J = 0 ) .We remark that the function Φ( x ) is defined globally but is not, in general, twice differentiableeverywhere. Indeed, discontinuities in its first derivatives are present if a) the Hamiltonian associ-ated with the Hamilton-Jacobi equation given in Eq. 5 is not integrable (non-integrability being thetypical situation), and b) if the system features more than one co-existing attractors. These latterdiscontinuities are of little practical relevance because they appear only for values of Φ larger thanthose at the saddles [82], for reasons that will become apparent below. Note that, in general, we can have an equilibrium state if and only if the the drift term has a gradient structure withrespect to the metric defined by the noise covariance tensor [30]. . Noise-induced Escape from the Attractor The quasi-potential Φ is key for determining the statistics of noise-induced escape from a givenattractor. Indeed, the probability that an orbit with initial condition in B j does not escape from itover a time p decays as: P j,σ ( p ) = 1 τ j,σ exp ( − p/τ j,σ ) , τ j,σ ∼ exp (cid:18) j σ (cid:19) , (8)where τ j,σ is the expected escape time and ∆Φ j = Φ(Π l ) − Φ(Ω j ) is the lowest quasi-potentialbarrier height [49], i.e. Φ has the lowest value in Π l compared to all the other saddles neighbouring Ω j . In general, one may need to add a correcting prefactor to P j,σ ( p ) [49].Note that τ j,σ given in Eq. 8 does not contain the pre-exponential factor. Ref. [66] provided anexpression for such pre-exponential factor for general non-equilibrium diffusion processes underthe assumption that attractors and saddles are simple points, thus generalising what is given in [83].As we aim at treating also a more general setting for the geometry of attractors and saddles, wepay below the price of having to take the pre-exponential factors as phenomenological parametersthat one can find from experiments or numerical simulations [84]. We also remark that, in thezero-noise limit, the transition paths outside a basin of attraction follow the instantons. Instantonsare defined as solutions of dx i /dt = ˜ F i ( x ) = R i ( x ) + C ij ( x ) ∂ j Φ( x ) (9)that connect a point in Ω j to a point in Π l . Instantonic trajectories have a reversed component ofthe gradient contribution to the vector field compared to regular - relaxation - trajectories.Let’s now take the simpler case of bistable systems, where we have two attractors Ω , Ω , andone saddle Π . We can then express the average transitions times as follows: τ → σ ∝ exp (cid:18) ) − Φ(Ω ) σ (cid:19) , (10) τ → σ ∝ exp (cid:18) ) − Φ(Ω ) σ (cid:19) , (11)so that τ → σ τ → σ ∝ exp (cid:18) ) − Φ(Ω ) σ (cid:19) . (12)This implies that, in the weak noise limit, both escape times diverge, but the escape time out ofthe attractor corresponding to the lower value of the quasi-potential diverges faster. Note that one15an expect the proportionality constant in Eq. 12 to be O (1) . Taking a maximally coarse-grainedview on the problem, where we consider the state as represented by the populations P ,σ , P ,σ ofthe neighbourhood of the two attractors, we can write the following master equation: ˙ P ,σ = − P ,σ τ → σ + P ,σ τ → σ (13) ˙ P ,σ = − P ,σ τ → σ + P ,σ τ → σ . (14)The master equation above makes sense if one assumes the presence of clear timescale separationbetween the relaxation motions near each attractor and those across the saddle, which dependscritically on the presence of weak noise [85, 86]. At steady state, we obtain that P ,σ P ,σ = τ → σ τ → σ ∝ exp (cid:18) ) − Φ(Ω ) σ (cid:19) . (15)Equation 15 could be obtained by integrating the invariant measure given in Eq. 3 in the neigh-borouhood of the attractors and taking a saddle point approximation. Additionally, Eq. 15 impliesthat in the weak-noise limit only one of the two deterministic attractors will be populated, andspecifically the one where the quasi-potential has lower value. We remark that two different noiselaws differing for the correlation matrix C acting on top of the same drift field will define twodifferent quasi-potentials, see Eq. 5. As a result of that, they will in general feature a differentselection of the dominating population in the zero noise limit. One can easily extend the mas-ter equation defined above to the case where multiple states and multiple paths of transitions arepresent. Finally, note that in [87] the mathematical framework described in this section has beenused to study stochastic resonance for general non-equilibrium systems. III. NUMERICAL MODELLING
The climate model considered here is constructed by coupling the primitive equations atmo-spheric model PUMA [88] with the Ghil-Sellers energy balance model [3], the latter describingsuccinctly the meridional oceanic heat transport. It has been already presented in [20] with thename of PUMA-GS, but we report here again its formulation in order to elucidate the role ofstochastic forcing, which was absent in the previous version. The stochastic forcing is added as afluctuating term modulating the value of the incoming radiation determining the energy input intothe system. 16 . The Atmospheric Component
The atmospheric component of the PUMA-GS model is provided by PUMA [88], which con-sists of a dynamical core: the dry hydrostatic primitive equations on the sphere (mapped laterallyby the latitude φ and longitude λ ), solved by a spectral transform method (only linear terms areevaluated in the spectral domain, nonlinear terms are evaluated in grid-point space). The equa-tions for the prognostic state variables, the vertical component (with respect to the local surface)of the absolute vorticity ζ = ξ + 2 ν Ω E (where ξ is the vertical component of the relative vortic-ity, ν = sin φ , and Ω E = 2 π /day is the angular frequency of the Earth rotation) the (horizontal)divergence of the velocity field D , the (atmospheric) temperature T a = ¯ T a + T (cid:48) a (separated intoa time-independent arbitrary reference part ¯ T a and anomalies T (cid:48) a ), and the logarithmic pressure(normalized by the surface pressure p s ) σ = ln p/p s , read as follows: ∂ t ζ = s ∂ λ F v − ∂ ν F u − τ − f ξ − K ∇ ξ, (16) ∂ t D = s ∂ λ F u + ∂ ν F v − ∇ [ s ( U + V ) / T a ln p s ] (17) − τ − f D − K ∇ D,∂ t T (cid:48) a = s ∂ λ ( U T (cid:48) a ) − ∂ ν ( V T (cid:48) a ) + DT (cid:48) a − ˙ σ∂ σ T a (18) + κT a ω/p + τ − c ( T R ( T S ) − T a ) − K ∇ T (cid:48) a ,∂ t ln p s = − s ∂ λ ln p s − V ∂ ν ln p s − D − ∂ σ ˙ σ, (19) ∂ ln σ Ψ = − T a , (20)where s = 1 / (1 − ν ) , F u = V ζ − ˙ σ∂ σ U − T (cid:48) a ∂ λ ln p s , F v = − U ζ − ˙ σ∂ σ V − T (cid:48) a s − ∂ ν ln p s , U = u cos φ , V = v cos φ , u , v being respectively the horizontal and vertical wind velocitycomponents, and Ψ is the geopotential height. Equations 16,17, and 19 express the conservationof momentum, Eq. 18 expresses the conservation of energy, and Eq. 20 is the equation of state.A number of simple parametrizations are adopted in order to improve the realism and the sta-bility of the model. Firstly, the hyperdiffusion operator K ∇ is added to the equations of vorticity,divergence and temperature, to represent subgrid-scale eddies. Secondly, large-scale dissipationof vorticity and divergence is facilitated by Rayleigh friction of time scale τ f . Thirdly, the physicsof diabatic heating due to radiative heat transport is parametrized by Newtonian cooling: the tem-perature field is relaxed (with a time scale τ c ) towards a reference or restoration temperature field T R , which can be considered a radiative-convective equilibrium solution. We adopt the following17imple expression for the restoration temperature [88]: T R = ( T R ) tp + (cid:113) [ L ( z tp − z ( σ )) / + S + L ( z tp − z ( σ )) / , (21) ( T R ) tp = (cid:104) T S (cid:105) − ¯ Lz tp , (22) L ( λ, φ ) = ∂ z T R = ( T S ( λ, φ ) − ( T R ) tp ) /z tp , (23)where ( T R ) tp and z tp are the temperature and height of the tropopause, respectively, L ( ¯ L ) is the(average) lapse rate, (cid:104) T S (cid:105) is the globally averaged surface temperature, and z ( σ ) is determinedby an iterative procedure [88]. The above expressions indicate that the restoration temperatureprofile is anchored to the surface temperature T S . However, as Eq. 22 indicates, T R at any onepoint on the sphere is determined by not only the local (dynamical) surface temperature, but alsothe global average (cid:104) T S (cid:105) . We note that T a ( σ = 1) is obtained by linear extrapolation, accordingto T a ( σ = 1) ≈ T a ( σ = 0 .
9) + η ( T S − T R ( σ = 0 . , < η < . With η = 1 the couplingterm is k ( T a ( σ = 1) − T S ) ≈ k ( T a ( σ = 0 . − T R ( σ = 0 . . Generally T a ( σ = 1) − T S (cid:54) = 0 (laterally inhomogeneous heating), but (cid:104) T a ( σ = 1) (cid:105) = (cid:104) T S (cid:105) , where the overbar denotes averagingwith respect to time.For our setup we choose: K − = 0 . days, τ c = 30 days, τ f = 1 day, ¯ L = 0 . K/m, z tp = 12000 m, and k = 10 − . We also adopt a coarse resolution of T21 (i.e., the series ofspherical harmonics are triangular-truncated at total wave number 21). This implies the optimalnumber of Gaussian grid points: N lon = 2 N lat = 64 . Finally, we consider N lev = 10 verticallayers and consider a vanishing orography, so that we have a zonally-symmetric configuration.The equations are integrated numerically using a ∆ t = 1 [hour] time step size. B. The Ocean Component and the Stochastic Forcing
The surface temperature T S is taken to be governed by the a 2D version of the GS EBM [3, 45].This model includes a simple yet effective representation of the ice-albedo feedback, and basi-cally defines the slow manifold of the coupled atmosphere-ocean system. The partial differentialequation describing the evolution of the ocean surface temperature field T S = T S ( t, φ, λ ) is: ∂ t T S ( t, φ, λ ) = µ I ( φ ) C ( φ ) S ∗ − α ( φ, T S )) − O ( T S ) C ( φ ) − ¯ D φ [ T S ] C ( φ ) + χ [ T S , T A ] C ( φ ) + s.f., (24)where S ∗ is the present solar irradiance , µ = S ∗ /S ∗ as introduced in Sect. I, while the heatcapacity C ( φ ) and the geometrical factor I ( φ ) are explicitly dependent on φ only, thus enforcing The factor 4 emerges as a result of the geometry of the Earth-Sun system [89] α depends on φ and, critically, on T S , witha rapid transition from strong albedo for low values of T S ( α max = 0 . ) to weak albedo for T S (cid:38) K ( α min = 0 . ), which fuels the positive ice-albedo feedback. Additionally, O is theoutgoing radiation per unit area, expressed as a monotonically increasing function of T S (this isresponsible for the negative Boltzmann feedback, taking into account also the greenhouse effect), ¯ D φ is a diffusion operator parametrizing the meridional heat transport, and χ describes the heatexchange with the atmosphere. See [20, 45] for further details.Finally, the last term on the right hand side s.f. is the stochastic forcing, which is introducedas a random modulation of the solar irradiance given by µS ∗ . Hence, we have: s.f. = σs ( T s , φ, λ ) dWdt = σµ I ( φ ) C ( φ ) S ∗ − α ( φ, T S )) dWdt , (25)where σ controls the intensity of the noise, s defines the noise law, and dW is the increment ofa one-dimensional Wiener process. Since s depends explicitly on T S via the term α ( φ, T S ) , weare dealing with a multiplicative noise law. We consider the Itˆo convention for noise, so that our(discretised) equations are in the form of Eq. 2. See a discussion on the relevance of the chosenconvention in Sect. IV E. Adding a Gaussian random variable of variance σ at each time step ∆ t ( hour ) of the model amounts to considering that, on the time scale τ = N × ∆ t , the relativefluctuation of the solar irradiance scales as σ τ = σ/ √ N .As mentioned in Sect. II B, our approach requires assuming the validity of the hypoelliptic-ity condition. In order to prove this, we should test the H¨ormander condition for the evolutionequations of the model. This is of great relevance but is beyond the specific scope of this paper,while indeed deserving a separate and accurate investigation. However, as discussed below wecan heuristically understand why, indeed, it is reasonable to assume that stochastic forcing actingon the oceanic surface temperature propagates to all degrees of freedom of the coupled system, asusually implicitly assumed in basically any numerical study of stochastically forced geophysicalflows .We can first approach this problem by looking at the structure of the evolution equations. Thestochastic forcing given in Eq. 25 impacts directly the T S field, as shown in Eq. 24. The T S fieldsdetermines the restoration temperature T R , see Eq. 21. In turn, the restoration temperature impactsthe anomaly of the atmospheric temperature field T (cid:48) a (see Eq. 18), which in turn affects the vorticityfield ζ - Eq. 16 - and divergence field D - Eq. 17. Finally, anomalies in D impact the surface Obviously, numerical truncation introduces some additional noise on all degrees of freedom of the system. p s , as clear from Eq. 19. The nonlinear terms corresponding to advective processes onthe right hand side of each Eq. 16-19, which contain two more of the above mentioned fields,make sure that noise propagates across all scales in each field. This latter point could be betterunderstood by taking a truncated Fourier representation of Eqs. 16-19, which, in fact, closelycorresponds to the actual formulation of the numerical model implemented here. Concluding, nodynamical or thermodynamical field and no scale within each field is insulated from the noise,even if the covariance matrix of the noise law is extremely singular, as noise impacts directly onlya small fraction of the degrees of freedom of the system.On more physical grounds, one can observe that the ocean surface temperature drives therestoration profile of the atmospheric temperature, and that fluctuations of such a profile mod-ulate the atmospheric instabilities, whose energy cascades down to the smallest scales resolved bythe model; see [44, 90, 91] for a detailed treatment of the energetics of the climate system. IV. RESULTS
We first treat in detail three cases inside the region of bistability depicted in Fig. 1, namely µ = 0 . (close to the tipping point µ W → SB ), µ = 1 . (corresponding to present-day solar irra-diance), and µ = 1 . (in the parametric region where the M state undergoes a symmetry-breakbifurcation). We then construct the weak-noise limit of the invariant measures for all the values of µ in the region of bistability. We remark that the results shown in Figs. 2, 3, and 5a) have alreadybeen reported in the short communication [56], but we deem extremely useful to present them hereas well because they are now part of a much more complex, coherent, and detailed narrative. A. Escapes from Basins of Attraction and Instantons
In the case of µ = 0 . , we first perform a set of simulations with noise of different intensityrangiW ng from σ τ = 0 . to σ τ = 1 . , with τ = 100 years ( y ). For each value of thenoise intensity, we initialise 50 trajectories in the basin of attraction of the W climate and studythe statistics of the escape time to the SB attractor. When the transition takes place, we stop theintegrations. We observe (not shown) that for each value of σ τ the escape times are to a goodapproximation exponentially distributed, thus obeying Eq. 8; the process of transition behaveslike a Poisson process. The results on the expectation value of the transition times are presented20 IG. 2. Statistics of the escape times p for the noise-induced W → SB transitions for various noisestrengths. Each coloured dot corresponds to an observed escape time p . The estimate of the expectedescape times ¯ τ σ are indicated by the black dots connected by the red line. The slope of the straight linefit gives the potential difference described in Eq. 10. See [84] for an optimal algorithm for estimating thepotential difference. Reproduced from [56]. in Fig. 2, where we show that, indeed, τ σ agrees with the prediction of Eq. 8. Hence, it is possibleto define the difference between the value of the quasi-potential Φ at the M state and at the Wattractor as the slope of the straight line. For reference, we have that for σ y = 0 . the averageescape time is about . × y . We can predict that the escape rate increases to about . × y when σ y ∼ . . We have that the slope of the straight line in Fig. 2 gives ∼ M ) − Φ( W )) .Therefore, we show that it is indeed possible to estimate quantitatively the properties of the quasi-potential also in a very high-dimensional dynamical system like the one considered here. Theoperation can be repeated for all the other values of µ in the range of multistability and for theprocesses of escape from the SB attractor, but we do not pursue here a systematic study of this.We then wish to look at the paths corresponding to the transitions. Following the discussion in[20, 45], we choose to consider the reduced phase space spanned by the globally averaged surfaceocean temperature (cid:104) T S (cid:105) and by the meridional temperature difference ∆ T S , defined as the differ-ence between the spatially averaged ocean temperature field between the Equator and ◦ N and21 IG. 3. Main graph: Logarithm of the transient density ˜ ρ in the reduced space (∆ T S , (cid:104) T S (cid:105) ) (in units of K in both axes), with indication of the actual position of the W attractor (red dot) and the M state (green dot)for µ = 0 . . We have used σ y = 1% . The W → SB approximate instanton is indicated. Top left inset:pdf along the path of the instanton ( (cid:104) T S (cid:105) on the x-axis). Reproduced from [56] . between ◦ N and the North Pole. This reduced phase space provides a minimal yet physically in-formative viewpoint on the problem, because it is directly linked with the main physical processesoccurring in the climate model: • The average surface temperature (cid:104) T S (cid:105) is directly associated to the positive ice-albedo feed-back and the negative Bolzmann radiative feedback; • The meridional temperature difference ∆ T controls the meridional heat transport performedby the ocean, as a result of the diffusive law we insert into its evolution equation; • The meridional temperature gradient ∆ T also controls the meridional heat transport per-formed by the atmosphere, as a result of the mechanism of baroclinic instability [92].Figure 3 depicts, for the case σ τ = 1 . , the transient two-dimensional distribution function ˜ ρ constructed using a frequentist approach using the 50 simulations described above, where thestatistics is collected only until the W → SB transition is realised. The distribution we obtaincannot be interpreted as an approximation of the invariant measure, because the integrations are22 IG. 4. Main graph: density in the projected phase space (∆ T S , (cid:104) T S (cid:105) ) (in units of K in both axes), withindication of the actual position of the W attractor (red dot), SB attractor (blue dot), M state (green dot)for µ = 0 . . The W → SB and SB → W approximate instantons are also indicated. We have used σ y = 1 . . Top left inset: marginal pdf with respect to (cid:104) T S (cid:105) . Bottom right inset: marginal pdf withrespect to ∆ T S . Center left inset: pdf along the path of the two instantons. stopped after the transitions. Nonetheless, it is apparent that the transitions take prominentlyplace in a very narrow band linking the W attractor and the M state . In order to obtain a betterunderstanding of the transition paths, we construct an estimate of the instanton linking the Wattractor to the M state and associated to the W → SB transitions by conditionally averaging thetrajectories according to the value of (cid:104) T S (cid:105) . To a good approximation, the instanton connects theW attractor to the M state, and follows a path of decreasing density. We do not find evidence ofdifferent paths for the trajectories leading to an escape and the relaxation trajectories, which is,instead, a typical signature of non-equilibrium [93]. This can be explained by considering [45],where it is shown that the ocean model evolve to a good approximation in an energy landscape.See also Section V. Note that, even if the W attractor and the M state look like dots, they have, in fact, a finite (yet very small) size,because they are both chaotic (see caption of Fig. 1). Here we are considering oceanic variables, which feature avery small variability in the deterministic chaotic case. )b)FIG. 5. Panel a) Main graph: density in the projected phase space (∆ T S , (cid:104) T S (cid:105) ) (in units of K in both axes),with indication of the actual position of the W attractor (red dot), SB attractor (blue dot), M state (greendot) for µ = 1 . The W → SB and SB → W approximate instantons are also indicated. We have used σ y = 1 . . Top left inset: marginal pdf with respect to (cid:104) T S (cid:105) . Bottom right inset: marginal pdf withrespect to ∆ T S . Center left inset: pdf along the path of the two instantons. Reproduced from [56]. Panelb): same as panel a), with σ y = 1 . . IG. 6. Estimate of the instantons for µ = 1 obtained using σ y = 1 . (dashed lines), σ y = 1 . (dash-dotted lines), and σ y = 2 . (dotted lines) (units of K in both axes). The red (blue) lines show theestimates for the W → SB ( SB → W ) instanton. The dots indicate the actual position of the W attractor(red dot), SB attractor (blue dot), and M state (green dot). The estimate of the instanton improves as theintensity of the noise is reduced. B. Construction of the Invariant Measure
The information contained in Fig. 3 is limited because we are studying only W → SB escapeprocesses, and we do not allow for the establishment of an invariant measure. The problem lies inthe fact that the quasi-potential minimum associated to the SB attractor is much deeper than theone associated to the W attractor, so that the average escape time associated to SB → W tran-sitions is prohibitively long for the range of (rather weak) noise intensities used for constructingFig. 3. In fact, in order to be able to construct the invariant measure of the system, we need toobserve a sufficient number of W → SB and SB → W transitions, in order to be sure that wehave collected a satisfactory statistics; see also the master equations for the populations given inEqs. 13-14. We next increase the noise intensity by setting σ y = 1 . , so that, in an integrationlasting about . × y , we observe 10 SB → W and W → SB transitions. The number oftransitions is low because it is extremely hard to escape from the SB state.Our results are shown in Fig. 4. We portray the logarithm of the projection of the invariant mea-25 IG. 7. Main graph: density in the projected phase space (∆ T S , (cid:104) T S (cid:105) ) , with indication of the position ofthe W attractor (red dot), SB attractor (blue dot), transient M state (green dot) for µ = 1 . (units of K inboth axes). The squares indicate the properties of the asymptotic M state: W sector (red square); cold sector(blue square); average (magenta square). The estimates of the W → SB and SB → W instantons are alsoindicated as the red and the blue line, respectively. We have used σ y = 1 . . Top left inset: marginalpdf with respect to (cid:104) T S (cid:105) . Bottom right inset: marginal pdf with respect to ∆ T S . Center left inset: pdf alongthe path of the two instantons. sure on the (∆ T S , (cid:104) T S (cid:105) ) plane; we refer to this also as the two-dimensional probability distributionfunction (pdf). We have that the peaks of the pdf are in good agreement with the position of the W and SB attractors, as implied by the large deviation result presented in Eq. 3. The agreementis even clearer when considering the two marginal pdfs, constructed by projecting the invariantmeasure on the one dimensional spaces defined by ∆ T S and (cid:104) T S (cid:105) . Note that the peak over the W state is hardly noticeable because the occupation of the state is extremely low (less than 5% of thetotal); compare with the cases studied below where higher values of µ are considered. We are alsoable to construct both the W → SB and the SB → W instantons, whose starting points agreeremarkably well with the position of the W and the SB attractors, respectively, while their finalpoints are in good agreement with the position of the M state. By constructing the pdf along theinstantons, we find that they follow a path of monotonic descent (indeed, they follow closely the26rests of the pdf), with the minimum located at the M state. Again, this last property can be hardlyvisualised in Fig. 3 for the W → SB instanton (see inset), because the population near the W state is quite small.Next, we repeat the analysis for µ = 1 ; results are shown in Figs. 5a)-b). In Panel a) weconsider an integration with σ y = 1 . lasting about . × y and characterised by 41 SB → W and W → SB transitions. We have an occupation of about for the W basin ofattraction, and of about for the SB basin of attraction; additionally, we have τ W → SBσ ∼ y and τ SB → Wσ ∼ y . Also in this case, the projection of the invariant measure in the (∆ T S , (cid:104) T S (cid:105) ) plane shows that there is good agreement between the position of the peaks of the pdfs and theattractors, and that the estimates of the instantons connect attractors and M states with a goodprecision. It is also clear that the instantons follow a path of descent in terms of probability, asshown by the central inset.In Panel b) we show the results of repeating the analysis for σ y = 1 . . In this case thesimulation lasts about . × y and we obtain 73 SB → W and W → SB transitions, andwe can draw similar conclusions as in Panel a) regarding the relative position of the attractors,of the M state, and of the instantons. The marginal pdfs are clearly less peaked than in Panel a);the occupancy rate changes slightly with respect to the previous case: it is about for the W basin of attraction, and of about for the SB . Instead, the average escape times change moresubstantially, and can be estimated as τ W → SBσ ∼ y and τ SB → Wσ ∼ y . These last two resultsindicate that the difference between the value of the quasi-potential at the two competing attractorsis relatively small. We will explore this matter in Sect. IV D, where we will try to deduce wherethe quasi-potential reaches its absolute minimum for each value of µ in the bistable region.It is worth looking more in detail at how the estimate of the instantons is impacted by theintensity of the noise used in the simulation. As instantons are defined in the weak-noise limit,we would expect that one achieves higher precision when weaker noise is used. This is confirmedby the results shown in Fig. 6: we have that the estimates of the instantons obtained using lowernoise intensity come closer to the attractors and to the M state. Nonetheless, also the instantonsobtained for very strong noise are still relatively accurate.27 . Instantons and Transitions across the Symmetry-broken Melancholia State We next examine the noise-induced transitions for µ = 1 . . This case is quite interestingbecause, as discussed in [20] and reported in Fig. 1, the longitudinally-symmetric M state istransient with a very long life time, and slowly evolves into a symmetry broken M state featuring arelatively cold and a relatively warm region, separated by two small regions with large longitudinaltemperature gradient at all latitudes. It seems relevant to test whether noise-induced transitionstake place through the transient, symmetric M state or the true, symmetry-broken one. Results areshown in Fig. 7, where we use data from a simulation lasting about y with σ y = 1 . . Weobserve only 6 transitions in both directions. This time, as opposed to the case of µ = 0 . , thefigure is so low because it is extremely hard to escape from the W state. We estimate the escapetimes as τ W → SBσ ∼ y and τ SB → Wσ ∼ y . We portray the logarithm of the invariant measureof the system in the usual projected space, and the estimates of the instantons. We first observethat in this case most of the density is concentrated around the W attractor, and a nontrivial relationexists between the paths of the instantons and the dynamical structures on the basin boundary. It isclear that the transient M state plays the role of the gateway of transitions as seen in the previouscases, despite its transient nature. The noise-induced transitions do not go through the actual Mstate (magenta square), while they seem to go thorough states resembling the properties of the W(red square) and cold (blue square) sectors in the asymptotic M state. This feature might resultfrom the consideration of noise with finite (and not infinitesimal) strength: the quasi-potential nearthe transient M state might be just barely higher than that of the asymptotic M state (and possiblywith a more favourable pre-exponential factor), so that a small but finite noise perturbation mightpush an orbit near the transient M state into the other basin of attraction. In order to address thispoint and find instantonic paths connecting the attractors with the true M states, one might need toresort to using more sophisticated numerical techniques. In particular, one should consider usingrare events algorithms [94, 95] to rigorously construct instantonic trajectories [96]. This is beyondthe current abilities of the authors but definitely deserves attention in future studies. D. Selection of the Limit Measure in the Weak-noise Limit and First-order Phase Transition
Equation 15 indicates that, in the weak-noise limit, all of the measure will be concentrated onthe attractor featuring the lowest value for the quasi-potential Φ . Since for µ < S ∗ W → SB /S ∗ only28he SB state is realised, physical intuition suggests that for low values of µ within the range ofbistability, the SB attractor should contain all the mass in the limit of weak noise, as one can alsoanticipate by looking at Fig. 3. Conversely, one expects that for high values of µ within the rangeof bistability, the W should be dominant in the weak noise limit. It is reasonable (yet far fromobvious or rigorous) to expect that there should be a critical value of µ = µ crit separating the tworegimes. Following [20], we consider 18 equally spaced values of µ ( ∆ µ = 0 . ) within themultistable regime. For each of these values of µ (excluding the case of µ = 1 . , where threestable states are realised) the fraction of the population residing within the basin of attraction ofthe deterministic W attractor P W,σ ( µ ) and its complement, residing in the basin of attraction ofthe SB attractor P SB,σ ( µ ) . Related results are shown in Fig. 8. We show how the probabilitydistribution of the variable (cid:104) T S (cid:105) depends on µ for three different noise levels: σ y = 1 . (Panela), σ y = 1 . (Panel b), and σ y = 2 . (Panel c). In these panels we superimpose thebifurcation diagram reported in Fig. 1a). We observe that as the noise is reduced, for all values of µ the distributions are a) more peaked around the attractors; and b) one of the attractors becomesclearly dominant.In Panel d) the plot for each value of µ the integral of the pdfs reported in the three panelsa), b), c) up to the values of (cid:104) T S (cid:105) corresponding to the M state (green continuous and dottedlines). To a very good degree of approximation, this corresponds to the integral of the invariantmeasure over the support of the deterministic basin of attraction of the SB climate. We obtainthat for decreasing values of the noise intensity, the emerging invariant measure converges to thedeterministic measure supported on the SB attractor for µ ≤ µ crit ≈ . , while the invariantmeasure converges to the deterministic measure supported on the W attractor for µ ≥ µ crit ≈ . . The absolute minimum of the quasi-potential Φ is realised in the W attractor for µ ≥ µ crit and in the SB attractor for µ ≤ µ crit . The changeover is, curiously, quite close to the referencecase µ = 1 , for which the weak-noise limit of the measure is given by the SB state.We add a note on the uncertainty associated to the figures reported in Fig. 8d). We remark thatfor all values of µ and σ we have used simulations lasting at least y . For very low ( ≤ . )and very large ( ≥ . ) values of µ in the considered range, the simulation length does not allowfor observing more than few transitions for the two lowest considered noise levels. Therefore,according to the fact that the transitions occur following a Poisson law, one expects in this rangean uncertainty on the figures reported in Fig. 8 of the order of the values of the smaller between P SB,σ ( µ ) and P W,σ ( µ ) . This corresponds, in fact, to a low uncertainty, because most of the mass29s concentrated near one of the two deterministic attractors. The uncertainty is quite small for . ≤ µ ≤ . , because, in all cases, we observe relatively many transitions. The uncertaintyin this range can be safely estimated to be below 5%. Summarising, while the values reported inFig. 8d) might have some non-negligible uncertainties, it seems that the estimate of µ crit is quiterobust.Finally, we have verified that for all values of µ the escape time τ W → SBσ and τ SB → Wσ growrapidly with decreasing values of the intensity of the noise for all values of µ . In agreement withwhat is shown in Fig. 8, we have that τ SB → Wσ grows more rapidly than τ W → SBσ for µ ≤ µ crit (andviceversa for µ ≥ µ crit ). As a result, in the weak-noise limit an individual trajectory might betrapped for a very long time in the metastable, non-asymptotic state. We remark that we have notperformed here a systematic evaluation of the exponential relationship between the escape timesand σ , as instead done for µ = 0 . and shown in Fig. 2, also because, as for the reasons explainedbefore, this would require computational resources that are beyond what has been allocated for thisstudy. We remark that this would allow for evaluating for each value of µ the difference betweenthe value of the quasi-potential realised at the attractors and at the M state. The algorithm proposedin [84] can be very useful to reduce the computational burden.We can say that for µ = µ crit our system exhibits a behaviour that is reminiscent of a first-orderphase transition for near-equilibrium statistical mechanical systems, like a liquid-gaseous transi-tion. In our case, µ is the control parameter (corresponding to the temperature in the equilibriumcase), the quasi-potential Φ is the equivalent of a thermodynamic potential, the scaling factor forthe noise intensity σ is the equivalent of (square root of) the temperature, and the globally averagedsurface temperature (cid:104) T S (cid:105) is the natural order parameter, e.g. density. The discontinuous change inthe properties of the system for µ = µ crit is associated to the change in the amount of absorbedand emitted radiation, as a result of the macroscopic change in the albedo of the planet due to thediscontinuity in the position of the ice-line. We remark that choosing a different noise law wouldin general lead to a different value of µ crit , as a result of the fact that the functional form of Φ would be different. E. Relevance of the Choice of the Itˆo Convention for the Noise
A reasonable question to ask concerns to what extent our results are sensitive to the fact thatwe have chosen the It ˆo convention for the noise, which provides the starting point of the results30 ) b)c) d)FIG. 8. Projection of the measure on the variable (cid:104) T S (cid:105) (units of K ) for σ y = 1 . (Panel a), σ y =1 . (Panel b), and σ y = 2 . (Panel c). The spacing of the isolines is the same in the three panels.Panel d): Fraction of the measure supported in the basin of attraction of the SB state as a function of µ andof the noise intensity. As the noise decreases, we observe a fast transition between SB - and W -dominatedpopulations for µ = µ crit ≈ . , which corresponds to a first-order phase transition. discussed in Sect. II B and II C. We argue that choosing other conventions would not alter essen-tially our findings because, to a first approximation, the stochastic forcing we have introduced canbe treated as corresponding to perturbing the system with additive noise of different strengths nearthe cold and W attractors, plus a transition region between the two attractors (which is evidentlyvery sparsely populated by the system), where the effective intensity of the noise decreases withthe globally averaged surface temperature (cid:104) T S (cid:105) and the multiplicative nature of the noise is moreevident.In the phase space region near the cold attractor, we have that − α ( φ, T S ) ∼ . since thetemperature T S is extremely low and the planet is fully glaciated (or almost entirely so), so that α ( φ, T S ) is virtually constant, with α ( φ, T S ) ∼ α min . Near the W attractor, the properties of the31eld α ( φ, T S ) are slightly more complex, because part of the planet is glaciated and part of it isice-free. Nonetheless, to a first approximation, the ratio of the variance of the noise in the SBattractor vs. W attractor is of the order ((1 − α min ) / (1 − α W )) ∼ . Loosely speaking, thecompeting W and SB climate states have different statistical mechanical, microscopic - as well asthermodynamical, macroscopic - temperatures. V. AN ALTERNATIVE CONSTRUCTION OF THE M STATES USING STOCHASTIC PERTUR-BATIONS
As discussed above, the construction of saddles for multistable systems is far from being atrivial task, and requires the use of the edge tracking algorithm introduced in [50, 51] and usedalso by the authors in [20, 45]. We wish to provide here a proof of concept of an alternativeprocedure for constructing the saddles - especially relevant when they are complex M states -without resorting to such an algorithm, but rather using only direct numerical simulations. Theprocedure discussed below might be useful when the edge tracking algorithm is hard to implement.As an example, this could be the case when the presence of similar time scales associated tothe instability along and across the basin boundary might hinder an accurate computation of thesaddles. Alternatively, it can be seen as a way to test the results obtained from the study of thedeterministic dynamics.The idea is to exploit the fact that, as discussed in Section II, under rather general conditionson the noise law, the saddle, in the weak noise limit, acts as the gate for noise-induced transitionsbetween the competing attractors. We then propose to proceed as follows. Let’s consider thefollowing stochastic differential equations: d x ( t ) = F ( x ( t )) dt + σ s k ( x ) d W , k = 1 , . . . , K. (26)We consider the possibility of perturbing the deterministic flow with K different noise laws, de-fined by the K functions s k ( x ) , each leading to a noise with different covariance matrix C ij,k ( x ) .We assume, for simplicity, that the deterministic system defined by ˙ x ( t ) = F ( x ( t )) is bistable.We choose K noise laws such that they obey the hypotheses discussed in Sects. II B and II C.In the weak noise limit σ → , the invariant measure of the k th SDE can be written as Π σ,k ∼ exp (cid:16) − k ( x ) σ (cid:17) . Clearly, for each noise law the quasi-potential Φ k ( x ) is different. Nonetheless, asdiscussed above, in all cases Φ k ( x ) has a local minimum (and is constant) on the support of the32wo deterministic attractors, and it is a saddle with constant value on the saddle separating the twoattractors. We here have to assume that either the saddle is unique, or that all the quasi-potentials Φ k ( x ) select the same saddle as the one with the lowest quasi-potential.Additionally, for each of the K SDEs the drift flow is split differently between the gradient-likecomponent and the rest (see Eq. 6). Therefore, as suggested by Eq. 9, the instantonic paths are alsodifferent; yet, they connect the same attractors to the same saddle. Assuming that the attractorsand the saddles are points or at least small sets given in a coarse-grained description of the phasespace, we can identify the saddle as the only point in space where the instantons corresponding toall the K noise laws will intersect.In order to show that this approach does indeed work, we investigate the properties of K = 3 variants of the Ghil-Sellers diffusive model we studied in [45], differing for the law of the stochas-tic perturbation impacting the energy balance of the climate system. The Ghil-Sellers diffusivemodel can be written by removing the atmosphere-ocean coupling term in Eq. 24. In order to con-form to Eq. 2 (note that we treat below the numerical discretization of a stochastically perturbedpartial different equation), we express the three stochastically perturbed models as follows: ∂T S ( φ, t ) ∂t = µ S ∗ I ( φ ) C ( φ ) (1 − α ( φ, T S )) − O ( T S ) C ( φ ) − D φ [ T S ] C ( φ ) + σs k ( φ, T S ) dWdt , k = 1 , , , (27)where dW is the increment of a one-dimensional Wiener process, and we have: s ( φ, T S ) = µ S ∗ I ( φ ) C ( φ ) , (28) s ( φ, T S ) = 1 , (29) s ( φ, T S ) = µ S ∗ I ( φ ) C ( φ ) , (1 − α ( φ, T S )) . (30)Specifically, we have that s ( φ, T S ) = s ( φ ) and s ( φ, T S ) = s ( φ ) correspond to additive noiselaws, which feature different diffusion matrices. Instead, s ( φ, T S ) is a multiplicative noise law asa result of the temperature-dependence of the albedo and is closely related to what was studied inthe rest of the paper; see Eq. 24. We construct for these three SDEs the invariant measure and,by stochastically averaging, we estimate the instantons connecting the attractors and the saddles.These sets are simple points. We make sure that the instantons are estimated using very weak noiseamplitudes. Results are presented in Fig. 9, where we show the projections on the (∆ T S , (cid:104) T S (cid:105) ) space. Panels a), b), and c) show the invariant measures and the instantons constructed for thenoise laws s (using σ = 0 . ) , s (using σ = 2 . ), and s (using σ = 1 . ), respectively. Note that33 ) b)c) d)FIG. 9. Projected measure in the (∆ T S , (cid:104) T S (cid:105) ) space (units of K on both axes) for selected stochasticallyperturbed simulations of the Ghil-Sellers models and related instantons. All results have been obtainedwith simulations lasting in . × model years. a) Additive noise s : Logarithm of the pdf; W → SB instanton (red dashed line); SB → W instanton (blue dashed line); W attractor (red dot); SB attractor (bluedot); saddle (green dot). The results have been obtained with σ = 0 . . b) Same as in a), for additive noise s (results obtained setting σ = 2 . c) Same as in a), for the multiplicative noise used in the rest of the paper(results obtained setting σ = 1 . . d) The instantons from a) (black line), b) (cyan line) and c) (magentaline) are plotted together. Dashed (dotted) lines correspond to the W → SB ( SB → W ) instantons. Theycross at the saddle (green dot) and at the W attractor (red dot) and at the SB attractor (blue dot). the instanton constructed with the multiplicative noise law looks qualitatively different from whatis show in Fig. 5-6, as a result of the lack of atmospheric motions in the simpler model discussedhere. Panel d) portrays the instantons constructed for the three noise laws. Indeed, we have aconfirmation that all of them are different, as a result of the different noise laws of the three SDEs,and intersect at the attractors and at the saddle. The position of the saddle in the projected phasespace can be identified through this geometric procedure, which is based exclusively on direct34umerical simulations. Considering additional noise laws can be helpful in resolving possiblegeometrical degeneracies due to the use of projections. Projecting in more than two dimensionscould also serve a similar scope and provide a better understanding of the alternative transitionpaths. VI. CONCLUSIONS
The goal of this paper has been the investigation of the properties of the noise-induced transi-tions across the multiple basins of attractions in an intermediate complexity climate model with O (10 ) degrees of freedom, describing the coupled evolution of atmospheric (fast) and oceanic(slow) variables. The model features the co-esistence of W and SB attractors for a fairly broadrange of values of the solar irradiance. In a previous investigation we had been able to constructthe full phase portrait of the deterministic version of the climate model considered here, and hadconstructed, beside the attractors, the M states of the climate system in the region of bistability[20].The stochastic forcing is introduced here as a random modulation for the incoming solar radi-ation, and leads to a nontrivial multiplicative noise law, because the radiative forcing is affectedby the albedo of the surface, which, in turn, depends on the surface temperature. The noise,by allowing transitions between the deterministic basins of attraction, allows for establishing an(apparently) ergodic invariant measure of the system. The theory of stochastic differential equa-tions indicates that for systems obeying the hypoellipticity condition, and for a suitable class ofnoise laws including additive laws and special multiplicative laws like the one considered here,one can write fairly generally the invariant measures in terms of a large deviation law, where therate function can be identified with the quasi-potential. We have clarified how to compute thequasi-potential from the drift and volatility fields of the SDE and explained its property of being aLyapunov function. The quasi-potential has local minima on the deterministic attractors, and has asaddle behaviour at the M states. Additionally, in the weak noise limit, transitions take place alongspecial paths called instantons, which link the deterministic attractors and the M states. While, fora given deterministic dynamics, instantons corresponding to different noise laws follow differentpaths, they all link the same deterministic attractors to the same M states. We have shown howthis property can be exploited to geometrically construct M states directly from direct numericalsimulations of stochastic systems. Indeed, while the edge tracking algorithm applied to the deter-35inistic system is a priori the preferred choice for finding M states, it might become nontrivialto implement in complex numerical models where the intermediate states constructed by bisec-tion might correspond to regions where the model is numerically unstable, possibly because therealised physical fields are extremely exotic or non-realisable.We have studied in detail the noise-induced transitions between the deterministic basins ofattraction in the range of multistability, extending the results presented in a short communication[56]. We have shown that by studying how the average escape time depends on the intensity ofthe noise it is possible to estimate the difference between the value of the quasi-potential at the Mstate and at the attractor that the trajectories are escaping from. The estimates of the instantons areshown to become more precise as we consider simulations where weaker noise is considered. Theinstanton, in a case of special interest where the M state was shown to undergo a symmetry-breakprocess, selects as optimal point of passage between the SB and the W climate the transient Mstate instead of the asymptotic one, possibly as a result of the finite amplitude of the noise.Finally, by studying how the populations of the W and SB climate change as a function of theintensity of the noise, and using the large deviation law for the measure predicted by the theory,we find an estimate of a critical value of µ = µ crit ≈ . such that for µ ≥ µ crit the zero-noiselimit of the invariant measure is supported on the W deterministic attractor, while for µ ≤ µ crit theweak noise limit of the invariant measure is supported on the SB attractor. The asymptotic statecorresponds to the attractor featuring the lowest value of the quasi-potential.These results obtained here indicate that, as soon as noise in some form is added to the system,multistability is factually lost in the weak-noise limit, as the noise law is responsible for selecting,for each value of the control parameter (here µ ), a specific asymptotic state. Changing the valueof the control parameter, one will find one or more abrupt transitions in the statistical properties ofthe system (here realised at µ crit ), i.e., in other terms discontinuities in the response of the systemto changes in the control parameter. What happens in our model at µ = µ crit is mathematicallyanalogous of a first-order phase transition occurring in a near-equilibrium statistical mechanicalsystem. We remark that, since the escape time away from either attractor grows exponentially withthe inverse of the parameter controlling the variance of the noise, an individual trajectory might betrapped for very long time in a metastable state.In collaboration with C. Kilic (Bern) and F. Lunkeit (Hamburg) the authors have started somepreliminary simulations where stochastic forcing is added to PLASIM [97], a much more complexclimate model than the one used in this study (yet missing some essential ocean dynamical pro-36
100 200 300 400 500 600 700 800 900 1000
Time ( G l o b a l M ea n S u r f a c e T e m p e r a t u r e ( K ) Time ( G l o b a l M ea n S ea I c e C o v e r ( ) σ = 1367x0.2 Wm -2 , μ = 1367 Wm -2 FIG. 10. Example of noise-induced transitions between SB and W climate state for present-day solar con-stant ( σ y = 2% ) for the climate model PLASIM [97] in a simulation of y . From top to bottom: zonalaverages of the temperature field and of the sea ice cover; globally averaged surface temperature; fractionof ice-covered ocean. The characteristic escape times is comparable with what is obtained in the simplermodel PUMA-GS used in this work (courtesy of F. Lunkeit). cesses). As shown in Fig. 10, the first results we have obtained are encouraging in indicating thatthe findings of this paper might be relevant for more realistic model configurations. For the future,we aim at obtaining detailed information on the large scale properties of the flow configurationsleading to the noise-induced transitions, taking the thermodynamic lens we originally explored in[9]. An important question of specific relevance for paleoclimatic and planetary science studiesis to understand whether the third climatic state found in [20] for µ = 1 . is recovered also for37ore realistic model configurations. Additionally, the complexity of the dynamical landscape ofthe climate system discussed in [21] suggests the existence of a possibly topologically non-trivialnetwork of transition paths between the many competing attractors, each crossing an M state.Maybe the itinerancy between possibly many competing attractors might be a way to explainingthe ultralow frequency variability of the climate. The computational needs of a naive approach tothese issues seem prohibitive, so that one should definitely take advantage of rare events algorithmsto construct instantonic trajectories [94–96].The approach presented here, based upon combining the knowledge of the dynamical landscapeof a multistable deterministic dynamical system with the analysis of the impacts of stochasticperturbations, seems of more general interested than the specific problem we have studied. Alongthese lines, in Appendix A we speculate on the the possible relevance of M states in the contextof the theory of biological evolution and of synthetic models of evolution. Specifically, the ideais that the presence of qualitatively diverse historical paths of evolution might result from the thepresence of M states in suitable defined dynamical landscapes, exactly because M states hinderpredictability of the second kind in the sense of Lorenz.In the specific case of geosciences, we believe that it can be key for addressing the challenge ofunderstanding tipping points in the Earth system [69] as well as providing insights to a large classof multistable systems [68]. We will dedicate future efforts exactly to exploring these researchlines, in particularly looking at the Atlantic meridional overturning circulation, an element of theglobal ocean circulation that is well-know to have more than one competing modes of operation .The occurrence of one or of the other state has important implications on the global climate, andrather dramatic ones for the regional climate of the North Atlantic sector; see [98] for a review ofthis topic. Multistability has been recently reported in the von Karman turbulent flow in [99]: themethods proposed in this paper could elucidate paths and mechanisms underlying the transitionsbetween the asymptotic states. ACKNOWLEDGMENTS
VL and TB acknowledge the support received by the EU Horizon2020 projects Blue-Action(grant No. 727852) and CRESCENDO (grant No. 641816). VL acknowledges the support theEU Horizon2020 project TiPES (grant No. 820970), of the DFG SFB/Transregio project TRR181,and of the Royal Society (grant No. IEC/R2/170001). TB acknowledges the support received38rom the Institute for Basic Science (IBS), Republic of Korea, under IBS-R028-D1. The authorswish to thank T. T´el for inspiring exchanges. VL wishes to thank S. Brusatte for suggesting somefascinating literature on evolutionary biology; M. Ghil and J. Yorke for great encouragement in thecourse of several cheerful discussions; T. Kuna, T. Lelievre, X.-M. Li, G. Pavliotis, and N. Zaglifor providing insights on stochastic dynamics. This paper has been partly prepared during theAdvanced Workshop on
Nonequilibrium Systems in Physics, Geosciences, and Life Sciences heldat ICTP, Trieste, on May 15-24 2018 and during the Scientific Programme on the
Mathematics ofClimate and the Environment held at the Institut H. Poincar´e, Paris on September-December 2019.The authors wish to thank F. Lunkeit for kindly sharing Fig. 10. The authors wish to dedicate thispaper to the memory of Bruno Eckhardt (1960-2019).
Appendix A: An Interdisciplinary Outlook: Evolutionary Biology
The concepts discussed in this paper might provide a conceptual and mathematical frameworkof possible interest for thinking at evolutive processes in biology and for constructing syntheticmodels of evolution. In the speculative discussion below, specifically, we want to highlight therole of M states in making possible the existence of multiple possible yet vastly different paths ofhistorical development of biological systems.In the influential book discussing the Cambrian fossils from the Burgess shale and their impor-tance for explaining the mechanisms of evolution, Gould [67] presents some ideas on the scientificmethodology inherent to historical sciences, and, specifically, to evolutionary biology. He clarifiesthat historical scientific explanations take the form of a narrative, whereby subsequent phenomenafollow in a specific order : the time-ordering and causal links between such phenomena can bediscovered and convincingly explained when multiple, independent sources provide indication forthe same historical pattern of change. Gould argues that this method is fundamentally differentfrom the classic - which he calls stereotypical - scientific method ´a la Galileo whereby one specificexperiment can be repeated in idealised conditions, and the final outcomes of the experiment canbe predicted using the fundamental laws of nature. The latter is the - indeed too simplified - ver-sion he gives of how hard sciences work. An obvious difference between the two methods comesfrom the fact that certain experiments - like the evolution of the climate and of the biosphere of aspecific planet - cannot be repeated. Another difference comes from the fact - Gould argues - that The related storyline approach is being currently proposed for studying weather and climate phenomena [100]. contingency : while we observe a specific path of historicaldevelopment, many others of such paths are indeed compatible with the laws of nature, and couldhave been realised had the system been forced in a slightly different way in the past. We are ableto understand the mechanisms of historical development, but not to predict accurately the specificpath. Gould argues that the evolution in our planet could have led to fundamentally different formsof life, had the contingencies been different in the distant past. The existence of our own speciesis the - a priori very unlikely - result of such contingencies. In general, in systems dominatedby contingency, if we could run again the movie the outcome would be vastly different . Gould’sviews on evolution have been criticised by authors, as Conway Morris, proposing that conver-gence , rather than contingency, is the main mechanism of evolution, so that evolution is seen as a -mostly - deterministic path of change, of which nowadays we see the unavoidable outcome [101];see the debate in [102]. Some authors have proposed that both mechanisms are indeed in action,yet dominant at different scales of diversification of the organisms [103].We argue that Gould’s view can be put in the context of the mathematical framework discussedin this paper. Stochastically forced complex systems evolve in a phase space where an individualtrajectory corresponds to a historical realisation of the system. The historical realisations, evenstarting from the same initial condition, can be vastly - even qualitatively - different if the de-terministic dynamics supports the existence of multiple attractors, because different realisationscan be trapped for very long times in very distant regions of the phase space. As discussed here,stochastic forcing allows for the system to jump across the various basins of attraction, and themechanism defining the evolution of the trajectory are defined by the differential equations. Thepredictability of the system, both of the first and second kind in the sense of Lorenz, is finite butnon-zero. Finally, we can interpret the M states as the true agents of the contingency discussed byGould. It is not the presence of a forcing, however strong, that changes radically the future historyof the system, but rather the existence of special regions of the phase space - near the M states- where even small perturbations can force two nearby trajectories towards qualitatively differentfuture histories. What Conway Morris proposes is associated with a scenario where M states areeither absent - so that the system is not multistable - or the noise is so weak that the probability ofgetting close to an M state is exceedingly low, despite the presence of stochastic perturbations, thesystem will be (almost) always quite close to a specific deterministic attractor, in correspondencewith which basin of attraction the initial condition belongs to. In other terms. even if we could runagain the movie , i.e., run another simulation, the outcome would be very similar.40he interpretation given above receives some supports from recent results obtained on syntheticmodels of evolution. Using the the Tangled Nature model, which is conjectured to be a prototyp-ical example of evolution, evolutive processes are interpreted as orbits of stochastic systems ina complex dynamical landscape featuring two or more competing metastable states [104, 105].This viewpoint, which has been openly inspired by Gould’s ideas, is in close correspondence withwhat has been discussed in this paper. While in these works the language and methodology areeminently of statistical mechanical nature and they aim at detecting and studying the metastablesstates, the viewpoint proposed in this paper has a stronger emphasis on understanding the transi-tions paths between such competing attractors through the M states. [1] M.I. Budyko. The effect of solar radiation variations on the climate of the earth.
Tellus , 21:611–619,1969.[2] WD Sellers. A global climatic model based on the energy balance of the earthatmosphere system.
J.Appl. Meteorol. , 8:392–400, 1969.[3] M. Ghil. Climate stability for a Sellers-type model.
J. Atmos. Sci. , 33:3–20, 1976.[4] J. L. Kirschvink. Late Proterozoic low-latitude global glaciation: The snowball Earth. In J. W. Schopfand C. Klein, editors,
The Proterozoic Biosphere: A Multidisciplinary Study , chapter 8, pages 91–92.Cambridge University Press, 1992.[5] Paul F. Hoffman and Daniel P. Schrag. The snowball Earth hypothesis: testing the limits of globalchange.
Terra Nova , 14(3):129–155, 2002.[6] R.T. Pierrehumbert, D.S. Abbot, Aiko Voigt, and D. Koll. Climate of the Neoproterozoic.
AnnualReview of Earth and Planetary Sciences , 39(1):417–460, 2011.[7] William T. Hyde, Thomas J. Crowley, Steven K. Baum, and W. Richard Peltier. Neoproterozoicsnowball earth simulations with a coupled climate/ice-sheet model.
Nature , 405:425–429, 2000.[8] Aiko Voigt and Jochem Marotzke. The transition from the present-day climate to a modern SnowballEarth.
Climate Dynamics , 35(5):887–905, 2010.[9] Valerio Lucarini, Klaus Fraedrich, and Frank Lunkeit. Thermodynamic analysis of snowball earthhysteresis experiment: Efficiency, entropy production and irreversibility.
Quarterly Journal of theRoyal Meteorological Society , 136(646):2–11, 2010.[10] R. Boschi, V. Lucarini, and S. Pascale. Bistability of the climate around the habitable zone: a hermodynamic investigation. Icarus , 227:1724–1742, 2013.[11] Thomas J. Crowley, William T. Hyde, and W. Richard Peltier. Co2 levels required for deglaciationof a near-snowball earth.
Geophysical Research Letters , 28(2):283–286, 2001.[12] Manuel Linsenmeier, Salvatore Pascale, and Valerio Lucarini. Climate of earth-like planets with highobliquity and eccentric orbits: Implications for habitability conditions.
Planetary and Space Science ,105:43 – 59, 2015.[13] C. Kilic, C. C. Raible, and T. F. Stocker. Multiple climate states of habitable exoplanets: The role ofobliquity and irradiance.
The Astrophysical Journal , 844(2):147, 2017.[14] Cevahir Kilic, Frank Lunkeit, Christoph C. Raible, and Thomas F. Stocker. Stable equatorial ice beltsat high obliquity in a coupled atmosphere–ocean model.
The Astrophysical Journal , 864(2):106, sep2018.[15] V Lucarini, S Pascale, R Boschi, E Kirk, and N Iro. Habitability and multistablility in earth-likeplantets.
Astr. Nach. , 334(6):576–588, 2013.[16] Dorian S. Abbot, Jonah Bloch-Johnson, Jade Checlair, Navah X. Farahat, R. J. Graham, DavidPlotkin, Predrag Popovic, and Francisco Spaulding-Astudillo. Decrease in hysteresis of planetaryclimate for planets with long solar days.
The Astrophysical Journal , 854(1):3, feb 2018.[17] Jade Checlair, Kristen Menou, and Dorian S. Abbot. No snowball on habitable tidally locked planets.
The Astrophysical Journal , 845(2):132, 2017.[18] J. P. Lewis, A. J. Weaver, and M. Eby. Snowball versus slushball earth: Dynamic versus nondynamicsea ice?
Journal of Geophysical Research: Oceans , 112(C11):C11014, 2007.[19] Dorian S. Abbot, Aiko Voigt, and Daniel Koll. The jormungand global climate state and implicationsfor neoproterozoic glaciations.
Journal of Geophysical Research: Atmospheres , 116(D18), 2011.[20] Valerio Lucarini and Tam´as B´odai. Edge states in the climate system: exploring global instabilitiesand critical transitions.
Nonlinearity , 30(7):R32, 2017.[21] M. Brunetti, J. Kasparian, and C. V´erard. Co-existing climate attractors in a coupled aquaplanet.
Climate Dynamics , 53(9):6293–6308, 2019.[22] Illeana G´omez-Leal, Lisa Kaltenegger, Valerio Lucarini, and Frank Lunkeit. Climate sensitivity tocarbon dioxide and the moist greenhouse threshold of earth-like planets under an increasing solarforcing.
The Astrophysical Journal , 869(2):129, dec 2018.[23] Illeana G´omez-Leal, Lisa Kaltenegger, Valerio Lucarini, and Frank Lunkeit. Climate sensitivity toozone and its relevance on the habitability of earth-like planets.
Icarus , 321:608 – 618, 2019.
24] James F. Kasting, Daniel P. Whitmire, and Ray T. Reynolds. Habitable zones around main sequencestars.
Icarus , 101(1):108 – 128, 1993.[25] Viviane Baladi.
Positive Transfer Operators and Decay of Correlations . World Scientific, Singapore,2000.[26] Mark Pollicott. On the rate of mixing of Axiom A flows.
Inventiones Mathematicae , 81(3):413–426,October 1985.[27] David Ruelle. Resonances of Chaotic Dynamical Systems.
Physical review letters , 56(5):405–407,1986.[28] Alexis Tantet, Valerio Lucarini, Frank Lunkeit, and Henk A Dijkstra. Crisis of the chaotic attractorof a climate model: a transfer operator approach.
Nonlinearity , 31(5):2221, 2018.[29] Masatoshi Shiino. Dynamical behavior of stochastic systems of infinitely many coupled nonlinearoscillators exhibiting phase transitions of mean-field type: H theorem on asymptotic approach toequilibrium and critical slowing down of order-parameter fluctuations.
Phys. Rev. A , 36:2393–2412,Sep 1987.[30] G.A. Pavliotis.
Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck andLangevin Equations . Texts in Applied Mathematics. Springer New York, 2014.[31] F. Ragone, V. Lucarini, and F. Lunkeit. A new framework for climate sensitivity and prediction: amodelling perspective.
Climate Dynamics , 46(5):1459–1471, March 2016.[32] V. Lucarini, F. Ragone, and F. Lunkeit. Predicting climate change using response theory: Globalaverages and spatial patterns.
Journal of Statistical Physics , 166(3):1036–1064, February 2017.[33] Valerio Lembo, Valerio Lucarini, and Francesco Ragone. Beyond Forcing Scenarios: PredictingClimate Change through Response Operators in a Coupled General Circulation Model. arXiv e-prints , page arXiv:1912.03996, Dec 2019.[34] David Ruelle. A review of linear response theory for general differentiable dynamical systems.
Nonlinearity , 22:855–870, 2009.[35] Michael Ghil, Micka¨el David Chekroun, and Eric Simonnet. Climate dynamics and fluid mechanics:Natural variability and related uncertainties.
Physica D: Nonlinear Phenomena , 237(14-17):2111–2126, August 2008.[36] Micka¨el David Chekroun, Eric Simonnet, and Michael Ghil. Stochastic climate dynamics: Randomattractors and time-dependent invariant measures.
Physica D: Nonlinear Phenomena , 240(21):1685–1700, October 2011.
37] A. N. Carvalho, J. Langa, and J. C. Robinson. The pullback attractor. In
Attractors for infinite-dimensional non-autonomous dynamical systems , volume 182 of
Applied Mathematical Sciences ,pages 3–22. Springer New York, 2013.[38] F. J. Romeiras, C. Grebogi, and E. Ott. Multifractal properties of snapshot attractors of random maps.
Phys. Rev. A , 41(2):784, 1990.[39] G. Dr´otos, T. B´odai, and T. T´el. Probabilistic concepts in a changing climate: A snapshot attractorpicture.
Journal of Climate , 28(8):3275–3288, 2015.[40] Amir Dembo and Jean-Dominique Deuschel. Markovian perturbation, response and fluctuation dis-sipation theorem.
Ann. Inst. H. Poincar Probab. Statist. , 46(3):822–852, 08 2010.[41] Roland Assaraf, Benjamin Jourdain, Tony Leli`evre, and Rapha¨el Roux. Computation of sensitivitiesfor the invariant measure of a parameter dependent diffusion.
Stochastics and Partial DifferentialEquations: Analysis and Computations , 6(2):125–183, Jun 2018.[42] Valerio Lucarini. Response operators for markov processes in a finite state space: Radius of conver-gence and link to the response theory for axiom a systems.
Journal of Statistical Physics , 162(2):312–333, Jan 2016.[43] Micka¨el David Chekroun, J. David Neelin, Dmitri Kondrashov, J. C. McWilliams, and Michael Ghil.Rough parameter dependence in climate models and the role of Ruelle-Pollicott resonances.
Pro-ceedings of the National Academy of Sciences of the United States of America , 111(5):1684–1690,March 2014.[44] Michael Ghil and Valerio Lucarini. The Physics of Climate Variability and Climate Change. arXive-prints , page arXiv:1910.00583, Oct 2019.[45] Tam´as B´odai, Valerio Lucarini, Frank Lunkeit, and Robert Boschi. Global instability in the ghil–sellers model.
Clim. Dyn. , 44(11):3361–3381, 2014.[46] Celso Grebogi, Edward Ott, and James A. Yorke. Fractal basin boundaries, long-lived chaotic tran-sients, and unstable-unstable pair bifurcation.
Phys. Rev. Lett. , 50:935–938, Mar 1983.[47] Carl Robert, Kathleen T. Alligood, Edward Ott, and James A. Yorke. Explosions of chaotic sets.
Physica D: Nonlinear Phenomena , 144(1):44 – 61, 2000.[48] E. Ott.
Chaos in Dynamical Systems . Cambridge University Press, 2002.[49] Y.-C. Lai and T. T´el.
Transient Chaos . Springer, New York, 2011.[50] Joseph D. Skufca, James A. Yorke, and Bruno Eckhardt. Edge of chaos in a parallel shear flow.
Phys.Rev. Lett. , 96:174101, May 2006.
51] Tobias M. Schneider, Bruno Eckhardt, and James A. Yorke. Turbulence transition and the edge ofchaos in pipe flow.
Phys. Rev. Lett. , 99:034502, Jul 2007.[52] D. A. W. Barton and J. Sieber. Systematic experimental exploration of bifurcations with noninvasivecontrol.
Phys. Rev. E , 87:052916, May 2013.[53] Jan Sieber, Oleh E. Omel’chenko, and Matthias Wolfrum. Controlling unstable chaos: Stabilizingchimera states by feedback.
Phys. Rev. Lett. , 112:054102, Feb 2014.[54] Daniel M. Abrams and Steven H. Strogatz. Chimera states for coupled oscillators.
Phys. Rev. Lett. ,93:174102, Oct 2004.[55] O. E. Omel’chenko. The mathematics behind chimera states.
Nonlinearity , 31(5):R121, 2018.[56] Valerio Lucarini and Tam´as B´odai. Transitions across melancholia states in a climate model: Recon-ciling the deterministic and stochastic points of view.
Phys. Rev. Lett. , 122:158701, Apr 2019.[57] Denis R. Bell.
Stochastic Differential Equations and Hypoelliptic Operators , pages 9–42. Birkh¨auserBoston, Boston, MA, 2004.[58] Peter Hanggi. Escape from a metastable state.
Journal of Statistical Physics , 42(1):105–148, Jan1986.[59] R.L. Kautz. Activation energy for thermally induced escape from a basin of attraction.
PhysicsLetters A , 125(6):315 – 319, 1987.[60] P Grassberger. Noise-induced escape from attractors.
Journal of Physics A: Mathematical andGeneral , 22(16):3283, 1989.[61] M. I. Freidlin and A.D. Wentzell.
Random Perturbations of Dynamical Systems . Springer, New York,1984.[62] R. Graham, A. Hamm, and T. T´el. Nonequilibrium potentials for dynamical systems with fractalattractors or repellers.
Phys. Rev. Lett. , 66:3089–3092, Jun 1991.[63] A. Hamm, T. T´el, and R. Graham. Noise-induced attractor explosions near tangent bifurcations.
Physics Letters A , 185(3):313 – 320, 1994.[64] Suso Kraut and Ulrike Feudel. Multistability, noise, and attractor hopping: The crucial role of chaoticsaddles.
Phys. Rev. E , 66:015207, Jul 2002.[65] S. Beri, R. Mannella, D. G. Luchinsky, A. N. Silchenko, and P. V. E. McClintock. Solution of theboundary value problem for optimal escape in continuous stochastic systems and maps.
Phys. Rev.E , 72:036131, Sep 2005.[66] Freddy Bouchet and Julien Reygner. Generalisation of the eyring–kramers transition rate formula to rreversible diffusion processes. Annales Henri Poincar´e , 17(12):3499–3532, Dec 2016.[67] S. J. Gould.
Wonderful Life: The Burgess shale and the Nature of History . W.W. Norton, New York,1989.[68] Ulrike Feudel, Alexander N. Pisarchik, and Kenneth Showalter. Multistability and tipping: Frommathematics and physics to climate and brainminireview and preface to the focus issue.
Chaos: AnInterdisciplinary Journal of Nonlinear Science , 28(3):033501, 2018.[69] Timothy M. Lenton, Hermann Held, Elmar Kriegler, Jim W. Hall, Wolfgang Lucht, Stefan Rahm-storf, and Hans Joachim Schellnhuber. Tipping elements in the earth’s climate system.
Proceedingsof the National Academy of Sciences , 105(6):1786–1793, 2008.[70] Juergen Vollmer, Tobias M Schneider, and Bruno Eckhardt. Basin boundary, edge of chaos and edgestate in a two-dimensional model.
New Journal of Physics , 11(1):013040, 2009.[71] Tam´as B´odai and Valerio Lucarini. Rough basin boundaries in high dimension: Can we classify themexperimentally. arXiv e-prints , page arXiv:XX.YY, Jan 2020.[72] E. N. Lorenz. The physical bases of climate and climate modelling. climate predictability. In
GARPPublication Series , pages 132–136. WMO, 1975.[73] Martin Hairer. On malliavin’s proof of h omander’s theorem.
Bulletin des Sciences Math´ematiques ,135(6):650 – 666, 2011. Special issue in memory of Paul Malliavin.[74] Pierre Gaspard. Trace formula for noisy flows.
Journal of Statistical Physics , 106(1):57–96, Jan2002.[75] Freddy Bouchet, Krzysztof Gawedzki, and Cesare Nardini. Perturbative calculation of quasi-potential in non-equilibrium diffusions: A mean-field example.
Journal of Statistical Physics ,163(5):1157–1210, Jun 2016.[76] P Ao. Potential in stochastic differential equations: novel construction.
Journal of Physics A: Math-ematical and General , 37(3):L25–L30, jan 2004.[77] L Yin and P Ao. Existence and construction of dynamical potential in nonequilibrium processeswithout detailed balance.
Journal of Physics A: Mathematical and General , 39(27):8593–8601, jun2006.[78] Joseph Xu Zhou, M. D. S. Aliyu, Erik Aurell, and Sui Huang. Quasi-potential landscape in complexmulti-stable systems.
Journal of The Royal Society Interface , 9(77):3539–3553, 2012.[79] R. D. Brackston, A. Wynn, and M. P. H. Stumpf. Construction of quasipotentials for stochasticdynamical systems: An optimization approach.
Phys. Rev. E , 98:022136, Aug 2018.
80] Ying Tang, Ruoshi Yuan, Gaowei Wang, Xiaomei Zhu, and Ping Ao. Potential landscape of highdimensional nonlinear stochastic dynamics with large noise.
Scientific Reports , 7(1):15762, 2017.[81] R. Graham. Macroscopic potentials, bifurcations and noise in dissipative systems. In L. Garrido, ed-itor,
Fluctuations and Stochastic Phenomena in Condensed Matter , pages 1–34, Berlin, Heidelberg,1987. Springer Berlin Heidelberg.[82] R. Graham and T. T´el. Nonequilibrium potential for coexisting attractors.
Phys. Rev. A , 33:1322–1337, Feb 1986.[83] Anton Bovier, Michael Eckhoff, Vronique Gayrard, and Markus Klein. Metastability in reversiblediffusion processes i. sharp asymptotics for capacities and exit times.
Journal of the European Math-ematical Society , 6:399–424, 10 2004.[84] Tam´as B´odai. An efficient algorithm to estimate the potential barrier height from noise-inducedescape time data. arXiv e-prints , page arXiv:1808.06903, Aug 2018.[85] T. Leli`evre. Accelerated dynamics: Mathematical foundations and algorithmic improvements.
TheEuropean Physical Journal Special Topics , 224(12):2429–2444, Sep 2015.[86] Giacomo Di Ges`u, Tony Leli`evre, Dorian Le Peutrec, and Boris Nectoux. Sharp asymptotics of thefirst exit point density.
Annals of PDE , 5(1):5, Mar 2019.[87] Valerio Lucarini. Stochastic resonance for nonequilibrium systems.
Phys. Rev. E , 100:062124, Dec2019.[88] Thomas Frisius, Frank Lunkeit, Klaus Fraedrich, and Ian N. James. Storm-track organization andvariability in a simplified atmospheric global circulation model.
Quarterly Journal of the RoyalMeteorological Society , 124(548):1019–1043, 1998.[89] Barry Saltzman.
Dynamical Paleoclimatology . Academic Press, New York, 2001.[90] Jos P Peixoto and Abraham H Oort.
Physics of Climate . American Institute of Physics, 1992.[91] Valerio Lucarini, Richard Blender, Corentin Herbert, Francesco Ragone, Salvatore Pascale, andJeroen Wouters. Mathematical and Physical Ideas for Climate Science.
Reviews of Geophysics ,52(4):1 – 51, 2014.[92] James R. Holton.
An introduction to dynamic meteorology . International Geophysics Series. ElsevierAcademic Press,, Burlington, MA, 4 edition, 2004.[93] J. Zinn-Justin.
Quantum Field Theory and Critical Phenomena . Oxford University Press, Oxford,1996.[94] G. Rubino and B. Tuffin.
Rare Event Simulation using Monte Carlo Methods . Wiley, 2009.
95] F. Ragone, J. Wouters, and F. Bouchet. Computation of extreme heat waves in climate models usinga large deviation algorithm.
Proceedings of the National Academy of Sciences , 115(1):24–29, 2018.[96] Tobias Grafke, Rainer Grauer, and Tobias Schfer. Instanton filtering for the stochastic burgers equa-tion.
Journal of Physics A: Mathematical and Theoretical , 46(6):062002, jan 2013.[97] Klaus Fraedrich, Heiko Jansen, Edilbert Kirk, Ute Luksch, and Frank Lunkeit. The Planet Simulator:Towards a user friendly model.
Meteorologische Zeitschrift , 14(3):299–304, 2005.[98] T Kuhlbrodt, A Griesel, M Montoya, A Levermann, M Hofmann, and S Rahmstorf. On the drivingprocesses of the Atlantic meridional overturning circulation.
Atlantic , 45(2004):RG2001, 2007.[99] D. Faranda, Y. Sato, B. Saint-Michel, C. Wiertel, V. Padilla, B. Dubrulle, and F. Daviaud. Stochasticchaos in a turbulent swirling flow.
Phys. Rev. Lett. , 119:014502, Jul 2017.[100] Theodore G. Shepherd, Emily Boyd, Raphael A. Calel, Sandra C. Chapman, Suraje Dessai, Ioana M.Dima-West, Hayley J. Fowler, Rachel James, Douglas Maraun, Olivia Martius, Catherine A. Senior,Adam H. Sobel, David A. Stainforth, Simon F. B. Tett, Kevin E. Trenberth, Bart J. J. M. van den Hurk,Nicholas W. Watkins, Robert L. Wilby, and Dimitri A. Zenghelis. Storylines: an alternative approachto representing uncertainty in physical aspects of climate change.
Climatic Change , 151(3):555–571,Dec 2018.[101] S. Conway Morris.
The crucible of creation: the Burgess Shale and the rise of animals . OxfordUniversity Press, Oxford, 1998.[102] Simon Conway Morris and Stephen Jay Gould. Showdown on the burgess shale.
Natural Historymagazine , 107(10):48–55, 1998.[103] J. B. Losos.
Improbable destinies: Fate, chance, and the future of evolution . Riverhead Books, NewYork, 2017.[104] Dominic Jones, Henrik Jeldtoft Jensen, and Paolo Sibani. Tempo and mode of evolution in thetangled nature model.
Phys. Rev. E , 82:036121, Sep 2010.[105] Henrik Jeldtoft Jensen. Tangled nature: a model of emergent structure and temporal mode amongco-evolving agents.
European Journal of Physics , 40(1):014005, dec 2018., 40(1):014005, dec 2018.