Early warnings of the transition to a superrotating atmospheric state

Several general circulation models (GCMs) show bifurcations of their atmospheric state under a broad range of warm climates. These include some of the more extreme global warming scenarios. This bifurcation can cause the transition to a superrotating state, a state where its angular momentum exceeds the solid body rotation of the planet. Here we use an idealised GCM to simulate this transition by altering a single non-dimensional control parameter, the thermal Rossby number. For a bifurcation-induced transition there is potential for early warnings, and we look for these here. Typically used early warning indicators, variance and lag-1 autocorrelation, calculated for the mean zonal equatorial wind speed, increase and peak just before the transition. The full autocorrelation function taken at multiple lags is also oscillatory, with a period of 25 d preceding the transition. This oscillatory behaviour is reminiscent of a local supercritical Hopf bifurcation. Motivated by this extra structure, we use a generalised early warning vector technique based on principal oscillation patterns (POPs) to diagnose the dominant spatial modes of the horizontal wind field fluctuations. We find a zonal-wavenumber-0 pattern that we call the “precursor” mode that appears shortly before and disappears soon after the transition. We attribute the increase in the early warning indicators to this spatial precursor mode. This mode is correlated to oscillations in strength of the Hadley cells. Following the transition, an eastward-propagating zonal-wavenumber-1 mode of period ∼4 d dominates the dynamics. This mode appears to be representative of the Kelvin–Rossby instability others have previously identified. Although the control parameter used to simulate the transition is unlikely to be relevant to future climate change, the Kelvin–Rossby transition mechanism may well be relevant, and the simulations reported here do show early warnings and serve as a test bed for whether we can detect this transition before it happens.

1 Introduction

Abrupt transitions in the climate often go by the name “tipping points”, and early warnings of such tipping points are currently an active area of research (Scheffer et al., 2009; Lenton, 2011; Scheffer et al., 2012). The IPCC has defined a tipping point as a “critical threshold beyond which a system reorganises, often abruptly and/or irreversibly” (IPCC, 2023) and “for the climate system, the term refers to a critical threshold at which global or regional climate changes from one stable state to another stable state” (IPCC, 2019). One such example is the Atlantic meridional overturning circulation (AMOC), which can abruptly change from an “on” to “off” state as freshwater flux into the North Atlantic is increased (Stommel, 1961; Vellinga and Wood, 2002; Mecking et al., 2016). There are a variety of other postulated climate tipping elements, including the Amazon rainforest, Greenland ice sheet and West Antarctic ice sheet. Changes in atmospheric circulation patterns under anthropogenic climate change are also interesting possible candidates for tipping points, with the 3D nature of atmospheric circulation presenting possibilities of novel early warning signals (EWSs) (Tantet et al., 2015). For example, the West African monsoon is sometimes classified as a tipping system (Loriani et al., 2023), leading to very different vegetation states, e.g. the “green Sahara” (Brovkin et al., 1998). The transition to an atmospheric superrotating state in warm climates is another possible tipping point. Here we look at whether there are EWSs of the tipping point to superrotation in an idealised scenario. However, the extension of EWSs we use is more generally applicable as an approach to detecting precursors of bifurcation-induced changes in state in any spatially extended or coupled system.

Mathematically, tipping points can be described by several kinds of transitions. (a) Noise-induced transitions are probabilistic transitions between multiple steady states induced by large-amplitude noise, with large here meaning that the probability of a perturbation large enough to push the system into a different attractor over the timescale of interest is non-negligible. These events have limited predictability (Ditlevsen and Johnsen, 2010). (b) Rate-induced transitions occur when the disturbance is too fast for the system to remain sitting in its current attractor, resulting in a jump to a new stable state (Ashwin et al., 2012; Wieczorek et al., 2010; Feudel, 2023). (c) Bifurcation-induced transitions occur when the current system state becomes unstable at the bifurcation as the disturbance (or control parameter) is slowly changed. The system will then transition to a new stable state. It is the third type, tipping due to a bifurcation, which has received the most attention, as there is potential for early warnings of this sort of transition (Wiesenfeld, 1985; Held and Kleinen, 2004; Livina and Lenton, 2007; Dakos et al., 2008; Lenton et al., 2012). A bifurcation of states has been reported for the atmosphere when it transitions to a superrotating state (Suarez and Duffy, 1992; Saravanan, 1993; Shell and Held, 2004; Arnold et al., 2012; Herbert et al., 2020).

Figure 1 Schematic of the variables used to describe a superrotating atmospheric state on a planet of radius r rotating with angular velocity Ω. θ is the latitude taken relative to the Equator, and u is the zonal horizontal velocity of the atmosphere in the planet's rotating frame.

An atmospheric superrotating state is said to occur when the atmosphere rotates quicker than the planet beneath it or (more precisely) when the zonal angular momentum of the atmosphere exceeds the angular momentum from the solid-body rotation of the planet. The solid-body angular momentum per unit mass L_E at latitude θ of a planet rotating at angular velocity Ω with radius r is (see Fig. 1)

$$\begin{array}{}\text{(1)}& L_{\mathrm{E}}=\mathrm{\Omega }{r}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }.\end{array}$$

For a shallow atmosphere moving with zonal velocity u at latitude θ measured in the planet's rotating frame, the angular momentum per unit mass L is given by

$$\begin{array}{}\text{(2)}& L=\left(\mathrm{\Omega }+{\displaystyle \frac{u}{r\cos \mathit{\theta }}}\right)r^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }.\end{array}$$

A superrotating state occurs when $L>max\left(L_{\mathrm{E}}\right)$ for some θ. L_E is maximal at the Equator (θ=0), giving $max\left(L_{\mathrm{E}}\right)=\mathrm{\Omega }{r}^{\mathrm{2}}$. For superrotation the zonal velocity of the atmosphere u must be greater than the right-hand side of the following equation:

$$\begin{array}{}\text{(3)}& u>\mathrm{\Omega }r{\displaystyle \frac{{\sin }^{\mathrm{2}}\mathit{\theta }}{\cos \mathit{\theta }}}.\end{array}$$

This condition has a minimum u at the Equator (θ=0) where zonal velocity need only be positive (u>0). Moving away from the Equator ($\left|\mathit{\theta }\right|>\mathrm{0}$) the zonal velocity required for superrotation increases steeply. If it occurs, superrotation will almost certainly be present at the Equator as superrotation at a higher latitudes would result in inertially unstable atmospheric configurations. What this means for Earth rotating west to east is that the prevailing easterlies at the Equator (u<0) change to prevailing westerlies. We label the bound for superrotation u_SR, which is found using the following equation:

$$\begin{array}{}\text{(4)}& u_{\mathrm{SR}}=\mathrm{\Omega }r{\displaystyle \frac{{\sin }^{\mathrm{2}}\mathit{\theta }}{\cos \mathit{\theta }}}.\end{array}$$

Superrotating atmospheres exist on Venus and Titan (Read and Lebonnois, 2018), and they may have occurred on Earth in warmer, past climates. Caballero and Huber (2010) showed a full-complexity general circulation model (GCM) that transitioned to a superrotating state when equatorial surface temperatures got sufficiently high, within the range of temperatures thought to have existed in the early Cenozoic period of Earth's history (65 million years ago). This was a period where it was around 10 °C warmer on average than today. The atmospheres of other warmer periods such as the early Pliocene and Eocene may also have been superrotating (Tziperman and Farrell, 2009). Simulations of warm climates that generate enhanced tropical convective variability seem particularly prone to superrotation (Arnold et al., 2012), as do those where convective heating strengthens (Laraia and Schneider, 2015). This has led to the hypothesis that it may be relevant to future climate projections (Held, 1999; Pierrehumbert, 2000) and account for the superrotation seen in extreme global warming simulations (Huang et al., 2001). Even though the probability of transitioning to a surface superrotating state is presently thought to be low as the Earth warms into the future (Caballero and Carlson, 2018), the possibility is a high risk. For example, if the prevailing Pacific equatorial westerlies switched to the easterlies of an atmosphere with superrotation at the surface, a permanent El Niño state would result, along with all of the global impacts that this would entail (Pierrehumbert, 2000). Tziperman and Farrell (2009) used this mechanism to explain the permanent El Niño state that is thought to have existed in the Pliocene. Although surface superrotation seems unlikely, superrotation happens readily in the upper troposphere under high-end global warming scenarios. The impact on surface climate would, however, be less severe (Caballero and Carlson, 2018).

Superrotation implies the atmospheric angular momentum locally exceeds that provided by the planet surface. This can only be achieved through an up-gradient momentum flux that is in excess or balances the removal of momentum by the poleward large-scale meridional circulation and the zonal-mean drag due to viscosity (Caballero and Carlson, 2018). This up-gradient flux can be supplied by non-axisymmetric eddies (Hide, 1969; Gierasch, 1975; Read, 1986). There are many studies that examine the eddy dynamics responsible for superrotation, of which there are two broad categories, although horizontal momentum transport by off-equatorial waves feature in both (i) generation through non-axisymmetric forcing (Scott and Polvani, 2008; Saravanan, 1993; Suarez and Duffy, 1992; Arnold et al., 2012; Shell and Held, 2004) and (ii) generation through barotropic, baroclinic or Kelvin–Rossby instability (Williams, 2003; Wang and Mitchell, 2014; Mitchell and Vallis, 2010; Zurita-Gotor and Held, 2018; Zurita-Gotor et al., 2022). It is the second category that is relevant to this paper, in particular the Kelvin–Rossby instability mechanism. Theoretical work on this mechanism (Wang and Mitchell, 2014; Zurita-Gotor and Held, 2018; Zurita-Gotor et al., 2022) shows the transition is via a loss of stability and so there should be early warning signals. In this mechanism, superrotation arises from the wave–mean-flow interaction, specifically the interaction of equatorial Kelvin waves and off-equatorial Rossby waves that mutually amplify each other and generate a planetary-scale mode that converges zonal momentum onto the Equator. Zurita-Gotor et al. (2022) show that the Kelvin–Rossby coupling is a viable mechanism for producing terrestrial superrotation and may explain the warm-climate superrotation found by Caballero and Huber (2010).

Transitions to superrotating states may be smooth (Huang et al., 2001; Caballero and Huber, 2010; Mitchell and Vallis, 2010) or abrupt (Suarez and Duffy, 1992; Saravanan, 1993; Arnold et al., 2012; Zurita-Gotor et al., 2022). The smooth, gradual type seems to be more common in simulations of global warming scenarios or warm paleoclimates with comprehensive, multilevel GCMs (Arnold et al., 2012). Abrupt transitions are seen more commonly in models with some degree of idealisation; however, that does not rule them out in more realistic scenarios (Held, 2016). The original example of an abrupt transition to a superrotating state was found by Suarez and Duffy (1992) in a two-layer GCM as the tropical heating was increased. This transition also showed hysteresis, meaning that as the tropical heating was reversed below the transition's critical value, the superrotating state persisted to a new and lower critical value. Arnold et al. (2012) have found similar behaviour in a higher-resolution and multilevel (albeit still idealised) GCM. This bistability is something we will also investigate in a different setting in this paper.

Tipping points and bifurcation-induced transitions are often, but not necessarily, abrupt (Armstrong McKay et al., 2022). Transcritical, supercritical pitchfork and supercritical Hopf bifurcations (Strogatz, 2001) are examples that result in smooth transitions but have early warnings due to the loss of stability of the local state. The transition to superrotation reported here is smooth. Nevertheless, early warnings of the transition are shown, and therefore these simulations serve as a test bed for whether we can detect this transition before it happens.

In this paper the focus will be the fluctuations around the mean wind field state and how their properties change as the atmosphere gets closer to a superrotating state, the idea being that these properties can serve as early warnings of the transition. We first introduce the idealised GCM used for these numerical simulations and the control parameter we vary to span Earth-like to fully superrotating atmospheric states in Sect. 2. We then describe this range of simulated mean atmospheric states in Sect. 3 before going to the main focus of the paper: early warnings of the transition. We start by calculating the widely used early warning indicators in Sect. 4. Motivated by the bistability in the superrotating states observed in the experiments of Suarez and Duffy (1992) and others, we look for multiple co-existing states with our control parameter in Sect. 5. In Sect. 6, we generalise the usual scalar early warning indicators to the vector setting by calculating the principal oscillation patterns (POPs). This allows us to diagnose the dominant spatial patterns in the wind field fluctuations preceding the superrotation transition and the dominant patterns following the transition. We discuss our results and make conclusions in Sect. 7.

2 Methods

