EQUATORIAL SUPERROTATION ON TIDALLY LOCKED EXOPLANETS

The past few years have witnessed major strides in our efforts to understand the atmospheric circulation of short-period exoplanets—both gas giants (hot Jupiters) and smaller terrestrial planets. Infrared photometry, spectra, and light curves from the Spitzer Space Telescope and Hubble Space Telescope now provide constraints on the three-dimensional (3D) temperature structure of several hot Jupiters, which hint at a vigorous atmospheric circulation on these bodies (e.g., Knutson et al. 2007, 2009; Charbonneau et al. 2008; Harrington et al. 2006; Cowan et al. 2007; Swain et al. 2009; Crossfield et al. 2010). These observations have motivated a growing effort to model the atmospheric circulation on these objects: to date, many 3D atmospheric circulation models of hot Jupiters have been published (Showman & Guillot 2002; Cooper & Showman 2005, 2006; Showman et al. 2008, 2009; Dobbs-Dixon & Lin 2008; Menou & Rauscher 2009; Rauscher & Menou 2010; Dobbs-Dixon et al. 2010; Thrastarson & Cho 2010; Lewis et al. 2010; Perna et al. 2010; Heng et al. 2010). These models have emphasized synchronously rotating hot Jupiters in circular, approximately 2–5 day orbits.

Just as the last decade witnessed the first characterization of hot Jupiters, the next decade will see a shift toward characterizing "super Earths" (planets of 1–10 Earth masses) and terrestrial planets. To date, roughly 30 super Earths have been discovered, including several that transit their host stars (Charbonneau et al. 2009; Léger et al. 2009; Batalha et al. 2011) with hundreds of additional candidates recently announced from the NASA Kepler mission (Borucki et al. 2011). Attempts to observationally characterize their atmospheres have already begun (Bean et al. 2010). In anticipation of this vanguard, several 3D circulation models of tidally locked, short-period terrestrial exoplanets have been published (Joshi et al. 1997; Joshi 2003; Merlis & Schneider 2010; Heng & Vogt 2010).

Intriguingly, the flows in most of these 3D models—both hot Jupiters and terrestrial planets—develop a fast eastward, or superrotating, jet stream at the equator, with westward flow typically occurring at deeper levels and/or higher latitudes. In hot-Jupiter models, the superrotating jet extends from the equator to latitudes of typically 20°–60° and is perhaps the dominant dynamical feature of the modeled flows. In some cases (depending on the strength of the imposed stellar heating and other factors), this jet causes an eastward displacement of the hottest regions from the substellar point by ∼10°–60° longitude. Showman & Guillot (2002) first predicted this feature and suggested that, if it existed on hot Jupiters, it would have important implications for infrared spectra and light curves. This prediction has been confirmed in Spitzer infrared observations of HD 189733b (Knutson et al. 2007, 2009), suggesting that this planet may indeed exhibit a superrotating jet.

However, despite the ubiquity of equatorial superrotation in 3D models of synchronously rotating exoplanets—and its relevance for observations—the mechanisms that produce this superrotation have yet to be identified. As demonstrated in a theorem due to Hide (1969), such superrotation cannot result from atmospheric circulations that are longitudinally symmetric or that conserve angular momentum per unit mass about the planetary rotation axis. The equatorial atmosphere is the region of the planet farthest from the planetary rotation axis, and a superrotating equatorial jet therefore corresponds to a local maximum in the angular momentum per unit mass about the planetary rotation axis. Thus, any angular-momentum-conserving circulation that moves air to the equatorial atmosphere from higher latitudes or deeper levels tends to produce westward equatorial flow. Equivalently, Coriolis forces always induce westward acceleration for air moving equatorward or upward, so an eastward equatorial jet cannot result from Coriolis forces acting on air that moves into the equatorial atmosphere from higher latitudes or deeper levels. In Earth's equatorial troposphere, for example, the flow is westward, which results from the tropospheric Hadley cell on Earth (a regime where a mean overturning circulation and its Coriolis accelerations plays the defining role; Held & Hou 1980). To maintain equatorial superrotation, a mechanism is needed that pumps angular momentum per unit mass from regions where it is small (outside the jet) to regions where it is large (within the jet)—a so-called up-gradient momentum transport. According to Hide's theorem, this transport can only be accomplished by waves or eddies.

Equatorial superrotation exists in several atmospheres of the solar system—the equatorial atmospheres of Venus, Titan, Jupiter, and Saturn all superrotate. Even localized layers within Earth's equatorial stratosphere exhibit superrotation, part of a phenomenon called the "Quasi-Biennial Oscillation" or QBO (Andrews et al. 1987). The mechanisms for driving the equatorial superrotation on these planets are diverse. Possible mechanisms include eddy transport associated with baroclinic instabilities, barotropic instabilities, turbulence, and the absorption/radiation of various types of atmospheric waves (e.g., Williams 2003a, 2003b; Lian & Showman 2008, 2010; Schneider & Liu 2009; Del Genio et al. 1993; Del Genio & Zhou 1996; Andrews et al. 1987; Mitchell & Vallis 2010).

Here, we demonstrate how the equatorial superrotation in 3D models of synchronously rotating, short-period exoplanets can result from the existence of standing, planetary-scale Rossby waves; such waves are naturally excited by the longitudinally dependent heating patterns—dayside heating and nightside cooling—that accompany the photospheric regions of short-period, synchronously rotating planets. Section 2 provides background. In Section 3, we present an analytic theory demonstrating the mechanism, and we systematically explore its behavior in idealized, nonlinear one-layer models. In Section 4, we extend the analysis to 3D circulation models. In Section 5, we summarize and discuss implications.

The ability of Rossby waves to accelerate jets can be schematically illustrated using the two-dimensional nondivergent barotropic vorticity equation, which is the simplest model for the global-scale flow of a planetary atmosphere (see discussion in Vallis 2006; Showman et al. 2010). The equation reads

$$\begin{equation} {d(\zeta + f)\over dt} = F, \end{equation} \tag{ 1 }$$

where ζ ≡ k · ∇ × v is the relative vorticity, v is the horizontal wind velocity, k is the local upward unit vector, f ≡ 2Ωsin ϕ is the Coriolis parameter, Ω is the planetary rotation rate (2π over the rotation period), ϕ is the latitude, and d/dt = ∂/∂t + v · ∇ is the material derivative (i.e., the derivative following the flow). The equation states that individual fluid parcels conserve the absolute vorticity, ζ + f, save for vorticity sources/sinks, which are represented by the term F. The equation can equivalently be written as

$$\begin{equation} {\partial \zeta \over \partial t} + {\bf v}\cdot \nabla \zeta + v\beta = F, \end{equation} \tag{ 2 }$$

where v is the meridional (northward) wind speed and β = df/dy is the gradient of the Coriolis parameter with northward distance y. Because the flow in this simple model is horizontally nondivergent, we can define a streamfunction ψ such that u = −∂ψ/∂y and v = ∂ψ/∂x, where x is the eastward coordinate and u is the zonal (eastward) wind speed. This allows the equation to be written as a function of one variable, ψ.

For purposes of illustration, consider the linearized version of Equation (2) with no sources and sinks. The solutions to this linearized, unforced equation are Rossby waves. For simplicity, we consider Cartesian geometry with constant β, representing a local region on the sphere. Decomposing variables into zonal means (denoted with overbars) and deviations therefrom (denoted with primes), and assuming that the mean flow is zero, leads to a solution given by $\psi ^{\prime }=\hat{\psi }\exp [i(kx+ ly)]$, where i is the imaginary number and k and l are the zonal and meridional wavenumbers. The dispersion relation is

$$\begin{equation} \omega = -{\beta k\over k^2 + l^2}. \end{equation} \tag{ 3 }$$

These waves propagate meridionally with a group velocity given by

$$\begin{equation} {\partial \omega \over \partial l}={2\beta kl\over (k^2+l^2)^2}. \end{equation} \tag{ 4 }$$

A simple argument, first clearly presented by Thompson (1971) and reviewed in Held (2000) and Vallis (2006), shows how these waves can produce an east–west acceleration of the zonal-mean flow. The latitudinal transport of eastward eddy momentum per unit mass is $\overline{u^{\prime }v^{\prime }}$, where u' and v' are the deviations of the zonal and meridional winds from their zonal means, respectively, and the overbar denotes a zonal average. Using the solutions for the wave-induced zonal and meridional wind, $u^{\prime }=-il\hat{\psi }\exp [i(kx+ly)]$ and $v^{\prime }=ik\hat{\psi }\exp [i(kx +ly)]$, yields a momentum flux

$$\begin{equation} \overline{u^{\prime }v^{\prime }}=-{1\over 2}\hat{\psi }^2 kl. \end{equation} \tag{ 5 }$$

Since the group velocity must point away from the region of wave generation (which we call the "wave source"), and since β is positive, we must have kl > 0 north of the source and kl < 0 south of the source. Thus, north of the source, $\overline{u^{\prime }v^{\prime }}$ is negative, implying southward transport of eastward momentum. But south of the source, $\overline{u^{\prime }v^{\prime }}$ is positive, implying northward transport of eastward momentum. Rossby waves therefore transport eastward momentum into the wave-source region. An eastward acceleration must therefore occur in the wave-source region and a westward acceleration must occur in the region of wave breaking or dissipation. This can lead to the formation of zonal (east–west) jet streams.⁵

Rossby waves correspond to latitudinal oscillations in surfaces of constant potential vorticity⁶; thus, any process that triggers such oscillations at large scales will tend to excite Rossby waves. In Earth's atmosphere, one of the predominant sources is baroclinic instability, which occurs in the midlatitude troposphere where latitudinal temperature gradients are large. Spatially varying tropospheric heating and cooling (e.g., due to land–sea contrasts) or flow over topography also perturb the potential vorticity contours and can therefore trigger Rossby waves. In the atmospheres of tidally locked, hot exoplanets, on the other hand, the day–night heating pattern constitutes the overriding dynamical forcing. For such planets, we expect this heating/cooling pattern to trigger Rossby waves at low latitudes (Figure 1).

Zoom In Zoom Out Reset image size

The above theory is for free waves. Consider now the extension to an atmosphere forced by vorticity sources/sinks and damped by frictional drag. The zonal-mean zonal momentum equation of the barotropic system reads

$$\begin{equation} {\partial \overline{u}\over \partial t}=-{\partial (\overline{u^{\prime }v^{\prime }})\over \partial y} - {\overline{u}\over \tau _{\rm drag}}, \end{equation} \tag{ 6 }$$

where overbars denote zonal means and primes denote deviations therefrom. The equation states that accelerations of the zonal-mean zonal flow result from convergences of the latitudinal eddy momentum flux and from drag, which we have parameterized as a term that relaxes the zonal-mean zonal wind toward zero over a drag time constant τ_drag. The relationship between the eddy acceleration in Equation (6) and the vorticity sources/sinks can be made in two steps. First, we note that the definition of vorticity implies that $\overline{v^{\prime }\zeta ^{\prime }} = - \partial (\overline{u^{\prime }v^{\prime }})/\partial y$. Second, we multiply the linearized version of Equation (2) by ζ' and zonally average. This leads to an equation for the budget of the so-called pseudomomentum (Vallis 2006, p. 493):

$$\begin{equation} {\partial {\cal A}\over \partial t} + \overline{v^{\prime }\zeta ^{\prime }} = {\overline{\zeta ^{\prime }F^{\prime }}\over 2(\beta -{\partial ^2 \overline{u}\over \partial y^2})}. \end{equation} \tag{ 7 }$$

For the two-dimensional nondivergent model, ${\cal A}=(\beta - \partial ^2 \overline{u}/\partial y^2)^{-1} \overline{\zeta ^{\prime 2}}/2$ is the pseudomomentum, which is a measure of wave activity. By combining Equations (6) and (7) and supposing that the wave amplitudes and zonal-mean zonal wind are statistically steady, i.e., $\partial {\cal A}/\partial t \approx 0$ and $\partial \overline{u}/\partial t \approx 0$, we obtain

$$\begin{equation} {\overline{u}\over \tau _{\rm drag}}= {\overline{\zeta ^{\prime }F^{\prime }}\over 2(\beta -{\partial ^2 \overline{u}\over \partial y^2})}. \end{equation} \tag{ 8 }$$

This equation relates the vorticity sources/sinks and drag to the zonal-mean zonal wind, $\overline{u}$. When eddy sources/sinks of relative vorticity on average exhibit the same sign as the vorticity itself (i.e., $\overline{\zeta ^{\prime }F^{\prime }}>0$), the eddy acceleration is eastward, and in steady state results in an eastward zonal-mean zonal wind. When sources/sink of relative vorticity tend to exhibit the opposite sign as the vorticity ($\overline{\zeta ^{\prime }F^{\prime }}<0$), the eddy acceleration is westward, and in steady state results in westward zonal-mean zonal wind.⁷ In analogy with the free solutions, this behavior is typically interpreted in terms of the generation, latitudinal propagation, and dissipation of Rossby waves.

