Edge States in the Climate System: Exploring Global Instabilities and Critical Transitions

Multistability is a ubiquitous feature in systems of geophysical relevance and provides key challenges for our ability to predict a system's response to perturbations. Near critical transitions small causes can lead to large effects and - for all practical purposes - irreversible changes in the properties of the system. The Earth climate is multistable: present astronomical and astrophysical conditions support two stable regimes, the warm climate we live in, and a snowball climate, characterized by global glaciation. We first provide an overview of methods and ideas relevant for studying the climate response to forcings and focus on the properties of critical transitions. Following an idea developed by Eckhardt and co. for the investigation of multistable turbulent flows, we study the global instability giving rise to the snowball/warm multistability in the climate system by identifying the climatic edge state, a saddle embedded in the boundary between the two basins of attraction of the stable climates. The edge state attracts initial conditions belonging to such a boundary and is the gate facilitating noise-induced transitions between competing attractors. We use a simplified yet Earth-like climate model constructed by coupling a primitive equations model of the atmosphere with a simple diffusive ocean. We refer to the climatic edge states as Melancholia states. We study their dynamics, their symmetry properties, and we follow a complex set of bifurcations. We find situations where the Melancholia state has chaotic dynamics. In these cases, the basin boundary between the two basins of attraction is a strange geometric set with a nearly zero codimension, and relate this feature to the time scale separation between instabilities occurring on weather and climatic time scales. We also discover a new stable climatic state characterized by non-trivial symmetry properties.


May 2018
Abstract. Multistability is a ubiquitous feature in systems of geophysical relevance and provides key challenges for our ability to predict a system's response to perturbations. Near critical transitions small causes can lead to large effects andfor all practical purposes -irreversible changes in the properties of the system. As is well known, the Earth climate is multistable: present astronomical and astrophysical conditions support two stable regimes, the warm climate we live in, and a snowball climate, characterized by global glaciation. We first provide an overview of methods and ideas relevant for studying the climate response to forcings and focus on the properties of critical transitions in the context of both stochastic and deterministic dynamics, and assess strengths and weaknesses of simplified approaches to the problem. Following an idea developed by Eckhardt and collaborators for the investigation of multistable turbulent fluid dynamical systems, we study the global instability giving rise to the snowball/warm multistability in the climate system by identifying the climatic edge state, a saddle embedded in the boundary between the two basins of attraction of the stable climates. The edge state attracts initial conditions belonging to such a boundary and, while being defined by the deterministic dynamics, is the gate facilitating noise-induced transitions between competing attractors. We use a simplified yet Earth-like intermediate complexity climate model constructed by coupling a primitive equations model of the atmosphere with a simple diffusive ocean. We refer to the climatic edge states as Melancholia states and provide an extensive analysis of their features. We study their dynamics, their symmetry properties, and we follow a complex set of bifurcations. We find situations where the Melancholia state has chaotic dynamics. In these cases, we have that the basin boundary between the two basins of attraction is a strange geometric set with a nearly zero codimension, and relate this feature to the time scale separation between instabilities occurring on weather and climatic time scales. We also discover a new stable climatic state that is similar to a Melancholia state and is characterized by non-trivial symmetry properties.

