Global Stability Properties of the Climate: Melancholia States, Invariant Measures, and Phase Transitions

TL;DR

This paper investigates the stability of Earth's climate states, focusing on transitions between warm and snowball conditions, using stochastic modeling, invariant measures, and Melancholia states to understand phase transitions and noise-induced shifts.

Contribution

It introduces a novel method for constructing Melancholia states from simulations and analyzes their role as gateways in noise-induced climate transitions.

Findings

Melancholia states sit at the boundary between competing climates.

Noise induces transitions between warm and snowball states.

Phase transition occurs at a critical solar irradiance level.

Abstract

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 snowball climate features global glaciation and conditions that can hardly support life forms. Paleoclimatic evidences suggest that in the past our planet flipped between these two states. The main physical mechanism responsible for such an instability is the ice-albedo feedback. In a previous work, we defined the Melancholia states that sit between the two climates. Such states are embedded in the boundaries between the two basins of attraction and feature extensive glaciation down to relatively low latitudes. Here, we explore the global stability properties of the system by introducing random perturbations as modulations to the intensity of the incoming 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 the Melancholia states are the gateways for the noise-induced transitions. In the region of multistability, in the zero-noise limit, the measure is supported only on one of the competing attractors. For low (high) values of the solar irradiance, the limit measure is the snowball (warm) climate. The changeover between the two regimes corresponds to a first-order phase transition in the system. The framework we propose seems of general relevance for the study of complex multistable systems. Finally, we put forward a new method for constructing Melancholia states from direct numerical simulations, black which provides a possible alternative with respect to the edge-tracking algorithm.

Figures

Equations

$${\dot{x}}_i=F_i(\mathbf{x}),i=1,\dots ,N,$$

$$d{x}_i=F_i(\mathbf{x})dt+\sigma s(\mathbf{x}{)}_{ij}d{W}_j,$$

$${\Pi }_{\sigma }(\mathbf{x})\sim \exp \left(-\frac{2\Phi (\mathbf{x})}{{\sigma }^2}\right),$$

$${\partial }_j{J}_j(\mathbf{x})=0,J_j(\mathbf{x})=-F_j(\mathbf{x}){\Pi }_{\sigma }(\mathbf{x})+{\sigma }^2{\partial }_i\left(C_{ij}(\mathbf{x}){\Pi }_{\sigma }(\mathbf{x})\right),$$

$$F_i(\mathbf{x}){\partial }_i\Phi (\mathbf{x})+C_{ij}(\mathbf{x}){\partial }_i\Phi (\mathbf{x}){\partial }_j\Phi (\mathbf{x})=0.$$

$$F_i(\mathbf{x})=R_i(\mathbf{x})-C_{ij}(\mathbf{x}){\partial }_j\Phi (\mathbf{x})$$

$$d\Phi (\mathbf{x})/dt=-C_{ij}(\mathbf{x}){\partial }_i\Phi (\mathbf{x}){\partial }_j\Phi (\mathbf{x})+R_i(\mathbf{x}){\partial }_i\Phi (\mathbf{x})=-C_{ij}(\mathbf{x}){\partial }_i\Phi (\mathbf{x}){\partial }_j\Phi (\mathbf{x}).$$

$$P_{j,\sigma }(p)=\frac{1}{{\tau }_{j,\sigma }}\exp \left(-p/{\tau }_{j,\sigma }\right),{\tau }_{j,\sigma }\sim \exp \left(\frac{2\Delta {\Phi }_j}{{\sigma }^2}\right),$$

$$d{x}_i/dt={\tilde{F}}_i(\mathbf{x})=R_i(\mathbf{x})+C_{ij}(\mathbf{x}){\partial }_j\Phi (\mathbf{x})$$