This mechanism is thought to be responsible for the eddy-driven jet streams (and the associated eastward surface winds) in Earth's midlatitudes: baroclinic instability generates Rossby waves that radiate away from the midlatitudes, causing eastward eddy acceleration there and leading to eastward surface flow in steady state (Held 2000; Vallis 2006). At the equator, Earth's troposphere does not superrotate; nevertheless, idealized Earth general circulation models (GCMs) have shown that the presence of strong zonally varying heating and cooling in the tropics can cause equatorial superrotation to emerge (Suarez & Duffy 1992; Saravanan 1993; Kraucunas & Hartmann 2005; Norton 2006). In analogy with the theory described above, Held (1999) suggested heuristically that the superrotation in these models results from the generation and poleward propagation of Rossby waves by the tropical heating sources.

Still, application of this barotropic theory to tidally locked exoplanets is problematic. Most models of atmospheric circulation on synchronously rotating, zero-eccentricity hot Jupiters exhibit relatively steady circulation patterns whose velocity and temperature patterns are approximately symmetric about the equator (Showman & Guillot 2002; Cooper & Showman 2005, 2006; Showman et al. 2008, 2009; Dobbs-Dixon & Lin 2008; Rauscher & Menou 2010). In such models, the relative vorticity is approximately antisymmetric about, and zero at, the equator. Equation (8) predicts that $\overline{u}=0$ at the equator for this situation. Thus, while attractive, this theory fails to explain the equatorial superrotation in these hot-Jupiter models.

There are other difficulties. First, the theory presented here assumes the flow is barotropic—i.e., it can be described by Equation (1)—and thereby ignores the role of finite Rossby deformation radius in shaping the wave properties. Second, vertical motions normally accompany regions of heating and cooling, and these vertical motions lead to nonzero horizontal divergence. The low-latitude flow is thus inherently divergent in the presence of heating/cooling, in conflict with the stated assumptions. Moreover, the theory cannot account in any rigorous way for the generation of the Rossby waves by thermal forcing. The above equations include no thermodynamics and only the effect of such forcing in producing relative vorticity is represented. Away from the equator, one expects that heating produces rising motion and horizontal divergence aloft, which generates anticyclonic relative vorticity by the action of Coriolis forces. At the equator, however, this vorticity source is small, leaving unclear the applicability of this barotropic theory in producing eastward flow at the equator.⁸ When the geostrophic assumption is relaxed and a finite deformation radius is included, analytic solutions of freely propagating equatorial waves show that such waves tend to be trapped within an equatorial waveguide and cannot propagate away from the equator (Holton 2004; Andrews et al. 1987); such wave solutions involve no net meridional (north–south) momentum transports, thus raising the question of whether the above mechanism is viable at the equator in the presence of finite deformation radius. Finally, for hot Jupiters, the Rossby waves are expected to be global in scale, and it is not clear a priori whether there is room for them to propagate poleward from the equatorial regions.

We present a sequence of models in the following sections that overcome these obstacles and provide a theoretical foundation for understanding equatorial superrotation on tidally locked exoplanets.

Full GCM solutions, although useful, involve so many interacting processes that it is often difficult to cleanly identify specific dynamical mechanisms from such solutions (see, e.g., Showman et al. 2010). Simplified models therefore play an important role in atmospheric dynamics. Here, we present a highly idealized model intended to capture the mechanism in the simplest possible context.

We adopt a two-layer model, with constant densities in each layer, where the upper layer represents the meteorologically active atmosphere and the lower layer represents the quiescent deep atmosphere and interior. In the limit where the lower layer becomes infinitely deep and the lower-layer winds and pressure gradients remain steady in time (which requires the upper layer to be isostatically balanced), this two-layer system reduced to the shallow-water equations for the flow in the upper layer (e.g., Vallis 2006, p. 129–130):

$$\begin{equation} {d{\bf v}\over dt}+g\nabla h + f{\bf k}\times {\bf v} = {\bf R} - {{\bf v}\over {\tau _{\rm drag}}} \end{equation} \tag{ 9 }$$

$$\begin{equation} {\partial h\over \partial t} + \nabla \cdot ({\bf v}h) = {h_{\rm eq}(\lambda,\phi) - h\over \tau _{\rm rad}} \equiv Q, \end{equation} \tag{ 10 }$$

where v(λ, ϕ, t) is the horizontal velocity, h(λ, ϕ, t) is the upper layer thickness, t is the time, g is the (reduced) gravity,⁹ f = 2Ωsin ϕ is the Coriolis parameter, k is the upward unit vector, Ω is the planetary rotation rate, and ϕ is the latitude. Again, d/dt = ∂/∂t + v · ∇ is the material derivative. The boundary between the layers represents an atmospheric isentrope, across which mass is exchanged in the presence of heating or cooling. Heating and cooling are therefore represented in the shallow-water system using mass sources and sinks, represented here as a Newtonian relaxation of the height toward a specified radiative-equilibrium height h_eq—thick on the dayside and thin on the nightside—over a radiative timescale τ_rad. The momentum Equations (9) include drag with a timescale τ_drag, which could represent the potentially important effects of magnetohydrodynamic friction (Perna et al. 2010), vertical turbulent mixing (e.g., Li & Goodman 2010), or momentum transport by breaking gravity waves (Lindzen 1981; Watkins & Cho 2010).

The term R in Equation (9) represents the effect on the upper layer of momentum advection from the lower layer, and takes the same form as in Shell & Held (2004) and Showman & Polvani (2010):

$$\begin{equation} {\bf R}(\lambda,\phi,t)= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}-{Q{\bf v}\over h}, &Q>0;\\ 0, &Q<0, \end{array}\right. \end{equation} \tag{ 11 }$$

where λ and ϕ are longitude and latitude. Air moving out of the upper layer (Q < 0) does not locally affect the upper layer's specific angular momentum or wind speed, hence R = 0 for that case. But air transported into the upper layer carries lower-layer momentum with it and thus alters the local specific angular momentum and zonal wind in the upper layer. For the simplest case where the lower-layer winds are assumed to be zero, this process preserves the column-integrated vh of the upper layer, leading to the expression in Equation (11) for Q > 0. Importantly, the expression for R follows directly from the momentum budget and contains no free parameters.¹⁰

The relative roles of the two terms on the right-hand side of the momentum equation can be clarified by rewriting it as follows:

$$\begin{equation} {d{\bf v}\over dt}+g\nabla h + f{\bf k}\times {\bf v}\!=\! -{\bf v}\left[{1\over {\tau _{\rm drag}}} + {1\over \tau _{\rm rad}}\left({h_{\rm eq}-h\over h}\right) {\cal H}(h_{\rm eq}-h)\right], \end{equation} \tag{ 12 }$$

where ${\cal H}(h_{\rm eq}-h)$ is the Heaviside step function, defined as 1 when h_eq − h > 0 and 0 otherwise. Dynamically, it should be clear that both terms on the right-hand side play a role analogous to drag; one can define the entire quantity in square braces as one over an effective drag time constant. Still, the second term (R) is spatially heterogeneous and only exists in regions of heating, and we will show that its effect on the zonal-mean flow is qualitatively different than that of the first term (frictional drag). For a strongly irradiated hot Jupiter, we might expect (h_eq − h)/h ∼ 0.01–1, and if so then the first term would dominate if τ_drag ≪ τ_rad whereas the second term would dominate if τ_drag ≫ τ_rad.

We present linear, analytic solutions and fully nonlinear, numerically determined solutions of Equations (9)–(11) in the next two subsections.

3.1. Linear Solutions

To enable analytic solutions, we solve Equations (9)–(11) in Cartesian geometry assuming that the Coriolis parameter can be approximated as f = βy, where y is the northward distance from the equator and β (the gradient of Coriolis parameter with northward distance) is assumed constant. This approximation, called the "equatorial β-plane," is strictly valid only at low latitudes, but we will see in Section 3.2 that the qualitative features of these solutions are recovered by the full solutions in spherical geometry.

We now linearize Equations (9)–(11) about a state of rest. By definition, all the terms in the linearized equations have magnitudes that scale with the forcing amplitude. Note that the term R involves the product of the velocity and forcing amplitude, and therefore is quadratic in the forcing amplitude and does not appear in the linearized equations. (We will come back to it when evaluating the implications of the solutions for the zonal-mean zonal wind.) The linearized equations read

$$\begin{eqnarray} &&{\partial u\over \partial t} + g {\partial \eta \over \partial x} -\beta y v = -{u\over \tau _{\rm drag}} \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&{\partial v\over \partial t} + g{\partial \eta \over \partial y} + \beta y u = -{v \over \tau _{\rm drag}} \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&{\partial \eta \over \partial t} + H\left({\partial u\over \partial x}+{\partial v\over \partial y}\right) = S(x,y) - {\eta \over \tau _{\rm rad}}, \end{eqnarray} \tag{ 15 }$$

where x is the eastward distance and η is the deviation of the thickness from its constant reference value H, such that h = H + η. The quantity S ≡ (h_eq − H)/τ_rad is the forcing, which can also be expressed as S ≡ η_eq/τ_rad, where η_eq ≡ h_eq − H is the deviation of the radiative-equilibrium height from H.

The free solutions to these equations (i.e., when the right-hand sides are set to zero) are the well-known equatorially trapped wave modes, described for example in Holton (2004) and Andrews et al. (1987). Given the intense heating and cooling experienced by hot, tidally locked exoplanets, however, we seek solutions to the forced problem. Most 3D dynamical models of hot Jupiters exhibit relatively steady circulation patterns (Showman & Guillot 2002; Cooper & Showman 2005, 2006; Dobbs-Dixon & Lin 2008; Dobbs-Dixon et al. 2010; Showman et al. 2008, 2009; Rauscher & Menou 2010), and so we seek steady solutions in the presence of forcing and damping.

We nondimensionalize Equations (13)–(15) with a length scale $L=(\sqrt{gH}/\beta)^{1/2}$, a velocity scale $U=\sqrt{gH}$, and a timescale ${\cal T}=(\sqrt{gH}\beta)^{-1/2}$, which correspond respectively to the equatorial Rossby deformation radius, the gravity wave speed, and the time for a gravity wave to cross a deformation radius in the shallow-water system. The thickness is nondimensionalized with H, the drag and thermal time constants with ${\cal T}$, and the forcing with $H/{\cal T}$. This yields, for steady flows,

$$\begin{eqnarray} &&{\partial \eta \over \partial x} -yv = - {u\over \tau _{\rm drag}} \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&{\partial \eta \over \partial y} + yu = - {v\over \tau _{\rm drag}} \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &&\left({\partial u\over \partial x} + {\partial v\over \partial y}\right) = S(x,y) - {\eta \over \tau _{\rm rad}}, \end{eqnarray} \tag{ 18 }$$

where all quantities, including τ_rad and τ_drag, are now nondimensional.

In pioneering investigations, Matsuno (1966) and Gill (1980) obtained analytic solutions to Equations (16)–(18) for the special case where the drag and radiative time constants are equal and drag is neglected from the meridional momentum Equation (17); the latter assumption is called the "longwave approximation" because it is valid in the limit where zonal length scales greatly exceed meridional ones. On tidally locked exoplanets, however, the drag and radiative timescales can differ greatly, and the longwave approximation may not apply, because the flow exhibits comparable zonal and meridional scales. We therefore retain the full form of Equations (16)–(18).

Equations (16)–(18) can be combined to yield a single differential equation for the meridional velocity v (e.g., Wu et al. 2001):

$$\begin{eqnarray} &&\nonumber {1\over \tau _{\rm drag}} \left({\partial ^2v \over \partial y^2} + {\partial ^2 v\over \partial x^2}\right) + {\partial v\over \partial x} - {1\over \tau _{\rm rad}}\left(y^2 + {1\over \tau _{\rm drag}^2}\right)v\\ &&\quad=\left(-y{\partial S\over \partial x} +{1\over \tau _{\rm drag}}{\partial S\over \partial y}\right)\!{.} \end{eqnarray} \tag{ 19 }$$

If one seeks separable solutions, then, as described by Gill (1980) and Wu et al. (2001), the meridional structure of the solutions to this equation with finite τ_drag are the parabolic cylinder functions ψ_n(y), which are simply Gaussians times Hermite polynomials¹¹:

$$\begin{equation} \psi _n(y)=\exp \left(-{y^2\over 2 {\cal P}^2}\right) H_n\left({y\over {\cal P}}\right), \end{equation} \tag{ 20 }$$

where ${\cal P}\equiv (\tau _{\rm rad}/\tau _{\rm drag})^{1/4}$ is the fourth root of a Prandtl number.

Our goal is to specify the thermal forcing, S(x, y), and solve for the unknowns u, v, and η. In general, any desired pattern of thermal forcing can be represented as

$$\begin{equation} S(x,y) = \sum _{n=0}^{\infty } S_n(x) \psi _n(y). \end{equation} \tag{ 21 }$$

For tidally locked exoplanets, we expect this pattern to consist of a day–night variation in heating/cooling whose amplitude peaks at low latitudes and diminishes near the poles. We take the forcing to be symmetric about the equator (appropriate for a planet with zero obliquity) and, to keep the mathematics tractable, retain solely the term S₀, corresponding to the pattern of heating and cooling that is a Gaussian, centered about the equator, with a latitudinal half-width of the equatorial Rossby radius of deformation modified by frictional and radiative effects. While the full solution would require consideration of S_n for all n ⩾ 0, the first term, S₀, will be the dominant term for cases where the deformation radius is similar to a planetary radius, as is the case on typical hot Jupiters. Consideration of this term alone will therefore suffice to illustrate the qualitative features relevant for inducing an equatorially superrotating jet on tidally locked exoplanets.