We use an idealised general circulation model (GCM) to simulate the transition to superrotation. The model framework we use is Isca (Vallis et al., 2018), which is based on a spectral dry dynamical core obeying the primitive equations of motion of Gordon and Stern (1982). Isca is designed to be highly configurable for use in exoplanet research in terms of both complexity (multiple radiation schemes and moist or dry atmospheres are options) and parameter values (planetary radius, rotation rate and fluid density are among the many values that can be altered). The configuration we use here is as simplified as possible and follows the approach of Mitchell and Vallis (2010) closely. Simulations are performed at T42 resolution (about 2.8° resolution at the Equator) and with 25 vertical sigma levels (equal pressure levels $\mathit{\sigma }=p/p_{\mathrm{s}}$, where p is pressure and p_s is surface pressure).

The surface of the planet is prescribed to be isotropically smooth with no topography or continents present. No radiation scheme is used either. Instead temperature is linearly relaxed to a prescribed mean annual average state (Newtonian cooling) using the scheme of Held and Suarez (1994). In particular, the surface temperature is relaxed to the zonally symmetric profile

$$\begin{array}{}\text{(5)}& T_{\mathrm{0}}=\stackrel{\mathrm{‾}}{T}[\mathrm{1}+{\displaystyle \frac{{\mathrm{\Delta }}_H}{\mathrm{3}}}(\mathrm{1}-\mathrm{3}{\sin }^{\mathrm{2}}\mathit{\theta }\left)\right],\end{array}$$

with T₀ being the lowest-level temperature and $\stackrel{\mathrm{‾}}{T}$ the global average surface temperature. Δ_H is a non-dimensional parameter specifying the pole–Equator temperature difference. The vertical temperature is assumed to relax to a moist adiabat with lapse rate 6 K km⁻¹ capped at a minimum of 200 K in the stratosphere. The radiative relaxation time is 40 d in the free troposphere and approximately 4 d in the boundary layer, the top of which is fixed at σ=0.7 (see Held and Suarez, 1994, for the full functional forms). As previously mentioned, we follow Mitchell and Vallis (2010) closely, although there are some minor differences: Mitchell and Vallis (2010) use $\stackrel{\mathrm{‾}}{T}=\mathrm{315}$ K and Δ_H=0.2, here we use $\stackrel{\mathrm{‾}}{T}=\mathrm{314}$ K and Δ_H=0.215. As both of these parameters combine together in a non-dimensional number that will be used as the control parameter in numerical experiments (the thermal Rossby number; see below), these small differences should be negligible in comparing our results with theirs.

Friction from the planet's surface on the atmosphere relaxes horizontal velocities back towards zero in the atmospheric boundary layer, and this is specified as an extra, linear and vertical-height-dependent (Rayleigh) drag term in the equation for the horizontal velocity. This extra drag term is applied as

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial \mathit{v}}{\partial t}}=\mathrm{\cdots }-k_v\left(\mathit{\sigma }\right)\mathit{v},\text{(7)}& {\displaystyle }{k}_v\left(\mathit{\sigma }\right)=k_{f}max\left(\mathrm{0},{\displaystyle \frac{\mathit{\sigma }-{\mathit{\sigma }}_{\mathrm{b}}}{\mathrm{1}-{\mathit{\sigma }}_{\mathrm{b}}}}\right),\end{array}$$

where v is a vector of the horizontal velocities $:=(u,v)$. For the full horizontal momentum equations, see Mitchell and Vallis (2010). The reciprocal of the Rayleigh relaxation time, k_f, is taken to be 0.5 d⁻¹ at the planet surface (σ=1) and linearly decreases to zero at the top of the boundary layer at σ_b=0.7. This is slightly smaller than the value of 1 d⁻¹ used in Mitchell and Vallis (2010).

The horizontal momentum equations governing the atmospheric dynamics simplify in non-dimensional form when subject to Rayleigh friction and Newtonian cooling. These equations are then governed by just a few non-dimensional parameters. The relevant non-dimensional parameter that is varied in the experiments here and in Mitchell and Vallis (2010) is the thermal Rossby number:

$$\begin{array}{}\text{(8)}& R{o}_T={\displaystyle \frac{R{T}_{\mathrm{0}}{\mathrm{\Delta }}_H}{(\mathrm{2}\mathrm{\Omega }r{)}^{\mathrm{2}}}}.\end{array}$$

R here is the specific gas constant for dry air (287 J K⁻¹ kg⁻¹). However, as the numerical model code is dimensional, a dimensional parameter must be varied, and here the planetary radius r has been chosen for this purpose. By making this choice, the other non-dimensional parameters, the Ekman number and thermal relaxation time, remain unaltered.

Strictly Ro_T as given by Eq. (8) varies with latitude. However, in the following experiments we attribute a single value of Ro_T to each simulation by replacing T₀ with the average temperature $\stackrel{\mathrm{‾}}{T}$, i.e. $T_{\mathrm{0}}\to \stackrel{\mathrm{‾}}{T}$ in Eq. (8), which is the average thermal Rossby number over the whole planet.

To summarise, in the following experiments we vary the thermal Rossby number to go from Earth-like atmospheric states through the transition to a fully superrotating state by changing the planetary radius r for each simulation. All other parameters will remain fixed. Each simulation is initialised from a homogenous atmospheric state and integrated for either 3 or 10 years (1080 or 3600 d, respectively). A stationary state is reached after approximately 200 d. In the following results we use the last 2 years of the 3-year simulation unless stated otherwise.

3 Mean atmospheric state with varying Ro_T

In Fig. 2, we plot the mean zonal and vertical structure of the atmosphere at selected values of Ro_T. The coloured line contours show the zonal time mean wind speed u with latitude θ and height $\mathit{\sigma }=p/p_{\mathrm{s}}$ after the atmosphere has reached a stationary state (the last 2 years of a 3-year integration). The red line contours show regions of u>0, the blue line contours show regions of u<0 and the black line contour marks u=0. Line contours are spaced at 5 m s⁻¹ intervals. The region enclosed by the dotted black line is superrotating, meaning that it marks the region where u=u_SR (Eq. 4) and is centred at the Equator. The filled contours show the mean meridional overturning.

We start simulations at Earth-like values of the control parameter (Ro_T=0.02, Fig. 2a) and incrementally increase it through the transition. No superrotation at any vertical level is present in the Earth-like simulation. The transition occurs first in the upper troposphere near the tropopause. This first transition is around Ro_T∼0.07 (Fig. 2c, dotted black region at the Equator around σ=0.2). We will refer to this as the first transition to an equatorial superrotating state. This small superrotating region then extends further downwards into the troposphere, upwards in the stratosphere and out into higher latitudes as Ro_T is increased (Fig. 2d–f).

Superrotation is present at Ro_T=0.87 (Fig. 2e) in the majority of the equatorial troposphere and lower stratosphere above the height of the boundary layer (σ=0.7) but absent within it. For Ro_T>0.87 (Fig. 2f), the third-lowest vertical level in the model also starts to superrotate at the Equator. For all simulations we have investigated (Ro_T≤11.4), this is the lowest vertical level that superrotates and marks the final vertical transition. We will refer to this as the final transition to an equatorial superrotating state. All of the troposphere at σ<0.8 (and some of the stratosphere) is superrotating after this point. We study this final transition at the lowest altitude that becomes superrotating in some detail in this paper. This lowest vertical level, the third vertical level is at σ=(0.7 0.79], with a mid-point at σ=0.74, and this level lies at a vertical height of ∼1–2 km above sea level.

As previously mentioned, none of the simulations show superrotation in the two lowest vertical levels (σ>0.8). This could be due to Rayleigh drag within the boundary layer (σ≥0.7) decelerating the horizontal flow and acting as a negative angular momentum flux for any u>0. However, the final transition occurs at a level within the atmospheric boundary layer, and thus it still experiences a small amount of drag. However, surface superrotation has been shown to occur in small parameter regimes (Caballero and Carlson, 2018) and is not commonly seen in GCM simulations.

Other notable features are the midlatitude jets (the largest-value red line contours in Fig. 2) in the Earth-like simulation centred at σ∼0.3 and at latitude $\mathit{\theta }\sim \pm \mathrm{50}$°. These jets tend to move poleward as Ro_T increases. The mean meridional overturning circulation (filled contours in Fig. 2) also tends to reduce in intensity with increasing Ro_T.

Figure 2 Mean vertical and zonal structure of the atmosphere at various values of the control parameter Ro_T. Zonal mean u with latitude and vertical height $\mathit{\sigma }=p/p_{\mathrm{s}}$ is shown as coloured line contours: u>0 contours are shown as red lines, the u=0 contour is the black line and u<0 contours are shown as blue lines. Contour lines are spaced at 5 m s⁻¹ intervals. The region enclosed by the dotted black line is superrotating (c–f). This means that it marks the region where u=u_SR (Eq. 4). The filled contours show the mean meridional overturning as given by the mass streamfunction, with blue and red regions having anti-clockwise and clockwise circulation, respectively. An Earth-like state is shown in (a) (Ro_T=0.02). (b) A state close to the first superrotation transition. (c) A small region in the equatorial upper troposphere that becomes superrotating. This region expands further down into the troposphere and up into the stratosphere as Ro_T increases (d–f). Panel (e) shows the zonal-mean atmospheric state just before the third-lowest vertical level becomes superrotating. This is the lowest level that superrotates, at least for Ro_T∼O(1), and it marks the final transition. In (f) all of the troposphere σ<0.8 is superrotating.

We also show snapshots t=900 d into the simulation of the horizontal wind field ($\mathit{v}=(u,v)$) at the two vertical levels where the first and final transitions occur (σ=0.19 and σ=0.74, respectively) for various key values of the control parameter Ro_T in Fig. 3. Earth-like simulations are shown in Fig. 3a and b. The midlatitude instabilities are evident, as is the lack of superrotation at the Equator (blue shading). In Fig. 3c, the wind field just before the first transition is shown at the vertical level it occurs at (Ro_T=0.06, σ=0.19). There is a roughly equal mixture of red and blue shading giving a zero mean zonal u. Lower levels are not showing equatorial superrotation at this value of Ro_T (Fig. 3d). Figure 3f shows the wind field just before the final transition at the vertical level it occurs at (Ro_T=0.87, σ=0.74). Again, there is a roughly equal mixture of red and blue shading, giving zero mean zonal u. Figure 3e, g and h are all fully superrotating equatorial states (mainly red shading).

Figure 3 Snapshots (t=900 d into each simulation) of the horizontal $\mathit{v}=(u,v)$ wind field (shown as arrows) at the two vertical levels σ=0.19 and σ=0.74. These vertical levels correspond to the first transition to superrotation at Ro_T∼0.07 and the final transition at Ro_T≳0.87, respectively. The underlying colour plot shows the u component of this field, with reds corresponding to u>0 and blues to u<0 (u=0 is white). Panels (a) and (b) show an Earth-like state. Panel (c) shows the wind field just before the first transition at the vertical level it occurs at. Panel (f) shows the wind field just before the final transition at the vertical level it occurs at. Panels (e), (g) and (h) show fully superrotating equatorial wind field states (predominately red colours at the Equator).

More information on the mean states and the mechanisms of superrotation in these simulations is reported in Mitchell and Vallis (2010). From here on we get to the main focus of this paper, namely that we study the fluctuations around the mean states and their properties as precursors to the transitions.

4 Temporal early warning signals

In this section we calculate early warning signals (EWSs) for the varying Ro_T simulations with the aim of early detection of the transition to a superrotating state. These indicators are widely used to detect abrupt transitions resulting from local bifurcations (Lenton, 2011; Thompson and Sieber, 2011). Typically, these are designed to detect the increasing recovery time and amplitude of small fluctuations δx(t) around the mean stationary state, denoted as $\stackrel{\mathrm{‾}}{x}$, of some scalar state parameter x(t) that is a function of time t. The state of the system is therefore given by

$$\begin{array}{}\text{(9)}& x\left(t\right)=\stackrel{\mathrm{‾}}{x}+\mathit{\delta }x\left(t\right),\end{array}$$

and it is often assumed to evolve autonomously according to some function f(x) as $\dot{x}=f\left(x\right)$ (overdots denote ordinary time derivatives, i.e. $\dot{x}:=\frac{\mathrm{d}x}{\mathrm{d}t}$). f(x) is a non-linear function of x in general. Provided the fluctuations are not too large, one can approximate the dynamics well by performing Taylor expansion on f(x) to the first order around the mean state, i.e. $f\left(x\right)\approx f\left(\stackrel{\mathrm{‾}}{x}\right)+J\left(\stackrel{\mathrm{‾}}{x}\right)\mathit{\delta }x$, where the Jacobian is given by $J\left(\stackrel{\mathrm{‾}}{x}\right)=\frac{\partial f\left(x\right)}{\partial x}{|}_{x=\stackrel{\mathrm{‾}}{x}}$. In this case the dynamics of the fluctuations can be approximated with a linear first-order ordinary differential equation (ODE), i.e.

