Wave-mean Flow Interactions in the Atmospheric Circulation of Tidally Locked Planets

The dispersion relation for the free modes is:

$$\begin{eqnarray}&&{\omega }^{2}-{k}_{x}^{2}-\displaystyle \frac{{k}_{x}}{\omega }=2n+1\end{eqnarray} \tag{ 3 }$$

and the free modes ξ are:

$$\begin{eqnarray}&&\left(\begin{array}{l}v\\ u\\ h\end{array}\right)=\left(\begin{array}{c}i(-{\omega }_{{nl}}^{2}-{k}_{x}^{2}){\psi }_{n}\\ \displaystyle \frac{1}{2}(-{\omega }_{{nl}}-{k}_{x}){\psi }_{n+1}+n(-{\omega }_{{nl}}+{k}_{x}){\psi }_{n-1}\\ \displaystyle \frac{1}{2}(-{\omega }_{{nl}}-{k}_{x}){\psi }_{n+1}-n(-{\omega }_{{nl}}+{k}_{x}){\psi }_{n-1}\end{array}\right)\end{eqnarray} \tag{ 4 }$$

where ${\psi }_{n}={e}^{-\tfrac{1}{2}{y}^{2}}{H}_{n}(y)$ (the parabolic cylinder functions, see Appendix A.3). In these solutions, l marks the three roots of each mode denoted by n (Matsuno 1966).

Matsuno (1966) derives the stationary response to forcing Q with a damping rate α. We introduce forcing and damping terms and then impose a sinusoidal dependence on x as before:

$$\begin{eqnarray}&&\begin{array}{l}{\alpha }_{\mathrm{dyn}}u-{yv}+{{ik}}_{x}h={F}_{x}\\ {\alpha }_{\mathrm{dyn}}v+{yu}+\displaystyle \frac{\partial h}{\partial y}={F}_{y}\\ {\alpha }_{\mathrm{rad}}h+{{ik}}_{x}u+\displaystyle \frac{\partial v}{\partial y}=Q\end{array}\end{eqnarray} \tag{ 5 }$$

for constant non-dimensional forcing $\sigma =({F}_{x},{F}_{y},Q)$, dynamical damping α_dyn, and radiative damping α_rad. Matsuno (1966) sets α_dyn = α_rad = α and imposes a Gaussian forcing on the height h to give an exact solution:

$$\begin{eqnarray}&&\sigma =\left(\begin{array}{c}0\\ 0\\ {Q}_{0}{e}^{-\tfrac{1}{2}{y}^{2}}\end{array}\right).\end{eqnarray} \tag{ 6 }$$

The boundary conditions are:

$$\begin{eqnarray}&&u,v,h\to 0\quad \mathrm{for}\quad y\to \pm \infty .\end{eqnarray} \tag{ 7 }$$

Matsuno (1966) shows that a forced solution χ can be expressed as a sum of the free solutions ξ:

$$\begin{eqnarray}&&\chi ={\rm{\Sigma }}{a}_{m}{\xi }_{m}\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{a}_{m} & = & \displaystyle \frac{1}{\alpha -i{\omega }_{m}}{b}_{m}\\ {b}_{m} & = & \displaystyle \frac{\displaystyle \int {\xi }_{m}(y)\sigma (y){dy}}{\displaystyle \int | {\xi }_{m}(y){| }^{2}{dy}}.\end{array}\end{eqnarray} \tag{ 9 }$$

The second plot in Figure 1 shows this forced response, given the sum of the free modes weighted by their response coefficients a_m. We used an equal radiative and dynamical damping rate α = 0.2 (this equality is relaxed later on), and a forcing magnitude Q₀ = 1 to match those of Matsuno (1966). This is also consistent with the range of values used in Showman & Polvani (2011) to represent a hot Jupiter.

The even parity forcing Q only excites even modes in u and h, and odd modes in v, which is π out of phase with u and h in the shallow-water equations. The Rossby and Kelvin modes dominate the forced response as they have the lowest frequencies, so are quasi-resonant with the stationary forcing with zero frequency.

Equation (12) shows that a uniform background flow can Doppler shift the frequency ω_m, changing the imaginary part of the wave response, which shifts its longitude. The third plot in Figure 1 shows how the forced response in the second plot changes when Doppler shifted by an eastward zonal flow.

These standing waves can transport momentum through the atmosphere. The zonal momentum equation predicts that the zonal acceleration depends on the gradient of the stationary wave response (Andrews & McIntyre 1976). Showman & Polvani (2011) show that this acceleration is:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \bar{U}(y)}{\partial t} & = & {\bar{v}}^{* }\left(f-\displaystyle \frac{\partial \bar{u}}{\partial y}\right)-\displaystyle \frac{1}{\bar{h}}\displaystyle \frac{\partial (\overline{{\left({hv}\right)}^{{\prime} }{u}^{{\prime} }})}{\partial y}\\ & & +\ \left(\displaystyle \frac{1}{\bar{h}}\overline{{u}^{{\prime} }{Q}^{{\prime} }}+{\overline{R}}_{u}^{* }\right)-\displaystyle \frac{{\bar{u}}^{* }}{{\tau }_{\mathrm{drag}}}-\displaystyle \frac{1}{\bar{h}}\displaystyle \frac{\partial \left(\overline{{h}^{{\prime} }{u}^{{\prime} }}\right)}{\partial t}\end{array}\end{eqnarray} \tag{ 10 }$$

where overbars are zonal means, primes are deviations from zonal means, and the "thickness-weighted zonal average of any quantity A" is ${\overline{A}}^{* }\equiv \overline{{hA}}/\overline{h}$ (Showman & Polvani 2011). R_u is the zonal component of the vertical momentum transport ${\boldsymbol{R}}$:

$$\begin{eqnarray}{\boldsymbol{R}}=\left\{\begin{array}{ll}-\displaystyle \frac{Q{\boldsymbol{v}}}{h}, & Q\gt 0\\ 0, & Q\lt 0\end{array}\right..\end{eqnarray} \tag{ 11 }$$

Showman & Polvani (2011) show why the tilted wave pattern (such as in the first plot in our Figure 6) transports eastward zonal momentum toward the equator. The tilts result in velocities "such that $u^{\prime} v^{\prime\prime} \lt 0$ in the northern hemisphere and $u^{\prime} v^{\prime} \gt 0$ in the southern hemisphere," so the term in $u^{\prime} v^{\prime}$ in Equation (10) produces an eastward acceleration.

Figure 2 shows the latitudinal profile of this jet acceleration, calculated from the forced response in the second plot of Figure 1. The black line is $\partial \bar{U}(y)/\partial t$, which is the sum of the first term (blue), second term (red), third term (purple), and fourth term (green) in Equation (10). Given our length and timescales (Matsuno 1966), the equatorial acceleration of 0.1 implies a dimensional acceleration of 0.01 ms⁻² on an Earth-like planet or 0.1 ms⁻² on a hot Jupiter. GCM simulations show that the jets of order 100 ms⁻¹ on Earth-like planets (see Section 7 and jets of order 1000 ms⁻¹ on hot Jupiters form from rest in about 20 days. The linear shallow-water model on the beta-plane therefore predicts an acceleration about 20 times too large. We will show later that this is because the beta-plane model requires an unrealistically large forcing in order to balance the height perturbation due to the jet (as the jet height perturbation scales linearly with the jet speed on the beta-plane). The model in spherical geometry gives a more realistic equatorial acceleration, as the jet height perturbation scales quadratically (correctly) with the jet speed, so is balanced by a smaller forcing Q.

Zoom In Zoom Out Reset image size

The net equatorial acceleration in Figure 2 is eastwards, but there is a westward acceleration further from the equator (Showman & Polvani 2011). Our GCM simulations do not show significant zonal-mean westward flows at the pressure level of the standing Rossby waves, although some simulations in different studies with different planetary parameters show strong westward jets (Showman et al. 2015). The westward acceleration predicted by the linear model may be opposed by an eastward acceleration from a Hadley circulation in the GCM, which is not represented in the linear model. Eastward flow may also take place in the lower or upper atmosphere of a real planet or a GCM simulation, which is not represented by our single-layer model.

The eastward flow from this acceleration affects the forced wave solutions. A zonally uniform jet modifies the projection coefficients in Equation (12) (Tsai et al. 2014), Doppler shifting the waves eastwards to a maximum of π/2 ahead of the forcing. The flow shifts the frequency ω of the free modes by ${k}_{x}\bar{U}$, which then affects the imaginary part of the coefficients projecting the free modes onto the forced solution:

$$\begin{eqnarray}&&{a}_{m}=\displaystyle \frac{1}{\alpha -i({\omega }_{m}-{k}_{x}\bar{U})}.\end{eqnarray} \tag{ 12 }$$

A positive eastward flow $\bar{U}$ shifts the wave response eastwards. The third plot in Figure 1 shows how a uniform flow shifts the maximum of the forced response toward +90° longitude. For a significant shift, the frequency shift ${k}_{x}\bar{U}$ must be on the order of ω_m, i.e., the jet speed $\bar{U}$ must be on the order of ω_m/k_x. The frequencies are the solutions of ${\omega }^{2}-{k}_{x}^{2}-{k}_{x}/\omega \,=2n+1$ (Matsuno 1966) for zonal wavenumber k_x (k_x = 1 in our periodic system) and mode n. So, the Doppler shift of the n = 1 Rossby wave is significant when ${k}_{x}\bar{U}\approx 0.25$. The scale of the zonal wavenumber of these planetary waves is set by the planet's radius: k ∼ 1/a ∼ 1. The dimensional velocity is $[U]=c\,={({gH})}^{1/2}$, so the critical jet speed at which ${\omega }_{m}\sim {k}_{x}\bar{U}$ and the Rossby mode is significantly shifted is ${\bar{U}}_{\mathrm{crit}}={({gH})}^{1/2}/4$.

For an Earth-like planet, g ∼ 10 ms⁻² and h_e ∼ 5 km, so ${\bar{U}}_{\mathrm{crit}}=50\,{\mathrm{ms}}^{-1}$. GCM simulations show equilibrium jet speeds of order 100 ms⁻¹, so the Doppler shift should be significant. A hot Jupiter has g ∼ 20 ms⁻² and h_e ∼ 200 km (Showman & Polvani 2011), so ${\bar{U}}_{\mathrm{crit}}=500\,{\mathrm{ms}}^{-1}$. GCMs and observations show equilibrium jet speeds of order 1000 ms⁻¹, so the Doppler shift should again be significant. In Section 4, we will linearize the shallow-water equations about a nonuniform jet.

In this section, we discuss the effect of the shear flow on the free modes of Equation (14), in order to understand changes to the forced response in the next section. In Section 3, we showed how the forced solution to the shallow-water equation in a uniform background flow (Equation (5)) can be expressed as a sum of the free modes using Equation (9). It is not possible to express the forced solution to Equation (15) in a shear background flow using such a simple sum of the free solutions to Equation (14). However, we can still use the free solutions of Equation (14) to qualitatively understand the forced solutions that we will find numerically in Section 4.2.

We write the free solutions to the shallow-water equations as a complex function of latitude A(y). Later, we will write and plot the forced solutions as functions of latitude and longitude, in the form $A(y){e}^{i\delta (y)x}$. The phase shift δ(y) in a shear flow is equivalent to the phase shift $({\omega }_{m}-{k}_{x}\bar{U})$ in a uniform flow (Equation (12)). The sign of the eigenvalue ω_m of any mode in shear determines where the free mode appears in the forced solution, as in the previous section for uniform background flow.

The free modes of the beta-plane shallow-water system (Equation (14)) depend on the background flow and height fields. For zero background flow, the free modes are the same as those discussed in Section 3. For an analytic solution with a uniform background flow ${\bar{U}}_{0}$, the free modes are linearly Doppler shifted as discussed in Section 3.3 (Tsai et al. 2014). For a shear background flow $\bar{U}(y)$, we found the free modes of Equation (14) using the method described in the Appendix, with the parameters listed at the start of Section 4 (and equal radiative and damping rates, α_rad = α_dyn).

Figure 3 shows the real parts of the eigenvalues of the free Kelvin mode and the symmetrical free Rossby modes of Equation (14), for a background flow $\bar{U}(y)={U}_{0}{e}^{-{y}^{2}/2}$ with a variable magnitude U₀. These are the lowest-order modes excited by the symmetric forcing. The Kelvin mode already has a positive eigenvalue for ${U}_{0}=0$, which increases and shifts further east as the flow increases, moving the Kelvin mode toward a maximum shift of +90° and contributing to the hotspot shift.

Zoom In Zoom Out Reset image size

Tsai et al. (2014) suggests that in a uniform background flow the n = 1 Rossby mode is shifted eastward of the substellar point, adding to the hotspot shift. Figure 3 shows that in this shear flow $\bar{U}(y)$, the n = 1 Rossby mode eigenvalue increases but does not become positive for U₀ = 1.0 (corresponding to the jet speed in the forced response plotted later in Figure 6). This means that it is shifted eastwards in the forced response, but does not reach the east of the substellar point. However, the higher-order Rossby mode eigenvalues do become positive, shifting to the east of the substellar point and contributing to the hotspot shift. So, the forced response and the hotspot shift are affected by the higher-order Rossby modes, and not just dominated by the Kelvin and n = 1 Rossby modes.

That is not to say that the n = 1 mode is never responsible for the hotspot shift—later, we will show that in a spherical geometry the n = 1 mode shifts close to +90° eastwards. It is also possible in the beta-plane system for different input parameters (flow speed, damping rates) to shift the n = 1 Rossby mode past the substellar point. But, our free mode expansion has shown that the n = 1 Rossby mode is not the only important mode, and that the higher-order modes are also important to the forced response.

For zero damping, half of these eigenvalues will have positive imaginary parts, and the modes corresponding to them will grow exponentially. Nonzero damping decreases the imaginary part of all the modes, so will make some or all of these modes stable. In general, the free linear system in Equation (14) will have some unstable modes unless the damping is very large. These unstable modes are similar to those discussed by Wang & Mitchell (2014), who show how similar modes can produce superrotation even on a planet without a permanent day–night heating difference.

These unstable modes technically make the linear forced wave problem ill-posed, since the result of any linear initial value problem will be eventually dominated by the most rapidly growing modes rather than the stationary response. Later comparison with nonlinear GCM simulations in Section 7 will show that the forced response still has considerable explanatory power. This may be because in reality the unstable modes equilibrate due to damping or nonlinear effects, at a sufficiently low amplitude that they take the form of mobile waves propagating across the forced stationary pattern without significantly disrupting its basic structure. Future work should investigate the exact nature of these instabilities, and the effect of damping and shear flow on their growth rates.

The shear flow also affects the latitudinal structure A(y) of the modes. The lowest-order free solutions of Equation (14) (the Kelvin and Rossby modes), plotted in Figures 4 and 5, resemble the free solutions with zero shear flow (Matsuno 1966), with their latitudinal structure slightly changed by a weak shear flow $\bar{U}=0.1{e}^{-{y}^{2}/2}$. The shear flow perturbs the solutions by adding higher-order meridional structure. For example, the meridional wind of the Rossby wave in Figure 5 resembles the n = 1 parabolic cylinder function added to the n = 3 function (see Figure 19 in Appendix A.3). Boyd (1978) discusses how a shear flow affects the meridional structure of these modes in more detail.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we discuss solutions of the forced linear shallow-water equations with ω = 0 (Equation (15)) as an inhomogeneous forced linear problem.

Figure 6 shows the main result of this paper: how linearizing the equations about a background jet changes the response to stationary forcing. The first plot in the figure is the case without a jet—the exact forced solution of Matsuno (1966). The case with a jet is our new solution to Equation (15), linearized about a background flow $\bar{U}(y)$ and $\bar{H}(y)$. Both plots were calculated with the method in the Appendix, but the first plot has identical form to the analytic solution in zero flow (the second plot in Figure 1), as discussed in Appendix A.3.

Zoom In Zoom Out Reset image size

We suggest that the background shear flow changes the solution in two key ways. First, the forced response is Doppler shifted eastwards as discussed previously. Second, the jet velocity and height adds to the forced solution, adding to the height along the equator, magnifying the cold Rossby lobes, and combining with the shifted stationary waves to form the hotspot.

In Section 4.1 we showed how the Doppler shift from the jet affects the Kelvin mode and the Rossby modes to different extents. Comparing Figure 6 to Figure 3 shows how the position of each mode in the forced response depends on its eigenvalue. The Kelvin mode has the most positive eigenvalue in Figure 3, so produces a large hotspot shift on the equator.

Unlike the Kelvin mode, the eigenvalue of the n = 1 Rossby mode is initially negative. It increases as the flow increases, but does not become positive for U₀ = 1, so does not shift past the substellar point. This means that in Figure 6 the hotspot shift is smaller at high latitudes. The higher-order Rossby modes are shifted further east, contributing more to the hotspot shift. However, the modes get progressively weaker at higher orders so affect the total response less. This shift of all the wave modes is dominated by the Kelvin mode and lowest-order Rossby mode. It makes the position of the hotspot and the cold Rossby lobes match our GCM output in Section 7.

The solutions in this section are limited by the linear approximation of the Coriolis parameter and non-dimensionalized y-coordinate on the beta-plane. This makes them inaccurate at high latitudes, and difficult to compare directly to real planets. Solutions on the beta-plane in other studies are similar to solutions in spherical coordinates (Showman & Polvani 2011; Heng & Workman 2014), which suggests that its limitations are not too problematic.

The beta-plane system also requires a forcing with the same scale L_R as the scale of the y dimension for a simple solution (Matsuno 1966). We found that the solution was not well represented by the parabolic cylinder functions if these scales were not similar. This limits our solutions to tidally locked planets where the scale of the forcing (the planetary radius) is comparable to the the equatorial Rossby radius L_R. This is not appropriate for exoplanets such as tidally locked planets with similar size and rotation rate as Earth. In Section 5, we will solve the linear shallow-water equations in a shear flow on a sphere, which relaxes some of the constraints of the beta-plane.

In this section, we show how the zonal acceleration of the linear shallow-water system, given by Equation (10), decreases as the zonal flow increases. This explains why the system should reach an equilibrium rather than accelerating the zonal flow indefinitely.

Figure 2 shows the different physical sources making up the zonal acceleration. The equatorial acceleration is primarily controlled by a balance of acceleration due to the convergence of zonal momentum (the second term in Equation (10)), deceleration due to vertical momentum transport (the third term), and deceleration due to drag on the zonal flow (the fourth term). When the mean zonal flow is zero, Figure 2 shows that there is a positive eastward acceleration on the equator that accelerates the flow initially (Showman & Polvani 2011).

As the mean eastward zonal flow $\bar{U}(y)$ increases (we model this as an eastward jet with Gaussian shape despite the westward acceleration at high latitudes, as discussed previously), the equatorial eastward acceleration in our linear model decreases. Showman & Polvani (2011) suggested that the shallow-water system should reach equilibrium when the deceleration due to drag on the zonal flow increases sufficiently to balance the other terms. Tsai et al. (2014) suggested that in addition to the increased drag, the acceleration from zonal momentum convergence decreases as the zonal flow shifts the forced response eastwards.

Our solutions also show this decrease in acceleration. Figure 6 shows that when the background flow is zero, the Rossby and Kelvin components of the forced response are out of phase (i.e., at different longitudes) (Matsuno 1966), so the part h'u' of the zonal momentum convergence term is nonzero. The magnitude of this term depends on the phase difference between the modes (Tsai et al. 2014), and corresponds to the tilt of the forced u, v, and h fields (Showman & Polvani 2011). As the background jet flow increases in the forced linear solutions, the Rossby and Kelvin components (and other higher-order terms) shift eastwards, tending toward +90°. So as the zonal flow increases, the waves become more in phase with each other, and the acceleration due to zonal momentum convergence decreases.

Figure 7 shows the acceleration due to zonal momentum convergence at three different background zonal flow speeds. As the background flow increases, the equatorial acceleration from this term decreases. This combines with the increased drag $-{\bar{u}}^{* }/{\tau }_{\mathrm{drag}}$ on the mean zonal flow to give a net equatorial acceleration of zero at high enough flow speeds, in agreement with Tsai et al. (2014).

Zoom In Zoom Out Reset image size

With the length and timescales on the beta-plane discussed previously, a non-dimensional velocity of U₀ = 1 corresponds to a velocity of order 100 ms⁻¹ on an Earth-like planet, and to a velocity of order 1000 ms⁻¹ on a hot Jupiter. These are comparable to the velocities seen in GCM simulations, so the suggestion of Tsai et al. (2014) that equilibrium is reached when the flow and Doppler shift becomes great enough is consistent with our linear shallow-water model.

In the previous section we assumed that the radiative damping ${\alpha }_{\mathrm{rad}}$ and dynamical damping term α_dyn were equal. This allows an analytic solution but is not physically justified. Radiative damping in the shallow-water equations has a clear physical basis, but linear dynamical damping does not, apart from effects like eddy viscosity, magnetohydrodynamic damping, or extra nonlinear terms (Heng & Workman 2014). In this section, we suggest that although the damping rates could be very different in reality, this does not alter the overall form of the solution.

Figure 8 shows the solution with zero dynamical damping. In this first plot, which has no background flow, this appears to weaken the Kelvin wave response relative to the Rossby wave response (compared to Figure 6). The Rossby lobes also tilt in the opposite direction to Figure 6, matching Figure 6(c) in Showman & Polvani (2011). In the second plot, which has a background jet, this effect is less significant as the jet velocity and height dominate on the equator where the Kelvin response would be. So, the dynamical damping is not crucial to the form of the solution, and it is not a large problem that it is not a well-defined quantity. In summary, once a strong jet forms, the form of our linear shallow-water solutions is the same regardless of the exact values of the damping α_rad and α_dyn.

Zoom In Zoom Out Reset image size

Understanding the mechanism behind the global circulation and eastward hotspot shift observed on many tidally locked planets (Parmentier & Crossfield 2017) is vital to interpreting observations of such planets, particularly why the magnitude of the hotspot shift varies. Showman & Polvani (2011) suggest that the hotspot shift is caused by "eastward group propagation of Kelvin waves," as the free Kelvin waves in the shallow-water model propagate eastwards (opposite to the Rossby waves) (Matsuno 1966). This is similar to the non-periodic Gill model, where the forced response propagates eastwards along the equator (Gill 1980). Other explanations have represented the wave dynamics indirectly. For example, Komacek & Showman (2016) and Zhang & Showman (2017) use a linear model based on balancing an advective timescale against a radiative timescale to explain the day–night contrast and hotspot shift on tidally locked planets. This model predicted a hotspot shift on the equator, but did not represent wave dynamics directly so did not include their effect on the hotspot shift, or predict the height field off the equator.

Tsai et al. (2014) focused on the effect of waves using a linear shallow-water model on a beta-plane with uniform flow ${\bar{U}}_{0}$, and suggested that the hotspot shift is instead caused by "the eastward shift of the Rossby wave with almost no zonal phase shift of the Kelvin wave." Our linear model builds on Tsai et al. (2014) by using a nonuniform flow $\bar{U}(y)$ and and a height perturbation $\bar{H}(y)$ in balance with this flow. The height perturbation and sheared flow are crucial to the overall form of the circulation seen in GCM results such as Figure 14.

The forced linear solution in shear in Figure 6 can be approximated as the analytic solution in Section 3, modified by the eastward equatorial jet in two ways:

• 1.  
  The flow Doppler shifts the wave response eastwards (see Section 3)
• 2.  
  The flow is balanced by a zonally uniform height perturbation $\bar{H}(y)$

Figure 9 shows a first-order model of these effects, using the analytic solutions of the shallow-water equations in a uniform flow. It shows how the jet height perturbation completes the global circulation pattern, particularly the hotspot pattern. The hotspot in the second plot in Figure 9 (without the jet height perturbation) resembles a Rossby wave. Summed with the height perturbation in the third plot, it becomes centered on the equator and matches the form of the full solution in Section 4 and the GCM simulations in Section 7. We therefore agree with Tsai et al. (2014) that the hotspot shift is caused by an eastward Doppler shift of the stationary Kelvin and Rossby waves toward +90°, rather than eastward propagating Kelvin waves or pure heat transport from air in the eastward flowing jet.

Zoom In Zoom Out Reset image size

The analytic solutions in Figure 9 are approximately the first-order terms in the pseudo-spectral expansion in Appendix. The pattern in Figure 6 is more complex than the sum of patterns in Figure 9 as the flow $\bar{U}(y)$ is not constant in y, but the eastward shift is still the most important effect. The shear flow introduces a tilt to the stationary Rossby waves, which we suggest is due to the faster flow at the equator producing a stronger Doppler shift there.

The linear model in Showman & Polvani (2011) does not include a background zonal flow, which may be why it matches their GCM in spin-up but not in equilibrium. Showman & Polvani (2011) also compare the linear model to numerical solutions of the nonlinear shallow-water equations, which should include the effect of any background flow. The nonlinear solutions in Figure 8 of Showman & Polvani (2011) do have a mean zonal flow but it appears that either the damping is too strong, or the jet is too weak or narrow, to produce a hotspot shift. In their case (a), the maximum zonal speed is approximately 100 ms⁻¹, which we suggest in Section 3 is not large enough for a significant phase shift on a hot Jupiter. Their case (b) has a maximum zonal wind of approximately 700 ms⁻¹, which is large enough for a phase shift, but the jet is so narrow that this shift only happens very close to the equator, shown by the narrow eastward "bump" in the hotspot in the top right panel of their Figure 8. If the cases are sufficiently strongly damped, then a zonal flow may still not produce a hotspot shift. Equation (12) shows that the projection coefficient will be mostly real if α is much larger than ${k}_{x}\bar{U}$, so the maximum of the wave response will be in phase with the forcing (at the substellar point, with no hotspot shift).

The pseudo-spectral solutions to Equation (15) such as in Figure 6 show the forced response of a shallow-water system with particular parameters such as damping rate and jet speed, but do not provide a direct relation between these parameters and observable quantities such as the position of the atmospheric hotspot. In this section, we calculate the on-equator height response on the beta-plane, where β = 0 so we can ignore the effects of waves in the linear shallow-water equations. This gives a simpler 1D model, which predicts a direct scaling relation between the equatorial height maximum and the planetary parameters.

For y = 0 (on the equator), retaining only damping and advection terms, the third line of Equation (13) simplifies from

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {H}^{{\prime} }u}{\partial x}+{H}^{{\prime} }\displaystyle \frac{\partial v}{\partial y}-y\bar{U}(y)v+\displaystyle \frac{\partial \bar{h}}{\partial t}\\ +\ \alpha h+\displaystyle \frac{\partial \bar{U}(y)h}{\partial x}={Q}_{0}{e}^{-{y}^{2}/2}\end{array}\end{eqnarray} \tag{ 20 }$$

to

$$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{h}}{\partial t}+\alpha h+\displaystyle \frac{\partial {U}_{0}h}{\partial x}={Q}_{0}.\end{eqnarray} \tag{ 21 }$$