Appendix B describes the solution method of Equations (16)–(18) and presents the solution for the specific case where the forcing consists solely of the S₀ term varying sinusoidally in longitude, i.e., $S(x,y)=\hat{S}_0 e^{ikx}\psi _0(y)$, where $\hat{S}_0$ is a constant. Figure 2 shows an example for parameter values typical of a hot Jupiter or hot super Earth (zonal wavelengths associated with the day–night heating contrast of a planetary circumference and radiative time constants of order 10⁵ s). For this example, the drag time constant is taken equal to the radiative time constant. Figure 2(a) shows the radiative-equilibrium height field and Figure 2(b) presents the steady-state height and velocity fields.

Zoom In Zoom Out Reset image size

The solutions exhibit several important features. Although the radiative-equilibrium height field is symmetric in longitude about the substellar point (Figure 2(a)), the actual height field deviates significantly from radiative equilibrium and exhibits considerable dynamical structure (Figure 2(b)). Two fundamental types of behavior are present. First, at mid-to-high latitudes (|y| ∼ 1–3 in the figure), the flow exhibits vortical behavior. The dayside contains an anticyclone in each hemisphere, manifesting as a pressure high (i.e., local maximum of the height) around which winds flow clockwise in the northern hemisphere and counterclockwise in the southern hemisphere; the nightside contains a cyclone in each hemisphere, manifesting as a pressure low around which winds flow counterclockwise in the northern hemisphere and clockwise in the southern hemisphere. Second, at low latitudes (|y| ≲ 1), the flows are nearly east–west; they diverge from a point east of the substellar longitude (marked with a cross in Figure 2(b)) and converge toward a point east of the antistellar longitude.

As discussed by Gill (1980), these features can be interpreted in terms of forced, damped, and steady equatorial wave modes. The mid-to-high latitude feature described above is dynamically analogous to that of an n = 1 equatorially trapped Rossby wave, which exhibits cyclones and anticyclones—alternating in longitude—that peak off the equator (see Matsuno 1966, Figure 4(c) for an example of the flow field in this mode). The low-latitude feature discussed above is dynamically analogous to a superposition of the n = 1 Rossby wave and the equatorial Kelvin wave, which is a fundamental equatorially trapped wave mode with strong zonal winds but very weak meridional winds and whose amplitude is symmetric about, and peaks at, the equator (see, e.g., Holton 2004 or Andrews et al. 1987). Both of these wave modes exhibit winds that are primarily east–west at the equator; in the example shown in Figure 2, the Kelvin component dominates over the n = 1 Rossby component at the equator.

Equations (16)–(18) indicate that this is a problem governed by two parameters: the radiative time constant and the drag time constant. We now examine how the behavior depends on their values. Figure 3 shows linear solutions, as presented in Appendix B, for dimensionless radiative time constants of 1, 10, and 100 (top, middle, and bottom rows, respectively) and drag time constants of 1, 10, and 100 (left, middle, and right columns, respectively). For parameters appropriate to hot Jupiters (rotation period of 3 Earth days and gH ≈ 4 × 10⁶ m² s⁻²), these dimensionless values correspond to dimensional time constants of ∼3 × 10⁴, 3 × 10⁵, and 3 × 10⁶ s, respectively (see Appendix A). When the radiative and drag time constants are short (upper left corner of Figure 3), the maximum and minimum thermal (height) perturbations lie on the equator and are close to the substellar and antistellar points; in this limit, the height field is close to radiative equilibrium (compare the top left of Figure 3 with Figure 2(a)),¹² and distinct Rossby-wave gyres do not appear. When the radiative and drag time constants have intermediate values (middle of Figure 3), cyclones and anticyclones become visible and—as in Figure 2(b)—exhibit height extrema that are phase shifted westward of the extrema in radiative equilibrium. Similarly, the height extrema along the equator become phase shifted eastward relative to radiative equilibrium; thickness variations along the equator become modest relative to those in midlatitudes. When the radiative and drag time constants are long (lower right corner of Figure 3), the height field becomes dominated by the off-equatorial anticyclones and cyclones, with minimal variation of height at the equator. In the limits τ_rad → ∞ and τ_drag → ∞, the solution becomes flat at the equator and is symmetric in longitude about the x = 0 axis (a point demonstrated explicitly in Appendix C); Figure 3 shows that this limit is almost reached even for τ_rad and τ_drag of 100.

Zoom In Zoom Out Reset image size

Much of the behavior in Figure 3 can be understood in terms of the zonal propagation of equatorially trapped Rossby and Kelvin modes. Kelvin waves exhibit eastward group propagation while long-wavelength, equatorially trapped Rossby waves exhibit westward group propagation. When τ_rad and τ_drag are very short (upper left corner of Figure 3), the damping is so strong that the waves are unable to propagate zonally. As a result, the height is close to the radiative equilibrium height field. When the two time constants have intermediate values, the propagation produces an eastward phase shift of the height field at the equator (the Kelvin component) and a westward phase shift of the height field in the off-equatorial cyclones and anticyclones (the Rossby component)—exactly as seen in Figure 2(b) and the middle of Figure 3. As the two time constants become very long, the westward phase offset of the off-equatorial cyclones and anticyclones achieves maximal values of 90°. At the equator, however, the height variations go to zero; this is explained by the fact that Coriolis forces are zero at the equator, so the linearized force balance is between pressure-gradient forces and drag. Weak drag requires weak pressure-gradient forces and hence a flat layer at the equator.

Now, the key point of our paper is that these linear solutions have major implications for the development of equatorial superrotation on tidally locked exoplanets. As can be seen in Figures 2(b) and 3, the wind vectors exhibit an overall tilt from northwest-to-southeast in the northern hemisphere and southwest-to-northeast in the southern hemisphere. This pattern, which resembles a chevron centered at the equator and pointing east, is particularly strong when the radiative and drag time constants are short, but occurs in all the cases shown. This structure implies that, on average, equatorward moving air has faster-than-average eastward wind speed while poleward moving air has slower-than-average eastward wind speed, so that $\overline{u^{\prime }v^{\prime }} < 0$ in the northern hemisphere and $\overline{u^{\prime }v^{\prime }} > 0$ in the southern hemisphere. As shown schematically in Figure 1, this is exactly the type of pattern that causes a flux of eastward eddy momentum to the equator and can induce equatorial superrotation. Since momentum is being removed from the midlatitudes, one would expect westward zonal-mean flow to develop there.

The physical mechanism responsible for producing these phase tilts is twofold. First, the differential wave propagation discussed above: this propagation causes an eastward phase shift of the height field in the Kelvin waves and a westward shift of the height field in the Rossby waves relative to the radiative-equilibrium height field. Because the Rossby wave lies on the poleward flanks of the Kelvin wave, the result is a chevron pattern where the height contours tilt northwest–southeast in the northern hemisphere and southwest–northeast in the southern hemisphere. To the extent that velocity vectors approximately parallel the geopotential contours (as they tend to do away from the equator when drag is weak or moderate), this will generate tilts in the velocities such that $\overline{u^{\prime }v^{\prime }}<0$ in the northern hemisphere and $\overline{u^{\prime }v^{\prime }}>0$ in the southern hemisphere.

The second mechanism for generating the velocity tilts needed for equatorial superrotation is simply the three-way force balance between Coriolis, drag, and pressure-gradient forces. Because drag acts opposite to the velocity, and Coriolis forces are perpendicular to the velocity, this three-way force balance requires the velocities to be rotated clockwise of −∇η in the northern hemisphere and counterclockwise of −∇η in the southern hemisphere. Given the expected day–night gradients in η, this balance implies that the velocities will tend to tilt northwest–southeast in the northern hemisphere and southwest–northeast in the northern hemisphere. We demonstrate this fact explicitly with an analytic solution in the limit of τ_rad → 0 in Appendix D; even when the height field is nearly in radiative equilibrium and hence exhibits no overall phase tilts, the velocities themselves develop tilts such that $\overline{u^{\prime }v^{\prime }}<0$ in the northern hemisphere and $\overline{u^{\prime }v^{\prime }}>0$ in the southern hemisphere (see Figure 15). The calculation in the limit τ_rad → 0 is particularly interesting because, in this limit, there is no zonal propagation of the Kelvin and Rossby waves: the radiative damping is infinitely strong and the zonal phase shift of the height field (relative to radiative equilibrium) is zero. This is the dominant mechanism for the velocity tilts in the top left panel of Figure 3.

To demonstrate explicitly how superrotation would emerge from these standing-wave patterns, we analyze the zonal accelerations associated with these linear solutions. Decomposing variables into their zonal means (denoted by overbars) and deviations therefrom (denoted with primes) and zonally averaging the zonal-momentum equation (Equation (9)) leads to (e.g., Thuburn & Lagneau 1999)

$$\begin{eqnarray} &&\nonumber {\partial \overline{u}\over \partial t}=\underbrace{\overline{v}^*\left[f - {\partial \overline{u}\over \partial y}\right]}_{{\rm I}} \nonumber \underbrace{-{1\over \overline{h}}{\partial \over \partial y} [\overline{(hv)^{\prime }u^{\prime }}]}_{{\rm II}}\\ &&+ \underbrace{\left[{1\over \overline{h}} \overline{u^{\prime }Q^{\prime }} + \overline{R_u}^*\right]}_{{\rm III}} \underbrace{-{\overline{u}^*\over \tau _{\rm drag}}}_{{\rm IV}} - {1\over \overline{h}}{\partial (\overline{h^{\prime }u^{\prime }})\over \partial t} \end{eqnarray} \tag{ 22 }$$

where a is the planetary radius and $\overline{A}^*\equiv \overline{hA}/ \overline{h}$ denotes the thickness-weighted zonal average of any quantity A. Equation (22) is the shallow-water version of the Transformed Eulerian Mean (TEM) momentum equation, analogous to that in the isentropic-coordinate form of the primitive equations (see Andrews et al. 1987, Section 3.9). On the right-hand side, terms I, II, and III represent accelerations due to (1) momentum advection by the mean-meridional circulation, (2) the convergence of the meridional flux of zonal eddy momentum, and (3) correlations between the regions of eddy zonal flow and eddy mass source (essentially vertical eddy-momentum transport). Within this term, the quantity R_u is the zonal component of R (equal to −Qu/h when Q > 0 and 0 when Q < 0). Term IV is the frictional drag. The final term represents the time rate of change of the eddy momentum. In the linear limit, all the terms on the right-hand side of Equation (22) have vanishingly small amplitude and, in this case, the solutions in Figures 2 and 3 represent true steady states. At any finite amplitude, however, terms I–IV are nonzero and would cause generation of a zonal-mean zonal flow.

Figure 2(c) depicts these terms for the example solution presented in Figure 2(b). As expected, horizontal convergence of eddy momentum, term II, causes a strong eastward acceleration at the equator and westward acceleration in the midlatitudes (black curve). On the other hand, the acceleration associated with vertical eddy-momentum transport, term III, is strong and westward at the equator (blue), implying downward transport of eddy momentum at the equator. The remaining terms—the mean-meridional circulation (term I, cyan) and mass-weighted friction (term IV, light green)—are small at the equator. The two eddy terms partially cancel at the equator, but the acceleration due to horizontal eddy momentum convergences exceeds that due to vertical eddy momentum convergences, leading to a net eastward acceleration at the equator and westward acceleration in midlatitudes (red curve).

Remarkably, despite the wide range of morphologies that occur when τ_rad and τ_drag are varied (Figure 3), all the solutions exhibit an equatorward flux of eddy momentum and a net eastward acceleration at the equator. This is shown in Figure 4, which presents the two eddy acceleration terms from Equation (22) for each of the cases shown in Figure 3. These solutions therefore suggest that superrotation at the equator and westward mean flow in the midlatitudes should occur at essentially any value of the control parameters.

Zoom In Zoom Out Reset image size

The patterns of spatial velocity and mass source/sink illuminate the physical origin of the westward equatorial acceleration caused by the vertical eddy exchange. The solutions show that the longitudes of zero zonal wind at the equator lie east of the mass-source extrema (Figure 2(b)), a feature also clearly visible in the steady, linear calculations of Matsuno (1966, Figure 9) and Gill (1980, Figure 1). Because of this shift, equatorial mass sources (sinks) occur predominantly in regions of westward (eastward) eddy zonal flow. On average, therefore, the mass sink regions transport air with eastward column-integrated eddy momentum out of the layer. The mass source regions transport air with no relative zonal momentum from the quiescent abyssal layer into the upper layer; this process conserves the local, column-integrated relative momentum vh of the upper layer. Thus, when zonally averaged, vertical exchange at the equator removes momentum from the layer, leading to $(\overline{u^{\prime }Q^{\prime }} + \overline{R_u}^*)<0$ and contributing a westward acceleration (blue curve in Figure 2(c)).

The above argument, however, does not determine which of the two eddy terms (II and III in Equation (22)) dominates. To determine which is larger—and hence whether the net equatorial eddy acceleration is eastward or westward—we write the zonally averaged zonal momentum equation in the form