$${\tau }_{\sigma }^{1\to 2}\propto \exp \left(\frac{2(\Phi ({\Pi }_1)-\Phi ({\Omega }_1)}{{\sigma }^2}\right),$$

$${\tau }_{\sigma }^{2\to 1}\propto \exp \left(\frac{2(\Phi ({\Pi }_1)-\Phi ({\Omega }_2)}{{\sigma }^2}\right),$$

$$\frac{{\tau }_{\sigma }^{1\to 2}}{{\tau }_{\sigma }^{2\to 1}}\propto \exp \left(\frac{2(\Phi ({\Omega }_2)-\Phi ({\Omega }_1)}{{\sigma }^2}\right).$$

$${\dot{P}}_{1,\sigma }$$

$${\dot{P}}_{2,\sigma }$$

$$\frac{P_{1,\sigma }}{P_{2,\sigma }}=\frac{{\tau }_{\sigma }^{1\to 2}}{{\tau }_{\sigma }^{2\to 1}}\propto \exp \left(\frac{2(\Phi ({\Omega }_2)-\Phi ({\Omega }_1)}{{\sigma }^2}\right).$$

$${\partial }_t\zeta$$

$${\partial }_{t}D$$

$${\partial }_t{T}_a^{\prime }$$

$${\partial }_t\ln p_s$$

$${\partial }_{\ln \sigma }\Psi$$

$$T_R=(T_R{)}_{tp}+\sqrt{[L(z_{tp}-z(\sigma ))/2{]}^2+S^2}+L(z_{tp}-z(\sigma ))/2,$$

$$(T_R{)}_{tp}=⟨T_S⟩-\bar{L}{z}_{tp},$$

$$L(\lambda ,\varphi )={\partial }_z{T}_R=(T_S(\lambda ,\varphi )-(T_R{)}_{tp})/z_{tp},$$

$${\partial }_t{T}_S(t,\varphi ,\lambda )=\mu \frac{I(\varphi )}{C(\varphi )}\frac{S_0^{\ast }}{4}(1-\alpha (\varphi ,T_S))-\frac{O(T_S)}{C(\varphi )}-\frac{{\bar{D}}_{\varphi }[T_S]}{C(\varphi )}+\frac{\chi [T_S,T_A]}{C(\varphi )}+s.f.,$$

$$s.f.=\sigma s(T_s,\varphi ,\lambda )\frac{dW}{dt}=\sigma \mu \frac{I(\varphi )}{C(\varphi )}\frac{S_0^{\ast }}{4}(1-\alpha (\varphi ,T_S))\frac{dW}{dt},$$

$$d\mathbf{x}(t)=\mathbf{F}(\mathbf{x}(\mathbf{t}))dt+\sigma {\mathbf{s}}_{\mathbf{k}}(\mathbf{x})d\mathbf{W},k=1,\dots ,K.$$

$$\frac{\partial T_S(\varphi ,t)}{\partial t}=\mu \frac{S_0^{\ast }}{4}\frac{I(\varphi )}{C(\varphi )}(1-\alpha (\varphi ,T_S))-\frac{O(T_S)}{C(\varphi )}-\frac{D_{\varphi }[T_S]}{C(\varphi )}+\sigma s_k(\varphi ,T_S)\frac{dW}{dt},k=1,2,3,$$

$$s_1(\varphi ,T_S)$$

$$s_2(\varphi ,T_S)$$

$$s_3(\varphi ,T_S)$$

Full text

MSC Classification: 76U05 Rotating fluids; 82C26 Dynamic and nonequilibrium phase transitions (general); 60Gxx Stochastic processes; 37D45 Strange attractors, chaotic dynamics; 85A20 Planetary atmospheres; 76E20 Stability and instability of geophysical and astrophysical flows

I Introduction

In the late 1960s and in the 1970s, Budyko, Selllers, and Ghil Budyko1969 ; Sellers1969 ; Ghil1976 proposed the idea that the Earth, in the current astrophysical and astronomical configuration, supports two co-existing attractors, the warm (W) state we live in, and the so-called snowball (SB) state, which is characterised by global glaciation and globally averaged surface temperature of about 200-220 K. Using parsimonious yet physically meaningful energy balance models, they indicated that the bistability of 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 Boltzmann radiative feedback (a warmer surface emits more radiation to space). The relevance of the theoretical ansatz became apparent when paleoclimatic data showed that, indeed, our planet has been flipping in and out of states of global glaciations corresponding to the predicated SB states during the Proterozoic, about 650 Mya Kirschvink ; Hoffman2002 ; Pierrehumbert2011a . W According to these energy balance models, the Earth’s climate is bistable for a substantial range of values of the solar irradiance $S^{*}$, which include the present day value. Below the critical value $S^{*}_{W\rightarrow SB}$, only the SB state is permitted, whereas above the critical value $S^{*}_{SB\rightarrow W}$, only the W state is permitted. Such critical values, which determine the boundaries of the region in parametric space where bistability is realised, are defined by bifurcations that occur when, roughly speaking, the strength of the positive, destabilising feedbacks becomes as strong as the negative, stabilizing feedbacks. Indeed, models of different levels of complexity ranging up to the state-of-the-art Earth System Models currently used for climate projections agree on predicting the existence of multistability in the climate system and point to the fundamental mechanisms described above as responsible for it, as well as providing values for $S^{*}_{W\rightarrow SB}$ that are in broad agreement with those obtained using simple models Hyde2000 ; Voigt2010 ; Lucarini2010a . We remark that both the concentration of greenhouse gases and the position of the continents have an impact on the values of $S^{*}_{W\rightarrow SB}$ and $S^{*}_{SB\rightarrow W}$ and on the properties of the W and SB states Boschi2013 . Extremely high values of the concentration of CO2 seem to be needed to deglaciate from a SB state crowley2001 .

Improving our understanding of the critical transitions associated to such a bistability is a key challenge of geosciences and has strong implications also in terms of the quest for understanding or anticipating planetary habitability. Planets in the habitable zone have astronomical and astrophysical configurations that allow, in principle, the presence of water at surface. Therefore, $S^{*}_{W\rightarrow SB}$ defines the cold boundary of the habitable zone. Clearly, an exoplanet in the habitable zone can be in the regime of bistability: if the planet is in the SB state, it will have very hard time supporting life111The project EDEN (http://project-eden.space/) combines ideas and methods in astrophysics, planetary sciences, and astrobiology to search for and characterize nearby habitable worlds.. Additionally, astronomical parameters such as the obliquity of the planet Linsenmeier2015 ; Kilic2017 ; Kilic2018 , eccentricity Linsenmeier2015 , or the length of the year LucAstr2013 ; Abbot2018 can have a dramatic effect on the properties 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-free areas at much larger distance from the host star than planets without seasonal variability, which leads to a substantial expansion of the outer edge of the habitable zone. Additionally, one expects that tidally locked planets with an active carbon cycle can never be found in a SB state Checlair2017 .

I.1 Multistability of the Climate System

Further investigations have proposed the possibility of the existence of alternative cold states with respect to the SB one. Such states are characterised by the existence of a thin strip of ice-free region near the equator. Clearly, this possibility has key implications in terms of habitability and the evolution of life on Earth. Different physical mechanisms have been proposed for explaining the existence of such a state, based either on the role of a dynamical ocean lewis2007 or the specific properties of the albedo of sea ice abbott2011 . In fact, in a previous work Lucarini2017 we have found that, indeed, a third co-existing stable state, intermediate between the SB and the W climate, can be found 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 rather robust property of the climate system.

Indeed, the dynamical landscape of the Earth system might be even more complex than what is usually expected. A recent study brunetti2019 , performed using a rather sophisticated climate model run using aquaplanet boundary conditions (without continents), indicates the existence of at least five competing 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 complexities when one considers even warmer conditions. For sufficiently large values of $S^{*}$, the W climate state loses its stability as a result of the dramatic strengthening of the positive feedback associated to the presence of water vapour in the atmosphere. Indeed, warm conditions favour, through the thermodynamic effect associated with the Clausius-Clayperon relation, the retaining capacity of water vapour of the atmosphere. The water vapour is a powerful greenhouse gas, as it is active in the infrared radiation. As a result, when the concentration of water vapour is sufficiently high, the planet performs a transition to either the so-called moist greenhouse state or the so-called runaway greenhouse state (associated to a complete evaporation of the oceans) GomezLeal2018a ; GomezLeal2018b , which define the warm (or inner) edge of the habitable zone Kastings1993 .

In what follows, we consider the simpler - yet extremely relevant - scenario where the only relevant co-existing climates are the SB state and the W state. As discussed in Lucarini2010a , the physics of the system is especially interesting when the critical transitions are approached. As the solar irradiance $S^{*}$ nears the critical value $S^{*}_{W\rightarrow SB}$ with the system being in the W state, the climatic engine becomes more efficient, because larger temperature gradients are realised inside the domain. Such an increased efficiency leads to a stronger atmospheric circulation, which is fuelled by temperature gradients and tends to reduce them by transporting heat from warm to cold regions, acting as a non-trivial diffusion process. Such a non-linear equilibration mechanism acts as a negative feedback and, broadly speaking, is a macroscopic manifestation of the second law of thermodynamics. One of the results of the heat transport performed by the atmospheric circulation is the stabilization of the ice-line. When $S^{*}=S^{*}_{W\rightarrow SB}$, the ice-albedo feedback becomes as strong as the negative feedbacks of the system, and the system flips to SB state with the ice-line reaching the equator. Similar mechanisms are in place when the system is in the SB state and $S^{*}$ nears, instead, $S^{*}_{SB\rightarrow W}$. When $S^{*}=S^{*}_{SB\rightarrow W}$, the ice begins to melt near the equator, leading to a rapid poleward retreat of the ice line.

At the critical transitions the climate system is not anymore able to dampen fluctuations due to (infinitesimal) external forcings. Using methods borrowed from transfer operator theory Baladi2000 , the investigation of the behaviour of the same model used in Lucarini2010a has indeed shown that when $S^{*}$ nears $S^{*}_{W\rightarrow SB}$, the spectral gap - defined as the absolute value of the (negative) imaginary part of the subdominant Ruelle-Pollicott pole Pollicott1985 ; Ruelle1986 - of the transfer operator constructed in a suitably defined reduced phase space vanishes. As a result, exponential decay of correlation is lost and the system experiences what is often referred to as critical slowing down Tantet2018 ; see also Refs. Shiino1987 ; Pavliotis2014 .

Far from critical transitions, it has been shown Ragone2016 ; Lucarini2017b ; Lembo2019 that it is possible to perform climate change projections resulting from a time-dependent CO2 forcing using Ruelle response theory Ruelle2009 . In the case of perturbations not depending explicitly on time, response theory allows one to describe how the measure of the system changes differentiably with respect to small changes in the dynamics of the system. In the case of time-dependent perturbations, response theory makes it possible to reconstruct the measure supported on the pullback attractor Ghil2008 ; Chekroun2011 ; CLR13 (see also the closely related concept of snapshot attractor Ott.ea.1990 ; DBT15 ) the non-autonomous system through a perturbative approach around a reference state, which, in the case of the climate studies referred to here, corresponds to the pre-industrial conditions. When we are nearing a critical transition, it is reasonable to expect a monotonic decrease of the spectral gap Dembo2010 ; Assaraf2018 . Hence, in the vicinity of the critical transitions the presence of a vanishing spectral gap leads to having a vanishing radius of expansion for response theory Lucarini2016 . Indeed, it is expected that the (near) closure of the spectral gap is associated to a strongly enhanced sensitivity of the system’s statistics to perturbations Chekroun2014 . See Ghil2019 for a thorough discussion on the various regimes of climatic response to forcings and of the relationship between climate change and climate variability across various temporal scales.

I.2 Melancholia States of the Climate System

A key question is what lies in-between the stable co-existing climates within the region of the parameters space where bistability is found. In simple models, it is often possible to identify unstable solutions sitting in-between the two stable climates. Such unstable solutions are embedded in the boundary between the two basins of attraction and, are roughly speaking, ice-covered up to the mid-latitudes. These solutions are saddles because they attract orbits starting from initial conditions on the basin boundary but are not stable, as small generic perturbations pushing the orbit outside the basin boundary with probability one and then lead to the system falling eventually into one of the competing attractors Ghil1976 ; Bodai2014 .

When studying more comprehensive climate models featuring chaotic dynamics, things are, as described in the next section, considerably more complex from a mathematical point of view, and the individuation of the unstable saddle solutions is much harder Grebogi1983 ; Robert2000 ; Ott2002 ; LT:2011 . Since these solutions are unstable, they cannot be found by direct numerical simulation. In a previous investigation Lucarini2017 , we adapted the edge tracking algorithm Skufca2006 ; Schneider2007 introduced for constructing the edge states, i.e. the special solutions separating laminar from long-lived turbulent regimes of motion in a fluid dynamical setting222Note that in such systems, technically, there are no competing attractors, because the only true attractor is the laminar state.. We used a recursive technique of bisection on the initial conditions333A different approach based on control theory aims at finding unstable solutions by a feedback loop involving changes in the value of the control parameter defining the region of bistability Sieber2013 ; Sieber2014 . for shadowing trajectories on the basin boundary separating the two co-existing W and SB states, and managed to populate the corresponding saddles.

The analysis was performed using an intermediate complexity climate model with $O(10^{4})$ degrees of freedom. The saddles we found had the remarkable property of featuring chaos (we found evidence that the first Lyapunov exponent was positive), and were named as Melancholia (M) states. Figures 1a)-b) summarize the main properties of the system, by showing the long term averages of the globally averaged surface temperature $\langle T_{S}\rangle$ and of the temperature difference between low and high latitudes $\Delta T_{S}$ as a function of the relative solar irradiance $\mu=S^{*}/S^{*}_{0}$, where $S^{*}_{0}$ is the present day value. By focusing on the M states we have been able to show the existence of much richer than previously thought dynamical landscape.

As discussed in Lucarini2017 , up to $\mu\sim 1.01$, the M state is characterised by longitudinal symmetry in its statistical properties, just as the boundary conditions of the system, are, indeed, longitudinally symmetric. The chaotic dynamics manifests itself as weather variability in a form not too dissimilar from the usual one observed in stable climates. Nonetheless, on long time scales, orbits initialised near the M states drift to either the W or the cold SB state, as a result of the dominating positive ice-albedo feedback. For $\mu\sim 1.01$, we observed that the symmetric M state becomes transient, evolving very slowly (on a time scale scale much longer than the other ones typical of the 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 feature rather different dynamical behaviour and the boundary between them rotates very slowly in time. The nontrivial bifurcation associated to such a symmetry break leads to dynamical regimes that resembles chimera states in extensive systems Abrams2004 ; Omelchenko2018Nonlinearity . The third climate state mentioned before exists in a small parametric window near $\mu\sim 1.045$ Lucarini2017 .

We remark that the dynamical systems viewpoint clarifies that the critical transitions for the W (snowball) state occurring when $S^{*}$ approaches the critical value $S^{*}_{W\rightarrow SB}$ ($S^{*}_{SB\rightarrow W}$) are associated to the collision between the M state and one of the W (snowball) climates, according to the dynamical scenario of boundary crisis Ott2002 . The system’s reduced ability to dampen fluctuations near the tipping points and the associated shrinking of the spectral gap described above can be seen, dynamically, as the result of the fact that the attractor attracts orbits less effectively in its immediate neighbourhood because of the presence of a nearby M state state Tantet2018 .

I.3 This Paper: Goals and Main Results

Studying multistable systems in general, and the climate system in particular, using deterministic autonomous dynamical systems faces two important issues, both resulting from the fact that the phase space is partitioned into disjoint invariant sets - the various basins of attractions and their respective boundaries. First, it is not possible to account for transitions between the co-existing basins of attraction. Instead, transitions between distinct regimes of motion are observed in many systems of interest. Additionally, one cannot establish an ergodic, physically relevant invariant measure, 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 the system of interest to be found at a given time in a specific regime of motion. Instead, all the statistical 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 LB2018 and giving them a much broader scope and setting them a in much more robust mathematical framework, we want to address these points. For the benefit of the reader, we report below the main goals of our investigation and the main findings presented in this paper.

We will introduce a stochastic forcing to the model studied in Lucarini2017 as a Gaussian perturbation with variance $\sigma^{2}$ modulating the intensity of the incoming solar radiation. This impacts the radiative budget of the planet in a very nontrivial way and, thanks to the ice-albedo feedback, acts as a multiplicative noise. Additionally, we conjecture, supported by physical and qualitative mathematical arguments, that the combination of noise and nonlinear deterministic dynamics leads to a hypoelliptic diffusion process, i.e. the noise propagates to all degrees of freedom of the system Bell2004 . The presence of a random forcing allows the system to perform transitions between the neighbourhoods of the deterministic attractors, crossing the basin boundaries that, in the unperturbed case, are impenetrable; see some classical results in Hanggi1986 ; Kautz1987 ; Grassberger1989 . The consideration of a random forcing will allow us to construct the ergodic invariant measure of the system by performing long integrations, and investigate the properties of noise-induced transitions.

Following Freidlin1984 ; Graham1991 ; Hamm1994 ; LT:2011 , we show that, in the weak noise limit and under suitable conditions, the invariant 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 $\sigma^{2}/2$. We will show how to compute the quasi-potential from the drift and volatility fields, and show that the drift field can be 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 the system in the absence of noise. In the region of bistability, the quasi-potential has local minima associated 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-potential has 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 attraction increases linearly with the product of $2/\sigma^{2}$ times the difference between the value of the quasi-potential at the M state through which the transitions takes place and its value at the corresponding attractor. We also show that the stochastically averaged exit trajectories connect the attractor with the M state, which is indeed the most likely exit point from the deterministic basin of attraction. In the weak noise limit, such paths correspond to the instantons Kautz1987 ; Grassberger1989 ; Kraut2002 ; Beri2005 ; Bouchet2016 . In our simulations, we show that such an identification becomes more accurate as the intensity of the noise is decreased. Exploiting this property, we also propose a new method for constructing the M states via direct numerical integration of stochastic differential equations.

We discover that, since generally the local minima of the quasi-potential corresponding to the two attractors have different values, in the zero-noise limit only one attractor - the one corresponding to the global minimum of the quasi-potential - is populated. Nonetheless, an individual trajectory may in fact persist near the competing metastable state for a very long time, as the permanence time in the corresponding basin of attraction also diverges (but at a slower rate than for the asymptotic state). As indicated by physical intuition, for low values of $S^{*}$, the noise selects as 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 $\mu=\mu_{crit}$, where $\mu_{crit}\approx 1.005$. For such a value of $\mu$, the equivalent of a first order phase transition for equilibrium systems takes place in the system. The order parameter that more obviously captures the transition is the globally averaged surface temperature.

The paper is structured as follows. In Section II we summarise the main mathematical concepts and ideas we have used to frame our investigation, to perform the data analysis, and to interpret our results. In Section III we report the modelling suite we have used, describe the data we have produced, and give information on how they have been processed. In Section IV we present and discuss our results. In Section V we propose and test numerically a new method for constructing the M states by looking at the intersections of instantons constructed by taking into account three different noise laws superimposed on to the same deterministic dynamics. In Section VI we summarize and interpret our findings, describe the limitations of our work, and propose ideas for future investigations. In Appendix A, we speculate on the fact that the framework presented in the paper is potentially suitable for clarifying the role played in complex historical processes by contingency, as discussed by Gould G1989 in the context of evolutionary biology.

II Mathematical Background

II.1 Geometry of the Phase Space: Attractors and M States

The investigation of systems possessing multiple steady states is an extremely active area of interdisciplinary research, encompassing mathematics, natural sciences, and social sciences. The recent review by Feudel et al. Feudel2018 , introducing a special journal issue on the topic, gives a rather complete overview of the state-of-the-art of the ongoing activities on this topic and provides several interesting examples. In the context of Earth sciences, it is more common to refer to critical transitions in multistable systems using the expression tipping points, which has received a great popularity following a paper by Lenton et al. Lenton2008 .

A possible way to introduce the mathematical background for multistable systems is the following. We consider a smooth autonomous continuous-time deterministic dynamical system defined on a smooth finite-dimensional compact manifold $\mathcal{M}$ (in what follows, a subset of $\mathbb{R}^{N}$). We assume that the dynamical system is dissipative, so that the phase space continuously contracts and the Lebesgue measure is not conserved. The orbit evolves from an initial condition $\mathbf{x}_{0}\in\mathcal{M}$ at time $t=0$. We define $\mathbf{x}(t,\mathbf{x}_{0})=S^{t}(\mathbf{x}_{0})$ as the orbit at a generic time $t$, where $S^{t}$ indicates the evolution operator. The corresponding set of ordinary differential equations written componentwise is:

[TABLE]

where $\mathbf{F}(\mathbf{x})=\textrm{d}/\textrm{d}\tau S^{\tau}(\mathbf{x})|_{\tau=0}$ is a smooth $N$-dimensional vector field. We define such a dynamical system as multistable if it features more than one asymptotic states. Specifically, we mean that there are two or more attractors $\Omega_{j}$, $j=1,\ldots,J$, each a corresponding basin of attraction having a finite Lebesgue measure. Each of the attractors is an invariant set and is the support of an invariant measure of the system. The asymptotic state where the trajectory falls into is determined by its initial conditions, and the phase space is partitioned between the basins of attraction $B_{l}$ of the various attractors $\Omega_{j}$ and the boundaries $\partial B_{l}$, $l=1,\ldots,L$ of such basins.

We assume, for simplicity, that within each basin of attraction an ergodic measure can be defined as the limit of the empirical measure constructed by averaging over infinitely long forward trajectories 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,\ldots,J$ ergodic components, each supported on the corresponding attractor $\Omega_{j},j=1,\ldots,J$.

We also assume that orbits initialized on the basin boundaries $\partial B_{l}$, $l=1,\ldots,L$ are attracted towards invariant saddles. In each basin boundary $\partial B_{l}$ there can be one or more of such saddles. We assume that each saddle features only one unstable direction. Such instability repels trajectories 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 the competing attractors can be chaotic (meaning here that at least one Lyapunov exponent is positive and unstable periodic orbits are dense). In some systems, chaotic dynamics can be realised on one or more of the saddles embedded in the basin boundaries Grebogi1983 ; Robert2000 ; Ott2002 ; Vollmer2009 . These are what we refer to as M states. M states are non-trivial geometrical objects being the support of a typically nontrivial measure LT:2011 .

A fascinating aspect of multistable systems is the following. If the first Lyapunov exponent of a chaotic saddle is larger than the inverse of the life time of the saddle state itself (which measures the rate at which orbits initialised near the saddle state are repelled towards the two nearby asymptotic states), then the basin boundary separating the basin of attraction of the two asymptotic sets has co-dimension strictly smaller than one. In two-dimensional maps, it has been proved that in this case the basin boundary is not a manifold, but rather a rough fractal Grebogi1983 ; Vollmer2009 ; LT:2011 . The authors have recently proposed a generalisation of this result to the case of $N$-dimensional maps BodaiLucarini2020 . In Lucarini2017 , despite high-dimensionality of the system, one could detect that the co-dimension of the basin boundary is strictly smaller than one. In fact, such co-dimension was found to be very close to zero, as a result of time scale difference between the most relevant instabilities: the climatic instability due to the ice-albedo feedback acting near the M state is much slower than the fast weather-like instability associated to baroclinic processes occurring on the M state. This implies that near the basin boundary there is basically no predictability of the second kind in the sense of Lorenz Lorenz1975 . The basin boundary is, de facto, a grey zone, and it is difficult to assess where orbits initialised near the boundary will end up.

II.2 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 on the deterministic dynamics discussed above. Processes that can be described as a noise-induced escape from an attractor have long been studied in the natural sciences; see Hanggi1986 ; Kautz1987 ; Grassberger1989 .

We generalise Eq. 1 by adding a stochastic component. We then consider a stochastic differential equation (SDE) in Itô form written as

[TABLE]

where $\mathbf{x},\mathbf{F}\in\mathbb{R}^{N}$, $\mathbf{F}$ and $s_{ij}\in\mathbb{R}^{N\times N}$ are smooth, $\sigma\geq 0$, $d{W}_{j}$ is the increment of an $N-$dimensional brownian motion, and $C_{ij}(\mathbf{x})={s}_{ik}(\mathbf{x}){s}_{jk}(\mathbf{x})$ is the noise covariance matrix. We consider the case of a hypoelliptic diffusion process. This amounts to assuming that while the covariance matrix of noise can be singular, the drift term modified according to the Stratonovich convention and the columns of the matrix $s$ satisfy the so-called Hörmander condition, i.e., the Lie algebra generated by them has dimension $N$ everywhere Hairer2011 . As a result of this, a smooth invariant density with respect to Lebesgue is realised because the noise is propagated to all the coordinates through the drift term; see Bell2004 444Assuming an elliptic diffusion process is extremely restrictive because it requires all degrees of freedom to be driven by Gaussian noise.. We discuss below in Sect. III.2 why such a mathematically important assumption can be heuristically justified in the specific case studied here.

Taking inspiration from the Freidlin-Wentzell Freidlin1984 theory and modifications thereof Graham1991 ; Hamm1994 ; LT:2011 , in the weak-noise limit $\sigma\rightarrow 0$ we seek a special functional form for the invariant measure. Indeed, we look for a large deviation law:

[TABLE]

where the rate function $\Phi(\mathbf{x})$ is referred to as quasi-potential, and we have neglected the pre-exponential term. Specifically, the symbol $\sim$ in Eq. 3 implies that $\Phi(\mathbf{x})=-1/2\lim_{\sigma\rightarrow 0}\sigma^{2}\log\Pi_{\sigma}(\mathbf{x})$. The function $\Phi(\mathbf{x})$ can be obtained as follows. We solve the stationary Fokker-Planck equation corresponding to Eq. 2:

[TABLE]

where $\mathbf{J}$ is the current density. We then consider the weak noise limit, and use as ansatz the expression given in Eq. 3. We obtain the following Hamilton-Jacobi equation Gaspard2002 :

[TABLE]

This equation allows one to express $\Phi$ in terms of the drift and volatility fields introduced in Eq. 2. The quasi-potential $\Phi$ can also be computed by solving the variational problem associated with the Freidlin-Wentzell action Nardini2016 . Finally, alternative routes for computing $\Phi$ have been proposed in Ao2004 ; Yin2006 .

The explicit computation of $\Phi$ is far from trivial, yet of great interest in many applications; see e.g., Zhou2012 for the case of biological systems. Brackston et al. Brackston2018 have recently proposed an algorithm for estimating $\Phi$ in the case that the governing equations are polynomial and involves solving an optimization over the coefficients of a polynomial function. Instead, Tang et al. Tang2017 proposed a variational method for estimating in the populations corresponding to each deterministic attractor without resorting to computing the invariant measure.

Following Graham1991 ; Hamm1994 , we now describe the dynamical meaning of $\Phi$. Indeed, solving the previous Hamilton-Jacobi equation corresponds to the fact that it is possible to write the drift vector field as the sum of two vector fields:

[TABLE]

that are mutually orthogonal, so that ${R}_{i}(\mathbf{x})\partial_{i}\Phi(\mathbf{x})=0$. In the case Eq. 2 describes a thermodynamical system near equilibrium, $\mathbf{R}$ defines the time reversible dynamics, while $\mathbf{F}-\mathbf{R}$ defines the irreversible, dissipative dynamics Graham1987 . One finds that

[TABLE]

As a result of this, $\Phi$ has the role of a Lyapunov function whose decrease describes the convergence of an orbit to the attractor. Specifically, $\Phi(\mathbf{x})$ has local minima at the deterministic attractors $\Omega_{j}$, $j=1,\ldots,J$, and has a saddle behaviour at the saddles $\Pi_{l}$, $l=1,\ldots,L$. If an attractor (saddle) is chaotic, $\Phi$ has constant value over its support, which can then be a strange set Graham1991 ; Hamm1994 .

Note that in the standard case of dynamics taking place in an energy landscape defined by a (confining) potential $U(\mathbf{x})$ and noise correlation matrix proportional to the identity (obtained by setting ${F}_{i}(\mathbf{x})=-\partial_{i}U(\mathbf{x})$ and $C_{ij}(\mathbf{x})=\mathbf{1}$ in the previous equations) one has $\Phi=U$. Additionally, one derives $\dot{U}(\mathbf{x})=-\partial_{i}U(\mathbf{x})\partial_{i}U(\mathbf{x})<0$ and $U(\mathbf{x})$ is a Lyapunov function, and one recovers an equilibrium state, where detailed balance applies and, by definition, the current vanishes ($\mathbf{J}=0$)555Note that, in general, we can have an equilibrium state if and only if the the drift term has a gradient structure with respect to the metric defined by the noise covariance tensor Pavliotis2014 ..

We remark that the function $\Phi(\mathbf{x})$ is defined globally but is not, in general, twice differentiable everywhere. Indeed, discontinuities in its first derivatives are present if a) the Hamiltonian associated with the Hamilton-Jacobi equation given in Eq. 5 is not integrable (non-integrability being the typical situation), and b) if the system features more than one co-existing attractors. These latter discontinuities are of little practical relevance because they appear only for values of $\Phi$ larger than those at the saddles Graham1986 , for reasons that will become apparent below.

II.3 Noise-induced Escape from the Attractor

The quasi-potential $\Phi$ is key for determining the statistics of noise-induced escape from a given attractor. Indeed, the probability that an orbit with initial condition in $B_{j}$ does not escape from it over a time $p$ decays as:

[TABLE]

where $\tau_{j,\sigma}$ is the expected escape time and $\Delta\Phi_{j}=\Phi(\Pi_{l})-\Phi(\Omega_{j})$ is the lowest quasi-potential barrier height LT:2011 , i.e. $\Phi$ has the lowest value in $\Pi_{l}$ compared to all the other saddles neighbouring $\Omega_{j}$. In general, one may need to add a correcting prefactor to $P_{j,\sigma}(p)$ LT:2011 .

Note that $\tau_{j,\sigma}$ given in Eq. 8 does not contain the pre-exponential factor. Ref. Bouchet2016 provided an expression for such pre-exponential factor for general non-equilibrium diffusion processes under the assumption that attractors and saddles are simple points, thus generalising what is given in Bovier2004 . As we aim at treating also a more general setting for the geometry of attractors and saddles, we pay below the price of having to take the pre-exponential factors as phenomenological parameters that one can find from experiments or numerical simulations Bodai2018 . We also remark that, in the zero-noise limit, the transition paths outside a basin of attraction follow the instantons. Instantons are defined as solutions of

[TABLE]

that connect a point in $\Omega_{j}$ to a point in $\Pi_{l}$. Instantonic trajectories have a reversed component of the 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 $\Omega_{1}$, $\Omega_{2}$, and one saddle $\Pi_{1}$. We can then express the average transitions times as follows:

[TABLE]

so that

[TABLE]

This implies that, in the weak noise limit, both escape times diverge, but the escape time out of the attractor corresponding to the lower value of the quasi-potential diverges faster. Note that one can expect the proportionality constant in Eq. 12 to be $O(1)$. Taking a maximally coarse-grained view on the problem, where we consider the state as represented by the populations $P_{1,\sigma}$, $P_{2,\sigma}$ of the neighbourhood of the two attractors, we can write the following master equation:

[TABLE]

The master equation above makes sense if one assumes the presence of clear timescale separation between the relaxation motions near each attractor and those across the saddle, which depends critically on the presence of weak noise Lelievre2015 ; Gesu2019 . At steady state, we obtain that

[TABLE]

Equation 15 could be obtained by integrating the invariant measure given in Eq. 3 in the neighborouhood of the attractors and taking a saddle point approximation. Additionally, Eq. 15 implies that in the weak-noise limit only one of the two deterministic attractors will be populated, and specifically the one where the quasi-potential has lower value. We remark that two different noise laws differing for the correlation matrix $C$ acting on top of the same drift field will define two different quasi-potentials, see Eq. 5. As a result of that, they will in general feature a different selection of the dominating population in the zero noise limit. One can easily extend the master equation defined above to the case where multiple states and multiple paths of transitions are present. Finally, note that in Lucarini2019 the mathematical framework described in this section has been used to study stochastic resonance for general non-equilibrium systems.

III Numerical Modelling

The climate model considered here is constructed by coupling the primitive equations atmospheric model PUMA puma with the Ghil-Sellers energy balance model Ghil1976 , the latter describing succinctly the meridional oceanic heat transport. It has been already presented in Lucarini2017 with the name of PUMA-GS, but we report here again its formulation in order to elucidate the role of stochastic forcing, which was absent in the previous version. The stochastic forcing is added as a fluctuating term modulating the value of the incoming radiation determining the energy input into the system.

III.1 The Atmospheric Component

The atmospheric component of the PUMA-GS model is provided by PUMA puma , which consists of a dynamical core: the dry hydrostatic primitive equations on the sphere (mapped laterally by the latitude $\phi$ and longitude $\lambda$), solved by a spectral transform method (only linear terms are evaluated in the spectral domain, nonlinear terms are evaluated in grid-point space). The equations for the prognostic state variables, the vertical component (with respect to the local surface) of the absolute vorticity $\zeta=\xi+2\nu\Omega_{E}$ (where $\xi$ is the vertical component of the relative vorticity, $\nu=\sin\phi$, and $\Omega_{E}=2\pi$/day is the angular frequency of the Earth rotation) the (horizontal) divergence of the velocity field $D$, the (atmospheric) temperature $T_{a}=\bar{T_{a}}+T_{a}^{\prime}$ (separated into a time-independent arbitrary reference part $\bar{T_{a}}$ and anomalies $T_{a}^{\prime}$), and the logarithmic pressure (normalized by the surface pressure $p_{s}$) $\sigma=\ln p/p_{s}$, read as follows:

[TABLE]

where $s^{2}=1/(1-\nu^{2})$, $F_{u}=V\zeta-\dot{\sigma}\partial_{\sigma}U-T_{a}^{\prime}\partial_{\lambda}\ln p_{s}$, $F_{v}=-U\zeta-\dot{\sigma}\partial_{\sigma}V-T_{a}^{\prime}s^{-2}\partial_{\nu}\ln p_{s}$, $U=u\cos\phi$, $V=v\cos\phi$, $u$, $v$ being respectively the horizontal and vertical wind velocity components, and $\Psi$ is the geopotential height. Equations 16,III.1, and 19 express the conservation of momentum, Eq. III.1 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 stability of the model. Firstly, the hyperdiffusion operator $K\nabla^{8}$ is added to the equations of vorticity, divergence and temperature, to represent subgrid-scale eddies. Secondly, large-scale dissipation of vorticity and divergence is facilitated by Rayleigh friction of time scale $\tau_{f}$. Thirdly, the physics of diabatic heating due to radiative heat transport is parametrized by Newtonian cooling: the temperature field is relaxed (with a time scale $\tau_{c}$) towards a reference or restoration temperature field $T_{R}$, which can be considered a radiative-convective equilibrium solution. We adopt the following simple expression for the restoration temperature puma :

[TABLE]

where $(T_{R})_{tp}$ and $z_{tp}$ are the temperature and height of the tropopause, respectively, $L$ ($\bar{L}$) is the (average) lapse rate, $\langle T_{S}\rangle$ is the globally averaged surface temperature, and $z(\sigma)$ is determined by an iterative procedure puma . The above expressions indicate that the restoration temperature profile is anchored to the surface temperature $T_{S}$. However, as Eq. 22 indicates, $T_{R}$ at any one point on the sphere is determined by not only the local (dynamical) surface temperature, but also the global average $\langle T_{S}\rangle$. We note that $T_{a}(\sigma=1)$ is obtained by linear extrapolation, according to $T_{a}(\sigma=1)\approx T_{a}(\sigma=0.9)+\eta(T_{S}-T_{R}(\sigma=0.9))$, $0<\eta<1$. With $\eta=1$ the coupling term is $k_{3}(T_{a}(\sigma=1)-T_{S})\approx k_{3}(T_{a}(\sigma=0.9)-T_{R}(\sigma=0.9))$. Generally $\overline{T_{a}(\sigma=1)-T_{S}}\neq 0$ (laterally inhomogeneous heating), but $\overline{\langle T_{a}(\sigma=1)\rangle}=\overline{\langle T_{S}\rangle}$, where the overbar denotes averaging with respect to time.

For our setup we choose: $K^{-1}=0.25$ days, $\tau_{c}=30$ days, $\tau_{f}=1$ day, $\bar{L}=0.0065$ K/m, $z_{tp}=12000$ m, and $k_{3}=10^{-4}$. We also adopt a coarse resolution of T21 (i.e., the series of spherical harmonics are triangular-truncated at total wave number 21). This implies the optimal number of Gaussian grid points: $N_{lon}=2N_{lat}=64$. Finally, we consider $N_{lev}=10$ vertical layers and consider a vanishing orography, so that we have a zonally-symmetric configuration. The equations are integrated numerically using a $\Delta t=1$ [hour] time step size.

III.2 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 Ghil1976 ; Bodai2014 . This model includes a simple yet effective representation of the ice-albedo feedback, and basically defines the slow manifold of the coupled atmosphere-ocean system. The partial differential equation describing the evolution of the ocean surface temperature field $T_{S}=T_{S}(t,\phi,\lambda)$ is:

[TABLE]

where $S_{0}^{*}$ is the present solar irradiance666The factor 4 emerges as a result of the geometry of the Earth-Sun system saltzman_dynamical , $\mu=S^{*}/S^{*}_{0}$ as introduced in Sect. I, while the heat capacity $C(\phi)$ and the geometrical factor $I(\phi)$ are explicitly dependent on $\phi$ only, thus enforcing zonally-symmetric boundary conditions. The albedo $\alpha$ depends on $\phi$ and, critically, on $T_{S}$, with a rapid transition from strong albedo for low values of $T_{S}$ ($\alpha_{max}=0.6$) to weak albedo for $T_{S}\gtrsim 260$ K ($\alpha_{min}=0.2$), which fuels the positive ice-albedo feedback. Additionally, $O$ is the outgoing radiation per unit area, expressed as a monotonically increasing function of $T_{S}$ (this is responsible for the negative Boltzmann feedback, taking into account also the greenhouse effect), $\bar{D}_{\phi}$ is a diffusion operator parametrizing the meridional heat transport, and $\chi$ describes the heat exchange with the atmosphere. See Bodai2014 ; Lucarini2017 for further details.

Finally, the last term on the right hand side $s.f.$ is the stochastic forcing, which is introduced as a random modulation of the solar irradiance given by $\mu S_{0}^{*}$. Hence, we have:

[TABLE]

where $\sigma$ controls the intensity of the noise, $s$ defines the noise law, and $dW$ is the increment of a one-dimensional Wiener process. Since $s$ depends explicitly on $T_{S}$ via the term $\alpha(\phi,T_{S})$, we are dealing with a multiplicative noise law. We consider the Itô convention for noise, so that our (discretised) equations are in the form of Eq. 2. See a discussion on the relevance of the chosen convention in Sect. IV.5. Adding a Gaussian random variable of variance $\sigma$ at each time step $\Delta t$ ($1$ hour ) of the model amounts to considering that, on the time scale $\tau=N\times\Delta t$, the relative fluctuation of the solar irradiance scales as $\sigma_{\tau}=\sigma/\sqrt{N}$.

As mentioned in Sect. II.2, our approach requires assuming the validity of the hypoellipticity condition. In order to prove this, we should test the Hörmander condition for the evolution equations 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 we can heuristically understand why, indeed, it is reasonable to assume that stochastic forcing acting on the oceanic surface temperature propagates to all degrees of freedom of the coupled system, as usually implicitly assumed in basically any numerical study of stochastically forced geophysical flows 777Obviously, numerical truncation introduces some additional noise on all degrees of freedom of the system..

We can first approach this problem by looking at the structure of the evolution equations. The stochastic forcing given in Eq. 25 impacts directly the $T_{S}$ field, as shown in Eq. 24. The $T_{S}$ fields determines the restoration temperature $T_{R}$, see Eq. 21. In turn, the restoration temperature impacts the anomaly of the atmospheric temperature field $T_{a}^{\prime}$ (see Eq. III.1), which in turn affects the vorticity field $\zeta$ - Eq. 16 - and divergence field $D$ - Eq. III.1. Finally, anomalies in $D$ impact the surface pressure $p_{s}$, as clear from Eq. 19. The nonlinear terms corresponding to advective processes on the 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 better understood by taking a truncated Fourier representation of Eqs. 16-19, which, in fact, closely corresponds to the actual formulation of the numerical model implemented here. Concluding, no dynamical 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 only a small fraction of the degrees of freedom of the system.

On more physical grounds, one can observe that the ocean surface temperature drives the restoration profile of the atmospheric temperature, and that fluctuations of such a profile modulate the atmospheric instabilities, whose energy cascades down to the smallest scales resolved by the model; see Peixoto1992 ; Lucarini2014b ; Ghil2019 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 $\mu=0.98$ (close to the tipping point $\mu_{W\rightarrow SB}$), $\mu=1.0$ (corresponding to present-day solar irradiance), and $\mu=1.02$ (in the parametric region where the M state undergoes a symmetry-break bifurcation). We then construct the weak-noise limit of the invariant measures for all the values of $\mu$ in the region of bistability. We remark that the results shown in Figs. 2, 3, and 5a) have already been reported in the short communication LB2018 , but we deem extremely useful to present them here as well because they are now part of a much more complex, coherent, and detailed narrative.

IV.1 Escapes from Basins of Attraction and Instantons

In the case of $\mu=0.98$, we first perform a set of simulations with noise of different intensity rangiW ng from $\sigma_{\tau}=0.5\%$ to $\sigma_{\tau}=1.4\%$, with $\tau=100$ years ($y$). For each value of the noise intensity, we initialise 50 trajectories in the basin of attraction of the W climate and study the statistics of the escape time to the SB attractor. When the transition takes place, we stop the integrations. We observe (not shown) that for each value of $\sigma_{\tau}$ the escape times are to a good approximation exponentially distributed, thus obeying Eq. 8; the process of transition behaves like a Poisson process. The results on the expectation value of the transition times are presented in Fig. 2, where we show that, indeed, $\tau_{\sigma}$ agrees with the prediction of Eq. 8. Hence, it is possible to define the difference between the value of the quasi-potential $\Phi$ at the M state and at the W attractor as the slope of the straight line. For reference, we have that for $\sigma_{100y}=0.5\%$ the average escape time is about $5.2\times 10^{3}$ $y$. We can predict that the escape rate increases to about $1.5\times 10^{7}$ $y$ when $\sigma_{100y}\sim 0.3\%$. We have that the slope of the straight line in Fig. 2 gives $\sim 2(\Phi(M)-\Phi(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. The operation can be repeated for all the other values of $\mu$ in the range of multistability and for the processes 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 Bodai2014 ; Lucarini2017 , we choose to consider the reduced phase space spanned by the globally averaged surface ocean temperature $\langle T_{S}\rangle$ and by the meridional temperature difference $\Delta T_{S}$, defined as the difference between the spatially averaged ocean temperature field between the Equator and $30^{\circ}N$ and between $30^{\circ}N$ and the North Pole. This reduced phase space provides a minimal yet physically informative viewpoint on the problem, because it is directly linked with the main physical processes occurring in the climate model:

• •

The average surface temperature $\langle T_{S}\rangle$ is directly associated to the positive ice-albedo feedback and the negative Bolzmann radiative feedback;

• •

The meridional temperature difference $\Delta T$ controls the meridional heat transport performed by the ocean, as a result of the diffusive law we insert into its evolution equation;

• •

The meridional temperature gradient $\Delta T$ also controls the meridional heat transport performed by the atmosphere, as a result of the mechanism of baroclinic instability Holton .

Figure 3 depicts, for the case $\sigma_{\tau}=1.0\%$, the transient two-dimensional distribution function $\tilde{\rho}$ constructed using a frequentist approach using the 50 simulations described above, where the statistics is collected only until the $W\rightarrow SB$ transition is realised. The distribution we obtain cannot be interpreted as an approximation of the invariant measure, because the integrations are stopped after the transitions. Nonetheless, it is apparent that the transitions take prominently place in a very narrow band linking the W attractor and the M state888Note 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 a very small variability in the deterministic chaotic case.. In order to obtain a better understanding of the transition paths, we construct an estimate of the instanton linking the W attractor to the M state and associated to the $W\rightarrow SB$ transitions by conditionally averaging the trajectories according to the value of $\langle T_{S}\rangle$. To a good approximation, the instanton connects the W attractor to the M state, and follows a path of decreasing density. We do not find evidence of different paths for the trajectories leading to an escape and the relaxation trajectories, which is, instead, a typical signature of non-equilibrium ZinnJustin1996 . This can be explained by considering Bodai2014 , where it is shown that the ocean model evolve to a good approximation in an energy landscape. See also Section V.

IV.2 Construction of the Invariant Measure

The information contained in Fig. 3 is limited because we are studying only $W\rightarrow SB$ escape processes, and we do not allow for the establishment of an invariant measure. The problem lies in the fact that the quasi-potential minimum associated to the $SB$ attractor is much deeper than the one associated to the $W$ attractor, so that the average escape time associated to $SB\rightarrow W$ transitions is prohibitively long for the range of (rather weak) noise intensities used for constructing Fig. 3. In fact, in order to be able to construct the invariant measure of the system, we need to observe a sufficient number of $W\rightarrow SB$ and $SB\rightarrow W$ transitions, in order to be sure that we have collected a satisfactory statistics; see also the master equations for the populations given in Eqs. 13-14. We next increase the noise intensity by setting $\sigma_{100y}=1.8\%$, so that, in an integration lasting about $1.0\times 10^{4}$ $y$, we observe 10 $SB\rightarrow W$ and $W\rightarrow SB$ transitions. The number of transitions 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 measure on the $(\Delta T_{S},\langle T_{S}\rangle)$ plane; we refer to this also as the two-dimensional probability distribution function (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 agreement is even clearer when considering the two marginal pdfs, constructed by projecting the invariant measure on the one dimensional spaces defined by $\Delta T_{S}$ and $\langle T_{S}\rangle$. Note that the peak over the $W$ state is hardly noticeable because the occupation of the state is extremely low (less than 5% of the total); compare with the cases studied below where higher values of $\mu$ are considered. We are also able to construct both the $W\rightarrow SB$ and the $SB\rightarrow W$ instantons, whose starting points agree remarkably well with the position of the $W$ and the $SB$ attractors, respectively, while their final points are in good agreement with the position of the M state. By constructing the pdf along the instantons, we find that they follow a path of monotonic descent (indeed, they follow closely the crests of the pdf), with the minimum located at the M state. Again, this last property can be hardly visualised in Fig. 3 for the $W\rightarrow SB$ instanton (see inset), because the population near the $W$ state is quite small.

Next, we repeat the analysis for $\mu=1$; results are shown in Figs. 5a)-b). In Panel a) we consider an integration with $\sigma_{100y}=1.5\%$ lasting about $2.9\times 10^{4}$ $y$ and characterised by 41 $SB\rightarrow W$ and $W\rightarrow SB$ transitions. We have an occupation of about $30\%$ for the $W$ basin of attraction, and of about $70\%$ for the $SB$ basin of attraction; additionally, we have $\tau_{\sigma}^{W\rightarrow SB}\sim 210y$ and $\tau_{\sigma}^{SB\rightarrow W}\sim 460y$. Also in this case, the projection of the invariant measure in the $(\Delta T_{S},\langle T_{S}\rangle)$ plane shows that there is good agreement between the position of the peaks of the pdfs and the attractors, and that the estimates of the instantons connect attractors and M states with a good precision. It is also clear that the instantons follow a path of descent in terms of probability, as shown by the central inset.

In Panel b) we show the results of repeating the analysis for $\sigma_{100y}=1.8\%$. In this case the simulation lasts about $2.7\times 10^{4}$ $y$ and we obtain 73 $SB\rightarrow W$ and $W\rightarrow SB$ transitions, and we 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 $35\%$ for the $W$ basin of attraction, and of about $65\%$ for the $SB$. Instead, the average escape times change more substantially, and can be estimated as $\tau_{\sigma}^{W\rightarrow SB}\sim 160y$ and $\tau_{\sigma}^{SB\rightarrow W}\sim 300y$. These last two results indicate that the difference between the value of the quasi-potential at the two competing attractors is relatively small. We will explore this matter in Sect. IV.4, where we will try to deduce where the quasi-potential reaches its absolute minimum for each value of $\mu$ in the bistable region.