Introductory Remarks
Providing a comprehensive framework able to bring together the understanding of climate dynamics, able to encompass its variability and its responses to forcings, is one of the grand challenges of contemporary science [1,2,3,4]. While classical approaches to the definition of climate sensitivity are based on concepts borrowed from the analysis of systems at or near equilibrium, a more general point of view is indeed needed [5,6]. The issue has relevance both for the understanding of a) the history of our planet, the relationship between climatic conditions and the development of terrestrial life across the geological epochs, and the pressing issue of anthropogenic climate change; and b) for the possibility of having an encompassing view of planetary circulations, a much coveted goal in view of the extraordinary achievements of contemporary planetary science. Many planets are being discovered by the day, and the quest for discovering atmospheres able to support life (as we know it) is ongoing; see the NASA exoplanet archive (http://exoplanetarchive.ipac.caltech.edu) for some useful information at this regard.
The climate is a complex system featuring variability of a vast range of spatial and temporal scales, resulting from instabilities fundamentally fuelled by the inhomogeneous absorption of solar radiation, and from the corresponding fluxes of mass, chemical species, and energy that, acting as negative feedbacks, allow for the establishment of steady state conditions. Interestingly, one can interpret the presence of organized motions of the geophysical fluids as a result of mechanical work, transforming available potential energy into kinetic energy, which is eventually dissipated by friction [7,8]. Altogether, the climate can be seen as a thermal engine able to transform heat into mechanical energy with a given efficiency, and featuring many different irreversible processes that make it non-ideal [9,10,3]. An overarching theory of climate dynamics, able to encompass instabilities, re-equilibrating feedbacks, and large scale dynamical structures is still far from being achieved, despite lots of progress in the last decades.
Nonequilibrium systems like the climate can be essentially characterized as being in contact with at least two thermostats with different temperatures (or, e.g. chemical potentials) [11]. Nonequilibrium systems exhibit an extremely complex phenomenology: while a lot is known for (near)equilibrium systems, our understanding of nonequilibrium systems is comparatively poor and limited, in spite of the wealth of phenomena occurring out of equilibrium.
Obviously, much of the interest on the understanding of the response of the climate system to forcings comes from the evidence of the complex interplay between climate conditions and life. Since the early 1990s, the Intergovernmental Panel on Climate Change (IPCC) has started collecting and analyzing in a systematic way the scientific literature on the topic of climate dynamics, looking at the distant and near paleoclimatological conditions, but with a special focus on global warming and on the socio-economic impacts of anthropogenic climate change [12,13,14]. An enormous scientific effort is aimed at improving our understanding of climate change, and climate-related investigations account for some of the most extensive data collection campaigns and computational exercises across all disciplines.
The response of the climate system to perturbations can be qualitatively divided in two different yet coexisting forms. In some cases, such a response is smooth with respect to the perturbation parameters. In other cases, the smoothness is lost and abrupt changes can result from small perturbations. In some previous works, we have studied using a statistical mechanical approach the regime of smooth response [3,15,4,16] and the conditions under which such smoothness is lost [17]. In other previous works, we have, instead, used thermodynamical concepts to better understand under which conditions small perturbations can lead to qualitative changes in the properties of the climate system [18,19,20]. In this work we wish to further advance the understanding of the global instabilities of the climate system by taking advantage of some powerful concepts and methods borrowed from the theory of dynamical systems.
In Sect. 2 we review the two scenarios of smooth vs. nonsmooth climatic response to perturbations, in order to provide a meaningful context for the results contained in this paper and for the motivations behind our work. We introduce the concept of critical transitions (tipping points) in geosciences [21,22], and provide basic concepts and models for understanding the multistability of the Earth climate coming from the coexistence of the so-called snowball and warm climate states in a wide range of boundary conditions [23,24,25,26,27,28].
In Sect. 3 we present an overview of some methods discussed in the literature for approaching the problem of investigating the critical transitions. We briefly summarize the rationale behind using stochastic models for studying climate dynamics [29,30,27] and discuss how large deviations theory [31] can be relevant, in the context of the Freidlin-Wentzell theory [32] for explaining the transitions of the system between coexisting regimes [33,34]. Following the work of Eckhardt and collaborators in the context of the fluid dynamics of turbulent flows, and sticking to the paradigm of deterministic dynamics, we introduce the concept of edge state, the unstable saddle separating different regimes of motions and discuss how it is possible to find it using the edge tracking algorithm [35,36,37]. We also briefly discuss the application we have previously presented [38] on a simple geophysical model [25].
In Sect. 4 we introduce PUMA-GS, a new simple yet potentially relevant climate model constructed by coupling a modified version of the model introduced in [25] with PUMA [39], a simple open-source primitive equations model of the atmospheric circulation. We then clarify how using a suitably adapted version of the edge tracking algorithm it is possible to identify the edge states of the system, i.e. unstable climatic states living on the boundary between the basins of attraction of the snowball and warm climate states ‡. ‡ We have decided to use for the edge states of the climate system the name of Melancholia states, as a homage to the 2011 movie Melancholia directed by L. von Trier, where a masterful and accurate portrait of the challenges of prediction in the vicinity of critical transitions is given. The movie itself has provided a fundamental inspiration for this investigation. Just to give a -possibly In Sect. 5 we study the geometrical properties of the boundary between competing basins of attraction, along the lines of [40,41,42,43,37], the properties of the climatic edge states, investigating whether they correspond to fixed points, periodic points or chaotic saddles featuring a limited horizon of predictability, and the bifurcations between such regimes, including a process of symmetry breaking. We also show how the methodology proposed here allows for identifying a new stable climate regime that would hardly be discoverable through standard direct numerical integration.
In Sect. 6 we summarize the main ideas presented in this paper and present our conclusions and perspectives for future work.

Response Theory and Climate Change
In (near-)equilibrium systems, we have formidable tools for relating forced and free fluctuations of a system. Kubo [44,45] constructed a mathematically nonrigorous but physically extremely powerful theory able to predict how a system whose statistics is described by the canonical ensemble responds to external perturbations. Kubo's response theory also includes the fluctuation-dissipation theorem, which basically shows that the linear response of a system to forcings can be derived from its natural fluctuation through a suitably defined linear operator. This result has been extremely influential in many areas of the physical and chemical sciences. See a comprehensive review in [46] and a recent extension to the nonlinear case in [47].
The lack of mathematical rigour of the response theory was later addressed by Ruelle [48,49,50], who showed that if one considers Axiom A dynamical systems, and limits oneself to considering sufficiently smooth observables, it is possible to prove the differentiability of the invariant measure of the system with respect to small parameters describing the change in the dynamics, and to provide explicit formulas for computing the change in the measure. Axiom A systems are chaotic and uniformly hyperbolic on the attractor, where the asymptotic dynamics takes place.
Moreover, modern methods of spectral theory have provided different and elegant proofs of Ruelle's results. The response theory can be developed by comparing the Perron-Frobenius transfer operator [51] of the unperturbed and of the perturbed system. The spectrum of the transfer operator describes, among other things, how an initial probability distribution converges to -in the case of mixing systems -a unique invariant measure. In this way, one can study directly the change in the invariant measure resulting from the the applied perturbation [52]. Such a point of view has also allowed the extension of Ruelle's results to the study of the response for smooth observables in more general dynamical systems [53,54] or of non-smooth observables in Axiom A systems [55].
Axiom A systems are indeed far from being typical dynamical systems, but, according to the chaotic hypothesis of Gallavotti and Cohen [56], they can be taken as effective models for chaotic physical systems with many degrees of freedom. Specifically, this means that when looking at macroscopic observables in sufficiently chaotic (to be intended in a qualitative sense) high-dimensional systems, it is extremely hard to distinguish their properties from those of an Axiom A system, including some degree of structural stability. One can interpret the chaotic hypothesis as the possibility of constructing robust physical properties for the system under investigation. Therefore, providing results for Axiom A systems can be thought of as being of rather general physical relevance.
Axiom A systems corresponding to equilibrium physical systems possess an invariant measure that is absolutely continuous with respect to Lebesgue because the phase space does not contract nor expand, as the flow is nondivergent. Indeed, in this case the fluctuation-dissipation theorem applies.
Axiom A systems corresponding to nonequilibrium conditions feature a phase space that contracts on the average, so that the asymptotic dynamics takes place on a strange attractor with dimension smaller than the dimension of the phase space. As a result of the singularity of the invariant measure, the usual form of the fluctuation dissipation theorem does not hold. Ruelle [48,49,50] provides a mathematical explanation of this property, while a physical interpretation is given in e.g. [57,58,59]. This implies that for nonequilibrium systems there is no obvious relationship between free and forced fluctuations of the system, as already suggested by Lorenz [60].
Various proposals have been put forward to avoid the problem of the lack of smoothness of the invariant measure in nonequilibrium statistical mechanical systems. Some authors consider adding some stochastic forcing on top of the deterministic dynamics, in order to take into account the effect of unresolved scales [61], having in mind the Mori-Zwanzig theory [62,63].
Interestingly, while on one side there have been positive examples of applications of the fluctuation-dissipation theorem for studying the climate, it is clear that the quality of the obtained response operator depends substantially on the chosen observable [64,65,66]. The difficulties in constructing ab-initio the response operator using Ruelle's explicit formulas come from the extremely different mathematical properties of the contributions coming from the unstable and stable manifold [61] plus from the possible presence, in the case of high-dimensional chaos, of rather diverse time scales. Algorithms based on adjoint methods seem to partially ease these issues [67,68]. A possible way to sort out the problem of estimating the response operator relies on projecting the perturbation flow on the covariant Lyapunov vectors, which allow for constructing a covariant base spanning the stable and unstable manifolds [69], or exploiting the reconstruction of the invariant measure obtained using unstable periodic orbits [70,71]. In fact, through a combined use of the formalism of covariant Lyapunov vectors and of unstable periodic orbits, it has been recently possible to elucidate interesting mathematical and physical aspects behind the violation of the fluctuation-dissipation theorem in a simple atmospheric model [16].
Convincingly good results in terms of climate prediction performed using the linear response theory have instead been obtained by bypassing the problem of constructing the response operator and using, instead, the formal properties of the Green's function [59,3,15]. Note that the response theory is able to predict not only the response of globally averaged quantities, but also the spatial patterns of the change [4].

Where the Response Theory Fails: Critical Transitions of the Climate System
The response theory is based on a perturbative expansion of the invariant measure, so its intrinsic limitation comes from the fact that the radius within which the expansion is valid could be extremely small. Of course, one might think of covering an extended range of parametric changes by repeating the expansion in neighbouring intervals. However, the procedure will fail in the vicinity of critical transitions, which correspond to qualitative change in the dynamics of the system resulting from a crisis of the attractor.
Bifurcations involving non-chaotic attractors have been very accurately studied, usually taking the point of view of analysing how the existence and the stability properties of invariant sets of the considered system change when one parameter is changed [72]. The resulting bifurcations can be classified according to universal classes [73], which has proved of immense relevance for studying a large variety of physical systems, and indeed successfully in the context of climate dynamics [1,74].
We remark that Ashwin et al. [75] have recently emphasized the need for widening the usual treatment of bifurcations in dynamical systems in order to accommodate for more general scenarios of critical transitions. In addition to the usual bifurcation scenario described above (B-tipping), they have studied the transitions between competing basins of attraction depending on the presence of a finite rate of change of one or more parameters of the system (R-tipping), or on the presence of stochastic forcing of finite amplitude (N-tipping). See in [76,77] an earlier discussion of the role of time-dependent forcings in the context of another geophysically relevant critical transition, i.e., the collapse of the large scale oceanic circulation [78].
The situation is much more complex in the case of high-dimensional and/or chaotic attractors, where critical transitions are associated to a large variety of mathematical mechanisms, including boundary crises, where an attractor is destroyed by contact with its basin; interior crises, where the shape of the attractor undergoes a sudden change; synchronization and symmetry breaking (and restoring) processes, whereby different subsystems exhibit (or lose) special temporal and/or spatial correlations; and many others; see comprehensive treatments of these fascinating issues in, e.g., [79], and a survey of recent advances on synchronization phenomena (which have been initially studied in the simpler context of non-chaotic attractors) in [80].
For reasons that will be clear in Sec. 3, in the rest of the paper our main (but not exclusive) focus will be on the boundary crises, which are accompanied by the disappearance of an attractor. The vicinity of such critical transitions is accompanied, on the one hand, by a rough dependence of the system's properties on the parameter to be tuned, as a result of the Ruelle-Pollicott resonances, and, on the other hand, by the presence of a slow decay of correlations (critical slowing-down) of the system's observables, with the decorrelation times diverging when the control parameter reaches its critical value [81,82,17]. Looking at the Ruelle response formulas, it is clear that the lack of a slow decay of correlations leads to a lack of convergence in the estimate of the Green's function, even though the relationship between the presence of sufficiently fast decay of correlations and the validity of response theory can be clarified more easily using spectral theory [83].
The climate system is characterized by such critical transitions, which in the geophysical literature are usually dubbed as tipping points. Relevant examples of tipping points include the dieback of the Amazon forest, the shut-down of the thermohaline circulation of the Atlantic ocean, the methane release resulting from the melting of the permafrost, and the collapse of the atmospheric circulation regime associated to Indian monsoon; see [21,22] for a useful perspective. Obviously, not all tipping points are equally critical, both in terms of global climatic impacts, and, more technically, in the sense that the lack of smoothness in the response might be practically detectable only when looking at specific observables. As an example, one can expect that nonsmoothness in the response associated to the dieback of the Amazon forest will be easier to detect when looking at climatic fields in South America rather than in Siberia or Australia. Understanding the qualitative and quantitative properties of climatic tipping points is another great challenge of climate science.

Snowball and Warm Climate States
Probably the most impressive example of critical transition in the Earth system can be found when considering that the current astronomical and astrophysical conditions support at least two possible steady climate states, one being the one we are experiencing, and the other one characterized by much lower temperatures, almost total absence of water vapour in the atmosphere, and global ice-cover, referred to as snowball Earth. Models of different levels of complexity support such a statement, and, more importantly, geological evidence suggests that during the Neoproterozoic (the period spanning from about 1000 million to about 540 million years ago), the Earth suffered two of its most severe periods of glaciation [26,84,85] and entered into a snowball climate state.
The onset and decay of snowball conditions are associated to changes in the value of the solar constant and in the composition of the atmosphere, with the ensuing variation of the greenhouse effect. For a fixed atmospheric composition one finds a range of values of the solar constant S * where the system features multistability: the warm climate can exist for all values of S * > S * W →SB and the snowball state can exist for all values of S * < S * SB→W , with S * W →SB < S * SB→W . The critical value S * = S * W →SB (S * = S * SB→W ) is such that if we are in the warm (snowball) state and slowly reduce (increase) the value of S * , for S * = S * W →SB (S * = S * SB→W ) the stability of the warm (snowball) climate is lost and the system tips to the snowball (warm) state.
The critical transitions associated with the boundary crisis happen in conditions where the stabilizing negative feedbacks, to be described later, are overcome by the positive ice-albedo feedback. The positive feedback works as follows: colder (warmer) temperatures lead to an increase (a decrease) in the ice cover, which, by increasing (decreasing) the albedo, leads in turn to lower (higher) temperatures [27]. When the ice-albedo feedback dominates, the temperature then decreases (increases) substantially until a new balance sets in, corresponding to the cold (warm) stable state realized after the critical transition.
We wish to remark that some theoretical investigations and observational evidence point at the possible existence of a third possible stable climate state, the so called slushball Earth. Such a state is characterized by cold conditions (yet not as extreme as in the snowball), where the equatorial belt is ice-free and characterized by vigorous oceanic and atmpsoheric circulations and a non-trivial hydrological cycle; see, e.g., [86,87,88]. In the rest of the paper, we will only comment briefly on this -indeed relevant -geophysical problem.
A simple mathematical formulation of the physical processes responsible for the snowball-warm multistability of the Earth system can be given by energy balance models (EBMs), which provide a simplified yet relevant representation of the energy exchanges in the climate system [89,27]. Reducing the Earth to a 0-dimensional (0-D) system (in physical space) where the only relevant quantity is the global temperature indicator T (which can be heuristically interpreted as a globally averaged effective temperature), one can construct an EBM by modelling the energy budget as follows [27]: where t is time, C is an effective average heat capacity per unit area, I = S * /4 is the average incoming solar radiation per unit area, S * is the solar constant (the factor 4 emerges looking at the geometry of the Earth-Sun system, see [27]), α is the albedo, which is a non-increasing function of T , and O is the outgoing radiation per unit area. O(T ) is a monotonically increasing function of T ,i.e. an increasing surface temperature leads to an increase in the outgoing radiation, which is the basic mechanism of the Boltzmann (negative) radiative feedback. By suitable (and reasonable) choices of α and O, one indeed finds bistability as a function of S * ; see the discussion in [27] and Fig. 1a, where the two stable (warm and snowball) solutions are separated by the unstable solution II.
Note that the time derivative of the temperature in Eq. 1 can be written as minus the derivative of a potential V (T ), and in the case of bistability the two local minima of V correspond to the stable solutions, and the local maximum of V corresponds to the unstable solution, see Fig. 2 for a qualitative description.
EBMs can be extended in such a way as to include a latitudinal dependence of the Earth's temperature, as in the case of the celebrated models by Budyko [23], Sellers Bifurcation diagram obtained using a 1-D EBM as described in Eq. 2, where γ controls the strength of the diffusion operator. The y-axis shows the ice-line latitude ϕ i , the control parameter in the abscissa is S * /S * 0 . Adapted from [27]. [24], and Ghil [25]. The time evolution of the temperature for one-dimensional EBMs (1D-EBMs) can be written in terms of the following partial differential equation: where, as opposed to Eq. 1, the heat capacity C, the incoming radiation I, and the albedo α are explicitly dependent on the latitude φ, and we have to include in the energy budget the divergence of energy transport performed by the geophysical fluids, Snowball Warm Figure 2: Scalar double well potential function V (T ); the warm and the snowball states correspond to the two attractors of the system, separated by the saddle point. The noise triggers the transition from one basin of attraction to the other with a mean exit time described by Eq. 7.
represented as the diffusion operator D[T, φ]. It is important to note that the diffusion operator describes in a very simple yet efficient way the action of the negative feedbacks regarding the meridional temperature gradient, through the meridional heat transport facilitated by atmospheric and oceanic motions. On weather time scales (few days), the feedback is powered by the baroclinic instability: the stronger the large scale meridional temperature gradient, the higher the intensity of the midlatitude atmospheric eddies and of their related counter-gradient heat transport [8,90]. On climatic time scales (several years) the ocean acts in a similar way (though by different specific physical mechanisms), by dampening the large scale temperature differences through the transport of heat from warm to cold regions [28]. With a suitable choice for the parameters controlling the terms responsible for the energy budget in Eq. 2, we obtain a range of values of S * where bistability is found; see The existence of multistability for the Earth's climatic conditions is a robust feature across a modelling hierarchy, ranging from the simple EBMs described above to global climate models: see a detailed analysis of acting feedbacks as well as glaciation/deglaciation scenarios and mechanisms in [28]. Global climate models feature a range of bistability [S * W →SB , S * SB→W ] that includes the present value of the solar constant S * 0 , see Fig. 3 a) for an example. While in simple EBMs the stable solutions correspond to fixed points in the phase space, in global climate models, the attracting solutions correspond to states featuring, in general, chaotic behaviour, where the § It is historically remarkable to note that, after the Budyko and Sellers models were introduced, it became apparent that a Nuclear Winter, by reducing the incoming solar radiation as a result of increased albedo due to a dramatic increases in the particulate matter in the atmosphere, could potentially trigger an even greater disaster for life on Earth than the nuclear war itself. This has been influential in reducing the size of the nuclear arsenals at the end of the Cold War. It is comforting to note that, somewhat ironically, the two contributions came almost simultaneously from scientists belonging to the two opposing geopolitical blocks. invariant measure is supported on a strange attractor [18,3]. The temperature profiles of the stable climates obtained with the 1D-EBM given in Eq. 2 with a suitable choice of the value of the parameters are in broad agreement with what obtained with complex general circulation models, see [38]. Figure 3: Multistability of the climate system corresponding to the coexistence of warm and snowball conditions. a) Globally averaged surface temperature T S . b) Average intensity of the Lorenz energy cycle C(P, K). c) Efficiency η of the climate system. W-SB (SB-W) indicates the occurrence of transitions from warm to snowball conditions (from snowball to warm conditions). These results are obtained using the open source climate model PlaSim [91]. Adapted from [18].
For a given value of S * within the range of bistability, the warm climate is characterized by a very active atmosphere whose variability is sustained by the presence of a large meridional temperature gradient, which supports the existence of eddies as result of the baroclinic conversion between potential and kinetic energy, whereas the snowball climate has a much lower level of atmospheric variability, with a correspondingly low meridional temperature gradient. Figure 3 b) shows that the intensity of the Lorenz energy cycle C(P, K), given by the average rate of conversion of potential into kinetic energy, which, at steady state, is also equal to the average rate of dissipation of kinetic energy due to friction [7,8,10], is substantially higher in the warm than in the snowball state. See [28] for a comprehensive discussion of the atmospheric circulation in the two regimes.
Thermodynamics provides additional insight to the critical transitions. One can define a Carnot-like climate efficiency η = (Θ + − Θ − )/Θ + in terms of the average temperatures inside the system where heating (Θ + ) and cooling (Θ − ) processes take place on the average, respectively [10]. The efficiency is higher when there is stronger correlation between temperature fluctuations and heating rates. The tipping point is accompanied by a sudden decrease in the efficiency η, which can be interpreted as a result of the system getting closer to thermodynamic equilibrium, see Fig. 3c) and discussion in [3,18,19,20].