Then as before, we set $\partial /\partial t=0$ and $\partial /\partial x={{ik}}_{x}$:

$$\begin{eqnarray}&&-\alpha h(y=0)+{{ik}}_{x}{U}_{0}h(y=0)={Q}_{0}.\end{eqnarray} \tag{ 22 }$$

So then the real part of the full solution $h(x,y=0)\,=h(y=0){e}^{{{ik}}_{x}x}$ is:

$$\begin{eqnarray}&&h(x,y=0)=\displaystyle \frac{{Q}_{0}}{{\alpha }^{2}+{k}_{x}^{2}{U}_{0}^{2}}(-\alpha \cos ({k}_{x}x)+{k}_{x}{U}_{0}\sin {k}_{x}x).\end{eqnarray} \tag{ 23 }$$

This corresponds to a sinusoidal height perturbation on the equator. The magnitude of this height perturbation is approximately:

$$\begin{eqnarray}&&{h}_{0}=\displaystyle \frac{\alpha +{k}_{x}{U}_{0}}{{\alpha }^{2}+{k}_{x}^{2}{U}^{2}}{Q}_{0}.\end{eqnarray} \tag{ 24 }$$

The position of the hotspot x₀ is the maximum of this height perturbation, where $\partial h/\partial x=0$:

$$\begin{eqnarray}&&{x}_{0}=\displaystyle \frac{1}{{k}_{x}}{\tan }^{-1}({k}_{x}\displaystyle \frac{{U}_{0}}{\alpha }).\end{eqnarray} \tag{ 25 }$$

