Issue 
A&A
Volume 633, January 2020



Article Number  A2  
Number of page(s)  26  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201936110  
Published online  19 December 2019 
Acceleration of superrotation in simulated hot Jupiter atmospheres
^{1}
IRAP, Université de Toulouse, CNRS, UPS, Toulouse, France
email: florian_debras@hotmail.com
^{2}
Ecole normale supérieure de Lyon, CRAL, UMR CNRS 5574,
69364 Lyon Cedex 07, France
^{3}
Physics and Astronomy, College of Engineering, Mathematics and Physical Sciences, University of Exeter,
Exeter, EX4 4QL, UK
^{4}
Centre for Fusion, Space & Astrophysics, Department of Physics, University of Warwick,
Coventry, CV4 7AL, UK
^{5}
Mathematics, College of Engineering, Mathematics and Physical Sciences, University of Exeter,
Exeter, EX4 4QL, UK
Received:
14
June
2019
Accepted:
7
November
2019
Context. Atmospheric superrotating flows at the equator are a nearly ubiquitous result when conducting simulations of hot Jupiters. One theory explaining how this zonallycoherent flow reaches equilibrium has already been developed in the literature. This understanding, however, relies on the existence of either an initial superrotating flow or a sheared flow, coupled with a slow evolution that permits a linear steady state to be reached.
Aims. A consistent physical understanding of superrotation is needed for arbitrary drag and radiative timescales, along with the relevance of taking linear steady states into account, needs to be assessed.
Methods. We obtained an analytical expression for the structure, frequency, and decay rate of propagating waves in hot Jupiter atmospheres around a state at rest in the 2D shallowwater βplane limit. We solved this expression numerically and confirmed the robustness of our results with a 3D linear wave algorithm. We then compared it with 3D simulations of hot Jupiter atmospheres and studied the nonlinear momentum fluxes.
Results. We show that under strong daynight heating, the dynamics do not transit through a linear steady state when starting from an initial atmosphere in solid body rotation. We further demonstrate that nonlinear effects favor the initial spinup of superrotation and that acceleration due to the vertical component of the eddymomentum flux is critical to the initial development of superrotation.
Conclusions. We describe the initial phases of the acceleration of superrotation, including the consideration of differing radiative and drag timescales, and we conclude that eddymomentumdriven superrotating equatorial jets are robust, physical phenomena in simulations of hot Jupiter atmospheres.
Key words: planets and satellites: gaseous planets / planets and satellites: atmospheres / hydrodynamics / waves / methods: analytical / methods: numerical
© F. Debras et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Understanding the atmospheric dynamics of hot Jupiters, which are Jovian planets in shortperiod orbits, has been a major challenge since their discovery (Mayor & Queloz 1995). Due to their proximity to the host star, hot Jupiters are expected to be tidallylocked (for a review, see Baraffe et al. 2010), resulting in a permanent dayandnightside driving atmospheric circulations with no equivalent in our solar system. This, in turn, is likely to lead to the mixing of material between the two hemispheres (Drummond et al. 2018a,b).
Cooper & Showman (2005) performed the first study of the atmosphere of HD 209458b (e.g., Charbonneau et al. 2002; Sing et al. 2008; Snellen et al. 2008) using a general circulation model (GCM). Such GCMs have been used extensively to characterize hot Jupiters (e.g., Showman et al. 2008; Heng et al. 2011; Rauscher & Menou 2012; DobbsDixon & Agol 2013; Mayne et al. 2014a; Helling et al. 2016). The physical complexity, or completeness, of these GCMs varies greatly; for example, treatments of the dynamics and radiative transfer range from those adopting the primitive equations of dynamics, along with simple Newtonian cooling, to those solving the full Navier–Stokes equations and more accurate radiative transfer (see, notably, Amundsen et al. 2014, 2016). Recent advances have also been made in the treatment of chemistry, regarding both the gas phase (see Drummond et al. 2016, 2018c,a,b; Tsai et al. 2018) and the condensates, or clouds (Lee et al. 2016; Lines et al. 2018a,b; Roman & Rauscher 2019).
There is a common feature that has emerged from almost all GCM studies of hot Jupiters: the atmosphere exhibits equatorial superrotation, a prograde atmospheric wind velocity, greater than that which arises from the rotation of the planet alone, over a range of pressures. Observations have detected an eastward shift of the peak infrared flux from the substellar point in the atmosphere of hot Jupiters (Knutson et al. 2007; Zellem et al. 2014), which is consistent with that found in simulations caused by the advection of heat by the superrotating jet. Mayne et al. (2017) attempted to suppress the formation of the equatorial jet in simulated hot Jupiter atmospheres by forcing the deep atmospheric flow or by altering the model parameters. They found superrotation to be a very robust feature in numerical simulations. However, a recent measurement has inferred an opposite, westward shift for CoRoT2b (Dang et al. 2018). Armstrong et al. (2016) previously obtained variability in the position of the hot spot over time, suggesting an additional layer of complexity (a potential link to magnetic fields has recently been investigated by Hindle et al. 2019). Superrotation, therefore, must be explained on the basis of sound physical arguments.
Showman & Polvani (2011) were the first to study the formation of a superrotating jet in simulated hot Jupiter atmospheres using a simplified twolayer model. In exploring the linear steady state of the atmosphere, Showman & Polvani (2011) highlight the formation of a Matsuno–Gill (hereafter MG, see Matsuno 1966; Gill 1980) pattern, whereby the atmospheric perturbations are “tilted” in the latitudelongitude plane driving momentum transport tothe equator and accelerating the jet. Showman & Polvani (2011) posit that the nonlinear equilibrium is reached when the transport of meridional and vertical eddy momentum into the region, acting to accelerate and decelerate the jet, respectively, are balanced by the atmospheric drag. Tsai et al. (2014) extend the study to a full 3D dynamical model and include a consideration of the resonance of the atmospheric wave response, as well as the “tilt” of the vertical component which acts to drive the vertical eddymomentum transport, under the assumption of equal drag and radiative timescales. This is followed by Hammond & Pierrehumbert (2018), who explore superrotation in 2D with the addition of a shearing flow. PerezBecker & Showman (2013) consider the propagation of waves and the resulting balance for the equilibrated jet and propose diagnostics for predicting the daytonighttemperature contrast, controlled by the efficiency of the zonal advection. This analysis is improved upon by Komacek & Showman (2016) across an extensive range of dissipation timescales.
There is, however, an inherent discrepancy between the works of Komacek & Showman (2016) and Showman & Polvani (2011): when the atmospheric drag timescale is large, superior to a few 10^{5} s, the linear steady states obtained in Komacek & Showman (2016) tend to decelerate the equator although the associated nonlinear steady states exhibit equatorial superrotation. This raises the question of whether superrotation can be properly explained through the transition from a linear steady state.
Specifically, the study of Tsai et al. (2014) is only valid through the moderate to strong dissipation limit. The study of Hammond & Pierrehumbert (2018) requires an initial sheared superrotation. However, Komacek & Showman (2016) show that superrotation develops only if the dissipation is sufficiently low (see their Fig. 4). Current theories are, therefore, applicable only once an initial flow has been set up and its evolution is slow compared to the wave propagation time.
In this study, we address the issue of what drives the initial spinup of superrotation in simulated hot Jupiter atmospheres. In order to do this we develop a description of the timedependent waves supported by our simulations of hotJupiter atmospheres with arbitrary drag and radiative timescales and we determine which are responsible for driving the evolution of the jet. Firstly, in Sect. 2, we state our main assumptions and develop the mathematical framework we have adopted throughout this work, before finding the form of the timedependent linear solution to the betaplane equations. Additionally, in this section we summarize the main results of Showman & Polvani (2011), Tsai et al. (2014) and Komacek & Showman (2016), which we have based our study upon. Obtaining the form of the solution to the timedependent case is not sufficient as the controlling parameters remain unconstrained and are not easily accessible analytically. Therefore, in Sect. 3 we numerically explore the sensitivity of the steady linear solution to the shape of the forcing, or heating, showing that the linear steady state requires a daynight heating contrast but that it is insensitive to the exact shape of the forcing itself. This confirms that the limitations of the current theory do not come from the simplified form of the forcing but the timedependent linear considerations must be included. Therefore, we conduct a numerical study of the propagating waves, apart from the special case of Kelvin waves (which have an analytical expression) to build a more complete picture of the physical process of the acceleration of superrotation. In Sect. 4 we determine the characteristic decay timescale for different waves and arbitrary dissipative timescales which are used to explain the structure of both the linear steady states presented in Fig. 5 of Komacek & Showman (2016) and the timedependent linear evolution of simulated hot Jupiters. In Sect. 5, we combine the understanding developed throughout this study to detail the transition to superrotation in 3D GCM simulations through eddymean flow interaction under different conditions, revealing the importance of timedependent linear considerations as well as vertical momentum transport across different drag and forcing regimes. Finally, we summarize our conclusions in Sect. 6. Overall, our study shows that the paradigm of equatorial superrotation in hot Jupiters is robust: superrotation is accelerated by an eddymean flow interaction (i.e., atmospheric waves interacting with the background flow) and it is strongly influenced by the wavedissipation timescales and vertical momentum convergence.
2 Solution to 2D shallow water equations
2.1 Theoretical framework
For this study, we adopted the 2D shallow water equations under the equatorial βplane approximation. As in Showman & Polvani (2011), the shallow water approximation consists of considering that the planet can be decomposed into an upper, constant density but dynamically active layer with a free surface at the bottom exchanging heat and momentum with a lower, quiescent layer of a much higher density. The equatorial βplane simplifies the view of the spherical planet as a local Cartesian plane at the equator, with the Coriolis parameter depending linearly on the Cartesian meridional coordinate, y, with a factor of proportionality β = 2Ω∕R, where Ω is the rotation rate of the planet and R its radius. The further away from the equator, the less valid this approximation, but it allows for analytic solutions, making it particularly suited for the study of equatorial superrotation. Wu et al. (2000) show that the 3D structure of solutions to the βplane system can be decomposed onto an infinite sum of solutions of 2D βplanes with different characteristic heights. The importance of this decomposition regarding hot Jupiter atmospheric dynamics is emphasized by Tsai et al. (2014), where a vertical shift of the wave response is presented when the mean background velocity is changed (their Fig. 10). We begin bysummarising the main results of Matsuno (1966), Gill (1980), and Showman & Polvani (2011), all of which solve the 2D βplane equations.
Following Showman & Polvani (2011), the nondimensional, linearized equations of motion for a forced 2D, equatorial βplane can be written as $$\begin{array}{l}\frac{\partial u}{\partial t}yv+\frac{\partial h}{\partial x}+\frac{u}{{\tau}_{\text{drag}}}=0,\\ \frac{\partial v}{\partial t}+yu+\frac{\partial h}{\partial y}+\frac{v}{{\tau}_{\text{drag}}}=0,\\ \frac{\partial h}{\partial t}+\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{h}{{\tau}_{\text{rad}}}=Q,\end{array}$$
where x and y are the horizontal coordinates, t is time, u is the zonal velocity (x direction), v the meridional (y direction), h the height (H) of the shallow water fluid minus the initial, horizontallyconstant height (H_{0}), that is, h = H − H_{0}, τ_{drag} the drag timescale, τ_{rad} the radiative timescale, and Q the heating function. The characteristic length, speed and time, corresponding to the Rossby deformation radius, the gravity wave speed, and the time for a gravity wave to cross a deformation radius in the shallow water system are: $$\begin{array}{l}L={\left({\beta}^{1}\sqrt{g{H}_{0}}\right)}^{1/2},\\ U=\sqrt{g{H}_{0}},\\ T={\left(\beta \sqrt{g{H}_{0}}\right)}^{1/2},\end{array}$$
respectively, where g is the gravitational acceleration assumed constant, and β = 2Ωcosϕ∕R or the derivative of the Rossby parameter with ϕ the latitude of the βplane. For the rest of this paper, we only consider ϕ = 0.
Equations (1)–(3) form a linear differential equation of the form $\partial X/\partial t=\mathcal{L}X+\mathcal{Q}$, where X = (u, v, h) is a vector of solutions, $\mathcal{L}$ a linear operator, and $\mathcal{Q}=(0,0,Q)$ is the vector form of the forcing. Hence, the solution is the sum of a homogeneous and a particular solution. The spatial part of the 3D solution can be expressed as an infinite sum of modes indexed by m with equivalent depth H_{m} instead of H_{0}. The orthogonal base functions are sinusoidal in z, the vertical coordinate, and of the form e^{i(mz)} with $m\in \mathbb{N}$ and the heating must be decomposed onto these functions (see Tsai et al. 2014 Sect. 2 for the rescaling of z and Wu et al. 2000 Sect. 2 for a discussion on boundary conditions).
When neither drag nor heating are considered, Matsuno (1966) expresses the analytic solutions to the homogeneous equations in the form of $\{u,v,h\}=\{\tilde{u},\tilde{v},\tilde{h}\}\text{\hspace{0.17em}exp}(\text{i}\omega t+\text{i}kx)$, where ω is the complex frequency, k the longitudinal wavenumber and a tilde denotes a function of y only. Dropping the tilde for simplicity, the homogeneous equations can be expressed as: $$\begin{array}{l}\text{i}\omega uyv+\text{i}kh=0,\\ \text{i}\omega v+yu+\frac{\partial h}{\partial y}=0,\\ \text{i}\omega h+\text{i}ku+\frac{\partial v}{\partial y}=0.\end{array}$$
Matsuno (1966) showed that this system can be reduced to a single equation for v, namely, $$\frac{{\partial}^{2}v}{\partial {y}^{2}}+\left({\omega}^{2}{k}^{2}+\frac{k}{\omega}{y}^{2}\right)v=0.$$(10)
Through an analogy with the Schrödinger equation of a simple harmonic oscillator, the boundary condition v → 0 when y→∞ requires $${\omega}^{2}{k}^{2}+\frac{k}{\omega}=2n+1,$$(11)
with $n\in \mathbb{N}$. As this is a third order equation, the eigenvalues for the frequency are labelled ω_{n,l} with l = 0, 1, 2, and the corresponding eigenvectors are labelled ${\tilde{X}}_{n,l}=({u}_{n,l},{v}_{n,l},{h}_{n,l})$. Finally, the case where n = 0 is treated separately and the important case where v is identically null is similar to a coastal Kelvin wave, with ω = −k (Matsuno 1966). The form of the solutions in the y direction are expressed through the use of the parabolic cylinder functions ψ_{n}, given by $${\psi}_{n}(y)={H}_{n}(y){e}^{{y}^{2}/2},$$(12)
where H_{n} is the nth Hermite polynomial. Finally, simple mathematical arguments show that the eigenvalue ω is always real: the homogeneous solutions are only neutral modes. Matsuno (1966) also shows that the eigenvectors of Eqs. (7)–(9) form a complete, orthogonal set of the 2D betaplane: at a given time, any function on the beta plane can be written as a linear combination of the ψ_{n} (y)exp(ikx) functions.
Matsuno (1966) and Gill (1980) obtained the steady state solution to Eqs. (1)–(3) under the inclusion of heating (and cooling) and drag. The completeness and orthogonality of the above functions allows us to write: $$\mathcal{Q}={\sum}_{n,l}{q}_{n,l}{\tilde{X}}_{n,l},$$(13)
where q_{n,l} is the projection of $\mathcal{Q}$ onto ${\tilde{X}}_{n,l}$. In their Eq. (34), Matsuno (1966) showed that a steady solution X to the forced problem with τ_{drag} = τ_{rad} is given by $$X={\sum}_{n,l}\frac{{q}_{n,l}}{{\tau}_{\text{drag}}^{1}\text{i}{\omega}_{n,l}}{\tilde{X}}_{n,l}.$$(14)
Showman & Polvani (2011) show that for a horizontal wavenumber one representing the asymmetric heating of hot Jupiters and τ_{drag} = τ_{rad} = 10^{5} s, the steady linear solution exhibits a “chevronshaped” pattern (in pressure, density, or temperature) which has been denominated the Matsuno–Gill (MG) solution, leading to a net acceleration of the equator at the nonlinear order. However, it is not clear whether a linear steady state is relevant in a case where nonlinearities are likely dominant, such as in the case of hot Jupiters, where the extreme forcing is likely to trigger nonlinear effects over short timescales. It also is not clear whether it is appropriate to choose equal values for both dissipation timescales. Therefore, we require the timedependent solution of Eqs. (1)–(3) in the general case, which are expressed in Sect. 2.3.
2.2 Nonlinear accelerations from the linear steady state
Now that we have reviewed the main assumptions and equations for our basic framework, we move on in this section to a summary of the key results from Showman & Polvani (2011) and Tsai et al. (2014). A key conclusion drawn by Showman & Polvani (2011) is that the MG pattern is a linear steady state but that the nonlinear accelerations from this circulation trigger an equatorial superrotation. In considering a linear perturbation (a wave or a steady linear circulation) associated with velocities u′, v′, w′ in the longitudinal, latitudinal, and vertical directions, respectively, the nonlinear momentum fluxes per unit mass from this perturbation scale as $$\begin{array}{l}{\varphi}_{l}\propto \overline{{u}^{\prime}{v}^{\prime}},\\ {\varphi}_{v}\propto \overline{{u}^{\prime}{w}^{\prime}},\end{array}$$
where ϕ_{l} is the latitudinal flux of momentum, ϕ_{v} the vertical and an overline denotes a longitudinal average. When ϕ_{l} is positive, there is a net transport of eastward momentum to the North, with a negative value resulting in a net transport to the South. When ϕ_{v} is positive, there is a net transport of eastward momentum upward or downwards when it is negative. Therefore, if ϕ_{l} is negative in the midlatitudes of the Northern hemisphere and positive in the Southern hemisphere, there is a net meridional convergence of eastward momentum (a similar argument applies in the vertical coordinate for a 3D systems). For the shallowwater βplane system used in Showman & Polvani (2011), the vertical momentum flux is accounted for by the addition of a coupling term between the deeper (highpressure), quiescent atmosphere and the dynamically active (lowerpressure) atmosphere (R term in Eqs. (9) and (10) of Showman & Polvani 2011).
In Fig. 1a, we present the temperature (color scale, K) and wind vectors (vector arrows) as a function of latitude and longitude for a typical MG circulation. This 3D linear steady state was obtained using ECLIPS3D (see Debras et al. 2019 for the details and benchmarking of ECLIPS3D), a linear solver for waves, instabilities,and linear steady states of an atmosphere under prescribed heating and drag. The initial state around which the equations are linearized is an axisymmetric, hydrostaticallybalanced state at rest. The bottom pressure is set to 220 bars and the temperature profile follows that of Iro et al. (2005), with the polynomial fit in the log of pressure of Heng et al. (2011), assuming anideal gas equation of state. Due to the very high inner boundary pressure, the choice of the inner boundary condition does not impact our results. The physical parameters relevant for HD 209458b used in this setup, namely the radius R_{p}, the rotation rate Ω, the depth of the atmosphere R_{top}, the surface gravity g_{p}, the inner boundary pressure p_{max}, the specific heatcapacity c_{p}, and the ideal gas constant R are given in Table 1. Finally, the heating as well as drag and radiative timescales were prescribed as in Komacek & Showman (2016) with the addition of an exponential decay in the heating in the upper, lowerpressure part of the atmosphere. The exponential damping acts to mimic the damping of vertical velocities close to the outer boundary, or “sponge layer” applied in 3D GCMs (see e.g., Mayne et al. 2014b). In turn, this damping layer allows the adoption of a “no escape” or reflective outer boundary condition, which would otherwise reflect waves back into the numerical domain. The equilibrium temperature towards which the atmosphere relaxes is a sinusoidal function of latitude and longitude, and ΔT_{eq}, the equilibrium day sidenight side temperature difference, decreases logarithmically in pressure between 10^{−3} bar where ΔT_{eq} = ΔT_{eq,top} (the value at the top of the atmosphere which is set to 10 K in this case) to 10 bars where ΔT_{eq} = 0. The drag timescale, τ_{drag}, is constant throughout the atmosphere at 10^{5} s and the radiative timescale, τ_{rad}, is a logarithmically increasing function of pressure (see Eq. (7) of Komacek & Showman 2016) between 10^{−2} bar where τ_{rad} = τ_{rad,top} = 10^{5} s and 10 bar where τ_{rad} = 10^{7} s.
As shown in Fig. 1a, the maximum temperature at the equator is shifted eastward from the substellar point (the substellar point is set at a longitude of 180 degrees) in our results consistent with observations (Knutson et al. 2007; Zellem et al. 2014). The meridional circulation exhibits a Rossby wavetype circulation at mid latitudes, with clockwise or counterclockwise rotation around the pressure maxima, and a Kelvin wave type circul ation at theequator, with no meridional velocities. The combination of both of these circulations brings eastward momentum to the equator to the East of the substellar point, and advects westward momentum to the midlatitudes to the West of the substellar point. Globally, it is easily shown that ϕ_{l} is indeed negative in the Northern hemisphere and positive in the Southern hemisphere: there is a net convergence of eastward momentum at the equator. According to Showman & Polvani (2011), this convergence is associated with divergence of vertical momentum flux, and equilibration occurs when the vertical and meridional terms balance.
Tsai et al. (2014) further extend this understanding to include vertical transport more completely. Tsai et al. (2014) project the heating function onto equivalent height βplane solutions (see Sect. 2.1), showing that the vertical behavior of the waves can be linked to the equilibration of the jet (a synonym used here and throughout this work for the term “equatorial superrotation”). More precisely, Tsai et al. (2014) show that in the limit of slow evolution (or strong dissipation), their linear development around a steady flow with constant background zonal velocity reproduces the wave processes occurring in 3D simulations extremely well. Tsai et al. (2014) show that the wave response of the atmosphere is shifted from west to east when the background zonal velocity is increased (their Fig. 10): this is interpreted as a convergence toward a single equilibrium state, where the nonlinear acceleration from the linear processes are canceled out. While they are very detailed and physically relevant, the results of Tsai et al. (2014) are, as they state, only applicable to the strong or modest damping scenario, dictated by the fact that the waves must have the time to reach a stationary state before nonlinearities become significant. Throughout this work, we define the drag regime relative to the initial acceleration of the jet: in the “weak” drag regime, nonlinear terms become nonnegligible before a linear steady state (MG) could have been reached; whereas in the “modest” drag regime, the time required to reach the MG state is comparable with the time to depart from this linear steady state; and in the “strong” drag regime, we can decouple the linear and the nonlinear evolution of the planet, as considered by Showman & Polvani (2011). Once an initial jet has been accelerated, the evolution of the atmosphere is much slower than its initial acceleration and the results of Tsai et al. (2014) apply, even in a weak drag regime, explaining the consistency of their work for the evolution of the jet toward an equilibrated state.
Komacek & Showman (2016) compare the steady states from various 3D GCM simulations across a range of τ_{rad} and τ_{drag} values (their Figs. 4 and 5). The simulations of Komacek & Showman (2016) extend from lowforcing, hence a linear steady state, to strongforcing, hence a nonlinear steady state. Contrary to the conclusions of Showman & Polvani (2011), Komacek & Showman (2016) also show that when the linear steady state resembles that of Fig. 1a, the associated nonlinear steady state is not (or only weakly) superrotating. This can be understood by considering the fact that τ_{drag} is smaller than the characteristic timescale of advection by the superrotating jet over the whole planet, hence the jet is dissipated before it can reach a steady state. In Fig. 1b, we present the results from ECLIPS3D that are obtained when reproducing a particular setup from Komacek & Showman (2016), namely, that with τ_{drag} = 10^{6} s and τ_{rad,top} = 10^{4} s. According to the analysis of Komacek & Showman (2016), the nonlinear steady state associated with the linear steady state of Fig. 1b does indeed exhibit strong superrotation, although the tilt of the wave in Fig. 1b would lead to a removal of the momentum at the equator. Komacek & Showman (2016) acknowledge this in writing, “these phase tilts are the exact opposite of those that are needed to drive superrotation”. As for the explanation of Showman & Polvani (2011), the stationary wave pattern obtained from the heating is postulated to accelerate superrotation but Komacek & Showman (2016) present results which are in opposition to this scenario: when the linear steady state accelerates the equator (Fig. 1a, with τ_{drag} = 10^{5} s) the associated nonlinear steady state is not superrotating. However, when the linear steady state takes momentum away from the equator (Fig. 1b, with τ_{drag} = 10^{6} s), the nonlinear steady state is superrotating. Thus, there is an inherent discrepancy between Showman & Polvani (2011) and Komacek & Showman (2016). In order to understand this discrepancy, we need to go a step beyond the sole consideration of a linear steady state and study the evolution of the linear solution over time. This is the objective of Sects. 2.3 and 4.
Fig. 1 Temperature (colorscale in K) and horizontal wind (arrows) as a function of longitude (x axis) and latitude (y axis) at the 40 mbar pressure level of the linear steady state (denominated MG circulation) obtained using ECLIPS3D (Debras et al. 2019) with heating function, drag and radiative timescales following the definitions of Komacek & Showman (2016). Following the notation of Komacek & Showman (2016): (a) ΔT_{eq,top}= 100 K, τ_{drag,top} = 10^{5} s and τ_{rad,top} = 10^{5} s. The maximum speed at this pressure range is 10 m s^{−1}. (b) ΔT_{eq,top} = 100 K, τ_{drag,top} = 10^{6} s and τ_{rad,top} = 10^{4} s. The maximum speed at this pressure range is 100 m s^{−1}. Note that the maximum speed has been multiplied by ten as the drag timescale has been multiplied by ten. 
Value of the standard parameters for HD 209458b, following Mayne et al. (2014b).
2.3 Timedependent solutions
Now that we have established the basic mathematical system and summarized the current picture of superrotation in hot Jupiter atmospheres, we move on to expressing the timedependent solution to the forced problem which provides us with the shape of the atmospheric wave response. Our main assumption is that the heating function can be decomposed onto the homogeneous solutions of Eqs. (17)–(19). When τ_{rad} = τ_{drag}, the horizontal part (defined as X_{n,l}(x, y, t = 0) in Eq. (20)) of the eigenvectors still form a complete set of the equatorial βplane because of the orthogonality and completeness of the Hermite functions, as in Matsuno (1966). However, when τ_{rad} ≠τ_{drag}, the eigenvectors are no longer orthogonal, as shown in Appendix A, and a rigorous proof would be needed to show that they still form a complete set of solutions^{1}. From a physical perspective, it is expected that the heating function would trigger linear waves which are solutions of the homogeneous equation and such a decomposition of the heating function onto these waves therefore probably exists, although it is no longer simply given by a scalar product. Finally, it is worth stating that solving $\partial X/\partial t=aX+\mathcal{Q}$ is straightforward (except that we don’t know the eigenvalues and eigenvectors of the homogeneous equation), however, employing a Green’s function to solve this equation provides a more physically intuitive result in terms of wave propagation and dissipation.
With the addition of a drag timescale, τ_{drag}, and a radiative timescale, τ_{rad}, Eqs. (7)–(9) can be modified to yield $$\begin{array}{l}\left(\text{i}\omega +\frac{1}{{\tau}_{\text{drag}}}\right)uyv+\text{i}kh=0,\\ \left(\text{i}\omega +\frac{1}{{\tau}_{\text{drag}}}\right)v+yu+\frac{\partial h}{\partial y}=0,\\ \left(\text{i}\omega +\frac{1}{{\tau}_{\text{rad}}}\right)h+\text{i}ku+\frac{\partial v}{\partial y}=0.\end{array}$$
Indexing again the solutions by n and l as in Matsuno (1966), we define $${X}_{n,l}=({u}_{n,l},{v}_{n,l},{h}_{n,l})={\tilde{X}}_{n,l}(y){e}^{\text{i}kx+(\text{i}{\nu}_{n,l}{\sigma}_{n,l})t}$$(20)
as an eigenvector of Eqs. (17), (18), and (19), with ${\tilde{X}}_{n,l}(y)$ the amplitude of the wave, k the horizontal wave number, ν_{n,l} its frequency, and σ_{n,l} its damping (or growing/growth) rate (note that ω_{n,l} = ν_{n,l} + iσ_{n,l}). We also define $\mathcal{L}$ as the operator of the same equations, such that $\mathcal{L}{X}_{n,l}=0$ for all n and l. The general equation can then be written as $\mathcal{L}{X}_{\text{F}}=\mathcal{Q}$, where $\mathcal{Q}$ is the forcing which, as it is only present in the third individual equation is given by $\mathcal{Q}=(0,0,Q)$, and X_{F} is the forced, timedependent solution. A homogeneous solution X_{H} can be written in its general form as $${X}_{\text{H}}={\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l},$$(21)
where α_{n,l} are scalars.
When τ_{drag} = τ_{rad}, ν_{n,l} are similar to the ω_{n,l} of Matsuno (1966) and ${\sigma}_{n,l}={\tau}_{\text{drag}}^{1}$ for all (n, l). When τ_{drag}≠τ_{rad}, the analytical expressions for ν_{n,l} and σ_{n,l} are not known a priori. In order to solve the general equation, we seek the causal Green function X_{G} that represents the solution at time t due to switching on the forcing at time t′ only. Therefore, for all t and t′: $$\mathcal{L}{X}_{\text{G}}(x,y,t,{t}^{\prime})=\delta (t{t}^{\prime})F(x,y,{t}^{\prime}),$$(22)
where δ(t) is the Dirac distribution and F is the heating function. In the case of simulated hotJupiter atmospheres, the star is effectively “switched on” at t = 0, after which the heating is constant with time (in the linear limit). F can then simply be expressed as $F(x,y,{t}^{\prime})=\Theta ({t}^{\prime})\mathcal{Q}(x,y)$, where Θ(t) is the Heavyside function (null when t < 0 and equal to1 otherwise). The forced solution of Eqs. (17)–(19) is simply the integral over t′ of the causal Green function, hence the sum of the responses of the atmosphere excited by a Dirac function of the forcing at time t′: $${X}_{\text{F}}={\displaystyle \underset{\infty}{\overset{\infty}{\int}}{X}_{\text{G}}}(t{t}^{\prime})\Theta ({t}^{\prime})\text{d}{t}^{\prime}={\displaystyle \underset{0}{\overset{t}{\int}}{X}_{\text{G}}}(t{t}^{\prime})\Theta ({t}^{\prime})\text{d}{t}^{\prime},$$(23)
where the change in the upper limit of integration can be made due to the fact that the Green function is causal and is, therefore, zero when t − t′ is negative. From the definition of the Green function (Eq. (22)), for (t− t′) > 0 we have $\mathcal{L}{X}_{\text{G}}=0$. Hence, X_{G} is a homogeneous solution of Eqs. (17)–(19) when t > t′. It is then logical to choose for X_{G}: $${X}_{\text{G}}(t,{t}^{\prime})=\Theta (t{t}^{\prime}){\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t{t}^{\prime}),$$(24)
and it is easily verified that $$\begin{array}{ll}\mathcal{L}{X}_{\text{G}}=\hfill & \text{\hspace{0.17em}}\delta (t{t}^{\prime}){\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t{t}^{\prime})\hfill \\ \hfill & +\text{\hspace{0.17em}}\Theta (t{t}^{\prime})\frac{\partial}{\partial t}{\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t{t}^{\prime})\hfill \\ \hfill & +\text{\hspace{0.17em}}{\mathcal{L}}_{h}\Theta (t{t}^{\prime}){\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t{t}^{\prime}),\hfill \\ =\hfill & \text{\hspace{0.17em}}\delta (t{t}^{\prime}){\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t{t}^{\prime})+\Theta (t{t}^{\prime})\mathcal{L}{\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t{t}^{\prime}),\hfill \\ =\hfill & \text{\hspace{0.17em}}\delta (t{t}^{\prime}){\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t{t}^{\prime}),\hfill \end{array}$$(25)
where we have used the fact that the derivative of the Heavyside function is the Dirac distribution, $\mathcal{L}=\partial /\partial t+{\mathcal{L}}_{h}$ where ${\mathcal{L}}_{h}$ is an operator acting on the horizontal coordinates only and $\mathcal{L}{X}_{n,l}=0$. Therefore, in order to solve the forced problem we can project the forcing onto the homogeneous solutions and write for t = t′: $${\sum}_{n,l}{\alpha}_{n,l}{X}_{n,l}(t={t}^{\prime})={\sum}_{n,l}{\alpha}_{n,l}{\tilde{X}}_{n,l}=\mathcal{Q}(x,y).\text{\hspace{0.17em}}$$(26)
The first equality simply arises from the definition of ${\tilde{X}}_{n,l}$, whereas the second uses Eq. (22). Therefore, if we know the projection of $\mathcal{Q}$ onto the ${\tilde{X}}_{n,l}$, the α_{n,l} quantities, the final solution can be obtained. By setting α_{n,l} = q_{n,l}, we recover the results of the previous section (these are termed b_{m} and b_{m,n,l} in Matsuno (1966) and Tsai et al. (2014), respectively, where the latter is in 3D). The solution to the forced problem is given by injecting the Green function (Eq. (24)) into Eq. (23): $$\begin{array}{ll}{X}_{\text{F}}\hfill & ={\displaystyle \underset{\infty}{\overset{\infty}{\int}}{\sum}_{n,l}}{q}_{n,l}{X}_{n,l}(t{t}^{\prime})\Theta (t{t}^{\prime})\Theta ({t}^{\prime})\text{d}{t}^{\prime},\hfill \\ \hfill & ={\sum}_{n,l}{\displaystyle \underset{0}{\overset{t}{\int}}{q}_{n,l}}{\tilde{X}}_{n,l}(x,y){e}^{(i{\nu}_{n,l}{\sigma}_{n,l})(t{t}^{\prime})}\Theta ({t}^{\prime})\text{d}{t}^{\prime}.\hfill \end{array}$$(27)
Under this integral form, it is clear that the solution is the continual excitation (the Θ term) of waves with a characteristic frequency ν_{n,l} and time of decay of ${\sigma}_{n,l}^{1}$. The amplitude of the excited waves is proportional to their projection onto the forcing function Q, as explained by Gill (1980) and Tsai et al. (2014). Solving this integral yields: $${X}_{\text{F}}={\sum}_{n,l}\frac{{q}_{n,l}{\tilde{X}}_{n,l}}{{\sigma}_{n,l}i{\nu}_{n,l}}\left(1{e}^{(i{\nu}_{n,l}{\sigma}_{n,l})(t)}\right)\text{}.$$(28)
In this form, we simply recover the results of Matsuno (1966) and, notably, their Eq. (34) or our Eq. (14), $${X}_{\text{MG}}={\sum}_{n,l}\frac{{q}_{n,l}{\tilde{X}}_{n,l}}{{\sigma}_{n,l}i{\nu}_{n,l}},$$(29)
where X_{MG} is the MG solution, hence the steady solution to the forced problem. The timedependent part of the solution could have been obtained from a simple firstorder equation solution. However, the Green function formalism allows us to determine that the solution consists of permanentlyforced waves that are damped with time and that the shape of the atmosphere is given by the interactions between these waves. Using Eq. (27), we can write $${X}_{\text{MG}}=\underset{t\to \infty}{\mathrm{lim}}{\sum}_{n,l}{\displaystyle \underset{0}{\overset{t}{\int}}{q}_{n,l}}{\tilde{X}}_{n,l}{e}^{(i{\nu}_{n,l}{\sigma}_{n,l})(t{t}^{\prime})}\Theta ({t}^{\prime})\text{d}{t}^{\prime},$$(30)
the form of which confirms the interpretation of the stationary solution as an infinite interaction of waves. Additionally, it shows that for a given heating function, changing the value of τ_{drag} (but keeping τ_{rad} = τ_{drag}), thereby not altering the $\tilde{{X}_{n,l}}$ and ν_{n,l} but only ${\sigma}_{n,l}={\tau}_{\text{drag}}^{1}$, will change the linear steady solution. This is because the excited waves will not propagate to the same length before being damped. This was first realized by Wu et al. (2001), where they show that the zonal wave decay length is of the order of $\sqrt{{\tau}_{\text{rad}}^{1}{\tau}_{\text{drag}}^{1}}$ for arbitrary τ_{rad} and τ_{drag}.
From Eq. (27), we see that the linear solution is controlled by three main parameters: the shape of the forcing function (q_{n,l}), the global behavior of the waves (horizontal shape of ${\tilde{X}}_{n,l}$ and ν_{n,l}), and the dissipation of the waves (σ_{n,l}). Although Eq. (27) provides the form of the solution to the timedependent problem, the three main parameters we have detailed remain unknown. Therefore, we move on to quantifying the sensitivity of the solution to the forcing function and, hence, the influence of the q_{n,l} in the next section.
3 Insensitivity of the Matsuno–Gill solution to differential heating
The interpretation of Showman & Polvani (2011) relies on a simplified treatment of the forcing, the impact of which we must first assess before discussing the impact of the linear evolution of the atmosphere. Firstly, in order to derive analytical results, Showman & Polvani (2011) impose an antisymmetric (sinusoidal) heating where the night side of the planet is cooled as much as the day side is heated. In this case, the linear steady solution gives rise to the chevronshaped pattern they denominate the “MatsunoGill” circulation. However, the actual structure of the heating process is not just a sinusoidal function. Moreover, we know from Mayne et al. (2017) and Amundsen et al. (2016) that there are qualitative differences between the steady state of GCM simulations of hot Jupiters calculated using either a Newtonian heating or with a more sophisticated radiative transfer scheme. In that regard, the first intuitive idea to test is whether the MG pattern is robust when the heating function is changed. With the addition of the vertical dimension, Tsai et al. (2014) have shown that the linear solution strongly resembles the MG circulation at low pressures in the atmosphere. It is not clear whether this holds with realistic, threedimensional heating functions.
From Eqs. (14) or (30), as ${X}_{n,l}(x,y)={\tilde{X}}_{n,l}(y)\mathrm{exp}(\text{i}kx)$, we know that the projection of the heating function onto different wavenumbers, k, can alter the resulting Matsuno–Gill circulation from that obtained at wavenumber one. To verify this point, using ECLIPS3D we solve for linear steady circulations across a range of prescribed heating rates, adopting the drag and radiative timescales of Komacek & Showman (2016), as used previously for the data presented in Fig. 1a. We use three forms of heating, the first two modified from that used for Fig. 1a, and a third one, inspired from 3D GCM simulations: (i) he same dayside forcing as Fig. 1a, but no nightside cooling ; (ii) the same forcing as Fig. 1a but with the nightside cooling enhanced by a factor of 3; (iii) a heating profile matching that obtained from 3D GCM simulations (taken from Amundsen et al. 2016), including net day side heating and night side cooling.
The first two cases allow us to test the robustness of the pattern under extreme situations, with the third mimicking the GCM simulations. ECLIPS3D calculates the linear steady states by inverting the linear matrix obtained with the full 3D equations (see description in Debras et al. 2019). The outer boundary condition is a “no escape” condition (zero vertical velocity w), and we have applied an exponential decay of the forcing in the high atmosphere (low pressure) to mimic the damping of vertical velocities, or sponge layer, employed in the Unified Model (Mayne et al. 2014a) and other GCMs. For the inner boundary, we adopt a solid boundary condition for the vertical velocity (w = 0), and a free slip (no vertical derivatives of the horizontal velocities, u and v) or no slip (no horizontal velocities) condition on u and v give qualitatively similar results because of the very high inner pressure. These assumptions are obviously simplifications of the real physical situation, but we have assessed their impact by changing the range of pressures at the inner and outer boundaries and we find no significant change in the qualitative results. The physical parameters adopted for HD 209458b are the same as in Mayne et al. (2017), which are also given in Table 1.
Figure 2 shows the perturbed temperature (steady temperature minus initial temperature, colorscale), and winds (arrows) as a function of longitude and latitude of the MG solutions for the three different heating profiles described previously, at a height where the initial pressure is 50 mbar. Figure 2 shows that the shape of the linear steady circulation does not qualitatively depend on the forcing: we always recover an eastward shift of the hot spot, associated with a tilt of the eddy patterns leading to a net acceleration of the equator for the nonlinear order. As long as there is a differential heating between the day and the night side, the linear steady solution of the atmosphere exhibits the chevronshaped pattern of the MG circulation (we have verified that a constant heating across the whole planet or just a cooling on the night side does not lead to solutions of a similar form as the MG solutions). There is, however, a change in the quantitative values, which does not affect the qualitative aspects of the nonlinear momentum transfer (see Sect. 5). Therefore, relaxing the approximation of a wavenumber one (e.g., daynight antisymmetric forcing) heating function has no influence on the nonlinear acceleration around the linear steady state: the projection onto different zonal wavenumbers changes the zonal amplitude of the MG circulation, but not its qualitative shape. We only present the results with a characteristic drag and radiative timescale of 10^{5}s but we haveverified that changing these values affects the shape of the solutions in all cases in a similar way.
Globally, even for simulations with a proper treatment of the radiative transfer, as long as the planet is tidally locked, we expect the linear steady circulation to be MGlike. In the paradigm proposed by Showman & Polvani (2011) and Tsai et al. (2014), the propagation of planetary waves impose this global MG circulation, exhibiting no superrotation, with the acceleration on to superrotation being due to eddy mean flow interactions around this primordial state. Therefore, the time to reach the linear steady state, which is by no means a nonlinear steady state, must be small relatively to other dynamical times in the system. In the case where the drag and radiative timescales are equal, this led Tsai et al. (2014) to conclude that their work could not apply to long diffusion timescales. We now seek to assess if similar conclusions can be made in cases where the drag and radiative timescales differ, in order to determine how the atmosphere reacts at first order. To develop this argument, we analytically restrict ourselves to the quasi geostrophic set of equations (2D Cartesian, shallow water βplane), as detailed in Sect. 2.1.
Fig. 2 Temperature (colorscale in K) and winds (arrows) as a function of longitude and latitude with different heating functions, with drag and radiative timescales defined as in Komacek & Showman (2016) with nondimensional values of 1 (almost half an Earth day). (a–c): the initial pressure at this height is 50 mbar which corresponds to an initial temperature of 1326 K and the forcing is ΔT_{eq,top} = 50 K. (d): the initial pressure is 40 mbar and initial temperature 1322 K. (a): heating as in Komacek & Showman (2016) (b): same as in Komacek & Showman (2016) but with a cooling at the night side three times more efficient than the heating on the day side. (c): same as in Komacek & Showman (2016) with no cooling at the night side. (d): heating rate extracted from the full radiative transfer calculations of the GCM (divided by ten to obtain comparable values). 
4 Wave propagation and dissipation
We have established, using a Green’s function, that the linear solution to the 2D shallow water, βplane equations, can be expressed as the continual excitation of waves with a characteristic frequency and decay timescale (Eq. (27), Sect. 2.3). The decay timescales themselves are crucial as they can be used to determine the qualitative form of the linear steady state itself, and provide insight into the response of the atmosphere over short timescales. In this section, we focus on the mathematical derivation of the characteristic decay timescales reaching an expression (Sect. 4.1.1), which we then solve numerically for different types of atmospheric waves (gravity; Rossby and Kelvin waves). A short summary is then provided in Sect. 4.1.5. Following this, we extend our arguments to 3D in Sect. 4.2. The entirety of this section is focused on the mathematical nature of the supported waves, with a physical interpretation provided later in Sect. 5.
4.1 Characteristics of waves in 2D
4.1.1 Decay time of damped waves
The first task is to derive an equation for the decay timescale, or growth rate, for a general damped wave in our framework. If τ_{drag} = τ_{rad}, using the complex frequency $\omega \in \u2102$ (contrary to Matsuno (1966) where $\omega \in \mathbb{R}$), we can define iλ = iω + 1∕τ_{drag} to transform Eqs. (17)–(19) into Eqs. (7)–(9), equations of which we know the eigenvalues from Matsuno (1966). Therefore, λ is real which implies that the realpart of ω, the frequency, is unaltered from the original equations of Matsuno (1966), but the imaginary part of ω becomes, $$\Im (\omega )=\frac{1}{{\tau}_{\text{drag}}}.$$(31)
One can then express e.g., u_{n,l} as: $${u}_{n,l}={\tilde{u}}_{n,l}{\text{e}}^{t/{\tau}_{\text{drag}}}{\text{e}}^{\text{i}({\nu}_{n,l}t+kx)}.$$(32)
This shows that all modes decay over a characteristic timescale, namely, the drag (or radiative) timescale. For the case where τ_{rad} = τ_{drag} the time to converge to the Matsuno–Gill circulation is the decay timescale, as one would naively expect, and all waves have the same exponential decay in time.
If τ_{drag}≠τ_{rad}, Eq. (10) must be modified to obtain: $$\begin{array}{l}\frac{{\partial}^{2}v}{\partial {y}^{2}}(\left(\text{i}\omega +{\tau}_{\text{drag}}^{1}\right)\left(\text{i}\omega +{\tau}_{\text{rad}}^{1}\right)+{k}^{2}\\ \text{\hspace{1em}}+{y}^{2}\frac{\text{i}\omega +{\tau}_{\text{rad}}^{1}}{\text{i}\omega +{\tau}_{\text{drag}}^{1}}\frac{ik}{\text{i}\omega +{\tau}_{\text{drag}}^{1}})v=0.\end{array}$$(33)
It is easy to verify that setting ${\tau}_{\text{drag}}^{1}={\tau}_{\text{rad}}^{1}=0$ gives back Eq. (10).
To go one step further, we define a complex number c such that: $${c}^{4}=\frac{\text{i}\omega +{\tau}_{\text{rad}}^{1}}{\text{i}\omega +{\tau}_{\text{drag}}^{1}},$$(34)
and choose for c the only root with a positive real and imaginary part. This definition ensures that the real part of c^{2} is always positive, which allows for the solutions to decay at infinity, see Appendix A. We then choose as a variable z = cy (the cases c = 0 and c = ∞ are of no physical interest). Using this new variable, ∂^{2}v∕∂y^{2} = c^{2}∂^{2}v∕∂z^{2}, Eq. (33) is simplified to $${c}^{2}\frac{{\partial}^{2}v}{\partial {z}^{2}}\left(\left(\text{i}\omega +{\tau}_{\text{drag}}^{1}\right)\left(\text{i}\omega +{\tau}_{\text{rad}}^{1}\right)+{k}^{2}+{y}^{2}{c}^{4}\frac{\text{i}k}{\text{i}\omega +{\tau}_{\text{drag}}^{1}}\right)v=0.$$(35)
Dividing this equation by c^{2} (and recalling the definition of z, z = cy), we obtain $$\frac{{\partial}^{2}v}{\partial {z}^{2}}\frac{\left(\text{i}\omega +{\tau}_{\text{drag}}^{1}\right)\left(\text{i}\omega +{\tau}_{\text{rad}}^{1}\right)+{k}^{2}\text{i}k/\left(\text{i}\omega +{\tau}_{\text{drag}}^{1}\right)}{{c}^{2}}v{z}^{2}v=0.$$(36)
Defining themultiplier of v in the second term as m, the equation can be expressed as $$\frac{{\partial}^{2}v}{\partial {z}^{2}}+(m{z}^{2})v=0.$$(37)
The major difference with the Matsuno case, Eq. (10), is that now $z\in \u2102$ and so the boundary conditions are altered. As z→∞ when y →±∞, we have to solve this equation with the following boundary condition: $$v\to 0\text{\hspace{0.17em}\hspace{0.17em}}\text{\hspace{0.17em}when}\text{\hspace{0.17em}\hspace{0.17em}}\text{\hspace{0.17em}}\leftz\right\to \infty .$$(38)
As in the case where m is real, one could show (see, e.g., Abramowitz & Stegun 1965) that the only solutions are the parabolic cylinder functions, Eq. (12): ${H}_{n}(cy){\text{e}}^{{c}^{2}{y}^{2}/2}$ where $n\in \mathbb{N}$, hence the need for ℜ(c^{2}) > 0 so that the solutions decay when y →±∞, and provided that: $$\frac{\left(\text{i}\omega +{\tau}_{\text{drag}}^{1}\right)\left(\text{i}\omega +{\tau}_{\text{rad}}^{1}\right){k}^{2}+\text{i}k/\left(\text{i}\omega +{\tau}_{\text{drag}}^{1}\right)}{{c}^{2}}=2n+1,$$(39)
Defining $x=\text{i}\omega +{\tau}_{\text{drag}}^{1}$ and $\gamma ={\tau}_{\text{rad}}^{1}{\tau}_{\text{drag}}^{1}$, followed by taking the square of Eq. (39) in order to obtain c^{4} yields, $$\begin{array}{l}{x}^{4}{(x+\gamma )}^{2}+2{k}^{2}{x}^{3}(x+\gamma )2\text{i}k{x}^{2}(x+\gamma )\\ \text{\hspace{1em}}{(2n+1)}^{2}x(x+\gamma )+{k}^{4}{x}^{2}2\text{i}{k}^{3}x{k}^{2}=0.\end{array}$$(40)
Equation (40) has already been obtained by Heng & Workman (2014) (their Eq. (121) with different notations: their F_{0} is ${\tau}_{\text{rad}}^{1}$ in our work, ω_{R} is our ω and ω_{I} is ${\tau}_{\text{drag}}^{1}$), in order to derive steady state solutions as performed in Wu et al. (2001) and Showman & Polvani (2011). Equation (40) is a polynomial of order six, but a thorough study of Eq. (39), as reported in Appendix B.1, shows that only three different waves propagate, and two for n = 0, as in Matsuno (1966).
The horizontal shape of the solutions to Eqs. (17)–(19) in the general case are given by Eq. (A.7), but we would also require the solutions of Eq. (40) to obtain a fully analytic expression for the waves. Therefore, we have solved Eqs. (17)–(19) numerically over an extensive range of n, k, τ_{drag} and τ_{rad} values (we have verified that our numerical solutions recover the limits ${\tau}_{\text{drag}}^{1}={\tau}_{\text{rad}}^{1}=0$ and ${\tau}_{\text{drag}}^{1}={\tau}_{\text{rad}}^{1}$). Firstly, as expected, no mode can exponentially grow given a background state at rest. The cases of k ~ 1 and n = 1, 3, 5 and 7 are the most important for hot Jupiters as the heating function is dominated by wavenumber 1, that is, a dayandnight side (nondimensional value around 0.7, see Sect. 5). Here the n number represents the order of the Hermite polynomial, hence the number of zero nodes in latitude (note that as c^{2} can be complex, the number of zeros in latitude is no more solely defined by the Hermite polynomials as in the c = 1 case). If the heating function is a Gaussian function, as chosen by Matsuno (1966), Showman & Polvani (2011) and Tsai et al. (2014), then the projection of $\mathcal{Q}$ on to the parabolic cylinder function stops at the third order (this is not necessarily true when τ_{drag}≠τ_{rad} but we don’t expect large amplitude in the n >3 waves as the forcing exhibits no zeros in latitude).
Using our numerical solutions, we explore the behavior of gravity (Sect. 4.1.2) and Rossby (Sect. 4.1.3) and Kelvin waves (Sect. 4.1.4) before summarising our results (Sect. 4.1.5). For all cases, the frequency of the waves remains within an order of magnitude of the free wave frequency (when ${\tau}_{\text{drag}}^{1}={\tau}_{\text{rad}}^{1}=0$). Hence the major changes between cases are obtained for the decay rate and horizontal shapes. The shapes of the waves with nonzero τ_{drag} and τ_{rad} are shown in Appendix C. Regarding the decay rates, we express them as power laws fitting reasonably well with the numerical values in the next section. We also implemented the 2Dshallow water equations in ECLIPS3D (detailed in Debras et al. 2019)to verify our numerical results discussed in this section. For all cases for matching parameters (characteristic values, τ_{rad}, τ_{drag}), the agreement between ECLIPS3D and the results discussed here (obtained using the Mathematica software) for both the decay timescales or growth rates and frequency of waves is excellent. Furthermore, inserting the numerical values obtained here in Eq. (A.7) recovers the exact modes as obtained using the shallow water version of ECLIPS3D.
4.1.2 Gravity waves
In this work we define “gravity waves” as the solutions to Eqs. (17)–(19) which tend to the standard definition of a gravity wave in the limit τ_{drag} → τ_{rad} (we have verified that the identification is unchanged for τ_{rad} → τ_{drag}). As there are only three solutions to these equations (see, e.g., Matsuno 1966), the Rossby wave is, therefore, the last mode.
Our numerical results give a characteristic time of decay for gravity waves of ~ τ_{drag}, when τ_{drag} ~ τ_{rad}, as expected. For cases where the drag is dominant over the radiative forcing, τ_{drag} ≪ τ_{rad}, the decay timescale obtained numerically is ~τ_{drag}, which indicates that the drag controls the timescale of the convergence of the atmosphere. Physically, this is expected as drag will preventthe wave from propagating and damp the perturbed velocity efficiently, preventing the temperature and pressure to depart significantly from the forced equilibrium profile.
However, when the radiative forcing is dominant over the drag, τ_{drag} ≫ τ_{rad}, we find two cases. Firstly, when τ_{rad} ≪ 1 the numerically obtained decay rate for the gravity waves (σ_{g1}) is given by $${\sigma}_{\text{g1}}\approx {\tau}_{\text{drag}}^{1}+{\tau}_{\text{rad}}^{1/(n+2)}.$$(41)
For this case, if τ_{rad} → 0, ${\sigma}_{\text{g1}}^{1}$ (the decay timescale) is given by the drag timescale, as for the case of dominant drag. However, if τ_{rad} is larger the behavior is more complex and includes a dependence on the order of the Hermite polynomial n, although this still results in the decay timescale being the same order of magnitude as the drag timescale. We interpret it based on the fact that although the temperature and pressure perturbation will be dictated by the forcing, the time to damp the wave is still governed by the time for the velocities to be damped, hence, the drag timescale.
Secondly, for the case where τ_{rad} ≫ 1 the numerically obtained limit for the decay rate (σ_{g2}) is given by $${\sigma}_{\text{g2}}\approx \frac{{\tau}_{\text{rad}}^{1}}{3}.$$(42)
In this case, the radiative timescale is long enough to be imposed as the characteristic time of decay even for the velocities, and the decay of the wave is only controlled by this timescale.
4.1.3 Rossby waves
The behavior of the Rossby wave decay timescale is more complex than that of gravity waves. When ${\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1}\gtrsim 0.5$, for all individual values of ${\tau}_{\text{drag}}^{1}$ or ${\tau}_{\text{rad}}^{1}$, the absolute value of the imaginary part of c^{2} is much larger than that of the real component. This means that the Rossby wave modes oscillate in the y direction several times before being damped, in these conditions. Additionally, for increasing values of ${\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1}$, the amplitude at the equator becomes negligible, and the mode’s peak amplitude moves to higher latitudes, where the equatorial βplane approximation begins to break down. Therefore, our numerical results show that the simple shallowwater, equatorial βplane framework adopted in this work is not usefully applicable to the case of Rossby waves where ${\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1}\gtrsim 0.5$. This is confirmed by the graphical representation of these waves in Appendix C. We will therefore rely on numerical results of Sect. 4.2 for this region of the parameter space. However, for ${\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1}\lesssim 0.5$ (which is the case for all ${\tau}_{\text{drag}}^{1},{\tau}_{\text{rad}}^{1}<0.5$), the decay rate for Rossby waves (σ_{R}) we have derived numerically can then be approximated by, $${\sigma}_{\text{R}}\approx \frac{1}{2}\left({\tau}_{\text{rad}}^{1}+{\tau}_{\text{drag}}^{1}\right)\text{\hspace{0.17em}}.$$(43)
Therefore, for long radiative and drag timescales, Rossby waves are equally sensitive to the damping of velocities and temperature. Such a result is expected from the conservation of potential vorticity, which gives rise to Rossby waves, and is defined in the shallow water system as (ξ + f)∕h where ξ =∇∧u is the vorticity and f =2Ωcosϕ the Coriolis parameter. When neither ξ and h are strongly damped, we might therefore expect a combination of the drag and radiative damping to return to equilibrium. Comparing the decay timescale for Rossby waves with that obtained for gravity waves in Sect. 4.1.2 shows that the decay rates differ between these two cases.
4.1.4 Kelvin waves
Kelvin waves are a particular solution of the homogeneous Eqs. (17)–(19), as first pointed out by Matsuno (1966), where the meridional velocity is zero, and can be characterized analytically. Combining Eqs. (17) and (19) with v =0 yields, $$\frac{\partial u}{\partial y}=\frac{\text{i}k}{\text{i}\omega +1/{\tau}_{\text{drag}}}yu,$$(44)
and hence $$u=A\text{exp}\left(\frac{\text{i}k}{\text{i}\omega +1/{\tau}_{\text{drag}}}\frac{{y}^{2}}{2}\right),$$(45)
where A is a constant. If the boundary condition u = 0 for y →±∞ is to be satisfied, the factor of y^{2}∕2 must have a negative real component. Additionally, in order for u and h not to be identically zero, ω must be a root of a second order polynomial given by $${\omega}^{2}\text{i}\omega \left(\frac{1}{{\tau}_{\text{rad}}}+\frac{1}{{\tau}_{\text{drag}}}\right){k}^{2}\frac{1}{{\tau}_{\text{rad}}{\tau}_{\text{drag}}}=0,$$(46)
that is, $$\omega =\frac{1}{2}\left(\text{i}\left(\frac{1}{{\tau}_{\text{rad}}}+\frac{1}{{\tau}_{\text{drag}}}\right)\pm \sqrt{4{k}^{2}{\left(\frac{1}{{\tau}_{\text{rad}}}\frac{1}{{\tau}_{\text{drag}}}\right)}^{2}}\right),$$(47)
where the term under the square root can be negative and, therefore, it can provide an imaginary component. Further algebraic manipulation then yields $$\frac{\text{i}k}{\text{i}\omega +\frac{1}{{\tau}_{\text{drag}}}}=\frac{2\text{i}k}{\left(\frac{1}{{\tau}_{\text{drag}}}\frac{1}{{\tau}_{\text{rad}}}\right)\pm \text{i}\sqrt{4{k}^{2}{\left(\frac{1}{{\tau}_{\text{rad}}}\frac{1}{{\tau}_{\text{drag}}}\right)}^{2}}}.$$(48)
In order to satisfy the boundary conditions, the term under the square root in this equation must be positive or the result is a pure imaginary number. In other words, Kelvin waves are able to propagate only when the condition $$4{k}^{2}>{\left(\frac{1}{{\tau}_{\text{rad}}}\frac{1}{{\tau}_{\text{drag}}}\right)}^{2}\text{\hspace{0.17em}}$$(49)
is met. Additionally, this simple estimation of the regimes where Kelvin waves can be supported in the atmosphere may well be an overestimate for the 3D, spherical coordinate case as the characteristic scale of the damping of the Kelvin wave must be smaller than the scale of the planet’s atmosphere itself. The real part of $\text{i}k/(\text{i}\omega +{\tau}_{\text{drag}}{}^{1})$ must therefore be large enough (and negative). Finally, when Kelvin waves propagate their characteristic decay rate (σ_{K}) is given by $${\sigma}_{\text{K}}=\frac{1}{2}\left(\frac{1}{{\tau}_{\text{rad}}}+\frac{1}{{\tau}_{\text{drag}}}\right).$$(50)
This result is similar to the decay timescale for Rossby waves (compare Eqs. (43) and (50)). τ_{drag} and τ_{rad} therefore have a symmetric contribution for Kelvin waves, as expected when considering Eqs. (17) and (19) for v = 0: they have a symmetric effect on u and h.
4.1.5 Decay timescale summary
We obtained expressions for the asymptotic values of the decay timescales for damped waves under the 2D shallowwater, βplane framework (see Sect. 2.1). In particular, for the case of Kelvin waves we obtained an analytical expression for the decay rate, Eq. (50). We also show that for the regime where the analytical calculations are valid, Rossby waves exhibit the same decay rate as found for Kelvin waves. For the more general case, aside from considerations of whether the waves can be supported by the atmosphere, we have two limits: (i) for τ_{drag} ~ τ_{rad} and τ_{drag} ≪ τ_{rad}, simply, ${\sigma}_{\text{R}}~{\sigma}_{\text{K}}~{\sigma}_{\text{g}}~{\tau}_{\text{drag}}^{1}$ within a factor of ~2; (ii) for τ_{drag} ≫ τ_{rad}, the Kelvin decay rate is the invert of the radiative timescale. The Rossby decay rate is the invert of the drag timescale if ${\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1}\lesssim 0.5$ and we obtained no semianalytical solution if ${\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1}\gtrsim 0.5$. Finally, the gravity waves decay rate becomes ${\sigma}_{\text{g}}~{\tau}_{\text{drag}}^{1}+{\tau}_{\text{rad}}^{1/(n+2)}$ for τ_{rad} ≪ 1 and ${\sigma}_{\text{g}}~{\tau}_{\text{rad}}^{1}/3$ for τ_{rad} ≫ 1.
Some of the numerical values we used to derive these asymptotic evaluations are reported in plain lines on Fig. 3. As our results have been derived for the 2D shallowwater, βplane system, for the case of Rossby waves in particular, the behavior of the waves may not be captured correctly. Therefore, we move on to verifying and extending our approach into 3D coordinates using ECLIPS3D.
4.2 Extension to 3D with ECLIPS3D
So far we have determined the characteristic frequencies and decay timescales for various atmospheric waves in our 2D framework, introduced in Sect. 2.1. In this section we extend our analysis to full 3D spherical coordinates using ECLIPS3D. In 3D, spherical coordinates, the dependency of the waves on the stratification, and the value of the drag and radiative timescales is difficult to predict theoretically. Therefore, we approach the problem numerically by using ECLIPS3D, studying the modes which propagate in a stratified atmosphere at rest, as is used for the initial condition when simulating hot Jupiter atmospheres. The background temperature–pressure profile is set to that of Iro et al. (2005), employing the polynomial fit of Heng et al. (2011). The pressure at the bottom of the atmosphere is set to 10 bars (the depth of the atmosphere is now 7 × 10^{6} m), capturing the dynamically active region of the atmosphere, driven by the forcing in the first phase of simulation but without detailing the innermost regions where the density is much higher. As before, we have varied the inner boundary condition totest its impact on our results, and find our conclusions to be robust to this choice. The selection of the modes of interest is performed by first excluding modes with unrealistic amplitudes at the pole or boundary, and then selecting modes with only one oscillation in longitude, that is, wavenumber 1, matching the heating function. Additionally, we restrict the selection to those with at most two nodes in the latitude direction which are the dominant modes (see discussion earlier in this section). This selection process is also described in Debras et al. (2019).
For the first study, we verified that all modes are supported when τ_{drag} = τ_{rad} = 10^{6} s is adopted throughout the atmosphere and are similar in form to when τ_{drag} = τ_{rad} = 0 s but exhibit a characteristic decay frequency of 10^{−6} s^{−1}. Although Rossby and gravity waves are supported with characteristic heights on the order of the height of the atmosphere itself, Kelvin waves seem to only be supported at the pressure scale height or smaller.
We subsequently applied ECLIPS3D with τ_{rad} prescribed as a function of pressure following Iro et al. (2005) and τ_{drag} set as a constant between 10^{5} and 10^{6} s. For this setup, we recover the usual (see, e.g., Thuburn et al. 2002) Rossby and gravity modes over different vertical wavelengths. Therefore, although our 2D analysis breaks for Rossby modes with large τ_{drag} (see Sect. 4.1.3), we recover them in the full 3D spherical coordinate treatment using ECLIPS3D, with, notably, their decay rate always comparable to the reciprocal of the drag timescale. However, in this setup, with a pressure dependent radiative timescale, we do not find Kelvin modes with atmosphericscale characteristic heights (we use 40 points in the z coordinates, therefore we are unable to resolve modes with H ≲ 10^{5} m). However, we do obtain Kelvin modes with smaller characteristic heights in the deepest, highestpressure region of the atmosphere where the radiative timescale is long, in agreement with our previous estimations (see Sect. 4.1.4, Eq. (49)). For this setup, we also obtain mixed Kelvingravity modes, which are concentrated at the equator, as well as mixed Rossbygravity modes, with the Rossby component dominating in the high latitudes and the gravity component dominating near the equator. In the case of the pure gravity modes, the resulting frequencies and decay rates from ECLIPS3D in 3D spherical coordinates are in good agreement with the 2D estimations (Sect. 4.1.2). However, for Kelvin modes, although the decay rates obtained from ECLIPS3D are in good agreement with our 2D analytical expressions (Sect. 4.1.4), the obtained frequencies are a little larger than our analytical analysis would suggest. Finally, for pure Rossby modes, for the range where our 2D analysis is valid, the decay timescales are again in good agreement between our 2D estimations (Sect. 4.1.3) and the numerical results in 3D but similarly to the case of the Kelvin modes, the frequencies are slightly underestimated.
From our ECLIPS3D outputs we have calculated the value of $\sqrt{{\sigma}^{2}+{\nu}^{2}}$ for all modes supported in the simulated atmosphere (keeping τ_{rad} prescribed following Iro et al. (2005) and τ_{drag} = 10^{6} s), as the amplitude of a given mode in the linear steady state of the atmosphere is inversely proportional to $\sqrt{{\sigma}^{2}+{\nu}^{2}}$ (Eq. (28)). These results show that the value of $\sqrt{{\sigma}^{2}+{\nu}^{2}}$ is an order of magnitude smaller for Rossby modes, compared with Kelvin or gravity modes. As discussed in Sect. 4.1.3 the frequencies of the modes are not significantly altered by the drag timescale, and Matsuno (1966) has shown that Rossby waves have frequencies ~ 10 times smaller than gravity waves and ~ 3 times smaller than Kelvin waves in this regime. Additionally, when τ_{drag} is large but τ_{rad} is modest, the decay rate of Kelvin waves is controlled by the radiative timescale, whereas the decay rate of the gravity waves and Rossby waves is conversely controlled by the drag timescale. Therefore, the value of $\sqrt{{\sigma}^{2}+{\nu}^{2}}$ will be bigger for Kelvin waves over that obtained for Rossby and gravity waves, for both of which this quantity will be on the order of the frequency, which is ten times smaller for Rossby waves compared to gravity waves. This indicates that the Rossby waves will propagate over greater timescales and lengths and dissipate, globally, more energy in the linear steady state (see Eq. (27)). The influence of the Rossby waves in the steady linear circulation regime will be dominant over Kelvin and gravity modes, explaining the qualitative shape of Fig. 1b.
As a summary, Fig. 3, shows the decay timescales for gravity and Rossby and Kelvin waves, obtained from Eqs. (40), (50) as a function of the drag timescale for τ_{rad}= 0.3 (in dimensional units, a few 10^{4} s), a characteristic value in the superrotating regions of HD 209458b, or as a function of the radiative timescale for τ_{drag} = 20 (about 10^{6} s), a value which allows for superrotation in the nonlinear limit. We also plot values extracted from numerical exploration with ECLIPS3D, although the radiative timescale is not constant in these numerical results and the characteristic height might differ.We recover the fact that Kelvin waves are more damped than other waves for the timescales used in this section,thought to be representative of hot Jupiter atmospheres, but not for all timescales (notably when τ_{drag} = τ_{rad}). Additionally, there are many regions of the parameter space where Kelvin waves don’t actually propagate, as evidenced in Sect. 4.1.4. This highlights the need to constrain the timescales to understand the spinup of superrotation and to understand the wave behavior across different timescales.
Figure 4 shows the pressure perturbations (color scale, total pressure minus initial pressure) and horizontal winds (vectors) as a function of longitude and latitude, for four Rossby modes. Two of the modes in Fig. 4 are from ECLIPS3D, resulting in 3D spherical coordinate calculations, and two are from the analytical studies (i.e., derived using equations in Appendix A; with Matsuno (1966) notations, they both have n =1 and l = 3), shown as the left and right columns, respectively. The locations in longitude are arbitrary as the initial state is axisymmetric and at rest. The Rossby modes shown in Fig. 4 from ECLIPS3D have been chosen such that the maximum amplitude was present in the deeper, highpressure regions where drag is dominant (top panel) in one case and for the other case, the amplitude was maximum in the upper, lowpressure part of the atmosphere, where the radiation timescale is shorter than that of the drag (bottom panel). Figure 4 shows that the ECLIPS3D results and those from our 2D analytical treatments are in good agreement. Specifically, the “tilt” of the modes in the latitude–longitude plane and the location of the maximum perturbation in pressure are broadly consistent between the analytical 2D and numerical 3D results. This agreement is comforting given that one case is a simplification of an atmosphere on a 2D shallow water βplane and the other is a restriction onto one height of a fully 3D, spherical mode. There are, however, discrepancies, notably at mid and high latitudes.
Interestingly, with ECLIPS3D we also recover Rossby waves with the opposite tilt in latitude than the results of the shallow water equations (right panel of Fig. 4), as well as Rossby waves with no tilt (their horizontal shape being similar to the Rossby waves from Matsuno 1966). All of these waves have comparable frequencies and decay rates, no matter the tilt. We attribute the existence of these waves to the density stratification and the dependence of the radiative timescale with height, not considered in the shallow water equations. Therefore, it is no easy task to predict what will be the shape of the linear steady state as it depends on the projection of the heating function on all these waves with different tilts, but also different characteristic heights. The fact that the waves are no longer orthogonal (Appendix A) further prevents an easy evaluation of the projection of the heating function on each wave. However, it is primordial to note that Rossby waves always exhibit a much lower amplitude in the pressure and temperature perturbation at the equator rather than in mid latitudes. This provides additional confidence in the assertion that Fig. 1b really is dominated by the Rossby wave component.
It is, therefore, difficult to draw conclusions about the behavior of the linear solutions solely based on the results of the ECLIPS3D calculations in 3D as we also require the projection of the heating function onto the waves. However, we clearly recover tilted Rossby waves and an absence of Kelvin waves in the upper part of the atmosphere where superrotation develops (or they must have a very small characteristic height). This is in contradiction to the findings of Showman & Polvani (2011), where they link the acceleration of superrotation to the interaction between a standing equatorial Kelvin wave and mid latitude Rossby wave; but it is in agreement with our analytical estimations for the domain of existence of Kelvin waves. As already stated in this paper, this does not refute the theories of Showman & Polvani (2011) and Tsai et al. (2014) regarding the equilibration of superrotation but, rather, it shows that the spinup of an initial jet is not due to a linear, chevronshaped steady state and that timedependent linear considerations must be taken into account. Globally, our 2D semianalytical arguments seem to be a good approximation of the 3D linear evolution of the atmosphere of hot Jupiters. We devote the next section to an application of these estimations to provide a physical understanding to the acceleration of superrotation.
Fig. 3 Typical decay rate σ for gravity and Rossby and Kelvin waves as a function of (top) τ_{drag} for τ_{rad}= 0.3 and (bottom) τ_{rad} for τ_{drag} = 20. The lines are values obtained using Eqs. (40), (50) while crosses are results from seven ECLIPS3D calculations.As τ_{rad} is not constant with depth in the ECLIPS3D results, we have chosen to use an arbitrary value of τ_{rad} = 0.3 for comparison. However, when τ_{drag} increasesthe location where the wave exhibits its maximum perturbation, it moves to higher pressures which should correspond to an increase in the equivalent τ_{rad}. For the Rossby waves, the low τ_{drag} limit has not been studied numerically as it is irrelevant for superrotation. Kelvin waves of comparable height with Rossby and gravity waves are only clearly identified in four ECLIPS3D calculations. 
Fig. 4 Pressure perturbations (colorscale) and winds (arrows) as a function of longitude and latitude for four Rossby waves for an axisymmetric hot Jupiter atmosphere at rest. Units are arbitrary. The values of the drag and radiative timescales are described in Sect. 4.2. Left column: results from ECLIPS3D 3D spherical coordinate calculations and right column: for the analytical results to the equation developed in Appendix A. For the top row, the drag is dominant, and the bottom row radiation is dominant. We note that as the initial state is at rest and axisymmetric, there is an uncontrolled phase in longitude, meaning the longitudinal location is abitrary. 
5 Transition to superrotation
In Sect. 3, we show that the linear steady state of our atmosphere is not significantly altered when moving to a more realistic heating profile taken from a 3D GCM simulation. Therefore, to gain further insight into the acceleration of the superrotation, we turn to timedependent, linear effects. This leads us to develop expressions for the timedependent linear solution to the problem in Sect. 4. In the current section, we use these solutions to understand the transition to superrotation in simulated hot Jupiter atmospheres. We first assess the validity of our own results by comparing them with the simulations of Komacek & Showman (2016) in Sect. 5.1. This is followed in Sect. 5.2 by an orderofmagnitude analysis, which allows us to conclude that we can not explain the acceleration to superrotation under the simplifications employed for the linear steady state (although such simplifications are wellsuited for studying the equilibration of the superrotation, see, e.g., Tsai et al. 2014). Finally, we test our interpretation against the results of 3D GCM simulations in Sect. 5.3, revealing that vertical accelerations are vital to the process.
5.1 Qualitative structure of linear steady state
Before discussing the transition to superrotation in the nonlinear limit, we first apply our findings from Sect. 4 to interpret the form of the various linear steady states presented in Komacek & Showman (2016), as seen in their Fig. 5. As shown by Wu et al. (2001), the zonal damping rate is proportional to $\sqrt{{\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1}}$. Therefore, if the two timescales (drag and radiative) are small, or one of them is vanishingly small, the zonal propagation of waves will be extremely limited in longitude. This is clearly seen in Komacek & Showman (2016) when one or two of the timescales are short and the temperature gradient is huge between the day and night side, while the zonal flows are restricted largely to the day side. Essentially, in this case, the waves excited on the day side have been damped before reaching the night side and, therefore, they cannot lead to efficient wind generation and heat redistribution (in the linear limit). This is also discussed in Komacek & Showman (2016).
Due to the strong asymmetric, steady forcing in the atmospheres of hot Jupiters, the dominant wavenumber k ensures one oscillation around the planet, that is, 2πrk = 2π. Typical conditions for hot Jupiter atmospheres yield gH ~ 4 × 10^{6} m^{2} s^{−2} (Showman & Polvani 2011) and therefore L ~ 7 × 10^{7} m and T ~ 3.5 × 10^{4} s (see Eqs. (6), (4)).
Let us suppose that one of the timescales (i.e., radiative or drag) is much shorter than the other, hence it is the dominant timescale, which we denote simply as τ. In order for Kelvin waves to propagate we require, from the dimensional form of Eq. (49): $$\tau \gtrsim \frac{r}{2\sqrt{gH}}~2.5\times {10}^{4}\text{\hspace{0.17em}s}.$$(51)
Therefore, if the dominant (shorter) timescale is inferior to ~ 2.5 × 10^{4} s, Kelvin waves cannot propagate unless both timescales are equal. However, if both these timescales are ≲ 2.5 × 10^{4} s, the dissipation time for Kelvin waves will be very short. Assuming a simple characteristic speed of waves to be $\sqrt{gH}$, which in our case is roughly 2 × 10^{3} m s^{−1}, the time for a wave to travel around the whole planet is ~2 × 10^{5} s. Therefore, for drag or radiative timescales of ≲2.5 × 10^{4} s, even in the cases where Kelvin waves exist they can not propagate around the whole planet and, thereby, generate the stationary chevronshaped MG pattern of Showman & Polvani (2011). The chevronshaped pattern of the linear steady state can, thus, only exist when both timescales are ≳10^{5} s and comparable in value. This is seen in Fig. 5 of Komacek & Showman (2016), where this pattern is clearly not ubiquitous, and this puts a strong restriction on the cases where the explanation of Showman & Polvani (2011) is applicable. In other words, acceleration of superrotation in hot Jupiters from the MG (or chevronshaped) tilted linear steady state is only possible over a restricted parameter space, which is not likely, therefore, to provide an explanation for all exoplanet cases.
Additionally, Eq. (51) shows that for a given τ but varying H, there is a minimum height for the propagation of Kelvin waves. As shown by Wu et al. (2000) and later by Tsai et al. (2014), the 3D structure of the propagating waves can be decomposed onto waves in a shallow water system but with differing equivalent depths. The modes are solutions to the homogeneous equations but featuring a characteristic height which varies between modes. Tsai et al. (2014) also show in their Fig. 2 that the projection of the heating function of the vertical modes has a high amplitude for modes with equivalent height between 5H_{P} and 0.2H_{P}, where H_{P} is the pressure scale height, roughly of the order of 4 × 10^{5} m in hot Jupiter atmospheres. Therefore, by adopting these limits, we have H = 5H_{P} ~ 2 × 10^{6} m and H = 0.2H_{P} ~ 8 × 10^{4} m, and obtain $${10}^{4}\text{s}<\tau <6\times {10}^{4}\text{s}.$$(52)
In our case, this implies that if τ < 10^{4} s, Kelvin waves are unable to propagate with characteristic height less than 5H_{P}: the linear steady state will have an almost null projection onto Kelvin waves. However, if τ > 6 × 10^{4} s, the majority of the Kelvin waves excited by the forcing can propagate (we recall that if both timescale are equal, all the wave can propagate, as in the neutral case, but we expect the radiative timescale to be at least an order of magnitude smaller in superrotating regions). Additionally, as introduced in Sect. 4, our estimates for the regimes where Kelvin waves are supported by the atmosphere is likely to be wider than the real situation, meaning that the criteria for Kelvin wave propagation are also likely to be stricter.
The behavior outlined in this section is readily apparent in Fig. 5 of Komacek & Showman (2016). In the limit where the drag is strong, the waves are damped efficiently, the thermal structure strongly resembles the thermal forcing, and there is no planetary Kelvin wave structure evident at the equator. However, in the case of weak drag (τ_{drag} > 10^{5} s), the Kelvin wave component is visible in the temperature structure only when the radiative timescale is comparable (i.e., > a few 10^{4} s, in other cases the temperature is almost uniform at the equator). Finally, in the limit of short radiative timescale, the Kelvin waves do not propagate, the dynamical shape of the atmosphere is dominated by other components, and the linear steady state strongly resembles Fig. 1b. Interestingly, it appears that the cases that superrotate in the nonlinear limit all have a negligible Kelvin wave contribution in their linear steady state. In other words, either the equator is not dominated by Kelvintype circulation or the high latitudes are dominant.
To conclude this section, we further detail the case of τ_{drag} = 10^{5} s and τ_{rad,top} = 10^{4} s. Following Eq. (9) of Komacek & Showman (2016), this choice for τ_{rad,top} implies a radiative timescale of 8 × 10^{4} s at P = 80 mbar (the isobaric surface presented in their Fig. 5) and the drag timescale 10^{5} s. Using (Eq. (49)), Kelvin waves are able to propagate in this scenario. Additionally, the difference between the drag and radiative timescales allows us to properly consider the behavior of the Rossby wave component. Solving Eq. (40) for these prescribed radiative and drag timescales yields a decay timescale of ~ 0.38 (in nondimensional units) for both Rossby and Kelvin waves, while the decay timescale of the gravity waves is ~ 0.36. Therefore, when τ_{drag} = 10^{5} s and τ_{rad,top} = 10^{4} s, we have Kelvin wave propagation, with a decay timescale long enough for the waves to traverse the entire planet and comparable lifetimes for all three wave types considered (Kelvin, Rossby, and gravity). In this instance we expect the steady state of the atmosphere to be a combination of planetary waves all with roughly the same magnitude (depending on the projection of the heating function), which leads to the chevronshaped pattern predicted by Showman & Polvani (2011) in the Matsuno–Gill circulation. This is confirmed by Fig. 5 of Komacek & Showman (2016) where, in the limit discussed, the linear steady state is similar to that shown in our Fig. 1a.
These comparisons of our estimations with results from this work and previous studies show that our semianalytical analysis is actually rather powerful in its explanation of the resulting linear steady state response of a hotJupiterlike atmosphere. The natural nextstep is to explore the implications for numerical simulations of a hot Jupiter atmosphere from the initial condition to the final superrotating state.
5.2 Orderofmagnitude analysis
In this section, we estimate the maximum forcing under which the consideration of a linear steady state is relevant, then we estimate thetimedependent wave response in the weak drag regime.
As we are performing a linear study, for constant τ_{rad} and τ_{drag}, the value of the maximum velocity is proportional to the amplitude of the forcing, which is well represented by the daysidenightside equilibrium temperature difference that we apply at the top of the atmosphere, ΔT_{eq,top}. Assuming that the radiative timescale as a function of pressure within the atmosphere is appropriately represented by the polynomial fit of Heng et al. (2011), adapted from Iro et al. (2005), the amplitude of the MG circulation will then only depend on ΔT_{eq,top} and the drag timescale.
Denoting the zonal velocity of the linear steady state as u_{MG}, the linear steady state can only be reached if the nonlinear terms can be neglected when u = u_{MG}, that is, once the linear steady state is formed. The nonlinear terms scale with the zonal advection, ${u}_{\text{MG}}\partial {u}_{\text{MG}}/\partial x~{u}_{\text{MG}}^{2}/L$, where L is a characteristic length. Whereas, the linear terms are of the order of u_{MG}∕τ_{drag}. Therefore, equating these two estimates provides a maximum zonal speed for which the nonlinear terms may be accurately neglected, u_{max}, where $${u}_{\text{max}}~\frac{L}{{\tau}_{\text{drag}}},$$(53)
and above which a linear steady state will not be reached by the 3D GCM. For the case of hot Jupiters, L ranges from half the planetary circumference in the MG case to the full circumference in the superrotating case, that is, L ~ πR. If we denote the maximum zonal, equatorial wind by u_{MG,1}, for the MG solution with ΔT_{eq,top} = 1 K, using the linear relationship of u_{MG} to the forcing we have $${u}_{\text{MG}}\approx {u}_{\text{MG},1}\times \left(\frac{\Delta {T}_{\text{eq},\text{top}}}{1\text{K}}\right)\text{\hspace{0.17em}},$$(54)
for any ΔT_{eq,top} at a constant τ_{drag}. Combining Eqs. (53) and (54) then yields a maximum forcing temperature difference value for which the linear steady state can be reached: $$\left(\frac{\Delta {T}_{\text{max}}}{1\text{\hspace{0.17em}}K}\right)\approx \frac{1}{{u}_{\text{MG},1}}\frac{L}{{\tau}_{\text{drag}}}.$$(55)
For equilibrium temperature contrast at the top of the atmosphere greater than the value in Eq. (55) nonlinear effects can no longer be neglected during the acceleration to the linear steady state, which would not be reached by a GCM. This has already been noted in Sect. 3.3.2 of Tsai et al. (2014), where they acknowledge that their analysis is strictly valid only in the strong or moderate damping scenario. As we see in Eq. (55), if τ_{drag} is too large, indicating a low damping scenario, the maximum forcing will be very small and the linear approximation becomes invalid for forcing amplitudes relevant to hot Jupiters. This analysis allows us to more rigorously define the weak, modest, and strong damping regimes we had mention in Sect. 2.
In the strong or moderate damping scenario, ΔT_{eq,top} ≲ ΔT_{max}, the atmosphere first reaches the linear steady solution and then the subsequent evolution is controlled by nonlinear acceleration. In this regime, the nonlinear evolution from the Matsuno–Gill linear steady state is given by $\partial {u}_{\text{MG}}/\partial t~{u}_{\text{MG}}^{2}/L$, where the ${u}_{\text{MG}}^{2}/L$ term comes from advection, then the characteristic time τ_{evol} for the atmosphere to significantly depart from the MG state is $${\tau}_{\text{evol}}=\frac{L}{{u}_{\text{MG}}}.$$(56)
This limit is that studied in Tsai et al. (2014), where the waves change the mean flow in a quasistatic way leading to the emergence and equilibration of superrotation.
In the low damping scenario however, the atmosphere never reaches the MG steady state, and nonlinear considerations must be taken into account when the characteristic speed exceeds u_{max} (hence Eq. (29) is irrelevant and only Eq. (27) and (28) can be used to understand the transition to superrotation). Let us suppose for example, that the atmosphere is composed of two waves: a slowly oscillating, slowly decaying Rossby wave, hence ν_{1}, σ_{1} ≪ 1, and a quickly oscillating, strongly damped Kelvin or gravity wave with ν_{2}, σ_{2} ≫ 1. Our analysis of Sect. 4 shows that when τ_{rad} ~ 10^{4} s and τ_{drag} ~ 10^{6} s, Rossby waves indeed have small frequency and decay rate, whereas Kelvin and gravity waves have orderofmagnitude higher frequencies and decay rates (provided they can propagate). In this simplified case, Eq. (28) is simplified to $$\begin{array}{ll}{X}_{\text{F}}=\hfill & \text{\hspace{0.17em}}\frac{{q}_{1}}{{\sigma}_{1}\text{i}{\nu}_{1}}{\tilde{X}}_{1}(x,y)\left(1{e}^{(\text{i}{\nu}_{1}{\sigma}_{1})(t)}\right)\hfill \\ \hfill & +\text{\hspace{0.17em}}\frac{{q}_{2}}{{\sigma}_{2}\text{i}{\nu}_{2}}{\tilde{X}}_{2}(x,y)\left(1{e}^{(\text{i}{\nu}_{2}{\sigma}_{2})(t)}\right),\hfill \end{array}$$(57)
where X_{F} is the timedependent solution vector. We know that in the linear steady state, assuming q_{1} ~ q_{2}, the Rossby wave component will hold much more power than the Kelvin or gravity wave as i ν_{1} − σ_{1}≪iν_{2} − σ_{2}. However, if we select a time t such that iν_{2} − σ_{2}t ≪ 1 and (necessarily) iν_{1} − σ_{1}t ≪ 1, Eq. (57) can be linearized to first order, yielding $${X}_{\text{F}}({\sigma}_{1}t,{\sigma}_{2}t\ll 1)\approx {q}_{1}{\tilde{X}}_{1}t+{q}_{2}{\tilde{X}}_{2}t.$$(58)
Therefore, in the limit of very early times in the evolution the two wave components in the timedependent solution of the atmosphere are comparable (provided q_{1} ~ q_{2}). The wave components remain comparable even when iν_{2} − σ_{2}t ~ 1 but iν_{1} − σ_{1}t ≪ 1, where the exponential term for the Kelvin or gravity wave is almost zero, but the linearisation holds for the Rossby wave, hence, $${X}_{\text{F}}({\sigma}_{1}t\ll 1,{\sigma}_{2}t\gtrsim 1)\approx {q}_{1}{\tilde{X}}_{1}t+\frac{{q}_{2}}{{\sigma}_{2}\text{i}{\nu}_{2}}{\tilde{X}}_{2}.$$(59)
Dividing the amplitude of our first, Rossby wave, in this linear timedependent state, α_{1}, by the amplitude of wave 2, α_{2} yields, $$\left\frac{{\alpha}_{1}}{{\alpha}_{2}}\right=\left\frac{{q}_{1}}{{q}_{2}}({\sigma}_{2}\text{i}{\nu}_{2})t\right~\frac{{q}_{1}}{{q}_{2}}.$$(60)
Provided that the projection of the heating function on the two wave components are similar, that is, q_{1} ~ q_{2}, the timedependent solution will exhibit comparable amplitudes for both waves before the (slowest) Rossby wave has grown much larger than the asymptotic amplitude of the Kelvin or gravity wave. However, the steady state will be dominated by the Rossby wave component. Thus, although the linear steady state strongly depends on τ_{rad} and τ_{drag}, in the limits of short timescales the structure of the atmosphere only depends on the projection of the heating function on the propagating waves. In Fig. 5, we show the modulus of $\left(1{e}^{(\text{i}\nu \sigma )(t)}\right)/(\sigma \text{i}\nu )$ for a Rossby, Kelvin and gravity waves obtained by ECLIPS3D (see Fig. 3) when τ_{drag} = 2 × 10^{5} s and for a Rossby and gravity waves when τ_{drag} = 10^{6} s. In both cases, after a few 10^{4} s (a day being ~9 × 10^{4} s), the amplitude of the different waves is of the same order of magnitude whereas Rossby wave contribution dominates at later times.
This analysis suggests that the linear steady state is not responsible for accelerating superrotation in hot Jupiter atmospheresfor the low drag limit and it allows us to resolve the problem explained in Sect. 2.2. As discussed in Sect. 2.2, superrotating atmospheres were found by Komacek & Showman (2016) despite structures implying negative convergence of momentum at the equator in the linear limit (a “reverse” MG structure), in contradiction with the mechanism of Showman & Polvani (2011) which invokes the linear MG steady state solution.
This problem persists for physically plausible choices on the equilibrium temperature contrast and the drag and radiative timescale, typically, ΔT_{eq,top} = 500 K, τ_{drag} ~ 5 = 10^{5} s and τ_{rad}= 5 × 10^{4} s which yields u_{max} ~ 200 m.s^{−1}. When solved numerically, we obtain u_{MG,1} ~ 5 m s^{−1} hence ΔT_{max} ~ 150 K: the linear steady state cannot be reached. Specifically, after one day of simulation, our analysis shows that the linear steady state will not have been reached; however, numerically, we already have u > u_{max}, hence nonlinear effects can no longer be neglected. Regarding the dissipation timescales for the waves (Sect. 4), after one day we have σ_{K} t, ν_{K}t ~ 1, σ_{g}t, ν_{g}t ~ 1 and σ_{R}t, ν_{R}t ≪ 1, where ν_{K}, ν_{g}, ν_{R} are the Kelvin, gravity, and Rossby waves frequencies, respectively, and σ_{K}, σ_{g} σ_{R} the Kelvin, gravity and Rossby dissipation rate, respectively. Our estimates of this section, Eq. (60), show that this leads to an equivalent contribution of Rossby, gravity and Kelvin waves in the circulation.
In summary, after 1 day of simulation: (i) the structure of the atmosphere exhibits similar contribution in Rossby, gravity and Kelvin waves, which is characteristic of the chevronshaped pattern of Fig. 1a. As shown by Showman & Polvani (2011), the eddies from such a circulation favour the meridional convergence of eastward momentum at the equator; (ii) in addition, simple orders of magnitude show that the nonlinear terms (hence the contribution from the eddies) can no longer be neglected as they are of the same order of magnitude as the linear terms.
Therefore, when the nonlinear terms become dominant, they lead to a net acceleration of the equator. Even though the eddy acceleration from the linear steady state (our Fig. 1b) would tend to decelerate the equator, the equator is accelerated nonlinearly after one day of simulation because of this chevronshaped, nonsteady linear circulation. This initiates superrotation and the later, slower evolution is well described by Tsai et al. (2014). We study this further in the next section using our own general circulation model (GCM) based on the unified model (UM, presented in Mayne et al. 2014b, 2017).
We must emphasize that this separation of scales between linear and nonlinear behavior is obviously very simplified. As already noted by Showman & Polvani (2011), the nonlinear accelerations are mostly due to the wavemean flow interactions (rather than wavewave or mean flowmean flow, as we confirm in the next section). Therefore, a quasilinear analysis or statistical studies of momentum transfer might allow for more rigorous insight into the jet acceleration (see, e.g., Srinivasan & Young 2012; Tobias & Marston 2013; Bakas & Ioannou 2013; Bouchet et al. 2013; Bakas et al. 2015; Herbert et al. 2019).
Fig. 5 Modulus of $\left(1{e}^{(\text{i}\nu \sigma )(t)}\right)/(\sigma \text{i}\nu )$ normalized by the final value for the Rossby wave as a function of time for gravity, Rossby and Kelvin waves with the numerical frequencies and decay rates obtained with ECLIPS3D. τ_{rad} follows the prescription of Iro et al. (2005) whereas τ_{drag} is set constant either to 2 × 10^{5} s (plain lines) or 10^{6} s (dashed lines). 
5.3 3D GCM simulations
In order to assess the applicability of the linear shallow water result developed in this work to a full 3D calculation, we performed simulations using the UM across a range of forcing scenarios. The simulations employ Newtonian heating as discussed in Mayne et al. (2014b) and adopt the baseline hot Jupiter setup presented in Mayne et al. (2014a) which follows that of Heng et al. (2011), and for the radiative timescale, Iro et al. (2005). For our simulations, we then varied the day to night temperature contrast, ΔT from 1 to 100 K (see Eqs. (B2) and (B3) of Heng et al. 2011). Obviously, this is only a toy model as we do not expect to find a tidallylocked planet with an effective temperatureof 1300 K and a day night contrast of 0.1 K, but this allows us to study the physical mechanisms controlling the atmospheric structure. Atmospheric drag has not been explicitly implemented in the UM but it is provided by a diffusion scheme as detailed in Mayne et al. (2014a). We verified that all of our simulations conserve mass and angular momentum to an order of 10^{−6}. We used these simulations to first explore the resulting qualitative structure of the atmosphere and then the accelerations within it, detailed in Sects. 5.3.1 and 5.3.2, respectively.
5.3.1 Qualitative structure of the atmosphere
As long as the linear effects are dominant, we expect all our simulations to have qualitatively similar features but with quantitativevalues that are scaled with the forcing. After one day of simulation, this is indeed the case, whereupon all of our simulations have the same qualitative structure matching Fig. 1a, although the magnitude of the temperature differences and wind velocities vary between simulations (increasing with larger temperature contrast). The structure matches the chevronshaped pattern of Showman & Polvani (2011) but, again, we state that this is not the steady MG solution. It is a specific time in the evolution of the atmosphere when all waves have the same projection in the circulation, as discussed in the previous section.
Comparing the very low temperature contrast case with linear steady states from ECLIPS3D, the dissipation within our simulations is equivalent to τ_{drag} ~ a few 10^{5} s. At P ~ 80 mbar the radiative timescale of Heng et al. (2011), adapted from Iro et al. (2005), is of the order of 2.5 × 10^{4} s. For these parameters the linear MG state can only be considered as reached after ~ 10 days, the time for the gravity and Rossby waves (which are the most long lived components) to be completely dissipated (see Sect. 4.1.5). Figure 6 shows the temperature and wind structure for three calculations, the first one using ECLIPS3D with parameters set to those matching the GCM simulations (Fig. 6a), and the next two from GCM simulations after 10 days of simulations, at 80 mbar, for a small and large temperature contrast at the top of the atmosphere ΔT_{eq,top} = 1 K (Fig. 6b) and ΔT_{eq,top} = 100 K, (Fig. 6c).
The ECLIPS3D calculation, hence the linear steady state, in Fig. 6a recovers the dominant midlatitude Rossby gyres associated with equatorial winds but it clearly differs from Fig. 1a in the sense that the equatorial circulation is not impacting the Rossby gyres significantly. We have, therefore, an intermediate between Figs. 1a and b. For completeness, we also show the evolution of the atmosphere for the same simulations but with τ_{drag} = 10^{6} s in Appendix D, which leads to a linear steady circulation equivalent to Fig. 1b. For the GCM simulations adopting ΔT_{eq,top} = 1 K, Fig. 6b, the circulation of Fig. 6a is recovered after 10 days of simulation. For the GCM simulations adopting ΔT_{eq,top} = 100 K in Fig. 6c, respectively, the longitudinal extent of the westward wind is reduced after 10 days compared to the weak forcing regime. Hence, for the simulation with the larger temperature contrast, after ten days the atmosphere has already diverged from the linear evolution of the atmosphere. Although both simulations were qualitatively identical after one day of simulation, the lowforcing simulations then reach the linear steady state, whereas higherforcing simulations never go through this state. Using the steady state wind velocities for the smaller temperature contrast simulation, ~ 2 m s^{−1} and Eq. (56), alongside the fact that we expect the linear steady wind to scale with the forcing, we can estimate the timescale to depart from the MG state in the stronger forcing simulation: ~ 10^{6} s, which is about 10 days. This timescale matches the timescale estimated above from the analysis of the atmospheric waves for convergence to the MG state; hence, the MG state is never actually reached.
In our simulations, the case when the linear steady state is not reached occurs for daynight temperature contrasts of ~ 100 K or greater at drag timescales of τ_{drag} ~ a few 10^{5} s (as presented in Fig. 6). If the drag timescale is further increased, the temperature contrast for which the linear steady state could not be reached would decrease further^{2} and, rather, HD 209458b is expected to experience a ~500 K temperaturecontrast. Hence, it does not seem possible to reconciliate the initial acceleration of superrotation with the consideration of steady linear effect. The linear considerations can only be used in the first day or so of simulation and the linear steady stateis irrelevant.
Therefore, if superrotation does exist in hot Jupiter atmospheres, the steady linear considerations are not likely to sufficiently explain the initial acceleration as nonlinear effects quickly dominate. Notably, the study of both Tsai et al. (2014) and Showman & Polvani (2011) (and more recently Hammond & Pierrehumbert 2018) would only apply to the limit of slow evolution, hence, once the atmosphere is already close to a nonlinear steady state. This explains why Fig. 16 of Tsai et al. (2014), which represents their linear consideration, compares so well with their Fig. 15 taken from 3D numerical simulations: when an initial superrotation is already settled, the further evolution is slow and can be understood in the linear limit. However, Tsai et al. (2014) provide no comparison between the linear expectation and the 3D simulations during the original acceleration phase of superrotation.
As we see in this section, linear considerations apply as long as we consider the timedependent solution. Although we cannot use our linear prescriptions to predict the evolution of the atmosphere in the nonlinear phase (as stressed in Sect. 5.2), we can estimate the duration of validity of the linear approximation and show that the atmosphere does not go through a linear steady state during the acceleration phase. Therefore, it is worth noting that the westward shift of the hot spot in the steady linear limit studied by Hindle et al. (2019), with the addition of a magnetic field, is not a robust enough diagnostic to predict whether the atmosphere is superrotating.
Interestingly, our low forcing simulation also converged to a superrotating state after a much longer time (scaling as one would expect as the inverse of the forcing). This was unexpected as Fig. 6b is not associated with a strong deposition of eastward momentum at the equator. More surprisingly, we also recover a superrotating jet for lowforcing simulations when we further increase the drag timescale and the atmosphere goes through a linear steady state resembling Fig. 1b. To understand this phenomenon, we conclude our study with the considerations of 3D accelerations in the spinup and equilibration of superrotation for GCM results.
Fig. 6 Temperature (color) and winds (arrows) as a function of longitude and latitude for (a) linear steady state calculated with ECLIPS3D with ΔT_{eq,top} = 100 K, τ_{drag} =2 × 10^{5} s and τ_{rad} following Iro et al. (2005). Maximum speed is 400 m s^{−1} (b): 3D GCM result at the 80 mbar level after 10 days of simulation, with ΔT_{eq,top} = 1 K and the radiative timescale of Iro et al. (2005). The maximum speed is 2.5 m s^{−1}; (c): 3D GCM result at the 80 mbar level after 10 days of simulation, with ΔT_{eq,top} = 100 K and the radiative timescale of Iro et al. (2005). The maximum speed is 350 m s^{−1} 
5.3.2 Accelerations
The final step in this process is to study the acceleration of the mean flow in the initial stages of the acceleration of superrotation within the 3D GCM simulations and assess the relevance of the 2D studies. As discussed, once the jet is settled and evolving slowly, the studies of Tsai et al. (2014) and Hammond & Pierrehumbert (2018) describe the evolution of superrotation but the initial acceleration is less clear, as we have explained in this paper. For this purpose, we study the acceleration of the jet in a simulation with ΔT_{eq,top} = 100 K and τ_{drag} ~ 10^{5} s, where the first phases of the development of superrotation can be captured over about 30 days. Following Mayne et al. (2017) (who adapted the treatment of Hardiman et al. 2010), the acceleration of the zonal mean flow can be written as $$\begin{array}{ll}{(\overline{\rho}\text{\hspace{0.17em}}\overline{u})}_{,t}=\hfill & \frac{{(\overline{\rho v}\text{\hspace{0.17em}}\overline{u}{\mathrm{cos}}^{2}\varphi )}_{,\varphi}}{r{\mathrm{cos}}^{2}\varphi}\frac{{(\overline{\rho w}\text{\hspace{0.17em}}\overline{u}{r}^{3})}_{,r}}{{r}^{3}}+2\Omega \overline{\rho v}\mathrm{sin}\varphi \hfill \\ \hfill & 2\Omega \overline{\rho w}\mathrm{cos}\varphi {(\overline{{\rho}^{\prime}{u}^{\prime}})}_{,t}\frac{{\left[\overline{{(\rho v)}^{\prime}{u}^{\prime}}{\mathrm{cos}}^{2}\varphi \right]}_{,\varphi}}{r{\mathrm{cos}}^{2}\varphi}\hfill \\ \hfill & \frac{{\left[\overline{{(\rho w)}^{\prime}{u}^{\prime}}{r}^{3}\right]}_{,r}}{{r}^{3}}+\overline{\rho {G}_{\lambda}},\hfill \end{array}$$(61)
where G_{λ} denotes the body forces acting in the longitudinal direction (not considered here), the subscripts denote partial derivatives and every quantity X is defined as $X=\overline{X}+{X}^{\prime}$ where a bar denotes an average on longitude. In this section, we do not consider the mean flowmean flow accelerations as they are negligible during the initial acceleration within our simulations. However, once the superrotating jet has formed, these accelerations should be taken into account as they balance the eddy accelerations and eventually lead to a nonlinear steady state (see, notably, the conclusions of Tsai et al. 2014).
Following Showman & Polvani (2011), the meridional eddy accelerations, involving v′ and u′, should lead to momentum convergence at the equator from the MG steady state, whereas the vertical component acts to decelerate the equatorial region. In Fig. 7, we show the value of $(\overline{\rho}\text{\hspace{0.17em}}\overline{u})$ for ΔT_{eq,top} = 100 K after 50 days of simulations, as well as the vertical and meridional accelerations. After 50 days, the jet extends from roughly 1 mbar to 1bar with the maximum of $(\overline{\rho}\text{\hspace{0.17em}}\overline{u})$ at around 0.2 bar. As in Fig. 15 of Tsai et al. (2014), we observe that the vertical accelerations are slowing down the upper part of the jet, while extending the jet to deeper pressures. The meridional accelerations on the other hand compensate the vertical component in agreement with both Tsai et al. (2014) and Showman & Polvani (2011). It is interesting to note that the explanation for radius inflation presented in Tremblin et al. (2017) relies on the vertical wind in the deep atmosphere due to the equatorial jet and that the spinup of the jet shows that the vertical accelerations are pushing the equatorial jet downwards. This seems to point towards a circulation in depth between the jet and the vertical velocities that gets deeper over time.
We also used these simulations to study the jet acceleration in more detail during the earlier phases. Firstly, we observed that the jet sets up initially between 0.08 and 0.1 bar in about 15 days (not shown) and then extends upwards and downwards. We show in Fig. 8 the meridional and vertical accelerations after 1, 5, and 20 days in the ΔT_{eq,top} = 100 K case.
Figure 8 shows that in the region where superrotation is the strongest (see Fig. 7), the vertical eddy acceleration always provides the maximum momentum convergence in the first 20 days, although the spatial extent of vertical acceleration decreases with time. Additionally, the location where the vertical motions accelerate the jet gets deeper with time, that is, it moves to higher pressures. We believe that this can be understood in the following way: the superrotation does not affect the midlatitude eddy circulation, which keeps acting to converge eastward momentum at the equator. However, as radiation penetrates deeper and the jet extends, the vertical circulation is changed and the vertical winds carry momentum away from the jet. This inhibits the deposition of eastward momentum at the equator, which is accompanied by a deceleration of westward winds on the night side of the planet. Globally, it appears that both meridional and vertical accelerations set up the initial superrotation (with the vertical accelerations being dominant), which then tends to decrease and even change the sign of the vertical acceleration. Subsequently, the global meridional motions act to sustain the jet once the vertical eddy acceleration is negative.
A key question is whether one reaches a limiting level of the daynight temperature contrast as a proxy for the radiative forcing at which the superrotation would transition. Such a transition would occur once the nonlinear terms become significant and it would depend on the state of the atmosphere at that point. If the atmosphere is similar to that of Fig. 1a, then superrotation would be favoured, but for a state such as 1b, the superrotation would be impeded. In our simulations, it seems that there is no threshold to superrotation because the vertical accelerations spin up the initial jet in all cases. The only limiting aspect is the time required to reach a nonlinear superrotating state as the forcing gets lower.
In the case of long drag and short radiative timescales, Komacek & Showman (2016) have already noted that the meridional motion of the linear steady state is opposite to what is necessary to drive superrotation, although they do observe that the nonlinear steady state is actually superrotating. We have explained this by considering the timedependent linear state in Sect. 5.3.1. More precisely, Fig. 9 shows the meridional and vertical accelerations as a function of time for two simulations adopting ΔT_{eq,top} = 1 K and ΔT_{eq,top} = 100 K and the same dissipation. Figure 10 shows the evolution of the zonally averaged zonal wind for comparison.
For both simulations, Figs. 9a and b, in the first two days meridional and vertical accelerations are both positive and of comparable magnitude. This was expected as over the course of the first few days, the atmosphere is comparable to the MG steady state explored by Showman & Polvani (2011). The zonallyaveraged zonal speed in Fig. 10 is almost zero, as expected from the propagation of waves with zonal wavenumber m = 1. In the lowforcing case, the vertical accelerations eventually dominate by almost one order of magnitude and the meridional terms can even end up opposing the creation of a jet (for a given pressure level, which is not seen on the vertically averaged Fig. 9). This is due to the fact that the atmosphere has reached the linear steady state of Fig. 6a, associated with strong vertical accelerations but with almost zero (or negative) meridional accelerations. For the highforcing case, although vertical accelerations contribute to the initiation of superrotation, they get smaller with time as the jet is being created, and eventually become negative (after ~60 days, not shown). This confirms that the vertical accelerations initiate the jet but then tend to extend it to deeper pressure, while meridional momentum convergence allows for an equilibrated state.
Globally, we can now resolve the discrepancy between Figs. 4 and 5 of Komacek & Showman (2016), which is discussed throughout this paper (notably Sect. 2.2 and this section): nonlinearly superrotating atmospheres can have linear steady states that seem to oppose the triggering superrotation, which is in contradiction with Showman & Polvani (2011). The study of Tsai et al. (2014) in the limit of strong dissipation does not offer explanation for the apparent paradox. Firstly, as we have shown,the treatment of the timedependent linear solution shows that the linear steady state is not relevant in the highforcing case. After one day, the nonlinear terms become dominant, whereas the MG steady state would require linear effects to dominate for at least 10 days. The shape of the atmosphere after one day is again given by our Fig. 1a: it is the usual chevronshaped pattern presented by Showman & Polvani (2011). Therefore, when nonlinear effects become dominant, they tend to accelerate the equator, whereas if the linear steady state had been reached, the deceleration by meridional motion could have been dominant: for adequate forcing in hot Jupiter conditions, the state of the atmosphere after one day always leads to meridional momentum convergence at the equator. Then, as seen in Figs. 8 and 9, the vertical accelerations also need to be taken into account: at the 80 mbar pressure range, vertical motion triggers the emergence of superrotation contrary to what is proposed by Showman & Polvani (2011). Only when superrotation is settled (Fig. 7b) do the vertical accelerations tend to decelerate the jet, and extend it deeper, that is, to higher pressures.
Later on, once the superrotation is settled, the study of Tsai et al. (2014) can be applied to the slower evolution of the atmosphere, which leads to the possible existence of a unique steady state. This steady state is permitted by both meridional and vertical accelerations, as we have seen throughout this work. This explains why superrotation is reached even when the dissipation is very low, although the initial acceleration is not included in the explanation given in Tsai et al. (2014).
Fig. 7 (a) $(\overline{\rho}\text{\hspace{0.17em}}\overline{u})$ (kg m^{−2} s^{−1}) as a functionof pressure and latitude after 50 days of simulation with ΔT_{eq,top} = 100 K. (b) Vertical eddy acceleration (kg m^{−2} s^{−2}) for the same simulation. (c) Meridional eddy acceleration (kg m^{−2} s^{−2}) for the samesimulation. 
Fig. 8 Left column: meridional eddy accelerations as a function of pressure and latitude in the ΔT_{eq,top} = 100 K case after 1 (top), 5 (middle) and 20 (bottom) days. Right column: same for vertical eddy acceleration. Units in kg m^{−2} s^{−2}. 
6 Conclusion
In this study, we explore the initial acceleration of superrotation in the context of a hot Jupiter atmosphere. We also focuse on the inherent discrepancy between the works of Showman & Polvani (2011) and Komacek & Showman (2016). Showman & Polvani (2011) propose that the superrotation is triggered by nonlinear accelerations around the linear steady state which converge momentum to the equator. On the other hand, Komacek & Showman (2016) show that certain configurations that exhibit superrotation are also associated with momentum divergence at the equator in the linear steady state limit.
In order to resolve this apparent contradiction, we studied the general form of the timedependent linear response of the atmosphere to a constant, asymmetrical heating. This response depends on the shape of the forcing, the global shape and frequency of the waves it generates, and the decay rate of these waves. Our first conclusion, based on the use of ECLIPS3D, is that changing the longitudinal form of the forcing is not of prime importance with regard to the qualitative understanding of superrotation, although it does affect the results in a quantitative sense. The use of a Newtonian cooling scheme with a wavenumber of one in longitude is, therefore, a reasonable approximation.
We have also obtained an equation for the frequency and decay rates of the propagating waves, as in Heng & Workman (2014). We could not solve this equation analytically for Rossby and gravity waves; so we have, therefore, estimated the asymptotic behavior of the waves numerically. For Kelvin waves, on the other hand, the analytical solution is obtained. The estimated decay rates are reported in Sect. 4.1.5 and Fig. 3.
From there, we explain in qualitative terms the structure of the linear solutions with different drag and radiative timescales, as presentedin Fig. 5 of Komacek & Showman (2016) and our Fig. 1. The zonal dependency had also already been estimated by Wu et al. (2001) by other means. A major result of this present work is that in the limit of short times (compared to the damping rates), the waves present in the decomposition of the heating function contribute almost equally to the timedependent linear solution. This tends to create a Matsuno–Gill like circulation (Fig. 1a) in the first day of the evolution of a hot Jupiter atmosphere in a GCM, although the actual linear steady state would be a “reverse” Matsuno–Gill, exhibiting eastward momentum divergence at the equator (Fig. 1b). With orderofmagnitude analysis, we concluded that in simulations representative of hot Jupiters, the linear steady state could not be reached, however, nonlinear terms were dominant after approximately one day, hence, when the atmosphere resembles that of the Matsuno–Gill circulation. As a consequence, the equator is accelerated even though the linear steady state would tend to decelerate the equator, resolving part of the discrepancy between Komacek & Showman (2016) and Showman & Polvani (2011).
Finally, we considered the nonlinear accelerations taking place during the spinup of superrotation from 3D GCM simulations with different contrasts in temperature between the day and night or in strengths of forcing to assess the importance of the vertical accelerations. Once the jet is formed, the vertical acceleration tends to decelerate the upper part of the jet while extending it to deeper pressure. The meridional accelerations oppose this deceleration and a steady state can be reached, as already shown by Tsai et al. (2014) and Showman & Polvani (2011). On the other hand, during the acceleration, the vertical component contributes equally to the meridional component to form an initial superrotation. This is in disagreement with Fig. 11 of Showman & Polvani (2011); however, the data for this figure are averaged across the upper atmosphere (above 30 mbar) where superrotation does not develop or it is weaker. Numerically, it seems that as a jet is initiated, the vertical circulation is altered, preventing the vertical deposition of eastward momentum at the equator, whereas the meridional circulation is roughly unchanged.
Overall, in this work, we study the acceleration of superrotation based on theoretical and semi analytical grounds. We have provide a complement to previous studies to provide a coherent understanding of the initial acceleration of the equator of hot Jupiters. Combined with the works of Showman & Polvani (2011); Tsai et al. (2014); Komacek & Showman (2016) and Hammond & Pierrehumbert (2018), a relatively complete picture of the initial phase of the atmospheric dynamics of simulated hot Jupiterscan now be drawn.
Our simulations suggest that there are regions of parameter space for which the linear steady state does not accelerate superrotation but the early spinup from the rest does so. This suggests that multiple long term nonlinear states might be possible (as explored by, e.g., Thrastarson & Cho 2010; Liu & Showman 2013) depending on initial conditions and it might be a track towards understanding unusual observations in the field, such those of Dang et al. (2018).
Fig. 9 Eddy accelerations at the equator at the 80 mbar pressure level (plain lines) or averaged from 40 mbar to 0.4 bar (dashed lines) as a function of time for (a): ΔT_{eq,top} = 1 K, (b): ΔT_{eq,top} = 100 K. 
Fig. 10 Zonallyaveraged zonal speed at the 80 mbar pressure level (plain lines) or averaged from 40 mbar to 0.4 bar (dashed lines) as a function of time for (a): ΔT_{eq,top} = 1 K, (b): ΔT_{eq,top} = 100 K. 
Acknowledgements
We wish to thank an anonymous referee for a very detailed report that greatly improved the argumentation of the paper. F.D. wishes to thank Antoine Venaille for his valuable insights on many questions regarding atmospheric dynamics. F.D. thanks the European Research Council (ERC) for funding under the H2020 research and innovation programme (grant agreement #740651 NewWorlds). Some of the calculations for this work were performed using Met Office software. Additionally, some of calculations used the Dirac Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC Dirac HPC Facility (www.Dirac.ac.uk). This equipment is funded by BIS National EInfrastructure capital grant ST/K000373/1 and STFC Dirac Operations grant ST/K0003259/1. Dirac is part of the National EInfrastructure. This research also made use of the ISCA High Performance Computing Service at the University of Exeter. N.J.M. is part funded by a Leverhulme Trust Research Project Grant and partly supported by a Science and Technology Facilities Council Consolidated Grant (ST/R000395/1). This work is also partly supported by the ERC grant 787361COBOM.
Appendix A Nonorthogonality in the case τ_{rad}≠τ_{drag}
In the case of τ_{rad} = τ_{drag}, Matsuno (1966) obtained a simple equation linking u_{n,l} and v_{n,l}, the horizontal speeds of the eigenvector (n, l), from Eqs. (7) and (9): $${u}_{n,l}=\frac{{\omega}_{n,l}y{v}_{n,l}+k\partial {v}_{n,l}/\partial y}{i({\omega}_{n,l}k)({\omega}_{n,l}+k)}.$$(A.1)
A similar equation is obtained for h_{n,l}. As ${v}_{n,l}\propto {\psi}_{n}(y)={H}_{n}(y){e}^{{y}^{2}/2}$, with H_{n} the nth Hermite polynomial, the use of the recurrence relations of the Hermite polynomials: $$\begin{array}{lll}\frac{\text{d}{H}_{n}(y)}{\text{d}y}\hfill & =\hfill & 2n{H}_{n1}(y),\hfill \\ {H}_{n+1}(y)\hfill & =\hfill & 2y{H}_{n}(y)2n{H}_{n1}(y),\hfill \end{array}$$
implies that the eigenvector (n,l) is simply written as: $$\left(\begin{array}{c}{v}_{n,l}\\ {u}_{n,l}\\ {h}_{n,l}\end{array}\right)=\left(\begin{array}{c}\text{i}({\omega}_{n,l}^{2}{k}^{2}){\psi}_{n}(y)\\ \frac{1}{2}({\omega}_{n,l}k){\psi}_{n+1}(y)+n({\omega}_{n,l}+k){\psi}_{n1}(y)\\ \frac{1}{2}({\omega}_{n,l}k){\psi}_{n+1}(y)n({\omega}_{n,l}+k){\psi}_{n1}(y)\end{array}\right).$$(A.4)
The orthogonality of the eigenvectors is easily proven, and the completeness is proved by use of the completeness of the Hermite polynomial functions.
In the case τ_{rad}≠τ_{drag}, Eq. (A.1) is slightly changed to get: $${u}_{n,l}=\frac{\left(\text{i}{\omega}_{n,l}+{\tau}_{\text{rad}}^{1}\right)y{v}_{n,l}+\text{i}k\partial {v}_{n,l}/\partial y}{{k}^{2}+\left(\text{i}{\omega}_{n,l}+{\tau}_{\text{rad}}^{1}\right)\left(\text{i}{\omega}_{n,l}+{\tau}_{\text{drag}}^{1}\right)},$$(A.5)
and we recall that ${v}_{n,l}\propto {H}_{n}({c}_{n,l}y){e}^{{c}_{n,l}^{2}{y}^{2}/2}$ where $${c}_{n,l}^{4}=\frac{\text{i}{\omega}_{n,l}+{\tau}_{\text{rad}}^{1}}{\text{i}{\omega}_{n,l}+{\tau}_{\text{drag}}^{1}}.$$(A.6)
Therefore, the expression of the eigenvector (n,l) is: $$\begin{array}{l}\left(\begin{array}{c}{v}_{n,l}\\ {u}_{n,l}\\ {h}_{n,l}\end{array}\right)=\\ \left(\begin{array}{c}\left({k}^{2}+\left(\text{i}{\omega}_{n,l}+{\tau}_{\text{rad}}^{1}\right)\left(\text{i}{\omega}_{n,l}+{\tau}_{\text{drag}}^{1}\right)\right){\psi}_{n}({c}_{n,l}y)\\ \frac{1}{2}\left(\frac{\text{i}{\omega}_{n,l}+{\tau}_{\text{rad}}^{1}}{{c}_{n,l}}\text{i}k{c}_{n,l}\right){\psi}_{n+1}({c}_{n,l}y)+n\left(\frac{\text{i}{\omega}_{n,l}+{\tau}_{\text{rad}}^{1}}{{c}_{n,l}}+\text{i}k{c}_{n,l}\right){\psi}_{n1}({c}_{n,l}y)\\ \frac{1}{2}\left({c}_{n,l}\left(\text{i}{\omega}_{n,l}+{\tau}_{\text{drag}}^{1}\right)\frac{\text{i}k}{{c}_{n,l}}\right){\psi}_{n+1}({c}_{n,l}y)n\left({c}_{n,l}\left(\text{i}{\omega}_{n,l}+{\tau}_{\text{drag}}^{1}\right)+\frac{\text{i}k}{{c}_{n,l}}\right){\psi}_{n1}({c}_{n,l}y)\end{array}\right).\end{array}$$(A.7)
Due to the dependency with c_{n,l} in the parabolic cylinder function, the eigenvectors of Eq. (A.7) are not orthogonal anymore (the calculation is cumbersome but not difficult). On the other hand, it is easily shown that the set of eigenvectors (n,l) are linearly independent (thanks to the Hermite polynomials). A rigorous proof would be needed to assess that they form a complete set, allowing for a projection of the heating function onto these eigenvectors. A graphical representation of these waves is provided in Appendix C.
Appendix B Solutions to Eq. (40)
B.1 The argument principle
The lefthand side of Eq. (39) may be confused with a second order polynomial but the dependency of c with ω actually leads to a polynomial of order 6. Based on Eqs. (39) and (34), we obtain an equation for X = c^{2}: $$\begin{array}{l}{X}^{6}+\left(\frac{2\text{i}(2n+1)\Delta \tau}{k}\right){X}^{5}\text{}+(3+2\text{i}k\Delta \tau ){X}^{4}\left(\frac{4\text{i}(2n+1)\Delta \tau}{k}\right){X}^{3}\\ \text{\hspace{1em}}+\left(3+\frac{8\text{i}\Delta {\tau}^{3}}{k}4\text{i}k\Delta \tau \right){X}^{2}+\left(\frac{2\text{i}(2n+1)\Delta \tau}{k}\right)X\\ \text{\hspace{1em}}+(1+2\text{i}k\Delta \tau )=0,\end{array}$$(B.1)
with $\Delta \tau =({\tau}_{\text{drag}}^{1}{\tau}_{\text{rad}}^{1})/2$. Determiningthe number of propagating waves consists, therefore, of determining the number of roots of Eq. (B.1) with a positive real part, as explained in Sect. 4. We make use of the argument principle: denoting P(X) the polynomial in X of Eq. (B.1), the number N of roots of P in a domain K is given by: $$N(\gamma )=\frac{1}{2\text{i}\pi}{\displaystyle {\oint}_{\gamma}\frac{{P}^{\prime}(z)}{P(z)}\text{d}z},$$(B.2)
where γ is a positively oriented contour encompassing K. Denoting C_{1∕2,r} the semicircle of radius r, centered on the origin and cut by the pure imaginary line (Fig. B.1), the number of zeros of P with positive real part is given by Eq. (B.2) with γ = C_{1∕2,r} in the limit r →∞.
The right hand side of Eq. (B.2) can be estimated by decomposing the integral on the pure imaginary line and on the circle of radius r. Namely: $$2\text{i}\pi N({C}_{1/2,r})={\displaystyle {\int}_{r}^{r}\frac{{P}^{\prime}(\text{i}t)}{P(\text{i}t)}}\text{id}t+{\displaystyle {\int}_{\pi /2}^{\pi /2}\frac{{P}^{\prime}(r\text{\hspace{0.17em}}\mathrm{exp}(\text{i}\theta ))}{P(r\text{\hspace{0.17em}}\mathrm{exp}(\text{i}\theta ))}}r\text{i}\mathrm{exp}(\text{i}\theta )\text{d}\theta .$$(B.3)
For large r, the second term of this expression can be calculated thanks to the asymptotical expansion of a polynomial of order 6: $$\frac{{P}^{\prime}(r\text{\hspace{0.17em}}\mathrm{exp}(\text{i}\theta ))}{P(r\text{\hspace{0.17em}}\mathrm{exp}(\text{i}\theta ))}=\frac{6}{r\text{\hspace{0.17em}}\mathrm{exp}(\text{i}\theta )}+\mathcal{O}\left(\frac{1}{{r}^{2}}\right).$$(B.4)
Hence the second term in the right hand side of Eq. (B.3) tends to 6i π when r →∞.
The first term deserves further consideration. If the image of the imaginary numbers by P does not cross a given half line $\mathbb{D}$ with origin at z = 0, said otherwise $P(\text{i}\mathbb{R})\subset \u2102\setminus \mathbb{D}$, we can define a holomorphicfunction Ln such that Ln(exp(z)) = exp(Ln(z)) = z and in that case: $$\begin{array}{ll}{\displaystyle {\int}_{r}^{r}\frac{{P}^{\prime}(\text{i}t)}{P(\text{i}t)}}\text{id}t\hfill & ={\displaystyle {\int}_{r}^{r}\frac{\text{d}}{\text{d}t}}\text{Ln}(P(\text{i}t))\text{d}t\hfill \\ \hfill & =\text{Ln}\left(\frac{P(\text{i}r)}{P(\text{i}r)}\right){\to}_{r\to \infty}\text{Ln}(1)=0,\hfill \end{array}$$(B.5)
where the limit comes from the fact that P is an even degree polynomial. In that case, our two expressions for the terms of the right hand side of Eq. (B.3) would yield: $$N({C}_{1/2,\infty})=3,$$(B.6)
confirming that Eq. (40) has exactly three roots of positive real part. We go on to show, therefore, that (except for n = 0) the image of the imaginary numbers by P is always included in $\u2102\setminus i{\mathbb{R}}^{}$ or $\u2102\setminus i{\mathbb{R}}^{+}$, where $i{\mathbb{R}}^{}$ and $i{\mathbb{R}}^{+}$ are the imaginary numbers with negative and positive imaginary part, respectively.
Fig. B.1 Graphical representation of the C_{1∕2,r} contour, encompassing the complex numbers with positive real part when r →∞. 
B.2 Image of the polynomial
Based on further calculation, it can be shown that the polynomial P from Eq. (B.1) may be written as: $$P(\text{i}Y)={(1+{Y}^{2})}^{2}({Y}^{2}2{\alpha}_{n}Y+1)+2\text{i}k\Delta \tau ({Y}^{4}2(2{\alpha}_{0}^{2}1){Y}^{2}+1),$$(B.7)
for $Y\in \mathbb{R}$ and with α_{n} = (2n + 1)Δτ∕k. Hence we can write the imaginary part of P(iY) as: $$\Im (P(\text{i}Y))=2k\Delta \tau ({Y}^{2}2{\alpha}_{0}Y+1)({Y}^{2}+2{\alpha}_{0}Y+1).$$(B.8)
We have two cases: (i) if α_{0} = Δτ∕k < 1, the imaginary part of P(iY) has no roots for $Y\in \mathbb{R}$. Hence for all n, $P(\text{i}\mathbb{R})\subset \u2102\setminus i{\mathbb{R}}^{}$ or $P(\text{i}\mathbb{R})\subset \u2102\setminus i{\mathbb{R}}^{+}$ depending on the sign of Δτ, and we have the result: for all n, there are exactly three roots of positive real part of Eq. (40), hence only three waves propagate. (ii) For α_{0}  = Δτ∕k > 1, the real part of the polynomial is: $$\Re (P(\text{i}Y))={(1+{Y}^{2})}^{2}({Y}^{2}2{\alpha}_{n}Y+1),$$(B.9)
which has two real roots: $${y}_{n}^{\pm}={\alpha}_{n}\pm \sqrt{{\alpha}_{n}^{2}1}.$$(B.10)
Calculating the imaginary part of P for $Y={y}_{n}^{\pm}$ yields: $$\Im (P(\text{i}{y}_{n}^{\pm}))=(2k\Delta \tau )({\alpha}_{n}^{2}{\alpha}_{0}^{2})(4{y}_{n}^{\pm}{}^{2}).$$(B.11)
Hence, for n > 0 the imaginary part of $P(\text{i}{y}_{n}^{\pm})$ is always of the same sign as Δτ: for all y, P(i y) cannot be an imaginary number with a negative (resp. positive) imaginary part if Δτ is positive (resp. negative). Therefore $P(\text{i}\mathbb{R})\subset \u2102\setminus i{\mathbb{R}}^{}$ (resp. $P(\text{i}\mathbb{R})\subset \u2102\setminus i{\mathbb{R}}^{+}$) and we have the same result: for n > 0, there are exactly three roots of positive real part of Eq. (40).
We confirm, therefore, that for Δτ∕k < 1, only three waves propagate for all n and this result holds for Δτ∕k > 1 and n > 0. The case of Δτ∕k > 1 and n = 0 is more complicated as when the real part of P(iy) cancels out its imaginary part as well and we can’t apply the same argument as for n >0. However, it means that the polynomial P for n = 0, P_{0}, can be easily factorized: $${P}_{0}(\text{i}Y)=({Y}^{2}2{\alpha}_{0}Y+1)\{{(1+{Y}^{2})}^{2}+2\text{i}k\Delta \tau ({Y}^{2}+2{\alpha}_{0}Y+1)\},$$(B.12)
which yields: $$\begin{array}{ll}{P}_{0}(Y)=\hfill & \text{\hspace{0.17em}}({Y}^{2}+2\text{i}{\alpha}_{0}Y+1)\hfill \\ \hfill & \times \{{(1{Y}^{2})}^{2}+2\text{i}k\Delta \tau ({Y}^{2}2\text{i}{\alpha}_{0}Y+1)\}\hfill \\ \equiv \hfill & \text{\hspace{0.17em}}({Y}^{2}+2\text{i}{\alpha}_{0}Y+1)Q(Y).\hfill \end{array}$$(B.13)
Looking for the roots of P_{0} with a positive real part is, therefore, equivalent to looking for the roots of Q with positive real part. If we can show that $Q(\text{i}\mathbb{R})\subset \u2102\setminus i{\mathbb{R}}^{}$ or $Q(\text{i}\mathbb{R})\subset \u2102\setminus i{\mathbb{R}}^{+}$, then we can use the argument principle on Q, a polynomial of order 4, and Eq. (B.3) ensures that P has only two roots with positive real part. This is easily proven by looking at the real part of Q: $$\Re (Q(\text{i}Y))={(1+{Y}^{2})}^{2}>0,$$(B.14)
and, hence, for all $y\in \mathbb{R},$ the real partof Q(iy) never cancels out, so that $Q(\text{i}\mathbb{R})\subset \u2102\setminus i\mathbb{R}$: P has only two roots of positive real part and only two waves can propagate.
It might seem surprising to obtain a different behavior for n = 0 but this was already found by Matsuno (1966), where only two of the three solutions of the n = 0 case were actual solutions of the linearized equations of motion. When considering that τ_{drag}≠τ_{rad}, this degeneracy is removed when Δτ∕k > 1, where only two roots of the polynomial have a positive real part. These findings have been tested and confirmed numerically.
Appendix C Waves when τ_{rad}, τ_{drag} ≠ 0
This appendix aims to demonstrate the change in the structure of the waves when τ_{rad}, τ_{drag}≠0 based on the numerical solutions of Eqs. (40) and (A.7). Figure C.1 shows a Rossby, gravity, and Kelvin wave when τ_{rad} = τ_{drag} = 0.35, which correspondsto ~10^{5} s in dimensional units. Their shape is similar to the waves studied by Matsuno (1966) as the only difference is in the decay timescale, which is not zero but equal to τ_{drag}.
Fig. C.1 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = τ_{drag} = 2.8 which is ~ 10^{5} s in dimensional units. Presenting (a) westward propagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 
Fig. C.2 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = 2.8 and τ_{drag} = 28. Presenting (a) westward propagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 
Subsequently, Figs. C.2 and C.3 show the same waves when either the drag or radiative timescales has been multiplied by ten, respectively. The shape of the waves is almost unaltered for gravity and Kelvin waves, apart from a tilt of the perturbation with latitude. Rossby waves, on the other hand, are more affected by the change in the drag and radiative timescale.
Fig. C.3 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = 28 and τ_{drag} = 2.8. Presenting (a) westward propagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 
Fig. C.4 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = 1 (~ 3.5 × 10^{4} s in dimensional units) and τ_{drag} = 28 (~ 10^{6} s in dimensional units). Presenting: (a) westwardpropagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 
Finally,Fig. C.4 shows the same waves when τ_{rad}= 1 and τ_{drag} = 28, which is about 3.5 × 10^{4} s and 10^{6} s, respectively.The tilt in the gravity and Kelvin modes is amplified in comparison to Fig. C.2 and the shape of the Rossby wave is extremely altered. Because the majority of the energy of the Rossby wave is in very high latitudes, where the equatorial βplane framework breaks, we excluded these modes in our semianalytical treatment. Physical Rossby waves are, nonetheless, recovered in 3D by ECLIPS3D with such timescales. For more, see Sect. 4.2.
Appendix D Numerical development of the linear steady state
In this appendix, we show the development of the linear steady state for a simulation with ΔT_{eq,top} = 1 K, τ_{drag} = 10^{6} s and τ_{rad} following the prescription of Iro et al. (2005) and the same simulation with ΔT_{eq,top} = 100 K. Our estimates show that such a drag timescale should lead to a linear steady circulation similar to Fig. 1b after ~ 20 days of evolution in the lowforcing case; whereas the time to depart from the linear evolution in the highly forcing case should be about approximately one day. This is clearly recovered in Fig. D.1.
Figure D.1 shows the temperature and winds of both simulations after one, three, and 50 days of evolution. We see that the lowforcing simulation resembles Fig. 1b after 50 days while the highly forced simulation is superrotating.
After one day of simulation, a MGlike circulation is recovered in both cases. The maximum speed in the highforcing case is a hundred times the maximum speed of the lowforcing case. After just three days, the maximum speed of the high forcing case is more than 100 times that of the lowforced case. After 50 days, there is a factor of 200 in the difference between the two cases, highlighting the influence of nonlinear terms.
In the lowforcing case, the propagation and dissipation of midlatitude Rossby waves shifts the circulation obtained after one day towards a reverseMG state after 50 days, whereas in the highlyforced case the eddies from the circulation after one day accelerate the equator towards a superrotating state. Notably, we see on the middle panel of Fig. D.1 that the winds in the lowforced simulation are slowly shifted westward by the midlatitude Rossby waves, whereas the highly forced simulation leads to a shrink of the westward winds that disappear after ~ 30 days.
Fig. D.1 Temperature (colors) and winds at a height corresponding to 80 mbar pressure at t = 0. τ_{drag} =10^{6} s and τ_{rad} follow the prescription of Iro et al. (2005) for all simulations. The left column has ΔT_{eq,top} = 1 K. while the right column has ΔT_{eq,top} = 100 K. Finally, the top panel shows the atmosphere after one day of evolution, middle panel after three days, andbottom panel after 50 days. Maximum speeds are (a) 1 m s^{−1} ; (b) 100 m s^{−1}; (c) 2 m s^{−1} ; (d) 300 m s^{−1}; (e) 5 m s^{−1} and (f) 1000 m s^{−1}. 
References
 Abramowitz, M., & Stegun, I. A. 1965, Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables (New York: Dover Publications), 55 [Google Scholar]
 Amundsen, D., Baraffe, I., Tremblin, P., et al. 2014, A&A, 564, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Amundsen, D. S., Mayne, N. J., Baraffe, I., et al. 2016, A&A, 595, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Armstrong, D. J., de Mooij, E., Barstow, J., et al. 2016, Nat. Astron., 1, 0004 [NASA ADS] [CrossRef] [Google Scholar]
 Bakas, N. A., & Ioannou, P. J. 2013, J. Atm. Sci., 70, 2251 [NASA ADS] [CrossRef] [Google Scholar]
 Bakas, N. A., Constantinou, N. C., & Ioannou, P. J. 2015, J. Atm. Sci., 72, 1689 [NASA ADS] [CrossRef] [Google Scholar]
 Baraffe, I., Chabrier, G., & Barman, T. 2010, Rep. Prog. Phys., 73, 16901 [Google Scholar]
 Bouchet, F., Nardini, C., & Tangarife, T. 2013, J. Stat. Phys., 153, 572 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Cooper, C. S., & Showman, A. P. 2005, ApJ, 629, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Dang, L., Cowan, N. B., Schwartz, J. C., et al. 2018, Nat. Astron., 2, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Debras, F., Mayne, N., Baraffe, I., Goffrey, T., & Thuburn, J. 2019, A&A, 631, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DobbsDixon, I., & Agol, E. 2013, MNRAS, 435, 3159 [NASA ADS] [CrossRef] [Google Scholar]
 Drummond, B., Tremblin, P., Baraffe, I., et al. 2016, A&A, 594, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Drummond, B., Mayne, N. J., Manners, J., et al. 2018a, ApJ, 869, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Drummond, B., Mayne, N. J., Manners, J., et al. 2018b, ApJ, 855, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Drummond, B., Mayne, N. J., Baraffe, I., et al. 2018c, A&A, 612, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gill, A. E. 1980, Q. J. Roy. Meteorol. Soc, 106, 447 [Google Scholar]
 Hammond, M., & Pierrehumbert, R. T. 2018, ApJ, 869, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Hardiman, S. C., Andrews, D. G., White, A. A., Butchart, N., & Edmond, I. 2010, J. Atm. Sci., 67, 1983 [NASA ADS] [CrossRef] [Google Scholar]
 Helling, C., Lee, G., DobbsDixon, I., et al. 2016, MNRAS, 460, 855 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., & Workman, J. 2014, ApJS, 213, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Menou, K., & Phillipps, P. 2011, MNRAS, 413, 2380 [NASA ADS] [CrossRef] [Google Scholar]
 Herbert, C., Caballero, R., & Bouchet, F. 2019, J. Atm. Sci., 77, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Hindle, A. W., Bushby, P. J., & Rogers, T. M. 2019, ApJ, 872, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Iro, N., Bézard, B., & Guillot, T. 2005, A&A, 436, 719 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Komacek, T., & Showman, A. 2016, ApJ, 821, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, G., DobbsDixon, I., Helling, C., Bognar, K., & Woitke, P. 2016, A&A, 594, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lines, S., Mayne, N. J., Boutle, I. A., et al. 2018a, A&A, 615, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lines, S., Manners, J., Mayne, N. J., et al. 2018b, MNRAS, 481, 194 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, B., & Showman, A. P. 2013, ApJ, 770, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuno, T. 1966, J. Meteorol. Soc. Jpn, 44, 25 [Google Scholar]
 Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014a, A&A, 561, A1 [Google Scholar]
 Mayne, N., Baraffe, I., Acreman, D., et al. 2014b, Geosci. Model Dev., 7, 3059 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mayne, N. J., Debras, F., Baraffe, I., et al. 2017, A&A, 604, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
 PerezBecker, D., & Showman, A. P. 2013, ApJ, 776, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Rauscher, E., & Menou, K. 2012, ApJ, 750, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Roman, M., & Rauscher, E. 2019, ApJ, 872, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A., Cooper, C., Fortney, J., & Marley, M. 2008, ApJ, 682, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Sing, D. K., VidalMadjar, A., Lecavelier des Etangs, A., et al. 2008, ApJ, 686, 667 [NASA ADS] [CrossRef] [Google Scholar]
 Snellen, I. A. G., Albrecht, S., de Mooij, E. J. W., & Le Poole, R. S. 2008, A&A, 487, 357 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Srinivasan, K., & Young, W. R. 2012, J. Atm. Sci., 69, 1633 [NASA ADS] [CrossRef] [Google Scholar]
 Thrastarson, H. T., & Cho, J. Y.K. 2010, ApJ, 716, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Thuburn, J., Wood, N., & Staniforth, A. 2002, Q. J. Roy. Meteorol. Soc., 128, 1771 [NASA ADS] [CrossRef] [Google Scholar]
 Tobias, S. M., & Marston, J. B. 2013, Phys. Rev. Lett., 110, 104502 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Chabrier, G., Mayne, N. J., et al. 2017, ApJ, 841, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Tsai, S.M., DobbsDixon, I., & Gu, P.G. 2014, ApJ, 793, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Tsai, S.M., Kitzmann, D., Lyons, J. R., et al. 2018, ApJ, 862, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Z., Sarachik, E. S., & Battisti, D. S. 2000, J. Atm. Sci., 57, 2169 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Z., Sarachik, E. S., & Battisti, D. S. 2001, J. Atm. Sci., 58, 724 [NASA ADS] [CrossRef] [Google Scholar]
 Zellem, R. T., Lewis, N. K., Knutson, H. A., et al. 2014, ApJ, 790, 53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Temperature (colorscale in K) and horizontal wind (arrows) as a function of longitude (x axis) and latitude (y axis) at the 40 mbar pressure level of the linear steady state (denominated MG circulation) obtained using ECLIPS3D (Debras et al. 2019) with heating function, drag and radiative timescales following the definitions of Komacek & Showman (2016). Following the notation of Komacek & Showman (2016): (a) ΔT_{eq,top}= 100 K, τ_{drag,top} = 10^{5} s and τ_{rad,top} = 10^{5} s. The maximum speed at this pressure range is 10 m s^{−1}. (b) ΔT_{eq,top} = 100 K, τ_{drag,top} = 10^{6} s and τ_{rad,top} = 10^{4} s. The maximum speed at this pressure range is 100 m s^{−1}. Note that the maximum speed has been multiplied by ten as the drag timescale has been multiplied by ten. 

In the text 
Fig. 2 Temperature (colorscale in K) and winds (arrows) as a function of longitude and latitude with different heating functions, with drag and radiative timescales defined as in Komacek & Showman (2016) with nondimensional values of 1 (almost half an Earth day). (a–c): the initial pressure at this height is 50 mbar which corresponds to an initial temperature of 1326 K and the forcing is ΔT_{eq,top} = 50 K. (d): the initial pressure is 40 mbar and initial temperature 1322 K. (a): heating as in Komacek & Showman (2016) (b): same as in Komacek & Showman (2016) but with a cooling at the night side three times more efficient than the heating on the day side. (c): same as in Komacek & Showman (2016) with no cooling at the night side. (d): heating rate extracted from the full radiative transfer calculations of the GCM (divided by ten to obtain comparable values). 

In the text 
Fig. 3 Typical decay rate σ for gravity and Rossby and Kelvin waves as a function of (top) τ_{drag} for τ_{rad}= 0.3 and (bottom) τ_{rad} for τ_{drag} = 20. The lines are values obtained using Eqs. (40), (50) while crosses are results from seven ECLIPS3D calculations.As τ_{rad} is not constant with depth in the ECLIPS3D results, we have chosen to use an arbitrary value of τ_{rad} = 0.3 for comparison. However, when τ_{drag} increasesthe location where the wave exhibits its maximum perturbation, it moves to higher pressures which should correspond to an increase in the equivalent τ_{rad}. For the Rossby waves, the low τ_{drag} limit has not been studied numerically as it is irrelevant for superrotation. Kelvin waves of comparable height with Rossby and gravity waves are only clearly identified in four ECLIPS3D calculations. 

In the text 
Fig. 4 Pressure perturbations (colorscale) and winds (arrows) as a function of longitude and latitude for four Rossby waves for an axisymmetric hot Jupiter atmosphere at rest. Units are arbitrary. The values of the drag and radiative timescales are described in Sect. 4.2. Left column: results from ECLIPS3D 3D spherical coordinate calculations and right column: for the analytical results to the equation developed in Appendix A. For the top row, the drag is dominant, and the bottom row radiation is dominant. We note that as the initial state is at rest and axisymmetric, there is an uncontrolled phase in longitude, meaning the longitudinal location is abitrary. 

In the text 
Fig. 5 Modulus of $\left(1{e}^{(\text{i}\nu \sigma )(t)}\right)/(\sigma \text{i}\nu )$ normalized by the final value for the Rossby wave as a function of time for gravity, Rossby and Kelvin waves with the numerical frequencies and decay rates obtained with ECLIPS3D. τ_{rad} follows the prescription of Iro et al. (2005) whereas τ_{drag} is set constant either to 2 × 10^{5} s (plain lines) or 10^{6} s (dashed lines). 

In the text 
Fig. 6 Temperature (color) and winds (arrows) as a function of longitude and latitude for (a) linear steady state calculated with ECLIPS3D with ΔT_{eq,top} = 100 K, τ_{drag} =2 × 10^{5} s and τ_{rad} following Iro et al. (2005). Maximum speed is 400 m s^{−1} (b): 3D GCM result at the 80 mbar level after 10 days of simulation, with ΔT_{eq,top} = 1 K and the radiative timescale of Iro et al. (2005). The maximum speed is 2.5 m s^{−1}; (c): 3D GCM result at the 80 mbar level after 10 days of simulation, with ΔT_{eq,top} = 100 K and the radiative timescale of Iro et al. (2005). The maximum speed is 350 m s^{−1} 

In the text 
Fig. 7 (a) $(\overline{\rho}\text{\hspace{0.17em}}\overline{u})$ (kg m^{−2} s^{−1}) as a functionof pressure and latitude after 50 days of simulation with ΔT_{eq,top} = 100 K. (b) Vertical eddy acceleration (kg m^{−2} s^{−2}) for the same simulation. (c) Meridional eddy acceleration (kg m^{−2} s^{−2}) for the samesimulation. 

In the text 
Fig. 8 Left column: meridional eddy accelerations as a function of pressure and latitude in the ΔT_{eq,top} = 100 K case after 1 (top), 5 (middle) and 20 (bottom) days. Right column: same for vertical eddy acceleration. Units in kg m^{−2} s^{−2}. 

In the text 
Fig. 9 Eddy accelerations at the equator at the 80 mbar pressure level (plain lines) or averaged from 40 mbar to 0.4 bar (dashed lines) as a function of time for (a): ΔT_{eq,top} = 1 K, (b): ΔT_{eq,top} = 100 K. 

In the text 
Fig. 10 Zonallyaveraged zonal speed at the 80 mbar pressure level (plain lines) or averaged from 40 mbar to 0.4 bar (dashed lines) as a function of time for (a): ΔT_{eq,top} = 1 K, (b): ΔT_{eq,top} = 100 K. 

In the text 
Fig. B.1 Graphical representation of the C_{1∕2,r} contour, encompassing the complex numbers with positive real part when r →∞. 

In the text 
Fig. C.1 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = τ_{drag} = 2.8 which is ~ 10^{5} s in dimensional units. Presenting (a) westward propagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 

In the text 
Fig. C.2 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = 2.8 and τ_{drag} = 28. Presenting (a) westward propagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 

In the text 
Fig. C.3 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = 28 and τ_{drag} = 2.8. Presenting (a) westward propagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 

In the text 
Fig. C.4 Temperature (colors) and winds (arrows), both in arbitrary unit, as a function of longitude and latitude for the shallow water solutions of Eqs. (40) and (A.7) considering τ_{rad} = 1 (~ 3.5 × 10^{4} s in dimensional units) and τ_{drag} = 28 (~ 10^{6} s in dimensional units). Presenting: (a) westwardpropagating gravity wave; (b) Rossby wave; and (c) Kelvin wave. 

In the text 
Fig. D.1 Temperature (colors) and winds at a height corresponding to 80 mbar pressure at t = 0. τ_{drag} =10^{6} s and τ_{rad} follow the prescription of Iro et al. (2005) for all simulations. The left column has ΔT_{eq,top} = 1 K. while the right column has ΔT_{eq,top} = 100 K. Finally, the top panel shows the atmosphere after one day of evolution, middle panel after three days, andbottom panel after 50 days. Maximum speeds are (a) 1 m s^{−1} ; (b) 100 m s^{−1}; (c) 2 m s^{−1} ; (d) 300 m s^{−1}; (e) 5 m s^{−1} and (f) 1000 m s^{−1}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.