$$\begin{equation} {\partial \overline{u}\over \partial t}=\overline{v^{\prime }\zeta ^{\prime }} + \overline{v}(f + \overline{\zeta }) - {\overline{u}\over \tau _{\rm drag}} + \overline{R_u}, \end{equation} \tag{ 23 }$$

where ζ is the relative vorticity. For the case where the forcing is symmetric about the equator, the solutions are symmetric about the equator in u and h but antisymmetric about the equator in v and ζ. As a result, the meridional velocity and relative vorticity are zero at the equator, so the terms $\overline{v}(f+\zeta)$ and $\overline{v^{\prime }\zeta ^{\prime }}$ vanish there. Therefore,

$$\begin{equation} {\partial \overline{u}\over \partial t}=-{\overline{u}\over \tau _{\rm drag}} +\overline{R_u}\qquad \qquad {\rm at}\;y=0. \end{equation} \tag{ 24 }$$

Essentially, at the equator, $\overline{R_u}$ is the mismatch between the accelerations caused by horizontal and vertical eddy-momentum fluxes. The analytic solutions, which assume $\overline{u}=0$, show that u is predominantly westward in regions where Q > 0, which therefore implies that $\overline{R_u}>0$. From Equation (24), the net eddy-induced acceleration is therefore eastward. This explains, in a general way, the sign of the net eddy accelerations at the equator in Figure 4. Of course, once a zonal-mean flow ($\overline{u}\ne 0$) develops, the magnitude of $\overline{R_u}$ changes and the friction term becomes important in Equation (24); eventually these terms balance and allow a steady state to be achieved. We discuss the possible steady states in light of this equation in Section 3.2.

We have so far emphasized the spatial patterns of the circulation, but it is also interesting to examine the magnitudes of the velocities predicted by our linear solutions. When the day–night difference in the radiative-equilibrium height is comparable to the mean value and the radiative time constant is a few days or less (as expected for the strongly forced conditions on hot Jupiters), the winds shown in Figures 2(b) and 3 reach nondimensional speeds of order unity. For a hot Jupiter, with typical g = 20 m s⁻² and H = 200 km, this corresponds to dimensional speeds of ∼2 km s⁻¹. To within a factor of a few, this is similar to the speeds obtained in fully nonlinear 3D atmospheric circulation models of hot Jupiters (Showman & Guillot 2002; Cooper & Showman 2005; Showman et al. 2008, 2009; Dobbs-Dixon & Lin 2008; Dobbs-Dixon et al. 2010; Menou & Rauscher 2009; Rauscher & Menou 2010; Thrastarson & Cho 2010). For a tidally locked, Earth-like planet in the habitable zone of an M-dwarf, with g = 10 m s⁻², H = 10 km, and an Earth-like radiative time constant of ∼10 days (corresponding to dimensionless time constants of 10–100), the solutions then yield nondimensional speeds of ∼0.02–0.1. This corresponds to dimensional speeds of up to a few tens of  m s⁻¹, similar to speeds obtained in models of tidally locked terrestrial planets (Joshi et al. 1997; Heng & Vogt 2010; Merlis & Schneider 2010).

3.2. Nonlinear Solutions

Next, we relax the small-amplitude and Cartesian constraints to demonstrate how nonlinearity and full spherical geometry affect the solutions, and we show how the wave-induced accelerations interact with the mean flow to generate an equilibrated state exhibiting equatorial superrotation. To do so, we solve the fully nonlinear forms of Equations (9)–(11) in global, spherical geometry, using a radiative-equilibrium thickness given by

$$\begin{equation} h_{\rm eq}=H + \Delta h_{\rm eq}\cos \lambda \cos \phi, \end{equation} \tag{ 25 }$$

where H is the mean thickness, Δh_eq is the day–night contrast in radiative-equilibrium thickness, and the substellar point is at longitude 0° and latitude 0°. The planet is assumed to be synchronously rotating, so that the pattern of h_eq(λ, ϕ) remains fixed in time. For concreteness, we adopt planetary parameters appropriate to a hot Jupiter, although we expect qualitatively similar solutions to apply to super Earths. For a typical gravity of 20 m s⁻² and scale height of 200 km appropriate to hot Jupiters, we might expect gH = 4 × 10⁶ m² s⁻², and we adopt this value for all our runs. (Note that g and H do not need to be specified independently.) We also take Ω = 3.2 × 10⁻⁵ s⁻¹ and a = 8.2 × 10⁷ m, corresponding to rotation period and planetary radius of 2.3 Earth days and 1.15 Jupiter radii, respectively, similar to the values for HD 189733b.

We reiterate that the equations represent a two-layer system with an active layer overlying a quiescent, infinitely deep lower layer. Because of coupling between the layers (specifically, mass exchange in the presence of heating/cooling), the solutions readily reach a steady state for any value of the drag time constant, including the limit where drag is excluded entirely in the upper layer (τ_drag → ∞). This in fact is a simple representation of the situation in many full 3D GCMs of solar system atmospheres, including Earth, which often have strong frictional drag near the surface, little-to-no friction in the upper layers, and yet easily reach a steady configuration throughout all the model layers. In our case, we find that, when drag is strong, the solutions reach steady states in runtimes ≲ 10τ_drag. In the case where drag is turned off, the time to reach steady state is determined by the magnitude of momentum and energy exchange between the layers (e.g., by the magnitude of the R term), and is generally ≲ 10τ_rad|H/Δh|, where |Δh|/H is a characteristic value of the fractional height variations in the active layer. All solutions shown here are equilibrated and steady.

We solve Equations (9)–(11) using the Spectral Transform Shallow Water Model (STSWM) of Hack & Jakob (1992). Rather than integrating the equations for u and v, the code solves the momentum equations in a vorticity-divergence form. The initial condition is a flat layer of geopotential gH at rest; the equations are integrated using a spectral truncation of T170, corresponding to a resolution of 07 in longitude and latitude (i.e., a global grid of 512 × 256 in longitude and latitude). A ∇⁶ hyperviscosity is applied to each of the dynamical variables to maintain numerical stability. The code adopts the leapfrog time-stepping scheme and applies an Asselin filter at each time step to suppress the computational mode. These methods are standard practice; for further details, the reader is referred to Hack & Jakob (1992).

To facilitate comparison with the analytic theory in Section 3.1, we first describe the solutions at very low amplitude where the behavior is linear. Figure 5 shows the geopotential (i.e., gh) for equal radiative and drag time constants of 0.1, 1, and 10 (Earth) days, respectively. Qualitatively, the numerical solutions in spherical geometry bear a striking resemblance to the analytic solutions on a β plane. At time constants of a fraction of day, the geopotential maxima occur on the equator, and for time constants of 0.1 days (panel (a)), the geopotential resembles the radiative-equilibrium solution, with wind flowing from the substellar point to the antistellar point. Longer time constants (1 day, panel (b)) allow zonal energy propagation of the Kelvin and Rossby waves, leading to an eastward phase shift of the geopotential at the equator and a westward phase shift at high latitudes (∼40°–90°). The result is contours of geopotential that develop northwest–southeast tilts in the northern hemisphere and southwest–northeast tilts in the southern hemisphere. When the time constants are long (10 days, panel (c)) off-equatorial cyclones and anticyclones dominate the geopotential, with only weak geopotential variations along the equator. These vortices are oval in shape, exhibiting no overall phase tilt, though the regions close to the equator do develop phase tilts (westward–poleward to eastward–equatorward). The momentum fluxes cause a prograde eddy acceleration (and superrotation) at the equator for all these cases. All of these features are also shared by the analytic solutions (Figure 3).

Zoom In Zoom Out Reset image size

We now explore how the solutions change when the radiative and frictional timescales are different. In the linear limit, the latitudinal width of the region exhibiting prograde phase tilts (i.e., northwest-to-southeast in the northern hemisphere and southwest-to-northeast in the southern hemisphere) contracts toward the equator when the drag time constant greatly exceeds the radiative time constant. This is illustrated in Figure 6, which shows the equilibrated (steady-state) solutions for τ_rad = 1 day and τ_drag/τ_rad = 1 (top), 10 (middle), and infinite (bottom). When the time constants are equal, the entire northern (southern) hemisphere exhibits northwest-to-southeast (southwest-to-northeast) phase tilts. When τ_drag/τ_rad = 10, these phase tilts are confined within ∼20° latitude of the equator, and for τ_drag → ∞, the width shrinks toward zero. This behavior is explained by the analytic theory in Section 3.1. As shown in Equation (20), the parabolic cylinder functions comprising the latitudinal structure exhibit a characteristic latitudinal width of L(τ_rad/τ_drag)^1/4, where $L=(\sqrt{gH}/\beta)^{1/2}$ is the equatorial Rossby deformation radius; these functions thus collapse toward the equator as τ_drag/τ_rad becomes infinite.¹³ Poleward of this region, the solutions exhibit phase tilts of the opposite direction (northeast-to-southwest in the northern hemisphere and southeast-to-northwest in the southern hemisphere). Appendix C gives the explanation for this reversal in phase tilts; the low-amplitude, full spherical numerical solutions at τ_drag/τ_rad ≫ 1 strongly resemble analytic solutions in the absence of drag, presented in Appendix C (compare Figures 6(c) and 14).

Zoom In Zoom Out Reset image size

Nonlinearity alters the solutions in several important ways, which we illustrate in Figure 7, showing a sequence of solutions for τ_rad = 0.1 days, τ_drag = 10 days and, from top to bottom, amplitudes Δh_eq/H of 0.01, 0.1, and 0.5, respectively. This choice of time constants is representative of the regime of strong radiative forcing and weak drag that may be appropriate to typical hot Jupiters. Increasing the forcing amplitude (i.e., increasing Δh_eq/H while holding τ_rad and τ_drag constant) of course leads to increased wind speeds and day–night geopotential variations; for the parameters in Figure 7, the zonal-mean zonal wind speed at the equator ranges from ∼10 m s⁻¹ at the lowest amplitude to almost 1000 m s⁻¹ for the highest amplitude shown. Moreover, beyond a critical value of Δh_eq/H (depending on the values of τ_rad and τ_drag), the solutions begin to deviate qualitatively from the linear solutions.

Zoom In Zoom Out Reset image size

First, nonlinearity allows greater geopotential variations to occur along the equator, such that at extreme forcing amplitude the geopotential extrema can in some cases occur along the equator when they otherwise would not. In the linear limit, the zonal force balance at the equator in the steady state is between the pressure-gradient force and drag (see Equation (13)); therefore, when drag is weak, the pressure-gradient force must likewise be small, implying that minimal variations of geopotential occur along the equator. This restriction does not apply at higher latitudes (where the Coriolis force can balance the pressure-gradient force), so for very weak drag the thickness extrema generally occur off the equator (as can be seen in the lower right portion of Figure 3; Figures 5(b) and (c); Figure 6; and Figure 7(a)). At large forcing amplitude, however, the momentum advection term v · ∇v and the R terms become important and can balance the pressure-gradient force, allowing significant zonal pressure gradients—and hence significant variations in thickness—to occur along the equator. For the parameters in Figure 7, the thickness variations peak at the equator when the forcing amplitude is sufficiently large (bottom panel).

Second, at high amplitude, the phase tilts of wind and geopotential tend to be from northwest-to-southeast (southwest-to-northeast) throughout much of the northern (southern) hemisphere—as in Figures 7(b) and (c)—even when the phase tilts are in the opposite direction at low amplitude (as in Figure 7(a)). This effect can be directly attributed to the term R in the momentum equations. As shown in Equation (12), R plays a role analogous to drag. When true drag is weak or absent, the effective drag time constant (one over the quantity in square brackets in Equation (12)) has a characteristic magnitude τ_radh/(h_eq − h). At large forcing amplitude, h/(h_eq − h) ∼ 1, and in that case the effective drag time constant is comparable to τ_rad. The linear solutions show that prograde phase tilts dominate over much of the globe when the radiative and drag time constants are comparable, but when the drag time constant greatly exceeds the radiative time constant, the phase tilts are in the opposite direction (see Figure 6). In Figure 7, τ_drag/τ_rad = 100, but the ratio of the effective drag time constant to the radiative time constant decreases from top to bottom and reaches ∼1 in the bottom panel, explaining the transition in the phase tilts from Figures 7(a) to (c). Through the momentum fluxes that accompany these phase tilts, the equatorial jet becomes broader and more dominant with increasing nonlinearity.

The momentum balance in the equatorial jet can achieve steady state in two ways. In steady state, Equation (24) becomes

$$\begin{equation} {\overline{u}\over \tau _{\rm drag}}=\overline{R_u} \qquad \qquad ({\rm at \;the\; equator}). \end{equation} \tag{ 26 }$$