It is worth looking more in detail at how the estimate of the instantons is impacted by the intensity 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 confirmed by the results shown in Fig. 6: we have that the estimates of the instantons obtained using lower noise intensity come closer to the attractors and to the M state. Nonetheless, also the instantons obtained for very strong noise are still relatively accurate.

IV.3 Instantons and Transitions across the Symmetry-broken Melancholia State

We next examine the noise-induced transitions for $\mu=1.02$. This case is quite interesting because, as discussed in Lucarini2017 and reported in Fig. 1, the longitudinally-symmetric M state is transient with a very long life time, and slowly evolves into a symmetry broken M state featuring a relatively cold and a relatively warm region, separated by two small regions with large longitudinal temperature gradient at all latitudes. It seems relevant to test whether noise-induced transitions take place through the transient, symmetric M state or the true, symmetry-broken one. Results are shown in Fig. 7, where we use data from a simulation lasting about $10^{4}$ $y$ with $\sigma_{100y}=1.5\%$. We observe only 6 transitions in both directions. This time, as opposed to the case of $\mu=0.98$, the figure is so low because it is extremely hard to escape from the $W$ state. We estimate the escape times as $\tau_{\sigma}^{W\rightarrow SB}\sim 1400y$ and $\tau_{\sigma}^{SB\rightarrow W}\sim 140y$. We portray the logarithm of the invariant measure of the system in the usual projected space, and the estimates of the instantons. We first observe that in this case most of the density is concentrated around the W attractor, and a nontrivial relation exists between the paths of the instantons and the dynamical structures on the basin boundary. It is clear that the transient M state plays the role of the gateway of transitions as seen in the previous cases, despite its transient nature. The noise-induced transitions do not go through the actual M state (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 result from the consideration of noise with finite (and not infinitesimal) strength: the quasi-potential near the transient M state might be just barely higher than that of the asymptotic M state (and possibly with a more favourable pre-exponential factor), so that a small but finite noise perturbation might push an orbit near the transient M state into the other basin of attraction. In order to address this point and find instantonic paths connecting the attractors with the true M states, one might need to resort to using more sophisticated numerical techniques. In particular, one should consider using rare events algorithms rubino2009rare ; Ragone2017 to rigorously construct instantonic trajectories Grafke_2013 . This is beyond the current abilities of the authors but definitely deserves attention in future studies.

IV.4 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 on the attractor featuring the lowest value for the quasi-potential $\Phi$. Since for $\mu<S^{*}_{W\rightarrow SB}/S^{*}_{0}$ only the SB state is realised, physical intuition suggests that for low values of $\mu$ within the range of bistability, the SB attractor should contain all the mass in the limit of weak noise, as one can also anticipate by looking at Fig. 3. Conversely, one expects that for high values of $\mu$ within the range of bistability, the $W$ should be dominant in the weak noise limit. It is reasonable (yet far from obvious or rigorous) to expect that there should be a critical value of $\mu=\mu_{crit}$ separating the two regimes. Following Lucarini2017 , we consider 18 equally spaced values of $\mu$ ($\Delta\mu=0.005$) within the multistable regime. For each of these values of $\mu$ (excluding the case of $\mu=1.045$, where three stable states are realised) the fraction of the population residing within the basin of attraction of the deterministic W attractor $P_{W,\sigma}(\mu)$ and its complement, residing in the basin of attraction of the SB attractor $P_{SB,\sigma}(\mu)$. Related results are shown in Fig. 8. We show how the probability distribution of the variable $\langle T_{S}\rangle$ depends on $\mu$ for three different noise levels: $\sigma_{100y}=1.5\%$ (Panel a), $\sigma_{100y}=1.8\%$ (Panel b), and $\sigma_{100y}=2.5\%$ (Panel c). In these panels we superimpose the bifurcation diagram reported in Fig. 1a). We observe that as the noise is reduced, for all values of $\mu$ the distributions are a) more peaked around the attractors; and b) one of the attractors becomes clearly dominant.