Energy Landscapes, Large Deviations, and Transitions
Note that, far from the critical transitions indicated in Fig. 3, it is possible to jump from the warm attractor to the snowball attractor and vice versa by perturbing the system through suitably defined time-dependent deterministic and/or stochastic forcing. The transition from one basin of attraction to another one is most easily obtained through a Dirac's δ-like perturbation applied to the orbit (this might correspond, in physical terms, to the impact of an asteroid on our planet).
Nonetheless, arguably the most interesting scenario of forcing is given by the presence of stochastic perturbations. Climate science has indeed been one of the first areas of natural sciences where the paradigm of stochastically perturbed dynamical systems has been extensively used, following the Hasselmann programme [29]; see also later discussions in [30,27,92]. The basic idea is the following: if one considers separately slow and fast modes of variability of the climate system, in the limit of infinite time scale separation, the impact of the fast modes on the slow modes can be essentially treated as a stochastic correction to the deterministic dynamics of the slow modes. This point of view is closely related to the mathematical theory of the averaging method for performing the coarse-graining of the dynamics of multiscale dynamical systems [93,94].
The Freidlin-Wentzell theory [32] provides a comprehensive framework in multistable systems for studying the probability of transitions outside one of the basins of attraction triggered by stochastic forcing. The Freidlin-Wentzell theory is most easily understood in the case of a system whose dynamics takes place in an energy landscape A larger (smaller) concentration of greenhouse gases leads to a shift of the range of multistability towards lower (higher) values of S * . The presence of gigantic concentration of CO 2 is deemed responsible for the exit of the Earth from the snowball state experienced in the Neoproterozoic [28]. Note also that the range of multistability is altered when one considers variations in the rotation rate of the planet, with multistability being eventually lost and a unique climate emerging when very slow rotation rates are considered [20,19].
(plus additive white noise), such that: where y ∈ R n , V : R n → R is sufficiently smooth, and dW is a vector of n independent increments of Brownian motions (the Itô or Stratonovich conventions are equivalent, in this case), and determines the strength of the noise. As is well known, one can associate to the stochastic differential equation given by Eq. 3 a Fokker-Planck equation describing the evolution of the probability density function (pdf) of an ensemble of trajectories obeying the stochastic differential equation as follows: where, under suitable conditions of integrability and in the weak noise limit → 0, the stationary solution corresponding to the invariant measure is given by lim t→∞ p(y, t) = p(y), where at leading order p(y) ∝ exp[−2V (y)/ 2 ]. The local minima of the potential V (fixed-point attractors in the case of deterministic dynamics with = 0) correspond to the local maxima of p, so that, e.g., a double-well potential corresponds to a bimodal pdf. In the limit of weak noise, trajectories starting near a local minimum of V typically a wait long time before moving to the neighborhood of a different local minimum of V , and the transitions take place most likely through the lowest energy saddle linking the initial basin of attraction to any other basin. A comprehensive treatment of the critical transitions accounting for the effect of noise can be found in [95]. Let's now consider the simple case of y ∈ R and focus again on the problem of the transitions between the snowball and warm climates. If one includes stochastic forcing in the form of additive white noise on the right hand side of Eq. 1, we have where dW is the increment of a Brownian motion and > 0. In the limit of weak noise, the stationary distribution can be expressed at leading order as Using large deviations theory, one can derive that the mean exit time corresponding to the transition from the basin of attraction of the stable solution A to the basin of attraction of the stable solution B through the unstable saddle U can be written in first approximation as: while the constant in Eq. 7 is specified by the celebrated Kramer's escape rate formula [96,97]. See Fig. 2, where the states A and B can be though of as corresponding to the warm and snowball state, and the saddle U coincides with the boundary of the two basins of attraction.
In the geophysical literature, tipping points are mostly studied by using Eq. 6 to construct an equivalent one-dimensional pseudo-potential from the pdf derived starting from the time series of a selected climate observable. Subsequently, Kramer's theory is used to estimate the probability of transition from one basin of attraction to another one within a given time frame [21,22]. Nonetheless, inconsistencies in this method emerge from the fact that the operation of projecting the dynamics onto only one observable is very likely to cause errors in the evaluation of the time scales of the dynamics, thus breaking the simple and powerful relations existing between invariant density and escape probability from a local mimimum of the effective potential, see the discussion in [98].
When considering dynamical systems featuring multiple attractors and undergoing gaussian stochastic forcing, large deviations theory [31] provides methods for constructing the so-called instantons, which are the most likely noise-driven trajectories leading to the transitions from one basin of attraction to the other one. An instanton is constructed as a minimizer of the related Freidlin-Wentzell action and can be interpreted as the most efficient way to exit a local potential minimum. See [99] for an introduction in the context of fluid dynamics.
Some geophysical fluid dynamical systems feature extremely rare transitions between different regimes of motions. These are rare excursions between regions of the attractor of the system that take place only on ultralong time scales, so that, performing an operation of coarse graining to the dynamics, one can describe them as noise-induced transitions between separate attracting sets for the deterministic part of the dynamics [33,34].
The instantons method allows for treating more general cases than those of systems whose dynamics is governed by an energy landscape. What we find a bit unsatisfactory in this otherwise extremely powerful framework of noise-driven transitions is the fact that one is left with the question of what is the mathematical nature and the physical justification of the noise. Note that, given the lack of time-scale separation in the climate variability [8], it is a bit of a stretch to use standard arguments to justify the presence of white noise as resulting from fast atmospheric motions. In other terms, despite the huge insight of the Hasselmann program mentioned above, climate dynamics is in fact not the best setting for using the averaging method. As discussed in [100,101] in the spirit of [62,63] and confirmed in the mathematically more rigorous treatment in [102,103], when no scale separation is present between the scales of motions we want to resolve and those we want to parameterize, the effect of the latter on the former entails considering time-correlated noise and non-markovian effects.
As a final note we wish to remark that, even in the context of systems of the form given in Eq 3, if one assumes that the noise possesses instead slow decay of correlations, such as Lévy processes, the standard Freidlin-Wentzell picture is dramatically altered, with an entirely different dynamical scenario to be considered. While in the case of white noise the transition from one basin of attraction to the other occurs in the very unlikely case that many subsequent stochastic perturbations conjure to push the orbit across the saddle, in the case of Lévy noise one or few large kicks do the job [104,105], so that the expression for mean exit time is rather different from what reported in Eq. 7. This scenario seems to be relevant in a climate context as well [106,107,108].

