## Abstract

Strong eastward jets at the equator have been observed in many planetary atmospheres and simulated in numerical models of varying complexity. However, the nature of the transition from a conventional state of the general circulation, with easterlies or weak westerlies in the tropics, to such a superrotating state remains unclear. Is it abrupt or continuous? This question may have far-reaching consequences, as it may provide a mechanism for abrupt climate change in a planetary atmosphere, both through the loss of stability of the conventional circulation and through potential noise-induced transitions in the bistability range. We study two previously suggested feedbacks that may lead to bistability between a conventional and a superrotating state: the Hadley cell feedback and a wave–jet resonance feedback. We delineate the regime of applicability of these two mechanisms in a simple model of zonal acceleration budget at the equator. Then we show using numerical simulations of the axisymmetric primitive equations that the wave–jet resonance feedback indeed leads to robust bistability, while the bistability governed by the Hadley cell feedback, although observed in our numerical simulations, is much more fragile in a multilevel model.

## 1. Introduction

A long-standing question in the study of the general circulation of the atmosphere, formulated early on by Lorenz (1967), is the uniqueness of the solution, for fixed boundary conditions. This question is an important one, because it may have deep consequences on climate dynamics. Indeed, in the presence of multiple attractors, the system may exhibit abrupt transitions from one to the other, induced either by internal variability or by an external forcing. Paleoclimatic records provide evidence for such abrupt climate changes (e.g., Dansgaard–Oeschger events; Dansgaard et al. 1993). Events of this type have so far been linked to nonlinear behavior of the oceanic circulation. For instance, it is well understood from conceptual models how the Atlantic meridional overturning circulation can be bistable (Dijkstra and Ghil 2005), and some full complexity, high-resolution ocean models show evidence of bistability (Jackson and Wood 2018). Some feedback mechanisms rely solely on internal ocean dynamics, while others invoke a coupling with ice sheets and sea ice [see Boers et al. (2018) for a recent example]. The atmosphere itself may admit multiple equilibria. As a matter of fact, turbulent flows often exhibit coexisting steady states for given external parameters, as well as spontaneous transitions between the two stable states, as has been reported in both numerical studies (Bouchet and Simonnet 2009; Cortet et al. 2010; Bouchet et al. 2019) and laboratory experiments (Berhanu et al. 2007; Cortet et al. 2010; Saint-Michel et al. 2013; Michel et al. 2016). Some of these experiments (Weeks et al. 1997; Tian et al. 2001) are actually inspired by geophysical flows (Charney and DeVore 1979). However, the question remains if such phenomena could occur at the level of the general circulation of the atmosphere.

An interesting candidate for bistability of the general circulation of the atmosphere is *superrotation* (Held 1999): this refers to an atmospheric flow for which there exists a region carrying a larger angular momentum than the one associated to solid-body rotation at the equator. While the *conventional* circulation of the atmosphere of Earth has midlatitude westerly jets and weak easterlies in the tropics (Lee 1999; Dima et al. 2005) (and everywhere smaller angular momentum than the surface at the equator), a superrotating atmosphere exhibits westerlies in the tropics. This is actually observed on other planets of the solar system, such as Jupiter, Saturn (and its moon Titan), or Venus (see, e.g., Read and Lebonnois 2018). On Earth, superrotation may have played a role in the climate of the past: it was observed in numerical simulations of warm climates such as the Eocene (Caballero and Huber 2010), and it has been suggested that it could explain the permanent El Niño conditions indicated by paleoclimatic proxies during the Pliocene (Tziperman and Farrell 2009). Another indicator of the robustness of superrotation is that it has been observed in numerical experiments with models of varying complexity: shallow-water models (Scott and Polvani 2008; Showman and Polvani 2010, 2011; Suhas et al. 2017), two-level primitive equations (Suarez and Duffy 1992; Saravanan 1993), and multilevel comprehensive GCMs (Kraucunas and Hartmann 2005; Schneider and Liu 2009; Caballero and Huber 2010; Showman and Polvani 2011; Arnold et al. 2012; Potter et al. 2014).

A natural question to ask first is, how is superrotation maintained at a dynamical level? Because axisymmetric dynamics, in the absence of forcing and dissipation, conserve angular momentum, the mean meridional circulation cannot generate superrotation. Momentum diffusion opposes superrotation, since upgradient fluxes of angular momentum are required. Hence, superrotation can only be achieved by eddy fluxes. This is often referred to as Hide’s theorem (Hide 1969). There are potentially many ways eddies could accelerate the zonal wind toward the east in the tropics, and several routes to superrotation have already been found. In a first series of studies, the basic physical parameters of the planet, such as the planetary rotation rate (Dias Pinto and Mitchell 2014) or the radius of the planet (Mitchell and Vallis 2010; Potter et al. 2014) are modified. The emerging scenario in this type of setup is that a hydrodynamic instability known as the *Kelvin–Rossby instability* (Iga and Matsuda 2005; Wang and Mitchell 2014; Zurita-Gotor and Held 2018) generates the eddies that converge momentum in the tropics. Instead of relying on an instability, a second thread of works has explored the possibility of stimulating wave emission from the tropics to account for equatorial momentum convergence, akin to the classical picture for midlatitude jets (Vallis 2006). Enhanced wave activity in the tropics can be the result of several physical processes: convection, day–night contrast in tidally locked exoplanets (Merlis and Schneider 2010; Showman and Polvani 2011), etc. Broadly speaking, such processes can be modeled as nonzonal heating of the tropics: idealized GCM studies including such an additional forcing term have led to abrupt transitions to superrotation once the forcing amplitude reaches a certain threshold (Suarez and Duffy 1992; Saravanan 1993; Kraucunas and Hartmann 2005; Arnold et al. 2012).

In fact, coexistence of the superrotating state with the conventional circulation for some range of parameters requires more than just eddy momentum flux convergence onto the equator. Indeed, some *positive feedback* mechanism is needed, so that the zonal-mean zonal wind budget may admit several solutions. Such a feedback mechanism may come directly from the eddy forcing or, alternatively, from the mean meridional circulation. The first possibility has been explored in particular by Arnold et al. (2012), who suggested a resonant feedback mechanism based on the properties of equatorial Rossby waves on a background mean flow. Relying on an explicit computation of the linear response of a shallow-water atmosphere to nonzonal tropical heating, in the spirit of the pioneering work of Matsuno (1966) and Gill (1980), they have argued that the amplitude of the response depends on the background zonal wind in such a way that a resonance appears close to the opposite of the phase velocity of free Rossby waves.