As described previously, when the zonal-mean zonal winds are weak, the zonal wind u is predominantly westward in regions where Q > 0, so that $\overline{R_u}>0$ (Section 3.1 and Figure 4). This implies an eastward eddy acceleration of the zonal-mean zonal winds at the equator, which induces equatorial superrotation. In the first type of steady state, corresponding to a regime of strong friction (short drag time constant), this superrotation implies a strong westward acceleration due to friction ($-\overline{u}/\tau _{\rm drag}$). Steady state occurs when the zonal-mean equatorial jet becomes strong enough for the friction to balance the eastward eddy-induced acceleration at the equator. We call this the "high Prandtl number" regime. In the second type of steady state, which we call the "low Prandtl number" regime, the friction is sufficiently weak that the term $-\overline{u}/\tau _{\rm drag}$ is unimportant in the momentum balance. Because of the eastward eddy acceleration, the zonal-mean zonal winds can build to high speed. Once they do, they change the nature of $\overline{R_u}$—the larger $\overline{u}$ becomes, the smaller the extent to which u < 0 in the region Q > 0, as necessary for $\overline{R_u}>0$. Eventually, for sufficiently large $\overline{u}$, the quantity $\overline{R_u}$ goes to zero at the equator. The equatorial jet thus achieves steady state.

Figure 8 shows examples of each of these regimes illustrating how the momentum balance occurs. The left column presents an example with strong drag (τ_rad = τ_drag = 1 day) and the right column presents an example with weak drag (τ_rad = 0.1 days and τ_drag = ∞); both are equilibrated and steady. These are high-amplitude cases, so the thickness (top row, a and e) exhibits large fractional variations, and the phase tilts exhibit an overall trend of northwest-to-southeast (southwest-to-northeast) in the northern (southern) hemisphere, as explained in previous discussion (Figure 7). These tilts indicate transport of eddy momentum from midlatitudes to the equator; as a result, the zonal-mean zonal winds are eastward at the equator and westward in the midlatitudes (b and f). Interestingly, however, the relative strengths of the equatorial and midlatitude jets differ and reflect the range of possible variation. Panels (c) and (g) show the terms in the zonal-mean momentum equation, just the spherical equivalent of Equation (22):

$$\begin{eqnarray} &&\nonumber {\partial \overline{u}\over \partial t}=\underbrace{\overline{v}^*\left[f - {1\over a\cos \phi }{\partial (\overline{u}\cos \phi)\over \partial \phi }\right]}_{{\rm I}}\\ &&\nonumber \underbrace{-{1\over \overline{h}a\cos ^2\phi }{\partial \over \partial \phi } [\overline{(hv)^{\prime }u^{\prime }}\cos ^2\phi ]}_{{\rm II}}\\ &&+ \underbrace{\left[{1\over \overline{h}} \overline{u^{\prime }Q^{\prime }} + \overline{R_u}^*\right]}_{{\rm III}} \underbrace{-{\overline{u}^*\over \tau _{\rm drag}}}_{{\rm IV}} - {1\over \overline{h}}{\partial (\overline{h^{\prime }u^{\prime }})\over \partial t}. \end{eqnarray} \tag{ 27 }$$

As expected, horizontal convergence of eddy momentum, term II, causes a strong eastward acceleration at the equator and westward acceleration in midlatitudes (black curves). The vertical eddy-momentum transport, term III (dark blue curves), causes a westward acceleration at the equator that counteracts the eastward acceleration due to horizontal eddy-momentum convergence. In the case of strong drag (Figure 8(c)), the cancellation is imperfect, leading to a net eddy-induced acceleration that is eastward at the equator—as predicted by the linear, analytic theory in Section 3.1 (compare to Figure 2(c)). A superrotating equatorial jet therefore emerges and only reaches steady state when the jet becomes sufficiently strong that the zonal-mean drag on the jet, $-\overline{u}^*/ \tau _{\rm drag}$, balances the eastward acceleration at the equator (Figure 8(c)). On the other hand, when drag is absent, the superrotation induced by the eddy fluxes becomes quite strong (Figure 8(f)). This mean flow alters the eddy fluxes, causing them to self-adjust to an equilibrium where the accelerations at the equator due to horizontal and vertical momentum fluxes cancel, leading to no net eddy-induced acceleration at the equator in steady state (Figure 8(g)).

Zoom In Zoom Out Reset image size

We find that equatorial superrotation occurs at all forcing amplitudes, even arbitrarily small amplitudes where the solutions behave linearly. This is illustrated in Figure 9, which shows the equilibrated equatorial zonal-mean zonal wind versus forcing amplitude for solutions with a range of τ_rad and τ_drag combinations. We therefore conclude that the mechanism for generating equatorial superrotation described here has no inherent threshold. Nevertheless, other processes—not included in the shallow-water model—can in some cases overwhelm the desire of the day–night forcing to trigger superrotation, particularly when the day–night forcing is weak. This occurs for example in the cases examined by Suarez & Duffy (1992) and Saravanan (1993), where superrotation only developed for forcing amplitudes exceeding a threshold value. In their case, the tropical wave forcing only triggers superrotation when it attains sufficiently great amplitudes to overcome the westward torques provided by midlatitude eddies propagating into the tropics. These issues are discussed further in Section 5.

Zoom In Zoom Out Reset image size

In all the cases shown in Figure 9 where τ_drag is finite, the zonal-mean speed of equatorial superrotation scales with the square of the forcing amplitude when the forcing amplitude is sufficiently small. This in fact is the expected low-amplitude behavior in the high-Prandtl-number regime described above: at low amplitude, the solutions become linear, such that the velocities, height perturbations, and mass source/sink scale with the forcing amplitude. Because $\overline{R_u}$ scales as the product of the mass source/sink and the velocities at low amplitude, it is therefore quadratic in the forcing amplitude. In the frictional regime, Equation (26) implies that $\overline{u}$ at the equator is simply $\tau _{\rm rad}\overline{R_u}$, and therefore $\overline{u}$ itself is quadratic in the forcing amplitude. This behavior breaks down when the solutions become sufficiently high amplitude, as can be seen in Figure 9. The low-Prandtl-number regime is more complex and can lead to a variety of scaling behaviors depending on the parameters.

The flow in the shallow-water models differs from that in 3D models in one major respect. In many 3D models of hot Jupiters, eastward equatorial flow occurs not only in the zonal mean but at all longitudes, at least over some range of pressures. In contrast, although the shallow-water models described here all exhibit eastward zonal-mean flow at the equator, the zonal wind at the equator is always westward over some range of longitudes. This can be seen as follows: $\overline{R_u}$ is essentially the mismatch in equatorial zonal acceleration between horizontal and vertical eddy-momentum transport and in steady state, when $\overline{u}>0$, will be greater than or equal to zero. From the definition of R_u, this implies westward equatorial flow at some longitudes. This trait probably arises because the meteorologically active atmosphere has here been resolved with only one layer overlying a deep interior; in future work, it would be interesting to explore models that represent the flow with two or more layers overlying a quiescent interior to see whether they can develop equatorial flow that is eastward at all longitudes.

Here, we show how the basic mechanism for generating equatorial superrotation identified in Section 3 occurs also in 3D under realistic conditions. To do so, we analyze the 3D model of HD 189733b presented in Showman et al. (2009). Showman et al. (2009) coupled the dynamical core of the MITgcm (Adcroft et al. 2004), which solves the primitive equations of meteorology in global, spherical geometry, using pressure as a vertical coordinate, to the state-of-the-art, nongray radiative transfer scheme of Marley & McKay (1999), which solves the multi-stream radiative transfer equations using the correlated-k method to treat the wavelength dependence of the opacities. This coupled model, dubbed the Substellar and Planetary Atmospheric Circulation and Radiation (SPARC) model, is to date the only GCM to include realistic radiative transfer for hot Jupiters. Gaseous opacities were calculated assuming local chemical equilibrium for a specified atmospheric metallicity, assuming rainout of any condensates (i.e., ignoring cloud opacity). Showman et al. (2009) presented synchronously rotating models of HD 189733b with 1, 5, and 10 times solar metallicity and of HD 209458b with solar metallicity, along with several models with nonsynchronous rotation. Their HD 189733b models in particular compare favorably with a variety of observational constraints (Showman et al. 2009; Agol et al. 2010), and here we focus on their solar metallicity, synchronously rotating HD 189733b case. This model adopts planetary radius and gravity of 8.24 × 10⁷ m and 9.36 m s⁻². The rotation rate is 3.3 × 10⁻⁵ s, corresponding to a rotation period of 2.2 Earth days.

Figure 10 shows the velocity and temperature structure at the 30 mbar level during the spin-up phase of this model—after the forcing has had sufficient time to trigger a global wave response but before the equatorial jet has spun up to high speed. The velocity pattern in the 3D model (Figure 10) strongly resembles the standing Kelvin and Rossby-wave pattern described in Section 3. The flow clearly exhibits the east–west divergence along the equator, emanating from a point near the substellar longitude, identified in Section 3 as the standing Kelvin-wave response. The longitude of peak divergence (i.e., the longitude at the equator where the zonal velocity changes sign) lies east of the substellar longitude, as expected from the analytic theory and nonlinear shallow-water runs in Section 3. Moreover, the flow exhibits the broad gyres in each hemisphere, anticyclonic on the dayside and cyclonic on the nightside, identified in Section 3 as the standing Rossby-wave response. As predicted analytically, the velocities in these gyres exhibit a northwest-to-southeast (southwest-to-northeast) phase tilt in the northern (southern) hemisphere. These phase tilts imply that $\overline{u^{\prime }v^{\prime }}$ is negative in the northern hemisphere and positive in the southern hemisphere. Eddy momentum therefore fluxes from the midlatitudes to the equator, and it is this flux that produces the superrotating equatorial jet (see Figure 1). The overall qualitative resemblance to the analytic calculation in Figure 2 is striking.

Zoom In Zoom Out Reset image size

As in the shallow-water solutions, the 3D models exhibit a net downward eddy momentum flux at the equator throughout the upper atmosphere where the radiative heating/cooling is strong. This momentum flux results from the fact that, at the equator, (1) the Matsuno–Gill-type standing-wave patterns lead to net zonal eddy velocities that are predominantly westward on the dayside and eastward on the nightside (see Figure 10) and (2) net radiative heating occurs on much of the dayside, leading to net upward velocities, whereas net radiative cooling occurs on the nightside, leading to net downward velocities. Thus, at the equator, upward velocities tend to be correlated with westward eddy velocities and vice versa. This transports eastward momentum downward and causes a westward acceleration at the equator throughout the upper atmosphere, which counteracts the eastward equatorial acceleration caused by latitudinal eddy-momentum transport—just as predicted by the analytic and numerical shallow-water solutions in Section 3 (see Figures 2, 4, and 8).

To quantify the accelerations resulting from these momentum fluxes, we consider the Eulerian-mean zonal-momentum equation in pressure coordinates. By expanding the dynamical variables into zonal-mean and deviation (eddy) components, and zonally averaging the zonal-momentum equation, and adopting pressure as the vertical coordinate, we obtain¹⁴

$$\begin{eqnarray} &&\nonumber {\partial \overline{u}\over \partial t}=\overline{v}\left[f - {1\over a\cos \phi }{\partial (\overline{u}\cos \phi)\over \partial \phi }\right] - \overline{\omega }{\partial \overline{u}\over \partial p}\\ &&-{1\over a\cos ^2\phi }{\partial (\overline{u^{\prime }v^{\prime }}\cos ^2\phi)\over \partial \phi } -{\partial (\overline{u^{\prime }\omega ^{\prime }})\over \partial p} + \overline{X}. \end{eqnarray} \tag{ 28 }$$

On the right-hand side, the terms describe the meridional momentum advection by the zonal-mean circulation, vertical momentum advection by the zonal-mean circulation, the meridional eddy-momentum convergence, the vertical eddy-momentum convergence, and friction (represented generically by $\overline{X}$), respectively. At the equator, the Coriolis term is zero. Because of the approximate symmetry of the flow about the equator, $\overline{v}$ and the meridional gradient of $\overline{u}$ are small there, so the mean-meridional advection term is small at the equator. The mean vertical-advection term also tends to be weak for the flow considered here, and the net zonal acceleration at the equator is then determined primarily by a competition between the horizontal and vertical eddy-momentum convergence terms (the analogs of terms II and III in Equation (27) for the shallow-water system).

Figure 11 shows the time evolution of the zonal-mean zonal wind and the two eddy acceleration terms at the equator for the solar-metallicity model of HD 189733b from Showman et al. (2009). These are vertical averages through the top portion of the atmosphere where the radiative heating/cooling is strong. The zonal-mean zonal wind accelerates rapidly from the initial rest state and approaches an equilibrium within ∼100 days (top). As expected, the acceleration due to horizontal eddy transport is eastward, while that due to vertical eddy transport is westward (bottom). Moreover, as suggested by the linear and nonlinear shallow-water calculations, the magnitude of the horizontal momentum convergence exceeds that of the vertical momentum convergence during spin-up, so the net acceleration is eastward at early times. A superrotating equatorial jet therefore develops. As the jet speed builds, the two acceleration terms weaken significantly, and the ratio of their magnitudes approaches one. As a result, the net acceleration drops to zero, allowing the jet to equilibrate to a constant speed (top). This model is in the same regime as the shallow-water calculation presented in the right column of Figure 8.

Zoom In Zoom Out Reset image size