In Panel d) the plot for each value of $\mu$ the integral of the pdfs reported in the three panels a), b), c) up to the values of $\langle T_{S}\rangle$ corresponding to the M state (green continuous and dotted lines). To a very good degree of approximation, this corresponds to the integral of the invariant measure over the support of the deterministic basin of attraction of the $SB$ climate. We obtain that for decreasing values of the noise intensity, the emerging invariant measure converges to the deterministic measure supported on the SB attractor for $\mu\leq\mu_{crit}\approx 1.005$, while the invariant measure converges to the deterministic measure supported on the W attractor for $\mu\geq\mu_{crit}\approx 1.005$. The absolute minimum of the quasi-potential $\Phi$ is realised in the W attractor for $\mu\geq\mu_{crit}$ and in the SB attractor for $\mu\leq\mu_{crit}$. The changeover is, curiously, quite close to the reference case $\mu=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 that for all values of $\mu$ and $\sigma$ we have used simulations lasting at least $10^{4}$ $y$. For very low ($\leq 0.98$) and very large ($\geq 1.02$) values of $\mu$ in the considered range, the simulation length does not allow for 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 range an uncertainty on the figures reported in Fig. 8 of the order of the values of the smaller between $P_{SB,\sigma}(\mu)$ and $P_{W,\sigma}(\mu)$. This corresponds, in fact, to a low uncertainty, because most of the mass is concentrated near one of the two deterministic attractors. The uncertainty is quite small for $0.98\leq\mu\leq 1.03$, because, in all cases, we observe relatively many transitions. The uncertainty in this range can be safely estimated to be below 5%. Summarising, while the values reported in Fig. 8d) might have some non-negligible uncertainties, it seems that the estimate of $\mu_{crit}$ is quite robust.