$$\begin{array}{}\text{(10)}& \dot{\mathit{\delta }x}\approx J\left(\stackrel{\mathrm{‾}}{x}\right)\mathit{\delta }x.\end{array}$$

For a time-invariant mean state $\stackrel{\mathrm{‾}}{x}$ (and time-invariant J(x)), this simple ODE has the solution

$$\begin{array}{}\text{(11)}& \mathit{\delta }x\left(t\right)=e^{J\left(\stackrel{\mathrm{‾}}{x}\right)t}\mathit{\delta }x\left(\mathrm{0}\right).\end{array}$$

There are three possible fates of an initial fluctuation in amplitude δx(0) depending on the sign of the real part of $J\left(\stackrel{\mathrm{‾}}{x}\right)$. (i) The real part of $J\left(\stackrel{\mathrm{‾}}{x}\right)$ is negative ($\mathcal{R}\left(J\right(\stackrel{\mathrm{‾}}{x}\left)\right)<\mathrm{0}$), the initial fluctuation δx(0) decays and the system state will be stable to perturbations eventually recovering to $\stackrel{\mathrm{‾}}{x}$. In this case a recovery or e-folding time can be defined as the time taken, τ, for the fluctuation to reach $\mathrm{1}/e$th of its initial value as $\mathit{\tau }=-\mathrm{1}/\mathcal{R}\left(J\right(\stackrel{\mathrm{‾}}{x}\left)\right)$. (ii) The real part of $J\left(\stackrel{\mathrm{‾}}{x}\right)$ is positive ($\mathcal{R}\left(J\right(\stackrel{\mathrm{‾}}{x}\left)\right)>\mathrm{0}$), any initial fluctuation will get exponentially larger with time and the system will not recover to $\stackrel{\mathrm{‾}}{x}$, i.e. $\stackrel{\mathrm{‾}}{x}$ is unstable. (iii) The real part of $J\left(\stackrel{\mathrm{‾}}{x}\right)$ is zero, i.e. $\mathcal{R}\left(J\right(\stackrel{\mathrm{‾}}{x}\left)\right)\to \mathrm{0}$. In this case, perturbations will not grow or decay. This neutral stability occurs at a bifurcation.

Early warning indicators are designed to detect the increasingly less negative value of $\mathcal{R}\left(J\right(\stackrel{\mathrm{‾}}{x}\left)\right)$ or equivalently the increasing system recovery time τ as it approaches the bifurcation. This phenomenon is termed “critical slowing down”. One popular indicator is the variance Var(x) of x with time

$$\begin{array}{}\text{(12)}& {\displaystyle }\text{Var}\left(x\right)& {\displaystyle }=E\left[\right(x\left(t\right)-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}],\text{(13)}& {\displaystyle }& {\displaystyle }=E\left[\mathit{\delta }x\right(t{)}^{\mathrm{2}}],\end{array}$$

where E[.] is the expectation value. For the case where x is a scalar variable with first-order dynamics subject to random (white-noise) forcing (an Ornstein–Uhlenbeck process), J(x) is also a real scalar, and the variance is given by

$$\begin{array}{}\text{(14)}& \text{Var}\left(x\right)={\displaystyle \frac{{\mathit{\sigma }}_Q^{\mathrm{2}}}{\mathrm{2}\left|J\right(\stackrel{\mathrm{‾}}{x}\left)\right|}}.\end{array}$$

This measures the mean amplitude of the system fluctuations ${\mathit{\sigma }}_x=\sqrt{\text{Var}\left(x\right)}$ to a random forcing of mean amplitude σ_Q (σ_Q is generally assumed to be constant as the bifurcation is approached). The state $\stackrel{\mathrm{‾}}{x}$ becomes less stable closer to the bifurcation, and $J\left(\stackrel{\mathrm{‾}}{x}\right)$ gets less negative, showing up in the indicator as an increase in variance.

White-noise forcing is often assumed to approximate internally generated fast-timescale chaotic weather variability when applied to climate applications, such as early warning indicators, empirical orthogonal functions and linear regressions (von Storch and Zwiers, 1999), and this is also what we do throughout this paper.

Another popular early warning indicator is time-lagged autocorrelation

$$\begin{array}{}\text{(15)}& \mathit{\alpha }\left(t_{\mathrm{lag}}\right)=e^{J\left(\stackrel{\mathrm{‾}}{x}\right)\left|t_{\mathrm{lag}}\right|},\end{array}$$

which measures the recovery time of the system to perturbations. Again, for the case where x is a real scalar variable with first-order dynamics subject to random (white-noise) forcing, $J\left(\stackrel{\mathrm{‾}}{x}\right)$ is a negative real number if the state is stable, and α(t_lag) is just an exponentially decaying function of t_lag. As a bifurcation is approached $\mathcal{R}\left(J\right(\stackrel{\mathrm{‾}}{x}\left)\right)\to \mathrm{0}$ and the indicators σ_x→∞ and α(t_lag)→1.

The temporal evolution of a scalar state variable x described by a linear first-order ODE is constrained to only permit exponential decay from a perturbation. However, a real-valued state variable x may have oscillatory behaviour and therefore must be described by higher-order ODEs. In this case α(t_lag) may have a more complicated functional form. For example, a system with second-order dynamics can have an autocorrelation function with an oscillatory term. Relevant to our case is the approach to a Hopf bifurcation, whose α(t_lag) has the generic form (Bury et al., 2020)

$$\begin{array}{}\text{(16)}& \mathit{\alpha }\left(t_{\mathrm{lag}}\right)=e^{-\frac{\left|t_{\mathrm{lag}}\right|}{\mathit{\tau }}}\cos \left(\mathit{\omega }{t}_{\mathrm{lag}}\right).\end{array}$$

For a real-valued J(x) to describe a Hopf bifurcation, it must be at least a two by two (non-symmetric) matrix with complex eigenvalues, λ and λ^*, that appear in conjugate pairs (the * denotes complex conjugation). The real part of the complex eigenvalues of J(x) ($\mathit{\tau }=-\mathrm{1}/\mathcal{R}\left(\mathit{\lambda }\right)$) determines the exponential decay envelope of the autocorrelation function. This envelope decays over a longer timescale as the bifurcation is approached. Within this envelope there may also be oscillations, whose frequency is given by the imaginary part of J(x), ω=ℐ(λ).

In most studies, the full autocorrelation function at all values of t_lag is generally not used. It is usually the value of the autocorrelation function at t_lag=1 (the lag-1 autocorrelation) that is used. This often works well as it is more robust to sample size and signal to noise provided the oscillation period, $P=\frac{\mathrm{2}\mathit{\pi }}{\mathit{\omega }}$, is large compared to the value of t_lag that is used to evaluate α(t_lag). Formally, if ωt_lag≪1, then Eq. 16 approximates to $\mathit{\alpha }\left(t_{\mathrm{lag}}\right)\sim e^{-\frac{\left|t_{\mathrm{lag}}\right|}{\mathit{\tau }}}$; i.e. it is a good indicator of the real part of the Jacobian eigenvalues and therefore of critical slowing down.

We will also use the autocovariance function, R(t_lag), which encodes information of both variance and autocorrelation, giving

$$\begin{array}{}\text{(17)}& R\left(t_{\mathrm{lag}}\right)=E\left[\right(x(t+t_{\mathrm{lag}})-\stackrel{\mathrm{‾}}{x}\left)\right(x\left(t\right)-\stackrel{\mathrm{‾}}{x}\left)\right],\end{array}$$

or $R\left(t_{\mathrm{lag}}\right)=\mathrm{Var}\left(x\right)e^{J\left(\stackrel{\mathrm{‾}}{x}\right)\left|t_{\mathrm{lag}}\right|}$. The variance is just given by R(0) and $\mathit{\alpha }\left(t_{\mathrm{lag}}\right)=\frac{R\left(t_{\mathrm{lag}}\right)}{R\left(\mathrm{0}\right)}$.

4.1 Temporal EWSs of the superrotation transition

We apply these early warning indicators to study the transition to superrotation at vertical level σ as follows. We look at the statistics of the fluctuations, δu(σ,t), around the mean zonal wind field state, $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$, where $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ is the time- and area-weighted spatial mean of a spherical segment centred on the Equator at vertical height σ with edges at ${\mathit{\theta }}_{\pm }=\pm \mathrm{10}$° latitude. An equatorial band centred on the Equator is chosen to analyse the zonal wind field u because, as previously discussed, this is the most likely place to observe superrotation if it occurs. Formally, we compute R(σ,t_lag) for $\mathit{\delta }u(\mathit{\sigma },t)=u(\mathit{\sigma },t)-\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ where

$$\begin{array}{}\text{(18)}& u(\mathit{\sigma },t)={\displaystyle \frac{r}{A}}\underset{\mathit{\varphi }=\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}\mathit{\varphi }\underset{\mathit{\theta }={\mathit{\theta }}_{-}}{\overset{{\mathit{\theta }}_{+}}{\int }}u(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}$$

where ϕ is the longitude and $A=r{\int }_{\mathit{\varphi }=\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}\mathit{\varphi }{\int }_{\mathit{\theta }={\mathit{\theta }}_{-}}^{{\mathit{\theta }}_{+}}\cos \mathit{\theta }\mathrm{d}\mathit{\theta }$ is the area of the spherical segment. We choose a section of the time series u(σ,t) starting at a time t_s that is sufficiently after the simulation spinup so that the atmosphere is in an approximately stationary state. In the model configuration used, the atmosphere reaches a steady state after roughly 200 d, although we use t_s=361 d. The finish time t_f is chosen to be large for accurate determination of the statistical estimators. Ideally this would be t_f→∞; however, due to computational time we mostly limit this to t_f=1080 d. However, when studying the first transition we have performed longer simulations (10 years) to better resolve the long timescales present around the tropopause and set t_f=3600 d; u(σ,t) is then averaged over a period $T=t_{\mathrm{f}}-t_{\mathrm{s}}$ to give

$$\begin{array}{}\text{(19)}& \stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}={\displaystyle \frac{\mathrm{1}}{T}}\underset{t_{\mathrm{s}}}{\overset{t_{\mathrm{f}}}{\int }}u(\mathit{\sigma },t)\mathrm{d}t.\end{array}$$

In Fig. 4, $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ is shown as it varies with vertical height σ and the control parameter Ro_T.

Figure 4 Here we show $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ in an equatorial band with edges at ±10° latitude as it varies with vertical height σ and control parameter Ro_T. Red contours show positive $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$. Blue contours show negative $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$. Both are spaced at 5 m s⁻¹ intervals. The black line denotes the zero contour and approximately marks the boundary between superrotating and non-superrotating regions.

Figure 5 Mean u around an equatorial band at vertical height σ=0.74. (a) Mean zonal wind speed u(σ,t) in a spherical segment centred on the Equator and with edges at ±10° latitude during a 3-year simulation (1080 d). Different colours correspond to a simulation with a different value of Ro_T. We also show highlighted values of Ro_T∈{0.02 (blue) 0.87 (red) 1.93 (green)}. In (c) we show $\mathit{\delta }u(\mathit{\sigma },t)=u\left(t\right)-\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$, which shows the fluctuations in u(σ,t) around the temporal mean $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ over the last 2 years of the simulation as a function of Ro_T. Dotted lines indicate $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$. (b) The autocovariance function R(σ,t_lag) as a function of t_lag. This function is oscillatory for RoT=0.87 (red line) just before the final transition.

4.1.1 Final transition: Ro_T∼0.87 at σ=0.74

The final transition where the majority of the equatorial troposphere (σ<0.8) starts to superrotate occurs at Ro_T∼0.87. In this subsection the statistics of δu(σ,t) at σ=0.74 are analysed. This corresponds to the third-lowest vertical level in the model. We look for early warnings of the first transition in the following subsection.

In Fig. 5a, we plot the $u(\mathit{\sigma }=\mathrm{0.74},t)$ time series for selected values of Ro_T (0.02 is blue (Earth-like), 0.87 is red (close to the final transition) and 1.93 is green (fully superrotating)) over the 3-year simulation period. The mean u(σ,t) over the last 2 years ($\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$) is plotted as a corresponding dotted coloured line. In Fig. 5c we show the fluctuations around the mean, δu(σ,t), as a function of Ro_T. The larger, oscillatory response of u(σ,t) close to the final transition is clearly visible in the time series. The time series looks much more random and akin to first-order noisy dynamics (such as an Ornstein–Uhlenbeck process) on either side of the final transition (blue and green).