Edge States
Following a classical point of view within statistical mechanics, we would like to be able to understand how noise-induced transitions between separate climatic attractors may take place by gaining insight on the underlying purely deterministic evolution laws. In this regard, we take advantage of the edge tracking method recently developed by Eckhardt and collaborators for studying fluid dynamical systems possessing multiple (quasi-)steady states [35,109,36].
The plane Couette flow or the pipe flow, in the case of sufficiently high Reynolds numbers, feature the coexistence of a (quasi-)attractor corresponding to the turbulent state and of an attracting fixed point corresponding to laminar flow. We remark that the turbulent state cannot be properly described as corresponding to a separate attractor, because turbulence is transient -even if with exceedingly long life time, which grows superexponentially fast with the Reynolds number [110,111,112,109]. But, using arguments of time-scale separation, we can say that the scenario is (almost) indistinguishable from that of two separate attractors, where the turbulent state is steady, so we will use an abuse of terminology in describing these results, and refer to the turbulent state as corresponding to a proper attractor, possessing its own basin of attraction.
The basic goal is to understand what lies in-between the two attractors in order to a) have a global view on the properties of the phase space of the system, beyond the paradigm of looking at steady states; and b) understanding which combinations of external forcings and internal processes might more easily lead to transitions between the steady states.
Clearly, in the case of a bistable system, there must be a separatrix, a repelling set which acts as boundary between the two basins of attraction. The natural question one might ask is: what is the evolution of an orbit starting nearby or on such a separatrix? In the case of an energy landscape, such a basin boundary is the mountain crest separating the two basins of attraction, and the evolution law would lead an orbit starting from a point on the boundary to the nearby local minimum restricted to the subset of the phase space corresponding to the basin boundary), where the evolution would stop. Clearly, the end point is in fact a saddle, because there is an (additional) unstable direction. Orbits departing near the basin boundary but not exactly on it, would end up in the attracting set corresponding to the basin they start from.
Eckhardt and collaborators showed that trajectories initialized on the boundary between the two basins of attraction converge to a unstable saddle, which is a relative attractor with respect to the basin boundary (but obviously a repelling set when viewed globally). The scenario has been confirmed in the case of truly separate attractors in a toy model discussed in [37].
Such an unstable saddle, the so-called edge state, has been found to be either a fixed point, or a periodic orbit, or, in some cases, a chaotic solution taking place on a strange geometrical set. The relevance of the edge state lies in the fact that it is crucial a) b) Figure 4: Edge state in the pipe flow. a) Time evolution of the energy of the turbulent flow for two trajectories starting near the edge state but belonging to two different basins of attraction: one trajectory collapses to zero (laminar state), and another one increases its turbulent kinetic energy until the quasi-steady turbulent state is reached. b) Turbulent kinetic energy of the steady turbulent state (square dots) and of the edge state (triangles) for various values of the Reynolds number. The turbulent energy of the other steady state, the laminar state, is zero by definition. Note that we are using an abuse of language in defining the turbulent state as corresponding to a true attractor (it is a very long transient, instead, see text). Adapted from [36].
for understanding the global properties of the system and that it provides the gate between the two attractors: trajectories resulting from (weak) external forcings that are successful in achieving the critical transition have to pass nearby the edge state.
Of course, given its (extreme) instability, one cannot find the edge state by direct forward numerical integration; instead, one can find it by a suitable algorithmic procedure named edge tracking. The basic idea is to find by bisections (and long forward integrations) two nearby points in the phase space that are on the two sides of the separatrix, follow their forward evolution for a given time (typically short), and then repeat the bisection in order to reduce the divergence between the two trajectories due to the global instability connected with the bistability of the system. The bisection is obtained by interpolation to the midpoint (more general constructions of convex combinations are possible) between the initial conditions belonging to the two sides of the basin boundary, checking, using a long run, whether the new initial condition leads to one or to the other attracting set, and repeating until two nearby (according to a prescribed criterion) initial conditions leading to different asymptotic dynamics are obtained. One retains short segments of the pair of control trajectories belonging to these most tightly bracketing initial conditions. The length of these segments are limited by a requirement on the maximal separation of the control trajectories. After a transient, any of the segments obtained iterating this procedure can be taken to represent the edge state, approximated by many discontinuous segments of the trajectory.
A key mathematical fact to be considered here is that the basins of attraction are invariant sets, which ensures that the procedure is well-defined. It is a crucial -in terms of robustness and efficiency -aspect of the procedure to choose a useful observable able to tell us unambiguously and soon enough whether the trajectory is ending up on one or on the other attractor, see Fig. 4a. An efficient observable has a relatively small variability in the transient runs compared with the actual difference between the expectation values computed on the two separate attractors.
Additionally, it seems quite natural to investigate the geometrical properties of the boundary separating the the basins of attraction, as it is not a-priori clear whether one should expect a manifold of co-dimension one or a more complex geometrical object. We will come back to this aspect in Sect. 5.3.
In the case of the problems studied by Eckhardt and collaborators, the natural choice is to use the total turbulent kinetic energy of the flow. In the limit of infinitely small separation between the two bracketing trajectories, one actually obtains the transient behaviour leading to the the edge state, and, then, is able to reconstruct the dynamics on the edge state. While this procedure is rather intuitive, one might reasonably expect lots of technical difficulties related to the need of controlling strong instabilities in a high-dimensional dynamical system. Eckhardt and collaborators have convincingly shown the robustness of the edge tracking method [35,36,37].
In a previous paper [38] we have introduced in the geophysical community the edge tracking method to provide a fresh outlook on the properties of the Ghil-Sellers (GS) 1D-EBM) [25], which can be written in the form of Eq. 2. Ghil [25] solved the appropriate boundary value problem and found for a wide range of values of the solar constant three coexisting stationary solutions, and carried out a stability analysis able to identify the stable warm and snowball climates, plus an unstable state. Such a state was identified by Ghil as the saddle point of a suitably constructed potential, and we proved that it coincides with the edge state found using the tracking procedure, see Fig. 5. Our investigation made it apparent that the most efficient (but definitely not unique) observable to be used for characterizing during the edge tracking whether  [25] model defining the constitutive relations between the normalized intensity of the incoming radiation (µ = S * /S * 0 = 1 corresponds to the current conditions), the globally averaged temperature [T ] and the low-to-high latitudes temperature difference ∆T . The thick red, green, and blue lines correspond to the warm, edge, and snowball state, respectively. The thin blue lines depict heteroclinic trajectories -the instantons. Adapted from [38].
the trajectory ends up in the snowball or warm state is the globally averaged surface temperature, because warming and cooling are mainly controlled by the anomalies in radiation emission and changes in albedo, which are both to a good extent controlled by changes in the globally averaged temperature. We also studied the physics of the three states, relating their instabilities to relevant macroscopic thermodynamical properties such as large scale temperature gradients and entropy production, using the conceptual framework introduced in [18,19,20].