Finally, we have verified that for all values of $\mu$ the escape time $\tau_{\sigma}^{W\rightarrow SB}$ and $\tau_{\sigma}^{SB\rightarrow W}$ grow rapidly with decreasing values of the intensity of the noise for all values of $\mu$. In agreement with what is shown in Fig. 8, we have that $\tau_{\sigma}^{SB\rightarrow W}$ grows more rapidly than $\tau_{\sigma}^{W\rightarrow SB}$ for $\mu\leq\mu_{crit}$ (and viceversa for $\mu\geq\mu_{crit}$). As a result, in the weak-noise limit an individual trajectory might be trapped for a very long time in the metastable, non-asymptotic state. We remark that we have not performed here a systematic evaluation of the exponential relationship between the escape times and $\sigma$, as instead done for $\mu=0.98$ and shown in Fig. 2, also because, as for the reasons explained before, this would require computational resources that are beyond what has been allocated for this study. We remark that this would allow for evaluating for each value of $\mu$ the difference between the value of the quasi-potential realised at the attractors and at the M state. The algorithm proposed in Bodai2018 can be very useful to reduce the computational burden.

We can say that for $\mu=\mu_{crit}$ our system exhibits a behaviour that is reminiscent of a first-order phase transition for near-equilibrium statistical mechanical systems, like a liquid-gaseous transition. In our case, $\mu$ is the control parameter (corresponding to the temperature in the equilibrium case), the quasi-potential $\Phi$ is the equivalent of a thermodynamic potential, the scaling factor for the noise intensity $\sigma$ is the equivalent of (square root of) the temperature, and the globally averaged surface temperature $\langle T_{S}\rangle$ is the natural order parameter, e.g. density. The discontinuous change in the properties of the system for $\mu=\mu_{crit}$ is associated to the change in the amount of absorbed and emitted radiation, as a result of the macroscopic change in the albedo of the planet due to the discontinuity in the position of the ice-line. We remark that choosing a different noise law would in general lead to a different value of $\mu_{crit}$, as a result of the fact that the functional form of $\Phi$ would be different.