Note that this is the same as the approximation to the hotspot shift from Zhang & Showman (2017), ${\lambda }_{s}={\tan }^{-1}(\tfrac{{\tau }_{\mathrm{rad}}}{{\tau }_{\mathrm{adv}}})$, if we set τ_rad = 1/α and τ_adv = k_x/U₀. This is because their prediction for the hotspot shift is calculated from Equation (41) in Zhang & Showman (2017):

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}+\displaystyle \frac{1}{{\tau }_{\mathrm{adv}}}\displaystyle \frac{\partial T}{\partial \lambda }=\displaystyle \frac{{T}_{\mathrm{eq}}(\lambda )-T}{{\tau }_{\mathrm{rad}}}.\end{eqnarray} \tag{ 26 }$$

This is equivalent to our Equation (22), if we set $T\equiv h$, τ_adv = k_x/U₀, $\lambda \equiv x$, and $({T}_{\mathrm{eq}}(\lambda )-T)/{\tau }_{{rad}}={Q}_{0}$. Zhang & Showman (2017) use a more complex day–night forcing difference than our Q(x, y) to solve this equation, which is why our prediction for the hotspot shift only matches the approximate solution to their equation.

This 1D model predicts a hotspot shift varying from 0°–90° east of the substellar point. This is similar to the linear shallow-water model of Tsai et al. (2014), where the uniform background flow could shift the maximum of the wave response to a maximum of 90° east. With zero zonal flow, the 1D model predicts a hotspot shift at 0°, the substellar point. This is only the case in our 2D linear solutions if the damping is strong, i.e., if the damping rate α is much greater than the Rossby and Kelvin wave frequencies. At higher zonal flow U₀ the 1D model predicts a hotspot shift approaching 90°. This can agree with the 2D model—see Figure 12, where the second plot has very low damping—but in general, wave terms and damping terms prevent such a large hotspot shift (as in Figure 6).