Remark
Looking at the top right projection of the three-dimensional curve portrayed in Fig. 5, one sees that the edge state is characterized by the fact that lower average temperatures correspond to higher values of solar constant. One can interpret this feature as representative of the fact that the edge state has, effectively, a negative heat capacity, and is not thermodynamically stable. Making a Gedankenexperiment where the Earth system in the edge state is put in an isolated box with another body emitting and absorbing radiation and looking at the dynamics of fluctuations should clarify this point.
Interestingly, a qualitatively similar property can be found when looking at the case of turbulent flows analysed by Eckhardt and co. [36] and portrayed in Fig. 4 b): in the edge state, larger values of the Reynolds number, which is representative of the applied pressure gradient, correspond to lower values of the turbulent kinetic energy of the flow. This can be interpreted as corresponding to conditions of mechanical instability, where the viscosity is negative.

Melancholia States in a Climate Model
In this paper we want to show how the concept of the edge state can be used for characterizing globally the dynamical properties of a severely simplified yet Earth-like new climate model featuring multistable properties in a realistic range of values of the solar constant. As mentioned in Sect. 1, we propose to call the climatic edge states Melancholia states. Figure 6 clarifies how this work relates to previous analyses we have performed on the thermodynamics and statistical mechanics of climate. While most of the studies we consider here have been performed using the open source climate model PlaSim [91], in this specific investigation, given the complexity of the procedure of edge tracking, we resort to using a simpler model, later described in detail in Sect. 4.
The climate model PUMA-GS is constructed by coupling PUMA, an open source three-dimensional primitive equations spectral model of the atmosphere [39], with the surface described by a slightly modified version of the reaction-diffusion GS model [25], extended symmetrically along the longitudinal direction, thus making it possible to couple it with the atmosphere aloft. The coupling is realized by imposing that the relaxation temperature profile that provides the baroclinic forcing to the atmospheric component is enslaved, through a simple representation of radiative convective adjustment, to the surface temperature field of the GS model, which basically plays the role of the ocean. This provides a natural and coherent separation between fast variables (atmosphere) and slow variables (ocean). The ocean exchanges energy with the atmosphere aloft as well as featuring absorption of shortwave radiation and emission of longwave radiation. The model features multistability for reasonable values of its parameters and is described in detail later in the paper. This newly introduced model sits in terms of complexity between PUMA [39] and PlaSim [91].
PlaSim is comparatively simple with respect to the state-of-the-art climate models considered in the compilation of the latest IPCC report [14], but nonetheless includes a (simplified) treatment via parameterizations of important physical processes such as convection, sea ice formation, clouds formation, precipitation in liquid and solid forms, boundary layer processes, which are entirely absent in PUMA-GS, plus a much more advanced treatment of radiative processes. Additionally, PlaSim allows for a flexible configuration of the land surface mask and of the orography, with an ensuing large variety of configurations for the ocean. Nonetheless, the ocean component of the model, similarly to PUMA-GS, has merely the role of horizontally transporting heat via diffusion, while not featuring any description of actual ocean currents.
The challenge we address is threefold:  Figure 6: In [18,113,19,20] we have studied the macroscopic thermodynamic properties of the climate system in the W state and SB state. In [15,4] we have studied the differentiability of the invariant measure in the W state around present climate conditions (purple shading, see DIFF). In [17] we have studied the loss of differentiability of the invariant measure near the W → SB tipping point (orange shading, see NonDIFF). All of these analyses have been performed using the climate model PlaSim [91]. The goal of this paper is to charactacterize what is in-between the W and SB states, i.e., constructing and studying the edge states (green shading, see EDGE STATE). For technical reasons, we use a simplified yet Earth-like climate model, PUMA-GS. The 3D surface plot is adapted from [19].
• We want to show that the edge tracking algorithm works also in the context of a system characterized by multiscale dynamics, a variety of positive and negative feedbacks, inhomogeneous physical and mathematical properties, and featuring intense fluxes across its subdomains. These features potentially result in a timedependent edge state, possibly characterized by nontrivial dynamics and nontrivial geometrical properties, and that cannot be found as a solution to a boundary value problem. This suggests an increased potential and scope for the methodology originally proposed by Eckhardt and collaborators; • We want to specifically use the edge tracking algorithm for identifying the edge states of the climate system separating the warm climate states from the snowball states. We want to study what lies in-between the two curves corresponding to the stable climates, portrayed in, e.g., Fig. 3a, and derive the long term statistics of the climatic edge states, mirroring what has been shown in the case of edge states in turbulent flows in, e.g., Fig. 4.
• As we observe that edge states are relative (in the reduced space of the basin boundary) attractors, we want to complement the characterization of their climatology with the investigation of their dynamical properties, and in particular understand whether they feature trivial, quasi-periodic, or chaotic dynamics, in order to classify the properties of the weather of such special states. Note that in this way we separate the slow global climatic instability (driven by the ice-albedo feedback) leading to the presence of multiple steady states, from the fast baroclinic instability [90] responsible for the variability of the weather and acting as negative feedback.

Model Formulation
The atmospheric component of the PUMA-GS model is provided by PUMA [39], which consists of a dynamical core: the dry hydrostatic primitive equations on the sphere (mapped laterally by the latitude φ and longitude λ), 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 ζ = ξ + 2νΩ E (where ξ is the vertical component of the relative vorticity, ν = sin φ, and Ω E = 2π/day is the angular frequency of the Earth rotation) the (horizontal) divergence D, the (atmospheric) temperature T a =T a + T a (separated into a time-independent arbitrary reference partT a and anomalies T a ) and the logarithmic pressure (normalized by the surface pressure p s ) σ = ln p/p s , read as follows: where , v being respectively the horizontal and vertical wind velocity components, and Φ is the geopotential height. Equations 8,9, and 11 express the conservation of momentum, Eq. 10 expresses the conservation of energy, and Eq. 12 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∇ 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 τ 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 τ 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 [39]: where (T R ) tp and z tp are the temperature and height of the tropopause, respectively, L (L) is the (average) lapse rate, [T s ] is the globally averaged surface temperature, and z(σ) is determined by an iterative procedure [39]. The above expressions indicate that the restoration temperature profile is anchored to the surface temperature T s . However, as Eq. 14 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 [T s ]. The surface temperature is taken to be governed by the 2D version of the GS EBM [25], which is in the form of Eq. 2, while the precise expression can be found in [38]. Since land masses are absent in this configuration, we are in fact dealing with an aquaplanet climate model. Two changes are introduced to the GS EBM evolution equations with respect to our previous study [38]: • we introduce a longitudinal component in the model (even if diffusion takes places only along the meridional direction); adding a second dimension is important because we introduce the following atmosphere-to-surface coupling term on the right-hand-side of the equation for the field of the surface temperature: k 3 (T a (σ = 1) − T s ), where T a (σ = 1) depends on latitude, longitude, and time; • we change the value of the original meridional diffusion coefficient, in order to make sure it represents the transport of the ocean only: in the original model the diffusion mimics the effects of both atmospheric and oceanic heat transports, but in PUMA-GS we have now a dynamical model for the atmosphere. We choose to reduce the value of the diffusion coefficient by a factor of 4 in order to make sure that in conditions resembling the present-day climate the ocean heat transport is few times weaker than the atmospheric one, in agreeement with observations and models outputs [8,114,3].
We note that T a (σ = 1) is obtained by linear extrapolation, according to T a (σ = 1) ≈ T a (σ = 0.9) + η(T s − T R (σ = 0.9)), 0 < η < 1. With η = 1 the coupling term is k 3 (T a (σ = 1) − T s ) ≈ k 3 (T a (σ = 0.9) − T R (σ = 0.9)). Generally T a (σ = 1) − T s = 0 (laterally inhomogeneous heating), but [T a (σ = 1)] = [T s ], where the overbar denotes averaging with respect to time. For our setup we chose: K −1 = 0.25 days, τ c = 30 days, τ f = 1 day.L = 0.0065 K/m, z tp = 12000 m. 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; and a number N lev = 10 of levels are taken. No orography is defined, i.e., we have a zonally-symmetric configuration, empirical functions of the GS EBM, like C(φ), depend on the latitude only. The equations are integrated numerically using a ∆t = 1 [hour] time step size. With such a setup we find bistability in a range of the relative solar strength µ = S * /S * 0 : the warm-to-cold (cold-to-warm) tipping point is found at µ W →SB = S * W →SB /S 0 0 ≈ 0.97 (µ SB→W = S * SB→W /S 0 0 ≈ 1.055). We take 18 equally spaced (∆µ = 0.005) sample values in order to study how the properties of the system change when µ is varied.