In Fig. 5b we have plotted the autocovariance function (Eq. 17), $R(\mathit{\sigma }=\mathrm{0.74},t_{\mathrm{lag}})$, for each of the three selected Ro_T values during the last 2 years of each simulation. The value at t_lag=0 is equivalent to the variance of u, Var(u). One can see close to the transition (red) that variance is much higher than either side of the transition (green and blue). One can also see the functional form of R(σ,t_lag) is also very different close to the transition, showing an oscillation of around 25 d contained in an exponential decay envelope $\left(e^{-\frac{\left|t_{\mathrm{lag}}\right|}{\mathit{\tau }}}\right)$ with an e-folding time of approximately τ∼35 d (see Eq. 16). In contrast, the response of u(σ,t) to perturbations is neither oscillatory nor long lasting on either side of the final transition.

Figure 6 Early warning indicators of the final transition at σ=0.74. Each dot is an individual simulation. Coloured dots correspond to the same coloured simulations in Fig. 5. (a) Mean zonal wind speed $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ in a spherical segment centred on the Equator with edges at ±10° latitude during the last 2 years of the simulation against RoT at vertical height σ=0.74. The dotted red line is the threshold for superrotation to occur (u_SR). This is non-zero as the wind speed is averaged over a non-zero-thickness latitude band. (c) Typical early warning indicators of a local bifurcation: with the left-hand plot showing the variance of mean zonal u(σ,t) and the right-hand plot showing lag-1 autocorrelation of mean zonal u(σ,t). (b) Autocovariance function R(σ,t_lag) and how it changes with Ro_T (y axis) and t_lag. Oscillations of the 25 d period appear for values of RoT∼[0.5 1.2] and have a maximum in amplitude around RoT∼0.9.

In Fig. 6a we plot $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ in the equatorial spherical segment against Ro_T. Black dots are individual simulations, and black lines connect these values. The dotted red line is the threshold $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ required for superrotation, u_SR; see Eq. (4). This is larger than the zero bound required for superrotation at the Equator due to the finite latitude width of the spherical segment. Superrotation at this vertical level occurs approximately for Ro_T>1. Interestingly, this plot is not monotonically increasing with Ro_T – there is a local maximum at Ro_T∼0.54 followed by a plateau to Ro_T∼0.75 below the superrotation threshold. There is another local (and global) maximum in $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ around Ro_T∼1.9 when the atmosphere is in a fully superrotating state. In Fig. 6c we plot the usual early warning indicators, variance of u(σ,t) (left-hand plot, Var(u)) and lag-1 autocorrelation (right-hand plot, α(t_lag=1 d)). Both indicators peak at Ro_T=0.87, just below $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}\sim \mathrm{0}$, suggesting this is where the atmospheric state at this vertical level is least stable.

Figure 7 Mean u around an equatorial band at vertical height σ=0.19. (a) Mean zonal wind speed u(σ,t) in a spherical segment centred on the Equator and with edges at ±10° latitude during a 10-year simulation (3600 d). Different colours correspond to a simulation with a different value of Ro_T. We also show highlighted values of Ro_T∈{0.02 (blue) 0.08 (red) 0.12 (green)}. In (c) we show $\mathit{\delta }u(\mathit{\sigma },t)=u\left(t\right)-\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$, which gives the fluctuations in u(σ,t) around the temporal mean $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ over the last 9 years of the simulation as a function of Ro_T. Dotted lines indicate $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$. (b) The autocovariance function R(σ,t_lag) as a function of t_lag. All of these functions have slow decay timescales, indicating very lightly restored or unstable regions, as can be seen from the meandering time series in (a) and (c).

There is an initially high value of lag-1 autocorrelation when Ro_T=0.02 is at an Earth-like value, and this looks like an outlier when compared to the other simulations in Fig. 6c. The full (multiple lag) autocorrelation function of this simulation as shown in Fig. 5b in blue is not oscillatory or with high variance, but it does not look like a typical exponentially decaying function like the other simulations performed here. It has an approximate linear decay with increasing t_lag. As the autocorrelation is assumed to be an exponentially decaying function of t_lag, lag-1 autocorrelation may not be the best way to infer the decay timescale (and critical slowing down) for this particular value of Ro_T. The mean equatorial zonal flow may be in a different regime compared to the other simulations as Ro_T is increased.

As shown in Fig. 5b, the autocorrelation function is not necessarily only exponentially decaying with increasing t_lag as required for lag-1 autocorrelation to be a good indicator of the stability of the state. However, as argued in Sect. 4, it is a good indicator for an oscillatory autocorrelation function provided ωt_lag is much smaller than 1. In this case $\mathit{\omega }{t}_{\mathrm{lag}}\sim \mathrm{1}/\mathrm{4}$, meaning that taking α(t_lag=1 d) should be sufficient for those oscillatory functions.

Motivated by the extra information in multiple t_lag values of the autocovariance function close to the final transition, we plot R(σ,t_lag) in Fig. 6b as it varies with Ro_T. As previously mentioned, R(σ,t_lag) encodes both variance (R(σ,0)) and lag-1 autocorrelation ($R(\mathit{\sigma },\mathrm{1})/R(\mathit{\sigma },\mathrm{0})$) as limiting cases, making R(σ,t_lag) a more complete information source and early warning precursor. R(σ,t_lag) is an even function around t_lag with a maximal value at R(σ,0).

For values of RoT<0.57 and Ro_T>1.2, R(σ,t_lag) is monotonically decreasing with $\left|t_{\mathrm{lag}}\right|$ much like the exponential decay expected for first-order dynamics. Oscillations in R(σ,t_lag) start at Ro_T∼0.57 and remain until Ro_T∼1.2. The start of the oscillations coincides with the start of the plateau in $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$. The end of the oscillatory phase coincides with the time when the full spherical segment begins to superrotate ($\stackrel{\mathrm{‾}}{u}>u_{\mathrm{SR}}$ as shown by the dotted red line). Oscillations reach a maximum in amplitude at Ro_T∼0.87 just before $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}\sim \mathrm{0}$. Throughout the range of their existence, the oscillations have a constant period of approximately 25 d.

4.1.2 First transition: Ro_T∼0.07 at σ∼0.2

The atmosphere first becomes superrotating at the Equator just below the tropopause (σ=0.19) at Ro_T∼0.07. We look at the statistics of δu(σ,t) at σ=0.19 in this subsection. This corresponds to the 14th-lowest vertical level in the model. We look at the last 9 years of a 10-year integration. We have performed longer simulations around the point of the first transition because of the long timescales that are present around the tropopause.

Similar to Fig. 5, in Fig. 7a, $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ is shown near the tropopause (σ=0.19) with varying Ro_T. In Fig. 7b we show the autocovariance function. All time series are very meandering, appearing to be weakly guided random walks when compared to the time series at lower vertical levels.

In Fig. 8, $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ in the equatorial spherical segment against Ro_T is shown. Superrotation at this vertical level occurs approximately for Ro_T>0.08. In Fig. 8b and c we plot the usual early warning indicators, variance of u(σ,t) (upper RH, Var(u)) and lag-1 autocorrelation (lower RH, α(t_lag=1 d)). Both indicators peak just before Ro_T=0.08; however, the results are less convincing than the EWSs of the final transition in the preceding section. As mentioned in Sect. 1, previous work has shown this transition is due to a loss of stability, and therefore it should be detectable by an increase in the EWSs. However, the decay timescales at this vertical height in the atmosphere are long (O(100) d), and thus they may require long time series to detect changes in EWSs reliably. For example, the typical values of α(t_lag=1) corresponded to an e-folding time of O(100) d. Compare this to the e-folding timescales in the lower troposphere, which are an order of magnitude smaller (O(10) d). Changes may also occur over a small range of Ro_T that we have not resolved well with the simulations performed.

Figure 8 Early warning indicators of the first transition at σ=0.19. Each dot is an individual simulation. Coloured dots correspond to the same coloured simulations in Fig. 7. (a) Mean zonal wind speed $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ in a spherical segment centred on the Equator with edges at ±10° latitude during the last 9 years of the simulation against RoT at vertical height σ=0.19. The dotted red line is the threshold for superrotation to occur (u_SR). This is nonzero as the wind speed is averaged over a non-zero thickness latitude band. (b, c) Typical early warning indicators of a local bifurcation: (b) variance of mean zonal u(σ,t) and (c) lag-1 autocorrelation of mean zonal u(σ,t).

4.1.3 Full vertical structure

In Fig. 9, the full vertical structure of the autocovariance function of the zonal mean u in the equatorial band, R(σ,t_lag), is shown. The oscillations of 25 d period found in R(σ,t_lag) that peaked in amplitude around the final transition at σ=0.74 and Ro_T=0.87 are found throughout the lower troposphere (σ>0.3) between values of $\mathrm{0.54}<R{o}_T<\mathrm{1.29}$. These oscillations have their highest amplitude at Ro_T=0.87, σ∼0.45. They have roughly the same period as vertical height, and Ro_T varies throughout their existence. Oscillations are found up to roughly the height of the top of the Hadley cells (σ∼0.3). We have found the oscillations in zonal mean u are also correlated with the strength of the Hadley cells in particular. Larger zonal mean u is simultaneously correlated with a weaker streamfunction in both the Northern Hemisphere and Southern Hemisphere Hadley cells. Conventional Hadley cells with ascent at the Equator are known to flux momentum out of the equatorial region, and this appears to be happening here. The transition to this oscillatory phase of the lower troposphere at Ro_T>0.54 is accompanied by a large increase in decay timescale around the tropopause (σ=0.19) at Ro_T=0.54.

Figure 9 R(σ,t_lag) calculated for δu(σ,t) in a spherical segment centred on the Equator and with edges at ±10° latitude over the last 2 years of each Ro_T simulation. Subplot titles indicate the value of Ro_T.

5 Evidence of multiple states

As previously mentioned, Suarez and Duffy (1992) found an abrupt transition to a superrotating state in a two-layer GCM when the tropical heating was increased above a critical value. This tropical heating term was their control parameter and is different to the control parameter we alter in these experiments (the thermal Rossby number). Not only was the transition Suarez and Duffy (1992) found to be abrupt, unlike the transition in this paper that seems smooth, but the mean state also showed bistability, meaning that two stable states exist for the same value of the control parameter. By slowly changing their control parameter across the transition's critical value from non-superrotating to a superrotating state and then slowly reversing the control parameter back again below the critical value, the superrotating state persisted to a new and lower critical value. This showed that there was a bistable region in which stable superrotating and non-superrotating states can coexist. Unlike Suarez and Duffy (1992) we find no evidence of bistability in the mean equatorial zonal wind speed, although there is evidence of “flickering” (Dakos et al., 2013), meaning a range of Ro_T where there are noise-induced transitions between two stable states.

To investigate whether bistability like that found in Suarez and Duffy (1992) is present, we look at the final transition at vertical height σ=0.74. We slowly increase the value of Ro_T=0.11 through the final transition at Ro_T∼1 and continue up to Ro_T=1.64. We will refer to this part of the simulation as the “ramp up”. Ro_T is then slowly decreased from 1.64 back down to Ro_T=0.11 at the same rate. We will refer to this part of the simulation as the “ramp down”. As in the fixed Ro_T simulations, the magnitude of Ro_T is changed by altering the planetary radius r. Rather than a continuous change in Ro_T, we step Ro_T up and then down by a small increment of 0.11 for every 2 years of simulation due to ease of implementation in the GCM. This gives a total simulation time of 58 years. Because of the step transition every 2 years, there is a small transient effect in the model state before settling back to a stationary state after approximately 200 d. To exclude any spurious effects in the early warning indicators, we calculate them only for the last year (360 d) of each step change in Ro_T.

Figure 10 Ramped control parameter runs with time. (a) Control parameter Ro_T value with simulation year. We start at Ro_T=0.11 and make stepped increases of 0.11 every 2 years to mimic a slow continuous change in the control parameter as closely as possible. We ramp up Ro_T to 1.64 over a period of 30 years (solid black line) before slowly ramping down back to 0.11 over the next 28 years (dotted black line). (c) Zonal mean wind speed u(σ,t) at σ=0.74 in the equatorial spherical segment with time (grey lines) and $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ (black dots are the ramp-up values, black crosses are the ramp-down values) for the last year of each stepped change in Ro_T. (b) Variance of u(σ,t) in the last year of each stepped change in Ro_T. (d) Lag-1 autocorrelation α(t_lag=1) of u(σ,t) in the last year of each stepped change in Ro_T.