So, the 1D model balancing advection and damping is useful for intuition and basic scaling relations, and can predict the hotspot shift of a tidally locked planet to first order, but is not accurate when wave effects become important. The 2D model is more useful when considering the full disk-integrated phase curve, where the off-equator response is very different to the 1D one-equator model.

In the previous section, we used a 1D model to calculate simplified predictions of how observables such as hotspot shift and day–night contrast should scale with planetary parameters such as damping rate. In this section, we will test the predictions from the 1D model and our discussion in previous sections.

The solutions such as those in Figure 6 used a representative set of non-dimensional parameters, which roughly matched a typical Earth-sized tidally locked planet or hot Jupiter. In this section, we identify the free parameters of the beta-plane solution and the solution in spherical coordinates, and vary these parameters to test the predictions of the previous section.

The free parameters on the beta-plane are the background jet strength and width U₀ and y₀, the forcing strength Q₀, and the radiative and dynamical damping rates α_rad and α_dyn. The solution in spherical geometry also depends on these parameters. The jet and forcing strengths U₀ and Q₀ affect the relative magnitude of the jet and wave responses without changing the actual wave solutions itself, as we show below for the case in spherical geometry. Figure 11 shows how changing the relative magnitude of the forcing Q₀ and the jet U₀ changes the solution in spherical coordinates. For strong forcing, the wave response dominates and the height field is very asymmetric, resulting in a large day–night contrast. For weak forcing or a strong jet, the jet velocity and height field dominates, leading to a more zonally homogeneous circulation.