Algorithm for Tracking the Melancholia States
We adapt the edge tracking algorithm so that in the initializing phase of the procedure the basin boundary between the warm and snowball states is closely bracketed by two initial conditions. See Fig. 7 for a cartoon representing the procedure. The initial (j = 0) bracketing can be obtained as follows. One starts with two arbitrary points in phase space, x W,i,j and x SB,i,j , i = 0, j = 0 (small black disks in the schematics), belonging to the two different basins of attraction. These points need not to belong to the attractors themselves, which are marked by W and SB in Fig. 7 (red and blue objects). By iterative bisecting the basin boundary (thick green line -a simple geometrical object for convenience of visualization) is bracketed by points along a straight line in phase space a distance ε 1 apart. After each bisection, a control simulation (thick black line) reveals on which side of the boundary the midpoint (indicated by the horizontal arrows) time Warm climate Snowball Figure 8: Scheme of the tracking method used to construct the edge state in the climate model. We let the two initial conditions defined in Fig. 7 evolve until the two globally averaged surface temperatures differ by ε 2 . We then repeat the bisection as in Fig. 7 until a steady state is obtained. See details in the text.
is situated. This is realized by checking the vicinity of which attractor (dashed black circles) the trajectory enters starting from the midpoint, possibly after a long waiting time. It is easier to consider a suitably defined scalar indicator quantity, and measure ε 1 in this single dimension, and find scalar threshold values (indicated by horizontal dashed lines) that unmistakably indicate the outcome.
The best indicator to be used for all phases of the edge tracking algorithm is the globally averaged surface temperature [T s ], similarly to what was discussed in [38]. The basic reason for such a choice is simply the fact that at 0 th order the dynamics of the system can be approximated, as discussed above, by Eq. 1, where the energy landscape is shown in Fig. 2. In the case of such a simple model, for initial conditions near the saddle, the ice-albedo feedback pushes the system towards the W (SB) attractor through a monotonic increase (decrease) of temperature. In case of the more complex model, monotonicity does not strictly hold, but its violation is restricted to small temperature scales: compare Fig. 9 a) to Fig. 9 b), and Fig. 9 c) to Fig. 9 d), respectively.
After the initialization, we run the actual edge tracking algorithm; a cartoon is presented in Fig. 8. We launch the two simulations with initial conditions portrayed in Fig. 7 and stop them when the value of their globally averaged surface temperature is apart more than a given value ε 2 . We then repeat the bisection procedure outlined in Fig. 7 and define two nearby initial conditions whose globally averaged surface temperature differs by less than ε 1 . The two orbits originating from those two initial conditions end up in two different climates. But, since they are both initialized very close to the basin boundary, they remain near it for quite some time, quantified by the Lyapunov exponent corresponding to the dominating unstable dimension leading to the bistability. In the limit of ε 1 , ε 2 → 0, we construct the dashed line representing a trajectory on the basin boundary. We repeat the procedure until a statistically steady state is realized, corresponding to having reached the edge state indicated with E in Fig. 7. The procedure is then continued in order to be able to construct a trajectory long enough for reconstructing accurately enough the statistical properties of the edge state, and in particular for computing the maximum Lyapunov exponent (MLE).

Practical Implementation of the Edge Tracking Algorithm
Let's delve a bit into the details of how edge tracking works in two selected cases. Results are presented in Fig. 9, where we plot the time series of the globally averaged surface temperature [T s ] for the bracketing and control trajectories. The reason why atmospheric variables are not included in the choice of the observable has to do with the fact that they have much larger fluctuations than the surface temperature, because of the smaller heat capacity of the atmosphere with respect to the surface. As discussed above, in order to construct an efficient algorithm along the lines of what depicted in Fig. 8, we want the short-term natural fluctuations of the observable to be much smaller than the difference of their climatological values in the two attractors.
After some testing, we have chosen ε 1 = 0.1 K and, additionally, we have selected ε 2 = 2ε 1 , which implies that a single bisection is taken per edge tracking j−cycle, as done in [38]. Such a bracket size is much larger than the fast variability of [T s ], so that we have little risk of being influenced by spurious trends. The control trajectories are continued up to a threshold value chosen to be close to the mean value for the stable warm or cold climate.
In Fig. 9a) we show how the procedure works for µ = 1.03, with 9b) providing a zoom: we are able to conclude that the procedure converges already in the j = 1 cycle almost completely (bracketing trajectories shown in red), an edge state having [T s ] ≈ 241.8 K is identified, and we observe a (closely) monotonic (in fact, almost exponential) change of [T s ] as trajectories diverge from the edge state, pulled apart by the ice-albedo feedback.
Somewhat more interesting results are shown in Fig. 9c) and its zoomed-in version 9d), where the edge tracking procedure is applied for µ = 0.98. Figure 9d) clarifies that the time series of [T s ] of trajectories going to different attractors can, in fact, cross, so that in general no trivial monotonic increase/decrease of the chosen indicator should be expected when a complex model is considered. In other terms, the cartoon provided in Fig. 8 (and, a fortiori, ultra-simplied points of view as expressed by Eq. 1 and Fig. 2) should be taken with a grain of salt.
In the supplementary material, we include a movie showing the evolution of three orbits initialized near the Melancholia state realized for µ = 1.0 towards the W state, the SB state, plus one trajectory held on the edge state through the edge tracking algorithm described here. See the movie three trajs atm temp.mp4 -URL: https://youtu.be/mLYZiyzO8c4 (three trajs surf temp.mp4 -URL: https://youtu.be/OOYqUuG VUE) for the evolution of the atmospheric (surface) temperature field.