IV.5 Relevance of the Choice of the Itô Convention for the Noise

A reasonable question to ask concerns to what extent our results are sensitive to the fact that we have chosen the Itô convention for the noise, which provides the starting point of the results discussed in Sect. II.2 and II.3. We argue that choosing other conventions would not alter essentially our findings because, to a first approximation, the stochastic forcing we have introduced can be treated as corresponding to perturbing the system with additive noise of different strengths near the cold and W attractors, plus a transition region between the two attractors (which is evidently very sparsely populated by the system), where the effective intensity of the noise decreases with the globally averaged surface temperature $\langle T_{S}\rangle$ and the multiplicative nature of the noise is more evident.

In the phase space region near the cold attractor, we have that $1-\alpha(\phi,T_{S})\sim 0.4$ since the temperature $T_{S}$ is extremely low and the planet is fully glaciated (or almost entirely so), so that $\alpha(\phi,T_{S})$ is virtually constant, with $\alpha(\phi,T_{S})\sim\alpha_{min}$. Near the W attractor, the properties of the field $\alpha(\phi,T_{S})$ are slightly more complex, because part of the planet is glaciated and part of it is ice-free. Nonetheless, to a first approximation, the ratio of the variance of the noise in the SB attractor vs. W attractor is of the order $((1-\alpha_{min})/(1-\alpha_{W}))^{2}\sim 3$. Loosely speaking, the competing W and SB climate states have different statistical mechanical, microscopic - as well as thermodynamical, macroscopic - temperatures.