The second possibility was suggested by Shell and Held (2004, hereafter SH04), who showed that the Hadley cell itself could admit multiple equilibrium states. Indeed, a conventional Hadley cell with updraft on the equator advects low-momentum air into the upper troposphere, thereby inhibiting the onset of superrotation. However, the contribution to the zonal momentum budget is the product of two terms: *ω*∂_{p}*u*. The core idea behind the feedback structure of the Hadley cell is that when westerly winds increase in the tropical upper troposphere, vertical shear increases, but vertical velocity decreases. As a result of this nonlinear behavior, there exists a first regime, for weak westerly wind, where the feedback is negative, and a second regime where the feedback is positive. For even larger westerly wind, the feedback becomes negative again. Assuming that the Hadley cell and some frictional dissipation balance positive eddy momentum flux convergence at the equator, this feedback structure leads to multiple equilibria (SH04). These arguments are supported by numerical simulations in a simple framework (1D axisymmetric shallow-water equations with a constant imposed torque). A natural question to ask is whether this behavior remains in more realistic conditions.

In this paper, we explore the robustness of these two bistability mechanisms: Hadley cell feedback and resonant response to equatorial heating. First, we explicitly show in an analytical model how the resonant structure of the eddy momentum flux convergence can lead to bistability, and observe the corresponding hysteresis phenomenon in numerical simulations of the axisymmetric primitive equations. Second, we investigate whether the results of SH04 extend to a multilevel model. We find numerically that bistability may be obtained in this framework, but that it is relatively fragile as it depends sensitively on vertical viscosity. Finally, we investigate the interplay between the two mechanisms. We show that depending on the parameter characterizing the width of the wave–jet resonance, two types of superrotating states can be found. For wide resonances, the superrotating state has a weaker mean meridional circulation than the conventional state, and the range of forcing amplitudes for which both states coexist is quite small (Hadley cell–driven superrotation). On the other hand, for narrow resonances, the strength of the mean meridional circulation does not change much across the bifurcation point, and the coexistence range is much wider (resonance-driven superrotation).

To reach these conclusions, we combine theoretical arguments obtained in a simplified framework based on the shallow-water model (section 2), and numerical simulations of the axisymmetric primitive equations (section 3). More precisely, the theoretical part relies on the fact that a linear response computation of the Matsuno–Gill type, has been found to agree relatively well with GCM results (Arnold et al. 2012). In section 2c, after recalling this computation, we describe the wave–jet resonance mechanism that provides a positive feedback. Then we present analytical arguments to disentangle the effects of the two nonlinear mechanisms by studying the fixed points of the zonal momentum budget at the equator (sections 2d and 2e). Finally, we test the scenarios outlined through the analytical study of the shallow-water model in a more realistic model of the atmosphere, by carrying out numerical simulations of the 2D axisymmetric primitive equations (section 3).

## 2. Bistability in an analytical model of equatorial momentum balance

### a. The shallow-water model

We first consider the simplest possible model that can account for both feedback mechanisms: a thin layer of fluid, described by the shallow-water equations, exchanging mass and momentum with a quiescent underlying layer. The fluid is forced by diabatic heating *Q* and dissipates energy through a Rayleigh friction *ε*. In spherical coordinates, these equations may be written as

where *u* and *υ* are the zonal and meridional components of the wind, *h* the thickness of the fluid layer, *ϕ* the latitude, *λ* the longitude, $g\u22c6$ the reduced gravity, Ω the rotation rate, and *a* the planetary radius. The mass source/sink term *Q* accounts both for radiative forcing and for an additional nonzonal heating term, confined in the tropics, which represents in a rough manner convective effects or day–night contrast in tidally locked exoplanets. As a consequence, the fluid layer also exchanges momentum with the underlying “sponge” layer through the terms *R*_{u} = −*Qu*/*h*Θ(*Q*) and *R*_{υ} = −*Qυ*/*h*Θ(*Q*), where Θ is the Heaviside function. This mechanism provides a rudimentary representation of the Hadley cell in the shallow-water model. It is required to obtain superrotation in this setup (Showman and Polvani 2010).

We now decompose all the fields into their zonal average, denoted by an overbar, and their eddy component, denoted by a prime: $u=u\xaf+u\u2032$, $\upsilon =\upsilon \xaf+\upsilon \u2032$, $h=h\xaf+h\u2032$. In this context, the zonal mean wind profile $u\xaf\u2061(\varphi )$ satisfies the equation:

Our goal is to study the possibility of multiple equilibria in this equation. In general, this depends on the form of the eddy momentum flux convergence $F=\u2212\u2061(1/a\u2009cos\varphi )\upsilon \u2032\u2202\varphi \u2061(u\u2032\u2009cos\varphi )\xaf$. In a first step, we assume that it does not depend on $u\xaf$ and discuss the other feedback mechanism $R\xafu$, associated to the Hadley cell (section 2b). We shall discuss the wave–mean flow interaction in section 2c.

### b. Simplified zonal momentum balance at the equator

At the equator, in a perpetual equinox configuration, the steady-state zonal-mean zonal momentum budget, (4), reduces to a balance between the eddy forcing *F*, the vertical advection by the Hadley cell $R\xafu$ and the frictional dissipation. This balance writes $F+R\xafu\u2212\epsilon u\xaf=0$. In this section, we study the existence of multiple solutions to this balance equation. Reducing this way the problem to a zero-dimensional model allows for a qualitative understanding of the physical mechanisms, as we can easily draw the different terms of this balance relation as functions of the parameters of the problem. As a matter of fact, SH04 have used this simple zonal momentum balance model to show how, for a constant forcing *F*, the Hadley cell feedback leads to bistability. We recall their argument in this section. As we shall only be working with zonally averaged fields, we should drop the overbar from now on. The value of all fields at the equator will be denoted with a null subscript.

Modeling radiative forcing by a Newtonian relaxation *Q* = (*h*_{eq} − *h*)/*τ* to a prescribed radiative equilibrium profile *h*_{eq} with relaxation time *τ*, and recalling that *R* = −*Qu*/*h*Θ(*Q*), the balance relation becomes

A relation between the layer thickness and the zonal wind velocity can be obtained through a simple model of the Hadley cell (Held and Hou 1980; Vallis 2006; SH04). The idea is that the thickness *h* is in geostrophic equilibrium with the angular momentum–conserving wind *u*_{m} = (*u*_{0} + Ω*a* sin^{2}*ϕ*)/cos*ϕ* in the tropics, $(g\u22c6/a)\u2202\varphi h=\u22122\Omega um\u2009sin\varphi $, and in radiative equilibrium, *h* = *h*_{eq} outside. Integrating the geostrophic equilibrium equation and matching the resulting profile with the radiative equilibrium at a latitude determined by mass conservation yields

Note that the model only makes sense for *u*_{0} < *u*_{0eq}. Introducing the nondimensional variables *U* and *H* through *h*_{0} = *Hh*_{0eq}, *u*_{0} = *Uu*_{0eq}, the Eqs. (6) and (5) reduce to the simple algebraic system

where we have assumed *H* ≈ 1, with

Parameter values are given in Table 1. Hence, the balance between the forcing, Rayleigh friction, and the Hadley cell advecting low-momentum wind from the lower layer is governed by the equation