A Plethora of States beyond the Stable Climates
In what follows we will show how much additional information one can draw from the climate model studied here beyond the usual bifurcation diagrams of the form shown in Fig. 3. We will be able to construct the Melancholia states, characterize their properties in terms of symmetry, variability, and degree of chaoticity, and prove that, near each tipping point, their properties converge to those of the stable climate that is in the process of losing stability. Figure 10 provides a first comprehensive summary of our results and should be compared with Fig. 5 discussed in [38] and with Fig. 4b) published in [36]. We  Fig. 5 for the simple model studied in [38]. Figure 11 shows that the state corresponding to the warm attractor -#1 -features a positive MLE, i.e., it has a limited horizon of predictability, while the cold state -#2 -is a fixed point, so that its MLE vanishes. What we obtain in terms of qualitative properties of the snowball state is partially unsatisfactory. Previous studies performed using comprehensive atmospheric models feature much weaker (yet present) atmospheric variability in the snowball state than in the warm state; see e.g. [18]. In our case the atmosphere is too quiet, as a result of the very weak atmospheric meridional temperature gradient which is below (even if close to) the critical value needed for having sufficient baroclinicity supporting unstable motions. The main culprits for this behaviour we can indicate are the fact that a) the PUMA-GS has no continents and no orography, which  . Wheb the MLE is larger than zero we have a limited horizon of predictability due to chaotic dynamics. Note that while the MLE can be rigorously defined for the states #1, #3, and #5, a word of caution must be used in the case of states #4. Same number coding as in Fig. 10. Details in the text.
makes it impossible to have geographically limited regions where higher baroclinicity can be present, thus leading to generation of cyclones (as preferentially happens in our planet in the westward portions of the Atlantic and Pacific storm tracks [90]), and b) the simple ocean model we use is slightly too effective in transporting heat from the equator to the polar regions, thus reducing too much the meridional temperature gradient. Such -in our opinion, minor -limitation of our investigation could be overcome by repeating the analysis using the full PlaSim model with a setup as in [18]. Most interestingly, we are able to identify chaoticity in the Melancholia state (state #3): we are here describing the properties of the dynamics restricted to the basin boundary between the W and SB attractors, so that the global instability due to the ice-albedo feedback is filtered out, and weather-related processes are instead retained. In other terms, we have that the Melancholia state features a dynamics of weather qualitatively analogous to what is found in the W state, including full life-cycles for mid-latitude disturbances and a non-trivial active Lorenz energy cycle. We present additional evidences of this in the movie included in the supplementary material. See the evolution of the atmospheric temperature field in the edge state belonging to µ=1.0 in the movie three trajs atm temp.mp4 -URL: https://youtu.be/mLYZiyzO8c4.
One finds -compare with Fig. 12 -that, broadly speaking (one cannot expect a one-to-one correspondence), the value of the MLE of states #1, #2, and #3 is controlled by the intensity of the large scale meridional temperature gradient at surface, as suggested by the basics of baroclinic instability theory [90]. We note that regions of strong meridional temperature gradient for the surface temperature correspond by and large to where the variability of the atmospheric temperature near surface is more pronounced, as a result of baroclinic disturbances (not shown).
Let's go back to Fig. 10. As µ is increased between 1.010 and 1.015, a rather interesting phenomenon appears. The Melancholia state shown with marker #3 goes through a symmetry breaking according to the following pattern. After an extremely long transient of the order of 100 years or more (marker #4), a large longitudinal modulation of the atmospheric and surface fields suddenly appears (see in the supplementary material the movies edge symm break atm temp.mp4 -URL: https://youtu.be/u4--tRnBBS8 and edge symm break surf temp.mp4 -URL: https://youtu.be/Q4YAbU9O15U corresponding to µ = 1.025), so that about half of the planet warms up (marker #7, showing the meridional average at the warmest longitude, following the revolution of the field) and is in conditions relatively similar to those of the co-existing W state, while the rest of the planet cools down and reaches conditions relatively similar to those of the co-existing SB state (marker #9, showing the meridional average at the coldest longitude). The longitudinal modulation impacts all latitudes and rotates with a constant and extremely low angular velocity, and one can define the average properties of the the Melancholia state (marker #5) without the need for considering ultralong time scales (order of 10 3 years). In order to define the properties of the system in the states #7 and #9 one needs to compute averages in a reference frame co-rotating with the temperature wave. Note that even the direction of rotation is not the same for all values of µ.
In the supplementary material (see the movie edge symm break atm temp.mp4 -URL: https://youtu.be/u4--tRnBBS8) one can also see that the warm portion of the domain features a large temperature gradient and supports the growth of baroclinic disturances, which tend to decay as they reach the cold region, where the low temperature gradient does not provide suffiicient baroclinic forcing. In a loose sense, the Melancholia states realized as symmetry breaking of the transient state #4 resemble some sort of chimera states [115], yet unstable ones.
An additional aspect of our results should be highlighted. For µ = 1.045 we have a fundamentally different phase portrait for our model, which seems to be specific to a small neighborhood of this value of µ. We find that the system features three stable climates, which include, apart from the W and SB climate described elsewhere, the result of the bifurcation of the Melancholia state described by the markers #5, #7, and #9 into a stable state (see markers #6, #8, and #10, respectively), which displays similar features in terms of the presence of a semi-stationary pattern of strong longitudinal gradients of temperature. As far as we know, we are unaware of any other study reporting the existence of such an exotic climate state.
One is naturally bound to expect the presence of additional Melancholia states sitting in-between the three stable climates found for this value of µ = 1.045. We do not pursue this investigation, because it would lead us to what we see at this stage as unnecessary complication. It is important to note that it would have been extremely unlikely to find such an additional attracting climatic steady state, had we not used the edge tracking algorithm, because the width of parametric window where it exists is small, and rather limited is also the size of its basin of attraction.
The presence of a large time scale separation between the rotation and all the other dynamical processes allows for constructing the statistics embodied in the markers #7 and #9, which are computed following the slow rotating wave, in order to guarantee homogeneity in the results.  : Bifurcation diagram for a large scale meridional surface temperature gradient defined as the difference between low-latitude and high latitude 'boxes', as done in [38]. Same number coding as in Fig. 10. Details in the text.
Additionally, for µ ≈ µ W →SB , the properties of the Melancholia state get very close to those of the W state, in agreement with the fact that at the tipping point the edge state and the W attractor intersect, leading to the loss of stability. For µ ≈ µ SB→W , we have that the properties of the Melancholia state (marker #5), transient Melancholia state (marker #4), and SB state tend all to converge, as the tipping point is reached.
By considering the argument of time scale separation, it is also possible to construct a (pseudo-)MLE for state #4 (we cannot obviously take an infinite time horizon to compute it, because of the symmetry breaking process in action), whose value joins on with continuity with what is found for the Melancholia state denoted by marker #3, as shown in Fig. 11. It is also possible to compute the MLE for state #5. One finds that the values of the (pseudo-)MLE obtained for state #4 follow closely those of the actual Melancholia state #3 (and obeys the previously discussed relationship between value of the MLE and the strength of the meridional temperature gradient), while the MLE values obtained for state #5 are much larger (even larger than those realized in state #1), despite the presence of a weaker average meridional temperature gradient; see Fig. 12. This can be explained by considering that in state #5, in addition to the presence of a meridional temperature gradient, the system experiences a large longitudinal temperature gradient too, which also supports the growth of instabilities, separating the warm and cold portions of the planet (states #7 and #9, respectively). The strongly nonlinear dynamics of cyclones when transiting from the relatively cold to the relatively warm regions (see in the supplementary material the movie edge symm break atm temp.mp4 -URL: https://youtu.be/u4--tRnBBS8) is a manifestation of this effect.
A further characterization of the dynamics of the model can be obtained by looking at the power spectrum of an atmospheric variable, keeping in mind that a broadband spectrum is a signature of chaotic dynamics. The power spectra for the T a (σ = 0.9) in a grid point situated at 30 • N are presented in Fig. 13. We have that, in agreement with what is reported in Fig. 11, the states corresponding to the W attractor feature a broadband spectrum for all values of µ, while no variability is present for the SB states (not shown).
As for states #3 and #4, we have that for increasing µ the total intensity S and the breadth of the spectrum decreases as well as featuring an overall shift to lower frequencies, as a result of the decrease of the meridional temperature gradient (see Fig.  12), and, by thermal wind relation [90], of the intensity of the jet, which defines the time scales of the mid-latitude. For very large values of µ, the system reaches a quasi-periodic behaviour, where one or few frequencies are present.
When constructing the power spectra of states #7-#10 (note that no information on these states could be obtained by looking at the MLE), we make sure to discard the effect of the longitudinally slowly moving temperature wave. The power spectra of the cold regions of the longitudinally-modulated climate, corresponding to states #9 and #10, are extremely weak, corresponding to the fact that the variability is only inherited from the active dynamics taking place in the warm regions corresponding to the states #7 and #8. The variability is intense with broad spectrum for states #7 when µ < 1.04, in agreement with the fact that a large horizontal temperature gradient is found across the mid-latitudes.

The Geometry of the Basin Boundary
A very nontrivial aspect of the problem we are studying is the characterization of the geometrical properties of the boundary between the basins of attraction of the stable climates. The cartoons shown in Figs. 7-8 suggest that the two basins of attraction are separated by a simple manifold of co-dimension one. Clearly, the existence of two separate true attractors (as opposed to the case of transient turbulence studied by Eckhardt and collaborators) implies that the co-dimension of the boundary can be at most one, because otherwise one would automatically have holes allowing for transitions between the two basins of attraction, in contrast with their property of being invariant. Nonetheless, nothing prevents in principle the boundary from being a more complex geometrical object of co-dimension smaller than one, such that the intersection between such a boundary and a line connecting the two competing attractors has a non-trivial Cantor-like structure. This issue has been studied in simple mathematical models in, e.g., [ The difference between the globally averaged surface temperature of the 1 st and of the 1024 th is 0.5 K. Outcome 1 indicates that the trajectory ends in the warm state, and outcome 0 indicates that the trajectory ends in the snowball state. b) Same as a), but for a different set of initial conditions. Details in the text. 43,37], and has practical relevance as the presence of a fractal boundary makes the prediction of the final state given the knowledge of the initial conditions of the system (with finite precision) a very non-trivial task. In other terms, it controls what Lorenz called the predictability of the second kind [7,8].
As discussed in, e.g., [116], the presence of a chaotic edge state is a necessary (but not sufficient) condition for the existence of a geometrically complex basin boundary. In fact, the geometrical complexity of the boundary requires that that the MLE characterizing the separation of trajectories inside the basin boundary is larger than the Lyapunov exponent describing the separation of trajectories in the unstable direction (and leading to the orbits ending up in either attractor) [40,42,43,37].
Following a suggestion by A. Pikovsky during a public presentation of some preliminary results then included in the present paper, we have decided to investigate the geometrical properties of the basin boundary in which the Melancholia states are embedded. We have considered the case µ = 0.98, for which we have already found a chaotic Melancholia state, and have considered two nearby initial conditions belonging, respectively, to the basin of attraction of the warm and of the snowball state. These two initial conditions are represented by the first bracket seen in Fig. 9d). The two initial conditions are very similar in terms of dynamical and thermodynamical properties and differ by only 0.5 K in terms of globally averaged surface temperature. We have then considered additional 1022 initial conditions obtained by equispaced convex linear interpolations of the two original initial conditions, which results in the fact that their corresponding globally averaged surface temperatures are ordered and equispaced, so that the globally averaged surface temperature of two neighbouring initial conditions differs by about 5 × 10 −4 K.
We have then integrated forward each of these initial conditions and have labelled with 1 the orbits ending in the warm state and 0 the orbits ending in the snowball state. The results are shown in Fig. 14a). Against intuition, one finds an extremely complex geometric structure across the boundary. Very often, nearby initial conditions belong to different basins of attraction. Similar results are obtained by repeating the investigation in another region of the phase space: see Fig. 14b), where we follow the procedure detailed above starting from the fifth bracket in Fig. 9d). Clearly, at finite precision one cannot in principle exclude the possibility that the basin boundary is a regular -yet very tightly folded -manifold, rather than being a strange geometrical object. Further analysis is shown in Fig. 15, where we have computed the information dimension of the intersection set between the segment including the 1024 initial conditions above and the basin boundary using the data shown in Fig.14a). Entirely compatible results are obtained using, instead, the data in Fig.14b). We find an information dimension of about 0.98, which indicates that -at least within a certain region -the boundary has almost full dimension. As a result, our ability to predict the final state of the system (warm or snowball) given the initial conditions (which can be measured using the uncertainty exponent [79]) is extremely low. Going back to the mathematical conditions discussed in [40,42,43,37], we must have that the instability inside the basin boundary is (much) stronger than the instability in the unstable direction. And here physical intuition comes to our help: the instability occurring inside the basin boundary is associated to weather processes, and corresponds to the baroclinic processes (having characteristic time scales of tens of days) responsible for the generation of eddies in the Melancholia state, see Fig. 11. Instead, the instability in the transversal direction results from the ice-albedo feedback, and acts on time scales of years [25,38]. Concluding, the strangeness of the basin boundary comes from the large time scale separation between weather and climatic unstable processes.
A natural question one may ask is why, despite such a geometrical complexity, the edge tracking algorithm still works well. The reason is that the bisection algorithm removes the motions in the unstable direction associated to the ice-albedo feedback. Therefore, the constructed pseudo-trajectory along the basin boundary is constrained to explore the saddle set with no regard to its geometry and to be eventually attracted to the edge state.

