Multistability and rare spontaneous transitions in barotropic β -plane turbulence
aa r X i v : . [ phy s i c s . a o - ph ] S e p Multistability and rare spontaneous transitionsbetween climate and jet configurations in abarotropic model of the Jovian mid-latitudetroposphere
Eric Simonnet ∗† Joran Rolland & Freddy Bouchet, ‡ September 22, 2020
Abstract
We demonstrate that turbulent zonal jets, analogous to Jupiter ones,which are quasi-stationary are actually metastable. After extremely longtimes they randomly switch to new configurations with a different num-ber of jets. The genericity of this phenomena suggests that most quasi-stationary turbulent planetary atmospheres might have many climatesand attractors for fixed values of the external forcing parameters. A keymessage is that this situation will usually not be detected by simply run-ning the numerical models, because of the extremely long mean transitiontime to change from one climate to another. In order to study such phe-nomena, we need to use specific tools: rare event algorithms and largedeviation theory. With these tools, we make a full statistical mechanicsstudy of a classical barotropic beta-plane quasigeostrophic model. It ex-hibits robust bimodality with abrupt transitions. We show that new jetsspontaneously nucleate from westward jets, with an exponentially smallprobability as the Ekman dissipation decreases (Arrhenius law). This phe-nomenology is controlled by two different types of instantons in the limit ofvanishing noise. Moreover, we are able to compute the saddles of the cor-responding effective dynamics. We also investigate the detailed dynamicsof solutions having three alternating jets. We uncover an unexpectedlyrich dynamics governed by the symmetric group S of permutations, withtwo distinct families of instantons, which is a surprise for a system whereeverything seemed stationary in the hundreds of simulations of this modelperformed in the past. ∗ Eric Simonnet, [email protected] † INPHYNI, UMR7010 CNRS-UNS, 1361, route des Lucioles, 06560 Valbonne ‡ Univ Lyon, Ens de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, Lyon,France Introduction
Many complex physical and natural systems have more than one local attrac-tor. The dynamics can then settle close to one of these local attractors fora very long time, and randomly switch to another one. Well-known exam-ples are phase transitions in condensed matter or conformational change ofmolecules in biochemistry (Maragliano et al., 2010). Such situations of multi-stability are ubiquitous also in geophysics and for turbulent flows (Ravelet et al.,2004; Bouchet and Simonnet, 2009). The reversal of the Earth magnetic fieldon geological timescales, which is due to turbulent motion of the Earth’s metalcore (Berhanu et al., 2007), is a paradigmatic example.The climate system is no exception and global multistability and abrupttransitions existed from the global Neoproterozoic glaciations (snowball Earthevents) (Pierrehumbert et al., 2011) to glacial-interglacial cycles of the Pleis-tocene (Paillard, 1998). Very often, those transitions are due to the internal dy-namics and are not caused by changes of external parameters. This is the casefor instance for the fast Dansgaard–Oeschger events (Dansgaard et al., 1993;Ditlevsen et al., 2007; Rahmstorf, 2002). More local multistability also exist inpart of the climate system, for instance the Kuroshio bimodality (Schmeits and Dijkstra,2001; Qiu and Miao, 2000) in the North Pacific.If multistability and abrupt transitions between distinguishable states aregeneric features of complex dynamical systems, why should planetary atmo-sphere dynamics be an exception? Actually Lorenz was already rising thisquestion in 1967 (Lorenz, 1967). The hypothesis of a transition to a superrotat-ing Atmosphere with westward equatorial jets has been addressed many times(Held, 1999) and observed in many numerical experiments (Arnold et al., 2012).Recently the more specific question of an abrupt transition to superrotation hasbeen studied carefully (Herbert et al., 2020). On Earth, superrotation may haveplayed a role in the climate of the past: it was observed in numerical simula-tions of warm climates such as the Eocene (Caballero and Huber, 2010), andit has been suggested that it could explain the permanent El Ni˜n o conditionsindicated by paleoclimatic proxies during the Pliocene (Tziperman and Farrell,2009).The aim of this paper is to address the hypothesis of planetary atmospheremultistability and abrupt transitions not for equatorial jets but for midlatitudeeddy-driven jets. This is a natural hypothesis as multiple midlatitude jets areobserved on many planets, and that even with Earth conditions slight changesof parameters lead to a different number of jets (Lee, 1997, 2005) and jets ofdifferent nature (Kaspi and Flierl, 2007). Because it is difficult to change froma state with n zonal jets, to a state with n + 1 zonal jets continuously, it isnatural to expect multistability, abrupt transitions, and spontaneous transitionsbetween those different states to be generic. We demonstrate this in this paper.We will focus on Jupiter like atmosphere with its alternating zonal bandswith colors correlated with the troposphere vorticity. They correspond to suc-cessions of westward and eastward velocity jets. There is a strong asymmetrybetween eastward and westward jets: eastward jets form cusps at their maximum2 hite ovals appear in 1939−1940 Figure 1: Upper left panel: Jupiter’s zonally averaged zonal velocity profileobserved by Voyager, Cassini (Porco et al., 2004) and Hubble (Simon et al.,2015). Upper right panel: Great Red Spot contraction after one century (left:from (Clerke, 1893), right: NASA). Lower panel: three white ovals appearedin 1939-1940 (Rogers, 1995) and later merged to from the oval BA, also called”Red Spot Junior” in 2006 when it turned red.velocity whereas westward jets have an almost parabolic profile (Porco et al.,2004; S´anchez-Lavega et al., 2008). What is striking is the stability of thesejets: despite the strong turbulent activity at small scales, one hardly notice anydifference between the two separate measurements by Voyager and Cassini 20years apart (see Fig.1). The Jovian large-scale atmosphere, including the GreatRed Spot (GRS), in fact appears to be stationary for decades. However, a fan-tastic event occurred in 1939-1940 when Jupiter lost one of its jets (Rogers,1995). Three white anticyclones were created which started to merge in thelate 90s into a large white anticyclone (see Fig.1). As such, Marcus (2004)called the 1939–1940 event a global climate change on Jupiter. The mechanismresponsible for the disappearance of one of the jet remains completely unknown.A key scientific aspect of bistability situations is that the time scale for ob-serving a spontaneous transition from one attractor to another can be extremelylong. For instance on Jupiter, no such events as the 1939-1940 jet disappearanceoccurred since then. The typical timescale for such transitions is thus longerthan 80 years, while energy balance for jet-mean flow interactions, Rossby wave3ynamics, and radiative effects have time scales much smaller ranging from oneJupiter day to several years. For this reason, we argue in this paper that mostof time, scientists running numerical simulations are not aware that their modelmay have more than one attractor, just because they never simulate it longenough. Actually to simulate the model long enough and wait is most of thetime impossible with currently available computational capabilities. This time-scale separation issue seems an impassable barrier. Breaking this constraint andstudying unobserved phenomenologies for such transitions are two fascinatingscientific aims. For instance in this paper, we will demonstrate that a barotropicmodel of Jupiter jets has a huge number of climate (attractors) for a fixed set ofboundary conditions and forcing parameters, with an amazingly rich long-timedynamical structure for the transitions from one attractor to another. This wascompletely unexpected and never observed, even if analogous models have beenrun by hundreds of groups before.In order to face this time scale separation issue, we need to build com-pletely new theoretical and numerical tools. For this it is scientifically soundto go back to the basis and to start with the simplest possible model that con-tains the phenomenology of Jupiter troposphere dynamics. We will considerforced–dissipative barotropic turbulence on a beta plane, in which small-scaleturbulence self-organizes into large-scale coherent structures, with upscale fluxesof energy. Although such a barotropic description of Jupiter troposphere haslong been recognized to be relevant (Galperin et al., 2014), it has clearly somelimitations as the energy transport through baroclinic instability is describedonly roughly by the stochastic forcing, and models with baroclinic instabilitymay lead to dynamics with different qualitative aspects (Jougla and Dritschel,2017). We come back to baroclinicity and multiple attractors in the conclusion.The phenomenology of stochastically forced barotropic turbulence has beendescribed for a very long time in the literature (Vallis and Maltrud, 1993; Dritschel and McIntyre,2008; Galperin et al., 2001; Bakas and Ioannou, 2014; Tobias and Marston, 2013).Starting from rest, as one forces randomly at scales smaller than the Rhinesscale, a 2D inverse energy cascade develops during a transient stage. Thiscascade is stopped at a scale of the order of the Rhines scale, where Rossbywaves start to play a fundamental role and energy transfer is inhibited. It iswell-known that this barrier is anisotropic and take the form of a dumbbellshape which favor zonal structure (Vallis and Maltrud, 1993). Zonal jets thenform, with a typical width of order of the Rhines scale, while jet waves oftencalled modified Rossby waves, or Zonons (Galperin et al., 2001), play an im-portant dynamical role. (Galperin and Read, 2019) is a survey of the latesttheoretical, numerical and experimental advancements atmosphere and oceanjets dynamics and their interactions with turbulence. As illustrated by manycontributions in (Galperin and Read, 2019) the beta-plane barotropic model hasbeen a test-bed for many phenomenological and theoretical studies of jet dy-namics, for instance second-order closures (also called S3T or CE2 expansions)(Farrell and Ioannou, 2003; Srinivasan and Young, 2011) or analogies with pat-tern formation (Parker and Krommes, 2013). Second order order closures, orthe related quasilinear approach, has been argued to be exact in the limit of4ime scale separation (Bouchet et al., 2013) between the zonal jet time scaleon one hand, and the eddy relaxation time on the other hand. This is a limitrelevant for Jupiter. Moreover (Woillez and Bouchet, 2019) has given a fairlycomplete analytical description of the zonal jet structure, including a dynamicalexplanation of the westward-eastward jet asymmetry, and a precise descriptionof the westward jet cups and their regularization.Those statistical approaches, for instance S3T, actually allowed to predictthe non-uniqueness of solutions (multistability) (see Farrel and Ioannou (2007),or Parker and Krommes (2013) for the transition regime to zonal jets). More-over, wonderful zonal jet experiments have demonstrated bistability in fairlybarotropic regimes (Lemasquerier et al., 2020).However, no spontaneous transitions between attractors were observed be-fore the work Bouchet et al. (2019). In this paper, with the beta-plane barotropicmodel with stochastic forcing we were able to observe multistability betweenstates with two eddy-driven zonal jets and states with three eddy-driven zonaljets, and to observe spontaneous transitions between those two states. Be-cause of the huge numerical cost, it would have been impossible to observeseveral of these very rare transitions, and to study their dynamics and prob-ability. In order to face this critical practical problem, we had to use a rareevent algorithm in order to concentrate the computational cost on transitionstrajectories from one attractor to another, rather than on extremely long pe-riod when nothing dynamically interesting occurs. Rare event algorithms aimat this goal (Kahn and Harris, 1951; P. Del Moral, 2004). In a dynamical con-text, they have been first applied for complex and bio-molecules, see for in-stance Metzner et al. (2009); Hartmann et al. (2013). More recently, in or-der to progressively go towards genuine geophysical applications, they havebeen applied to Lorenz models (Wouters and Bouchet, 2016), partial differen-tial equations (Rolland et al., 2016), turbulence problems (Grafke et al., 2015;Laurie and Bouchet, 2015; Ebener et al., 2019; Bouchet et al., 2019; Lestang et al.,2020), geophysical fluid dynamics (Bouchet et al., 2019), and climate appli-cations (Ragone et al., 2018; Webber et al., 2019; Ragone and Bouchet, 2019;Plotkin et al., 2019). Some approaches through minimum action methods, re-lated to large deviation theory, are reviewed in Grafke and Vanden-Eijnden(2019). In Bouchet et al. (2019), and in this paper, we rather use the Adap-tive Multilevel Splitting (AMS) algorithm, a very simple rare event algorithmwhich is well suited to study rare transitions (C´erou and Guyader, 2007) (seeC´erou et al. (2019) for a short review, and Rolland et al. (2016); Bouchet et al.(2019) for application for partial differential equations and geostrophic turbu-lence problems). In Bouchet et al. (2019), thanks to the AMS algorithm, wewere able to demonstrate that rare transition path concentrate close to instan-tons and that transition rates follow an Arrhenius law. The concentration oftransition paths close to predictable trajectories, called instantons, is a fasci-nating property shared by many rare events. This has been now been observedin several turbulent flow applications (Laurie and Bouchet, 2015; Grafke et al.,2015; Bouchet et al., 2014; Dematteis et al., 2019) and other geophysical appli-cations, for instance the dynamics of the solar system (Woillez and Bouchet,5020). The idea that Jupiter’s jet transitions should follow instantons and bedescribed by Arrhenius laws, like condensed matter phase transitions, is reallystriking.The aim of this paper is to develop and extend those results in several ways,pushing the power of rare event algorithm for studying unobserved phenomenaso far, for instance rare transitions for eddy-driven jets. We clearly demonstratethe generic nature of multistability by performing bistability experiments be-tween the attractors with two and three jets respectively. New results include acomplete description of the three jet structure. Amazingly there are actually sixdifferent types of three-jet states with different jet spacing. The dynamics canremain quasi-stationary close to the two-jet state for a very long time, switchto the three-jet states; change several times its type of three-jet states, beforecoming back to the two-jet states. We confirm the instanton phenomenologyfor transition between the two-jet states and the three jet-states, but we alsodemonstrate it for transitions between different three-jet states. This fascinat-ing complex phase space structure for the slow dynamics, occurring on verylong time scales, has profound implication on our understanding of Jovian likeplanets. It suggests that the observed state of Jovian troposphere is most-probably one among many different ones, for a fixed set of forcing parameters.This profound fact suggests a very simple explanation of the observed asym-metry between the northern and southern hemisphere structure of Jupiter. Inthe conclusion we explain why we believe that this observation for a barotropicdynamics is expected to be robust for more comprehensive models. This opena completely new set of scientific questions for planetary atmosphere studies.Another new result of this paper concerns the hypothesis of barotropic ad-justment of eddy driven jets. An eddy driven jet with a time-scale separationbetween eddy dynamics and zonal jet time scale, should be unstable in orderto transfer meridionally zonally-averaged momentum, while it should be at thesame time stable if quasi-stationary zonal jet states are observed for very longtimes. This apparent contradiction leads to the hypothesis of barotropic ad-justment: the state of the system should be marginally stable (or unstable) inorder to fulfill these two seemingly contradictory requirements. The relevanceof this adjustment hypothesis has been recently discussed in the context of Jo-vian planets and for a hierarchy of models (Read et al., 2020b,a). Through anempirical analysis of the Rayleigh-Kuo criteria, we demonstrate in this paperthat the zonal jets actual constantly remain in a state of marginal stability.The striking new results is that this remains true also during the transitionsbetween different attractors: The nucleation of a new jet, or the merging of twojets, both occur within barotropically adjusted states.The paper is organized as follow: we first show some bimodality results in-volving transitions between two-jet and three-jet solutions in Section 3.1. Westudy more in details these rare transitions using an advanced rare event algo-rithm (AMS) in Section 3. Section 4 focuses on the dynamics of the three-jetstates only and a summary of the full effective dynamics is given. Section 5addresses the question of the Arrhenius law when the Ekman dissipation α de-6reases. We then conclude in Section 6. We consider in the following the barotropic quasi-geostrophic equations, witha beta plane approximation for the variation of the Coriolis parameter. Theequations in a doubly periodic domain D = [0 , πLl x ) × [0 , πL ) read ∂ t ω + v · ∇ ω + β d v y = − λω − ν n,d ( − ∆) n ω + √ ση, (1)where v = e z × ∇ ψ is the non-divergent velocity, v y the meridional velocitycomponent, ω = ∆ ψ and ψ are the vorticity and the stream function, respec-tively. λ is a linear friction coefficient, ν n,d is a (hyper-)viscosity coefficient, and β d is the mean gradient of potential vorticity. η is a white in time Gaussianrandom noise, with spatial correlations E [ η ( r , t ) η ( r , t )] = C ( r − r ) δ ( t − t )that parametrizes the curl of the forces (physically due, for example, to the effectof baroclinic instabilities or convection). The correlation function C is assumedto be normalised such that σ represents the average energy injection rate, sothat the average energy injection rate per unit of area (or equivalently per unitof mass taking into account density and the layer thickness) is ǫ = σ/ π L l x .For atmospheric flows, viscosity is often negligible in the global energy bal-ance and this is the regime that we will study in the following. Then the mainenergy dissipation mechanism is linear friction. The evolution of the averageenergy (averaged over the noise realisations) E is given by dEdt = − λE + σ. In a stationary state we have E = E stat = σ/ λ , expressing the balance be-tween forces and dissipation. This relation gives the typical velocity associatedwith the coherent structure U ∼ √ E stat /L ∼ p ǫ/ λ . As will be clear in thefollowing, we expect the non-zonal velocity perturbation to follow an inviscidrelaxation, on a typical time scale related to the inverse of the shear rate. As-suming that a typical vorticity or shear is of order s = U/L corresponding toa time τ = L/U , it is then natural to define a non-dimensional parameter α asthe ratio of the shear time scale over the dissipative time scale 1 /λ , α = λτ = L r λ ǫ . When β is large enough, several zonal jets can develop in the domain. Animportant scale is the so-called Rhines scale R which gives the typical size of7he meridional jet width: L R = ( U/β d ) / = (cid:0) ǫ/β d λ (cid:1) / . We write the non-dimensional barotropic equation using the box size L as alength unit and the inverse of a typical shear τ = L/U as a time unit. We thusobtain (with a slight abuse of notation, due to the fact that we use the samesymbols for the non-dimensional fields): ∂ t ω + v · ∇ ω + βv y = − αω − ν n ( − ∆) n ω + √ αη, (2)where, in terms of the dimensional parameters, we have ν n = ν n,d τ /L n , β = β d Lτ . Observe that the above equation is defined on a domain D = [0 , πl x ) × [0 , π ) and the averaged stationary energy for ν n ≪ α is of order one. Pleaseobserve that the non dimensional number β is equal to the square of the ratioof the domain size divided by the Rhines scale. As a consequence, according tothe common belief, the number of jets should scale like β / when β is changed.The phenomenology of stochastically forced barotropic turbulence has beendescribed for a very long time in the literature (Vallis and Maltrud, 1993; Dritschel and McIntyre,2008; Galperin et al., 2001; Bakas and Ioannou, 2014; Tobias and Marston, 2013).Let us remind briefly the main aspects. Starting from rest, as one forces ran-domly at scales smaller than the Rhines scale, a 2D inverse energy cascade devel-ops during a transient stage. This cascade is stopped at a scale of the order of theRhines scale, where Rossby waves start to play a fundamental role and energytransfer is inhibited. It is well-known that this barrier is anisotropic and takethe form of a dumbbell shape which favor zonal structure (Vallis and Maltrud,1993). Zonal jets then form, with a typical width of order of the Rhines scale,while jet waves often called modified Rossby waves, or Zonons (Galperin et al.,2001), play an important dynamical role. For values of α of order one or larger,the system settles in a statistically stationary state with such a phenomenology.For smaller values of α , the jets become very strong. Once such strong jets areformed, they progressively expel the waves. For smaller values of α no waveexist anymore, and both the inverse cascade phenomenology and the modifiedRossby wave phenomenology become irrelevant. The system then settles in astatistically stationary state, where the dynamics is dominated by the strongzonal jets, the average width of which is still approximately given by the Rhinesscale, and with fluctuations of order √ α . Most of the energy is then directlytransferred from the forcing scale to the jet scale through direct interaction withthe average flow. This process, which is non local in Fourrier space, is no morea cascade one (Bakas and Ioannou, 2014).In the following we will be interested in this strong jet regime, which isrelevant for Jupiter. Except if otherwise stated, all the computations of thispaper will be performed with the parameters α = 1 .
20 10 − , ν = 1 . − , andusing a stochastic force with a uniform spectrum in the wave number band | k | ∈ [14 , α t
50 100 150 2000246 −1−0.500.51 y α t
50 100 150 2000246 −101 − . . − − . . β = 5 β = 10 Figure 2: Jets sampled in direct numerical simulations. Left panels show localtime averages of the zonally averaged vorticity (red curve) and velocity (greencurve) as a function of latitude y ( y is the vertical axis), for α = 1 . · − .Middle panels show the corresponding Hovm¨oller diagrams of the vorticity ω (the color represent the value of the zonally averaged vorticity as a function oftime and latitude). Right panels show one of the corresponding instantaneoussnapshots of the vorticity field. The first line shows these figures at β = 5, thesecond line shows these figures at β = 10.9urnover time or more. If α would be decreased, these jets would become evenmore stable. The averaged vorticity has a saw-tooth profile: it is composed ofhomogenized potential vorticity areas where the potential vorticity ω + βy isnearly constant (and for which the vorticity decreases approximately like − βy across the basin), separated by small areas of abrupt increase of the vorticity(or potential vorticity). For the velocity profile, homogenized potential vorticityareas correspond to westward jets, with a local quadratic velocity profile witha curvature approximately equal to β , while the jumps in vorticity give riseto cusp for the velocity close to eastward jets. This phenomenology has beenobserved in many simulations for barotropic flows and is actually analogous tothe one observed on Jupiter (see figure 1). The top panels of figure 2 showa state with two alternating jets. There are actually four jets: two eastwardjets with local maxima of the zonally averaged velocity and two westward jetswith local minima of the zonally averaged velocity. We stress also that, onthe Hovm¨oller diagram, the jet extrema are located on the black area (close tozero vorticity lines). The eastward jets are the black areas with negative (blue)vorticity on the south and positive (red) vorticity on the north. Related to thePV staircase phenomenology the black areas for eastward jets are very thin line,while the black areas for westward jets have a broader extent. In the doublyperiodic domain of this study, the number of eastward jet has to be equal to thenumber of westward jets, and jets will come by pairs. We will call a state with K alternating jets, a state with K eastward and K westward jets. As stressed in Dritschel and McIntyre (2008), assuming the PV staircase phe-nomenology with a westward jet curvature of order β is sufficient to roughlydetermine the number of jets: its order of magnitude is then given by theratio of the domain size divided by the Rhine scale, or equivalently in our non-dimensional units √ β . Figure 2 indeed shows that increasing β increases thenumber of jets.Compatible with this phenomenology, one expects to see transitions from K alternating jets to a state with K + 1 ones when β is increased. As there is nosymmetry breaking in this process, one may expect these transitions to be firstorder ones with discontinuous jumps of some order parameters. In situationswith discontinuous transitions when an external parameter β is changed, oneexpects for each bifurcation multistability ranges ( β , β ) in which two (or more)possible states are observed depending on the initial conditions, and for whichextremely rare transition from one state to another may occur due to either ex-ternal or internal fluctuations. Such a bistability phenomenology has first beenobserved in Bouchet and Simonnet (2009). The aim of this paper is to go much10
00 1000 1500 2000 2500 3000 350000.10.20.3 | q | αt | q | Figure 3: Lower panel: timeseries the moduli of the first zonal Fourier com-ponents ( | q n | with q n = | Ω | R Ω ω ( x, y ) e iny d x d y ) for n = 2 (red) and n = 3(black). The corresponding Hovm¨oller diagram for the zonally averaged vortic-ity is shown in the upper panel ( α = 1 . · − and β = 5 . , n ), defined as q n = | Ω | R Ω ω ( x, y ) e iny dx dy , for n = 2 and n = 3. Their moduli | q n | are relevantorder-parameter, for instance | q | will be large for a two-alternating jet state andsmall otherwise. In figure 3, one observes five spontaneous 2 → → turnover times. Thefluctuations of the three alternating jet states are noticeably larger than theones for the two alternating jet states. This corresponds to internal dynamicsbetween several three alternating jet states, as will ne explained in section 4.We finally perform an hysteresis experiment by varying adiabatically (veryslowly) the parameter β at rate ˙ β = α , in the range β ∈ [4 , α = 5 · − . Different numerical experiments with independent real-izations of the noise are performed. The result is shown in figure 4. This figuredisplays the modulus of | q | as a function of β in three of these experiments:the upper branch of | q | corresponds to the two alternating jet state and thelower branch to the three alternating jet states. It shows bistability in the range β ∈ [6 , T . − : withthis slowly varying β , transition events are not uncommon. Performing direct11 β | q | Hysteresis α =5.10 −4 Figure 4: We performed four independent hysteresis experiments, by increasingand then decreasing slowly the value of β . For the four experiments, each ofthe four curves represents | q | , the modulus of zonal Fourier component withwavenumber two, as a function of β ( α = 5 · − and ˙ β = α/ α will not show transitions however. Thereason is simply that, as explained in the next sections, the probability of seeingsuch transitions becomes then much too small. Based on the hysteresis curve, wechoose β = 8 as a reference value for future studies, unless otherwise specified. We discuss now the spontaneous 2 → → → → | q | , | q | and | q | . The remarkable result is that these trajectories seem to concentratein the phase space. The blue tube in figure 3 contains 80% of the 2 → → → → → → . Note that it is a veryrough estimate based on a single AMS realization. We did not search for moreprecise estimates here, as we are interested mostly on transitions dynamics. Inprinciple, one should use more realizations with different reaction coordinates(see Appendix B).The middle-left panel of figure 5 shows that the nucleation of a new pairof alternating jets proceeds through an inversion of the velocity curvature ofone of the parabolic westward jet in the exact middle between the two adjacenteastward jets. It takes the form of a small bump in the zonal velocity profile.This structure indeed nucleates from two bands of positive and negative vorticityinside a narrow band of zero vorticity around the westward jet. Such a nucleationis highly improbable, as the new vortex bands are initially too narrow to bestable. Then just by chance, one can have these band surviving long enough togrow until it reaches a stable state. However, once a critical size is reached, theband becomes stable and grows in size. The three jets slowly equilibrate andseparate apart. The new jet can relax either to the north or to the south of theinitial parabolic jet where it has nucleated. We will describe this in more detailsin the next section, by studying saddle zonally averaged zonal velocity profiles.In the following, we will call a transition 2 →
3, a nucleation .The lower-right panel of figure 5 shows how two jets can merge, which canbe interpreted as the disappearance of one pair of alternating jets. The phe-nomenology is again rather simple. Interestingly, if one reverses both the timeand the sign of the velocity, one would indeed obtain a nucleation from a cuspy(westward) jet. Moreover, the resulting two alternating jet state just after thetwo eastward jets have merged, is slightly asymmetric (see zonally averagedzonal velocity for the 3 → /α (not shown). We call the jet merging 3 → coalescence . From a fluid mechanics point of view, it is very natural to understand if therare transitions between attractors are related to hydrodynamical instabilities.It is clear from all the Hovm¨oller diagrams, on figures 2 and 5 and on movies ofthe dynamics (not shown), that the time scale for the dynamics is of order 1 /α ,both during a nucleation of new jets or during a coalescence. This is also truefor all other transitions we have observed. This phenomenology thus excludes afast instability in the barotropic equation (2), developing on time scales of orderone. As a consequence, if some instability occurs in equation (2), the flow hasto be marginally unstable with an instability rate scaling at most like α when α ↓ jets2 jets → → q q q α t y α t y yZonal velocity U yZonal velocity U → → Figure 5: The upper panel shows the level set of the reactive trajectories in thespace | q | , | q | , | q | . The blue and red reactive tubes contain 80% of the 2 → → α = 10 − and β = 8. Middle left and right panels show the zonally averaged zonal velocity asa function of latitude, 3 curves at 3 different times (blue first, then red, thenblack), during the 2 → → RK = U ′′ − β is positive ( α = 2 . − and β = 5 . RK = U ′′ − β , where U ′′ is the second derivativewith respect to y of the zonally averaged zonal velocity U . It is well known that,for a zonal flow U ( y ), a necessary condition for an hydrodynamic instability tooccur is that the Rayleigh–Kuo criteria changes sign (H.L.Kuo, 1949). Figure 6is a Hovm¨oller diagram of the sign of the Rayleigh–Kuo criteria: we draw a blackpoint each time the Rayleigh criteria is negative. This picture is drawn from anexample of a 2 → < αt < .
3) followed by a free dynamics (Eq. (2) without selection)(1 . < αt < α = 2 . − and β = 5 .
5. The eastwardjets are localized inside the thin white bands delimited by the area with denseblack dots. The westward jets are localized in the area with lots of scattered15lack dots, but with a much lower density than at the flank of eastward jets.The picture then shows subsequently: a nucleation of new jets at the locationof a westward jet (0 < αt < . . < αt < < αt < . . < αt < . . < αt < . < αt < U looksvery much like the one of the two jet state. On this curve, we see a characteristicshape for eastward jet velocity profile, with a characteristic cusp shape. At theedge of the eastward jet, the cusp is associated with a very strong negativecurvature and thus a very low value of the RK = U ′′ − β , explaining the whiteband without black dots on figure 6. On each flank of the jet, the velocity profilehas a concavity change, where RK change sign very often, explaining the twodense black dot bands surrounding the eastward jet. On all the other parts ofthe flow RK remains much of the time positive, except close to the westward jet,where U ′′ is very close to β and RK changes sign rather often. This explainsthe less dense band of black dots close to the westward jet edge. We stressalso that the shape of eastward jet cusps seems to depend significantly on thedissipative mechanism: indeed figures 2 and 5, show that for larger values of α for which hyperviscosity effects are negligible the cusp is rather rounded, whileother simulations with smaller values of α (not shown), for which hyperviscosityeffects become important, the velocity profile has a strong concavity change onthe flank of the jet. In a recent theoretical work Woillez and Bouchet (2019), ithas been established that: i) eastward jets should form a cusp, with a width oforder 1 /K where K is the typical random force wavenumber, and which specificprofile depends on the dissipative mechanism, ii) the mechanism preventing thegrowth of westward jet is a marginal instability of the velocity profile appearingwhen the positive curvature of the jet is of order β . Our numerical observationsare thus in line with this theoretical work.Using figure 6, we now describe the jet nucleations and coalescences, fromthe point of view of the Rayleigh–Kuo criterium. First we observe that duringthe nucleation processes, we do not see any of sign change for the Rayleigh–Kuocriteria that would qualitatively differ from the stable jet situation. Second,during the coalescence, a black band appears at the level of the westward jetsqueezed between two eastward jets, just before the coalescence. In this area,although we cannot conclude from this figure, the Rayleigh–Kuo criteria maychange sign. As the Rayleigh–Kuo criteria is only a necessary condition forinstability, we can not conclude from this figure if a very localized and very lowrate instability may occur during the coalescence. But such an instability dooccur, its rate has to be the small as the typical time scale for the change of thestructure, is very long, of order 1 /α . 16e thus conclude that a barotropic adjustment mechanism actual also takesplace during the periods of ”marginal instability” that lead to coalescences andnucleations of jets. In the following, we will use the language of dynamical system theory plusweak noise, in order to describe the phenomenology of the attracting statesand the transition dynamics between them. This dynamical system discussionapplies to the slow effective dynamics of the zonally averaged zonal velocity U ( y, t ) = R π u ( x, y, t ) dx . Indeed, in the regime of small values of α , as arguedin Bouchet et al. (2013), we expect the slow relaxation dynamics of zonal jets tobe described by a deterministic effective equation that can be obtained througha second cumulant closure. Moreover, the transitions from one attractor toanother are expected to be triggered by a noise, due to the fluctuations ofthe Reynolds stress divergence, which properties can be computed using largedeviation theory as explained in Bouchet et al. (2018). This is equivalent to theeffective dynamics ∂ τ U = F ( U ) + √ ασ ( U, τ ) , (3)with the rescaled time τ = αt and the average Reynolds stress divergence F ( U )and a white noise σ ( U, τ ) related to the fluctuation of time averaged Reynoldsstress.In the following, we directly verify the consequences of such a picture ofa dynamical system plus weak noise for the effective dynamics of U ( y ). Theattractors and saddle points should be thought as attractors and saddle pointsof this dynamics, as directly observed in fully turbulent direct numerical simu-lations, or fully turbulent simulations using the AMS algorithm. A saddle pointis thus a zonally averaged zonal velocity profile U ( y ) which is stationary for theeffective dynamics ( F ( U ) = 0) and unstable. In a fully turbulent simulationthis state is characterized by the fact that it belongs to a transition trajectory,and that when adding small perturbations it can either relax to one attractoror to the other depending on the perturbation, as explained better in appendixB.c.For the value β = 8, the hysteresis curve in figure 4 suggests that the systemis multistable with at least two attractors. For dynamical systems with weaknoise, in the simplest cases the transitions from one attractor to another gothrough a saddle point which belongs to the common boundary of the twobasins of attraction. In more complex situations, the transitions could ratherproceed through more complex structures, for instance saddle limit cycles orsaddle strange attractors. As we describe in the following, we have detectedpoint saddles points only for the dynamics of the 2 → → y0 1 2 3 4 5 6−0.3−0.2−0.100.10.20.3 N N ⋆ Figure 7: Zonally averaged velocity as a function of latitude, U ( y ), for thenucleation saddles. Each black curve is obtained by the ensemble average of theoutput of about 800 AMS runs ( α = 10 − and β = 8). Please note the slightasymmetry of each of the saddles, close to the westward tip of the nucleatingjet. For each value of y , we have superimposed the histograms of the U values ofthe 800 AMS states (light blue), in order to feature the spread of the ensemblevalues.or to the south. This is illustrated by the presence of two different ensemblesin the AMS reactive trajectories giving two different saddles. The differenceis small however and it is possible that in other parameter regimes, these twosaddles collapse into a single one at the exact middle of the two westward jets.This is discussed in the next sections. Note that, in addition, one must add twoother saddles related to the nucleation at the other westward jets. Even more,it is possible that one observes direct transitions 2 → → β . This distance increases as β decreases and decreases for larger values of β . This is consistent with the idea that for a regime where only the two-jetstates dominate, say for smaller values of β , the three-jet states will be easily18 C Figure 8: Zonally averaged velocity as a function of latitude, U ( y ), for thecoalescence saddle. Each black curve is obtained by the ensemble average of theoutput of about 200 AMS runs ( α = 10 − and β = 8). For each value of y , wehave superimposed the histograms of the U values of the 800 AMS states (lightblue), in order to feature the spread of the ensemble values.destabilized though the coalescence process. One may ask, why only one saddleis found ? The reason is simply related to the initial conditions used in theAMS. In fact, there are as many coalescence saddles as nucleation saddles. Thispoint is clarified in section 4. Figure 5 strongly suggests the existence of an instanton-driven dynamics. Bythis, we mean that jet nucleation or coalescence events are concentrated arounda predictable path, called instanton. These paths are composed of a rare noisedriven fluctuation path from the two alternating jet attractor toward a saddlepoint (for instance), then of a typical relaxation path from the saddle point tothe three alternating jet attractor. If our hypothesis is correct, as α → α .The important point is that two different saddles are involved in the 2 → → J and J respectively. In figure 9, the nucleation19nstanton is the red curve starting from J and going to N together with itsrelaxation part (black curve) from N to J . Similarly, the coalescence instantonis the red curve starting from J to C and its relaxation part (black curve) goesfrom C to J . The relaxation curves (black curves) correspond to trajectorieshaving zero action in the asymptotic limit since they correspond to the naturaldynamics of the system. Figure 7 shows that the left schematic of figure 9 hasmoreover a mirror equivalent shown to the right. The mirror states are denotedhere with a star for convenience.Even more, to each instanton is associated an anti-instanton going throughthe same saddle. The pair instanton anti-instanton forms a so-called figure-eightpattern in phase space.The time asymmetry does not only translate into a geometrical asymmetryof the pair instanton anti-instanton, but the transition probability of corre-sponding reactive trajectories is different as well. It appears that in this model,the asymmetry is very strong especially in terms of transition probabilities: the“anti-transition” probabilities are several order of magnitude smaller than thetransition probabilities. For instance, the ”anti-coalescence” instanton corre-sponds to a path going from J to C (thin red curve) and then from C to J (relaxation black curve). It requires an inversion of the curvature on the top ofone of the eastern cuspy jets followed by a separation of the two small jets whichare created off the cusp. Such events are of course very rare but can be detectedby AMS. The probability of an inversion of the cusp curvature is estimated byAMS to be O (10 − ). However, the nucleation of an eastward jet leading to theconfiguration shown in Fig. 8 has a considerably smaller probability, estimatedby AMS to be less than O (10 − ).We describe in the next section the internal effective dynamics of the statewith three alternating jets. The previous results indicate the existence of at least two types of states havingthree alternating jets, due to different possible nucleations (see figure 7) orcoalescences (see figure 8). It was previously referred to as a mirror symmetrywithout being too precise on what symmetry was involved. We show next thatthe reality is more complicated, with the description of several different stateswith three alternating jets, and a slow dynamics between them. This explainsin particular the large amplitude of the fluctuations of the three alternating jetstates, as seen in figure 3.
In order to describe this structure we need a more precise description of thejet profiles (the zonally averaged zonal velocity). If one measures precisely the20 C C ⋆ J ⋆ C J J NN ⋆ J J Figure 9: Schematic representation of the instantons involved in the transitions2 → →
2. Left figure is a simplification of the right figure, ie without repre-senting the dual mirror pair of states. J is the stationary two-jet state, J and J ⋆ are the two mirror-symmetric three-jet states obtained from the two nucle-ation saddles called N and N ⋆ (see figure 7). The coalescence saddles are called C and C ⋆ (see figure 8). The black arrows represent the relaxation trajectoriesfrom the saddle to the attractors (with zero action; see Appendix A). The redcurves correspond to the reactive part of the instantons (with nonzero action).The colored area inside the bold curves is what is shown in instanton tubes offigure 5, that is half of the actual instanton structure.21istance between the jet tips, one observes some slight asymmetries in the threealternating jet states. Indeed, the attractors do not equilibrate to an equidistantjet solution, see for instance the three alternating jet states in figure 3. Thisasymmetry is robust when β is changed. Based on this remark, it is relevantto characterize these states by the distance between the eastward jet tips forinstance. These distances give three numbers we call σ , σ , σ . Our conventionis to divide these quantities by the dimensionless domain width 2 π so that0 ≤ σ i ≤
1. Due to the periodicity in y , one has the constraint σ + σ + σ = 1 . Given this constraint, in the following, we prefer to consider the most economicalrepresentation and use only two of these quantities, say σ and σ . Looking atthe Hovm¨oller diagram of Figure 3 for instance, one observes some slow drift ofthe solution in the meridional direction. We therefore choose one of the easternjet as a reference moving frame, and compute σ and σ from this referencejet. σ is the distance between this reference jet and the next eastward-jet tip,northward. σ is the distance between this new jet and the next eastward-jettip, northward again. This procedure is for measurements only and requires toalways track the same reference jet at any time. A state with three alternatingjets will be thus described by ( σ , σ , σ = 1 − σ − σ ).It appears that the values taken by σ , σ and σ are always very close tothree fixed values σ ∗ , σ ∗ and σ ∗ which are all different from 1 / σ ∗ , σ ∗ , σ ∗ ) thereare two other symmetric configurations by meridional translation: ( σ ∗ , σ ∗ , σ ∗ )and ( σ ∗ , σ ∗ , σ ∗ ). However there exists also three other states that can not berecovered by meridional translation from the first one ( σ ∗ , σ ∗ , σ ∗ ), ( σ ∗ , σ ∗ , σ ∗ )and ( σ ∗ , σ ∗ , σ ∗ ), but can be recovered if we use a mirror symmetry. These 6states are in correspondence with all the permutations of { , , } , three evenpermutations and three odd permutations.For clarity, figure 10 shows a schematic of these configurations. The group ofpermutation is called S and possesses 3! = 6 possible permutations: { , , , , , } .We now show in the next subsections a remarkable result: the symmetric group S is indeed realized by the natural dynamics of the system. ( σ , σ ) plane We perform a direct numerical simulation with some long time integration ofequation (2), and we measure the distances σ and σ . This probability densityfunction of ( σ , σ ) is shown in Fig. 11. The time averaged values obtainedfrom this dataset gives < σ > = σ ∗ ≈ .
340 and < σ > = σ ∗ ≈ .
226 withfluctuations for σ and σ being of the same order, of about ± . < σ > the value of σ averaged over a sliding time windows over a small enoughperiod of time. It reaches three distinct values < σ > ∈ { . , . , . } corre-sponding to the three values { σ ∗ , σ ∗ , σ ∗ } . A closer inspection of the quantities22 σ , σ , σ ) ( σ , σ , σ )( σ , σ , σ ) ( σ , σ , σ )( σ , σ , σ )( σ , σ , σ ) Figure 10: Sketch of the zonal velocity as a function of latitude, U ( y ), gen-erated by the symmetry group for the six different configurations, with three-alternating jet. We have amplified the meridional asymmetry of the jets forbetter visualisation. σ , σ , shows that the system is indeed wandering around four distinct config-urations denoted a, b, c, d in in figure 11. When one integrates for a sufficientlylong period of time, the number of possible configurations is in fact 6, in corre-spondance with the 3! permutations of the symmetric group S .The probability density function of the full set is shown in figure 13. Oneclearly sees in light blue some traces of transitions between pairs of adjacentstates. As explained in the next subsection, all states are indeed connected bytwo type of transitions, one toward each adjacent state. The transitions observed must involve saddles with two equally distant jets andare therefore located on the white dot lines of figure 13. Let a state be writtenas ( σ , σ , σ ) and let i, j, k ∈ { , , } such that σ i < σ j < σ k . Then type-I transitions correspond to transitions where the two smallest distances havebeen switched so that the new state has σ j < σ i < σ k . Indeed, it can beseen in figure 13 by looking at the visible blue traces. This is a consequence ofthe two possible nucleations of the form ( − ǫ, + ǫ, ) and ( + ǫ, − ǫ, ).Type-II transitions follow the same type of rule, where instead, the two largestdistances are involved. In the example above, the new state will be such that σ i < σ k < σ j . It is clear that there are no possible transitions where the smallestand largest distances are switched, as it would require to go through the perfectlysymmetric state ( , , ): the system prefers to visit two configurations before,through one type-I transition, and one type-II transition. By symmetry, we23 −3 σ σ Figure 11: Probability density function of ( σ , σ ) for one of the states withthree alternating jets ( α = 10 − and β = 8).24
00 1000 1500 2000 2500 30000.20.250.30.350.40.45 α t a b c b a d a d Figure 12: Timeseries of σ (black curve) together with σ (blue curve) and σ (red curve) with the time axis rescaled by α = 10 − .25igure 13: Probability density of ( σ , σ ) featuring the six different states withthree alternating jets. The three dotted lines correspond to states with twoequally distant jets, they intersect at ( σ , σ , σ ) = ( , , ). Figure 11 is azoom of the states in the lower-left sector corresponding to configuration a .The blue traces inbetween the states feature reactive trajectories statistics.26 Type I saddles
Figure 14: The three type-I saddles ( σ, σ, − σ ) , ( σ, − σ, σ ) and (1 − σ, σ, σ )are located at the intersections between the dot lines and the blue traces ofFig.13 with σ ≈ .
30. The corresponding zonally averaged zonal velocities as afunction of latitude, U ( y ), which were obtained by direct numerical simulation,are represented. Type II saddles
Figure 15: The three plots show the zonally averaged zonal velocities as afunction of latitude, U ( y ), for the three type-II saddles ( σ, σ, − σ ) , ( σ, − σ, σ )and (1 − σ, σ, σ ) with σ ≈ .
37. They have been obtained by direct numericalsimulations.conclude that there are a total of three type-I saddles and three type-II saddles.Those saddles are computed as discussed in appendix B.c. Type-I saddles havetwo equally-small distances and are shown in figure 14. Type-II saddles havetwo equally-large distances and are shown in figure 15.
We first mention that one can represent a two-jet state in our ( σ , σ ) parameterplane in different ways as there is one missing jet. We decide to choose the onewhich is consistent with the dynamics in the limit where one jet appears or dis-appears. One must identify first the four corners (0 , ≡ (0 , ) ≡ ( , ≡ ( , )as a single state. During the nucleation process, the initial perturbation beforereaching the nucleation saddle is a jet of the form ( , , ) (plus permutations)27nd the nucleation saddle is of the form ( + ǫ, − ǫ, ) and its five symmetricpartners making a total of six different nucleation saddles. One can measurethe numerical value of ǫ from figure 7, it gives ǫ ≈ .
01. The same rule appliesfor the coalescence saddles, there are a total of six coalescence saddles. Figure16 is a summary of all the results.Despite its simplicity, one should remark that there is, in addition, sub-internal low-frequency variability within a single of the six states with threealternating jets as observed in figure 11. This suggests the existence of additionalsaddles possibly within the white S shape visible in figure 11. We do not showthe anti-instantons 2 → → S symmetry-breaking bifurcations of the effective dynamics have occurred, one related totype-I transitions and one to type-II. Due to that, internal and complex low-frequency variability can induce rather large fluctuations of the states with threealternating jets. The scenario is in fact similar to the typical deterministic one,where global bifurcations like homoclinic and heteroclinic bifurcations force thesystem to have large-amplitude oscillations. One can readily anticipate thatthere are simpler regimes unfolding from the perfectly symmetric state. Chang-ing β does not seem to change radically the internal structures shown in figure16, except for the saddles involved in the transitions 2 → → β is changed. It is plausible that, changing for instancethe aspect ratio of the domain, could reveal situations where the state ( , , )is an attractor. We did not explore these interesting possibilities. We now address the important question of the dependency of the attractorstructure and the transition rate between attractors, with respect to the param-eter α . From equation 3, as a direct consequence of Freidlin-Wentzell theory(Freidlin and Wentzell, 1984), a large-deviation principle for the full system canbe inferred. In particular, an Arrhenius law must be present. One can rephrasethis result even more radically by saying that the effective dynamics describedin figure 16 and obtained for α = 10 − is indeed valid in the limit α → α . Values of ( σ ∗ , σ ∗ , σ ∗ ), for different values of α , are reported intable 1. We conclude that the three alternating jet attractors, corresponding tofigure 11, are slightly more asymmetric when α is decreased. In particular, the28 N ⋆ N CC ⋆ Figure 16: Effective phase space structure in the ( σ , σ ) plane, including all S symmetric states. The red squares are the three-alternating jet states, theblue squares correspond to the two alternating jet symmetric states. The bluediamonds are asymmetric two-alternating jet states obtained just after the co-alescence process. Red circles are the coalescence saddles, blue circles are thenucleation saddles, black circles are type-I internal saddles and dark blue cir-cles are type-II saddles. The middle square is the (unstable) three alternatingjet symmetric state. The states just before the nucleations are shown as smallwhite squares. The positions correspond to observed numerical values. The dotlines are the internal mirror transitions responsible for metastability and thesolid blue and red lines are the 2 → → σ ∗ and σ ∗ for the jet attractors shown in figure 11 asa function of α . α < σ > < σ > std8 · − · − · − · − · − · − · − · − closest jets become closer as α decreases, whereas σ stays almost unchanged.The effects of the resolution and viscosity are probably becoming more impor-tant so that one cannot investigate further the limiting behavior.We now discuss our last experiment. Using the AMS algorithm, we computefor different values of α , the mean time it takes to observe 2 → mean return time or mean first transitiontime and essentially scales like the inverse of the transition probability. Note thatwe have used a different value of β = 5 . β = 8. The dependency on β is in fact irrelevant, as the result we discuss here is mostly insensitive to β . Notealso that although the probability of 2 → β ,the nucleation phenomenology remains unchanged. The result is shown in figure17 and suggests the existence of an Arrhenius law of the form T ∝ exp (cid:0) ∆ Uα (cid:1) ,which is expected from large deviation theory.If one decreases α even further, the scaling observed in figure 17 does notcorrespond anymore to the simple Arrhenius law (not shown). This is mostprobably because viscosity ν is not negligible anymore for so small values of α . Then the working hypothesis ν ≪ α in order to obtain T ∝ exp (cid:0) ∆ Uα (cid:1) isno more valid. To support this hypothesis, one can observe that the meanenergy expected to be equal to one for negligible values of ν , drops to valuessignificantly smaller than one when one uses such small values of α . Usingsmaller value of ν to test further the Arrhenius law would require to makemuch more costly numerical simulation with a higher resolution, which wouldbe extremely difficult.It is important to understand that such scaling law cannot be revealed by di-rect numerical simulations as illustrated by the table 2: The amount of requiredCPU time to study the Arrhenius law with a direct numerical simulation in-creases exponentially as α decreases, whereas using the AMS thr required CPUtime increases only linearly (see Appendix B). Using the rare event AMS algorithm, we are able to obtain a complete statisti-cal description of the dynamics of planetary jets in a simple barotropic quasi-30
33 1111 1667 22222.1e+062.0e+083.5e+093.5e+11 1/ α A ve r a g e t r a n s i t i on t i m e Figure 17: Logarithm of the mean first transition time T as a function of α forthe 2 → N = 1000, and two to three independent algorithm realizations were performed.Error bars are a crude approximation of the AMS variance based on the differentrealizations.Table 2: CPU time (d:day,y:year) needed to obtain about 1000 2 → α . α AMS DNS1 . · − . · − . · − ∼
51 y0 . · − ∼ α reveals a scaling law clearly compatiblewith an Arrhenius law. This implies that the observed phenomenology is indeedinstanton-driven. As a by-product of our rare event algorithm, we are also ableto compute the saddles of the effective dynamics.We next focus on the internal dynamics of solutions having three eastwardand westward jets. Surprisingly, it shows a striking example of internal metasta-bility. The state with three alternating jets appears to have jets not at equaldistance to each other. We show that there are in fact six distinct configu-rations which are governed by the symmetric group of permutation called S .The consequence is a complex dynamics with low-frequency variability and largefluctuations. This behavior originates from at least two S symmetry-breakingbifurcations of the effective dynamics associated with two different types oftransitions depending on the distance between the jets.We stress that these results are generic. For instance, the observed internalmetastability of the three-jet states should appear as soon as there is someasymmetry in the jet positions and shapes. It comes from deeper algebraicconstraints link to the symmetries of the solutions. To this respect, changing β should not alter this picture but rather has an impact on the stability of thesestates. As discussed before, changing α to even smaller values will not changethis picture either.One expects each metastable state to concentrate on a well-defined singleconfiguration. One may also ask whether the use of a doubly-periodic domainhas an impact on these results or not? As far as transitions between states withtwo and three alternating jets are concerned, the nucleation and coalescencephenomenology is mostly local and should prevail in more general domains. Theinternal symmetry-breaking dynamics of the state with three alternating jets islikely to be modified in more visible ways. However, we expect this mechanismto be robust, since it is the relative distance between the jets which plays arole. As a basis for these assertions, we demonstrate that this phenomenologyis preserved when one uses lateral meridional walls in appendix C. As mentionedin section 4, it is plausible that changing the aspect ratio to a small value hasa more radical impact on the internal dynamics. Effective bifurcations, eithersupercritical or subcritical, are potentially controlled by the domain aspect ratio.For instance, one can possibly obtain states with symmetric three alternatingjets which are stable configurations. In fact, all the visible structures shown inthis work should unfold from a single symmetric state. It is remarkable that allthe classical bifurcation concepts including normal form theory are recoveredhere in a statistical sense through large-deviation principles.One consequence of this approach is the possibility to infer many results32n the dynamics involving n jets. It is reasonable to think that in the generalcase there are n ! distinct configurations controlled by the symmetric group S n .The situation becomes rather involved however as saddles can be of many moredifferent types instead of type-I and type-II for n = 3.We stress that it is very difficult to obtain such results using classical tools.The rare-event AMS algorithm is performing exponentially better than anynaive Monte-Carlo methods and is particularly fitted to study large deviations.The price to pay is the difficulty sometimes to control the algorithmic variance.It essentially depends on the choice of a good reaction coordinate.We have not discussed one important aspect of the dynamics, which is calledjet migration , and is observed in advanced primitive-equation models of Jupiter(Williams, 2003). It takes the form of slow drift of the whole jet system whichcan be either northward or southward. Our mechanism of metastability doesnot explain such behavior since the relative distance between the jets must bethe same. We do however observe jet migrations either to the north or to thesouth in our simple model (see e.g. Fig. 3). The AMS algorithm can be usedin such a context by simply changing the reaction coordinate and is expectedto provide powerful insight on jet migration as well.In the future, we would like to consider more realistic models and, in particu-lar, two-layer baroclinic models (Phillips, 1951) using a small aspect ratio. Thesemodels are rather popular for the study of the midlatitude Jovian atmosphere(Williams, 1979; Haidvogel and Held, 1981; Panetta, 1993; Kaspi and Flierl,2007). Interestingly, one anticipates the same transition phenomenology thanthe one studied here, at least in square domains. The Hovm¨oller diagram ofPanetta (1993), Fig. 10 for instance is very suggestive as one observes severalsmall coalescences at the beginning and a well-defined 3 → → acknowledgments Eric Simonnet acknowledges support from CICADA (centre de calculs interac-tifs) of University Nice-Sophia Antipolis. The computation of this work werepartially performed on the PSMN platform of ENS de Lyon. The research lead-ing to these results has received funding from the European Research Councilunder the European Union’s seventh Frame- work Programme (FP7/2007-2013Grant Agreement No. 616811). During its last stage, this publication wassupported by a Subagreement from the Johns Hopkins University with fundsprovided by Grant No. 663054 from Simons Foundation. Its contents are solelythe responsibility of the authors and do not necessarily represent the officialviews of Simons Foundation or the Johns Hopkins University.33
Instantons, large deviations
Let us consider a stochastic partial differential equation written as du = F ( u ) dt + √ αdW t , where u ≡ u ( x , t ) , x ∈ R d , u ( x ,
0) = u ( x ) is an initial condition and W t isa Wiener process. We define V = V ( u , u , T ), the set of trajectories startingfrom u and ending at u at time T . A path integral approach (initiated byOnsager and Machlup, 1953) allows to formally express the transition probabil-ity from u at time t = 0 to u at time T asPr ( u , T | u ,
0) = Z V exp (cid:18) − A T [ u ] α (cid:19) D [ u ] . (4)Here D [ u ] is a measure over paths which can be defined rigorously, and A T isthe so-called action: A T [ u ] = 12 Z T (cid:12)(cid:12)(cid:12)(cid:12) dudt − F ( u ) (cid:12)(cid:12)(cid:12)(cid:12) dt. The saddle-point approximation of (4) in the limit α → T → ∞ , the large-deviation result − lim α → α log Pr = inf u ∈ ˜ V A ∞ [ u ] , with ˜ V = { u ( t ) , lim t →−∞ u ( t ) = u , u (0) = u s } . Here, u s is a saddle point at the boundary of the basin of attraction of u .Generically, there is a unique minimizer of the action A = A ∞ which is calledan instanton path , often referred to as the most probable path going from u tothe saddle u s in the limit of vanishing noise. Remarkably, Freidlin and Wentzell(1984) have shown that the saddle approximation of the path integral is indeedvalid provided u is a fixed isolated attracting point of the deterministic systemwith α = 0. Let v be a solution in the basin of attraction of u , the instantonpath is in this case simply the trajectory of the deterministic dynamics whichsatisfies dudt − F ( u ) = 0 with u (0) = v . It is the unique minimizer of the infinitetime action A which is trivially zero. One calls these solutions, relaxation paths(Bouchet et al., 2014). Of particular interest is the case where the deterministicsystem is a gradient system with F ( u ) = − δVδu . It is possible to define the time-reversed path of the relaxation path going from u to v and it is called fluctuationpath. Fluctuation paths have of course the same geometrical support thanrelaxation paths (see also Vanden-Eijnden and Heymann, 2008). One can easilydeduce that a fluctuation path is indeed the unique minimizer of the infinite timeaction which becomes equal to V ( v ) − V ( u ). The classical Arrhenius law forgradient systems is thus obtained , namely − lim α → α log Pr = V s − V , where V s is the value of the potential V at the saddle u s .34 Adaptive Multilevel Splitting algorithm (AMS)
The idea of splitting indeed took its source during the Manhattan project withJ.Von Neumann (unpublished) and was sometimes referred to as Von Neumannsplitting (Kahn and Harris, 1951; Rosenbluth and Rosenbluth, 1955). Theseideas were more or less forgotten until the 90s where they start to be ratherpopular in chemistry, molecular dynamics and networks. These algorithms arenow better understood from a mathematical point of view and have becomean important area of research in probability (P. Del Moral, 2004). We use herethe modern adaptive version of these splitting techniques which is proposed byC´erou and Guyader (2007). It is much more robust and versatile than the clas-sical versions used previously. The general idea is to decompose a very smallprobability into a chain of products of (larger) conditional probabilities. Assuch, one mimics the evolution of species by performing Darwinian selectionson the system dynamics, in a controlled (unbiased) way. These types of algo-rithms essential perform a large number of mutations and selections (branch-ing) by cloning the system dynamics. This is the reason they are sometimescalled genetic algorithms although they bear several different names like multi-level splitting, go-with-the-winner, rare event or even large-deviation algorithms.Note that they must not be confused with the well-known family of importancesampling techniques (see Ebener et al., 2018).
B.1 Algorithm description
One needs a Markovian model ( X x t ) t ≥ with X = x ∈ H , where H is an ad-hocfunctional space. Typically, X t corresponds to a stochastic partial differentialequation (PDE) in H . Let us consider two arbitrary sets A and B such that A ∩ B = ∅ . The goal is to estimate the probability p to enter the set B ,starting from the initial condition x , before returning to A . Let τ C ( x ) = inf { t ≥ , X x t ∈ C} . The probability p translates to p = p ( x ) = Pr( τ B ( x ) < τ A ( x )).A reactive trajectory is a particular realisation of this probability. We notethat p ( x ) = 0 if x ∈ A and p ( x ) = 1 if x ∈ B . The function p : x p ( x )is called either committor or importance function or equilibrium potential inmathematics. It can be shown to solve the backward Fokker-Planck equationfor diffusive processes (e.g., E and Vanden-Eijnden, 2006). Solving this PDE isout of question in large (infinite) dimension. One then needs to define a quantitywhich measures how far a trajectory is escaping from A . There is no uniquechoice and this will depend mostly on the question asked. Let Φ : H → R a chosen function which is called reaction coordinate or observable . One thendefines Q ( x ) = sup t ∈ [0 ,τ A ] Φ( X x t ) . For convenience Φ is renormalized such that Φ( A ) = 0 and Φ( B ) = 1 althoughit is not necessary. We can now give the description of the algorithm in itssimplest version. 35 seudo-code (last particle) : Let N be a fixed number, Φ given and x somefixed initial condition. The quantity b p is an unbiased estimate of p . Initialization : set the counter K = 0 and b p = 1. Draw N i.i.d. trajectoriesindexed by 1 ≤ i ≤ N , all starting at x until they either reach A or B . Compute Q i ( x ) for 1 ≤ i ≤ N . In the following we denote Q i thesequantities since their initial restarting conditions might change.DO WHILE (min i Q i ≤ Selection : find i ⋆ = argmin i Q i and set b Q = Q i ⋆ . Branching : select an index I uniformaly in { , · · · , N } \ { i ⋆ } . Find t ⋆ such that t ⋆ = inf t n Φ( { X ( I ) t ≥ b Q o . Compute the new trajectory i ⋆ with initial condition x ⋆ ≡ X ( I ) t ⋆ until it reaches either A or B .Compute the new value Q i ⋆ ( x ⋆ ). DO K ← K +1 and b p ← (1 − N ) × b p .END WHILEFigure 18 gives an illustration of these algorithmic steps. The performanceof the algorithm depends crucially on the choice of Φ. One can show thatwhen Φ has the same isosurfaces than the committor p ( x ) (which is alwaystrue in 1-D), the number of iterations K ∼
Poisson( − N log p ) (Guyader et al.,2011; Simonnet, 2016). Central limit theorems can now be demonstrated ingeneral situations (Br´ehier et al., 2016; C´erou and Guyader, 2016) and typicallytake the form √ N p − b pp D −−−−→ N →∞ N (0 , | log p | ). The advantage of these algorithmsis the unbiased estimate they provide. One should be careful however as abad choice of Φ may lead to lognormal law with heavy tails and apparent bias(Rolland and Simonnet, 2015). In practice, the larger N the better. A too smallvalue of N not only give results with large variance but also trigger extinction ofspecies by loss of trajectory diversity. One should take advantage of algorithmicvariants which eliminate n ≪ N trajectories at each step instead of just one.The number of iterations K in this case scales like Nn | log p | . Doing so in a parallelenvironment provide much faster results. This is indeed the strategy used in thispaper, where we typically have N = O (10 ) and n = O (100) on n cores witha speed-up scaling linearly like O ( n ). In practice, one should also run severali.i.d. algorithmic realisations with different choices of reaction coordinates andconsider a set of initial conditions x rather than a single initial state (see e.g.,C´erou et al., 2011). B.2 Reactive coordinates
There are many possible options to construct reasonable reaction coordinates.The first option is to consider a low-dimensional projection of the system phasespace, say on zonal Fourier modes. One can also consider EOFs or wavelets basisto better represent the solutions. In our case, we have chosen the Fourier moduli | q | , | q | , | q | for Fig. 5. One then needs to construct a function f ( | q | , | q | , | q | ) ∈ Q < Q < Q new Q x ⋆ A B x Figure 18: A1Sketch of the adaptive multilevel splitting algorithm (last particle algorithm)with N = 3 trajectories. One eliminates the trajectory (trajectory labeled 2in the sketch) which has the smallest maximum value of Q (denoted Q ), andcreates a new one (sketched by a red curve) which is resampled from the initialcondition from one of the others trajectory chosen randomly (here trajectory1), and when it crosses the level Q . 37 such that it describes well the fluctuations of the system. Let M a projectionof the system in the low-dimensional phase space, a natural choice is then toconsider f ( M ) = d A d B if d A ≤ d B and f ( M ) = 1 − d B d A if d B ≤ d A , where d C = dist( M , C ). The function f is equal to zero in the set A and one in theset B . One must be careful as the Euclidian distance may not reflect well thedynamics of the transitions. Moreover if the dimension of the projected phasespace is too low, relevant fluctuations may not be seen at all.Another option is to work in the physical space instead. In our case, it isnatural to work in the 1-D zonal velocity space such as those shown in manyinstances of this work. For the nucleation transitions, the strategy is to detectloss of convexity locally near the parabolic jet. A typical relevant quantity isthen || u − u ⋆ || where u ⋆ is the local convexification of u , say between two easternjets. This quantity typically can detect nonlocal changes of curvature and hasa smoothing effect. In principle the second derivative of u at the parabolicjet would be the right quantity but the corresponding fluctuations can be verylarge and may reflect more subtil changes in the velocity profile even far fromthe jet. Moreover, as the nucleation proceeds, the fluctuations of the curvatureat the newly created jet saturate and become more or less irrelevant. If one hasenough control and understanding on such reaction coordinates, only one scalarquantity is necessary to perform AMS. If it is not the case, one should considera mix of both approaches. Note that one is indeed doing some simplified formof image processing so that deep neural networks can bring interesting newperspectives. B.3 Saddle detection
Once a reactive trajectory reaches a saddle, it should relax, with a probabilityclose to 1, either to the set A or the set B depending whether it is ”before” or”after” the saddle. This property is at the basis of the three methods presented.1 The simplest method requires no extra computations and is fairly accurate.Once AMS algorithm has converged, one should simply record the statesover which the last branching has occurred. The reason is that theseinitial conditions all relax (by the ”natural” dynamics, ie the one withprobability close to 1) to the target set B by the AMS definition whereasprevious states relax to A again by definition.2 The dichotomy method is more precise than method 1, but it requiresextra costly computations as one must reconstruct all AMS trajectories(see Lucarini and B´odai, 2017; Schneider et al., 2007; Willis and Kerswell,2009, for a similar idea). One proceeds by dichotomy on the reactivetrajectories ( X t ) t ∈ [0 ,τ ] by picking initial conditions along the trajectory.For each initial conditions one relaxes the dynamics to decide whether itgoes to A or B . The dichotomy provides a critical time, say t s , and thesaddle is close to X t s . Note that this is correct in probability only, so thatone should perform several realisations, and for many AMS trajectories38n practice, one realisation is enough since the return time associated tothe transition is much larger than the natural relaxation dynamics.3 This method has mainly a theoretical interest. The idea is simply tocompute pairs of instantons, anti-instantons and to look for their inter-sections since they must intersect at saddles and form a so-called figure-eight pattern (see Fig. 9). For reversible dynamics like gradient systems,this approach is irrelevant as the geometric support of an instanton andanti-instanton is the same. In practice, computing anti-instantons is verychallenging. C Permutations in a channel flow
We investigate here the robustness of the S n symmetric group w.r.t. the bound-ary conditions. We conduct an experiment in a zonally-periodic channel withfree-slip boundary conditions, by imposing q = 0 at the north and south wallsand and an additional constraint ∆ q = 0 for the hyperviscosity. We focus on aparameter regime having three jets ( β = 8), and use the same isotropic annularstochastic forcing than in the doubly-periodic case. We generate O (10) i.i.d.realisations starting from the same initial condition q = 0. We obtain threedistinct mirror pairs of solutions w.r.t. the symmetry y → − y (denoted by abar), ie a total of 6 different configurations. No additional states have beenobserved even by increasing the number of realisations. We define σ as thedistance between the first and second eastward jet and σ , between the secondand third eastward jet. The results are shown in Fig. 19 and we call A, B, C thethree states (right of Fig. 19) together with their mirror counterparts ¯ A, ¯ B, ¯ C (upper part of Fig. 19).A close inspection of the numerical values ( σ , σ , σ ) where σ = 1 − σ − σ ,shows that not all permutations are realized. Let us look at the state A withvalues ( σ ,A , σ ,A , σ ,A ) ≈ (0 . , . , . C is simply ( σ ,A , σ ,A , σ ,A ),¯ A = ( σ ,A , σ ,A , σ ,A ) and ¯ C = ( σ ,A , σ ,A , σ ,A ), giving a total of 4 states. Thestate B does not correspond to a permutation of the previous states. It hasvalues ( σ ,B , σ ,B , σ ,B ) ≈ (0 . , . , . σ and σ . Note that the boundary conditions impose the presence of eitherwestward or eastward jets on the walls since q = 0 yields dudy = 0, in particularwe observe that the states have either two westward jets ( A, B, ¯ A, ¯ B ) or oneeastward and one westward jet ( C, ¯ C ). The existence of two eastward jets atthe north and south walls seems to be forbidden for this value of β , it has neverbeen observed.We thus conclude that the symmetric group appears to be broken into twodifferent subgroups. It suggests a more complex general picture where the sym-metries of the problem play a key role in the interaction between the jets. It isplausible for instance that in a four-jet regime, the symmetric group S wouldemerge as a subgroup.One of the experiment is shown in the form of an Hovm¨oller diagram in39igure 19: Probability density function ( σ , σ ) for the 6 different three alternat-ing jet states, in a zonally-periodic channel with free-slip boundary conditions.The computation used about 10 independent noise realizations, each startingfrom rest (zero initial condition). The physical parameters are β = 8 and α = 5 · − . The number of grid points in each direction is N = 400 usingRK4 time scheme for an annular isotropic noise like in the doubly-periodic caseat k f = 12. The diagonal ( , , ) is shown with a white dotted line.Fig. 20. A striking transition from A → ¯ A takes place through a nucleationinducing a 4-jet state quickly followed by a merging back to the mirror 3-jetstate. The solution goes through the perfectly marginally stable symmetry state( , , ) (see middle of Fig. 20). Transitions involving internal saddles have notbeen directly observed and deserve further investigations using the rare eventAMS algorithm. References
Arnold, N. P., E. Tziperman, and B. Farrell, 2012: Abrupt transition to strongsuperrotation driven by equatorial wave resonance in an idealized gcm.
Jour-nal of the atmospheric sciences ,
69 (2) , 626–640.Bakas, N., and P. Ioannou, 2014: A theory for the emergence of coherent struc-tures in beta-plane turbulence.
Journal of Fluid Mechanics , , 312–341.Berhanu, M., and Coauthors, 2007: Magnetic field reversals in an experimentalturbulent dynamo. Eur. Phys. Lett. , , 59 001.Bouchet, F., J. Laurie, and O. Zaboronsky, 2014: Langevin dynamics, large de-viations and instantons for the quasi-geostrophic model and tzo-dimensionalEuler equation. J. Stat. Phys , , 1066–1092.40igure 20: Hovm¨oller diagram of the zonally averaged zonal velocity in a channelflow showing a transition between A and ¯ A .Bouchet, F., J. Marston, and T. Tangarife, 2018: Fluctuations and large devi-ations of reynolds stresses in zonal jet dynamics. Physics of Fluids ,
30 (1) ,015 110.Bouchet, F., C. Nardini, and T. Tangarife, 2013: Kinetic theory of jet dynamicsin the stochastic barotropic and 2d navier-stokes equations.
J. Stat. Phys. ,
153 (4) , 572–625.Bouchet, F., J. Rolland, and E. Simonnet, 2019: Rare event algorithm links tran-sitions in turbulent flows with activated nucleations.
Physical Review Letters ,
122 (7) , 074 502.Bouchet, F., and E. Simonnet, 2009: Random Changes of Flow Topologyin Two-Dimensional and Geophysical Turbulence.
Physical Review Letters ,
102 (9) , 094 504, doi:10.1103/PhysRevLett.102.094504.Br´ehier, C. E., M. Gazeau, L. Gouden`ege, T. Leli`evre, and M. Rousset, 2016:Unbiasedness of some Generalized Adaptive Multilevel Splitting algorithms.
Annals. of App.Prob.
Caballero, R., and M. Huber, 2010: Spontaneous transition to superrotationin warm climates simulated by CAM3.
Geophys. Res. Lett. , , L11 701, doi:10.1029/2010GL043468.C´erou, F., and A. Guyader, 2007: Adaptive multilevel splitting for rare eventsanalysis. Stoch. Anal. and Appl. , , 417–443.C´erou, F., and A. Guyader, 2016: Fluctuation Analysis of Adaptive MultilevelSplitting. Anal.App.Prob. , , 3319–3380.C´erou, F., A. Guyader, T. Leli`evre, and D. Pommier, 2011: A Multiple ReplicaApproach to Simulate Reactive Trajectories. J.Chem.Phys. , , 054 108.C´erou, F., A. Guyader, and M. Rousset, 2019: Adaptive multilevel split-ting: Historical perspective and recent results. Chaos: An Interdisciplinary ournal of Nonlinear Science ,
29 (4) , 043 108, doi:10.1063/1.5082247, URLhttp://aip.scitation.org/doi/10.1063/1.5082247.Clerke, A. M., 1893: A popular history of astronomy during the nineteenthcentury.Dansgaard, W., and Coauthors, 1993: Evidence for general instability of pastclimate from a 250-kyr ice-core record.
Nature , , 218–220, doi:10.1038/364218a0.Dematteis, G., T. Grafke, M. Onorato, and E. Vanden-Eijnden, 2019: Exper-imental evidence of hydrodynamic instantons: the universal route to roguewaves. Physical Review X , , 041 057.Ditlevsen, P., K. K. Andersen, and A. Svensson, 2007: The do-climate eventsare probably noise induced: statistical investigation of the claimed 1470 yearscycle. Climate of the Past , , 129–134.Dritschel, D., and M. McIntyre, 2008: Multiple jets as pv staircases: the phillipseffect and the resilience of eddy-transport barriers. Journal of the AtmosphericSciences ,
65 (3) , 855–874.E, W., and E. Vanden-Eijnden, 2006: Towards a Theory of Transition Paths.
J.Stat.Phys. , , 504–523.Ebener, L., G. Margazoglou, J. Friedrich, L. Biferale, and R. Grauer, 2018:Instanton based importance sampling for rare events in stochastic PDE. arXiv:1812.03543 .Ebener, L., G. Margazoglou, J. Friedrich, L. Biferale, and R. Grauer, 2019:Instanton base importance sampling for rare events in stochastic pdes. Chaos:An Interdisciplinary Journal of Nonlinear Science .Farrel, B. F., and P. J. Ioannou, 2007: Structure and Spacing of Jets inBarotropic Turbulence.
J. Atmos. Sci. , , 3652–3664.Farrell, B. F., and P. J. Ioannou, 2003: Structural stability of turbulent jets. Journal of Atmospheric Sciences , , 2101–2118.Freidlin, M. I., and A. D. Wentzell, 1984: Random Perturbations of DynamicalSystems, volume 260 of Grundlehren der Mathematischen Wissenschaften .Springer-Verlag, New York.Galperin, B., and P. Read, Eds., 2019:
Zonal Jets: Phenomenology, Genesis,and Physics . Cambridge University Press, doi:10.1017/9781107358225.Galperin, B., S. Sukoriansky, and H.-P. Huang, 2001: Universal n − spectrumof zonal flows on giant planets. Phys. Fluids , , 1545.42alperin, B., R. M. Young, S. Sukoriansky, N. Dikovskaya, P. L. Read, A. J.Lancaster, and D. Armstrong, 2014: Cassini observations reveal a regime ofzonostrophic macroturbulence on jupiter. Icarus , , 295–320.Grafke, T., R. Grauer, and T. Sch¨afer, 2015: The instanton method andits numerical implementation in fluid mechanics. J.Phys.A:Math.Theor. , ,333 001.Grafke, T., and E. Vanden-Eijnden, 2019: Numerical computation of rare eventsvia large deviation theory. Chaos: An Interdisciplinary Journal of NonlinearScience ,
29 (6) , 063 118.Guyader, A., N. Hengartner, and E. Matzner-Løber, 2011: Simulation and Esti-mation of Extreme Quantiles and Extreme Probabilities.
Appl. Math. Optim. , , 171–196.Haidvogel, D., and I. Held, 1981: Homogeneous quasi-geostrophic turbulencedriven by a uniform temperature gradient. J. Atmos. Sci. , , 2644–2660.Hartmann, C., R. Banisch, M. Sarich, T. Badowski, and C. Sch ˜A tte,2013: Characterization of Rare Events in Molecular Dynam-ics. Entropy ,
16 (1)
Bernhard Haurwitz Memorial Lecture .Herbert, C., R. Caballero, and F. Bouchet, 2020: Atmospheric bistability andabrupt transitions to superrotation: Wave–jet resonance and hadley cell feed-backs.
Journal of the Atmospheric Sciences ,
77 (1) , 31–49.H.L.Kuo, 1949: Dynamics instability of two-dimensional nondivergent flow in abarotropic atmosphere.
J.Fluid.Mech. , , 105–122.Jougla, T., and D. G. Dritschel, 2017: On the energetics of a two-layer baroclinicflow. Journal of Fluid Mechanics , , 586–618.Kahn, H., and T. Harris, 1951: Estimation of particle transmission by randomsampling. Natl. Bur. Stand., Appl. Math. Ser. , , 27–30.Kaspi, Y., and G. R. Flierl, 2007: Formation of Jets by Baroclinic Instabilityon Gas Planet Atmospheres. J. Atmos. Sci. , , 3177–3194.Laurie, J., and F. Bouchet, 2015: Computation of rare transitions in thebarotropic quasi-geostrophic equations. New J. of Phys. , , 015 009.Lee, S., 1997: Maintenance of multiple jets in a barotropic flow. J. Atmos. Sci. , , 1726–1738.Lee, S., 2005: Baroclinic multiple zonal jets on the sphere. Journal of the at-mospheric sciences ,
62 (7) , 2484–2498.43emasquerier, D., B. Favier, and M. L. Bars, 2020: Zonal jets at thelaboratory scale: hysteresis and rossby waves resonance. arXiv preprintarXiv:2008.10304 .Lestang, T., F. Bouchet, and E. L´evˆeque, 2020: Numerical study of extrememechanical force exerted by a turbulent flow on a bluff body by direct andrare-event sampling techniques.
Journal of Fluid Mechanics , .Lorenz, E. N., 1967: The nature and theory of the general circulation of theatmosphere . World Meteorological Organization.Lucarini, V., and T. B´odai, 2017: Edge states in the climate system: exploringglobal instabilities and critical transitions.
Nonlinearity , , 32–66.Maragliano, L., G. Cottone, G. Ciccotti, and E. Vanden-Eijnden, 2010: Mappingthe Network of Pathways of CO Diffusion in Myoglobin. J.Am.Chem.Soc. , , 1010–1017.Marcus, P., 2004: Prediction of a global climate change on Jupiter. Nature , ,828–831.Metzner, P., C. Sch¨utte, and E. Vanden-Eijnden, 2009: Transition path theoryfor markov jump processes. Multiscale Modeling & Simulation , , 1192–1219.Onsager, L., and S. Machlup, 1953: Fluctuations and irreversible processes. Phys. Rev. , , 1505–1512.P. Del Moral, 2004: Feynman-Kac formulae, Genealogical and interacting parti-cle systems with applications . Probability and Applications, Springer-Verlag,New York.Paillard, D., 1998: The timing of pleistocene glaciations from a simple multiple-state climate model.
Nature ,
391 (6665) , 378–381.Panetta, R. L., 1993: Zonal Jets in Wide baroclinically unstable regions: per-sistence and scale selection.
J. Atmos. Sci. , , 2073–2106.Parker, J. B., and J. A. Krommes, 2013: Zonal flow as pattern formation. Physics of Plasmas ,
20 (10) , 100 703.Phillips, N., 1951: A simple three-dimensional model for the study of large-scaleextratropical flow patterns.
J. Meteor. , , 381–394.Pierrehumbert, R., D. Abbot, A. Voigt, and D. Koll, 2011: Climate of theneoproterozoic. Annual Review of Earth and Planetary Sciences , .Plotkin, D. A., R. J. Webber, M. E. O’Neill, J. Weare, and D. S. Abbot, 2019:Maximizing simulated tropical cyclone intensity with action minimization. Journal of Advances in Modeling Earth Systems ,
11 (4) , 863–891.44orco, C. C., and Coauthors, 2004: Cassini imaging science: Instrument char-acteristics and anticipated scientific investigations at saturn.
Space ScienceReviews ,
115 (1-4) , 363–497.Qiu, B., and W. Miao, 2000: Kuroshio path variations south of Japan: Bimodal-ity as a self-sustained internal oscillation.
J. Phys. Oceanogr. , , 2124–2137.Ragone, F., and F. Bouchet, 2019: Computation of extreme values of time aver-aged observables in climate models with large deviation techniques. Journalof Statistical Physics , 1–29.Ragone, F., J. Wouters, and F. Bouchet, 2018: Computation of extremeheat waves in climate models using a large deviation algorithm.
PROCEED-INGS OF THE NATIONAL ACADEMY OF SCIENCES OF THE UNITEDSTATES OF AMERICA ,
115 (1) , 24–29, doi: { } .Rahmstorf, S., 2002: Ocean circulation and climate during the past 120,000years. Nature , , 207–214.Ravelet, F., L. Mari´e, A. Chiffaudel, and F. Daviaud, 2004: Multistability andmemory effect in a highly turbulent flow: experimental evidence for a globalbifurcation. Phys.Rev.Lett. , , 164 501.Read, P., D. Kennedy, N. Lewis, H. Scolan, F. Tabataba-Vakili, Y. Wang,S. Wright, and R. Young, 2020a: Baroclinic and barotropic instabilities inplanetary atmospheres: energetics, equilibration and adjustment. NonlinearProcesses in Geophysics ,
27 (2) , 147–173.Read, P. L., R. M. Young, and D. Kennedy, 2020b: The turbulent dynamics ofjupiter’s and saturn’s weather layers: order out of chaos?
Geoscience Letters , , 1–18.Rogers, J. H., 1995: The Giant Planet Jupiter . Cambridge University Press, 418pp.Rolland, J., F. Bouchet, and E. Simonnet, 2016: Computing transition ratesfor the 1-D stochastic Ginzburg–Landau–Allen–Cahn equation for finite-amplitude noise with a rare event algorithm .
J. Stat. Phys. , , 277–311.Rolland, J., and E. Simonnet, 2015: Statistical behavior of adaptive multilevelalgorithm in simple models. J.Comput.Phys. , , 541–557.Rosenbluth, M., and A. Rosenbluth, 1955: Monte Carlo calculation of the av-erage extension of molecular chains. J.Chem.Phys. , , 356–359.S´anchez-Lavega, A., and Coauthors, 2008: Depth of a strong jovian jet from aplanetary-scale disturbance driven by storms. Nature , , 437–40.Schmeits, M. J., and H. A. Dijkstra, 2001: Bimodal behavior of the Kuroshioand the Gulf Stream. J. Phys. Oceanogr. , , 3435–3456.45chneider, T. M., B. Eckhardt, and J. A. Yorke, 2007: Turbulence transitionand the edge of chaos in pipe flow. Physical review letters ,
99 (3) , 034 502.Simon, A. A., M. H. Wong, and G. S. Orton, 2015: First results from the hubbleopal program: Jupiter in 2015.
Astrophys. J. , , 55–63.Simonnet, E., 2016: Combinatorial analysis of the adaptive last particle method. Stat. Comp. , , 211–230.Srinivasan, K., and W. R. Young, 2011: Zonostrophic Instability. J. Atmos.Sci. ,
69 (5) , 1633–1656.Tobias, S., and J. Marston, 2013: Direct statistical simulation of out-of-equilibrium jets.
Phys. Rev. Lett. ,
110 (10) , 104 502.Tziperman, E., and B. Farrell, 2009: Pliocene equatorial temperature: Lessonsfrom atmospheric superrotation.
Paleoceanography , , PA1101, doi:10.1175/1520-0469(2000)057 h i J. Phys. Oceanogr. , , 1346–1362.Vanden-Eijnden, E., and M. Heymann, 2008: The geometric minimum actionmethod: A least action principle on the space of curves. Comm. on Pure andApplied Math. , , 1052–1117.Webber, R. J., D. A. Plotkin, M. E. O’Neill, D. S. Abbot, and J. Weare, 2019:Practical rare event sampling for extreme mesoscale weather. Chaos: An In-terdisciplinary Journal of Nonlinear Science ,
29 (5) , 053 109, doi:10.1063/1.5081461, URL https://aip.scitation.org/doi/full/10.1063/1.5081461.Williams, G. P., 1979: Planetary ciruclations. 2: The Jovian quasi-geostrophicregime.
J. Atmos. Sci. , , 932–968.Williams, G. P., 2003: Jovian Dynamics. Part III: Multiple, Migrating, andEquatorial Jets. J. Atmos. Sci. , , 1270–1296.Willis, A. P., and R. R. Kerswell, 2009: Turbulent dynamics of pipe flow cap-tured in a reduced model: puff relaminarization and localized “edge” states. Journal of fluid mechanics , , 213–233.Woillez, E., and F. Bouchet, 2019: Barotropic theory for the velocity profile ofjupiter’s turbulent jets: an example for an exact turbulent closure. Journalof Fluid Mechanics , , 577–607.Woillez, E., and F. Bouchet, 2020: Instantons for the destabilization of the innersolar system. Physical Review Letters ,
125 (2) , 021 101.Wouters, J., and F. Bouchet, 2016: Rare event computation in deterministicchaotic systems using genealogical particle analysis.
Journal of Physics A:Mathematical and Theoretical ,
49 (37)49 (37)