V An Alternative Construction of the M States using Stochastic Perturbations

As discussed above, the construction of saddles for multistable systems is far from being a trivial task, and requires the use of the edge tracking algorithm introduced in Skufca2006 ; Schneider2007 and used also by the authors in Bodai2014 ; Lucarini2017 . We wish to provide here a proof of concept of an alternative procedure 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. The procedure 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 to the instability along and across the basin boundary might hinder an accurate computation of the saddles. Alternatively, it can be seen as a way to test the results obtained from the study of the deterministic dynamics.

The idea is to exploit the fact that, as discussed in Section II, under rather general conditions on the noise law, the saddle, in the weak noise limit, acts as the gate for noise-induced transitions between the competing attractors. We then propose to proceed as follows. Let’s consider the following stochastic differential equations:

[TABLE]

We consider the possibility of perturbing the deterministic flow with $K$ different noise laws, defined by the $K$ functions $\mathbf{s_{k}}(\mathbf{x})$, each leading to a noise with different covariance matrix $C_{ij,k}(\mathbf{x})$. We assume, for simplicity, that the deterministic system defined by ${\dot{{\mathbf{x}}}(t)}=\mathbf{F}(\mathbf{x(t)})$ is bistable. We choose $K$ noise laws such that they obey the hypotheses discussed in Sects. II.2 and II.3.

In the weak noise limit $\sigma\rightarrow 0$, the invariant measure of the $k^{th}$ SDE can be written as $\Pi_{\sigma,k}\sim\exp\left(-\frac{2\Phi_{k}(\mathbf{x})}{\sigma^{2}}\right)$. Clearly, for each noise law the quasi-potential $\Phi_{k}(\mathbf{x})$ is different. Nonetheless, as discussed above, in all cases $\Phi_{k}(\mathbf{x})$ has a local minimum (and is constant) on the support of the two deterministic attractors, and it is a saddle with constant value on the saddle separating the two attractors. We here have to assume that either the saddle is unique, or that all the quasi-potentials $\Phi_{k}(\mathbf{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-like component and the rest (see Eq. 6). Therefore, as suggested by Eq. 9, the instantonic paths are also different; yet, they connect the same attractors to the same saddle. Assuming that the attractors and the saddles are points or at least small sets given in a coarse-grained description of the phase space, we can identify the saddle as the only point in space where the instantons corresponding to all 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 Bodai2014 , differing for the law of the stochastic perturbation impacting the energy balance of the climate system. The Ghil-Sellers diffusive model can be written by removing the atmosphere-ocean coupling term in Eq. 24. In order to conform to Eq. 2 (note that we treat below the numerical discretization of a stochastically perturbed partial different equation), we express the three stochastically perturbed models as follows:

[TABLE]

where $dW$ is the increment of a one-dimensional Wiener process, and we have:

[TABLE]

Specifically, we have that $s_{1}(\phi,T_{S})=s_{1}(\phi)$ and $s_{2}(\phi,T_{S})=s_{2}(\phi)$ correspond to additive noise laws, which feature different diffusion matrices. Instead, $s_{3}(\phi,T_{S})$ is a multiplicative noise law as a result of the temperature-dependence of the albedo and is closely related to what was studied in the 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 noise amplitudes. Results are presented in Fig. 9, where we show the projections on the $(\Delta T_{S},\langle T_{S}\rangle)$ space. Panels a), b), and c) show the invariant measures and the instantons constructed for the noise laws $s_{1}$ (using $\sigma=0.4$) , $s_{2}$ (using $\sigma=2.0$), and $s_{3}$ (using $\sigma=1.0$), respectively. Note that the instanton constructed with the multiplicative noise law looks qualitatively different from what is show in Fig. 5-6, as a result of the lack of atmospheric motions in the simpler model discussed here. Panel d) portrays the instantons constructed for the three noise laws. Indeed, we have a confirmation 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 phase space can be identified through this geometric procedure, which is based exclusively on direct numerical simulations. Considering additional noise laws can be helpful in resolving possible geometrical degeneracies due to the use of projections. Projecting in more than two dimensions could also serve a similar scope and provide a better understanding of the alternative transition paths.

VI Conclusions