Summary and Conclusions
In this paper we have attempted to provide a review of some ideas concerning the study of critical transitions in the context of climate science. We have first briefly recapitulated under which conditions one can expect a regular response of the climate system to forcings, and have clarified that the response theory introduced by Ruelle gives the mathematical framework for addressing comprehensively the problem of predicting how the statistical properties of the climate system change as a result of (possibly time-dependent) perturbations to its dynamics. We have then clarified that the lack of regularity in the response of the system to perturbations, associated to a slowing-down in the rate at which the system decorrelates, flags being near to conditions where critical transitions take place.
In order to have a more complete picture of climate dynamics, one needs to take a global point of view in the study of the phase space of the climate system and consider the possibility of the existence of multiple steady states. We have then briefly discussed the literature concerning the existence of multistability in the Earth system, with a special focus on the key scientific issue regarding the co-existence of snowball and warm climate conditions given the current astrophysical and astronomical conditions. In this regard, we have presented a hierarchy of conceptual to more realistic mathematical models commonly used for studying multistability in a geophysical (and not only) context and have underlined virtues and limitations of using a stochastic dynamical point of view.
While typically one wants to study the properties of the co-existing climate states and of, e.g., the noise-induced probability of jumping from one to the other, the question we find most interesting and indeed novel regarding the specific climatic problem we focus on is to define what -dynamically -lies in-between the multiple (say, two) coexisting steady states of the climate system. In other terms: is there a special dynamical configuration that generalizes the saddle point in an energy landscape and acts as gate facilitating the transitions between the two attractors?
The edge state introduced by Eckhardt and collaborators for studying turbulent fluid systems has proven rather useful for our study. The edge state is the attracting set embedded in and relative to the boundary between the two basins of attraction of the co-existing attractors. The edge state attracts, in fact, almost all the initial conditions belonging to the boundary of the basins of attraction. The edge state cannot be found by direct numerical integration, but can be constructed through a control algorithmthe edge tracking algorithm -that allows for bracketing nearby trajectories evolving on either side of the basins boundary.
In recent years we have followed a scientific programme aiming at a comprehensive view on climate response that has led us in the last years to study -using the same climate model -a) the macroscopic thermodynamical properties of the climate system associated to the critical transitions responsible for the snowball/warm bistability; b) the smooth climate response to perturbations using the Ruelle response theory; and c) the breakdown of the smooth response near the critical transitions due to slow decay of correlations using a the transfer operator method in a reduced space.
In this paper we have aimed at completing the puzzle by taking a global point of view with the analysis of dynamical regimes that do not correspond to attracting steady states of the climate, but rather constitute the edge states sitting in-between co-existing climate states. We refer to these states as Melancholia states. We have decided to stick to the point of view of statistical mechanics that proposes to investigate complex systems using a purely deterministic dynamics, with the hope of being able to understand more deeply the effect of adding noise.
By adapting the edge tracking algorithm to a climate model that is somewhat simpler than what was used in previous investigations (yet Earth-like), we have been able to identify a plethora of climate states and study the properties of the Melancholia states existing in the region of bistability where snowball and warm climate attractors coexist.
Our control parameter is µ, the ratio between the considered and the present-day solar constant. The tipping points are in fact reached when one of the climate attractors collides with the Melancholia state. We provide ample diagnostics for characterizing such states, including some movies that are available as supplementary material. The Melancholia states are obviously unstable with respect to the ice-albedo feedback, and, once such an instability is kept under control, they exhibit in some parametric range the kind of variability usually associated to regular climate conditions. The Melancholia states feature, e.g. mid-latitude cyclones growing as a result of baroclinic instability fuelled by the meridional temperature gradient, and have a chaotic behaviour associated to a limited horizon of predictability.
As µ is increased from the value corresponding to the tipping point from the warm to the snowball state, one finds that the Melancholia states undergo a symmetry breaking bifurcation, after which they are characterized by strongly zonally non-symmetric conditions, in the form of a slowly rotating temperature anomaly wave. This regime is a sort of mix of the co-existing warm and snowball attractors, and feature a complex dynamics at the interface between the cold and warm regions. In a small window of values of µ, one finds three stable climate attractors, but this regime is not further explored in this paper. After such a window, the asymmetric Melancholia state is recovered and converges to the cold attractor at the tipping point associated to the transition from the snowball to the warm state.
The Melancholia states discussed in this paper correspond to climatic regimes that cannot be realized in nature nor can be obtained by a forward integration of a numerical model. We would like to suggest two arguments supporting the idea that the Melancholia states have, nonetheless, great physical relevance.
First, they can be used to flag potentially dangerous forcings acting on systems originally living in one of the coexisting attractors. Obviously, forced transitions between the two basins of attraction can take place also bypassing the Melancholia state, but, surely, if the forced system resembles more and more a Melancholia state, the risk of catastrophic changes taking place if the forcing is not stopped is very real. Clearly, the risk is much higher when the unperturbed system lives in an attractor close to the Melancholia state -so in a parametric range near the critical transition.
Second, as we have discussed in the latter part of our paper, Melancholia states might morph into actual stable yet very nontrivial stable climate states for small changes in the value of the parameters. These exotic climate states might be extremely hard to discover through direct numerical integration, given the limited size of their basin of attraction and the small parametric window where they actually exist.
The results presented provide a first example of a successful reconstruction of the global properties of a climate model possessing multiple steady states and complex regimes of motion, and hope to provide a positive stimulation in the direction of having a more thorough understanding of the properties of the climate system. It seems then relevant to extend the study of Melancholia state for the Earth and for other planets, as they may be key to understanding the evolution and the statistics of observed planetary atmospheres.
Additionally, we have shown that the basin boundary between the basins of attraction has a complex geometrical structure, as a result of the very different time scales associated to weather-like (baroclinic processes) and climatic instabilities (icealbedo feedback). Obtaining a rather detailed geometrical description in such a highdimensional system as the one analysed here seems quite an achievement. The presence of a thick basin boundary with a complex geometry makes the prediction of the final state given a finite precision knowledge of the initial conditions extremely challenging, thus indicating a reduced predictability of the second kind.
As a future investigation, we aim at relating the dynamical properties of the edge state to the statistics of noise-induced transitions between competing attractors, starting from simple up to more complex multistable models. As the edge state is the gate controlling the noise-induced transitions between the two co-existing attractors, its dynamical properties are crucial in assessing how likely a random forcing might lead to a catastrophic transition to another basin of attraction. Therefore, one could hope to derive estimates or semi-quantitative statements relating the properties of the edge state to the escape rates one can study via the Freidlin-Wentzell theory. Additionally, we wish to analyse in detail the geometry of the basin boundary and investigate the relationship between the escape rate and Lyapunov exponent along the unstable direction, along the lines indicated in [116].
Additionally, we remark that the existence of more than two stable states leads to the existence of a potentially more complex partition of the phase space in different basins of attraction, and of an ensuing more complex population of edge states sitting in-between different pairs of attractors.
In terms of specific geophysical relevance, it seems relevant to go deeper in the properties of the Melancholia states found in this study and analyse in greater detail a) more specific physical reasons behind the presence of such a complex basin boundary (e.g. the role of localised instabilities in the physical space) ; b) the physical mechanisms behind the symmetry break found in the upper range of bistability, and, c) indeed, the apparent exchange of stability that morphs a Melancholia state into a new stable climate state and viceversa for small modulations of the solar constant.
The properties and relevance of the newly found stable climate state characterized by the break-up of the zonal symmetry and by intense horizonal temperature gradients seem indeed worth investigating. As far as we are aware, no previous studies have found any qualitatively analogous stable state.
A minor pitfall of the model used in this study is that it features a stable atmospheric circulation with no fluctuations in the snowball states. Instead, the snowball states obtained using more comprehensive models have weak yet non-vanishing atmospheric variability, as baroclinic instability is active in the atmosphere, despite the low meridional temperature gradient. Hence, repeating our analysis with a climate model like the one used in [18] seems a useful future exercise.
As discussed earlier in the paper, a more sophisticated analysis of the mechanisms behind the multistability of the Earth climate using models able to resolve accurately the oceanic time scales suggests the possibility of an additional climate state (the slushball state) co-existing with the warm and the snowball state. A complete analysis of the Melancholia states situated in-between the possible climate states would definitely provide a more complete picture of the global dynamical properties of the Earth system and possibly lead to identifying additional, nontrivial stable climates of relevance for paleoclimate and for the study of exoplanetary systems.
Kirk for the great support on many computational aspects of this paper. The authors acknowledge the support of the FP7-ERC StG Grant NAMASTE, of the Horizon2020 Project CRESCENDO, and of the DFG Cluster of Excellence CliSAP. VL acknowledges the support of the DFG project MERCI. VL wishes to thank GGG for many things, which include forcing him to watch Melancholia.