Figure 10a shows how Ro_T varies with time. The solid black line indicates the ramp-up of Ro_T, while the dotted black line indicates the ramp-down of Ro_T. In Fig. 10c, u(σ,t) (grey line) in the equatorial band is plotted against time. The periods of increased variance and positive skewness are clearly visible between years 6–18 and 42–50 (also in Fig. 10b). Some of this large-amplitude response is a short transient effect from adjustment to the step change in Ro_T (although this quickly dies away); however, these periods of increased variance and skewness are also present in the static Ro_T simulations. Black dots (ramp-up) and black crosses (ramp-down) in Fig. 10c are $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ in the last year of the step change, avoiding this transient effect. EWS variance and lag-1 autocorrelation are shown in Fig. 10b and d, respectively.

Figure 11 Ramped control parameter runs as a function of Ro_T similar to Fig. 10 but with Ro_T as the x axis variable. (a) Mean zonal wind speed $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ in the equatorial spherical segment for the last year of each stepped change in Ro_T. Black dots show values as Ro_T is ramped up (years 0–30), and black crosses show those values as Ro_T is ramped back down (years 30–58). (b) Variance of u(σ,t) in the last year of each stepped change in Ro_T. (c) Lag-1 autocorrelation α(t_lag=1) of u(σ,t) in the last year of each stepped change in Ro_T. Black dots indicate values for the ramping up of Ro_T, while black crosses show Ro_T as it is ramping back down.

Although $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ looks relatively symmetric around 30 years (the turning point of the control parameter), the EWSs do not. We plot these same quantities as a function of Ro_T in Fig. 11 to assess if there is any hysteresis when the control parameter is reversed and therefore whether bistability can be detected. In Fig. 11a we plot $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ as it is ramping up (black dots) and ramping down (black crosses). Apart from a small region around the plateau in $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ between Ro_T∼0.57–0.76, there does not appear to be any difference in the values of $\stackrel{\mathrm{‾}}{u\left(\mathit{\sigma }\right)}$ when ramping up and ramping down. In Fig. 11b and c EWS differences are larger; however, this is probably due to the short length of the time series used to calculate these statistical estimators and small transient effects. Also, as the standard error in the statistical estimators for a fixed sample size grows as autocorrelation increases, estimation should be worse in regions of large autocorrelation. Such a region exists for Ro_T∈[0.4 1.2], and this is where we find the largest discrepancies. Therefore, we do not expect the small differences to be statistically significant.

Unlike Suarez and Duffy (1992), we find no evidence of bistability in the mean equatorial zonal wind speed. However, our control parameter is markedly different from theirs. This, however, does not completely rule out bistability. There may be smaller regions of Ro_T that these simulations have not been able to resolve.

Although we have not found direct evidence of bistability, there is some evidence for flickering. In a system with multiple basins of attraction and noisy perturbations sufficiently strong to kick the system between these basins, a system may flicker between the attractors (Dakos et al., 2013). Flickering between multiple states can show up as asymmetry in the probability distribution of the system state variable. This asymmetry can be measured by the skewness. As previously mentioned, there are clearly visible periods of increased variance and positive skewness in δu(σ,t) between years 6–18 and 42–50 (Ro_T∼0.5 and Ro_T∼1.2) in Fig. 10c. We plot skewness with Ro_T at vertical height σ=0.74 in Fig. 12c (solid black line). Skewness peaks at Ro_T=0.44 during the ramp-up and Ro_T=0.55 during the ramp-down, although there is high positive skewness (Skew (u)>0.5) for $\mathrm{0.33}<R{o}_T<\mathrm{0.87}$. This peak in skewness occurs just before the oscillations start, although high positive skewness is present throughout the existence of the oscillations. We have also looked at how skewness changes with vertical height. Positive skewness peaks at σ=0.74, which is the level of the final transition. Negative skewness peaks at σ∼0.2 (dotted line in Fig. 12b), which is the vertical height where the first transition occurs. At σ=0.19 a similar graph of skewness with Ro_T is found to that of σ=0.74; however, positive skewness values are mapped to negative values. Little skewness is found around σ∼0.3 (dashed–dotted line); this is roughly the height of the top of the Hadley cells. This may indicate the presence of multiple states in the range $\mathrm{0.33}<R{o}_T<\mathrm{0.87}$. This is also approximately the range where the ramp-up and ramp-down values differ. The flickering states seem to be the oscillating Hadley cell state (typical of the final transition) and the superrotating upper troposphere state (typical of the first transition).

Figure 12 Ramped control parameter R(σ,t_lag) at vertical height σ=0.74 and skewness for δu(σ,t) at vertical heights σ=0.74 (largest positive skewness, final transition vertical height), σ=0.32 (smallest absolute skewness values, height of top of Hadley cells) and σ=0.19 (largest negative skewness values, first transition vertical height). The ramp-up starts from the bottom of the y axis, and the ramp-down starts from Ro_T=1.64 and continues up the y axis. (a) Although there are some small differences in ramp-up and ramp-down values, these are minor, and oscillations in R(σ,t_lag) appear in roughly the same place with the same amplitude and period. (b) High positive skewness (>0.5) and high negative skewness are present in δu(σ,t) for $\mathrm{0.33}<R{o}_T<\mathrm{0.87}$ at vertical heights of σ=0.74 and σ=0.19, respectively. This could indicate the presence of flickering, i.e. noise-induced transitions between attractors.

6 Spatial precursors to the superrotation transition

The oscillations in the fluctuations, δu(σ,t), around the mean zonal u(σ,t) in the equatorial band around the final transition suggest the minimal approximate description of the dynamics must be at a second-order or higher level. Motivated by this and the interesting structure in the variance and autocovariance, we next use a method to diagnose any higher-order dynamics of the fluctuations in the full global horizontal wind field $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)=\left(\mathit{\delta }u\right(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t),\mathit{\delta }v(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t\left)\right)$. We look at the dominant spatial vector modes of δv and their stability at vertical height σ.

In Sect. 4, the state parameter x(t) was a real-valued scalar, and the theory developed to investigate the stability of a stationary state, $\stackrel{\mathrm{‾}}{x}$, to small perturbations, δx(t), assumed first-order dynamics. We remove the restriction on scalar x and therefore make higher-order dynamics possible to describe and diagnose here. This is done by upgrading δx(t) to δx(t), an N component vector (vectors are denoted by bold type face). The technique we next describe goes by several names, such as principal oscillation pattern (POP) analysis (Hasselmann, 1988; von Storch and Zwiers, 1999), linear inverse modelling (LIM) (Penland, 1989) or empirical normal mode (ENM) analysis (Penland and Ghil, 1993; Brunet, 1994). Amongst others, this approach has been used to study the El Niño–Southern Oscillation (Penland and Magorian, 1993) and midlatitude jets (Farrell and Ioannou, 1995). The technique we use here approximates a (typically) non-linear dynamical system with a reduced number of linear modes. However, there are approaches that remove this linear approximation. One such example is the reconstruction of the Koopman operator (e.g. Budišić et al., 2012).

We can repeat the steps in Sect. 4 for the vector generalisation (see Williamson and Lenton, 2015, and Appendix A for details). We again want to study the response of small fluctuations δx(t) around some mean state $\stackrel{\mathrm{‾}}{\mathit{x}}$ to assess the stability of that mean state. As the state becomes less stable approaching the transition, fluctuations become larger and longer lived. Analogously to the scalar case, the dynamics of the fluctuations are determined by the Jacobian J(x); however, this is now a N×N matrix. J(x) can be decomposed into N modes (the eigenvectors of J) associated with N timescales (the eigenvalues of J), the largest of which will dominate the dynamics of the fluctuations. The early warning signals, i.e. autocorrelation and variance, are still functions of J(x) in the vector setting, but they are now also matrices (autocovariance and covariance matrices, respectively). We will analyse the properties of J in the following as this is what ultimately determines stability (and the EWSs, which are functions of J).

We identify the vector δx(t) at each time t in Eq. (A1) to have components δx_i(t) that are the $\mathit{\delta }u(\mathit{\sigma },{\mathit{\theta }}_i,{\mathit{\varphi }}_i,t)$ at each horizontal grid box labelled i in the GCM with coordinates (θ_i,ϕ_i) at vertical level σ. This two-dimensional horizontal field is area-weighted by grid box size, reshaped into a column vector and concatenated with the analogous field of $\mathit{\delta }v(\mathit{\sigma },{\mathit{\theta }}_i,{\mathit{\varphi }}_i,t)$. Time is also a discrete variable in the GCM output and has integer day values with the following interval: Δt=1 d. We label the full vector as δv_t.

As the full-resolution horizontal field in this configuration of Isca has $\mathrm{64}\times \mathrm{128}=\mathrm{8192}$ grid boxes for each field of u and v, the vector has 16 384 components. The resulting reconstructed J would therefore be a real-valued result, but it would be a non-symmetric 16 384×16 384 matrix, and calculating this matrix would be computationally costly. To reduce this cost, we project fields of δu and δv at each time t onto their n largest empirical orthogonal functions (EOFs), the dominant eigenvectors of the covariance matrix Σ, to capture the dominant dynamics of δv(t) (see von Storch and Zwiers, 1999, for a description of the EOF technique). We use n=60, which captures between 76 % and 88 % of the total variance depending on the value of Ro_T and σ. We have experimented with using larger n, and while more of the total variance is captured, the results presented below remain invariant. This is because we are interested in the dominant modes of the dynamics (largest eigenmodes of J), and these are captured well by the projections on to the leading EOFs.

Using the reduced EOF basis for δv_t, we calculate A and Σ (see Appendix A) using the last 2 years of each Ro_T simulation and therefore infer J using Eqs. (A21) and (A22) for the full global horizontal vector wind field fluctuations.

6.1 Final transition: Ro_T∼0.87 at σ=0.74

We now apply the vector technique described in Appendix A to the full spatial horizontal vector wind field $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$ at vertical height σ=0.74 and diagnose the dominant eigenvalues and eigenmodes as Ro_T is varied through the final superrotation transition. Performing an eigendecomposition to find the dominant modes and their stabilities on each reconstructed J results in Fig. 13 for the final transition.

At vertical height σ=0.74, we show that one particular eigenmode, referred to as the “precursor” mode (Fig. 13b and red lines in the middle column of Fig. 13), is responsible for early warnings of the final transition. The precursor mode is an oscillating mode with a zonal wavenumber of 0 and a 25 d period. It becomes important in the dynamics just before the final transition, dominates the dynamics at the final transition and decays quickly following the final transition. We can therefore categorise the precursor mode as a spatial signature heralding the final transition in addition to the usual widely used temporal EWSs. This mode is correlated to oscillations in the equatorial meridional overturning strength (the Hadley cells) seen in Sect. 4.1.3.

After the final transition, the eigenmode referred to as the “superrotation” mode (Fig. 13c and blue lines in the middle column of Fig. 13) dominates the fluctuation dynamics. This is has a zonal wavenumber of 1 and propagates eastward (Fig. 14b) with a period of ∼4 d. The spatial pattern of this mode appears very similar to the Kelvin–Rossby instability that others have previously identified (Mitchell and Vallis, 2010; Wang and Mitchell, 2014; Zurita-Gotor and Held, 2018; Zurita-Gotor et al., 2022). In this mechanism, off-equatorial Rossby waves (diagonal lobes in Figs. 13c and 14b) centred around ±40° latitude couple with an equatorial Kelvin wave (equatorial lobe in the same figures) to converge momentum onto the Equator and enable superrotation.