Zoom In Zoom Out Reset image size

Section 6.1 predicts that increasing the damping rate α should decrease the hotspot shift and decrease the day–night contrast. Figure 12 shows that varying the damping α has this effect in the full solutions (both α_dyn and α_rad are set to be equal here). We can interpret this effect using Equation (12)—increasing the damping increases the real part of the projection coefficient, so the maximum of the forced solution is more in phase with the maximum of the forcing. The damping also affects the relative strength of the Rossby and Kelvin parts of the forced response. For strong damping, the Kelvin part dominates, leading to the single maximum centered on the equator. For weak damping, the Rossby part dominates, leading to the two maxima off the equator.

Zoom In Zoom Out Reset image size

The linear shallow-water equations in a spherical geometry described in Section 5 also have the additional parameter $G={(c/a{\rm{\Omega }})}^{2}={gH}/{a}^{2}{{\rm{\Omega }}}^{2}$. This extra parameter is due to the extra degree of freedom added by constraining the latitudinal direction, unlike on the beta-plane.

The additional free parameter G in the spherical geometry describes the effect of rotation rate, planetary radius, and gravity wave speed. It is the square of the ratio of the Rossby radius of deformation to the planetary radius, and as such represents the relative latitudinal scale of the waves in the system to the size of the planet. It is also equivalent to the square of the Weak Temperature Gradient (WTG) parameter, which Pierrehumbert & Ding (2016) used to characterize the global circulation of tidally locked planets.

Equation (17) shows that G is inversely proportional to the background jet height perturbation. Figure 13 shows that for $G\gg 1$ (a low rotation rate Ω) waves dominate the forced response, giving a large day–night contrast. For $G\ll 1$ the jet velocity and height dominate, giving a more zonally symmetric system. Our description of the effect of G matches the qualitative predictions of previous studies, which interpreted changes in atmospheric dynamics on these planets by comparing the equatorial Rossby radius and the planetary radius (Carone et al. 2015; Koll & Abbot 2015; Showman et al. 2015).

Zoom In Zoom Out Reset image size

This 2D linear model does not give such a straightforward prediction of the hotspot shift and day–night contrast as the 1D model, but does capture the important off-equator behavior and the effect of the planetary waves. It would be ideal for comparison with observational methods that retrieve 2D maps of the planets (Rauscher et al. 2018).

The plots in Figure 15 show the time-averaged eddy streamfunction and eddy velocities of our simulations of the planet discussed above, with a 10 day orbital period. In Sections 3 and 4, we showed how the forced planetary waves are Doppler shifted eastwards in the linear shallow-water model. This is clearest in the pattern of the eddy velocities in the GCM results in Figure 15. The pattern of velocities is 180° out of phase with the second plot in Figure 1, but is in phase with the third plot, which has been Doppler shifted eastwards by a background flow.

Zoom In Zoom Out Reset image size

This explains why in Figure 14 the equatorial u velocity is lower at about +90°, as the local Rossby wave velocity opposes it, and higher at about −90°, as the local Rossby wave velocity adds to it. Studies such as Carone et al. (2015) and Pierrehumbert & Hammond (2019) have shown this pattern of Rossby wave velocities, which differs from the Rossby wave velocities in the non-shifted linear shallow-water solutions (Matsuno 1966; Showman & Polvani 2011).