The weakening in time of the eddy accelerations seen in Figure 11 indicates that the mean flow, once it forms, exerts a back-reaction on the eddies that alters their structure. The nature of these changes is illustrated in Figure 12. The top panel shows the temperatures and winds at 30 mbar pressure after the flow at this altitude has become steady; the superrotating equatorial jet, eastward offset of the hottest region from the substellar point, and other features are evident as detailed in Showman et al. (2009). The bottom panel depicts the eddy temperature and eddy winds for the same pressure and time—that is, T' in color scale and (u',v') as arrows. Several features are similar to those in Figure 10: the eddy flow near the equator is approximately zonal and exhibits a Kelvin-wave-like character, with predominantly eastward flow at some longitudes and westward flow at others; the midlatitudes contain broad Rossby-wave gyres in each hemisphere, anticyclonic on the dayside and cyclonic on the nightside. Interestingly, however, the Kelvin-wave structure is shifted eastward, and the midlatitude velocity structure differs significantly, relative to that with weak mean flow (compare Figure 12(b) to Figures 2(b) and 10). From longitudes of about −75° to +40°, the midlatitude velocity structure induces equatorward momentum flux (i.e., u'v' negative in the northern hemisphere and positive in the southern hemisphere), but at longitudes ∼80°–150° the flux is reversed (i.e., u'v' positive in the northern hemisphere and negative in the southern hemisphere). Due to this cancellation, the magnitude of the zonally averaged flux $\overline{u^{\prime }v^{\prime }}$ is significantly weaker in the equilibrated state than during the spin-up phase, when the signs of the midlatitude u'v' add coherently at most longitudes (see Figure 10).

Zoom In Zoom Out Reset image size

The latitudinal pattern of zonal-mean eddy accelerations in the upper atmosphere of the 3D model, shown in Figure 13, exhibit a strong relationship to those from the shallow-water calculations. The comparison is most apt to shallow-water calculations with short radiative time constant (τ_rad ∼ 0.1–1 day), weak frictional drag (τ_drag → ∞), and large amplitude, as depicted for example in the right column of Figure 8. In the 3D model (Figure 13), the acceleration due to horizontal eddy-momentum convergence is eastward at the equator and westward in midlatitudes, in agreement with analytic theory and nonlinear shallow-water solutions (compare black curve in Figure 13 to Figures 4 and 8(g)). Near the poles, the situation is more complex. Poleward of ∼50° latitude, the acceleration in the upper atmosphere of the 3D run is eastward. In the shallow-water solutions, the pattern of eddy accelerations at high latitude depend on τ_rad, τ_drag, and the forcing amplitude (e.g., compare Figures 8(c) and (g) poleward of ∼60° latitude), but for the shallow-water cases most relevant to the 3D run shown here—such as the right column in Figure 8—the acceleration due to horizontal eddy convergence becomes eastward at high latitudes (Figure 8(g)), like that in the 3D run. In the 3D model, the acceleration due to vertical eddy-momentum convergence (blue curve in Figure 13) is westward at the equator and eastward in the midlatitudes, again like that arising in the shallow-water solutions, although a significant difference is that the midlatitude eastward acceleration is weak in the shallow-water runs but strong in the 3D run (relative to the magnitude of acceleration at the equator).

Zoom In Zoom Out Reset image size

Of course, the standing eddy patterns and resulting zonal-wind accelerations in 3D models depend on the strength of radiative heating/cooling and drag, just as they do in the shallow-water models. For example, in the shallow-water solutions, the Kelvin-wave structure and Rossby gyres are spatially distinct when the radiative and/or drag time constants are long and the forcing amplitude is small but not when the time constants are short or the forcing amplitude is large (e.g., compare the upper left versus the lower right of Figure 3, the top versus the bottom of Figure 5, and the top versus the bottom of Figure 7). The 3D models shown here lie at an intermediate position along this continuum, with Rossby and Kelvin-wave structures that are visibly distinct, analogous for example to the shallow-water case in Figure 2(b). 3D models with very strong heating rates, however, seem to exhibit eddy patterns lacking distinct Rossby wave gyres, more analogous to the top left case in Figure 3 and the top case in Figure 5. Examples of models in this regime include the topmost part of the atmosphere in the models of Cooper & Showman (2005, 2006), Koskinen et al. (2007), Rauscher & Menou (2010), and the HD 209458b model of Showman et al. (2009). In contrast, cases in the literature with more modest heating rates tend to exhibit distinct standing Rossby and Kelvin-wave structures; examples include Showman & Guillot (2002, Figure 5), Heng & Vogt (2010, Figures 1 and 12), the lower portion of some of the models of Koskinen et al. (2007, Figure 3(b)), and several of the runs in Thrastarson & Cho (2010), which exhibit a planetary-scale cyclone and anticyclone in each hemisphere.

Despite differences of detail, the overall broad similarities described here between the 3D and shallow-water models argue strongly that the mechanism for equatorial jet maintenance that we have identified occurs in both the shallow-water and 3D models.

The development of an eastward equatorial jet—that is, equatorial superrotation—is a common feature emerging from 3D models of synchronously rotating hot Jupiters and extrasolar terrestrial planets (Showman & Guillot 2002; Cooper & Showman 2005, 2006; Showman et al. 2008, 2009; Dobbs-Dixon & Lin 2008; Menou & Rauscher 2009; Rauscher & Menou 2010; Perna et al. 2010; Heng et al. 2010; Joshi et al. 1997; Merlis & Schneider 2010; Heng & Vogt 2010). Showman & Guillot (2002) first pointed out that, when the radiative and advective time constants are similar, this superrotation causes an eastward displacement of the hottest regions from the substellar point—a phenomenon discovered on HD 189733b five years later (Knutson et al. 2007, 2009). Despite its relevance, however, the dynamical mechanisms responsible for generating the equatorial superrotation on tidally locked exoplanets have not been previously identified.

Here, we have shown that the equatorial superrotating jet results from an interaction of the mean flow with standing, planetary-scale Rossby and Kelvin waves generated by the day–night thermal forcing. The strong longitudinal variations in radiative heating—namely intense dayside heating and nightside cooling—trigger the formation of standing, planetary-scale equatorial Rossby and Kelvin waves; this is essentially a linear response when wind speeds are modest, although nonlinearities affect the wave structure at high amplitude. The Kelvin waves straddle the equator while the Rossby waves lie on their poleward flanks. As a result of the differential zonal propagation—Kelvin waves propagating to the east and long-wavelength Rossby waves to the west—as well as the multi-way force balance between pressure-gradient, Coriolis, advective, and drag forces, the velocities develop tilts that resemble an eastward-pointing chevron centered at the equator. These velocity tilts pump eastward momentum from high latitudes to the equator, thereby inducing equatorial superrotation. In steady state, the zonal-mean equatorial jet speed near the photosphere is determined by a balance between this eastward, wave-induced acceleration and westward equatorial acceleration resulting from vertical eddy-momentum transport and/or drag. We demonstrated the mechanism in a hierarchy of dynamical models—including linear, analytic shallow-water models, fully nonlinear shallow-water models, and state-of-the-art 3D GCMs. For conditions relevant to hot Jupiters, such equatorial superrotation occurs over a wide range of radiative heating rates and drag time constants. The consistency of the picture emerging from this sequence of models with widely varying complexity is encouraging and suggests that the mechanism is robust.

The mechanism identified here has several implications.

• 1.  
  It implies that the equatorial jet results from a direct, essentially weakly nonlinear interaction between the thermally forced waves and the mean flow at the planetary scale. Eddy–eddy interactions, including the possibility of inverse or forward energy cascades or other turbulent interactions, may occur but are not essential to the basic mechanism. (This is analogous to the situation suggested by O'Gorman & Schneider (2007) for interaction of baroclinic midlatitude eddies with the mean flow on Earth.)
• 2.  
  The wave–mean-flow interaction produces an equatorial jet whose latitudinal width is comparable to that of the Rossby waves, namely the equatorial Rossby deformation radius modified by radiative and frictional effects. For conditions typical of synchronously rotating hot Jupiters, this length is comparable to a planetary radius, explaining the broad scale of the equatorial jet obtained in most hot Jupiter models.
• 3.  
  When acting in isolation, this mechanism for generating superrotation has no inherent forcing-amplitude threshold; it operates at any forcing amplitude, unlike the behavior reported in some earlier studies in the terrestrial context (Suarez & Duffy 1992; Saravanan 1993).
• 4.  
  For parameter combinations appropriate to hot exoplanets, our models yield flows whose hottest regions along the equator lie to the east of the substellar point. This property results from the eastward group propagation of Kelvin waves. The development of a strong mean flow (equatorial superrotation) further enhances the offset by its eastward advection of the temperature pattern. Together, these effects can explain the offsets observed on HD 189733b (Knutson et al. 2007, 2009).

Despite the ubiquity of superrotation in our models, the flow could shift regimes to one with westward zonal-mean equatorial flow if westward equatorial accelerations caused by other processes—not considered here—become important. For example, if baroclinic instabilities occur in midlatitudes, they could cause Rossby wave radiation at midlatitudes, potentially allowing absorption of Rossby waves near the equator. This would contribute a westward wave-induced acceleration near the equator. If this westward acceleration dominates over the eastward equatorial acceleration caused by the day–night heating contrast, then the net wave-induced acceleration at the equator could be westward and equatorial superrotation would not occur. This could occur on hot Jupiters if the planetary rotation rate is sufficiently fast and heating rates are sufficiently low. Note that baroclinic instabilities cannot occur in a one-layer model such as the shallow-water model, helping to explain why no such transitions to westward zonal-mean equatorial flow were seen in the shallow-water models presented here. In some cases, Hadley cells may also force the circulation into a regime of westward equatorial flow, particularly if the planetary obliquity is nonzero; this is the case in Earth's troposphere. On the other hand, if the planetary rotation is sufficiently slow, additional mechanisms for generating equatorial superrotation become possible, even when the forcing is axisymmetric (Del Genio & Zhou 1996; Mitchell & Vallis 2010). Exploring the combinations of orbital semimajor axes, stellar fluxes, planetary rotation rates, and atmospheric compositions under which such regime transitions can occur in 3D models is an important goal for future work, as such transitions could have major implications for visible and infrared light curves.

It is worth mentioning that the presence of a physical surface—and entropy gradients on that surface—promote the existence of baroclinic instabilities in midlatitudes (see, for example, Vallis 2006, chap. 6), so tidally locked terrestrial planets may be more prone than hot Jupiters to experience a regime transition to a flow containing midlatitude eastward jets in addition to, or instead of, a flow dominated by equatorial superrotation.

In the geophysical-fluid-dynamics (GFD) literature, the generation of midlatitude eastward zonal jets is often described theoretically using simple barotropic considerations analogous to those summarized in Section 2, involving the propagation of waves. While such barotropic considerations seem to work reasonably well for midlatitude jets, the barotropic framework fundamentally fails to explain the emergence of equatorial superrotation in 3D models of synchronously rotating exoplanets, where the flow is often steady and symmetric about the equator. The theory presented here overcomes this obstacle and represents an extension of the barotropic framework to a thermally forced, stratified, ageostrophic flow with finite Rossby deformation radius.

Specifically, the generation of eastward jets is often interpreted in terms of the meridional propagation of Rossby waves away from a source region and their dissipation at other latitudes (e.g., Held 2000; Vallis 2006). As described in Section 2, this interpretation seems to flow naturally from barotropic theory, in which free Rossby waves, even at the equator, can propagate not only in longitude but also in latitude. In contrast, although our work clearly shows how superrotation can emerge on tidally locked planets, meridional wave propagation plays no obvious role in our theory. Unlike barotropic Rossby waves, the baroclinic Rossby waves in our theory are equatorially trapped, confined to an equatorial waveguide: they can propagate in longitude and potentially height but not latitude.¹⁵ Moreover, under conditions appropriate to typical hot Jupiters, these waves exhibit meridional scales typically stretching from the equator to the pole. It is not at all clear that the paradigm of waves propagating from one latitude to another applies in this context. Rather, the velocity tilts that allow the meridional momentum fluxes to generate superrotation appear to be explained by the differential zonal—rather than meridional—propagation of equatorially trapped Kelvin and Rossby waves. Future work may further clarify the issue.

This research was supported by NASA Origins grant NNX08AF27G and PATM grant NNX10AB91G to A.P.S.

The nondimensional solutions to Equations (16)–(18) involve three dimensionless parameters—k, τ_rad, and τ_drag. Here we relate these dimensionless parameters to physical properties for exoplanets. For synchronously locked exoplanets, we expect the forcing to correspond to a zonal planetary wavenumber 1, implying a dimensional wavenumber of a⁻¹, where a is the planetary radius. Thus, our dimensionless wavenumber k has the value $a^{-1} (\sqrt{gH}/\beta)^{1/2}$, that is

$$\begin{equation} k = \left({\sqrt{gH}\over 2 \Omega a}\right)^{1/2}. \end{equation} \tag{ A1 }$$

If we assume for purposes of illustration that the mean layer thickness H is a scale height, then gH is just RT, where R is the specific gas constant and T is the mean temperature of the atmosphere. For a hydrogen atmosphere where $R=3700\rm \,J\rm \,kg^{-1}\rm \,K^{-1}$, this gives

$$\begin{equation} k = 0.75 \left({T\over 1000\rm \,K}\right)^{1/4} \left({P\over 3\,\rm days}\right)^{1/2}\left({R_{\rm J}\over a}\right)^{1/2}, \end{equation} \tag{ A2 }$$

where P is the rotational period and R_J is Jupiter's radius. For a carbon dioxide atmosphere where $R=189\rm \,J\rm \,kg^{-1}\rm \,K^{-1}$

$$\begin{equation} k=1.2\left({T\over 1000\rm \,K}\right)^{1/4} \left({P\over 3\,\rm days}\right)^{1/2}\left({R_{\oplus }\over a}\right)^{1/2}, \end{equation} \tag{ A3 }$$

where R_⊕ is Earth's radius. Thus, relevant values of k for a hot Jupiter and a hot Earth are similar, in the range of ∼0.5–2. Hotter atmospheres, longer rotation (=orbital) periods, and smaller planetary radii would promote larger values of k.

Consider now τ_rad and τ_drag. Their nondimensional values are their dimensional values times $(\sqrt{gH}\beta)^{1/2}$. Considering τ_nondim to be either dimensionless time constant (τ_rad or τ_drag) and τ_dim to be its dimensional counterpart, we have

$$\begin{equation} \tau _{\rm nondim} = \tau _{\rm dim} \left({2 \Omega \sqrt{gH}\over a} \right)^{1/2}. \end{equation} \tag{ A4 }$$

Again equating gH with RT, for a hot Jupiter with a hydrogen atmosphere we obtain

$$\begin{equation} \tau _{\rm nondim} = 3.6 \left({\tau _{\rm dim} \over 10^5\;{\rm s}}\right)\left({R_J\over a}\right)^{1/2} \left(T\over {1000{\rm \,K}}\right)^{1/4} \left({3\,{\rm days}\over P}\right)^{1/2}, \end{equation} \tag{ A5 }$$

while for a hot Earth with a CO₂ atmosphere

$$\begin{equation} \tau _{\rm nondim} = 5.7 \left({\tau _{\rm dim} \over 10^5\;{\rm s}}\right)\left({R_{\oplus }\over a}\right)^{1/2} \left(T\over {1000{\rm \,K}}\right)^{1/4} \left({3\,{\rm days}\over P}\right)^{1/2}. \end{equation} \tag{ A6 }$$

Here we present solutions to the nondimensional, linearized shallow-water Equations (16)–(18), subject to thermal forcing and drag. We follow the solution method outlined by Gill (1980) and Wu et al. (2001). For notational brevity, we define α ≡ τ⁻¹_drag and γ ≡ τ⁻¹_rad. Defining

$$\begin{eqnarray} &&q=\sqrt{\gamma }\eta + \sqrt{\alpha }u \end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray} &&r=\sqrt{\gamma }\eta - \sqrt{\alpha }u, \end{eqnarray} \tag{ B2 }$$