Figure 13 Spatial precursors of the final superrotation transition at vertical height σ=0.74. The seven leading eigenvalues of the Jacobian of the global horizontal wind field fluctuations (δv(t)) around the mean state $\stackrel{\mathrm{‾}}{\mathit{v}}$ as Ro_T changes. The real part of the eigenvalues, determining the stability and the dominance in the dynamics of the associated mode, are shown in the upper-middle panel as a timescale $\mathit{\tau }=-\mathrm{1}/\mathcal{R}\left({\mathit{\lambda }}_i\right)$. The corresponding imaginary part of the eigenvalues, determining the oscillation frequency of the associated mode, is shown as an oscillation period $\left(P=\frac{\mathrm{2}\mathit{\pi }}{\mathcal{I}\left({\mathit{\lambda }}_i\right)}\right)$. Three of the largest eigenmodes are colour coded and tracked through different values of Ro_T to show how they evolve in the spectrum of J (middle column). In (c) and (d) the eigenmodes associated with the colour-coded eigenvalues are shown. Eigenmodes are dynamic but are shown at t=0 before they oscillate with period P and damping takes place (exponential decay with timescale τ). The filled coloured contours indicate the magnitude of the $\mathit{\delta }u(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$ component of the mode (eigenmodes are dimensionless) with δu<0 (blue) and δu>0 (red) shown. Arrows show the full vector field $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$. The modes are (a) Earth-like (black), (b) precursor (red), (c) superrotation (blue) and (d) swirly (magenta).

6.1.1 Eigenvalues of J: timescales and oscillation frequencies

The spectrum of each J (the first seven eigenvalues, λ_i $i\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\cdots },N\mathit{\}}$, where $\mathcal{R}\left({\mathit{\lambda }}_i\right)\ge \mathcal{R}\left({\mathit{\lambda }}_{i+\mathrm{1}}\right)$) as a function of Ro_T are shown in the middle column of Fig. 13. Generally each λ_i has a real and imaginary part. The real part, ℛ(λ_i), determines how dominant the associated eigenmode is in the dynamics, with less negative values being more dominant. Equivalently, the eigenmodes with the largest e-folding timescales, τ_i, dominate the dynamics (recall ${\mathit{\tau }}_i=-\mathrm{1}/\mathcal{R}\left({\mathit{\lambda }}_i\right)$). We choose to represent the real part of each eigenvalue as a timescale τ_i in the upper-middle panel of Fig. 13.

The imaginary part of each eigenvalue, ω_i=ℐ(λ_i), determines the oscillation frequency of the associated eigenmode. We choose to represent the imaginary part of each eigenvalue as the period of the oscillation $P_i=\frac{\mathrm{2}\mathit{\pi }}{{\mathit{\omega }}_i}$ in the lower-middle panel of Fig. 13.

Each black dot represents one of the seven leading eigenvalues. At low values of Ro_T, the timescales of the leading modes are very close or the same (degenerate) within calculated error bounds from finite-sample-length t_f−t_s. This effective degeneracy in J can result in some non-uniqueness in the decomposed eigenmodes. In contrast, leading τ_i become well separated for Ro_T>0.36 (around the start of the oscillations where skewness peaks) and eigenmodes can be reliably inferred. We have highlighted how the τ_i of three leading eigenmodes around the transition change with Ro_T in the upper-middle panel and the associated oscillation periods P_i in the lower-middle panel Fig. 13. These are colour coded to the associated spatial mode: black is the Earth-like mode (Fig. 13a), red is the precursor mode (Fig. 13b), blue is the superrotation mode (Fig. 13c) and magenta is the “swirly” mode (Fig. 13d).

6.1.2 Eigenvectors of J: spatial mode structure

The leading modes of J, the eigenvectors associated with the eigenvalues with the largest τ_i around the final transition, are discussed next. These are dynamic modes that change with time according to Eq. (A10). We have tried to show the evolution of the precursor and superrotation mode through one full oscillation in Fig. 14 (animations of these modes are also available; see the “Video supplement”). These modes show how the vector wind field adjusts when perturbed from its mean state. An analogy would be the most visible pattern of a drum skin vibration (the dominant eigenmode) from rest (the mean state) when hit with a drumstick (the perturbation). Eigenvectors for a real-valued Jacobian are unique up to a global sign provided eigenvalues are non-degenerate, i.e. $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\equiv -\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i$.

Earth-like mode. This spatial pattern can only be identified for Ro_T=0.02, which is an Earth-like simulation (black, Fig. 13a). Although it is the dominant mode, due to the close spacing of the τ_i at Ro_T=0.02, this pattern alone does not account for very much of the dynamics of δv. It has a timescale of τ∼20 d and an oscillation frequency of P∼20 d. This mode is like a jet stream, capturing the midlatitude instabilities.

Precursor mode. The precursor mode spatial pattern (red, Fig. 13b) appears at Ro_T∼0.5 and is the dominant mode in the fluctuation dynamics between Ro_T≈0.75 and Ro_T=0.87, although a second mode, the superrotation mode, is similar in size during this Ro_T interval. At Ro_T=0.87 the mode peaks in τ with a decay timescale of τ∼25 d, and this is when the early warning indicators of Sect. 4 all peak just before the transition. It is responsible for the shape of the variance of zonal mean u in Fig. 6 and the oscillating autocovariance function R(t_lag) in the same figure as it is a zonal-wavenumber-0 mode with P∼25 d. This is why we refer to this mode as the precursor mode. After the peak, this mode then reduces quickly in importance until it is no longer detectable at Ro_T∼1.3. It is the third-largest mode from Ro_T>1. The evolution of this mode through one full cycle is shown in Fig. 14b. Physically, this appears to represent the fluctuations in the strength of mean meridional overturning circulation, particularly in fluctuations in the strength of the Hadley cells.

Superrotation mode. This mode alone clearly dominates the dynamics (largest τ_i) when the atmosphere becomes fully superrotating for Ro_T≥1 (blue, Fig. 13c). It is also an important component of the fluctuation dynamics around the transition at Ro_T∼0.75 when it is a similar size to the precursor mode. This is also when it first appears as a mode of J. We call this mode the superrotation mode. The e-folding timescale of this mode gets as large as τ∼200 d past Ro_T>1 and has an oscillation frequency of P∼4 d. It is a mode with a zonal wavenumber of 1. The evolution of this mode through one full cycle is shown in the second row of Fig. 14, and this shows that it is propagating eastward. This mode is has a similar pattern to the Kelvin–Rossby instability identified by Mitchell and Vallis (2010) and others for a fully superrotating state, meaning that this is where an equatorial Kelvin wave (the equatorial lobe) couples with off-equatorial Rossby waves (the chevron-like lobes) that mutually amplify each other and converge zonal momentum on to the Equator.

Swirly mode. This mode is significant in the fluctuation dynamics after the superrotation transition (magenta, Fig. 13d), particularly for larger values of Ro_T, although it is never the dominant mode. It is a mode with a zonal wavenumber of 2 and oscillation frequency P∼2 d. We call this mode the swirly mode.

The start of the plateau in $\stackrel{\mathrm{‾}}{u}$ in Fig. 6 coincides with the appearance of the precursor mode, and the plateau ends with the appearance of the superrotation mode. One could argue the plateau in $\stackrel{\mathrm{‾}}{u}$ is also a precursor to the superrotation transition. However, it is the appearance of the precursor mode and its properties that produce the EWSs (rising variance, increasing α(t_lag), R(t_lag)) and heralds the transition to superrotation.