This theory makes the feedback structure of the Hadley cell very clear: the *pU*(*U* − 1)^{2} term acts as a positive feedback between the two roots of its derivative, 1/3 < *U* < 1, and as a negative feedback for *U* < 1/3 and for *U* > 1. When the forcing term is a constant imposed torque, the equation is a simple cubic equation, and the condition for bistability can be easily obtained. A necessary condition is 0 ≤ *r*/*p* ≤ 1/3: it is the condition for the function *pU*(*U* − 1)^{2} + *rU* to have a local maximum. With the default parameter values, *r*/*p* ≈ 0.1, and the above condition is fulfilled. An illustration is provided in Fig. 1: we plot separately *U*(*U* − 1)^{2} + *rU*/*p* (the sum of vertical advection by the Hadley cell and friction) and the constant forcing *q*/*p* for two values of the ratio *r*/*p*. When this ratio is small enough (0 ≤ *r*/*p* ≤ 1/3), the positive feedback of the Hadley cell leads to the existence of three solutions to (10) for some range of forcing amplitude *q*/*p* (indicated by the two dashed lines in Fig. 1, right), two stable ones (*U* ≈ 0.1 and *U* ≈ 1.2 in the figure) and an unstable one (*U* ≈ 0.6 in the figure). As the forcing amplitude sweeps the range of positive values, two (saddle node) bifurcations are encountered: we start from an equilibrium with weak equatorial wind (*U* ≈ 0) for low values of the forcing, which loses stability when the forcing amplitude increases past some value (the dashed line at *q*/*p* ≈ 0.16 in Fig. 1, right). The system then jumps abruptly to the equilibrium state with strong westerly wind (*U* ≈ 1.3) and remains on this branch if the forcing is further increased. Now, this superrotating equilibrium in turn loses stability when the forcing decreases below some value (the dashed line at *q*/*p* ≈ 0.02 in Fig. 1, right). We have just described a *hysteresis* phenomenon. When the ratio *r*/*p* becomes too large, the negative feedback of friction overcomes the positive feedback of the Hadley cell, and there is only one solution to (10) (*U* ≈ 0.05 in the figure) for the whole range of forcing amplitude *q*/*p* (Fig. 1, left).

In fact, the eddy forcing *F* should not be a constant. In the next section, we show that it may be modeled as a resonant function of *U*, and we discuss in sections 2d and 2e the consequences for the balance relation (10).

### c. The wave–jet resonance: Matsuno–Gill computation of the eddy momentum flux convergence

The goal here is to compute the eddy momentum flux convergence induced by a nonzonal tropical heating. We rely on a classical approach, pioneered by Matsuno (1966) and Gill (1980): we assume that the zonal-mean zonal wind evolves slowly compared to the eddies, and we compute the linear eddy response to the heating term with a constant background wind. A major advantage is that for the shallow-water equations on an equatorial beta plane, the linear response can be computed explicitly. Typically, the stationary response of the atmosphere to a localized heating consists in the superposition of an equatorially trapped Kelvin wave east of the source and a Rossby wave west of the source. The relative phases of the two standing waves depend on the parameters of the problem. In a wide range of parameter values, the Matsuno–Gill response converges westerly momentum at the equator (Showman and Polvani 2010, 2011); Arnold et al. (2012) further argued that the response exhibits a resonant structure. Here, after briefly recalling their result, we compute the associated eddy momentum flux convergence and discuss its resonant structure.

To make the problem analytically tractable, we rewrite Eqs. (1)–(3) linearized around a uniform zonal mean flow $u\xaf$, using the beta-plane approximation (Vallis 2006):

where *x* and *y* represent the zonal and meridional directions, respectively, and *β* = 2Ω/*a* is the beta effect at the equator (on Earth, *β* ≈ 2.289 × 10^{−11} m^{−1} s^{−1}). We have assumed that the background flow has no meridional component $(\upsilon \xaf=0)$ and no meridional shear $(\u2202yu\xaf=0,\u2202yh\xaf=0)$. We are also neglecting the momentum exchange with the underlying layer. In the rest of this section, we use as time and length units $T=1/2\beta cg$ and $L=cg/\u2061(2\beta )$, with $cg=g\u22c6h\xaf$ the velocity of pure gravity waves. For simplicity, we also absorb the $g\u22c6$ factor into *h* (so in this section *h* is actually a nondimensionalized geopotential) and *Q*.

In the absence of mean flow $(u\xaf=0)$, Matsuno (1966) found the normal modes of the linear system (11)–(13) without forcing and dissipation (*Q* = *ε* = 0), and computed the stationary solution to the forced-dissipative problem by projecting onto those normal modes. We refer to Vallis (2006, chapter 8) or Gill (1982, chapter 11) for details of the methods, including the dispersion relation and spatial structure of the modes. The uniform mean flow $u\xaf$ Doppler shifts the response without modifying the structure of the modes. For a stationary tropical heating of the form $Q=Q0\u2009cos\u2061(kx)e\u2212y2/4$, Arnold et al. (2012) computed the stationary response and separated it into two contributions, both with zonal wavenumber *k*: a Kelvin mode $(uK\u2032,\upsilon K\u2032,hK\u2032)$, with

and a Rossby mode $(uR\u2032,\upsilon R\u2032,hR\u2032)$ with meridional wavenumber *n* = 1:

with $\gamma X=\epsilon /k\u2061(u\xaf+cX)$ a nondimensional parameter defined for the two indices *X* = *K* and *X* = *R*, and *c*_{X} the phase velocity of the free waves: *c*_{R} = −1/(3 + 2*k*^{2}) and *c*_{K} = 1 in nondimensional units. The total response is given by $u\u2032=uR\u2032+uK\u2032$, $\upsilon \u2032=\upsilon R\u2032,$, and $h\u2032=hR\u2032+hK\u2032$.

From this point, an explicit formula can be obtained for the corresponding eddy momentum flux convergence:

where the first and second terms in the braces correspond, respectively, to the contribution from the Rossby mode only $(\u2212\u2202y\u2329uR\u2032\upsilon R\u2032\u232a)$ and to the interaction between the Kelvin and Rossby modes $(\u2212\u2202y\u2329uK\u2032\upsilon R\u2032\u232a)$. The spatial structure of the eddy momentum flux convergence $F\u2061(u\xaf,y)$ as a function of the background mean-flow velocity $u\xaf$, and its contribution from the Rossby mode only $(\u2212\u2202y\u2329uR\u2032\upsilon R\u2032\u232a)$, are shown in Fig. 2 in dimensional units. It is obtained using parameter values *c*_{g} = 49 m s^{−1}, *ε* = 1 day^{−1} and *ka* = 1. Going back to the dimensional expression for the phase velocity of the Rossby and Kelvin waves yields the corresponding numerical values:

With these parameters, the Rossby deformation radius is *L* ≈ 1000 km. As expected, the eddy momentum flux convergence is symmetric with respect to the equator. For all the values of the background mean flow $u\xaf$, the Rossby component (Fig. 2, right) is positive in the equatorial region (within one deformation radius of the equator, roughly speaking, i.e., about 10°), inducing eastward acceleration of the jet, then negative (between one and two deformation radii) and positive again in the extratropics. A similar spatial structure is found in the full eddy momentum convergence flux (Fig. 2, left), except when the background mean flow corresponds to strong easterly wind. In that case, the contour lines are distorted, up to a point where the eddy momentum flux convergence becomes negative in the equatorial region. Both the full eddy momentum flux convergence and its Rossby component exhibit local maxima and minima, corresponding to resonance and antiresonance structures. Let us further describe these mechanisms by focusing on the equatorial area.

Let us denote $FRK\u2061(u\xaf)$ the full eddy momentum flux convergence at the equator (*y* = 0) and $FR\u2061(u\xaf)$ the contribution from the Rossby mode:

It is easily seen from (21) that the eddy momentum flux convergence at the equator $FRK\u2061(u\xaf)$ is positive as long as $u\xaf>\u2061(3cR\u2212cK)/2$. On the other hand, $FR\u2061(u\xaf)$ is always positive. Besides, $FR\u2061(u\xaf)$ has the shape of a Lorentz curve. The curves $FRK\u2061(u\xaf)$ and $FR\u2061(u\xaf)$ are shown in Fig. 3, using the same parameter values as above. Both cases exhibit a resonance for background velocities $u\xaf\u2248\u2212cR$. When the Kelvin mode is taken into account, there is a secondary peak with opposite sign for $u\xaf\u2248\u2212cK$. For the existence of multiple steady states, a critical point is the sign of the feedback associated to the eddy momentum flux convergence, that is, the sign of the derivative with respect to $u\xaf$, $dFRK\u2061(u\xaf)/du\xaf$ or $dFR\u2061(u\xaf)/du\xaf$. From Fig. 3, it is clear that the feedback is positive below the resonance ($dFRK\u2061(u\xaf)/du\xaf>0$ for $\u2212cK<u\xaf<\u2212cR$) and negative above it ($dFRK\u2061(u\xaf)/du\xaf<0$ for $u\xaf>\u2212cR$). Ultimately, the existence of multiple steady states for the mean flow $u\xaf$ depends on the other acceleration terms: qualitatively, bistability with a superrotating steady state hinges on the positive feedback described above overcoming the negative feedbacks due to other effects, such as linear friction for instance (see section 2d).

Of course, it seems natural that the linear response framework should break down when the amplitude of the forcing becomes too large. Then the dynamical feedback of the eddies on the mean flow cannot be neglected anymore. The linear and nonlinear responses have been compared for instance analytically using perturbative expansion (Gill and Phlips 1986), or numerically using idealized models (Nobre 1983) and full GCM simulation (Lutsko 2018). In practice however, it has been found in many studies that the linear response computation provides a useful starting point for interpreting results from observations or full nonlinear GCMs (Moura and Shukla 1981; Gill and Rasmusson 1983; Neelin 1988; Jin and Hoskins 1995; Kraucunas and Hartmann 2005; Norton 2006; Sobel and Maloney 2012; Arnold et al. 2012). Here, it should be kept in mind that the spatial structure of the response may differ significantly from the linear response in the superrotating state (Lutsko 2018). However, most of our reasoning does not depend on the details of the spatial structure, but rather on the resonant behavior that has been reported to hold in a full nonlinear GCM (Arnold et al. 2012) for heating rates and spatial structure similar to those considered here. Hence, we shall consider that the eddy momentum flux convergence computed in this section is a reasonable working hypothesis, and we shall now study how it may lead to bistability.

### d. Qualitative behavior of the wave–jet resonance

As shown in Fig. 3, the eddy momentum flux convergence associated to the full response (i.e., including the projection on the Kelvin mode) is amplified compared to the Rossby mode response, but the overall structure remains qualitatively similar (if we except the negative tail for strong easterly background winds). The functional form of the Rossby wave forcing $FR\u2061(u\xaf)$ (it is a Lorentzian function) makes it simpler than the full resonant eddy forcing $FRK\u2061(u\xaf)$, and it also reduces the number of free parameters. In this section, we exploit this to obtain a qualitative understanding of the steady states of the momentum budget, (10).

Injecting (22), in dimensional units, into the normalized parameters, (9), we obtain the corresponding forcing term for the zonal momentum balance model:

with $Q\u02dc=\beta Q02\tau 2/\u2061(6ru0eq)$ and Λ = (*ku*_{0eq}/*ε*)^{2}, where *k* is the zonal wavenumber of the forcing and *ε* the friction coefficient. In addition to the parameter *r*/*p* discussed in section 2b, which governs the competition between the feedbacks of the two damping mechanisms, vertical advection by the Hadley cell and friction, there are two parameters characterizing the eddy forcing. First, the position of the resonance is governed by a purely dynamical quantity −*c*_{R}/*u*_{0eq}, the phase velocity of free Rossby waves, nondimensionalized by the velocity associated to the radiative forcing. Second, the width of the resonance peak is governed by the parameter Λ, which depends upon the wavenumber of the nonzonal forcing, but also the radiative forcing and friction. Together, these parameters select the range of background wind values that can be maintained by the eddy forcing.

Ideally, we would like to know when, in the 3D parameter space (*r*/*p*, Λ, *c*_{R}/*u*_{0eq}), (10) admits multiple solutions for some range of forcing amplitude $Q\u02dc$. Even within this simplified framework, it is difficult to obtain such a full classification (in general, it amounts to counting the real roots of a fifth-order polynomial), and we shall not attempt to do so. Instead, let us try to get some qualitative insight by discussing some cases of physical relevance. Let us first discriminate the possibilities based on the parameter *r*/*p*.

When *r*/*p* > 1/3, the negative feedback of friction overcomes the positive feedback of the Hadley cell. The sum of the two is a monotonously increasing function of *U*. For bistability to appear, we need the wave–jet resonance to be sufficiently peaked for the positive feedback due to the eddy forcing to overcome the negative feedback of friction close to the resonance peak. This requires that the region with a significant positive feedback (i.e., the bump of the Lorentzian) is entirely contained in the *U* > 0 range, which can be expressed as Λ ≫ (*u*_{0eq}/*c*_{R})^{2}. This can be checked explicitly by setting *R* = 0 in the simplified zonal momentum balance, which yields the equation *q*_{R}(*U*) = *rU*: this is a cubic equation, which can be solved exactly. In this case, a bistability range appears as soon as Λ > 3. Then, provided the forcing amplitude is large enough, there are three solutions to the balance equation: an unstable one and two stable ones. We refer to this case as *resonance-driven* bistability: it is illustrated in Fig. 4 (top left). One of the stable states corresponds to *U* ≈ 0—on the left flank of the resonance peak—and the other one is a superrotating state, with *u*_{0} ≈ −*c*_{R} (for an infinitely narrow resonance)—on the right flank of the resonance peak. A first saddle-node bifurcation occurs when $Q\u02dc$ increases and the resonance peak intersects the friction curve, corresponding to the appearance of the superrotating state. A second saddle-node bifurcation occurs when the forcing becomes significantly nonzero for *U* close to zero, corresponding to the loss of stability of the conventional circulation. However, this second bifurcation is expected to occur for very large forcing amplitudes: in other words, the range of forcing amplitude for which bistability occurs should be very wide in this scenario. Note that the stable superrotating state is very close to the unstable state.

When *r*/*p* < 1/3, the positive feedback of the Hadley cell acts in the region 1/3 < *U* < 1. Multiple steady states may also exist in this case. First, for an infinitely wide resonance (Λ ≪ 1), we should recover the case studied in section 2b, governed entirely by the Hadley cell feedback. Second, for a very narrow resonance (Λ ≫ 1), bistability should also be obtained similarly to the case *r*/*p* > 1/3 discussed in the previous paragraph (in this case, it might even be possible to obtain three coexisting stable states). Now, let us discuss the case of a resonant eddy forcing with finite width (for the figures, we choose Λ = 10). We distinguish three cases, based on the position of the resonance, for a fixed value of *r*/*p*:

When −

*c*_{R}/*u*_{0eq}< 1/3 (Fig. 4, top right), the same kind of scenario as in the previous paragraph unfolds, except that in the regime where three equilibria exist, they are all on the right flank of the resonance peak, that is, in the region where the eddy forcing feedback is negative. Hence, bistability relies on the Hadley cell feedback, like in section 2b.When 1/3 < −

*c*_{R}/*u*_{0eq}< 1 (Fig. 4, bottom left), bistability is again possible. This time, the two stable states are always on different flanks of the resonance peaks, while the unstable one moves from the right flank to the left flank as the forcing amplitude increases (until it annihilates with the low-wind stable state at the saddle-node bifurcation). In other words, the appearance of the superrotating state occurs because the positive feedback of the Hadley cell sets in, like in the previous case, but, on the other hand, the loss of stability of the conventional state is due to the positive wave–jet feedback prevailing over the negative feedback of the Hadley cell.When −

*c*_{R}/*u*_{0eq}> 1, (Fig. 4, bottom right), the first saddle-node bifurcation, corresponding to the appearance of the superrotating case, occurs on the left flank of the resonance peak. As the forcing amplitude keeps increasing, the superrotating state moves to the right flank of the Lorentzian. In this case both feedbacks contribute with the same sign.

### e. Quantitative discussion

In section 2d, we have considered independently the role of the three nondimensional parameters (*r*/*p*, Λ, and *c*_{R}/*u*_{0eq}) characterizing the balance between zonal acceleration due to resonant eddy forcing, vertical advection by the Hadley cell and friction. We have given simple criteria for bistability due to the Hadley cell feedback (*r*/*p* ≤ 1/3) and the wave–jet resonance [Λ ≫ (*u*_{0eq}/*c*_{R})^{2}; i.e., $k2cR2/\epsilon 2\u226b1$]. Let us now discuss the applicability of these regimes for some typical parameter values.

We first consider the parameter values from SH04, summarized in Table 1, supplemented with forcing parameters *ka* = 1 and *c*_{R} = −16 m s^{−1}. Such values fall under the scenario where there is bistability, governed by the wave–jet feedback because, although *r*/*p* < 1/3, the resonance is very strongly peaked ($k2cR2/\epsilon 2\u22486\xd7104$, Λ ≈ 10^{6}), like in the top-left panel of Fig. 4. However, the value used for friction is lower than typical values for the atmosphere of Earth, by several orders of magnitude (about 0.001 day^{−1}, instead of 0.1–1 day^{−1}; e.g., Held and Suarez 1994). As explained by SH04, this is essentially a consequence of the simplistic vertical structure of the model. In reality, dissipative processes modeled by linear friction have a more complex physical nature (eddy viscosity, wave breaking, etc.). Increasing *ε* and keeping the other parameters constant, one may easily transition to a case without bistability (because *r* also increases and the resonance becomes too broad) or a case where bistability is governed by the Hadley cell if we keep *r* constant by decreasing simultaneously the radiative cooling time *τ* (one could equivalently decrease the layer thickness at the equator *h*_{0eq}). We list in Table 2 estimates of parameter values for different planetary atmospheres, which indicate that the bistability regime governed by the wave–jet feedback seems relevant in most cases of interest, although perhaps marginally for Earthlike planets. However, this conclusion hinges crucially on the friction coefficient *ε*, which is difficult to estimate, as mentioned above. Investigations with a more realistic model would be necessary to better understand which physical parameters govern the resonance. In section 3, we adopt a refined description of the vertical structure of the atmosphere, replacing linear friction by a turbulent diffusion scheme, but prescribing the resonance width. Before doing so, let us comment on the differences between the two bistability regimes in the simple model.

The hysteresis curves obtained by tracking the solution of the balance equation, (10), as we ramp up and down the forcing amplitude, both for velocity *U* and vertical advection of zonal momentum *R*, are shown in Fig. 5. We show the same figure for two cases: one where bistability is governed by the Hadley cell feedback (*ε* = 0.01 *τ*^{−1} = 1 day^{−1}; Fig. 5, left), and one where bistability is governed by the resonant eddy forcing feedback (*ε* = 0.01 *τ*^{−1} = 0.1 day^{−1}, Fig. 5, right). As anticipated in the qualitative study, while both cases exhibit bistability, the bistability range is much wider in the case dominated by the resonant eddy forcing. Care should be taken with the terminology: eddy forcing with a narrow resonance (i.e., acting on a narrow range of *U*) corresponds to a wide bistability range (coexistence of two steady states on a wide range of $Q\u02dc$), and vice versa. The behavior of the Hadley cell is also quite different in the two cases: it collapses in the superrotating state governed by the Hadley cell feedback (*R* decreases sharply on the lower branch of the hysteresis cycle), but this is not necessarily the case in the superrotating case induced by the resonant eddy forcing (*R* remains larger than in the conventional circulation over a wide range of forcing amplitudes).

## 3. Bistability in the axisymmetric primitive equations

### a. Numerical setup

We now investigate the interplay between the resonant eddy forcing and the Hadley cell feedbacks in a more realistic context. Instead of the zonally averaged shallow-water equations [(4) for the zonal wind], we consider the axisymmetric primitive equations:

where the zonally averaged zonal and meridional wind *u* and *υ* are now 2D fields (depending on latitude *ϕ* and pressure *p*); *ω* = *Dp*/*Dt* is the zonally averaged vertical velocity in pressure coordinates; and Φ, *T*, and *θ* are the zonally averaged geopotential, temperature, and potential temperature. Dissipative effects are represented generically by the zonal and meridional components of the zonal-mean stress tensor *τ*_{u} and *τ*_{υ}, respectively; and *F*_{u} and *F*_{υ} represent the divergence of the Reynolds stress tensor (i.e., the eddy forcing). In our numerical simulations, we prescribe the eddy forcing to account for the wave–jet resonance in a simplified manner. The only diabatic heating term is a Newtonian relaxation term that drives the temperature field toward a prescribed radiative–convective equilibrium: $\theta e\u2061(p,\varphi )=max\u2061(200\u2061(p0/p)R/cp,\theta \u22c6\u2212\Delta h\u2009sin2\varphi \u2212\Delta \upsilon \u2009ln\u2061(p0/p)\u2009cos2\varphi )$. We use standard values for the coefficients (Held and Suarez 1994): $\theta \u22c6=315\u2009K$, Δ_{h} = 60 K, Δ_{υ} = 10 K. The relaxation time *τ* is as in Held and Suarez (1994).

The main difference with the shallow-water model considered in section 2 is a more accurate description of the vertical structure of the atmosphere, which allows to properly resolve vertical momentum transport by the Hadley cell, through the term *ω*∂_{p}*u*.

The model is solved numerically using the Climt framework (Caballero et al. 2008; Monteiro et al. 2018), which solves the above equations in flux form, using a simple upwind scheme (Smolarkiewicz 1983). We use 91 grid points in latitude (i.e., a resolution slightly smaller than 2°), and 45 vertical levels (see appendix). The initial condition is a state of rest (*u* = *υ* = *ω* = 0) with a constant temperature field *T* = 283.15 K. A turbulent diffusion scheme is used for the stress tensor *τ*; we shall denote *ν* the kinematic viscosity (in m^{2} s^{−1}) in the vertical direction. Our runs use the value *ν* = 0.5 m^{2} s^{−1} by default. Surface momentum drag is parameterized through a bulk aerodynamic formula (Caballero et al. 2008), akin to the linear friction considered above.

We carry out two series of numerical experiments, corresponding to two different kinds of prescribed eddy forcing:

A resonant eddy forcing

*F*_{u}=*F*_{RK}(*u*(*ϕ*= 0)) with spatial structure given by the Matsuno–Gill problem (section 2c) and with a varying amplitude given by (21) (see Fig. 2, left). These experiments are designed to reproduce the behavior observed in GCM studies with nonzonal tropical heating, such as Suarez and Duffy (1992), Saravanan (1993), Kraucunas and Hartmann (2005), and Arnold et al. (2012). Since the model considered here is axisymmetric, we need to parameterize the effect of the eddy forcing, for which we use the analytical computation of the linear response of a shallow-water atmosphere to a nonzonal tropical heating carried out in section 2c. This allows us to explore parameter space at a much lower computational cost.A constant eddy forcing

*F*_{u}with the same spatial structure as above and fixed amplitude*F*_{RK}(*U*= 0). These experiments amount to adding a vertical dimension to the setup of SH04 (the meridional structure is also slightly different).

In both cases, the vertical structure is arbitrarily chosen as a Gaussian profile $e\u2212\u2061(p\u2212p0)/2\sigma 2$ centered on the *p*_{0} = 300-hPa level, as in Caballero and Carlson (2018), with vertical extent *σ* = 200. While the details of the vertical structure of the forcing affect the vertical profile of the flow, they do not change the results presented here, as long as the forcing acts in the upper tropical troposphere. In these experiments, we always have *F*_{υ} = 0. Note also that we do not parameterize eddy heat transport. Below, we vary the forcing amplitude: this refers to the coefficient *Q*_{0} entering (19), nondimensionalized in the way explained in section 2. While this coefficient does not have a simple physical interpretation for the axisymmetric primitive equations, and may therefore be considered as arbitrary, values in the range 0.01–0.1 used below correspond to maximum accelerations on the order of 10^{−7}–10^{−6} m s^{−2}.

For both types of experiments, we shall be interested in steady-state solutions of the 2D axisymmetric primitive equations, (24)–(28). Specifically, we want to know whether superrotating solutions exists, and whether multiple solutions may coexist for some values of the forcing parameters. In particular, we shall vary the resonance width parameter *ε* (in the case of the resonant eddy forcing) to illustrate the occurrence of both kinds of bistability identified in section 2. Note that while *ε* was the friction coefficient in the shallow-water model of section 2, we treat it as a free parameter in the numerical experiments below. We shall also discuss the role of viscosity *ν* and the vertical resolution.

### b. Control run

Before investigating bistability, let us first show a control run without eddy forcing (*F*_{u} = 0). The equilibrium zonal wind field is shown in Fig. 6. Jets (with maximum wind speed ≈60 m s^{−1}) are obtained in each hemisphere at the poleward edge of the Hadley cell, which extends approximately to 20° in both hemispheres. Easterly winds prevail in the tropical regions; in particular at the equator, the wind is easterly at all levels. This control run does not exhibit superrotation.

A more realistic control run could be obtained by prescribing additional eddy momentum (or heat) forcing in the midlatitudes (Schneider 1984; Singh and Kuang 2016), as was done in Caballero and Carlson (2018). For simplicity, we prefer not to do so here.

### c. Resonance-driven bistability

In a first set of experiments, using the resonant eddy forcing, we illustrate the type of hysteresis identified in section 2 where bistability is driven by the resonant response to the forcing.

First, we integrate the axisymmetric primitive equations until a statistically stationary state is reached (typically about 1500 days). We show in Fig. 7 (left) the vertically averaged (all the vertical averages shown here are restricted to the region within 1*σ* of the level of maximum forcing *p*_{0}, i.e., to the region between 100 and 500 hPa) zonal wind profile at steady state for a narrow resonance (*ε* = 0.1 day^{−1}), as the forcing amplitude *Q*_{0} is varied. For low values of the forcing amplitude, the zonal wind profile is essentially fixed by angular momentum conservation in the tropics and radiative equilibrium outside. Generally speaking, this state has similar characteristics as the control run: full spatial structure of the zonal wind field, mean meridional circulation. In particular, it exhibits jets close to 20° latitude, as we have seen in the control run. As the forcing amplitude *Q*_{0} increases, these jets move equatorward and weak westerlies appear in the tropics. When *Q*_{0} further increases, there is a relatively sharp transition (between *Q*_{0} = 0.04 and *Q*_{0} = 0.05) to a different regime where a jet appears on the equator, which quickly becomes as strong as the subtropical jets. In this regime, the atmosphere is clearly in a state of equatorial superrotation. The full spatial structure of the wind field in the conventional state is similar to the control run, shown in Fig. 6. The circulation in the superrotating state (shown in the right panel of Fig. 10) will be discussed in more details in section 3e.

We now carry out hysteresis experiments to investigate the possibility that the conventional and superrotating states coexist in some range of forcing amplitude. The experiment consists in increasing the forcing amplitude step by step and letting the system relax to its new equilibrium state at each step. This introduces a discontinuity (in time) in the forcing, but it allows for clearer diagnostics of the response of the system. Typically, we observe a smooth relaxation to a new equilibrium state, possibly with an initial overshoot. As expected, relaxation to the new steady state upon application of the step forcing takes longer close to the bifurcation points. To ensure that the system has relaxed, we choose a time interval between two steps several times longer than the typical relaxation time observed in previous runs. We apply this procedure up to a given value of the forcing amplitude (larger than the amplitude threshold for which we observe the abrupt transition to superrotation in the steady-state experiments above), then we reverse the procedure by decreasing the forcing amplitude step by step until we reach the initial forcing amplitude. Any observable can then be computed as a function of time, or equivalently as a function of the forcing amplitude, with the only difference that in the latter case, it may take one value on the way up and a different one on the way down.

The results of the hysteresis experiments are shown in Fig. 7 (right), for different values of the parameter *ε*. The observable plotted in the figure is the zonal wind, averaged over a range of latitude around the equator (here between 5°S and 5°N) and over the upper atmosphere (between 100 and 500 hPa). For small values of *ε* (narrow resonance; e.g., *ε* = 0.1 day^{−1}), the averaged zonal wind, initially negative, first increases slowly when the forcing amplitude is increased, then abruptly switches to a positive value (above 10 m s^{−1}), characteristic of a superrotating state. This corresponds to the behavior observed with the steady-state experiments in the above paragraph, and suggests that the conventional circulation becomes unstable (saddle-node bifurcation). Once in the superrotating state, the averaged zonal wind again increases slowly with the forcing amplitude until the maximum value of the forcing amplitude is reached. When the forcing amplitude is decreased, the averaged zonal wind decreases slowly, down to a forcing amplitude below the critical point where the conventional circulation became unstable. Then a second bifurcation occurs: the superrotating state becomes unstable and the averaged zonal wind suddenly switches back to its value in the conventional circulation.

The forcing amplitudes corresponding to the bifurcation points depend on *ε*. More precisely, the bistability range decreases significantly as *ε* is increased (see the curves for *ε* = 0.3 day^{−1} and *ε* = 0.5 day^{−1}), that is, as the resonance broadens, as anticipated in section 2. When *ε* is sufficiently large (e.g., *ε* = 0.7 day^{−1}), the bifurcation points disappear entirely: the upper and lower branch of the hysteresis curves collapse onto a single curve, describing the smooth growth of the averaged zonal wind with the forcing amplitude.

### d. Hadley cell–driven bistability

We now turn to the second series of runs, with a constant eddy forcing. Since the resonance mechanism is manually switched off in this case, the only possibility for bistability to occur is through the Hadley cell feedback. The goal is to investigate whether the bistability due to this feedback mechanism, obtained analytically in the simple zonal wind balance model of section 2 and observed in numerical simulations of the 1.5-layer shallow-water equations (SH04), subsists in our multilayer configuration.

Like in section 3c, we start by studying the steady states of the axisymmetric equations: the typical relaxation time from an initial state of rest is similar to the resonant eddy forcing (about 1500 days), although the larger values of *Q*_{0} require longer integrations (up to about 4500 days). Figure 8 shows the vertically averaged zonal wind obtained in the steady state for many forcing amplitudes *Q*_{0}. For both kinds of forcings, the qualitative behavior is similar to the one described in section 3c. For low values of the forcing amplitude, the zonal wind profile is similar to the control run, with subtropical jets at the poleward edge of the Hadley cell. As the forcing amplitude increases, there is a sharp transition to a superrotating circulation. The value of the threshold amplitude is similar in both forcing cases (*Q*_{0} ≈ 0.04), although the magnitude of the equatorial wind in the superrotating state is much larger with the constant forcing.

We now carry out a hysteresis experiment following the same protocol as in section 3c [Fig. 8 (right) shows the averaged zonal wind for this hysteresis experiment]: we increase step by step the forcing amplitude, allowing the system to relax to its new steady state at each time, until an abrupt transition to a superrotating state is found. Then we further increase the forcing amplitude to show that the averaged zonal wind increases smoothly on the superrotating branch, before reverting the loop. We decrease the forcing step by step, until the superrotating state loses stability, and the averaged zonal wind abruptly goes back to the values found on the way up. This shows that the constant eddy forcing also exhibits bistability, as was found in the simple analytical model (section 2) and in a single-layer shallow-water model (SH04). However, it should be noted that the bistability range is much smaller than in the case of resonance-driven bistability studied in section 3c.

Like in the zero-dimensional model of zonal momentum balance studied analytically in section 2, we can diagnose the zonal acceleration budget in our axisymmetric simulations. We show in Fig. 9 the three dominant terms: the eddy forcing (blue curve) as well as the vertical advection of zonal momentum by the Hadley cell, *ω*∂_{p}*u* (orange curve) and the turbulent momentum diffusion term ∇ ⋅ *τ*_{u} (green curve). While the former term is prescribed, the latter terms are dynamically adjusted. These curves are constructed by plotting these terms as functions of the zonal wind, both quantities being averaged over the tropical upper atmosphere, in the hysteresis experiment shown in Fig. 8 (right). In particular, it contains points that correspond to transient states. The vertical transport by the Hadley cell exhibits the same kind of cubic behavior as in the analytical model shown in Fig. 1. Since the dissipative mechanism is not linear friction, the turbulent momentum diffusion curve is not just a straight line, but it is nevertheless an increasing function of the local zonal wind, apart from very low values of the wind. Both mechanisms act essentially as damping effects (again, except for the lower values of the zonal wind as far as turbulent diffusion is concerned). We also display the sum of the two effects as a separate curve (red curve): for these parameter values, there exists a range of wind velocities where the positive feedback of the Hadley cell prevails over the negative feedback of the eddy viscosity, and the net damping is not a monotonous function of the zonal wind. Hence, the qualitative behavior is the same as in Fig. 1: steady-state solutions of the zonal momentum budget should equilibrate this net damping by a prescribed eddy forcing, which is a straight horizontal line in this case. For a fixed forcing amplitude in a given range, bistability may occur. In the figure, we show the prescribed eddy forcing in the hysteresis experiment, where the forcing amplitude is time dependent. This shows that the hysteresis experiment explores successive steady states over the two increasing branches of the net damping curve. The decreasing branch of the net damping curve can only be seen because we have included values from transient states.

### e. Comparing the two types of bistability

There are two major differences between bistability driven by the wave–jet resonance and by the Hadley cell: the behavior of the Hadley cell on the superrotating branch, and the sensitivity to vertical diffusion.

Figure 10 shows the 2D zonal wind field as well as the mean meridional circulation streamfunction for the two kinds of superrotating states: one on the upper branch of the hysteresis loop obtained with a constant eddy forcing (with *Q*_{0} = 0.038), and another on the upper branch of the hysteresis loop obtained with a resonant eddy forcing (*ε* = 0.1 day^{−1}, *Q*_{0} = 0.03). It illustrates the fact that the Hadley cell is almost as strong as in the conventional circulation in the resonance-induced superrotating state, while it is reduced by a factor of 5 in the Hadley cell-induced superrotating state. While in both cases, the equatorial jet is essentially confined to the upper atmosphere, it is much sharper in the case of the constant forcing: both the vertical and the meridional wind shear are larger than in the resonant eddy forcing case. The maximum velocity is also larger with the constant eddy forcing. This is consistent with the behavior of the Hadley cell in the two cases. If we further increase the resonant eddy forcing amplitude (not shown), we recover a state very similar to the superrotating state obtained with the constant eddy forcing at lower forcing amplitude, such as illustrated in the left panel of Fig. 10, with a sharper jet and collapsed Hadley cell. It should be noted that the Hadley cell in the superrotating state may be modified by physical mechanisms not taken into account here, such as eddy heat transport for instance.

Note that in both cases, the zonal flow is maximum near the upper boundary. This is due to a combination of vertical momentum diffusion and the fact that the eddy forcing remains finite there because we used a relatively broad vertical structure (*σ* = 200).

Vertical momentum transport by the eddy viscosity is also expected to play an important role: we have seen in section 2b that, in the 0D model, when bistability is driven by the Hadley cell, it can be destroyed by increasing the strength of dissipative processes. To test whether it is also the case in the 2D axisymmetric primitive equations, we show in Fig. 11 hysteresis experiments for several values of the vertical viscosity *ν*, both for the Hadley-driven and the resonance-driven cases. It is found that the Hadley-driven case exhibits high sensitivity of the stability thresholds of both the conventional and superrotating states (Fig. 11, left). We find that bistability disappears beyond a critical viscosity *ν*_{c} ≈ 0.7 m^{2} s^{−1}. On the other hand, bistability governed by the wave–jet resonance (Fig. 11, right) is much more robust to variations of the vertical diffusivity than the Hadley-driven case: bistability subsists for vertical viscosities up to 2 m^{2} s^{−1}, with an unaffected range of coexistence of the two states. Again, this is in agreement with the theoretical analysis of section 2d, where we have found that in the 0D model resonance-driven bistability subsists when friction is the main damping mechanism.

## 4. Conclusions

In this paper, we have considered the question of atmospheric bistability at the planetary scale through the special case of equatorial superrotation. This case is particularly interesting because it is frequently encountered in planetary atmospheres, and is hypothesized to have played a role in warm climates of the past on Earth. A crucial point is the nature of the transition to superrotation: continuous (akin to second-order phase transitions in condensed matter physics) or abrupt (first-order phase transition). In the latter case, the system exhibits hysteresis. Besides, the transition may even occur spontaneously below the bifurcation point where the conventional state loses stability, driven by the fluctuations inherent to a turbulent atmosphere. The mechanisms determining the nature of the transition may, but need not coincide with those maintaining the equatorial jet by converging angular momentum at the equator. We have studied two such mechanisms corresponding to the two different cases. On the one hand, it was suggested that an abrupt transition to superrotation may be triggered by a resonant response to nonzonal equatorial heating, exciting tropical waves that in turn accelerate the mean flow toward the east. We have shown in a simple model of zonal momentum balance at the equator that such a phenomenon indeed resulted in the appearance of multiple equilibria and hysteresis. On the other hand, it was shown in an idealized framework that the Hadley cell itself could admit two different modes, which results in the coexistence of a conventional and a superrotating state when a constant torque is applied by an external operator. In the same framework, we have studied the interplay between the two mechanisms and showed that there were two main regimes: Hadley-driven bistability, with a small coexistence range, and resonance-driven bistability, with a larger coexistence range. On the other hand, the latter only occurs if the resonance is sufficiently peaked, which in the shallow-water model studied here amounts to sufficiently small linear friction. Parameter values corresponding to the atmosphere of Earth lie close to the boundary separating the two idealized regimes. It should also be noted that, while most existing studies report a significant weakening of the Hadley cell in the superrotating state, our results indicate that in the resonance-driven case it is possible to obtain a superrotating state while maintaining a strong meridional circulation. These findings are confirmed by numerical simulations of an axisymmetric primitive equations model, with an arbitrary number of vertical levels. Nevertheless, we find that Hadley-driven bistability is relatively fragile, in the sense that it depends sensitively on vertical viscosity, while the resonance-driven bistability is much more robust to changes in this parameter.

These results may help shedding light on bistability and hysteresis (or the lack thereof) in full GCM simulations of superrotation. Indeed, while abrupt transitions to superrotation have been reported before, it is often thought that such phenomena should be absent from state-of-the-art models. Here, in the simpler case of a 2D axisymmetric model, we have observed unambiguously the existence of hysteresis phenomena. We have also isolated some factors upon which bistability relies primarily: our simulations provide evidence for a very sensitive dependence on the vertical resolution, and on the damping mechanisms. Further research is still needed to investigate whether the mechanisms described here hold in full 3D GCMs. A critical factor differing from the framework considered here is that the eddy momentum convergence flux should depend dynamically on the full structure of the zonal wind field. Here, we have considered a prescribed forcing, which, although consistent with diagnostics from 3D GCMs, was computed in a linear approximation assuming a uniform background wind, while there may be strong meridional shear in reality, especially in the superrotating state. Understanding how that affects the results reported here would be a key step toward providing a definitive answer to the question of the nature of the transition to superrotation.

## Acknowledgments

This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie Grant Agreement 753021. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement 616811). Part of R. Caballero’s work was performed thanks to a visiting professor grant of ENS de Lyon. Development of the Climt framework (available at https://github.com/CliMT/climt) was supported by the Swedish e-Science Research Centre. Computer time was provided by the “Pôle Scientifique de Modélisation Numérique” in Lyon. We thank the three anonymous referees for useful comments that helped improve the manuscript.

### APPENDIX

#### Numerical Convergence with Vertical Resolution

Bistability driven by the Hadley cell feedback had previously been observed in an axisymmetric model representing only one vertical mode (SH04). In section 3d, we have found that it subsists only marginally in a multilevel setup: there is still bistability but the range of coexistence of the conventional and superrotating states is very narrow. To better understand this behavior, we have performed steady-state and hysteresis experiments with various vertical resolutions. We have found that the meridional profile of vertically averaged zonal wind computed with only two vertical levels is quite similar to the one obtained with full vertical resolution (see Fig. 7, left, and Fig. 8, left). The main differences are that the transition to the superrotating regime seems even sharper in the two-level case, and that stronger equatorial jets are obtained in the resonant case.

Our experiments indicate that the hysteresis loop is quite sensitive to the vertical resolution, and depends on it in a nonmonotonic manner. This holds for both kinds of eddy forcings: resonant and constant. Several hypotheses may be done to account for this sensitivity. First of all, different choices of vertical levels result in different samplings of the vertical profile of the forcing. Given the structure of the forcing, this amounts to multiplying the forcing amplitude by a constant factor. Besides, since the Hadley cell feedback is proportional to vertical shear, poorly resolved vertical gradients may have large effects. Finally, numerical modes of the discretized vertical diffusion operator may also play a part.

Ultimately, the hysteresis loop strongly depends on the number of vertical levels for low resolutions but converges for larger resolutions: for resolutions larger than 10 vertical levels (including a run with 90 levels), the hysteresis cycle does not change significantly. In particular, all the hysteresis loops shown in the above sections have converged.

## REFERENCES

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).