we convert the coupled equations for u, v, and η (Equations (16)–(18)) to equivalent equations for q, r, and v:

$$\begin{eqnarray} &&{\partial q\over \partial x} + \sqrt{\alpha }{\partial v\over \partial y} = \sqrt{\gamma }yv + \sqrt{\alpha } S(x,y) - \sqrt{\alpha \gamma } q \end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray} &&{\partial r\over \partial x} - \sqrt{\alpha }{\partial v\over \partial y}=\sqrt{\gamma }yv + \sqrt{\alpha \gamma } r - \sqrt{\alpha } S(x,y) \end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray} &&yq - yr + \sqrt{\alpha \over \gamma }{\partial q\over \partial y} + \sqrt{\alpha \over \gamma } {\partial r\over \partial y} = -2\alpha ^{3/2} v. \end{eqnarray} \tag{ B5 }$$

Gill (1980) and Wu et al. (2001) neglected the drag term in the meridional momentum equation, which is equivalent to dropping the term on the right-hand side of Equation (B5). However, we retain the full form of Equations (B3)–(B5). These equations are separable, and we adopt series solutions

$$\begin{equation} \lbrace q, r, v, S\rbrace = \sum _{n=0}^{\infty } \lbrace q_n(x), r_n(x), v_n(x), S_n(x)\rbrace \psi _n(y). \end{equation} \tag{ B6 }$$

Recursion relations for the parabolic cylinder functions are

$$\begin{eqnarray} &&{d\psi _n\over dy}={2n\psi _{n-1}\over {\cal P}} - {y \psi _n\over {\cal P}^2} \end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray} &&{d\psi _n\over dy}=-{\psi _{n+1}\over {\cal P}} + {y\psi _n\over {\cal P}^2}. \end{eqnarray} \tag{ B8 }$$

Inserting expressions (B6) into Equations (B3)–(B5), using the recursion relations (B7) and (B8), and invoking the orthogonality of the parabolic cylinder functions leads to the system

$$\begin{eqnarray} &&{dq_n\over dx} + \sqrt{\alpha \gamma }q_n - (\alpha \gamma)^{1/4} v_{n-1}=\sqrt{\alpha }S_n\quad n\ge 1 \end{eqnarray} \tag{ B9 }$$

$$\begin{eqnarray} &&{dr_n\over dx} - \sqrt{\alpha \gamma } r_n - 2(n+1)(\alpha \gamma)^{1/4} v_{n+1}= - \sqrt{\alpha }S_n\quad n\ge 0\qquad\; \end{eqnarray} \tag{ B10 }$$

$$\begin{eqnarray} &&2(n+1)q_{n+1} - r_{n-1}= -2\alpha ^{3/2} \left({\gamma \over \alpha }\right)^{1/4} v_n.\quad n\ge 1 \end{eqnarray} \tag{ B11 }$$

Equations (B9) and (B11) do not apply for n = 0, and Equations (B3)–(B5) instead yield for that case

$$\begin{eqnarray} &&{dq_0\over dx} + \sqrt{\alpha \gamma }q_0 = \sqrt{\alpha }S_0 \end{eqnarray} \tag{ B12 }$$

$$\begin{eqnarray} &&(\alpha /\gamma)^{1/4}q_1 = -\alpha ^{3/2} v_0. \end{eqnarray} \tag{ B13 }$$

Given a specified longitude and latitude dependence of the forcing (and hence given S_n(x) for all n ⩾ 0), our goal is to solve for q_n(x), r_n(x), and v_n(x). To determine q₀, we simply use Equation (B12). For n = 1, Equation (B9) relates q₁ and v₀ to S₁. To determine v₀, use Equation (B13). Inserting into Equation (B9), we obtain

$$\begin{equation} {dq_1\over dx} + \left(\sqrt{\alpha \gamma } + {1\over \alpha }\right)q_1 = \sqrt{\alpha }S_1. \end{equation} \tag{ B14 }$$

Obtaining an equation for q_n for n ⩾ 2 requires full use of Equations (B9)–(B11). First, rewrite Equation (B10) as an equation for dr_{n − 2}/dx in terms of r_{n − 2}, v_{n − 1}, and S_{n − 2}. Next, obtain an equation for r_{n − 2} from Equation (B11), and differentiate this expression to obtain an equation for dr_{n − 2}/dx. Inserting these two expressions into the equation derived from Equation (B10) yields

$$\begin{eqnarray} &&2n{dq_n\over dx} + 2\alpha ^{3/2}\left({\gamma \over \alpha }\right)^{1/4} {dv_{n-1}\over dx}\nonumber\\ &&- \sqrt{\alpha \gamma }\left[2nq_n + 2\alpha ^{3/2}\left({\gamma \over \alpha }\right)^{1/4}v_{n-1}\right]\nonumber\\ &&-2(n-1)(\alpha \gamma)^{1/4}v_{n-1}= -\sqrt{\alpha }S_{n-2}. \end{eqnarray} \tag{ B15 }$$

Equations (B9) and (B15) form two coupled differential equations for q_n and v_{n − 1} in terms of the known coefficients S_i. We solve for v_{n − 1} from Equation (B9) and insert this into Equation (B15) to obtain a single ordinary differential equation for q_n in terms of S_n, dS_n/dx, and S_{n − 2}:

$$\begin{eqnarray} &&\alpha {d^2q_n\over dx^2} + {dq_n\over dx} - [(2n-1)\alpha ^{1/2}\gamma ^{1/2} + \alpha ^2\gamma ]q_n\nonumber\\ &&= -[\alpha ^2\gamma ^{1/2} + (n-1)\alpha ^{1/2}]S_n +\alpha ^{3/2} {dS_n\over dx} - {\alpha ^{1/2}\over 2} S_{n-2}\nonumber\\ \end{eqnarray} \tag{ B16 }$$

for n ⩾ 2. With Equations (B12), (B14), and (B16), all possible q_n are determined from the specified forcing terms. To determine v_n from the q_n, we use Equation (B13) for n = 0, while for n ⩾ 1, we use Equation (B9):

$$\begin{equation} (\alpha \gamma)^{1/4}v_n = {dq_{n+1}\over dx} + \sqrt{\alpha \gamma }q_{n+1} - \sqrt{\alpha }S_{n+1}. \end{equation} \tag{ B17 }$$

To determine r_n from the q_n and v_n, we use Equation (B11):

$$\begin{equation} r_n = 2(n+2)q_{n+2} + 2\alpha ^{3/2}\left({\gamma \over \alpha }\right)^{1/4}v_{n+1} \end{equation} \tag{ B18 }$$

valid for all n ⩾ 0.

We now specify the forcing and solve for the response. For simplicity, consider a sinusoidal variation of all the variables in longitude:

$$\begin{equation} \lbrace q_n(x), r_n(x), v_n(x), S_n(x)\rbrace = \lbrace \hat{q}_n, \hat{r}_n, \hat{v}_n, \hat{S}_n\rbrace e^{ikx}, \end{equation} \tag{ B19 }$$

where $\hat{q}_n$, $\hat{r}_n$, $\hat{v}_n$, and $\hat{S}_n$ are complex constants and k is the dimensionless zonal wavenumber associated with the day–night heating contrast. We take the forcing to be symmetric about the equator (appropriate for a planet with zero obliquity) and, to keep the mathematics tractable, assume that the forcing is represented solely by the term S₀(x), corresponding to the pattern of heating and cooling that is a Gaussian, centered about the equator, with a half-width (in latitude) of the equatorial Rossby radius of deformation modified by frictional and radiative effects:

$$\begin{equation} S(x,y)=\hat{S}_0\psi _0(y)e^{ikx}. \end{equation} \tag{ B20 }$$

While the full solution would require consideration of S_n for all n ⩾ 0, the first term, S₀, will be the dominant term for cases where the deformation radius is similar to a planetary radius, as is the case on typical hot Jupiters. Consideration of this term alone will therefore suffice to illustrate the qualitative features relevant for pumping the equatorial jet on hot Jupiters.

With these assumptions, Equation (B12) implies

$$\begin{equation} \hat{q}_0 = {\sqrt{\alpha }(\sqrt{\alpha \gamma } - ik)\over \alpha \gamma + k^2}\hat{S}_0. \end{equation} \tag{ B21 }$$

Given that S₁ = 0, Equation (B14) implies that q₁ = 0. Similarly, Equation (B16) implies that

$$\begin{equation} \hat{q}_2 = {\alpha ^{3/2} k^2 + 3\alpha \gamma ^{1/2} + \alpha ^{5/2}\gamma + ik\alpha ^{1/2} \over 2[(\alpha k^2 + 3\alpha ^{1/2}\gamma ^{1/2} + \alpha ^2\gamma)^2 + k^2]} \hat{S}_0. \end{equation} \tag{ B22 }$$

All $\hat{q}_n$, with n ⩾ 3, equal zero. From Equation (B13), v₀ = 0, whereas

$$\begin{equation} \hat{v}_1 = {(ik + \sqrt{\alpha \gamma })\over (\alpha \gamma)^{1/4}}\hat{q}_2 \end{equation} \tag{ B23 }$$

which can be expressed in terms of the real and imaginary components of $\hat{q}_2=\hat{q}_{2_{\rm real}}+ i\hat{q}_{2_{\rm imag}}$, as

$$\begin{eqnarray} \hat{v}_1&=&\left[(\alpha \gamma)^{1/4}\hat{q}_{2_{\rm real}}- {k\over (\alpha \gamma)^{1/4}}\hat{q}_{2_{\rm imag}}\right]\nonumber\\ && + i\left[(\alpha \gamma)^{1/4}\hat{q}_{2_{\rm imag}}+ {k\over (\alpha \gamma)^{1/4}}\hat{q}_{2_{\rm real}}\right]. \end{eqnarray} \tag{ B24 }$$

All v_n, with n ⩾ 2, equal zero. From Equation (B18),

$$\begin{equation} \hat{r}_0 = (4 + 2\alpha ^{3/2}\sqrt{\gamma } + 2\alpha i k)\hat{q}_2 \end{equation} \tag{ B25 }$$

or equivalently

$$\begin{eqnarray} \hat{r}_0&=&\big[\big(4 + 2\alpha ^{3/2}\sqrt{\gamma }\big)\hat{q}_{2_{\rm real}}- 2\alpha k\hat{q}_{2_{\rm imag}}\big]\nonumber\\ &&+ i\big[\big(4 + 2\alpha ^{3/2} \sqrt{\gamma }\big)\hat{q}_{2_{\rm imag}}+ 2\alpha k\hat{q}_{2_{\rm real}}\big]. \end{eqnarray} \tag{ B26 }$$