Figure 14 Dynamic evolution of the eigenmodes of J at vertical height σ=0.74. As the eigenmodes of $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$ change with time according to Eq. (A10), we have tried to show this evolution through one full oscillation (please see the Video supplement for more information) of two of the dominant modes in Fig. 13. The filled colour indicates the $\mathit{\delta }u(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$ component of the mode (eigenmodes are dimensionless) with δu<0 (blue) and δu>0 (red) shown. Arrows show the full vector field $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$. The columns show the mode at the start of the cycle (t=0) and advance by a quarter of a cycle left to right. The rows of the figure show the modes: (a) precursor (red) and (b) superrotation (blue).

6.2 First transition: Ro_T∼0.07 at σ∼0.2

We next look at the eigenmodes around the first transition to superrotation at vertical height σ∼0.2. This first transition occurs near the tropopause. Because of the long timescales near the tropopause, we have run 10-year simulations and calculated J for the last 9 years of the simulation to help resolve the distinct modes (Fig. 15).

The first transition to superrotation occurs at vertical height σ∼0.2 at Ro_T∼0.07, and here we find a zonal-wavenumber-1 pattern very similar to the superrotation mode developing prior to the transition and dominating the dynamics. It is a long-period (P∼50–100 d) travelling mode; however, it is propagating westward rather than eastward (Fig. 15a). This mode propagates westward, increasing in period until at the transition the mode becomes effectively a standing mode (P→∞), at least up to the temporal resolution the simulations are able to resolve. It then becomes a long-period travelling-wave pattern again around Ro_T∼0.09 but now propagates eastward like the superrotation mode following the final transition.

Around the point of the first transition, the westward-propagating superrotation mode is replaced by a standing zonal-wavenumber-0 mode as the largest mode of J (Fig. 15b), rapidly dominating the dynamics with relatively large τ (green mode and line in Fig. 15). This remains the dominant mode until at least Ro_T=0.14.

It is harder to resolve modes well at vertical heights above σ<0.22 due to the increased timescales. Even at σ=0.22 and with 9 years of daily data to reconstruct J, only the “standing” mode (green mode) has a sizeable gap in between eigenvalues to resolve the mode well.

Figure 15 Spatial precursors of the first superrotation transition at vertical height σ=0.22. The seven leading eigenvalues of the Jacobian of the global horizontal wind field fluctuations (δv(t)) around the mean state $\stackrel{\mathrm{‾}}{\mathit{v}}$ as Ro_T changes. The real part of the eigenvalues, determining the stability and the dominance in the dynamics of the associated mode, are shown in the left-hand upper panel as a timescale $\mathit{\tau }=-\mathrm{1}/\mathcal{R}\left({\mathit{\lambda }}_i\right)$. The corresponding imaginary part of the eigenvalues, determining the oscillation frequency of the associated mode, is shown as an oscillation period $\left(P=\frac{\mathrm{2}\mathit{\pi }}{\mathcal{I}\left({\mathit{\lambda }}_i\right)}\right)$. The two largest eigenmodes are colour coded and tracked through different values of Ro_T to show how they evolve in the spectrum of J (left column). In the right column we show the eigenmodes associated with the colour-coded eigenvalues. The zonal-wavenumber-1 superrotation eigenmode (blue) (a) is dynamic over a range of Ro_T but is shown at t=0 before it oscillates with period P and damping takes place (exponential decay with timescale τ). The zonal-wavenumber-0 standing mode (green) is a standing wave pattern (P→∞) and dominates the spectrum of J as the vertical level becomes superrotating. The filled coloured contours indicate the magnitude of the $\mathit{\delta }u(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$ component of the mode (eigenmodes are dimensionless) with δu<0 (blue) and δu>0 (red) shown. Arrows show the full vector field $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$. The modes are (a) the superrotation mode (blue) and (b) the standing mode (green).

6.3 Full vertical structure

The leading eigenmode structure follows a typical sequence, like those displayed in Figs. 16 and 17, at vertical heights between the first transition (σ∼0.2) and the final transition (σ=0.74). We have chosen vertical height σ=0.66 to show this typical sequence. This height is above the boundary layer and thus does not experience Rayleigh damping. The transition to superrotation at this vertical height occurs at Ro_T∼0.36.

The sequence of leading eigenmode from low to high Ro_T is as follows. A zonal-wavenumber-1 pattern similar to the superrotation mode (blue in Figs. 16 and 17) develops just before the local transition to superrotation at that particular vertical level. This eigenmode is then replaced by a zonal-wavenumber-0 standing mode pattern (green) at the transition to superrotation at that particular vertical level. As Ro_T is further increased, the precursor mode dominates until the final transition before the zonal-wavenumber-1 superrotation mode finally comes back and dominates the fluctuation dynamics following the final transition.

Generically, the local transition to superrotation at a particular vertical level is preceded by a long-period zonal-wavenumber-1 superrotation mode that is replaced at the transition by a standing zonal-wavenumber-0 mode. The final transition to superrotation, however, is preceded by the precursor mode.

Figure 16 Eigenmodes of J at vertical height σ=0.66. The seven leading eigenvalues of the Jacobian of the global horizontal wind field fluctuations, δv(t), around the mean state $\stackrel{\mathrm{‾}}{\mathit{v}}$ as Ro_T changes. The three largest eigenmodes are colour coded and tracked through different values of Ro_T to show how they evolve in the spectrum of J (middle column). In (c) and (d) we show the eigenmodes associated with the colour-coded eigenvalues. The zonal-wavenumber-0 standing eigenmode (a) dominates the dynamics just after the vertical level σ=0.66 becomes superrotating at Ro_T∼0.36. The precursor mode (b) takes over as the 25 d oscillations become more prevalent (Ro_T>0.5). Following the final transition at Ro_T>0.87, the superrotation mode (c) dominates the dynamics. The filled coloured contours indicate the magnitude of the $\mathit{\delta }u(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$ component of the mode (eigenmodes are dimensionless) with δu<0 (blue) and δu>0 (red) shown. Arrows show the full vector field $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi },t)$.

Figure 17 Leading eigenmodes of J at vertical height σ=0.66 as Ro_T is varied. The e-folding timescale of the mode and period is given in each subtitle in units of days. Subtitles are colour coded by the representative mode. The superrotation-like zonal-wavenumber-1 mode (blue) first appears as the leading eigenmode at Ro_T=0.21. It is an eastward-propagating mode at this Ro_T and vertical level. This is replaced by the standing mode (green), a mode with a zonal wavenumber of 0 and infinite period, as this vertical level becomes superrotating. This is then replaced by another zonal-wavenumber-0 mode, the precursor mode (red), as the leading mode at Ro_T=0.75. This is a previously identified mode with period of P∼25 d. This is then replaced by the superrotation mode (blue) again at Ro_T≥0.97, which is also eastward propagating.

7 Discussion and conclusion

The motivation for this study came from Isaac Held's excellent blog (Held, 2016). One of the posts on that blog introduced the work of Suarez and Duffy (1992), who observed an abrupt transition to a superrotating state in an idealised two-vertical-level GCM as a longitude-dependent heat source in the tropics was increased. The atmospheric state in this model also showed bistability, i.e. both superrotating and regular, and present-day atmospheric states were stable for the same value of tropical heating. This example has many similarities with other climate tipping points (Lenton et al., 2019) – abrupt changes in state caused by relatively small changes in forcing parameters. These are generally regarded as low-probability, high-risk events. One such example is the Atlantic meridional overturning circulation (AMOC), which can abruptly change from an on state to an off state as freshwater flux into the North Atlantic is increased (Stommel, 1961; Vellinga and Wood, 2002; Mecking et al., 2016). The AMOC can also show bistability. The off state can persist as the freshwater flux is decreased back below the original threshold. There are a variety of other widely studied climate tipping points (e.g. the Amazon rainforest, Greenland ice sheet and Arctic sea ice); however, superrotation has not received the same attention. One of the motivations of this work was to begin to correct this deficiency.

Unlike Suarez and Duffy (1992), where the control parameter varied (the tropical heating) is arguably relevant to future anthropogenic climate change, our control parameter, the thermal Rossby number (Ro_T), a non-dimensional number that controls the global atmospheric dynamics, is probably not relevant to future climate change. Unlike the abrupt transition in equatorial mean zonal $\stackrel{\mathrm{‾}}{u}$ observed in Suarez and Duffy (1992), our transition is smooth as the control parameter is varied. The third difference is that we have not been able to detect bistability (Sect. 5), although evidence of noise-induced transitions between multiple states (flickering) was found within a subset of the control parameter range.

Although the thermal Rossby number is not likely to change very much as the climate warms, the mechanism for the transition may well be relevant, namely the Kelvin–Rossby instability. By changing the latitude of the baroclinic instability forcing, Zurita-Gotor et al. (2022) found that an abrupt transition to superrotation could be achieved via the Kelvin–Rossby instability mechanism. They also postulated this mechanism might be relevant for the transition found in warm climates by Caballero and Huber (2010). The simulations reported here do show early warnings of this transition and therefore serve as a test bed for the second motivation of this paper, namely whether we can detect this transition before it happens.

In the simulations reported here, the first transition to superrotation occurs just below the tropopause at vertical height σ∼0.2 and Ro_T∼0.07. As Ro_T increases, this superrotating region expands horizontally into the tropics, vertically up into the stratosphere and downwards into the troposphere. The superrotating region finally terminates at σ=0.83 at Ro_T=0.87 and does not extend any further, at least for values up to Ro_T=11.4. We term the transition of the lowest vertical level to superrotation at σ=0.74 as the final transition. We therefore do not see surface superrotation in these simulations, although this is known to occur in small regions of parameter space (Caballero and Carlson, 2018).

The usual early warnings of local bifurcation type tipping (rising variance and autocorrelation, Sect. 4) all increase before the both the first and final superrotation transition, indicative of critical slowing down near a local bifurcation. These EWSs are harder to detect near the tropopause around the first transition due to the difficulty resolving the differences between the long timescales that are naturally present in this part of the atmosphere.

In addition, the full autocovariance function, R(t_lag), showed oscillations of a period of P∼25 d that became increasingly less damped as the final transition was approached, reminiscent of a Hopf type bifurcation (see Eq. 16). In a Hopf bifurcation, a pair of complex conjugate eigenvalues of the Jacobian move in the complex plane from the negative real half and cross the imaginary axis at the point of bifurcation (see Strogatz, 2001, for details). Before the bifurcation, this results in critical slowing down of the dynamics (increasingly less negative ℛ(λ)) and oscillations of frequency ω=ℐ(λ). These are observed in the superrotation transition. These oscillations are not present all the way through the troposphere. They exist for σ>0.3 and peak at σ∼0.45. These oscillations were correlated with oscillations in the strength of the Hadley cells. Both Northern Hemisphere and Southern Hemisphere cells simultaneously increase and decrease in strength as δu(σ,t) decreases and increases in strength, respectively.

In Thompson et al. (1994), the authors classify all codimension-1 bifurcations via observations from a slow control parameter sweep in a series of tables. Using this classification, our observations of a smooth transition, a Hopf-like autocorrelation function and lack of hysteresis point to a local supercritical Hopf bifurcation being present at the transition when the control parameter is the thermal Rossby number.

The presence of oscillations in δu motivated us to use a more general technique to study the transitions. The usual theory of EWSs relies on scalar state variables and first-order dynamics. This does not permit oscillations or higher-order dynamics. We gave a pedagogical introduction to the more general theory based on POP analysis in Appendix A and calculated the Jacobian and normal modes for the full global vector wind field $\mathit{\delta }\mathit{v}(\mathit{\sigma },\mathit{\theta },\mathit{\varphi })$. EWSs are designed to detect changes in the stability of the mean state and therefore must be functions of the Jacobian's eigenvalues that entirely describe the linear stability of the mean state. We therefore chose to study the eigenvalues and eigenmodes of the Jacobian as they are the sole determinant of traditional EWSs based on linear stability analysis.

When a system approaches a transition resulting from a local bifurcation, often a single eigenmode of the Jacobian starts to dominate the dynamics. Close to the final transition, we found two eigenmodes (or spatial patterns) that the fluctuations, δv, preferentially oscillate in before slowly decaying back towards the mean state $\stackrel{\mathrm{‾}}{\mathit{v}}$. One of these modes, the precursor mode, was a zonal-wavenumber-0 mode that we attribute to the 25 d oscillations in the equatorial zonal mean δu and the rise in variance and autocorrelation seen in Sect. 4. There is also a second mode of roughly the same size around the final transition – the superrotation mode – a zonal-wavenumber-1 mode that comes to dominate the dynamics as the majority of the troposphere became fully superrotating. This is an eastward-propagating mode of a roughly 4 d period that looks very similar to the coupled Kelvin–Rossby mode identified by others.

Detecting the signature of the precursor mode in the spatial or temporal correlations would serve as the best EWS of the final transition. These modes were calculated for the global wind field. It is plausible that the precursor mode would be even more significant in the dynamics of the tropics if the spatial analysis was performed over just the tropical region rather than the entire globe. The attribution of a spatial signature as an early warning signal is something that has not been reported before, as far as we are aware, for a climate tipping point.

Just before the first transition, a zonal-wavenumber-1 mode very much like the superrotation mode appears, although it is propagating very slowly westward rather than eastward. At the first transition, this is replaced by a standing zonal-wavenumber-0 mode. As Ro_T increases, at vertical heights σ>0.3 this mode is replaced by the oscillating precursor mode, which peaks at the final transition at Ro_T=0.87 before the superrotation mode returns and dominates for Ro_T>0.87. Generically it appears the local transition to superrotation at a particular vertical level is preceded by a long-period zonal-wavenumber-1 superrotation mode that is replaced at the transition by a standing zonal-wavenumber-0 mode.

One may wonder why oscillations are not observed for Ro_T>1.3 after the precursor mode disappears from the dynamics in the scalar EWSs even though the dominant modes in this regime, the superrotation mode and the swirly mode, are oscillatory with periods of P∼4 and P∼2 d, respectively. This is due to their zonal wavenumbers. Remember that equatorial zonal u was averaged over a latitude band in Sect. 4. For non-zero zonal wavenumbers, the average will be always be zero, as there are an equal number of peaks and troughs. A zonal-wavenumber-0 pattern like the precursor mode. however, will have a non-zero average in a zonal mean.

There are still questions about the dynamical mechanisms giving rise to the observed eigenmodes around the transition. There are a lot of interesting dynamics going on. The superrotation mode appears to be the Kelvin–Rossby instability previously identified and studied by others (i.e. Williams, 2003; Wang and Mitchell, 2014; Mitchell and Vallis, 2010; Zurita-Gotor and Held, 2018; Zurita-Gotor et al., 2022). In this mechanism, superrotation arises from the wave–mean-flow interaction, specifically the interaction of equatorial Kelvin and off-equatorial Rossby waves that mutually amplify each other and generate a planetary-scale mode that converges zonal momentum onto the Equator. The global 25 d oscillations present around the transition are correlated to the strength of the Hadley cell intensity and are present up until the top of the meridional overturning cells. The banded zonal-wavenumber-0 structure also seems to hint that this mode represents the meridional overturning dynamics. It is also known that conventional Hadley cells with uplift at the Equator flux momentum out of this equatorial region and therefore should reduce the magnitude of superrotation based on how strong these cells are (which is observed here). These are hints at the underlying dynamics and details that have been left to a future study.

The degree of spatial correlation in the fluctuations, particularly in the precursor mode, is reminiscent of a phase transition (Anderson, 1984). Traditionally, phase transitions are studied via the statistics of the spatial fluctuations around the mean state. The hallmarks of a phase transition are diverging correlation length and correlation time (resulting in critical slowing down) and an abrupt change in the order parameter of the state. In particular, continuous phase transitions have spatial correlation functions that are described by a power law at the critical point giving rise to scale-invariant physics and often involve a broken symmetry. What are the similarities with the transition to superrotation? Diverging correlation time is present (τ), and diverging correlation length is qualitatively present (as seen in the pattern of the precursor mode, at least zonally), although we have not formally calculated spatial correlation functions. However, Dylewsky has found quantitative evidence via a neural network trained on the spatial information from a 2D Ising model phase transition. This trained neural network clearly predicts the superrotation transition before it occurs (Daniel Dylewsky, personal communication, 2023; see also Dylewsky et al., 2023).

Although the probability of transitioning to an equatorial surface superrotating atmospheric state under future climate change seems unlikely (Caballero and Carlson, 2018), it is not zero (see Held, 2016). There is (limited) evidence that it has happened in the past and can happen in the future (Huang et al., 2001; Caballero and Huber, 2010; Tziperman and Farrell, 2009). As with other tipping points, this makes superrotation a low-probability but high-risk event. Even though a surface superrotating atmospheric state seems unlikely, superrotation happens readily in the upper troposphere under high-end global warming scenarios. However, the impact on surface climate would probably be less severe (Caballero and Carlson, 2018). It would be interesting to repeat the experiment of Suarez and Duffy (1992) with the Isca framework, which would show a more relevant way of transitioning to superrotation under climate change, and see if the EWSs identified here are still present.

Appendix A: Vector-valued state parameter generalisation for early warnings of bifurcation-induced tipping

We can repeat the steps in Sect. 4 for the vector generalisation (see Williamson and Lenton, 2015, for further details). We again want to study the response of small fluctuations δx(t) around some mean state $\stackrel{\mathrm{‾}}{\mathit{x}}$ to assess the stability of that mean state. As the state becomes less stable approaching the transition, fluctuations become larger and longer lived. The state of the system in analogy to Eq. (9) is

$$\begin{array}{}\text{(A1)}& \mathit{x}\left(t\right)=\stackrel{\mathrm{‾}}{\mathit{x}}+\mathit{\delta }\mathit{x}\left(t\right).\end{array}$$

The N components of the vector x(t) can be written explicitly as

$$\begin{array}{}\text{(A2)}& \mathit{x}\left(t\right)=\left(x_{\mathrm{1}}\right(t),x_{\mathrm{2}}(t),\mathrm{\dots },x_{N-\mathrm{1}}(t),x_N(t){)}^T.\end{array}$$

Similar expressions can be written for δx and $\stackrel{\mathrm{‾}}{\mathit{x}}$. The superscript T denotes the matrix transpose operator. The individual scalar components of x(t) are indexed with i, which can take the values $i\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N-\mathrm{1},N\mathit{\}}$. Each of the x_i values is assumed to evolve autonomously, each according to a different and generally non-linear function f_i(x). Note that each f_i has the multi-variable input x. This means that f_i is a function of all $\mathit{\{}{x}_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_N\mathit{\}}$ in general, meaning that each x_i value is coupled, possibly non-linearly. Formally, each x_i evolves according to $\dot{x_i}=f_i\left(\mathit{x}\right)$. Provided the fluctuations are not too large, the dynamics are approximated well by making a multi-variable Taylor expansion of each f_i(x) to the first order around the mean state, i.e. $f_i\left(\mathit{x}\right)\approx f_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)+{\sum }_{j=\mathrm{1}}^{N}J(\stackrel{\mathrm{‾}}{\mathit{x}}{)}_{i,j}\mathit{\delta }{x}_j$, where instead of a single number describing the stability of the system, the Jacobian is now an N×N matrix encoding the linear coupling between the components of x. The element in the ith row and jth column of J is given by $J(\stackrel{\mathrm{‾}}{\mathit{x}}{)}_{i,j}=\frac{\partial f_i\left(\mathit{x}\right)}{\partial x_j}{|}_{\mathit{x}=\stackrel{\mathrm{‾}}{\mathit{x}}}$. The dynamics of the fluctuations can be now be approximated with a linear first-order ODE; however, it is now a matrix equation