The goal of this paper has been the investigation of the properties of the noise-induced transitions across the multiple basins of attractions in an intermediate complexity climate model with $O(10^{4})$ 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 broad range of values of the solar irradiance. In a previous investigation we had been able to construct the full phase portrait of the deterministic version of the climate model considered here, and had constructed, beside the attractors, the M states of the climate system in the region of bistability Lucarini2017 .

The stochastic forcing is introduced here as a random modulation for the incoming solar radiation, and leads to a nontrivial multiplicative noise law, because the radiative forcing is affected by 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 equations indicates that for systems obeying the hypoellipticity condition, and for a suitable class of noise 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 the rate function can be identified with the quasi-potential. We have clarified how to compute the quasi-potential from the drift and volatility fields of the SDE and explained its property of being a Lyapunov function. The quasi-potential has local minima on the deterministic attractors, and has a saddle behaviour at the M states. Additionally, in the weak noise limit, transitions take place along special paths called instantons, which link the deterministic attractors and the M states. While, for a given deterministic dynamics, instantons corresponding to different noise laws follow different paths, they all link the same deterministic attractors to the same M states. We have shown how this property can be exploited to geometrically construct M states directly from direct numerical simulations of stochastic systems. Indeed, while the edge tracking algorithm applied to the deterministic system is a priori the preferred choice for finding M states, it might become nontrivial to implement in complex numerical models where the intermediate states constructed by bisection might correspond to regions where the model is numerically unstable, possibly because the realised physical fields are extremely exotic or non-realisable.

We have studied in detail the noise-induced transitions between the deterministic basins of attraction in the range of multistability, extending the results presented in a short communication LB2018 . We have shown that by studying how the average escape time depends on the intensity of the noise it is possible to estimate the difference between the value of the quasi-potential at the M state and at the attractor that the trajectories are escaping from. The estimates of the instantons are shown to become more precise as we consider simulations where weaker noise is considered. The instanton, in a case of special interest where the M state was shown to undergo a symmetry-break process, selects as optimal point of passage between the SB and the W climate the transient M state 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 the intensity 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 $\mu=\mu_{crit}\approx 1.005$ such that for $\mu\geq\mu_{crit}$ the zero-noise limit of the invariant measure is supported on the W deterministic attractor, while for $\mu\leq\mu_{crit}$ the weak noise limit of the invariant measure is supported on the SB attractor. The asymptotic state corresponds 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 $\mu$), a specific asymptotic state. Changing the value of the control parameter, one will find one or more abrupt transitions in the statistical properties of the system (here realised at $\mu_{crit}$), i.e., in other terms discontinuities in the response of the system to changes in the control parameter. What happens in our model at $\mu=\mu_{crit}$ is mathematically analogous of a first-order phase transition occurring in a near-equilibrium statistical mechanical system. We remark that, since the escape time away from either attractor grows exponentially with the inverse of the parameter controlling the variance of the noise, an individual trajectory might be trapped for very long time in a metastable state.

In collaboration with C. Kilic (Bern) and F. Lunkeit (Hamburg) the authors have started some preliminary simulations where stochastic forcing is added to PLASIM Fraedrich2005b , a much more complex climate model than the one used in this study (yet missing some essential ocean dynamical processes). As shown in Fig. 10, the first results we have obtained are encouraging in indicating that the 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 configurations leading to the noise-induced transitions, taking the thermodynamic lens we originally explored in Lucarini2010a . An important question of specific relevance for paleoclimatic and planetary science studies is to understand whether the third climatic state found in Lucarini2017 for $\mu=1.045$ is recovered also for more realistic model configurations. Additionally, the complexity of the dynamical landscape of the climate system discussed in brunetti2019 suggests the existence of a possibly topologically non-trivial network 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 explaining the ultralow frequency variability of the climate. The computational needs of a naive approach to these issues seem prohibitive, so that one should definitely take advantage of rare events algorithms to construct instantonic trajectories rubino2009rare ; Grafke_2013 ; Ragone2017 .

The approach presented here, based upon combining the knowledge of the dynamical landscape of a multistable deterministic dynamical system with the analysis of the impacts of stochastic perturbations, seems of more general interested than the specific problem we have studied. Along these lines, in Appendix A we speculate on the the possible relevance of M states in the context of the theory of biological evolution and of synthetic models of evolution. Specifically, the idea is that the presence of qualitatively diverse historical paths of evolution might result from the the presence of M states in suitable defined dynamical landscapes, exactly because M states hinder predictability 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 of understanding tipping points in the Earth system Lenton2008 as well as providing insights to a large class of multistable systems Feudel2018 . We will dedicate future efforts exactly to exploring these research lines, in particularly looking at the Atlantic meridional overturning circulation, an element of the global 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, and rather dramatic ones for the regional climate of the North Atlantic sector; see Kuhlbrodt2007 for a review of this topic. Multistability has been recently reported in the von Karman turbulent flow in Faranda2017 : the methods proposed in this paper could elucidate paths and mechanisms underlying the transitions between 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 the EU 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 received from the Institute for Basic Science (IBS), Republic of Korea, under IBS-R028-D1. The authors wish to thank T. Tél for inspiring exchanges. VL wishes to thank S. Brusatte for suggesting some fascinating literature on evolutionary biology; M. Ghil and J. Yorke for great encouragement in the course of several cheerful discussions; T. Kuna, T. Lelievre, X.-M. Li, G. Pavliotis, and N. Zagli for providing insights on stochastic dynamics. This paper has been partly prepared during the Advanced Workshop on Nonequilibrium Systems in Physics, Geosciences, and Life Sciences held at ICTP, Trieste, on May 15-24 2018 and during the Scientific Programme on the Mathematics of Climate and the Environment held at the Institut H. Poincaré, Paris on September-December 2019. The authors wish to thank F. Lunkeit for kindly sharing Fig. 10. The authors wish to dedicate this paper 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 framework of possible interest for thinking at evolutive processes in biology and for constructing synthetic models of evolution. In the speculative discussion below, specifically, we want to highlight the role of M states in making possible the existence of multiple possible yet vastly different paths of historical development of biological systems.

In the influential book discussing the Cambrian fossils from the Burgess shale and their importance for explaining the mechanisms of evolution, Gould G1989 presents some ideas on the scientific methodology inherent to historical sciences, and, specifically, to evolutionary biology. He clarifies that historical scientific explanations take the form of a narrative, whereby subsequent phenomena follow in a specific order999The related storyline approach is being currently proposed for studying weather and climate phenomena Shepherd2018 .: the time-ordering and causal links between such phenomena can be discovered and convincingly explained when multiple, independent sources provide indication for the same historical pattern of change. Gould argues that this method is fundamentally different from the classic - which he calls stereotypical - scientific method á la Galileo whereby one specific experiment can be repeated in idealised conditions, and the final outcomes of the experiment can be predicted using the fundamental laws of nature. The latter is the - indeed too simplified - version he gives of how hard sciences work. An obvious difference between the two methods comes from the fact that certain experiments - like the evolution of the climate and of the biosphere of a specific planet - cannot be repeated. Another difference comes from the fact - Gould argues - that historical processes are dominated by contingency: while we observe a specific path of historical development, many others of such paths are indeed compatible with the laws of nature, and could have been realised had the system been forced in a slightly different way in the past. We are able to understand the mechanisms of historical development, but not to predict accurately the specific path. Gould argues that the evolution in our planet could have led to fundamentally different forms of life, had the contingencies been different in the distant past. The existence of our own species is the - a priori very unlikely - result of such contingencies. In general, in systems dominated by contingency, if we could run again the movie the outcome would be vastly different. Gould’s views on evolution have been criticised by authors, as Conway Morris, proposing that convergence, 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 Morris1998 ; see the debate in MorrisGould1998 . Some authors have proposed that both mechanisms are indeed in action, yet dominant at different scales of diversification of the organisms Losos2017 .

We argue that Gould’s view can be put in the context of the mathematical framework discussed in this paper. Stochastically forced complex systems evolve in a phase space where an individual trajectory corresponds to a historical realisation of the system. The historical realisations, even starting from the same initial condition, can be vastly - even qualitatively - different if the deterministic dynamics supports the existence of multiple attractors, because different realisations can 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 the mechanism defining the evolution of the trajectory are defined by the differential equations. The predictability of the system, both of the first and second kind in the sense of Lorenz, is finite but non-zero. Finally, we can interpret the M states as the true agents of the contingency discussed by Gould. It is not the presence of a forcing, however strong, that changes radically the future history of 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 different future histories. What Conway Morris proposes is associated with a scenario where M states are either absent - so that the system is not multistable - or the noise is so weak that the probability of getting close to an M state is exceedingly low, despite the presence of stochastic perturbations, the system will be (almost) always quite close to a specific deterministic attractor, in correspondence with which basin of attraction the initial condition belongs to. In other terms. even if we could run again the movie, i.e., run another simulation, the outcome would be very similar.

The interpretation given above receives some supports from recent results obtained on synthetic models of evolution. Using the the Tangled Nature model, which is conjectured to be a prototypical example of evolution, evolutive processes are interpreted as orbits of stochastic systems in a complex dynamical landscape featuring two or more competing metastable states Jones2010 ; Jensen2018 . This viewpoint, which has been openly inspired by Gould’s ideas, is in close correspondence with what has been discussed in this paper. While in these works the language and methodology are eminently of statistical mechanical nature and they aim at detecting and studying the metastables states, the viewpoint proposed in this paper has a stronger emphasis on understanding the transitions paths between such competing attractors through the M states.