The linear shallow-water model predicts that the Rossby and Kelvin components of the forced response should shift eastwards as the eastward zonal flow increases, up to a maximum of +90°. In this section, we test if this is the case in the spin-up of our GCM. We initialized a tidally locked planetary atmosphere from rest with the parameters listed previously (idealized Earth-like conditions, with a 10 days orbital period). To compare the GCM results to the linear shallow-water results, we decomposed its daily output into Kelvin and Rossby components, following the method of Boulanger & Menkes (1995).

Figure 16 shows the Kelvin and Rossby components of the GCM height field at 0, 30, and 60 days after spin-up, during which time the mean zonal flow had accelerated from rest to a maximum of approximately 90 ms⁻¹ eastwards on the equator. The initial Kelvin and Rossby waves are centered at the substellar point, rather in the western hemisphere as might be expected for a Matsuno-type pattern in zero background flow. We suggest this is because the radiative and dynamical damping is strong enough in the GCM that the Kelvin and Rossby waves are mostly in phase with the stellar forcing in the absence of background flow. This is shown by Figure 3 of Showman & Polvani (2011), where the strongly damped Matsuno-type solutions have their maxima close to the substellar point. Figure 10 of Showman & Polvani (2011) is a snapshot of their GCM during its spin-up phase before the jet forms, which also has a Rossby wave maximum at the substellar point, rather than to its west.

Zoom In Zoom Out Reset image size

Although the initial position of the waves is not exactly as expected, Figure 16 shows that as the jet forms, the Rossby and Kelvin waves shift eastwards and reach an equilibrium at about toward +60° east. This matches the prediction of our linear shallow-water model.

Figures 17 and 18 show the mean mid-atmosphere temperature of a suite of GCM tests with variable stellar forcing and planetary rotation rate. Studies such as Carone et al. (2015), Showman et al. (2015), Haqq-Misra et al. (2018), and Pierrehumbert & Hammond (2019) varied the parameters of simulated tidally locked planets and identified changes in the global circulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We suggest that we can qualitatively explain these circulation patterns with our linear shallow-water solutions, using the ideas discussed in Section 6. Figure 17 shows the effect of rotation rate. At high rotation rates, the jet height (and temperature perturbation) dominates as it is proportional to Ω via G, as shown in Figure 13. This gives a more zonally homogeneous circulation. At high jet speed the wave response will be more Doppler shifted and the hotspot shift will be larger. This gives a smaller day–night contrast and larger hotspot shift.

Second, the effect of stellar forcing in Figure 18. At high forcing (instellation) the wave response dominates as shown in Figure 11, so the wave response and the day–night contrast are large. The temperature is also higher, so the radiative damping is stronger and the wave response moves in phase with the forcing (see Section 3). This leads to a smaller hotspot shift as well as the large day–night contrast—i.e., a zonally inhomogeneous circulation in phase with the stellar forcing, with strong temperature gradients.

These results match the scaling relations of the 1D and 2D models discussed in Section 6. There, we predicted that planets with higher damping α would have a smaller hotspot shift and larger day–night contrast. This is the case in Figure 18, where the hotter planet has a higher radiative damping rate. We also predicted that planets with a higher rotation rate should have a more zonally homogeneous circulation pattern, which we see in Figure 17.

These results are similar to the work of Komacek & Showman (2016), Komacek et al. (2017), and Zhang & Showman (2017), who predicted scaling relations for observables such as hotspot shift and day–night contrast, based on balancing jet transport against radiative and dynamical damping. Our wave-based approach to the circulation predicts some of the same scaling relations for different reasons. Zhang & Showman (2017) predicted that an atmosphere with a low heat capacity and short radiative timescale has a large day–night contrast and small hotspot shift, which makes sense as the hot air transported in their jet model will cool quickly on the nightside.

We make the same prediction, but for a different reason. A short radiative timescale t_rad corresponds to a strong damping term α_rad = 1/t_rad. Figure 12 shows how this weakens the forced wave response and brings the response in phase with the forcing (as discussed in Section 3), reducing the hotspot shift. For a long radiative timescale and weak radiative damping α_rad, the Doppler shift due to the jet will dominate the effect of the forcing in Equation (12) and the hotspot shift will be large.

Some tidally locked planets may not match the linear shallow-water model. Very slowly rotating planets are dominated by an overturning day–night circulation without an eastward jet. They may be limited by nonlinear constraints on geostrophic adjustment (Pierrehumbert & Hammond 2019), and the linear theory in this paper will not apply.

Boyd (1978) solves the linearized shallow-water equations, by reducing them to a single equation for a single variable, and applying the pseudo-spectral method. In this paper, we solve the entire system of shallow-water equations at once with the method in Appendix A.2, but explain the method for a single equation here as it naturally leads to the second method (Boyd 2000).

The matrix elements H_ij in Equation (29) are evaluated using the operator L at the collocation points x_i and for every mode ϕ_j, and the vector elements f_i are the terms q evaluated at the collocation points x_i:

$$\begin{eqnarray}&&{H}_{{ij}}=L{\phi }_{j}({x}_{i})\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{f}_{i}=q({x}_{i}).\end{eqnarray} \tag{ 31 }$$

This is then solved using a standard linear algebra routine to find a_n, and the solution u(x) is reconstructed using Equation (28).

The pseudo-spectral method can also be applied to systems of linear ordinary differential equations. For a system of forced, time-independent equations:

$$\begin{eqnarray}&&{\boldsymbol{Lu}}={\boldsymbol{q}}.\end{eqnarray} \tag{ 32 }$$

The condition that the differential equation is satisfied at the collocation points gives the equivalent matrix equation to Equation (29):

$$\begin{eqnarray}&&{\boldsymbol{Ha}}={\boldsymbol{f}}.\end{eqnarray} \tag{ 33 }$$

${\boldsymbol{H}}$ is an M × N square matrix with elements:

$$\begin{eqnarray}&&{H}_{{ij}}^{{kl}}={L}^{{kl}}{\phi }_{j}({x}_{i}),\end{eqnarray} \tag{ 34 }$$

i.e., the operator L^kl which acts on the lth variable in the kth equation, applied to the jth basis function and evaluated at the ith collocation point. ${\boldsymbol{f}}$ is a vector made up of N subvectors f_i, which are the forcing terms in each equation evaluated at each collocation point.

$$\begin{eqnarray}{\boldsymbol{H}}=\left(\begin{array}{cc}{\left(\begin{array}{cc}{H}_{{ij}} & ...\\ \vdots & \ddots \end{array}\right)}^{{kl}} & ...\\ \vdots & \ddots \end{array}\right)\left(\begin{array}{c}\left(\begin{array}{c}{\alpha }_{i}\\ \vdots \end{array}\right)\\ \vdots \end{array}\right)=\left(\begin{array}{c}\left(\begin{array}{c}{f}_{i}\\ \vdots \end{array}\right)\\ \vdots \end{array}\right).\end{eqnarray} \tag{ 35 }$$

${\boldsymbol{H}}$ is the same as the matrix in Equation (30) with the elements H_ij replaced by submatrices H^kl_ij. Solving this system returns the coefficients of the basis functions, and the solutions are:

$$\begin{eqnarray}&&u(y)=\displaystyle \sum _{n=0}^{N}{a}_{n}{\phi }_{n};\quad v(y)=\displaystyle \sum _{n=0}^{N}{b}_{n}{\phi }_{n};\quad h(y)=\displaystyle \sum _{n=0}^{N}{c}_{n}{\phi }_{n}.\end{eqnarray} \tag{ 36 }$$

This gives a linear matrix equation with one solution corresponding to the coefficient vectors a_n, b_n, and c_n of the forced solution.

Without forcing, the shallow-water equations define an eigensystem where the eigenvalue is the frequency ω.

$$\begin{eqnarray}&&{\boldsymbol{Lu}}=\omega {\boldsymbol{Pu}}.\end{eqnarray} \tag{ 37 }$$

The pseudo-spectral equation is then:

$$\begin{eqnarray}&&{\bf{Ha}}=\omega {\bf{Ra}}.\end{eqnarray} \tag{ 38 }$$

${\boldsymbol{R}}$ is an M × N square matrix with elements:

$$\begin{eqnarray}&&{R}_{{ij}}^{{kl}}={P}^{{kl}}{\phi }_{j}({x}_{i}),\end{eqnarray} \tag{ 39 }$$

i.e., the eigenvalue operator P^kl acting on the lth variable in the kth equation, applied to the jth basis function and evaluated at the ith collocation point.

$$\begin{eqnarray}{\boldsymbol{H}}=\left(\begin{array}{cc}{\left(\begin{array}{cc}{H}_{{ij}} & ...\\ \vdots & \ddots \end{array}\right)}^{{kl}} & ...\\ \vdots & \ddots \end{array}\right)\left(\begin{array}{c}\left(\begin{array}{c}{\alpha }_{i}\\ \vdots \end{array}\right)\\ \vdots \end{array}\right)=\omega \left(\begin{array}{cc}{\left(\begin{array}{cc}{R}_{{ij}} & ...\\ \vdots & \ddots \end{array}\right)}^{{kl}} & ...\\ \vdots & \ddots \end{array}\right)\left(\begin{array}{c}\left(\begin{array}{c}{\alpha }_{i}\\ \vdots \end{array}\right)\\ \vdots \end{array}\right).\end{eqnarray} \tag{ 40 }$$

This gives an eigenvalue matrix equation, with N eigenvalues and eigenvectors, corresponding to the frequencies and coefficient vectors a_n, b_n, and c_n for each free mode. Not all N modes must be physically realistic, so we identify the spurious modes by inspecting the eigenvalues for different values of N.

We use the parabolic cylinder functions ψ_n(y) (Showman & Polvani 2011) as defined in Equation (41) as a basis set for the pseudo-spectral method on the beta-plane (Equation (15)), as they are the exact free solutions of Matsuno (1966) and Boyd (2000).

Their collocation points are at their zeros (which are just the zeros of the Hermite polynomials H_n). Figure 19 shows the first few parabolic cylinder functions.

$$\begin{eqnarray}&&{\psi }_{n}(y)={e}^{-{y}^{2}/2}{H}_{n}(y)\end{eqnarray} \tag{ 41 }$$

Zoom In Zoom Out Reset image size

Figure 20 shows the magnitude of the coefficients (Equation (36)) of the pseudo-spectral solution of the shallow-water equations linearized about a jet on a beta-plane (plotted in Figure 6). The first plot shows that when the background jet flow is zero, only modes up to n = 2 are nonzero. This is the analytic solution from Matsuno (1966), which the pseudo-spectral method identifies because we have used the free modes (the parabolic cylinder functions) as our basis functions.

Zoom In Zoom Out Reset image size

For nonzero jet speed (corresponding to Figure 6), the pseudo-spectral series solution does not terminate, but the coefficients for the 30th mode are about eight orders of magnitude smaller than the largest mode. The beta-plane solutions in this paper were all calculated with at least 30 modes.

We use the Legendre polynomials as a basis set for the pseudo-spectral method in a spherical geometry (Equation (16)). Figure 19 shows the first few Legendre polynomials. Our collocation points are the zeros of these functions.

As discussed in Section 5, Equation (16) has a singularity at the the poles, which we avoided by using a rescaled height γ, where $\gamma =h/\cos \phi$ (Iga & Matsuda 2005). We replaced h with $\gamma \cos \phi$ in Equation (16), solved as normal, then multiplied the solution for γ by $\cos \phi$ to recover the solution for h.

Figure 21 shows how rescaling the h variable made the solutions converge much more quickly. In fact, the solutions without a rescaled h variable never reached a smooth solution at the poles.

Zoom In Zoom Out Reset image size