$$\begin{array}{}\text{(A3)}& \dot{\mathit{\delta }\mathit{x}}\approx J\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)\mathit{\delta }\mathit{x}.\end{array}$$

Analogous to the scalar x case, for a time-invariant mean state $\stackrel{\mathrm{‾}}{\mathit{x}}$ (and time-invariant J(x)), Eq. (A3) has the solution

$$\begin{array}{}\text{(A4)}& \mathit{\delta }\mathit{x}\left(t\right)=e^{J\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)t}\mathit{\delta }\mathit{x}\left(\mathrm{0}\right).\end{array}$$

J(x) still determines the behaviour and fate of small perturbations from $\stackrel{\mathrm{‾}}{\mathit{x}}$ and therefore $\stackrel{\mathrm{‾}}{\mathit{x}}$'s stability. However, J(x) is a more complicated and general object than the single number of Sect. 4 due to the couplings between the various δx_i (non-zero off-diagonal entries of J). One can uncouple the dynamics into N uncoupled modes, $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i$, by transforming each of the δx_i values into a new basis (denoted by a tilde) that is a linear combination of the old δx_i, i.e.

$$\begin{array}{}\text{(A5)}& \mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)=\sum _{j=\mathrm{1}}^N{a}_j^{\left(i\right)}\mathit{\delta }{x}_j\left(t\right),\end{array}$$

where $a_j^{\left(i\right)}$ represents scalar (and possibly complex) coefficients. These new $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i$ values are vectors, although they behave like scalars and evolve independently of each other. Therefore, one can treat the dynamics as N scalar (although generally complex, rather than real valued) first-order ODEs like the one in Eq. (10). The problem of finding this new uncoupled basis, where the off-diagonal entries of J are zero, is equivalent to finding the eigendecomposition of the Jacobian

$$\begin{array}{}\text{(A6)}& J=Q\stackrel{\mathrm{̃}}{J}{Q}^{-\mathrm{1}},\end{array}$$

where $\stackrel{\mathrm{̃}}{J}$ is a diagonal matrix (only non-zero entries are on the diagonal) with N eigenvalues, each labelled by λ_i on the diagonal of $\stackrel{\mathrm{̃}}{J}$ associated with the N eigenvectors of J, denoted q_j, given by the jth column of matrix Q. The same Q also diagonalises the matrix $e^{J\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)t}$. Thus, the new basis in Eq. (A5) is simply given by $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i={\sum }_{j=\mathrm{1}}^N(Q^{-\mathrm{1}}{)}_{i,j}\mathit{\delta }{x}_j$, where $(Q^{-\mathrm{1}}{)}_{i,j}\equiv a_j^{\left(i\right)}$ gives the coefficients to uncouple the equations. The new basis given by the set of $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i$ is the eigenmode, but these also go by other names; i.e. they are often known as linear, critical or normal modes.

How each of the $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i$ evolve now reduces to the N equation

$$\begin{array}{}\text{(A7)}& \mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)=e^{{\mathit{\lambda }}_{i}t}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(\mathrm{0}\right),\end{array}$$

and therefore again there are three possible fates of an initial perturbation of amplitude $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(\mathrm{0}\right)$ on mode i depending on the sign of the real part of each of the λ_i. (i) The real part of λ_i is negative (ℛ(λ_i)<0), where the initial fluctuation in the normal mode $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(\mathrm{0}\right)$ decays, and this particular mode will be stable to perturbations. In this case a recovery time for this mode can be defined as ${\mathit{\tau }}_i=-\mathrm{1}/\mathcal{R}\left({\mathit{\lambda }}_i\right)$. (ii) The real part of λ_i is positive (ℛ(λ_i)>0), where any initial fluctuation in this mode will get exponentially larger with time, and this mode will be unstable. (iii) The real part of λ_i is zero, i.e. ℛ(λ_i)→0. In this case, perturbations will not grow or decay in this mode. This neutral stability occurs at a bifurcation.

The evolution of the fluctuations around the mean of the entire system is described by a linear superposition of all N modes, i.e.

$$\begin{array}{}\text{(A8)}& \mathit{\delta }\mathit{x}\left(t\right)=\sum _{i=\mathrm{1}}^N{e}^{{\mathit{\lambda }}_{i}t}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(\mathrm{0}\right).\end{array}$$

Therefore, all N modes must be stable for the overall system to be stable. If one or more of the modes are unstable, as time increases the unstable mode with the most positive real part will swamp the response of the other N−1 modes and cause the entire system to be unstable. This means that for a stable system the requirement must be ℛ(λ_i)≤0 ∀i. Provided all the modes are stable, the $\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i$ with the least negative ℛ(λ_i) will dominate the dynamics of the system simply because they decay more slowly (larger τ_i) than those modes associated with shorter τ_i (more negative ℛ(λ_i)). This means that the dynamics of the entire system may be described well by just a few (n) modes, where n<N, with the largest timescales such that

$$\begin{array}{}\text{(A9)}& \mathit{\delta }\mathit{x}\left(t\right)\approx \sum _{i=\mathrm{1}}^n{e}^{{\mathit{\lambda }}_{i}t}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(\mathrm{0}\right),\end{array}$$

where λ_i values are ordered in increasing size, i.e. $\mathcal{R}\left({\mathit{\lambda }}_i\right)\ge \mathcal{R}\left({\mathit{\lambda }}_{i+\mathrm{1}}\right)$. An extreme example is when a local bifurcation is approached and the timescale of one (or more) of the modes is τ_i→∞ (ℛ(λ_i)→0) and dominates the dynamics to such a degree that the system can be described by just that one (or more) mode.

The eigenvalues λ_i can also have imaginary components. These give rise to oscillations in the associated mode with angular frequency ω_i=ℐ(λ_i). For a real-valued J, complex λ_i values always occur in complex conjugate pairs with associated complex conjugate eigenvectors. If eigenvalues (and therefore eigenvectors) are complex, the superposition with its conjugate mode gives a real but time-evolving spatial mode. Explicitly, for a single dominant complex (oscillatory) mode (ℛ(λ₁)≫ℛ(λ_i), $\forall i\ne \mathrm{1}$), the dynamics are

$$\begin{array}{}\text{(A10)}& {\displaystyle }\mathit{\delta }\mathit{x}\left(t\right)\sim e^{{\mathit{\lambda }}_{\mathrm{1}}t}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}}+e^{{\mathit{\lambda }}_{\mathrm{1}}^{*}t}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}}^{*},\text{(A11)}& {\displaystyle }\sim e^{-\frac{t}{{\mathit{\tau }}_{\mathrm{1}}}}\left[\cos \left({\mathit{\omega }}_{\mathrm{1}}t\right)\mathcal{R}\left(\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}}\right)-\sin \left({\mathit{\omega }}_{\mathrm{1}}t\right)\mathcal{I}\left(\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}}\right)\right],\text{(A12)}& {\displaystyle }\sim e^{-\frac{t}{{\mathit{\tau }}_{\mathrm{1}}}}\cos ({\mathit{\omega }}_{\mathrm{1}}t+{\mathit{\varphi }}_{\mathrm{1}})\left|\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}}\right|.\end{array}$$

This means that the mode spatially oscillates between its real and imaginary parts at angular frequency ω₁. In the last line, we define the following vector phase: ${\mathit{\varphi }}_{\mathrm{1}}=\arg \left(\mathit{\delta }{\stackrel{\mathrm{̃}}{x}}_{\mathrm{1}}\right)$.

For a dominant real mode, the dynamics are

$$\begin{array}{}\text{(A13)}& {\displaystyle }\mathit{\delta }\mathit{x}\left(t\right)\sim e^{{\mathit{\lambda }}_{\mathrm{1}}t}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}},\text{(A14)}& {\displaystyle }\sim e^{-\frac{t}{{\mathit{\tau }}_{\mathrm{1}}}}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}}.\end{array}$$

Since ω₁=0 and $\mathcal{I}\left(\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}}\right)=\mathrm{0}$, there is no conjugate mode. A mode associated with a real eigenvalue has real eigenmode, and therefore the spatial pattern remains fixed, decaying exponentially and uniformly from the constant multiplier on the pattern $e^{-\frac{t}{{\mathit{\tau }}_{\mathrm{1}}}}$. This case is analogous with the first-order dynamics scalar case of Sect. 4, with the constant multiplier being the lag-1 autocorrelation and the spatial mode analogous to the scalar variable x.

One may wonder how Eq. (A3) describes higher-order dynamics since it is a first-order ODE. A Nth-order system can be written as a first-order vector ODE with N dimensions. For example, the third-order system

$$\begin{array}{}\text{(A15)}& \stackrel{\dots }{x}+a\ddot{x}+b\dot{x}+cx=\mathrm{0},\end{array}$$

can be written as a first-order vector ODE in three dimensions by making the identifications $y=\dot{x}$ and $z=\dot{y}$ so that Eq. (A15) can be written

$$\begin{array}{}\text{(A16)}& {\displaystyle }\dot{x}=y,\text{(A17)}& {\displaystyle }\dot{y}=z,\text{(A18)}& {\displaystyle }\dot{z}=-cx-by-az,\end{array}$$

or as a first-order vector ODE

$$\begin{array}{}\text{(A19)}& {\displaystyle }\left(\begin{array}{c}\dot{x}\\ \dot{y}\\ \dot{z}\end{array}\right){\displaystyle }=\left(\begin{array}{ccc}\mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\\ -c& -b& -a\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right),\end{array}$$

which can be written more compactly in the same form as Eq. (A3) as follows:

$$\begin{array}{}\text{(A20)}& \dot{\mathit{x}}=J\mathit{x}.\end{array}$$

Under the assumption of white-noise forcing (or equivalently independent errors in δx_i), one can reconstruct J from time series of the δx_i, with each time series sampled at an interval Δt between time steps as follows (see Williamson and Lenton, 2015, for details)

$$\begin{array}{}\text{(A21)}& e^{J\left(\mathit{x}\right)\mathrm{\Delta }t}=A{\mathrm{\Sigma }}^{\mathrm{-}\mathrm{1}},\end{array}$$

where Σ is the covariance matrix of δx and A is the lagged covariance matrix:

$$\begin{array}{}\text{(A22)}& {\displaystyle }\mathrm{\Sigma }=E\left[\mathit{\delta }{\mathit{x}}_t\mathit{\delta }{\mathit{x}}_t^T\right],\text{(A23)}& {\displaystyle }A=E\left[\mathit{\delta }{\mathit{x}}_{t+\mathrm{\Delta }t}\mathit{\delta }{\mathit{x}}_t^T\right].\end{array}$$

One can then find the eigenmodes and eigenvalues of J (and hence their stability). For vector-valued x, the analogous EWSs to the scalar case (Sect. 4) are now matrix valued:

$$\begin{array}{}\text{(A24)}& {\displaystyle }A{\mathrm{\Sigma }}^{-\mathrm{1}}\equiv \mathit{\alpha }\left(t_{\mathrm{lag}}\right),\text{(A25)}& {\displaystyle }\mathrm{\Sigma }\equiv \text{Var}\left(x\right).\end{array}$$

However, we now have N modes (eigenvectors of J) associated with N timescales (eigenvalues of J). The autocorrelation and variance of early warning signals are still functions of J(x), but they are now also matrices (autocovariance and covariance matrices, respectively). We analyse the properties of J in this paper as this is what ultimately determines stability (and the EWSs, which are functions of J).