All r_n, with n ⩾ 1 equal zero. This completes the solutions for q, r, and v. The solutions for u and η are then determined from (see Equations (B1) and (B2))

$$\begin{eqnarray} &&\eta ={q+r\over 2\sqrt{\gamma }} \end{eqnarray} \tag{ B27 }$$

$$\begin{eqnarray} &&u={q-r\over 2\sqrt{\alpha }}. \end{eqnarray} \tag{ B28 }$$

It is interesting to consider limits of these solutions as γ and α become infinitely large (implying that the radiative or drag time constants go to zero). The solution presented above implies that, in the limit α → ∞ at constant γ, or in the limit γ → ∞ at constant α, the quantities $\hat{q}_0 \rightarrow \hat{S}_o/\sqrt{\gamma }$, $\hat{q}_2\rightarrow \hat{S}_0/(2\alpha ^{3/2}\gamma)$, and $\hat{r}_0\rightarrow \hat{S}_0/\sqrt{\gamma }$. Noting that $\eta =(2\sqrt{\gamma })^{-1} [\hat{q}_0\psi _0(y) + \hat{q}_2\psi _2(y) + \hat{r}_0\psi _0(y)]e^{ikx}$, these limits imply that $\eta \rightarrow \hat{S}_0 \gamma ^{-1} e^{ikx}$, which, by definition, is simply the radiative-equilibrium height field. Thus, when either time constant goes to zero, the height field becomes the radiative-equilibrium height field—even if the other time constant is nonzero.

The steady, linearized, nondimensional shallow-water equations on an equatorial β plane (Equations (16)–(18)) have a particularly simple analytic solution in the case when frictional drag is set to zero. In this case, the nondimensional equation governing the meridional velocity (Equation (19)) becomes

$$\begin{equation} {y^2 v\over \tau _{\rm rad}} - {\partial v\over \partial x}= y{\partial S\over \partial x}. \end{equation} \tag{ C1 }$$

This equation differs vastly from Equation (19) because it no longer contains any y derivatives. As before, we seek separable solutions that are sinusoids in x. We specify the forcing amplitude $S(x,y)=\tilde{S}(y)e^{ikx}$ (where $\tilde{S}(y)$ is a specified function that describes the latitude dependence of the radiative heating/cooling), and search for solutions $u=\tilde{u}(y)e^{ikx}$, $v=\tilde{v}(y) e^{ikx}$, and $\eta =\tilde{\eta }(y)e^{ikx}$, $\tilde{u}(y)$, $\tilde{v}(y)$, and $\tilde{\eta }(y)$ are complex functions of latitude that we seek to determine. Inserting these expressions into Equation (C1), we obtain

$$\begin{equation} \left({y\over \tau _{\rm rad}}- {ik\over y}\right)\tilde{v} = ik\tilde{S} \end{equation} \tag{ C2 }$$

from which we have the solution

$$\begin{equation} \tilde{v} = {\tilde{S}(y)} k \tau _{\rm rad}{i y^3 - k\tau _{\rm rad}y \over y^4 + k^2\tau _{\rm rad}^2}. \end{equation} \tag{ C3 }$$

The height field is determined by the fact that ∂η/∂x = yv (Equation (16)), which yields

$$\begin{equation} \tilde{\eta }= \tilde{S}(y)\tau _{\rm rad}{y^4 + iy^2 k \tau _{\rm rad}\over y^4 + k^2\tau _{\rm rad}^2}. \end{equation} \tag{ C4 }$$

We can then determine u using yu = −∂η/∂y (Equation (17)), which yields

$$\begin{eqnarray} \tilde{u}&=&\tilde{S}(y)\tau _{\rm rad}\left[4{y^6+ iy^4k\tau _{\rm rad}\over (y^4 + k^2\tau _{\rm rad})^2} - {4y^2 + 2ik\tau _{\rm rad}\over y^4 + k^2\tau _{\rm rad}^2}\right]\nonumber\\ &&- {\partial \tilde{S}(y)\over \partial y} \tau _{\rm rad}\left[{y^3 + iyk\tau _{\rm rad}\over y^4 + k^2\tau _{\rm rad}^2} \right]. \end{eqnarray} \tag{ C5 }$$

Note that $\tilde{S}(y)$ can be any function; there is no need to decompose the solution into a summation over normal modes (e.g., parabolic cylinder functions), as in the case of finite τ_drag.

The solutions take on a particularly simple form in the limit τ_rad → ∞:

$$\begin{eqnarray} &&\tilde{v}=-\tilde{S}(y) \end{eqnarray} \tag{ C6 }$$

$$\begin{eqnarray} &&\tilde{\eta }={i\tilde{S}(y) y^2\over k} \end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray} &&\tilde{u} =4\tilde{S} i\left[{y^4\over k^3} - {2\over k}\right] - {\partial \tilde{S}\over \partial y}{iy\over k}. \end{eqnarray} \tag{ C8 }$$

It can be seen that, in this limit, the amplitudes of v are in phase with S. The peaks in η and u are shifted by 90° in longitude relative to S(x,y) and are zero at x = 0. The solution in this case is mirror symmetric about the y-axis, unlike the case with finite time constants.

Figure 14 illustrates the solutions for the case $\tilde{S}(y)=\hat{S}e^{-l^2y^2}$ with l = 0.6 and $\hat{S}$ being a real constant. As in Figures 2 and 3, the solutions exhibit midlatitude cyclones and anticyclones, with equatorial flow that is zonally aligned and diverges from a longitude near the substellar point. Because drag is zero, the layer is flat at the equator (see Equation (C4), which shows that η = 0 for y = 0). Interestingly, when τ_rad is finite, the westward phase shift of the Rossby waves is large at low latitudes and approaches zero at high latitudes (Figure 14(b)). This leads to phase tilts that are southwest–northeast in the northern hemisphere and northwest–southeast in the southern hemisphere—opposite to the cases shown in Figures 2 and 3. In the limit τ_rad → ∞, these phase tilts disappear and the solution strongly resembles that shown in the bottom right corner of Figure 3.

Zoom In Zoom Out Reset image size

Analysis of the continuity equation explains these behaviors. In addition to the mass source/sink caused by the forcing S, two processes affect the local layer thickness—radiative damping and horizontal convergence/divergence. Consider how their relative strengths depend on latitude. In the absence of drag, the linearized force balance is geostrophic, i.e., yv = ∂η/∂x and yu = −∂η/∂y (see Equations (16) and (17)). This means that, for a given velocity amplitude, the thickness gradients—and hence lateral thickness variations themselves—will be larger farther from the equator. As a result, for a given velocity amplitude, the radiative damping (which is −η/τ_rad) is stronger farther from the equator. On the other hand, in geostrophic balance, the wind flows parallel to contours of constant height, which severely limits the horizontal convergence/divergence; convergence can only come about due to variations of Coriolis parameter f with latitude, which are stronger near the equator. In geostrophic balance, the horizontal divergence is −βv/f, which is just −v/y for the equatorial beta plane considered here. Thus, for a given velocity amplitude, the amplitude of horizontal convergence is large near the equator but small at high latitudes. Given these latitude dependences, we thus expect that the thickness changes caused by the forcing (S) will predominantly be balanced by radiative relaxation at high latitude but horizontal convergence/divergence at low latitude.

Zoom In Zoom Out Reset image size

Equations (C1) and (C2), which represent the zonal derivative of the continuity equation, describe exactly this balance. The three terms correspond to local changes in layer thickness due to the forcing S (right side), mass source/sinks due to radiative damping (first term on left side), and changes in the layer thickness due to horizontal convergence or divergence of the fluid flow (second term on left side). For a given amplitude of v, the radiative damping increases with latitude (scaling with y), while the effect of horizontal convergence decreases with latitude (scaling with y⁻¹). The transition occurs at $y\sim \sqrt{k\tau _{\rm rad}}$ (or $y\sim \sqrt{gH k \tau _{\rm rad}/\beta }$ using dimensional quantities): radiative relaxation balances the forcing poleward of this latitude, whereas thickness changes due to horizontal convergence balance the forcing equatorward of this latitude.

When radiative relaxation balances the forcing ($y\gtrsim \sqrt{k\tau _{\rm rad}}$), the height field is in phase (in longitude) with the radiative equilibrium height field. When convergence/divergence balances the forcing ($y\lesssim \sqrt{k\tau _{\rm rad}}$), the convergence ∂u/∂x + ∂v/∂y, and therefore v itself, are phase with the radiative-equilibrium height field. Since ∂η/∂x = yv, this implies that the height field is phase shifted westward by 90° relative to the radiative equilibrium height field. In the transition zone ($y\sim \sqrt{k\tau _{\rm rad}}$), these arguments imply that the phase lines of thickness—and therefore the velocities themselves—tilt from southwest-to-northeast in the northern hemisphere and northwest-to-southeast in the southern hemisphere. This explains the phase tilts seen in Figure 14(b). In the limit τ_rad → ∞, the whole domain lies in the region where convergence balances forcing, explaining the lack of phase tilts of $\tilde{\eta }$ in Figure 14(c).

There also exist simple analytic solutions to the dimensionless, linearized shallow-water Equations (16)–(18), for general τ_drag, in the limit τ_rad → 0. This limit is particularly relevant for the hottest of tidally locked exoplanets, which, due to high temperature and/or significant visible-wavelength opacity, are expected to have short radiative time constants. Moreover, as a simplification of the full system it provides insights into the dynamical mechanisms.

In Appendix B, we showed that, in the limit τ_rad → 0, the height field converges toward the radiative-equilibrium height field. Thus, in the momentum equations, we can replace the height with radiative-equilibrium height. The continuity equation involves the difference between η and η_eq, however, and so we must retain η in that equation. This leads to the system:

$$\begin{eqnarray} &&{\partial \eta _{\rm eq}\over \partial x}-yv=-{u\over \tau _{\rm drag}} \end{eqnarray} \tag{ D1 }$$

$$\begin{eqnarray} &&{\partial \eta _{\rm eq}\over \partial y}+yu=-{v\over \tau _{\rm drag}} \end{eqnarray} \tag{ D2 }$$

$$\begin{eqnarray} &&\left({\partial u\over \partial x}+ {\partial v\over \partial y} \right)= {\eta _{\rm eq}-\eta \over \tau _{\rm rad}}. \end{eqnarray} \tag{ D3 }$$

Equations (D1) and (D2) constitute an algebraic system for u and v that can directly be solved to yield

$$\begin{eqnarray} &&u=-{{\partial \eta _{\rm eq}\over \partial x}+ y\tau _{\rm drag}{\partial \eta _{\rm eq}\over \partial y}\over {1\over \tau _{\rm drag}} + y^2\tau _{\rm drag}} \end{eqnarray} \tag{ D4 }$$

$$\begin{eqnarray} &&v={-{\partial \eta _{\rm eq}\over \partial y}+ y\tau _{\rm drag}{\partial \eta _{\rm eq}\over \partial x}\over {1\over \tau _{\rm drag}} + y^2\tau _{\rm drag}}. \end{eqnarray} \tag{ D5 }$$

Given u and v, Equation (D3) can then be solved to determine the (tiny) mismatch between η and η_eq.

Consider the limits of this solution for extreme values of the drag time constant. As τ_drag → 0, Equations (D1) and (D2) imply that u = −τ⁻¹_drag∂η_eq/∂x and v = −τ⁻¹_drag∂η_eq/∂y. In this strong-drag limit, the winds simply flow down the pressure gradient. On the other hand, in the limit of large τ_drag, then away from the equator, we obtain u = −y⁻¹∂η_eq/∂y and v = y⁻¹∂η_eq∂x—implying geostrophic balance. Note that the approximate solutions (D4) and (D5) are singular at the equator in the limit τ_drag → ∞; in this limit, the divergence becomes infinite at the equator, and the assumption that η = η_eq used to derive Equations (D4) and (D5) no longer holds.

Figure 15 displays this solution for a drag-dominated case (τ_drag = 0.01, left column) and case where drag is comparable to Coriolis forces over much of the domain (τ_drag = 1, right column). The top row shows the winds with the assumed height field, and the bottom row shows the winds and the horizontal divergence. As expected, when drag is strong (left column), the air flows directly from the substellar point toward the antistellar point and is almost parallel to the pressure gradients. When drag is less dominant, however (right column), the solution forms a Matsuno–Gill-type wind pattern which exhibits velocities that tilt northwest–southeast in the northern hemisphere and southwest–northeast in the southern hemisphere. As discussed in Section 3, this pattern of velocity tilts would generate equatorial superrotation.

The mechanism for generating these velocity tilts is, simply, the three-way force balance between pressure-gradient, Coriolis, and drag forces. Because drag points in the opposite direction of the velocities, and Coriolis forces are perpendicular to the velocities, this three-way balance must give a velocity direction that is rotated clockwise of −∇η in the northern hemisphere and counterclockwise of −∇η in the southern hemisphere. A visual inspection of Figure 15 makes clear that, given the expected pattern of ∇η_eq, this rotation forces the flow pattern to adopt velocity tilts that are northwest–southeast (southwest–northeast) in the northern (southern) hemisphere.